/* * Experimental data distribution table generator * Taken from the uncopyrighted NISTnet code (public domain). * * Read in a series of "random" data values, either * experimentally or generated from some probability distribution. * From this, create the inverse distribution table used to approximate * the distribution. */ #include #include #include #include #include #include #include double * readdoubles(FILE *fp, int *number) { struct stat info; double *x; int limit; int n=0, i; if (!fstat(fileno(fp), &info) && info.st_size > 0) { limit = 2*info.st_size/sizeof(double); /* @@ approximate */ } else { limit = 10000; } x = calloc(limit, sizeof(double)); if (!x) { perror("double alloc"); exit(3); } for (i=0; i= DISTTABLESIZE) index = DISTTABLESIZE-1; ++table[index]; if (index > last) last = index +1; if (index < first) first = index; } return table; } /* replace an array by its cumulative distribution */ static void cumulativedist(int *table, int limit, int *total) { int accum=0; while (--limit >= 0) { accum += *table; *table++ = accum; } *total = accum; } static short * inverttable(int *table, int inversesize, int tablesize, int cumulative) { int i, inverseindex, inversevalue; short *inverse; double findex, fvalue; inverse = (short *)malloc(inversesize*sizeof(short)); for (i=0; i < inversesize; ++i) { inverse[i] = MINSHORT; } for (i=0; i < tablesize; ++i) { findex = ((double)i/(double)DISTTABLEGRANULARITY) - DISTTABLEDOMAIN; fvalue = (double)table[i]/(double)cumulative; inverseindex = (int)rint(fvalue*inversesize); inversevalue = (int)rint(findex*TABLEFACTOR); if (inversevalue <= MINSHORT) inversevalue = MINSHORT+1; if (inversevalue > MAXSHORT) inversevalue = MAXSHORT; if (inverseindex >= inversesize) inverseindex = inversesize- 1; inverse[inverseindex] = inversevalue; } return inverse; } /* Run simple linear interpolation over the table to fill in missing entries */ static void interpolatetable(short *table, int limit) { int i, j, last, lasti = -1; last = MINSHORT; for (i=0; i < limit; ++i) { if (table[i] == MINSHORT) { for (j=i; j < limit; ++j) if (table[j] != MINSHORT) break; if (j < limit) { table[i] = last + (i-lasti)*(table[j]-last)/(j-lasti); } else { table[i] = last + (i-lasti)*(MAXSHORT-last)/(limit-lasti); } } else { last = table[i]; lasti = i; } } } static void printtable(const short *table, int limit) { int i; printf("# This is the distribution table for the experimental distribution.\n"); for (i=0 ; i < limit; ++i) { printf("%d%c", table[i], (i % 8) == 7 ? '\n' : ' '); } } int main(int argc, char **argv) { FILE *fp; double *x; double mu, sigma, rho; int limit; int *table; short *inverse; int total; if (argc > 1) { if (!(fp = fopen(argv[1], "r"))) { perror(argv[1]); exit(1); } } else { fp = stdin; } x = readdoubles(fp, &limit); if (limit <= 0) { fprintf(stderr, "Nothing much read!\n"); exit(2); } arraystats(x, limit, &mu, &sigma, &rho); #ifdef DEBUG fprintf(stderr, "%d values, mu %10.4f, sigma %10.4f, rho %10.4f\n", limit, mu, sigma, rho); #endif table = makedist(x, limit, mu, sigma); free((void *) x); cumulativedist(table, DISTTABLESIZE, &total); inverse = inverttable(table, TABLESIZE, DISTTABLESIZE, total); interpolatetable(inverse, TABLESIZE); printtable(inverse, TABLESIZE); return 0; }