diff options
Diffstat (limited to 'netem/maketable.c')
| -rw-r--r-- | netem/maketable.c | 235 |
1 files changed, 235 insertions, 0 deletions
diff --git a/netem/maketable.c b/netem/maketable.c new file mode 100644 index 0000000..ad8620a --- /dev/null +++ b/netem/maketable.c @@ -0,0 +1,235 @@ +/* SPDX-License-Identifier: NIST-PD */ +/* + * 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; +} |