/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "sphere.p" */ #include /* sphere: plot density of shannon spheres Dr. Thomas D. Schneider National Cancer Institute Laboratory of Experimental and Computational Biology Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu http://www.lecb.ncifcrf.gov/~toms/ */ /* end of program */ /* begin module version */ #define version 1.47 /* of sphere.p 2001 July 27 2001 Jul 25, 1.43: document spherep.fractal 2001 Jul 25, 1.43: noninteger D allowed 2001 Jul 25, 1.42: upgrade documentation. 2001 Jul 25, 1.41: merge documentation with published ccmm version 1998 Jun 1, 1.40: add sphere.xyplop pointer 1994 Sep 6, 1.39: last previous change (used for ccmm paper) 1988 Jan 24, 1.00: origin */ /* end module version */ /* begin module describe.sphere */ /* name sphere: plot density of shannon spheres synopsis sphere(spherep: in, sigma: out, xyin: out, output:out) files spherep: parameters. The first line is the step size interval (0.01 works well). the second line is the maximum radius to calculate out to (= maxr, 3.1 works well). Each following line is a dimension to plot. If the dimension number is negative, it must be followed on the same line by the coordinates of the position to place the dimension numeral. The absolute value of the dimension is used in the calculation. If the dimension is negative AND not an integer, the coordinates of the position must be followed by the number of decimal places to display the dimension. sigma: lists the estimates for Rmaximum +/- sigma, taken as the radius when the curve passes through exp(-1/2). xyin: input to xylop, the plot output: messages to the user description Create a graph of radius versus density of Shannon spheres at various given dimensions. The output is run through xyplo. The function is: pd(R) = R^(D-1) * exp( sqr(R)/ (2* sqr(sigma))) where '^' means to exponentiate and where sqr(sigma) * (D-1) - sqr(Rmaximum) so setting Rmaximum = 1 relates sigma and D. The graph is in the range (0,0) to (r=maxr,1)). The curve is normalized so that its maximum is at (1,1). (except when dimension = 1, where it is at (1,0). Since xyplo can't plot several separate curves, without being told each symbol, this program simply starts at (0,pd(r)), draws the curve to (maxr,pd(maxr)), then circles back by drawing lines to the x axis (2*maxr,0) and then the origin (0,0). By setting the region that xyplo plots below maxr, one gets nice, fully correct curves that do not appear to be connected. documentation [1988 jan 23,5] @article{Schneider.ccmm, author = "T. D. Schneider", title = "Theory of Molecular Machines. {I. Channel} Capacity of Molecular Machines", journal = "J. Theor. Biol.", volume = "148", thenumber = "1", pages = "83-123", comment = "{(Note: The figures were printed out of order! Fig. 1 is on p. 97)}", note = "http://www-lecb.ncifcrf.gov/$\sim$toms/paper/ccmm/", year = 1991} see also {Schneider.ccmm paper:} http://www.lecb.ncifcrf.gov/~toms/paper/ccmm {example parameters:} spherep {example plotting file for xyplo:} sphere.xyplop {plotting program:} xyplo.p {resulting graph, postscript ... :} sphereinteger.ps {resulting graph, jpg .......... :} sphereinteger.jpg {example of fractal graph parameters:} spherep.fractal {resulting graph, postscript ... :} spherefractal.ps {resulting graph, jpg .......... :} spherefractal.jpg {Related programs:} {compress D dimensional sphere to 2D: ...} ring.p {plot output of ring: ...................} riden.p {Match fdr curve: .......................} fdr.p author Thomas Dana Schneider bugs none known */ /* end module describe.sphere */ /* begin module sphere.const */ #define wid 10 /* width of numbers */ #define dec 8 /* number of decimal places of numbers */ #define Rmaximum 1.0 /* where the curves are to have their maximum */ /* end module sphere.const */ Static _TEXT spherep; /* parameters */ Static _TEXT sigma; /* list of sigma values */ Static _TEXT xyin; /* input to xyplo */ #define tolerance 5e-2 /* how close we have to be to be "near" */ Static boolean near_(n, spot) double n, spot; { /* return true if n is near spot */ return (n > spot - tolerance && n < spot + tolerance); } #undef tolerance 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) double d, 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); } /* begin module sphere.themain */ Static Void themain(spherep, sigma, xyin) _TEXT *spherep, *sigma, *xyin; { /* the main procedure of the program */ double D; /* dimension, fractal dimension allowed! */ double Dx; /* x coordinate to plot the value of D */ double Dy; /* y coordinate to plot the value of D */ long decimals; /* number of decimal places to display the dimension */ double exphalf = exp(-1.0 / 2); /* exponent of -1/2 */ long interval; /* what segment of the curve we are currently plotting. 1: before first intercept with exphalf, 2: before rmax has been reached, 3: before the second intercept with exphalf has been reached, 4: tail of the curve */ double lastvalue; /* the previous value of the curve */ double lnpdmax; /* ln of maximum value of pd curve, for normalizing */ double lnvalue; /* intermediate value: ln of the final result */ double maxr; /* maximum radius to calculate to */ double pdmax; /* maximum value of pd curve, for normalizing */ double r; /* current radius */ double sigmaplus, sigmaminus; /* estimates of sigma above and below rmax */ double step; /* the step between calculations */ double twosigma2; /* 2 times sigma squared, where sigma is the standard deviation of the gaussian */ double value; /* the final result */ printf("sphere %4.2f\n", version); if (*sigma->name != '\0') { if (sigma->f != NULL) sigma->f = freopen(sigma->name, "w", sigma->f); else sigma->f = fopen(sigma->name, "w"); } else { if (sigma->f != NULL) rewind(sigma->f); else sigma->f = tmpfile(); } if (sigma->f == NULL) _EscIO2(FileNotFound, sigma->name); SETUPBUF(sigma->f, Char); fprintf(sigma->f, "* sphere %4.2f\n", version); fprintf(sigma->f, "* D: dimension\n"); fprintf(sigma->f, "* sigma-: deviation to inside of sphere\n"); fprintf(sigma->f, "* sigma+: deviation to outside of sphere\n"); fprintf(sigma->f, "* average: average of absolute values of the sigmas\n"); fprintf(sigma->f, "* estimate: of sigma (1/sqrt(2*(D-1))\n"); fprintf(sigma->f, "*\n"); if (*spherep->name != '\0') { if (spherep->f != NULL) spherep->f = freopen(spherep->name, "r", spherep->f); else spherep->f = fopen(spherep->name, "r"); } else rewind(spherep->f); if (spherep->f == NULL) _EscIO2(FileNotFound, spherep->name); RESETBUF(spherep->f, Char); fscanf(spherep->f, "%lg%*[^\n]", &step); getc(spherep->f); if (step <= 0) { step = 0.2; printf("step must be positive, using %*.*f\n", wid, dec, step); } fscanf(spherep->f, "%lg%*[^\n]", &maxr); getc(spherep->f); if (*xyin->name != '\0') { if (xyin->f != NULL) xyin->f = freopen(xyin->name, "w", xyin->f); else xyin->f = fopen(xyin->name, "w"); } else { if (xyin->f != NULL) rewind(xyin->f); else xyin->f = tmpfile(); } if (xyin->f == NULL) _EscIO2(FileNotFound, xyin->name); SETUPBUF(xyin->f, Char); fprintf(xyin->f, "* sphere %4.2f\n", version); fprintf(xyin->f, "* step: %*.*f\n", wid, dec, step); fprintf(xyin->f, "* maximum r: %*.*f\n", wid, dec, maxr); fprintf(xyin->f, "* kind (d=dimension numeral, c = curve) r, pd(r)\n"); while (!BUFEOF(spherep->f)) { fscanf(spherep->f, "%lg", &D); if (D < 0) { /* plot the dimension number at the coordinates given */ fscanf(spherep->f, "%lg%lg", &Dx, &Dy); D = fabs(D); if (D == (long)floor(D + 0.5)) decimals = 0; else fscanf(spherep->f, "%ld", &decimals); fscanf(spherep->f, "%*[^\n]"); getc(spherep->f); if (D == (long)floor(D + 0.5)) { fprintf(xyin->f, "D=%ld %*.*f %*.*f\n", (long)floor(D + 0.5), wid, dec, Dx, wid, dec, Dy); printf("D = %ld plotted at %*.*f %*.*f\n", (long)floor(D + 0.5), wid, dec, Dx, wid, dec, Dy); } else { fprintf(xyin->f, "D=%1.*f %*.*f %*.*f\n", (int)decimals, D, wid, dec, Dx, wid, dec, Dy); printf("D = %1.*f plotted at %*.*f %*.*f\n", (int)decimals, D, wid, dec, Dx, wid, dec, Dy); } } else { /*zzz*/ /* writeln(output,'decimals =',decimals:1); writeln(output,'round(D) =',round(D):1); writeln(output,' D =',D:1:dec); if D = round(D) then writeln(output,'D = round(D) SAME') else writeln(output,'D <> round(D) DIFFERENT: fractal!'); writeln(output); */ fscanf(spherep->f, "%*[^\n]"); getc(spherep->f); if (D == (long)floor(D + 0.5)) printf("D = % .*E\n", P_max(wid - 7, 1), D); else printf("D = %*.*f\n", wid, dec, D); } fprintf(sigma->f, "D:% .1E ", D); 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(pdmax); } /* make gaussian */ sigmaplus = 0.0; sigmaminus = 0.0; interval = 1; lastvalue = -LONG_MAX; /* preset for the way up */ r = step; while (r <= maxr) { /* writeln(xyin,r:wid:dec,' ', (pd(D,r, twosigma2)/pdmax):wid:dec) */ lnvalue = lnpd(D, r, twosigma2) - lnpdmax; /* by testing for the log of this value, we can avoid exponentiation of very small values. the result is not be affected, for these values come up only when the final result is nearly zero anyway. Just don't even plot them! */ if (lnvalue >= -25) { value = exp(lnvalue); fprintf(xyin->f, "c %*.*f %*.*f\n", wid, dec, r, wid, dec, value); switch (interval) { case 1: /* before first intercept */ if (value > exphalf) { sigmaminus = r - step + step * (exphalf - lastvalue) / (value - lastvalue) - Rmaximum; lastvalue = LONG_MAX; /* reset for the way down */ interval++; } else lastvalue = value; break; case 2: /* before peak */ /* wait for the maximum of the curve */ if (r >= Rmaximum) interval++; break; case 3: /* before second intercept */ if (value < exphalf) { sigmaplus = r + step + step * (value - exphalf) / (lastvalue - value) - Rmaximum; fprintf(sigma->f, " sigma-: %5.3f sigma+: %5.3f average: %5.3f", sigmaminus, sigmaplus, (sigmaplus - sigmaminus) / 2.0); if (D > 1) fprintf(sigma->f, " estimate: %5.3f\n", 1 / sqrt(2 * (D - 1))); else fprintf(sigma->f, " exactly!: %5.3f\n", exphalf); interval++; } else lastvalue = value; break; case 4: /* after second intercept */ break; } } r += step; } /* make a line to the axis, off the edge of the plot: */ fprintf(xyin->f, "c %*.*f %*.*f", wid, dec, 2 * maxr, wid, dec, 0.0); fprintf(xyin->f, " move to axis\n"); /* move back to the origin along the axis: */ fprintf(xyin->f, "c %*.*f %*.*f", wid, dec, 0.0, wid, dec, 0.0); fprintf(xyin->f, " move to origin\n"); } } /* end module sphere.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); xyin.f = NULL; strcpy(xyin.name, "xyin"); sigma.f = NULL; strcpy(sigma.name, "sigma"); spherep.f = NULL; strcpy(spherep.name, "spherep"); themain(&spherep, &sigma, &xyin); _L1: if (spherep.f != NULL) fclose(spherep.f); if (sigma.f != NULL) fclose(sigma.f); if (xyin.f != NULL) fclose(xyin.f); exit(EXIT_SUCCESS); } /* End. */