1 /*
2 * Experimental data distribution table generator
3 * Taken from the uncopyrighted NISTnet code (public domain).
4 *
5 * Read in a series of "random" data values, either
6 * experimentally or generated from some probability distribution.
7 * From this, create the inverse distribution table used to approximate
8 * the distribution.
9 */
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <math.h>
13 #include <malloc.h>
14 #include <string.h>
15 #include <sys/types.h>
16 #include <sys/stat.h>
17
18
19 double *
readdoubles(FILE * fp,int * number)20 readdoubles(FILE *fp, int *number)
21 {
22 struct stat info;
23 double *x;
24 int limit;
25 int n=0, i;
26
27 fstat(fileno(fp), &info);
28 if (info.st_size > 0) {
29 limit = 2*info.st_size/sizeof(double); /* @@ approximate */
30 } else {
31 limit = 10000;
32 }
33
34 x = calloc(limit, sizeof(double));
35 if (!x) {
36 perror("double alloc");
37 exit(3);
38 }
39
40 for (i=0; i<limit; ++i){
41 fscanf(fp, "%lf", &x[i]);
42 if (feof(fp))
43 break;
44 ++n;
45 }
46 *number = n;
47 return x;
48 }
49
50 void
arraystats(double * x,int limit,double * mu,double * sigma,double * rho)51 arraystats(double *x, int limit, double *mu, double *sigma, double *rho)
52 {
53 int n=0, i;
54 double sumsquare=0.0, sum=0.0, top=0.0;
55 double sigma2=0.0;
56
57 for (i=0; i<limit; ++i){
58 sumsquare += x[i]*x[i];
59 sum += x[i];
60 ++n;
61 }
62 *mu = sum/(double)n;
63 *sigma = sqrt((sumsquare - (double)n*(*mu)*(*mu))/(double)(n-1));
64
65 for (i=1; i < n; ++i){
66 top += ((double)x[i]- *mu)*((double)x[i-1]- *mu);
67 sigma2 += ((double)x[i-1] - *mu)*((double)x[i-1] - *mu);
68
69 }
70 *rho = top/sigma2;
71 }
72
73 /* Create a (normalized) distribution table from a set of observed
74 * values. The table is fixed to run from (as it happens) -4 to +4,
75 * with granularity .00002.
76 */
77
78 #define TABLESIZE 16384/4
79 #define TABLEFACTOR 8192
80 #ifndef MINSHORT
81 #define MINSHORT -32768
82 #define MAXSHORT 32767
83 #endif
84
85 /* Since entries in the inverse are scaled by TABLEFACTOR, and can't be bigger
86 * than MAXSHORT, we don't bother looking at a larger domain than this:
87 */
88 #define DISTTABLEDOMAIN ((MAXSHORT/TABLEFACTOR)+1)
89 #define DISTTABLEGRANULARITY 50000
90 #define DISTTABLESIZE (DISTTABLEDOMAIN*DISTTABLEGRANULARITY*2)
91
92 static int *
makedist(double * x,int limit,double mu,double sigma)93 makedist(double *x, int limit, double mu, double sigma)
94 {
95 int *table;
96 int i, index, first=DISTTABLESIZE, last=0;
97 double input;
98
99 table = calloc(DISTTABLESIZE, sizeof(int));
100 if (!table) {
101 perror("table alloc");
102 exit(3);
103 }
104
105 for (i=0; i < limit; ++i) {
106 /* Normalize value */
107 input = (x[i]-mu)/sigma;
108
109 index = (int)rint((input+DISTTABLEDOMAIN)*DISTTABLEGRANULARITY);
110 if (index < 0) index = 0;
111 if (index >= DISTTABLESIZE) index = DISTTABLESIZE-1;
112 ++table[index];
113 if (index > last)
114 last = index +1;
115 if (index < first)
116 first = index;
117 }
118 return table;
119 }
120
121 /* replace an array by its cumulative distribution */
122 static void
cumulativedist(int * table,int limit,int * total)123 cumulativedist(int *table, int limit, int *total)
124 {
125 int accum=0;
126
127 while (--limit >= 0) {
128 accum += *table;
129 *table++ = accum;
130 }
131 *total = accum;
132 }
133
134 static short *
inverttable(int * table,int inversesize,int tablesize,int cumulative)135 inverttable(int *table, int inversesize, int tablesize, int cumulative)
136 {
137 int i, inverseindex, inversevalue;
138 short *inverse;
139 double findex, fvalue;
140
141 inverse = (short *)malloc(inversesize*sizeof(short));
142 for (i=0; i < inversesize; ++i) {
143 inverse[i] = MINSHORT;
144 }
145 for (i=0; i < tablesize; ++i) {
146 findex = ((double)i/(double)DISTTABLEGRANULARITY) - DISTTABLEDOMAIN;
147 fvalue = (double)table[i]/(double)cumulative;
148 inverseindex = (int)rint(fvalue*inversesize);
149 inversevalue = (int)rint(findex*TABLEFACTOR);
150 if (inversevalue <= MINSHORT) inversevalue = MINSHORT+1;
151 if (inversevalue > MAXSHORT) inversevalue = MAXSHORT;
152 inverse[inverseindex] = inversevalue;
153 }
154 return inverse;
155
156 }
157
158 /* Run simple linear interpolation over the table to fill in missing entries */
159 static void
interpolatetable(short * table,int limit)160 interpolatetable(short *table, int limit)
161 {
162 int i, j, last, lasti = -1;
163
164 last = MINSHORT;
165 for (i=0; i < limit; ++i) {
166 if (table[i] == MINSHORT) {
167 for (j=i; j < limit; ++j)
168 if (table[j] != MINSHORT)
169 break;
170 if (j < limit) {
171 table[i] = last + (i-lasti)*(table[j]-last)/(j-lasti);
172 } else {
173 table[i] = last + (i-lasti)*(MAXSHORT-last)/(limit-lasti);
174 }
175 } else {
176 last = table[i];
177 lasti = i;
178 }
179 }
180 }
181
182 static void
printtable(const short * table,int limit)183 printtable(const short *table, int limit)
184 {
185 int i;
186
187 printf("# This is the distribution table for the experimental distribution.\n");
188
189 for (i=0 ; i < limit; ++i) {
190 printf("%d%c", table[i],
191 (i % 8) == 7 ? '\n' : ' ');
192 }
193 }
194
195 int
main(int argc,char ** argv)196 main(int argc, char **argv)
197 {
198 FILE *fp;
199 double *x;
200 double mu, sigma, rho;
201 int limit;
202 int *table;
203 short *inverse;
204 int total;
205
206 if (argc > 1) {
207 if (!(fp = fopen(argv[1], "r"))) {
208 perror(argv[1]);
209 exit(1);
210 }
211 } else {
212 fp = stdin;
213 }
214 x = readdoubles(fp, &limit);
215 if (limit <= 0) {
216 fprintf(stderr, "Nothing much read!\n");
217 exit(2);
218 }
219 arraystats(x, limit, &mu, &sigma, &rho);
220 #ifdef DEBUG
221 fprintf(stderr, "%d values, mu %10.4f, sigma %10.4f, rho %10.4f\n",
222 limit, mu, sigma, rho);
223 #endif
224
225 table = makedist(x, limit, mu, sigma);
226 free((void *) x);
227 cumulativedist(table, DISTTABLESIZE, &total);
228 inverse = inverttable(table, TABLESIZE, DISTTABLESIZE, total);
229 interpolatetable(inverse, TABLESIZE);
230 printtable(inverse, TABLESIZE);
231 return 0;
232 }
233