Blame netem/maketable.c

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
}