#include #include #include #include // k-nearest neighbor Monte Carlo high confidence bounds --------------------- #define K 4 // number of nearest neighbors #define X 'B' // type of component U/O/B #define BND 'L' // upper/lower bound U/L #define R 50. // radius of central disk #define S 0. // buffer zone around central disk #define ITER 100 // number of iterations #define M 100000 // maximum number of points in rectangle #define SEED 1 // random number generator seed // For simplicity, mostly use global variables ------------------------------- double vx[M],vy[M]; // coordinates of points int st[M]; // status 0-normal 1-left 2-right 3-done int pc[M]; // vertex is in cluster or comes from boundary int qq[M],qs,qe; // queue of vertices int nb[M][6*K+1]; // neighbor list int nbb[M][K]; // neighbor backup list double ds[K+2]; // list of distances^2 from given vertex int nt[K+2]; // temp neighbor list int nv; // total number of vertices int cl,cr; // left/right vertex count int sl,sr; // left/right success count int erc; // error count (too many vertices in rectangle) // Random number generator --------------------------------------------------- unsigned char rs[256],ri,rj; // random number generator status // Random number, uniform in [0,1] (never =0) double rand_uniform(){ unsigned int s,t,i; double w; for(w=.5,i=0;i<7;i++){ // generates 56 bits of random data ri=(ri+1)&0xFF;rj=(rj+rs[ri])&0xFF; // &0xFF to fix compiler bug s=rs[ri];rs[ri]=t=rs[rj];rs[rj]=s; w=(rs[(s+t)&0xFF]+w)/256; // rs[(s+t)&0xFF]=random byte } return w; // w=0.{56 random bits}1 in binary } // Initialize random number generator // Call before first use of rand_uniform() void rand_setup(){ unsigned int s,t; int i; ri=1; rj=SEED; // initialize with seed for(i=0;i<256;i++)rs[i]=i; // setup permutation array for(i=0;i<64;i++)rand_uniform(); // run algorithm to scramble rs[] } // Generate instance --------------------------------------------------------- // generate graph and modify edges according to X = U/O/I/B // int generate(){ double x,y,d; int i,j,k,t,w; // first generate the points nv=0; x=-2*(R+S); cl=cr=0; for(;x<2*(R+S) && nvy)x=y; d=x*x; for(j=0;j=0;j--){ x=vx[i]-vx[j]; if(x*x>=ds[K-1])break; y=vy[i]-vy[j]; d=x*x+y*y; if(d0 && ds[k-1]>d;k--){ nb[i][k]=nb[i][k-1]; ds[k]=ds[k-1]; } nb[i][k]=j; ds[k]=d; } } for(j=i+1;j=ds[K-1])break; y=vy[i]-vy[j]; d=x*x+y*y; if(d0 && ds[k-1]>d;k--){ nb[i][k]=nb[i][k-1]; ds[k]=ds[k-1]; } nb[i][k]=j; ds[k]=d; } } } if(X=='U'||X=='I'){ // undirected / in-component case for(i=0;i=0;k++){ j=nbb[i][k]; for(t=0;nb[j][t]>=0 && nb[j][t]!=i;t++); if(nb[j][t]<0){ if(t>=6*K){printf("\nFatal error - too many neighbors\n");exit(1);} nb[j][t]=i; nb[j][t+1]=-1; } } } } if(X=='B'){ // bidirectional case for(i=0;i=0;k++){ j=nb[i][k]; for(t=0;nb[j][t]>=0 && nb[j][t]!=i;t++); if(nb[j][t]==i){nb[i][w]=j;w++;} } nb[i][w]=-1; } } // check for incoming edges from boundary // (only lower bound U/O/B cases) if(BND!='U' && X!='I'){ if(nv0 && ds[k-1]>d;k--){ nt[k]=nt[k-1]; ds[k]=ds[k-1]; } nt[k]=i; ds[k]=d; } } // we mark the nearest K+1 points // and move (d_{k+2}-d_K)/2 around boundary // this guarantees we get K nearest points of any // boundary point. for(k=0;k=2*(R+S)){y+=d+x-2*(R+S);x=2*(R+S);d=0;} // right if(y>=(R+S)){x-=d+y-(R+S);y=R+S;d=0;} // top if(x<=-2*(R+S)){y-=d-x-2*(R+S);x=-2*(R+S);d=0;} // left if(y<-(R+S))break; } if(X=='B'){ for(i=0;i most points // in central disk on right. // int ubtrial(){ int i,j,k,t,w,b; // generate graph if(!generate())return 0; // Now percolate from each vertex for(b=0;b<2;b++){ // both directions for(i=0;iqs){j=qq[qs++]; for(k=0;nb[j][k]>=0;k++){ if(pc[nb[j][k]]=nv){printf("Fatal queue error\n"); exit(1);} qq[qe++]=nb[j][k]; pc[nb[j][k]]=w; } } if(st[j]==2)sr++; } if(X=='U' || X=='B'){ t=0; for(j=0;jcr)sl+=t; } if(2*sl<=cl)return 0; for(i=0;iR if(st[i])st[i]--; } i=cr;cr=cl;cl=i; } return 1; } // Lower bound trial --------------------------------------------------------- // Test whether there exists path from boundary cutting // line between centers of squares // int lbtrial(){ int i,j,k; // generate graph if(!generate())return 0; // percolate from boundary qs=0;qe=0; for(i=0;iqs){i=qq[qs++]; for(k=0;nb[i][k]>=0;k++){ j=nb[i][k]; if(vy[i]*vy[j]<0){ // test if crosses center line if(fabs((vx[i]*vy[j]-vx[j]*vy[i])/(vy[j]-vy[i]))<=R+S)return 0; } if(!pc[j]){ if(qe>=nv){printf("Fatal queue error\n"); exit(1);} qq[qe++]=j; pc[j]=1; } } } return 1; } // Compute statistics -------------------------------------------------------- main(){ time_t start,end; double f,s; int c,t,i; time(&start); rand_setup(); erc=c=0; for(t=0;t1){f/=1000;s/=1000;t+=3;} } for(;i<=ITER;i++){ s*=0.8639; if(s<1){s*=1000;t-=3;} } while(s>=10-5e-7){s/=10;t++;} while(s<1 && s>0){s*=10;t--;} printf("Confidence=%fE%d\n",s,t); if(erc)printf("%d errors due to too many points\n",erc); printf("%d seconds\n",(int)(end-start)); }