program ring(data, ringp, color, output); (* ring: z space ring Tom Schneider, 1989 *) label 1; (* end of program *) const (* begin module version *) 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 *) wid = 6; (* width of numbers *) dec = 4; (* number of decimal places of numbers *) Rmaximum = 1.0; (* where the curves are to have their maximum *) (* end module sphere.const *) (* begin module ring.const *) overlap = 1.10; (* percentage of overlap between squares, 1.01 means 1%. to avoid gaps. See jump and boxwidth. *) rangemax = 2.5; (* range of the fD(r) plot *) (* end module ring.const *) var (* begin module ring.var *) data: text; (* Gaussian data *) ringp: text; (* parameters *) color: text; (* density output *) pi: real; (* circumfrence divided by diameter of circle *) (* end module ring.var *) (* begin module halt *) procedure 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. *) begin writeln(output,' program halt.'); goto 1 end; (* end module halt version = 1.33; (@ of domod, 1989 July 25 *) (* begin module copyaline *) procedure copyaline(var fin, fout: text); (* copy a line from file fin to file fout *) begin (* copyaline *) while not eoln(fin) do begin fout^ := fin^; put(fout); get(fin) end; readln(fin); writeln(fout); end; (* copyaline *) (* end module copyaline version = 1.33; (@ of domod, 1989 July 25 *) (* begin module ring.readparameters *) procedure readparameters(var f: text; var D: integer; var n: integer; var steps: integer; var partial: real; var idstring: char); (* read parameters *) begin reset(f); if eof(f) then begin writeln(output,'missing parameter: D'); halt end; readln(f,D); if eof(f) then begin writeln(output,'missing parameter: n'); halt end; readln(f,n); if eof(f) then begin writeln(output,'missing parameter: steps'); halt end; readln(f,steps); if eof(f) then begin writeln(output,'missing parameter: partial'); halt end; readln(f,partial); if eof(f) then begin writeln(output,'missing parameter: idstring'); halt end; readln(f,idstring); if not (idstring in ['n','a','d','p']) then idstring := 'n'; end; (* end module ring.readparameters *) (* begin module ring.writeparameters *) procedure writeparameters(var f: text; D: integer; n: integer; steps: integer; partial: real; idstring: char); (* write parameters *) begin writeln(f,'* '); writeln(f,'* Parameters:'); writeln(f,'* ',D:10,' total dimensionality, D'); writeln(f,'* ',n:10,' requested number of points to plot'); writeln(f,'* ',steps:10,' steps for plotting smooth fD(r) graph'); writeln(f,'* ',partial:10:2,' partial value used for plotting'); write (f,'* ',idstring:10,' parameter data on plot: '); if idstring = 'd' then writeln(f,'dimension'); if idstring = 'p' then writeln(f,'dimension + point'); if idstring = 'a' then writeln(f,'all'); if idstring = 'n' then writeln(f,'none'); writeln(f,'* '); end; (* end module ring.writeparameters *) (* begin module pic.polrec *) procedure polrec(r,theta: real; var x,y: real); (* convert polar to rectangular coordinates, theta is in radians *) begin x := r*cos(theta); y := r*sin(theta) end; (* end module pic.polrec version = 1.33; (@ of domod, 1989 July 25 *) (* begin module ring.getpoint *) function getpoint(var data: text; D: integer; var plotx,ploty: real ): boolean; (* 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. *) var dindex: integer; (* index for number of D *) partialsqd: real; (* sum of square of Gaussians, MISSING 2 deg. free. *) Gaussian: real; (* a Gaussianly distributed variable *) gotpoint: boolean; (* true if the point was gotten *) rz: real; (* radius to the point in Z space *) x, y: real; (* the last two coordinates *) Znormalize: real; (* the radius in Z space is normalized with this *) Zrsqd: real; (* the square of the radius to the point (x,y) in Z space *) begin partialsqd := 0.0; dindex := 0; gotpoint := false; (* If it is not 2 dimensional, pick up the first D-2 coordinates, and calculate the square distance *) if D > 2 then while (not eof(data)) and (dindex < D - 2) do begin readln(data,Gaussian); partialsqd := partialsqd + Gaussian * Gaussian; dindex := succ(dindex) end; (* pick up the last D coordinates *) if (dindex < D - 2) and (D > 2) then gotpoint := false else begin if eof(data) then gotpoint := false else begin readln(data, x); if eof(data) then gotpoint := false else begin readln(data, y); gotpoint := true end end end; if gotpoint then begin Zrsqd := x*x + y*y; rz := sqrt((partialsqd + Zrsqd) / (D-1)); Znormalize :=rz/sqrt(Zrsqd); plotx := x * Znormalize; ploty := y * Znormalize; getpoint := true; end else begin plotx := 0.0; ploty := 0.0; getpoint := false end end; (* end module ring.getpoint *) (* begin module ring.functions *) (* the following two functions were taken from sphere.p *) function pd(d: integer; r: real; twosigma2: real): real; (* probability at a given radius in a given dimension *) begin pd := exp(-r*r/twosigma2) * exp((d-1)*ln(r)) end; function lnpd(d: integer; r: real; twosigma2: real): real; (* natural log f the pd function, to make calculation of high dimensions possible *) (* probability at a given radius in a given dimension *) begin lnpd := (-sqr(r)/twosigma2) + ((d-1)*ln(r)) end; (* end module ring.functions *) (* begin module ring.themain *) procedure themain(var data, ringp, color: text); (* the main procedure of the program *) var boxwidth: real; (* the size to plot, overlap % bigger than jump, to avoid aliasing effects, and ensure that the squares overlap *) idstring: char; (* d = dimension only, p = dimension + points a = dimension + points + other parameters (all) n = none *) D: integer; (* total dimensionality *) halfbox: real; (* half of boxwidth, for positioning the squares correctly *) halfjump: real; (* half of jump *) jump: real; (* the interval to increment in the fD(r) plot *) n: integer; (* number of points to do *) nindex: integer; (* index for n *) partial: real; (* a column of the output data for partial tones *) steps: integer; (* number of steps for fD(r) color graph *) xplot, yplot: real; (* coordinate to plot *) (* from sphere.p *) lnpdmax: real; (* ln of maximum value of pd curve, for normalizing *) (* pdmax: real; (@ maximum value of pd curve, for normalizing *) twosigma2: real; (* 2 times sigma squared, where sigma is the standard deviation of the gaussian *) density: real; (* density of the fD(r) function at radius r *) lndensity: real; (* ln of density *) r: real; (* radius on the plot *) procedure dofdr(xplot, yplot: real); (* do the fdr box *) begin writeln(color,'f', (* for fD(r) data *) ' ', (xplot-halfbox):wid:dec, ' ', (yplot-halfbox):wid:dec, ' ', boxwidth:wid:dec, ' ', boxwidth:wid:dec, ' ',(density*partial):wid:dec, ' ',(1.0 - (density*partial)):wid:dec, ' 1', (* max *) ' 0', (* min *) ' 1', (* max *) ' 0', (* min *) ' ', partial:4:2 (* partial tones *) ); end; begin writeln(output,'ring ',version:4:2); pi := 4.0 * arctan(1.0); (*writeln(output,'pi = ',pi:20:18);*) reset(data); rewrite(color); writeln(color,'* ring ',version:4:2); writeln(color,'*'); (* read in the data *) writeln(color,'* Gaussians from:'); while data^='*' do copyaline(data,color); readparameters(ringp,D,n,steps,partial,idstring); writeparameters(output,D,n,steps,partial,idstring); if D < 2 then begin writeln(output,'dimensionality must be greater than 1'); halt end; writeparameters(color,D,n,steps,partial,idstring); (* Calculation of fD(r) maximum (from sphere.p) *) if D > 1 then begin twosigma2 := 2*sqr(Rmaximum)/(D-1); (* pdmax := pd(D,1.0, twosigma2); *) lnpdmax := lnpd(D,1.0, twosigma2); end else begin (* make gaussian *) twosigma2 := 2*sqr(Rmaximum); (* pdmax := 1.0; *) lnpdmax := ln(1.0); end; writeln(color,'* xyin columns are:'); writeln(color,'* symbols: b=background rectangle,', ' f=from fD(r), s = simulated point'); writeln(color,'* symbols| x| y|', ' xwidth| ywidth|', ' density| 1-density|', ' max| min| max| min|', ' partial'); writeln(color,'* 1| 2| 3|', ' 4| 5|', ' 6| 7|', ' 8| 9| 10| 11|', ' 12'); writeln(color,'*'); (* create background of the plot points *) writeln(color,'b', (* for background *) ' ', -rangemax:wid:dec, ' ', -rangemax:wid:dec, ' ',2*rangemax:wid:dec, (* x width of rectangle *) ' ',2*rangemax:wid:dec, (* y width of rectangle *) ' ',0.0:wid:dec, (* density *) ' ',1.0:wid:dec, (* 1 - density *) ' 1', (* max *) ' 0', (* min *) ' 1', (* max *) ' 0', (* min *) ' ', partial:4:2 (* partial tones *) ); (* 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)); boxwidth := overlap * jump; halfbox := 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 then while xplot < rangemax do begin r := sqrt(xplot*xplot + yplot*yplot); (* pd(D,r, twosigma2)/pdmax *); if r <> 0.0 then lndensity := (lnpd(D,r, twosigma2)-lnpdmax) else lndensity := -25; (* 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 then begin density := exp(lndensity); (* plot as 8 quadrants to save calculation time *) dofdr( xplot, yplot); dofdr( yplot, xplot); dofdr( yplot, -xplot); dofdr(-xplot, yplot); dofdr(-xplot, -yplot); dofdr(-yplot, -xplot); dofdr(-yplot, xplot); dofdr( xplot, -yplot); end; yplot := yplot + jump; if yplot > xplot then begin yplot := halfjump; xplot := xplot + jump end end; (* make a tiny black box to reset the colors *) writeln(color,'b', (* for background *) ' ', 0.0:wid:dec, ' ', 0.0:wid:dec, ' ', 0.0:wid:dec, ' ', 0.0:wid:dec, ' ', 0.0:wid:dec, ' ', 0.0:wid:dec, ' 0', (* max *) ' 1', (* min *) ' 0', (* max *) ' 1', (* min *) ' ', partial:4:2 (* partial tones *) ); (* create simulated points *) nindex := 0; while (nindex < n) and (not eof(data)) do begin if getpoint(data, D, xplot, yplot) then begin nindex := succ(nindex); writeln(color,'s', (* for simulated data *) ' ', xplot:wid:dec, ' ', yplot:wid:dec, ' ', boxwidth:wid:dec, ' ', boxwidth:wid:dec, ' 0', (* density *) ' 1', (* 1 - density *) ' 1', (* full *) ' 0', (* minimal *) ' 1', (* full *) ' 0', (* minimal *) ' ', partial:4:2 (* partial tones *) ); end; end; (* create descriptive string *) if idstring in ['d','p','a'] then begin write (color,'D=',D:1); if idstring in ['p','a'] then begin write (color, '__', nindex:1,'_points'); if idstring in ['a'] then begin write (color, '__', steps:1,'_steps', '__', partial:4:2,'_partial'); end; end; writeln(color, ' ', -1.3:wid:dec, (* x location *) ' ', +2.5:wid:dec, (* y location *) ' ', 1:1, (* width *) ' ', 1:1, (* width *) ' 0', (* density *) ' 1', (* 1 - density *) ' 1', (* full *) ' 0', (* minimal *) ' 1', (* full *) ' 0', (* minimal *) ' ', partial:4:2 (* partial tones *) ); end; if idstring in ['a'] then writeln(color,'ring_',version:4:2, ' ', +2.0:wid:dec, (* x location *) ' ', +2.5:wid:dec, (* y location *) ' ', 1:1, (* width *) ' ', 1:1, (* width *) ' 0', (* density *) ' 1', (* 1 - density *) ' 1', (* full *) ' 0', (* minimal *) ' 1', (* full *) ' 0', (* minimal *) ' ', partial:4:2 (* partial tones *) ); if nindex = n then writeln(output,' ', nindex:1, ' points done') else writeln(output,' Did not reach ',n:1,' data points,', ' Only ',nindex:1, ' points done'); end; (* end module ring.themain *) begin themain(data,ringp,color); 1: end.