netem/maketable.c

   1 /*
   2  * Experimental data  distribution table generator
   3  * Taken from the uncopyrighted NISTnet code.
   4  *
   5  * Rread 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 *
  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
  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 *
  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
 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 *
 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
 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
 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
 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 }