# R-code computing the triangulation probability by Monte-Carlo simulations # as described in # Triangulation Properties of the Target Area in Wireless Sensor # Networks # by Xiaoyun Li, David K. Hunter and Sergei Zuyev # Author: Sergei Zuyev # Last modified: 2008-09-24 # # Typical use: source this file in R-session: # source("triang-simul.R") # To get estimate for a sequence of values from 0.1 to 1 by, say, 0.2 # of the density lambda=1 of sensors, run # res=tri.prob(1) # reported results are saved in the list res: type # res # By default only N=10^3 simulations is performed. If you want, say N=10^4 # res=tri.prob(1,10^4) # Reported are: # res$prob - estimated triangulation probability # res$sim - number of simulations, i.e. N # res$hits - number of runs when 0 was triangulated, i.e. res$prob=res$hits/N # res$cl.n - number of runs when triangulation involved the closest to 0 node xi_0 # res$sm.n - number of runs when triangulation involved the node with the # smallest polar angle called xi_1 in the paper # res$b.n - number of runs when triangulation involved both xi_0 and xi_1 # This is simulated version of the analytic bound. # # Simulation is long, so to compute for a range of lambda and save results in a file, # say, res.dat , do: # # lam=seq(0.1,1,by=0.1) # N=10^4 # res=NULL # for (l in lam) # { # a=tri.prob(l,N) # res=rbind(res,c(l,a$prob,a$cl.n/N,a$sm.n/N,a$b.n/N)) # cat(l,a$prob,a$cl.n/N,a$sm.n/N,a$b.n/N,"\n") # } # colnames(res)=c("lambda","prob","pr.clos.n","pr.minang.n","pr.both") # dput(res,"res.dat") # # As this is finished, you can always recover the results with # res=dget("res.dat") # # Under *NUX, one can put the above commands in a file, say, run.R and run R in a batch mode: # > R CMD BATCH run.R # Or # > cat "R CMD BATCH run.R" > go.sh # and run simulations at a specified time: # > at -f go.sh now # and go to sleep. # # The code also can be run interactively, showing actual configurations and the triangles # tried to triangulate 0. Edges larger than 2 are shown red. So triangulation occurs only # if all edges are green and the triangle contains the origin. See this by running: # tri.prob(1,N=10,details=TRUE,draw=TRUE) # type CTRL-C to interrupt # Other relevant info - polar coordinates of the nodes, hit or miss is also reported now. # # Enjoy! SZ is.Oinside=function(p1,p2,p3,draw=FALSE) { # check whether the origin O is inside the triangle over p1,p2,p3 # not actually used in the below code. if(draw) { x1=min(p1[1],p2[1],p3[1],0) x2=max(p1[1],p2[1],p3[1],0) y1=min(p1[2],p2[2],p3[2],0) y2=max(p1[2],p2[2],p3[2],0) plot(c(x1,x2),c(y1,y2),ty="n",xlab="",ylab="") lines(c(p1[1],p2[1],p3[1],p1[1]),c(p1[2],p2[2],p3[2],p1[2]),col="red") points(0,0,pch=19,col="blue") text(p1[1],p1[2],"p1") text(p2[1],p2[2],"p2") text(p3[1],p3[2],"p3") } a12=p1[1]*p2[2]-p1[2]*p2[1] # +- twice the area of triang O,p1,p2 a23=p2[1]*p3[2]-p2[2]*p3[1] # +- twice the area of triang O,p2,p3 a31=p3[1]*p1[2]-p3[2]*p1[1] # +- twice the area of triang O,p1,p3 # if O is inside, these sum up to twice the area of the triangle: a=abs(a12+a23+a31) return(!(abs(a12)+abs(a23)+abs(a31)>a)) } triangulated=function(r1,r2,r3,phi1,phi2,phi3,R=2) { # check whether the origin O is inside the triangle over 3 points # with polar radii r1--r3 and angles phi1--phi3 AND all the triangle # edges' lengths do not exceed R OK=FALSE a12=r1*r2*sin(phi2-phi1) # +- twice the area of triang O,p1,p2 a23=r2*r2*sin(phi3-phi2) # +- twice the area of triang O,p2,p3 a31=r3*r1*sin(phi1-phi3) # +- twice the area of triang O,p1,p3 # if O is inside, these sum up to twice the area of the triangle: a=abs(a12+a23+a31) if(!(abs(a12)+abs(a23)+abs(a31)>a)) { # now check whether all sides of the triange # have lengths not exceeding R l12=r1^2+r2^2 -2*r1*r2*cos(phi2-phi1) l23=r2^2+r3^2 -2* r2*r3*cos(phi3-phi2) l31=r1^2+r3^2 -2*r1*r3*cos(phi1-phi3) if(all(l12<=R^2,l23<=R^2,l31<=R^2)) OK=TRUE } return(OK) } gen.nodes=function(lam) { # generate Poisson PP in the ball b(0,2) # return radii and polar angles of the nodes # If the closest to O node xi0 is at the distance greater # than sqrt(2) or the number of nodes is less than 3, # no triangulation is possible. Return NULL is this case # Otherwise, the first point returned will be xi0, the second # is xi1 which is the one within the distance 2 from xi0 with # the smallest angle. If there is no such point, xi1=FALSE r=NULL phi=NULL xi1=FALSE r0=sqrt(-log(runif(1))/lam/pi) # dist. to the closest node if(r0=2) { phi0=runif(n,0,2*pi) rr=sqrt(r0^2+(4-r0^2)*runif(n)) # angles of the nodes not father than 2 from the closest # node to O ff=phi0[rr^2+r0^2+2*rr*r0*cos(phi0)<=4] if(length(ff)>0) { # there IS a xi_1-node xi1=TRUE ind.min.phi=which(phi0==min(ff)) phi=c(pi,phi0[ind.min.phi],phi0[-ind.min.phi]) r=c(r0,rr[ind.min.phi],rr[-ind.min.phi]) } else { phi=c(pi,phi0) r=c(r0,rr) } } } return(list(r=r,phi=phi,xi1=xi1)) } tri.prob=function(lam,N=10^3,details=FALSE,draw=FALSE) { # Estimate triangulation probability of the origin O # given density of nodes lam by N simulations # The communication radius R is set to 1 if(N>10 & draw) { cat("draw should be used with small N. Setting N=10.\n") N=10 } hits=0 closest.node.used=0 smallest.angle.node.used=0 both.nodes.used=0 triv.hole=0 for (i0 in 1:N) { nodes=gen.nodes(lam) r=nodes$r phi=nodes$phi xi1=nodes$xi1 n=length(r) if(details) { if(n>0) { cat("n=",n,"\n") } else { cat("No of nodes is less than 3\n") } } hit=FALSE straight.hit=FALSE if(n>=3) { i1=0 i2=0 i3=0 hit=triangulated(r[1],r[2],r[3],phi[1],phi[2],phi[3]) if(hit) { straight.hit=TRUE } while(!hit & i14,"red","green") lines(c(r[i1]*cos(phi[i1]),r[i2]*cos(phi[i2])), c(r[i1]*sin(phi[i1]),r[i2]*sin(phi[i2])),col=cl,lw=2) cl=ifelse(l23>4,"red","green") lines(c(r[i2]*cos(phi[i2]),r[i3]*cos(phi[i3])), c(r[i2]*sin(phi[i2]),r[i3]*sin(phi[i3])),col=cl,lw=2) cl=ifelse(l31>4,"red","green") lines(c(r[i1]*cos(phi[i1]),r[i3]*cos(phi[i3])), c(r[i1]*sin(phi[i1]),r[i3]*sin(phi[i3])),col=cl,lw=2) #trr=c(r[i1],r[i2],r[i3],r[i1]) #trphi=c(phi[i1],phi[i2],phi[i3],phi[i1]) #lines(trr*cos(trphi),trr*sin(trphi),col="red",lw=2) points(0,0,col="blue",pch=19) lines(2*cos(t)-r[1],2*sin(t),col="green") abline(h=0,col="grey") abline(v=0,col="grey") answ=readline("Hit RETURN to continue...") } } } } } } if(hit) { hits=hits+1 if(i1==1 | straight.hit) closest.node.used=closest.node.used+1 if(xi1 & i2==2 | xi1 & i1==2 | straight.hit) smallest.angle.node.used=smallest.angle.node.used+1 if(i1==1 & i2==2 & xi1 | straight.hit) both.nodes.used=both.nodes.used+1 if(all(r>1)) triv.hole=triv.hole+1 } } if(details) { cat("No. of nodes generated:",length(r),"\n") cat("r=\n") print(r) cat("phi=\n") print(phi) return(list(prob=hits/N,sim=N,r=r,phi=phi,hits=hits,cl.n=closest.node.used, sm.n=smallest.angle.node.used,b.n=both.nodes.used,tr.h=triv.hole)) } return(list(prob=hits/N,sim=N,hits=hits,cl.n=closest.node.used, sm.n=smallest.angle.node.used,b.n=both.nodes.used,tr.h=triv.hole)) } draw.triang=function(lambda=1,R=1,a=10,n=NULL,Z=NULL,IDX=NULL,IEX=NULL,col="lightblue",border="blue", col.cov="snow2",border.cov="snow3",draw.cov=FALSE) { # this draw trianglated area shaded in `col' with all edges not exceeding 2R inside # window W=[0,a]x[0,a] # Unless the Poisson points Z are given, generate those # nodes in [-2R,a+2R]x[-2R,a+2R] (only those nodes # contribute to the geometry inside W) # Z must be a matrix with 2 columns containing coordinates of nodes in each row # Returns matrix Z of generated nodes, matrix IDX containing in each row indices of # the nodes forming triangles; and matrix IEX containing in each row indices of # the edges of connected pairs of nodes. # Example of usage: # a=draw.triang(0.5) # Now change the colours: # x11(); draw.triang(Z=a$Z,IDX=a$IDX,IEX=a$IEX,col="grey") # if draw.cov=TRUE, the sensing zones are also plotted in `col.cov': # these are circles of radii R # # No of points to generate: if(is.null(Z)) { IDX=NULL IEX=NULL if(is.null(n)) { n=0 while(n<=2) n=rpois(1,lambda*(a+4*R)*(a+4*R)) } Z=matrix(runif(2*n,-2*R,a+2*R),ncol=2) } else n=nrow(Z) if(is.null(IDX) | is.null(IEX)) { # Constructing indices of the triangles and edges: for (i1 in 1:(n-2)) { for (i2 in (i1+1):(n-1)) { l12=sum((Z[i1,]-Z[i2,])^2) if(l12<=4*R^2) { IEX=rbind(IEX,c(i1,i2)) for (i3 in (i2+1):n) { l23=sum((Z[i2,]-Z[i3,])^2) if(l23<=4*R^2) { IEX=rbind(IEX,c(i2,i3)) l31=sum((Z[i3,]-Z[i1,])^2) if(l31<=4*R^2) { IEX=rbind(IEX,c(i3,i1)) IDX=rbind(IDX,c(i1,i2,i3)) } } } } } } } # Now plotting: require(MASS) eqscplot(c(0,a),c(0,a),ty="n",xlab="",ylab="",main="Triangulated Areas") if(draw.cov) { symbols(Z,circles=rep(R,n),inch=FALSE,add=TRUE,bg=col.cov) symbols(Z,circles=rep(R,n),inch=FALSE,add=TRUE,fg=border.cov) } for (i in 1:nrow(IDX)) { Tr=rbind(Z[IDX[i,1],],Z[IDX[i,2],],Z[IDX[i,3],]) polygon(Tr,col=col) } for (i in 1:nrow(IEX)) segments(Z[IEX[i,1],1],Z[IEX[i,1],2],Z[IEX[i,2],1],Z[IEX[i,2],2],col=border) points(Z,pch=20) segments(0,0,R,0,col="red",lw=3) segments(0,-0.1,0,0.1,col="red",lw=3) segments(R,-0.1,R,0.1,col="red",lw=3) text(R/2,0,labels="R",pos=3,col="red",font=4) invisible(list(Z=Z,IDX=IDX,IEX=IEX)) }