/*
 * 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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>


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<limit; ++i){
                if (fscanf(fp, "%lf", &x[i]) != 1 ||
                    feof(fp))
                        break;
                ++n;
        }
        *number = n;
        return x;
}

void
arraystats(double *x, int limit, double *mu, double *sigma, double *rho)
{
        int n=0, i;
        double sumsquare=0.0, sum=0.0, top=0.0;
        double sigma2=0.0;

        for (i=0; i<limit; ++i){
                sumsquare += x[i]*x[i];
                sum += x[i];
                ++n;
        }
        *mu = sum/(double)n;
        *sigma = sqrt((sumsquare - (double)n*(*mu)*(*mu))/(double)(n-1));

        for (i=1; i < n; ++i){
                top += ((double)x[i]- *mu)*((double)x[i-1]- *mu);
                sigma2 += ((double)x[i-1] - *mu)*((double)x[i-1] - *mu);

        }
        *rho = top/sigma2;
}

/* Create a (normalized) distribution table from a set of observed
 * values.  The table is fixed to run from (as it happens) -4 to +4,
 * with granularity .00002.
 */

#define TABLESIZE       16384/4
#define TABLEFACTOR     8192
#ifndef MINSHORT
#define MINSHORT        -32768
#define MAXSHORT        32767
#endif

/* Since entries in the inverse are scaled by TABLEFACTOR, and can't be bigger
 * than MAXSHORT, we don't bother looking at a larger domain than this:
 */
#define DISTTABLEDOMAIN ((MAXSHORT/TABLEFACTOR)+1)
#define DISTTABLEGRANULARITY 50000
#define DISTTABLESIZE (DISTTABLEDOMAIN*DISTTABLEGRANULARITY*2)

static int *
makedist(double *x, int limit, double mu, double sigma)
{
        int *table;
        int i, index, first=DISTTABLESIZE, last=0;
        double input;

        table = calloc(DISTTABLESIZE, sizeof(int));
        if (!table) {
                perror("table alloc");
                exit(3);
        }

        for (i=0; i < limit; ++i) {
                /* Normalize value */
                input = (x[i]-mu)/sigma;

                index = (int)rint((input+DISTTABLEDOMAIN)*DISTTABLEGRANULARITY);
                if (index < 0) index = 0;
                if (index >= 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;
}