#include "misc.h"
#include "stats.h"

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

char *StringOfDashes(int n)
{
    char *s = Malloc(n+1);
    assert(n>=0);
    s[n] = '\0';
    while(n-- > 0)
        s[n] = '-';
    return s;
}
char USAGE[] = "USAGE: %s [-g|-gain] [-[h|p][c] numBins min max]\n"
    "read numbers on the standard input (no files!), output basic statistics\n"
    "g = geometric\n"
    "gain = geometric + total product\n"
    "h = do a histogram; p = normalize the histogram to probabilities per unit\n"
    "c = make the histogram cumulative\n";
/*
** Read real numbers, and print stats.
*/
int main(int argc, char *argv[])
{
    STAT *st;
    double sample, histMin, histMax;
    int numBins = 0, ss;
    Boolean histCumulative = false, histNormalize = false, geom=false, gain=false;

    if(argc > 1)
    {
	if(argc == 2 && strcmp(argv[1], "-g") == 0)
	    geom = true;
	else if(argc == 2 && strcmp(argv[1], "-gain") == 0)
	{
	    geom = true;
	    gain = true;
	}
	else if(argc == 5 &&
	    (strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "-hc") == 0 ||
	     strcmp(argv[1], "-p") == 0 || strcmp(argv[1], "-pc") == 0) &&
	    (numBins = atoi(argv[2])) > 0 &&
	    (histMin = atof(argv[3])) < (histMax = atof(argv[4]))
	)
	{
	    if(strcmp(argv[1], "-hc") == 0)
		histCumulative = true;
	    if(strcmp(argv[1], "-p") == 0)
		histNormalize = true;
	    if(strcmp(argv[1], "-pc") == 0)
		histCumulative = histNormalize = true;
	}
	else
	{
	    fprintf(stderr, USAGE, argv[0]);
	    exit(1);
	}

	st = StatAlloc(numBins, histMin, histMax, geom);
    }
    else
	st = StatAlloc(0, 0, 0, geom);

    while(scanf("%lf", &sample) == 1)
	StatAddSample(st, sample);

    if(gain)
	printf("# %d mean %.16g min %.16g max %.16g stdDev %.16g var %.16g gain %.16g\n",
	    StatSampleSize(st), StatGeomMean(st), StatMin(st), StatMax(st),
	    StatGeomStdDev(st), StatGeomVariance(st), exp(st->geomSum));
    else if(geom)
	printf("# %d mean %.16g min %.16g max %.16g stdDev %.16g var %.16g\n",
	    StatSampleSize(st), StatGeomMean(st), StatMin(st), StatMax(st),
	    StatGeomStdDev(st), StatGeomVariance(st));
    else
	printf("# %d mean %.16g min %.16g max %.16g stdDev %.16g var %.16g skew %.16g\n",
	    StatSampleSize(st), StatMean(st), StatMin(st), StatMax(st),
	    StatStdDev(st), StatVariance(st), StatSkew(st));

    if(numBins)
    {
	double binSize = st->histWidth/numBins;
	/*double norm = histNormalize ? StatSampleSize(st)*binSize : 1;*/

	if(histCumulative)
	    StatCumulativeHistogram(st);

	for(ss=-1; ss <= st->numHistBins; ss++)
	{
	    int ivalue = st->histogram[ss];
	    char *s = StringOfDashes(ivalue);
	    double value;

	    if(histNormalize)
	    {
		/* For example, say we have 1000 samples, and 12 bins
		 * of size 10 each, so the histogram spans 0 to 120.
		 * Then if the histogram is cumulative, then it doesn't
		 * matter what the binSize was to left of us, it just
		 * matters what they added up to.  So there's no
		 * binSize-dependent normalization required for
		 * cumulative.  For non-cumulative, however, if we
		 * are doing cumulative based on probablity rather than
		 * the integer histogram value, then the integral across
		 * bins must be 1.  So normalization requires us to
		 * additionally divide by the binSize.
		 */
		if(histCumulative)
		    value = ivalue/(double)StatSampleSize(st);
		else
		    value = ivalue/binSize/StatSampleSize(st);
	    }
	    else
		value = ivalue;
	    printf("%.16g %.16g\n", st->histMin + ss*binSize, value);
	    Free(s);
	}
    }

    StatFree(st);

    return 0;
}
