/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "riden.p" */ #include /* riden: ring density graph Tom Schneider National Cancer Institute Laboratory of Mathematical Biology Frederick, Maryland 21701-1013 toms@ncifcrf.gov 1989 */ /* end of program */ /* begin module version */ #define version 1.29 /* of riden.p 1996 November 18 origin 1989 November 19 */ /* end module version */ /* begin module describe.riden */ /* name riden: ring density graph synopsis riden(color: in, ridenp: in, xyin: out, output: out) files color: output of the ring program ridenp: parameter file for this program. Two lines: First line: largest radial distance recorded Second line: number of bins to store the data in xyin: histogram of the density output: messages to the user description This program converts the graph generated by the ring program into a form that allows one to see if the results are as expected. examples documentation see also ring.p author Thomas Dana Schneider bugs Program only works for D=2. The curves don't match for D=4, but do for higher dimensions. It is not obvious why. technical notes */ /* end module describe.riden */ /* more constants */ #define binsmax 1000 /* maximum number of storage locations */ #define Rmaximum 1.0 /* where the curves are to have their maximum */ /* an array for storing integers */ typedef long storage[binsmax + 1]; Static _TEXT color; /* output of the ring program */ Static _TEXT ridenp; /* output, input to xyplo */ Static _TEXT xyin; /* parameter file of this program */ Static jmp_buf _JL1; /* 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 = 2.46; (@ of ring.p 1989 Nov 19 */ /* 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 = 2.46; (@ of ring.p 1989 Nov 19 */ Static double pd(d, r, twosigma2) long d; double r, twosigma2; { /* probability at a given radius in a given dimension */ if (r == 0.0) return 0.0; else return (exp(-(r * r / twosigma2)) * exp((d - 1) * log(r))); } /* begin module riden.themain */ Static Void themain(color, ridenp, xyin) _TEXT *color, *ridenp, *xyin; { /* the main procedure of the program */ long biggest = 0; /* the largest point in the histogram */ long bins; /* number of divisions of rmax to make */ long D; /* the dimensionality */ double pdmax; /* the largest expected density value */ double r; /* radial distance from the origin coordinate of the point */ double rmax; /* maximum radial distance to record */ long s; /* index for finding the dimensionality, and indexing the store */ long skipped = 0; /* number of points lost */ long stored = 0; /* number of points captured */ storage store; /* the place to build the histogram */ double twosigma2; /* 2 times sigma squared, where sigma is the standard deviation of the gaussian */ double x; /* x coordinate of the point */ double y; /* y coordinate of the point */ printf("riden %4.2f\n", version); if (*ridenp->name != '\0') { if (ridenp->f != NULL) ridenp->f = freopen(ridenp->name, "r", ridenp->f); else ridenp->f = fopen(ridenp->name, "r"); } else rewind(ridenp->f); if (ridenp->f == NULL) _EscIO2(FileNotFound, ridenp->name); RESETBUF(ridenp->f, Char); fscanf(ridenp->f, "%lg%*[^\n]", &rmax); getc(ridenp->f); fscanf(ridenp->f, "%ld%*[^\n]", &bins); getc(ridenp->f); if (bins > binsmax) { printf("too many bins, increase binsmax\n"); halt(); } /* determine the dimensonality */ if (*color->name != '\0') { if (color->f != NULL) color->f = freopen(color->name, "r", color->f); else color->f = fopen(color->name, "r"); } else rewind(color->f); if (color->f == NULL) _EscIO2(FileNotFound, color->name); RESETBUF(color->f, Char); for (s = 1; s <= 9; s++) { fscanf(color->f, "%*[^\n]"); getc(color->f); } getc(color->f); fscanf(color->f, "%ld%*[^\n]", &D); getc(color->f); if (D < 2) { printf("D must be > 1\n"); halt(); } printf("dimensionality: D = %ld\n", D); 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, "* riden %4.2f\n", version); if (*color->name != '\0') { if (color->f != NULL) color->f = freopen(color->name, "r", color->f); else color->f = fopen(color->name, "r"); } else rewind(color->f); if (color->f == NULL) _EscIO2(FileNotFound, color->name); RESETBUF(color->f, Char); while (P_peek(color->f) == '*') copyaline(color, xyin); /* clear the store */ for (s = 0; s <= bins; s++) store[s] = 0; /* build the histogram */ while (!BUFEOF(color->f)) { if (P_peek(color->f) != 's') { /* pick up the simulated points only */ fscanf(color->f, "%*[^\n]"); getc(color->f); continue; } getc(color->f); /* skip the s */ fscanf(color->f, "%lg%lg%*[^\n]", &x, &y); getc(color->f); /* grab the coordinates */ r = sqrt(x * x + y * y); s = (long)floor(bins / rmax * r + 0.5); if (s <= bins) store[s]++; else { skipped++; printf("skipping r=%10.8f\n", r); } } /* determine the largest stored value, and the number stored */ for (s = 0; s <= bins; s++) { if (store[s] > biggest) biggest = store[s]; stored += store[s]; } /* (* #1 *) twosigma2 := 2*sqr(Rmaximum)/sqrt(D-1); (* this worked with #1! *) (* ie, the fdr curve (4d) matched the data, but both were shifted high *) */ /* #2 */ twosigma2 = 2.0 / (D - 1); /* this forces the fdr curve to pass through 1 */ pdmax = pd(D, 1.0, twosigma2); /* write it out */ fprintf(xyin->f, "* symbol, radius, normalized histogram, raw histogram, expected fD(r)\n"); /* identification tag */ fprintf(xyin->f, "D=%ld %10.8f %10.8f %5d %10d\n", D, 2.0, 0.9, 0, 0); for (s = 0; s <= bins; s++) { r = rmax * s / bins; fprintf(xyin->f, "s %10.8f %10.8f %5ld %10.8f\n", r, (double)store[s] / biggest, store[s], pd(D, r, twosigma2) / pdmax); } /* now give the expected data in the normalized column to allow plotting */ for (s = 0; s <= bins; s++) { r = rmax * s / bins; fprintf(xyin->f, "e %10.8f %10.8f %5ld %10.8f\n", r, pd(D, r, twosigma2) / pdmax, store[s], (double)store[s] / biggest); } printf("number of points skipped: %5ld\n", skipped); printf("number of points histogrammed: %5ld\n", stored); printf("number of points in input file: %5ld\n", skipped + stored); } /* end module riden.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; xyin.f = NULL; strcpy(xyin.name, "xyin"); ridenp.f = NULL; strcpy(ridenp.name, "ridenp"); color.f = NULL; strcpy(color.name, "color"); themain(&color, &ridenp, &xyin); _L1: if (color.f != NULL) fclose(color.f); if (ridenp.f != NULL) fclose(ridenp.f); if (xyin.f != NULL) fclose(xyin.f); exit(EXIT_SUCCESS); } /* End. */