|
Packit Service |
3880ab |
/*
|
|
Packit Service |
3880ab |
* Experimental data distribution table generator
|
|
Packit Service |
3880ab |
* Taken from the uncopyrighted NISTnet code (public domain).
|
|
Packit Service |
3880ab |
*
|
|
Packit Service |
3880ab |
* Read in a series of "random" data values, either
|
|
Packit Service |
3880ab |
* experimentally or generated from some probability distribution.
|
|
Packit Service |
3880ab |
* From this, create the inverse distribution table used to approximate
|
|
Packit Service |
3880ab |
* the distribution.
|
|
Packit Service |
3880ab |
*/
|
|
Packit Service |
3880ab |
#include <stdio.h>
|
|
Packit Service |
3880ab |
#include <stdlib.h>
|
|
Packit Service |
3880ab |
#include <math.h>
|
|
Packit Service |
3880ab |
#include <malloc.h>
|
|
Packit Service |
3880ab |
#include <string.h>
|
|
Packit Service |
3880ab |
#include <sys/types.h>
|
|
Packit Service |
3880ab |
#include <sys/stat.h>
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
double *
|
|
Packit Service |
3880ab |
readdoubles(FILE *fp, int *number)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
struct stat info;
|
|
Packit Service |
3880ab |
double *x;
|
|
Packit Service |
3880ab |
int limit;
|
|
Packit Service |
3880ab |
int n=0, i;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
if (!fstat(fileno(fp), &info) &&
|
|
Packit Service |
3880ab |
info.st_size > 0) {
|
|
Packit Service |
3880ab |
limit = 2*info.st_size/sizeof(double); /* @@ approximate */
|
|
Packit Service |
3880ab |
} else {
|
|
Packit Service |
3880ab |
limit = 10000;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
x = calloc(limit, sizeof(double));
|
|
Packit Service |
3880ab |
if (!x) {
|
|
Packit Service |
3880ab |
perror("double alloc");
|
|
Packit Service |
3880ab |
exit(3);
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
for (i=0; i
|
|
Packit Service |
3880ab |
if (fscanf(fp, "%lf", &x[i]) != 1 ||
|
|
Packit Service |
3880ab |
feof(fp))
|
|
Packit Service |
3880ab |
break;
|
|
Packit Service |
3880ab |
++n;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
*number = n;
|
|
Packit Service |
3880ab |
return x;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
void
|
|
Packit Service |
3880ab |
arraystats(double *x, int limit, double *mu, double *sigma, double *rho)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
int n=0, i;
|
|
Packit Service |
3880ab |
double sumsquare=0.0, sum=0.0, top=0.0;
|
|
Packit Service |
3880ab |
double sigma2=0.0;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
for (i=0; i
|
|
Packit Service |
3880ab |
sumsquare += x[i]*x[i];
|
|
Packit Service |
3880ab |
sum += x[i];
|
|
Packit Service |
3880ab |
++n;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
*mu = sum/(double)n;
|
|
Packit Service |
3880ab |
*sigma = sqrt((sumsquare - (double)n*(*mu)*(*mu))/(double)(n-1));
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
for (i=1; i < n; ++i){
|
|
Packit Service |
3880ab |
top += ((double)x[i]- *mu)*((double)x[i-1]- *mu);
|
|
Packit Service |
3880ab |
sigma2 += ((double)x[i-1] - *mu)*((double)x[i-1] - *mu);
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
*rho = top/sigma2;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
/* Create a (normalized) distribution table from a set of observed
|
|
Packit Service |
3880ab |
* values. The table is fixed to run from (as it happens) -4 to +4,
|
|
Packit Service |
3880ab |
* with granularity .00002.
|
|
Packit Service |
3880ab |
*/
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
#define TABLESIZE 16384/4
|
|
Packit Service |
3880ab |
#define TABLEFACTOR 8192
|
|
Packit Service |
3880ab |
#ifndef MINSHORT
|
|
Packit Service |
3880ab |
#define MINSHORT -32768
|
|
Packit Service |
3880ab |
#define MAXSHORT 32767
|
|
Packit Service |
3880ab |
#endif
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
/* Since entries in the inverse are scaled by TABLEFACTOR, and can't be bigger
|
|
Packit Service |
3880ab |
* than MAXSHORT, we don't bother looking at a larger domain than this:
|
|
Packit Service |
3880ab |
*/
|
|
Packit Service |
3880ab |
#define DISTTABLEDOMAIN ((MAXSHORT/TABLEFACTOR)+1)
|
|
Packit Service |
3880ab |
#define DISTTABLEGRANULARITY 50000
|
|
Packit Service |
3880ab |
#define DISTTABLESIZE (DISTTABLEDOMAIN*DISTTABLEGRANULARITY*2)
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
static int *
|
|
Packit Service |
3880ab |
makedist(double *x, int limit, double mu, double sigma)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
int *table;
|
|
Packit Service |
3880ab |
int i, index, first=DISTTABLESIZE, last=0;
|
|
Packit Service |
3880ab |
double input;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
table = calloc(DISTTABLESIZE, sizeof(int));
|
|
Packit Service |
3880ab |
if (!table) {
|
|
Packit Service |
3880ab |
perror("table alloc");
|
|
Packit Service |
3880ab |
exit(3);
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
for (i=0; i < limit; ++i) {
|
|
Packit Service |
3880ab |
/* Normalize value */
|
|
Packit Service |
3880ab |
input = (x[i]-mu)/sigma;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
index = (int)rint((input+DISTTABLEDOMAIN)*DISTTABLEGRANULARITY);
|
|
Packit Service |
3880ab |
if (index < 0) index = 0;
|
|
Packit Service |
3880ab |
if (index >= DISTTABLESIZE) index = DISTTABLESIZE-1;
|
|
Packit Service |
3880ab |
++table[index];
|
|
Packit Service |
3880ab |
if (index > last)
|
|
Packit Service |
3880ab |
last = index +1;
|
|
Packit Service |
3880ab |
if (index < first)
|
|
Packit Service |
3880ab |
first = index;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
return table;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
/* replace an array by its cumulative distribution */
|
|
Packit Service |
3880ab |
static void
|
|
Packit Service |
3880ab |
cumulativedist(int *table, int limit, int *total)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
int accum=0;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
while (--limit >= 0) {
|
|
Packit Service |
3880ab |
accum += *table;
|
|
Packit Service |
3880ab |
*table++ = accum;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
*total = accum;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
static short *
|
|
Packit Service |
3880ab |
inverttable(int *table, int inversesize, int tablesize, int cumulative)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
int i, inverseindex, inversevalue;
|
|
Packit Service |
3880ab |
short *inverse;
|
|
Packit Service |
3880ab |
double findex, fvalue;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
inverse = (short *)malloc(inversesize*sizeof(short));
|
|
Packit Service |
3880ab |
for (i=0; i < inversesize; ++i) {
|
|
Packit Service |
3880ab |
inverse[i] = MINSHORT;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
for (i=0; i < tablesize; ++i) {
|
|
Packit Service |
3880ab |
findex = ((double)i/(double)DISTTABLEGRANULARITY) - DISTTABLEDOMAIN;
|
|
Packit Service |
3880ab |
fvalue = (double)table[i]/(double)cumulative;
|
|
Packit Service |
3880ab |
inverseindex = (int)rint(fvalue*inversesize);
|
|
Packit Service |
3880ab |
inversevalue = (int)rint(findex*TABLEFACTOR);
|
|
Packit Service |
3880ab |
if (inversevalue <= MINSHORT) inversevalue = MINSHORT+1;
|
|
Packit Service |
3880ab |
if (inversevalue > MAXSHORT) inversevalue = MAXSHORT;
|
|
Packit Service |
3880ab |
if (inverseindex >= inversesize) inverseindex = inversesize- 1;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
inverse[inverseindex] = inversevalue;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
return inverse;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
/* Run simple linear interpolation over the table to fill in missing entries */
|
|
Packit Service |
3880ab |
static void
|
|
Packit Service |
3880ab |
interpolatetable(short *table, int limit)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
int i, j, last, lasti = -1;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
last = MINSHORT;
|
|
Packit Service |
3880ab |
for (i=0; i < limit; ++i) {
|
|
Packit Service |
3880ab |
if (table[i] == MINSHORT) {
|
|
Packit Service |
3880ab |
for (j=i; j < limit; ++j)
|
|
Packit Service |
3880ab |
if (table[j] != MINSHORT)
|
|
Packit Service |
3880ab |
break;
|
|
Packit Service |
3880ab |
if (j < limit) {
|
|
Packit Service |
3880ab |
table[i] = last + (i-lasti)*(table[j]-last)/(j-lasti);
|
|
Packit Service |
3880ab |
} else {
|
|
Packit Service |
3880ab |
table[i] = last + (i-lasti)*(MAXSHORT-last)/(limit-lasti);
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
} else {
|
|
Packit Service |
3880ab |
last = table[i];
|
|
Packit Service |
3880ab |
lasti = i;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
static void
|
|
Packit Service |
3880ab |
printtable(const short *table, int limit)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
int i;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
printf("# This is the distribution table for the experimental distribution.\n");
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
for (i=0 ; i < limit; ++i) {
|
|
Packit Service |
3880ab |
printf("%d%c", table[i],
|
|
Packit Service |
3880ab |
(i % 8) == 7 ? '\n' : ' ');
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
int
|
|
Packit Service |
3880ab |
main(int argc, char **argv)
|
|
Packit Service |
3880ab |
{
|
|
Packit Service |
3880ab |
FILE *fp;
|
|
Packit Service |
3880ab |
double *x;
|
|
Packit Service |
3880ab |
double mu, sigma, rho;
|
|
Packit Service |
3880ab |
int limit;
|
|
Packit Service |
3880ab |
int *table;
|
|
Packit Service |
3880ab |
short *inverse;
|
|
Packit Service |
3880ab |
int total;
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
if (argc > 1) {
|
|
Packit Service |
3880ab |
if (!(fp = fopen(argv[1], "r"))) {
|
|
Packit Service |
3880ab |
perror(argv[1]);
|
|
Packit Service |
3880ab |
exit(1);
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
} else {
|
|
Packit Service |
3880ab |
fp = stdin;
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
x = readdoubles(fp, &limit);
|
|
Packit Service |
3880ab |
if (limit <= 0) {
|
|
Packit Service |
3880ab |
fprintf(stderr, "Nothing much read!\n");
|
|
Packit Service |
3880ab |
exit(2);
|
|
Packit Service |
3880ab |
}
|
|
Packit Service |
3880ab |
arraystats(x, limit, &mu, &sigma, &rho);
|
|
Packit Service |
3880ab |
#ifdef DEBUG
|
|
Packit Service |
3880ab |
fprintf(stderr, "%d values, mu %10.4f, sigma %10.4f, rho %10.4f\n",
|
|
Packit Service |
3880ab |
limit, mu, sigma, rho);
|
|
Packit Service |
3880ab |
#endif
|
|
Packit Service |
3880ab |
|
|
Packit Service |
3880ab |
table = makedist(x, limit, mu, sigma);
|
|
Packit Service |
3880ab |
free((void *) x);
|
|
Packit Service |
3880ab |
cumulativedist(table, DISTTABLESIZE, &total);
|
|
Packit Service |
3880ab |
inverse = inverttable(table, TABLESIZE, DISTTABLESIZE, total);
|
|
Packit Service |
3880ab |
interpolatetable(inverse, TABLESIZE);
|
|
Packit Service |
3880ab |
printtable(inverse, TABLESIZE);
|
|
Packit Service |
3880ab |
return 0;
|
|
Packit Service |
3880ab |
}
|