#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>	/* For stat()  */
/* #include <float.h>
#include <limits.h>  */
#include "blockstat.h"

struct IMGBM_OLDSTYLE {
	int	index;
	int	z;
};

struct WEIGHT_WORK {
	double	p1;
	double	p2;
};

int	main (int argc, char **argv) {

	char	*fname_in = NULL;	/* Input file name */
	char	*fname_out = NULL;	/* Output file name */
	
	struct BM_DATA_IN *bmdi;
	struct WEIGHT_WORK *wtw;
	
	FILE	*fp;

	size_t npoints = 0;		/* Number of bmdi struct points incoming  */
	size_t	nwork = 0;		/* Number of bmdi points in most full box  */
	
	int	nerrs = 0;
	int	old_style = 0;
	int	sid_out = 0;
	int	karg;

	int compare_bmdi (struct BM_DATA_IN *p1, struct BM_DATA_IN *p2);
/*	int (*compare_bmdi) ();  */
	void	make_old_output (struct BM_DATA_IN *bmdi, size_t npoints, FILE *fp);
	void	make_sid_output (struct BM_DATA_IN *bmdi, size_t npoints, FILE *fp);
	void	make_bmdo_output (struct BM_DATA_IN *bmdi, struct WEIGHT_WORK *wtw, size_t npoints, FILE *fp);
	size_t	check_file_size (char *fname);
	size_t	n_in_most_full_box (struct BM_DATA_IN *bmdi, size_t npoints);

	for (karg = 1; karg < argc; karg++) {
		if (argv[karg][0] == '-') {
			switch (argv[karg][1]) {
			
				case 'B':
					old_style = 1;
					break;
			
				case 'I':
					fname_in = &argv[karg][2];
					if (strlen(fname_in) < 1) nerrs++;
					break;
			
				case 'O':
					fname_out = &argv[karg][2];
					if (strlen(fname_out) < 1) nerrs++;
					break;
			
				case 'S':
					sid_out = 1;
					break;
			
				default:
					nerrs++;
					break;
			}
		}
		else {
			nerrs++;
		}
	}
	
	if (nerrs || (argc == 1) ) {
	
		fprintf (stderr, "usage:  bmprep2bm -I<inputfile> -O<outputfile> [-B]\n\n");
		if (argc == 1) exit (EXIT_SUCCESS);
		fprintf (stderr, "\t-I<inputfile> gives path to input holding concatenation of all binary cm2bm_prep outputs.\n");
		fprintf (stderr, "\t-O<outputfile> gives path to output holding binary BM_DATA_OUT structures.\n");
		fprintf (stderr, "\t-B option creates old-style 'imgbm' with index and depth.\n");
		fprintf (stderr, "\t-S option creates old-style 'imgbm' with index and source_id.\n");
		exit (EXIT_FAILURE);
	}

	npoints = check_file_size (fname_in);
	
	if ( (bmdi = (struct BM_DATA_IN *)malloc((size_t)(npoints * sizeof(struct BM_DATA_IN)))) == NULL)  {
		fprintf (stderr, "bmprep2bm:  Failed to malloc %lg bytes to hold %lg input points.\n", 
			(double)npoints * sizeof(struct BM_DATA_IN), (double)npoints);
		exit (EXIT_FAILURE);
	}
	
	if ( (fp = fopen (fname_in, "r") ) == NULL) {
		fprintf (stderr, "bmprep2bm:  Failed to open file %s\n", fname_in);
		exit (EXIT_FAILURE);
	}
	
	if ( (fread ((void *)bmdi, sizeof(struct BM_DATA_IN), npoints, fp)) != npoints) {
		fprintf (stderr, "bmprep2bm:  Failed to read file %s\n", fname_in);
		fclose (fp);
		exit (EXIT_FAILURE);
	}
	fclose (fp);
	
		
	if ( (fp = fopen (fname_out, "w") ) == NULL) {
		fprintf (stderr, "bmprep2bm:  Failed to open file %s\n", fname_in);
		exit (EXIT_FAILURE);
	}

	fprintf (stderr, "bmprep2bm:  Read completed.  Sorting...");
	
	qsort ((void *)bmdi, npoints, sizeof(struct BM_DATA_IN), compare_bmdi);
	
	fprintf (stderr, "Done.  Creating output.\n");
	
        if(old_style == 1 || sid_out == 1) {
		if (old_style) {
			make_old_output (bmdi, npoints, fp);
		}
	
		if (sid_out) {
			make_sid_output (bmdi, npoints, fp);
		}
	}
	else {
		nwork = n_in_most_full_box (bmdi, npoints);
		if ( (wtw = (struct WEIGHT_WORK *)malloc((size_t)(nwork	* sizeof(struct WEIGHT_WORK)))) == NULL)  {
			fprintf (stderr, "bmprep2bm:  Failed to malloc %lg weight work points.\n", (double)nwork);
			exit (EXIT_FAILURE);
		}
		make_bmdo_output (bmdi, wtw, npoints, fp);
		free ((void *)wtw);
	}
	
	fclose (fp);
	free ((void *)bmdi);
	exit (0);
}

size_t	check_file_size (char *fname) {

	struct stat buf;
	double	p, q, t;
	
	if (stat (fname, &buf) ) {
		fprintf (stderr, "bmprep2bm:  Cannot stat filename %s\n", fname);
		exit (EXIT_FAILURE);
	}
	p = (double) buf.st_size;
	q = (double) sizeof (struct BM_DATA_IN);
	t = p/q;
	p = floor(t);
	q = ceil (t);
	if (p <= 0.0 || p != q) {
		fprintf (stderr, "bmprep2bm:  Bad file size for %s\n", fname);
		exit (EXIT_FAILURE);
	}
	return ((size_t)p);
}

int	compare_bmdi (struct BM_DATA_IN *p1, struct BM_DATA_IN *p2) {

	if ( (p1->index) < (p2->index) ) {
		return (-1);
	} else if ( (p1->index) > (p2->index) ) {
		return (1);
	} else {
		if ( (p1->z) < (p2->z) ) {
			return (-1);
		} else if ( (p1->z) > (p2->z) ) {
			return (1);
		} else {
			return (0);
		}
	}
}

void	make_old_output (struct BM_DATA_IN *bmdi, size_t npoints, FILE *fp) {

	size_t j, k, i;
	int	n;
	struct IMGBM_OLDSTYLE old_out;
	
	j = 0;
	while (j  < npoints) {
		k = j + 1;
		while (k < npoints && bmdi[k].index == bmdi[j].index) k++;
		n = k - j;
		i = j + (n/2);
		if (n%2) {
			old_out.z = bmdi[i].z;
		} else {
			old_out.z = (int)rint(0.5*(bmdi[i].z + bmdi[i-1].z));
		}
		old_out.index = bmdi[j].index;
   	fwrite ((void *)&old_out, sizeof(struct IMGBM_OLDSTYLE), (size_t)1, fp);
		j = k;
	}
}


void	make_sid_output (struct BM_DATA_IN *bmdi, size_t npoints, FILE *fp) {

	size_t j, k, i;
	int	n;
	struct IMGBM_OLDSTYLE old_out;
	
	j = 0;
	while (j  < npoints) {
		k = j + 1;
		while (k < npoints && bmdi[k].index == bmdi[j].index) k++;
		n = k - j;
		i = j + (n/2);
                old_out.z = bmdi[i].source_id;
		old_out.index = bmdi[j].index;
   	fwrite ((void *)&old_out, sizeof(struct IMGBM_OLDSTYLE), (size_t)1, fp);
		j = k;
	}
}

void	make_bmdo_output (struct BM_DATA_IN *bmdi, struct WEIGHT_WORK *wtw, size_t npoints, FILE *fp) {

	struct	BM_DATA_OUT	bmdo;
	size_t j, k;
	void	process_a_box (struct BM_DATA_IN *bmdi, struct WEIGHT_WORK *wtw, int n_in_box, struct BM_DATA_OUT *bmdo);

	j = 0;
	while (j  < npoints) {
		k = j + 1;
		while (k < npoints && bmdi[k].index == bmdi[j].index) k++;
		process_a_box (&bmdi[j], wtw, (int)(k-j), &bmdo);
		fwrite ((void *)&bmdo, sizeof(struct BM_DATA_OUT), (size_t)1, fp);
		j = k;
	}
}


size_t	n_in_most_full_box (struct BM_DATA_IN *bmdi, size_t npoints) {

	size_t j, k, n, m = 0;
	
	j = 0;
	while (j  < npoints) {
		k = j + 1;
		while (k < npoints && bmdi[k].index == bmdi[j].index) k++;
		n = k-j;
		if (n > m) m = n;
		j = k;
	}
	return (m);
}

void	process_a_box (struct BM_DATA_IN *bmdi, struct WEIGHT_WORK *wtw, int n, struct BM_DATA_OUT *bmdo) {

	double	sum1, sum2, p, r;
	int	k;
	short int smin;

	bmdo->z_min = bmdi[0].z;
	bmdo->id_min = bmdi[0].source_id;
	k = n - 1;
	bmdo->z_max = bmdi[k].z;
	bmdo->id_max = bmdi[k].source_id;
	
	if (n == 1) {
		bmdo->z_mean = bmdi[0].z;
		bmdo->z_std = bmdi[0].s;
		bmdo->z_median = bmdi[0].z;
		bmdo->z_spread = bmdi[0].s;
		bmdo->id_med = bmdi[0].source_id;
		return;
	}

	/* Set weight work parameters.  p1 is proportional to 1/sigma.
	p2 is proportional to 1/sigma^2.  Both are normalized so their sums are 1.  */
	sum1 = 0.0;
	sum2 = 0.0;
	for (k = 0; k < n; k++) {
		p = 1.0 / bmdi[k].s;
		wtw[k].p1 = p;
		wtw[k].p2 = p*p;
		sum1 += wtw[k].p1;
		sum2 += wtw[k].p2;
	}
	sum1 = 1.0 / sum1;
	sum2 = 1.0 / sum2;
	for (k = 0; k < n; k++) {
		wtw[k].p1 *= sum1;
		wtw[k].p2 *= sum2;
	}
	
	/* I think that maybe 1/sqrt(sum2) is now what we want for the s.d. around the mean,
	under the hypothesis that all the z in this box are drawn from a distribution
	having the same mean.
	But I don't have my notes here at SIO (7 Feb 2007) so I'll have to check this when
	I get back to Silver Spring.
	In any case, the true s.d. will be bigger due to variation of the mean across the box.  */
	
	sum1 = 0.0;
	sum2 = 0.0;
	for (k = 0; k < n; k++) {
		if (sum1 < 0.5 && sum1 + wtw[k].p1 > 0.5) {
			/* Then this value of k holds the weighted median  */
			bmdo->z_median = bmdi[k].z;
			bmdo->id_med = bmdi[k].source_id;
		}
		if (sum1 == 0.5) {
			/* Then the weighted median is either this one or the previous one  */
			if (wtw[k].p1 > wtw[k-1].p1) {
				/* Take this one  */
				bmdo->z_median = bmdi[k].z;
				bmdo->id_med = bmdi[k].source_id;
			} else if (wtw[k].p1 < wtw[k-1].p1) {
				/* Take the previous one */
				bmdo->z_median = bmdi[k-1].z;
				bmdo->id_med = bmdi[k-1].source_id;
			} else {
				/* The median should be defined as the average of these two.
					The source id is ambiguous in this case.  */
				bmdo->z_median = (short int) rint (0.5 * (bmdi[k].z + bmdi[k-1].z));
				bmdo->id_med = bmdi[k].source_id;	/* Take the shallower one (arbitrary choice)  */
			}
		}
		sum1 += wtw[k].p1;
		sum2 += wtw[k].p2 * bmdi[k].z;
	}
	bmdo->z_mean = (short int) rint (sum2);
	
	/* It isn't clear how to define the best estimate of the uncertainty.  
	For now I will try this:  */
	
	smin = 32000;
	sum1 = 0.0;
	sum2 = 0.0;
	for (k = 0; k < n; k++) {
		if (bmdi[k].s < smin) smin = bmdi[k].s;
		sum1 += (bmdi[k].z - bmdo->z_median) * wtw[k].p1;
		r = bmdi[k].z - bmdo->z_mean;
		sum2 += r*r * wtw[k].p2;
	}
	bmdo->z_spread = (short int) ceil (sum1);	/* This is not an IQR but instead a Mean Abs Dev  */
	if (bmdo->z_spread < smin) bmdo->z_spread = smin;
	bmdo->z_std = (short int) ceil (sqrt(sum2));
	if (bmdo->z_std < smin) bmdo->z_std = smin;
}