/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "fdr.p" */ #include /* computes f_D(r) */ /* begin module version */ #define version 1.04 /* of fdr.p 1995 Jan 14 */ /* Revision history: v. 1.00 started = 10 Jan 1995 Written at National Cancer Institute, Frederick, Maryland */ /* end module version */ /* begin module describe.fdr */ /* name fdr: computes best values of sigma/D for f_D(r) fit of real data. synopsis fdr(histog: in, fdrp:in, outfile: out, output: out) files histog: output file from genhis. fdrp: parameter file for fdr. sigma,sigmastep (reals): initial value and incremental increase in sigma. dim,dimstep (reals): initial value and incremental increase in D. temperature (real): Kelvin temperature of system. outfile: writes best value of sigma/D and error in f_D(r) fit. output: sends messages to the user. description D-1 2 2 r exp(-r /2sigma ) f (r) = ----------------------- D D D/2-1 Gamma(D/2) sigma 2 We compute the values of D, sigma that minimze the least-squares error between f_D(r) and the actual data read in from histog. examples See fdr.eg (not yet written ) documentation See Appendix 3 to TDS's CCMM paper. see also genhis.p author V. Jejjala (vishnu@ncifcrf.gov) bugs Unknown but existent (existential bugs, as it were... :-) ). If you find a bug, please forward your discovery to the author. disclaimer The author of this code is in no way responsible for any havoc that this program may cause to your computer system, your files, or your reputation. Use this program and the data it generates at your own risk. technical notes Use large data sets in histog when possible. */ /* end module describe.fdr */ /* begin module globals */ /* more constants */ #define k 1.380662e-23 /* Boltzmann's constant */ #define maxval 50 /* maximum number of points to consider on f_D(r) */ #define maxpass 100 /* maximum number of passes of program */ #define maxpassp 101 /* maxpass + 1 */ #define big 1000000.0 /* very big real number */ #define machep (-500.0) /* exp(machep) ~ machine epsilon */ typedef struct singlept { double start, histval; } singlept; typedef struct recdata { singlept allpts[maxval]; double sig, D, err; } recdata; typedef double rdata[maxpassp + 1]; Static _TEXT histog; /* output of genhis from which data are read */ Static _TEXT fdrp; /* parameter file */ Static _TEXT outfile; /* file to which data are written */ /* end module globals */ /* begin module expo */ Static double expo_(a, b) double a, b; { /* yields a^b */ if (a == 0.0) return 0.0; else if (b * log(a) < machep) return 0.0; else return exp(b * log(a)); } #define stp 2.50662827465 /* constant to compute gamma */ /* end module expo */ /* begin module gamma */ Static double gamma(xx) double xx; { /* yields gamma(xx) using expansion formula/recurrance relation; code modified from _Numerical Recipes_ */ double x, tmp, ser; /* temporary variables */ double gammaln; /* ln(gamma(x)) */ if (xx < 1) x = xx; else x = xx - 1.0; tmp = x + 5.5; tmp = (x + 0.5) * log(tmp) - tmp; ser = 1.0 + 76.18009173 / (x + 1.0) - 86.50532033 / (x + 2.0) + 24.01409822 / (x + 3.0) - 1.231739516 / (x + 4.0) + 0.120858003e-2 / (x + 5.0) - 0.536382e-5 / (x + 6.0); gammaln = tmp + log(stp * ser); if (xx < 1) return (exp(gammaln) / xx); else return exp(gammaln); } #undef stp /* end module gamma */ /* begin module f */ Static double f(r, sigma, dim) double r, sigma, dim; { /* computes f_D(r) */ double g, h; /* temporary variables */ if (-1 * expo_(r, 2.0) / (2 * expo_(sigma, 2.0)) < machep) g = 0.0; else g = exp(-1 * expo_(r, 2.0) / (2 * expo_(sigma, 2.0))); h = gamma(dim / 2) * expo_(sigma, dim) * expo_(2.0, dim / 2 - 1); return (expo_(r, dim - 1) * g / h); } /* end module f */ /* begin module readdata */ Static Void readdata(histog, datareal, Rseq, Rstart, Rstep) _TEXT *histog; recdata *datareal; double *Rseq, *Rstart, *Rstep; { /* reads in data from histog */ double total = 0.0; /* total number of sequences */ long dummy; /* dummy variable */ long count1 = 0; long count2; /* for loops */ if (*histog->name != '\0') { if (histog->f != NULL) histog->f = freopen(histog->name, "r", histog->f); else histog->f = fopen(histog->name, "r"); } else rewind(histog->f); if (histog->f == NULL) _EscIO2(FileNotFound, histog->name); RESETBUF(histog->f, Char); while (!BUFEOF(histog->f)) { if (P_peek(histog->f) == '*') { fscanf(histog->f, "%*[^\n]"); getc(histog->f); continue; } count1++; fscanf(histog->f, "%lg%ld%*[^\n]", &datareal->allpts[count1-1].start, &dummy); getc(histog->f); datareal->allpts[count1-1].histval = dummy; total = dummy + total; } for (count2 = 0; count2 < count1; count2++) { *Rseq = datareal->allpts[count2].histval * datareal->allpts[count2].start + *Rseq; datareal->allpts[count2].histval /= total; } *Rseq /= total; *Rstart = datareal->allpts[0].start; *Rstep = datareal->allpts[1].start - datareal->allpts[0].start; for (count2 = count1; count2 < maxval; count2++) { datareal->allpts[count2].start = *Rstart + count2 * *Rstep; datareal->allpts[count2].histval = 0.0; } } /* end module readdata */ /* begin module readpar */ Static Void readpar(fdrp, sigma, sigmastep, dim, dimstep, temp) _TEXT *fdrp; double *sigma, *sigmastep, *dim, *dimstep, *temp; { /* reads in parameters from fdrp */ if (*fdrp->name != '\0') { if (fdrp->f != NULL) fdrp->f = freopen(fdrp->name, "r", fdrp->f); else fdrp->f = fopen(fdrp->name, "r"); } else rewind(fdrp->f); if (fdrp->f == NULL) _EscIO2(FileNotFound, fdrp->name); RESETBUF(fdrp->f, Char); fscanf(fdrp->f, "%lg%lg%*[^\n]", sigma, sigmastep); getc(fdrp->f); fscanf(fdrp->f, "%lg%lg%*[^\n]", dim, dimstep); getc(fdrp->f); fscanf(fdrp->f, "%lg%*[^\n]", temp); getc(fdrp->f); } /* end module readpar */ /* begin module computeerror */ Static Void computeerror(databest, error, datareal, sigma, dim, kT, Rstart, Rstep, sigmastep, dimstep, pass) recdata *databest; double *error; recdata datareal; double sigma, dim, kT, Rstart, Rstep, sigmastep, dimstep; long pass; { /* computes least-squares error--(sigma_(f_D(r)))^2 */ recdata datameas; /* f_D(r) data for this particular value of D,sig,r */ double r; /* current value of r */ double diff; /* difference between f_D(r) and actual curve */ double err; /* least-squares error in f_D(r)--sum of diff over all r */ long count; /* for loops */ datameas.sig = sigma; datameas.D = dim; for (count = 0; count < maxval; count++) { r = Rstart + count * Rstep; datameas.allpts[count].start = r; datameas.allpts[count].histval = f(r, sigma, dim) / kT; diff = datareal.allpts[count].histval - datameas.allpts[count].histval; if (diff < 0) diff = -1 * diff; err = expo_(diff, 2.0) + err; } datameas.err = err / maxval; /*temp below*/ printf("% .5E% .5E% .5E%12d%12d% .5E% .5E\n", sigma, dim, datameas.err / 225 * 0.84 + 0.16, 1, 1, sigmastep, dimstep); /*temp above*/ if (datameas.err < error[pass]) error[pass] = datameas.err; if (datameas.err >= databest->err) return; for (count = 0; count < maxval; count++) { databest->allpts[count].start = datameas.allpts[count].start; databest->allpts[count].histval = datameas.allpts[count].histval; } databest->sig = datameas.sig; databest->D = datameas.D; databest->err = datameas.err; } /* end module computeerror */ /* begin module writefile */ Static Void writefile(outfile, databest, sigmastep, dimstep, pass) _TEXT *outfile; recdata databest; double sigmastep, dimstep; long pass; { /* writes out best sigma/D to outfile */ long count; /* for loops */ fprintf(outfile->f, "* Sigma: % .5E +/-% .5E\n", databest.sig, sigmastep); fprintf(outfile->f, "* Dimension: % .5E +/-% .5E\n", databest.D, dimstep); fprintf(outfile->f, "* Std. Deviation: % .5E\n", expo_(databest.err, 0.5)); fprintf(outfile->f, "* in %12ld passes\n", pass); for (count = 0; count < maxval; count++) { fprintf(outfile->f, "% .5E", databest.allpts[count].start); fprintf(outfile->f, "% .5E\n", databest.allpts[count].histval); } } /* end module writefile */ /* begin module themain */ Static Void themain(histog, fdrp, outfile) _TEXT *histog, *fdrp, *outfile; { /* main procedure of fdr.p */ recdata datareal, databest; /* actual and best-fit data */ rdata error; /* array of data at various values o dim/sigma */ double Rseq; /* rough estimate of R_seq */ double Rstart, Rstep; /* initial value and steps in gradation of r */ double sigma, sigmastep; /* initial value and steps in gradation of sigma */ double dim, dimstep; /* initial value and steps in gradation of D */ double dime; /* actual value of D used at any given time */ double temp, kT; /* temperature and k*temperature */ long pass = 1; /* number of tries done */ long count; /* for loops */ boolean result = false; /* decides when we're done */ long FORLIM; if (*outfile->name != '\0') { if (outfile->f != NULL) outfile->f = freopen(outfile->name, "w", outfile->f); else outfile->f = fopen(outfile->name, "w"); } else { if (outfile->f != NULL) rewind(outfile->f); else outfile->f = tmpfile(); } if (outfile->f == NULL) _EscIO2(FileNotFound, outfile->name); SETUPBUF(outfile->f, Char); readdata(histog, &datareal, &Rseq, &Rstart, &Rstep); databest.err = big; error[0] = big; do { if (pass == 1) { readpar(fdrp, &sigma, &sigmastep, &dim, &dimstep, &temp); kT = 1.0; } dime = dim; sigma += sigmastep; FORLIM = (long)floor(Rseq + 0.5) * 3; for (count = 0; count <= FORLIM; count++) { dime += dimstep; error[pass] = big; computeerror(&databest, error, datareal, sigma, dime, kT, Rstart, Rstep, sigmastep, dimstep, pass); } if (error[pass] < error[pass-1] && pass > 1 || pass > maxpass) { writefile(outfile, databest, sigmastep, dimstep, pass); result = true; } pass++; } while (result != true); } /* end module themain */ /* begin module mainroutine */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); outfile.f = NULL; strcpy(outfile.name, "outfile"); fdrp.f = NULL; strcpy(fdrp.name, "fdrp"); histog.f = NULL; strcpy(histog.name, "histog"); themain(&histog, &fdrp, &outfile); if (histog.f != NULL) fclose(histog.f); if (fdrp.f != NULL) fclose(fdrp.f); if (outfile.f != NULL) fclose(outfile.f); exit(EXIT_SUCCESS); } /* end module mainroutine */ /* End. */