/* Code to generate samples from the autonormal distribution, also * known as a Gaussian field or free field. * Uses the algorithm described in Wilson 1999. * Makes a massless field on a 2-D toroidal grid with wraparound * boundary conditions. */ /* Parameters: $\sigma =$ standard deviation of normals Somebody else does: $X_0 := \Normal(\mu_0,\sigma)$ We generate random variables: $X := (X0 - \mu_0)$ $Y := \exp(-(X/\sigma)^2/2) * \Uniform(0,1)$ If ($X>0$) Then $Y := 1-Y$ $L := -\sigma\sqrt{-2\log(Y)}$ $R := \sigma\sqrt{-2\log(1-Y)}$ The random map is: $f(\mu) = \phi(\mu,(L,R,\mu_0,X)) = \lfloor (\mu-\mu_0+R-X)/(R-L)\rfloor (R-L) +\mu_0+X$ */ #include #include #define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE))) #define CALLOC(NUMBER,TYPE) ((TYPE*)Calloc((NUMBER),sizeof(TYPE))) #define REALLOC(PTR,NUMBER,TYPE) ((TYPE*)Realloc((PTR),(NUMBER)*sizeof(TYPE))) #define FREE(PTR) free(PTR) void *Malloc(), *Calloc(), *Realloc(), Free(); #define pi 3.14159265358979323846264338327950288419716939937510 inline double uniform() { return (((((double)random()) / 2147483648e0) + (double)random()) / 2147483648e0); } /* Use the (basic form of) the Box-Muller (1958) to get Gaussian */ double normal () { return sqrt(-2*log(uniform())) * cos(2*pi*uniform()); } /* Update a lower, middle and upper configuration of the autonormal * by doing a heat bath at each site in sequence. Middle configuration * may be NULL, in which case it does not exist and is not updated. * The three (or two) configurations are updated in a coupled fashion, * using the layered multishift coupler for the normal distribution. */ HeatBathSweep (int width, int height, double force, double *lower, double *middle, double *upper) { int i, n, north, east, west, south; double x,y,l,r,mu; /* Four springs with given force collectively behave like 4*force */ double sigma = sqrt(1/(4*force)); n = width*height; /* site 0 tied to 0 */ lower[0] = upper[0] = 0; if (middle) middle[0] = 0; /* sweep through remaining sites */ for (i=1; i=n) south -= n; east = i+1; if (east%width==0) east -= width; west = i-1; if (i%width==0) west += width; /* For the three (or two) configurations, compute the average of the neighbors, and set site to be normal centered on this average, using the layered multishift coupler. */ mu = (lower[north] + lower[east] + lower[west] + lower[south])/4; lower[i] = floor((mu+r-x)/(r-l))*(r-l) + x; if (middle) { mu = (middle[north] + middle[east] + middle[west] + middle[south])/4; middle[i]= floor((mu+r-x)/(r-l))*(r-l) + x; } mu = (upper[north] + upper[east] + upper[west] + upper[south])/4; upper[i] = floor((mu+r-x)/(r-l))*(r-l) + x; } } /* Generate a proposal state for a Metropolis-Hastings update, * using the tree method for picking the proposal. * Part of the independence sampler. */ MHproposal (int width, int height, double force, double *proposal) { int i, j; /* For the proposal, act as if springs in tree have only half strength */ double sigma = sqrt(1/(force/2)); proposal[0] = 0; for (i=1; i #include #include /* Apply a random composite map to a given state, and output whether * or not the composite map is officially coalescent. * First do a timing experiment. * Then update the given state using the independence sampler, * while also obtaining upper and lower bounds on the updated value * of any state. * Then update the given state and upper and lower bounds using * sweeps for the heatbath, for as long as specified by the timing * experiment. */ #define Diagnose fprintf(stderr,"%d,%lf ",count,Gap(width,height,lower,upper)) ApplyCompositeMap (int width, int height, double force, double *state) { double *proposal, *upper, *lower; double emax; int count; int flag; time_t last, current; proposal = CALLOC(width*height,double); lower = CALLOC(width*height,double); upper = CALLOC(width*height,double); /* Timing experiment */ MHproposal(width,height,force,proposal); emax = Energy(width,height,force,proposal,max); ExtremalStates(width, height, force, lower, upper, emax); last = time(0); for (count=0,Diagnose; memcmp(lower,upper,width*height*sizeof(double)); ++count) { HeatBathSweep(width,height,force, lower, NULL, upper); if ((current=time(0))-last > 1) { Diagnose; fflush(stderr); last=current; } } Diagnose; fprintf(stderr,"\n"); fflush(stderr); /* Actual updates */ MHproposal(width,height,force,proposal); /* Maybe the state gets changed to the proposal. */ if (uniform() < exp(Energy(width,height,force,proposal,tt_m_ht) - Energy(width,height,force,state,tt_m_ht))) memcpy(state,proposal,width*height*sizeof(double)); emax = Energy(width,height,force,proposal,max); ExtremalStates(width,height,force, lower, upper, emax); last = time(0); for (Diagnose; count; --count) { HeatBathSweep(width,height,force, lower, state, upper); if ((current=time(0))-last > 1) { Diagnose; fflush(stderr); last=current; } } Diagnose; fprintf(stderr,"\n"); fflush(stderr); flag = !(memcmp(lower,upper,width*height*sizeof(double))); FREE(upper); FREE(lower); FREE(proposal); return flag; } /* Do Read-Once Coupling from the Past using the composite maps. */ /* Structure for creating a stream of autonormal configurations * with given width, height, and force. */ typedef struct { int width, height; double force; double *state, *oldstate; } AutonormalStream; /* Initialize Read-Once Coupling from the Past. */ AutonormalStream AutonormalInit (int width, int height, double force) { AutonormalStream str; str.width = width; str.height = height; str.force = force; str.oldstate = CALLOC(width*height,double); str.state = CALLOC(width*height,double); /* Keep going until a coalescent map is found. */ while (!ApplyCompositeMap(width, height, force, str.state)); return str; } /* Generate the next sample using Read-Once Coupling from the Past. * If called more than once, then previously returned values are overwritten. */ double *NextAutonormal (AutonormalStream str) { while (1) { memcpy(str.oldstate,str.state,str.width*str.height*sizeof(double)); if (ApplyCompositeMap(str.width, str.height, str.force, str.state)) return str.oldstate; } } main (int argc, char *argv[]) { int width, height, num_iter, seed; double force; AutonormalStream str; if (argc!=6) { printf("%s \n",argv[0]); exit(-1); } width = atoi(argv[1]); height = atoi(argv[2]); force = atof(argv[3]); num_iter = atoi(argv[4]); seed = atoi(argv[5]); if (width<=0 || height<=0 || force<=0 || num_iter<=0) { printf("Width, height, force, and iterations must be positive.\n"); exit(-1); } printf("Generating %d autonormal's on %dX%d grid with coupling constant %lf\nusing seed %d.\n", num_iter,width,height,force,seed); srandom(seed); str = AutonormalInit(width, height, force); for (; num_iter; --num_iter) ShowAutonormal(width,height,NextAutonormal(str)); } ShowAutonormal (int width, int height, double *state) { int i,j; for (j=0; j