/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "ring.p" */ #include /* ring: z space ring Tom Schneider, 1989 */ /* end of program */ /* begin module version */ #define version 3.00 /* of ring.p 1989 Nov 25 origin 1989 July 25 */ /* end module version */ /* begin module describe.ring */ /* name ring: z space ring synopsis ring(data: in, ringp: in, color: output, output: out) files data: set of Gaussianly distributed variables from the program gentst. ringp: parameters: first line: total dimensionality D. second line: number of points to do. If the end of the data is reached, the actual number of points generated is reported to output. third line: number of steps to generate the fD(r) graph. (The range for this is always -2.5 to +2.5 on both x and y axes.) If the number of steps is less than 1, then no smooth graph is done. fourth line: a real number, "0 <= partial <= 1" by which to multiply the actual fD(r) density by to obtain the density reported to the color file. This allows one to tone down the gray scale, or to avoid having the highest density of color equal the lowest (as when the hue is used and a hue of 1 is the same as a hue of 0). fifth line: printing of data on plot (one character): d=dimension,p=dimension+point,a=all,n=none color: a xyin file for input to the xyplo or riden program. The columns are: 1 symbols: f=from fD(r), s = simulated point'); 2 x: x coordinate 3 y: y coordinate 4 xwidth: width of symbol on x axis 5 ywidth: width of symbol on y axis 6 density: density 7 inverse: 1 - density (for inverse plotting) 8 maximum: MAXimum density 9 minimum: MINimum density 10 maximum: MAXimum density 11 minimum: MINimum density 12 partial: partial density for grey tones Partial is the largest density allowed. When plotted in color, hues come from a color wheel in which the highest color is almost identical to the lowest color. That is the color of hue=1 is almost identical to the color of hue = 0. To avoid this effect, make partial less than 1.0. A partial less than 1.0 also avoids completely black gray scale plots. output: messages to the user, number of points generated. description Simulate mapping from many-dimensional to 2-dimensional Z space. Sets of D Gaussian values are read from the data file, squared, summed and square rooted. The x and y value in Z space is determined from an angle and a radius. The angle is found from the last two Gaussian values, while the radius is determined by the noise (rms) for all dimensions. The statistical function fD(r) is to be graphed in color or gray scale using xyplo, while the simulated points are graphed as points on top of the smooth fD(r) function. The program output is ready to read into the xyplo plotting program. examples ringp used for generating figures: 16 total dimensionality 100 number of points to do 128 steps for plotting smooth fD(r) graph 0.50 partial d d=dimension,p=dimension+point,a=all,n=none xyplop used for generating figures: 2 2 zerox zeroy graph coordinate center x -2.5 2.5 zx min max (character, real, real) if zx='x' then set xaxis y -2.5 2.5 zy min max (character, real, real) if zy='y' then set yaxis 10 10 xinterval yinterval number of intervals on axes to plot 4 4 xwidth ywidth width of numbers in characters 1 1 xdecimal ydecimal number of decimal places 5 5 xsize ysize size of axes in inches x y c zc if zc='c' then a crosshairs put on zero of x and y n 2 zxl base if zxl='l' then make x axis log to the given base n 2 zyl base if zyl='l' then make y axis log to the given base ********************************************************************* 2 3 xcolumn ycolumn columns of xyin that determine plot location 1 symbol column the xyin column to read symbols from 4 5 xscolumn yscolumn columns of xyin that determine the symbol size 10 8 7 hue saturation brightness columns for color manipulation ********************************************************************* r symbol to plot c(circle)bd(dotted box)x+Ifgpr(rectangle) b symbol flag character in xyin that indicates that this symbol -1.0 symbol sizex side in inches on the x axis of the symbol. -1.0 symbol sizey as for the x axis, get size from yscolumn n no connection (example for connection is c- 0.05 for dashed 0.05 inch) n 0.05 linetype size linetype l.-in and size of dashes or dots ********************************************************************* r symbol to plot c(circle)bd(dotted box)x+Ifgpr(rectangle) f symbol flag character in xyin that indicates that this symbol -1.0 symbol sizex side in inches on the x axis of the symbol. -1.0 symbol sizey as for the x axis, get size from yscolumn n no connection (example for connection is c- 0.05 for dashed 0.05 inch) n 0.05 linetype size linetype l.-in and size of dashes or dots ********************************************************************* c symbol to plot c(circle)bd(dotted box)x+Ifgpr(rectangle) s symbol flag character in xyin that indicates that this symbol 0.0858 symbol sizex side in inches on the x axis of the symbol. 0.0858 symbol sizey as for the x axis, get size from yscolumn n no connection (example for connection is c- 0.05 for dashed 0.05 inch) n 0.05 linetype size linetype l.-in and size of dashes or dots ********************************************************************* g symbol to plot c(circle)bd(dotted box)x+Ifgpr(rectangle) g symbol flag character in xyin that indicates that this symbol -1.0 symbol sizex side in inches on the x axis of the symbol. -1.0 symbol sizey as for the x axis, get size from yscolumn n no connection (example for connection is c- 0.05 for dashed 0.05 inch) n 0.05 linetype size linetype l.-in and size of dashes or dots ********************************************************************* . ********************************************************************* Useful color parameters are: 8 6 10 Light density plot, printable on a black and white device (best). 8 7 10 Dark density plot, printable on a black and white device. 6 8 10 Color plot, red background. 7 8 10 Color plot, purple background (neat). 6 7 10 Color and density varying to make the simulated points easy to see. (red background) 7 6 10 Color and density varying to make the simulated points easy to see. (white background - lovely!) Warning: since the program has changed, these may no longer be correct. documentation ccmm see also gentst.p xyplo.p riden.p author Thomas Dana Schneider bugs none known. Confirm that the density distribution is correct by using program riden. technical notes */ /* end module describe.ring */ /* begin module sphere.const */ #define wid 6 /* width of numbers */ #define dec 4 /* number of decimal places of numbers */ #define Rmaximum 1.0 /* where the curves are to have their maximum */ /* end module sphere.const */ /* begin module ring.const */ #define overlap 1.10 /* percentage of overlap between squares, 1.01 means 1%. to avoid gaps. See jump and boxwidth. */ #define rangemax 2.5 /* range of the fD(r) plot */ /* end module ring.const */ /* begin module ring.var */ Static _TEXT data; /* Gaussian data */ Static _TEXT ringp; /* parameters */ Static _TEXT color; /* density output */ Static double pi; /* circumfrence divided by diameter of circle */ Static jmp_buf _JL1; /* end module ring.var */ /* begin module halt */ Static Void halt() { /* stop the program. the procedure performs a goto to the end of the program. you must have a label: label 1; declared, and also the end of the program must have this label: 1: end. examples are in the module libraries. this is the only goto in the delila system. */ printf(" program halt.\n"); longjmp(_JL1, 1); } /* end module halt version = 1.33; (@ of domod, 1989 July 25 */ /* begin module copyaline */ Static Void copyaline(fin, fout) _TEXT *fin, *fout; { /* copy a line from file fin to file fout */ while (!P_eoln(fin->f)) { putc(P_peek(fin->f), fout->f); getc(fin->f); } fscanf(fin->f, "%*[^\n]"); getc(fin->f); putc('\n', fout->f); } /* copyaline */ /* end module copyaline version = 1.33; (@ of domod, 1989 July 25 */ /* begin module ring.readparameters */ Static Void readparameters(f, D, n, steps, partial, idstring) _TEXT *f; long *D, *n, *steps; double *partial; Char *idstring; { /* read parameters */ if (*f->name != '\0') { if (f->f != NULL) f->f = freopen(f->name, "r", f->f); else f->f = fopen(f->name, "r"); } else rewind(f->f); if (f->f == NULL) _EscIO2(FileNotFound, f->name); RESETBUF(f->f, Char); if (BUFEOF(f->f)) { printf("missing parameter: D\n"); halt(); } fscanf(f->f, "%ld%*[^\n]", D); getc(f->f); if (BUFEOF(f->f)) { printf("missing parameter: n\n"); halt(); } fscanf(f->f, "%ld%*[^\n]", n); getc(f->f); if (BUFEOF(f->f)) { printf("missing parameter: steps\n"); halt(); } fscanf(f->f, "%ld%*[^\n]", steps); getc(f->f); if (BUFEOF(f->f)) { printf("missing parameter: partial\n"); halt(); } fscanf(f->f, "%lg%*[^\n]", partial); getc(f->f); if (BUFEOF(f->f)) { printf("missing parameter: idstring\n"); halt(); } fscanf(f->f, "%c%*[^\n]", idstring); getc(f->f); if (*idstring == '\n') *idstring = ' '; if (*idstring != 'p' && *idstring != 'd' && *idstring != 'a' && *idstring != 'n') *idstring = 'n'; } /* end module ring.readparameters */ /* begin module ring.writeparameters */ Static Void writeparameters(f, D, n, steps, partial, idstring) _TEXT *f; long D, n, steps; double partial; Char idstring; { /* write parameters */ fprintf(f->f, "* \n"); fprintf(f->f, "* Parameters:\n"); fprintf(f->f, "* %10ld total dimensionality, D\n", D); fprintf(f->f, "* %10ld requested number of points to plot\n", n); fprintf(f->f, "* %10ld steps for plotting smooth fD(r) graph\n", steps); fprintf(f->f, "* %10.2f partial value used for plotting\n", partial); fprintf(f->f, "* %10c parameter data on plot: ", idstring); if (idstring == 'd') fprintf(f->f, "dimension\n"); if (idstring == 'p') fprintf(f->f, "dimension + point\n"); if (idstring == 'a') fprintf(f->f, "all\n"); if (idstring == 'n') fprintf(f->f, "none\n"); fprintf(f->f, "* \n"); } /* end module ring.writeparameters */ /* begin module pic.polrec */ Static Void polrec(r, theta, x, y) double r, theta, *x, *y; { /* convert polar to rectangular coordinates, theta is in radians */ *x = r * cos(theta); *y = r * sin(theta); } /* end module pic.polrec version = 1.33; (@ of domod, 1989 July 25 */ /* begin module ring.getpoint */ Static boolean getpoint(data, D, plotx, ploty) _TEXT *data; long D; double *plotx, *ploty; { /* read D numbers from data file, create the root mean square of the Gaussians to form the radus r. The radius is normalized by dividing by sqrt(D-1). Find the angle from the last two Gaussians. If there was not enough data in the file, gotpoint is false, otherwise true. */ long dindex = 0; /* index for number of D */ double partialsqd = 0.0; /* sum of square of Gaussians, MISSING 2 deg. free. */ double Gaussian; /* a Gaussianly distributed variable */ boolean gotpoint; /* true if the point was gotten */ double rz; /* radius to the point in Z space */ double x, y; /* the last two coordinates */ double Znormalize; /* the radius in Z space is normalized with this */ double Zrsqd; /* the square of the radius to the point (x,y) in Z space */ /* If it is not 2 dimensional, pick up the first D-2 coordinates, and calculate the square distance */ /* p2c: ring.p: Note: Eliminated unused assignment statement [338] */ if (D > 2) { while (!BUFEOF(data->f) && dindex < D - 2) { fscanf(data->f, "%lg%*[^\n]", &Gaussian); getc(data->f); partialsqd += Gaussian * Gaussian; dindex++; } } /* pick up the last D coordinates */ if (dindex < D - 2 && D > 2) gotpoint = false; else { if (BUFEOF(data->f)) gotpoint = false; else { fscanf(data->f, "%lg%*[^\n]", &x); getc(data->f); if (BUFEOF(data->f)) gotpoint = false; else { fscanf(data->f, "%lg%*[^\n]", &y); getc(data->f); gotpoint = true; } } } if (gotpoint) { Zrsqd = x * x + y * y; rz = sqrt((partialsqd + Zrsqd) / (D - 1)); Znormalize = rz / sqrt(Zrsqd); *plotx = x * Znormalize; *ploty = y * Znormalize; return true; } else { *plotx = 0.0; *ploty = 0.0; return false; } } /* end module ring.getpoint */ /* begin module ring.functions */ /* the following two functions were taken from sphere.p */ Static double pd(d, r, twosigma2) long d; double r, twosigma2; { /* probability at a given radius in a given dimension */ return (exp(-(r * r / twosigma2)) * exp((d - 1) * log(r))); } Static double lnpd(d, r, twosigma2) long d; double r, twosigma2; { /* natural log f the pd function, to make calculation of high dimensions possible */ /* probability at a given radius in a given dimension */ return ((d - 1) * log(r) - r * r / twosigma2); } /* Local variables for themain: */ struct LOC_themain { _TEXT *color; double boxwidth, halfbox; /* half of boxwidth, for positioning the squares correctly */ double partial; /* a column of the output data for partial tones */ /* 2 times sigma squared, where sigma is the standard deviation of the gaussian */ double density; /* density of the fD(r) function at radius r */ } ; Local Void dofdr(xplot, yplot, LINK) double xplot, yplot; struct LOC_themain *LINK; { /* for fD(r) data */ /* do the fdr box */ /* max */ /* min */ /* max */ /* min */ /* partial tones */ fprintf(LINK->color->f, "f %*.*f %*.*f %*.*f %*.*f %*.*f %*.*f 1 0 1 0 %4.2f\n", wid, dec, xplot - LINK->halfbox, wid, dec, yplot - LINK->halfbox, wid, dec, LINK->boxwidth, wid, dec, LINK->boxwidth, wid, dec, LINK->density * LINK->partial, wid, dec, 1.0 - LINK->density * LINK->partial, LINK->partial); } /* end module ring.functions */ /* begin module ring.themain */ Static Void themain(data, ringp, color_) _TEXT *data, *ringp, *color_; { /* the main procedure of the program */ struct LOC_themain V; /* the size to plot, overlap % bigger than jump, to avoid aliasing effects, and ensure that the squares overlap */ Char idstring; /* d = dimension only, p = dimension + points a = dimension + points + other parameters (all) n = none */ long D; /* total dimensionality */ double halfjump; /* half of jump */ double jump; /* the interval to increment in the fD(r) plot */ long n; /* number of points to do */ long nindex = 0; /* index for n */ long steps; /* number of steps for fD(r) color graph */ double xplot, yplot; /* coordinate to plot */ /* from sphere.p */ double lnpdmax; /* ln of maximum value of pd curve, for normalizing */ /* pdmax: real; (@ maximum value of pd curve, for normalizing */ double twosigma2, lndensity; /* ln of density */ double r; /* radius on the plot */ _TEXT TEMP; V.color = color_; printf("ring %4.2f\n", version); pi = 4.0 * atan(1.0); /*writeln(output,'pi = ',pi:20:18);*/ if (*data->name != '\0') { if (data->f != NULL) data->f = freopen(data->name, "r", data->f); else data->f = fopen(data->name, "r"); } else rewind(data->f); if (data->f == NULL) _EscIO2(FileNotFound, data->name); RESETBUF(data->f, Char); if (*V.color->name != '\0') { if (V.color->f != NULL) V.color->f = freopen(V.color->name, "w", V.color->f); else V.color->f = fopen(V.color->name, "w"); } else { if (V.color->f != NULL) rewind(V.color->f); else V.color->f = tmpfile(); } if (V.color->f == NULL) _EscIO2(FileNotFound, V.color->name); SETUPBUF(V.color->f, Char); fprintf(V.color->f, "* ring %4.2f\n", version); fprintf(V.color->f, "*\n"); /* read in the data */ fprintf(V.color->f, "* Gaussians from:\n"); while (P_peek(data->f) == '*') copyaline(data, V.color); readparameters(ringp, &D, &n, &steps, &V.partial, &idstring); TEMP.f = stdout; *TEMP.name = '\0'; writeparameters(&TEMP, D, n, steps, V.partial, idstring); if (D < 2) { printf("dimensionality must be greater than 1\n"); halt(); } writeparameters(V.color, D, n, steps, V.partial, idstring); /* Calculation of fD(r) maximum (from sphere.p) */ if (D > 1) { twosigma2 = 2.0 / (D - 1); /* pdmax := pd(D,1.0, twosigma2); */ lnpdmax = lnpd(D, 1.0, twosigma2); } else { twosigma2 = 2.0; /* pdmax := 1.0; */ lnpdmax = log(1.0); } /* make gaussian */ fprintf(V.color->f, "* xyin columns are:\n"); fprintf(V.color->f, "* symbols: b=background rectangle, f=from fD(r), s = simulated point\n"); fprintf(V.color->f, "* symbols| x| y| xwidth| ywidth| density| 1-density| max| min| max| min| partial\n"); fprintf(V.color->f, "* 1| 2| 3| 4| 5| 6| 7| 8| 9| 10| 11| 12\n"); fprintf(V.color->f, "*\n"); /* for background */ /* create background of the plot points */ /* x width of rectangle */ /* y width of rectangle */ /* density */ /* 1 - density */ /* max */ /* min */ /* max */ /* min */ /* partial tones */ fprintf(V.color->f, "b %*.*f %*.*f %*.*f %*.*f %*.*f %*.*f 1 0 1 0 %4.2f\n", wid, dec, -rangemax, wid, dec, -rangemax, wid, dec, 2 * rangemax, wid, dec, 2 * rangemax, wid, dec, 0.0, wid, dec, 1.0, V.partial); /* figure out box sizes and jumps */ /* We wish to pack a number of squares (steps of them) into the region of width 0 to rangemax, leaving a small amount of space at the edge. Instead of printing the squares, we will print a square a few percent bigger, so that they overlap without gaps. (PostScript leaves gaps between some of the rectangles, probably due to roundoff errors.) The variable 'overlap' defines how much the box width should be expanded during printing so that they overlap. It is the fraction of expansion plus 1.0. For example, a 10% expansion is defined by using overlap = 1.10. Half of this overlap should lie around the edge of the plot, but within the bounds of rangemax. So we have steps * rawboxwidth + (overlap - 1)*rawboxwidth/2 = rangemax This is rearranged to find the "raw box width". This raw number is also the amount to "jump" between printings of boxes. */ jump = rangemax / (steps + (overlap - 1) / 2); V.boxwidth = overlap * jump; V.halfbox = V.boxwidth / 2.0; halfjump = jump / 2.0; /* create fD(r) points */ xplot = halfjump; yplot = halfjump; /* only do the smooth graph if the number of steps requested is positive */ if (steps > 0) { /* for background */ while (xplot < rangemax) { r = sqrt(xplot * xplot + yplot * yplot); /* pd(D,r, twosigma2)/pdmax */ if (r != 0.0) lndensity = lnpd(D, r, twosigma2) - lnpdmax; else lndensity = -25.0; /* By testing for the log of this density, we can avoid exponentiation of very small density. The result is not affected, for these density come up only when the final result is nearly zero anyway. Don't plot if density is zero. */ if (lndensity >= -25) { V.density = exp(lndensity); /* plot as 8 quadrants to save calculation time */ dofdr(xplot, yplot, &V); dofdr(yplot, xplot, &V); dofdr(yplot, -xplot, &V); dofdr(-xplot, yplot, &V); dofdr(-xplot, -yplot, &V); dofdr(-yplot, -xplot, &V); dofdr(-yplot, xplot, &V); dofdr(xplot, -yplot, &V); } yplot += jump; if (yplot > xplot) { yplot = halfjump; xplot += jump; } } } /* make a tiny black box to reset the colors */ /* max */ /* min */ /* max */ /* min */ /* partial tones */ fprintf(V.color->f, "b %*.*f %*.*f %*.*f %*.*f %*.*f %*.*f 0 1 0 1 %4.2f\n", wid, dec, 0.0, wid, dec, 0.0, wid, dec, 0.0, wid, dec, 0.0, wid, dec, 0.0, wid, dec, 0.0, V.partial); /* create simulated points */ while ((nindex < n) & (!BUFEOF(data->f))) { if (getpoint(data, D, &xplot, &yplot)) { nindex++; /* for simulated data */ /* density */ /* 1 - density */ /* full */ /* minimal */ /* full */ /* minimal */ /* partial tones */ fprintf(V.color->f, "s %*.*f %*.*f %*.*f %*.*f 0 1 1 0 1 0 %4.2f\n", wid, dec, xplot, wid, dec, yplot, wid, dec, V.boxwidth, wid, dec, V.boxwidth, V.partial); } } /* create descriptive string */ if (idstring == 'a' || idstring == 'p' || idstring == 'd') { fprintf(V.color->f, "D=%ld", D); if (idstring == 'a' || idstring == 'p') { /* x location */ fprintf(V.color->f, "__%ld_points", nindex); if (idstring == 'a') fprintf(V.color->f, "__%ld_steps__%4.2f_partial", steps, V.partial); } /* y location */ /* width */ /* width */ /* density */ /* 1 - density */ /* full */ /* minimal */ /* full */ /* minimal */ /* partial tones */ fprintf(V.color->f, " %*.*f %*.*f %d %d 0 1 1 0 1 0 %4.2f\n", wid, dec, -1.3, wid, dec, 2.5, 1, 1, V.partial); } if (idstring == 'a') /* x location */ fprintf(V.color->f, "ring_%4.2f %*.*f %*.*f %d %d 0 1 1 0 1 0 %4.2f\n", version, wid, dec, 2.0, wid, dec, 2.5, 1, 1, V.partial); /* y location */ /* width */ /* width */ /* density */ /* 1 - density */ /* full */ /* minimal */ /* full */ /* minimal */ /* partial tones */ if (nindex == n) printf(" %ld points done\n", nindex); else printf(" Did not reach %ld data points, Only %ld points done\n", n, nindex); } /* end module ring.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; color.f = NULL; strcpy(color.name, "color"); ringp.f = NULL; strcpy(ringp.name, "ringp"); data.f = NULL; strcpy(data.name, "data"); themain(&data, &ringp, &color); _L1: if (data.f != NULL) fclose(data.f); if (ringp.f != NULL) fclose(ringp.f); if (color.f != NULL) fclose(color.f); exit(EXIT_SUCCESS); } /* End. */