/*
  gbbin  (ver. 5.6)  -- A program to bin data
  Copyright (C) 1998-2018 Giulio Bottazzi

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  (version 2) as published by the Free Software Foundation;
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

#include "tools.h"

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

  double **data=NULL; /* array of values */
  size_t rows=0,columns=0;

  size_t steps=10;
  size_t binum;
  int numpres;

  char *splitstring = strdup(" \t");

  double wmin,wmax;


  /* otpions */
  char * const output_opts[] = {"xmean","xstd","xmin","xmax","xmedian","ymean","yadev","ystd","yvar","yskew","ykurt","ymin","ymax","ymedian","num",NULL};

  size_t outnum;
  size_t *outtype;

  int o_verbose=0;
  int  o_printbinned = 0;
  int o_window=0; /* if the windows to bin is provided */
  

  /* variable to handle data splitting */
  size_t printnum = 0;
  size_t *printcolumn=NULL;

  /* variables for reading command line options */
  int opt;

  /* set default output */

  outnum=2;
  outtype = (size_t *) my_alloc(2*sizeof(size_t));
  outtype[0]=0;
  outtype[1]=5;

  /* COMMAND LINE PROCESSING */
  
  /* read the command line */    
  while((opt=getopt_long(argc,argv,"vn:hF:O:xyc:w:",gb_long_options, &gb_option_index))!=EOF){
    if(opt==0){
      gbutils_header(argv[0],stdout);
      exit(0);
    }
    else if(opt=='?'){
      fprintf(stderr,"option %c not recognized\n",optopt);
      return(-1);
    }
    else if(opt=='F'){
      /*set the fields separator string*/
      free(splitstring);
      splitstring = strdup(optarg);
    }
    else if(opt=='n'){
      /*set the number of steps*/
      steps = (size_t) atoi(optarg);
    }
    else if(opt=='v'){
      /*increase verbosity*/
      o_verbose = 1;
    }
    else if(opt=='c'){
      /* set the columns whose elements are required */
      char *stmp1=strdup (optarg);
      char *stmp2;
      char *stmp3=stmp1;

      o_printbinned = 1;

      while( (stmp2=strsep (&stmp1,",")) != NULL && strlen(stmp2)>0 ){
	printnum++;
	printcolumn = (size_t *) my_realloc((void *) printcolumn,printnum*sizeof(size_t));
	printcolumn[printnum-1]= (size_t) atoi(stmp2);
      }

      free(stmp3);
    }
    else if(opt=='x'){
      /* print binned x values */
      o_printbinned = 1;

      printnum++;
      printcolumn = (size_t *) my_realloc((void *) printcolumn,printnum*sizeof(size_t));
      printcolumn[printnum-1]=1;
    }
    else if(opt=='y'){
      /* print binned y values */
      o_printbinned = 1;

      printnum++;
      printcolumn = (size_t *) my_realloc((void *) printcolumn,printnum*sizeof(size_t));
      printcolumn[printnum-1]=2;
    }
    else if(opt=='O'){
      /* set the output type */
      char *subopts,*subvalue=NULL;
      /* clear default settings */
      outnum=0;
      free(outtype); outtype=NULL;

      /* define new output */
      subopts = optarg;
      while (*subopts != '\0'){
	int  itmp1 =  getsubopt(&subopts, output_opts, &subvalue);
	outnum++;
	outtype = (size_t *) my_realloc((void *)outtype,outnum*sizeof(size_t));

	if(itmp1 == -1){/* Unknown suboption. */
	  fprintf (stderr,"ERROR (%s): Unknown output type `%s'\n", 
		   GB_PROGNAME,subvalue);
	  exit(1);
	}
	else
	  outtype[outnum-1]=(size_t) itmp1;
      }
    }
    else if(opt=='w'){
      /*set the binning window */

      char *stmp1=strdup (optarg);
      char *stmp2;
      char *stmp3=stmp1;

      stmp2=strsep (&stmp1,",");
      if(strlen(stmp2)>0 && stmp1 != NULL && strlen(stmp1)>0){
	wmin=atof(stmp2);	
	wmax=atof(stmp1);
      }
      else
	fprintf(stderr,"ERROR (%s): provide comma separated min and max. Try option -h.\n",GB_PROGNAME);

      free(stmp3);

      o_window=1;
    }
    else if(opt=='h'){
      /*print short help*/
      fprintf(stdout,"Compute binned statistics. Data are read from standard input as records    \n");
      fprintf(stdout,"(X,Y1,Y2,...). Equipopulated bins are built with respect to the first field\n");
      fprintf(stdout,"X. Option -O decide which statistics are printed for each bin. With options\n");
      fprintf(stdout,"-x, -y or -c elements in different bins are split in different columns.    \n");
      fprintf(stdout,"\nUsage: %s [options]\n\n",argv[0]);
      fprintf(stdout,"Options:\n");
      fprintf(stdout," -n  set the number of equipopulated bins (default 10)                     \n");
      fprintf(stdout," -w  min,max set manually the binning window. Ignored with -x,-y and -c    \n");
      fprintf(stdout," -O  set the output with a comma separated list of variables: xmean, xmin, \n");
      fprintf(stdout,"     xmax, xstd, xmedian, ymean, yadev, ystd, yvar, yskew, ykurt, ymin,    \n");
      fprintf(stdout,"     ymax, ymedian, num (default xmean, ymean)                             \n");
      fprintf(stdout," -c#  the values in column # are split in different columns; # can be a comma\n");
      fprintf(stdout,"      separated list of columns                           \n");
      fprintf(stdout," -x  equivalent to -c 1                                                    \n");
      fprintf(stdout," -y  equivalent to -c 2                                                    \n");
      fprintf(stdout," -F  specify the input fields separators (default \" \\t\")                \n");
      fprintf(stdout," -v  verbose mode\n");
      fprintf(stdout,"Examples:\n");
      fprintf(stdout," gbbin -n 20 < file  split the records (line) in 20 bins according to first \n");
      fprintf(stdout,"                     field. Print the average value of the bin entries for \n");
      fprintf(stdout,"                     each column. If 'file' has 3 columns, the output has  \n");
      fprintf(stdout,"                     twenty rows and three columns.\n");
      fprintf(stdout," gbbin -O ymedian < file  bin the records with respect to the first value and\n");
      fprintf(stdout,"                          print the median value of the other columns in each bin\n");
      fprintf(stdout," gbbin -c 3 < file   the binned values of the third columns are printed on \n");
      fprintf(stdout,"                     separate columns. The output has ten columns.\n");
      exit(0);
    }
  }
  /* END OF COMMAND LINE PROCESSING */

  /* initialize global variables */
  initialize_program(argv[0]);

  /* load data: data[column][row] */
  loadtable(&data,&rows,&columns,0,splitstring);

  if(columns==2)
    sort2(data[0],data[1],rows,2);
  else if(columns>2)
    sortn(data,columns,rows,0,2);
  else{
    fprintf(stdout,
	    "ERROR (%s): Provide at least two columns of input\n",GB_PROGNAME);
    exit(1);
  }
    

  /* compute the size of the bins */
  binum=rows/steps;   /* number of points per bin  */
  numpres=rows%steps; /* number of residual points */

  /* ++++++++++++++++++++++++++++ */
  if(o_verbose == 1){
    size_t i;
    fprintf(stdout,"#read rows       %zd\n",rows);
    fprintf(stdout,"#read columns    %zd\n",columns);
    fprintf(stdout,"#number of bins  %zd\n",steps);
    fprintf(stdout,"#x support       [%f,%f]\n",data[0][0],data[0][rows-1]);
    fprintf(stdout,"#Output: ");

    if(outnum==0){
      fprintf(stdout,"%s %s\n",output_opts[0],output_opts[4]);
    }
    else{
      for(i=0;i<outnum;i++)
	fprintf(stdout,"%s ",output_opts[outtype[i]]);
      fprintf(stdout,"\n");
    }
  }
  /* ++++++++++++++++++++++++++++ */



  if(o_printbinned == 0){ /* option -O */

    size_t i,h,j;
    double **stats;

    size_t startindex=0,actualbinum;

    /* set to 1 if the median has to be computed */
    char o_xmedian,o_ymedian=0;

    /* allocate space for the results of stat */
    stats = (double **) my_alloc(columns*sizeof(double *));
    for(i=0;i<columns;i++)
      stats[i] = (double *) my_alloc(9*sizeof(double));

    /* check if median for x and ys should be computed */
    for(i=0;i<outnum;i++)
      if(outtype[i]==4){
	o_xmedian=1;
	break;
      }

    for(i=0;i<outnum;i++)
      if(outtype[i]==13){
	o_ymedian=1;
	break;
      }
    
    for(i=0;i<steps;i++){

      /* chose the number of observations */
      if(o_window == 1){/* windows width preset */

	const double xmax= wmin+(wmax-wmin)*(i+1)/steps;

	/* initially remove too small observations */
	if(i==0){
	  for(j=startindex;j<rows;j++){
	    if(data[0][j] >= wmin)
	      break;
	  }
	  if(j>startindex) startindex=j;
	}
	
	/* chose the limit index of the observations in the windows */
	actualbinum=0;
	for(j=startindex;j<rows;j++){
	  if(data[0][j] <= xmax)
	    actualbinum++;
	  else
	    break;
	}
      }
      else{/* equipopulated bins */
	if(i==0){
	  actualbinum=binum+numpres/2;
	}
	else if(i==steps-1){
	  actualbinum=binum+numpres-numpres/2;
	}
	else{
	  actualbinum=binum;
	}
      }

      /* compute stats for x values (first column) */
      moment_short(data[0]+startindex,actualbinum,stats[0],stats[0]+1,stats[0]+2,stats[0]+3);
      if(o_xmedian==1){		/* need to compute the median */
	double* vals = (double *) my_alloc(actualbinum*sizeof(double));
	for (j=0;j<actualbinum;j++)
	  vals[j]=data[0][startindex+j];
	qsort(vals,actualbinum,sizeof(double),sort_by_value);
	stats[0][8] = (actualbinum%2 ==0 ? .5*(vals[actualbinum/2-1]+vals[actualbinum/2]) : vals[(actualbinum-1)/2]);
	free(vals);
      }


      /* compute stats for y values (from the second columns) */
      for(h=1;h<columns;h++){
	moment(data[h]+startindex,actualbinum,
	       stats[h],stats[h]+1,stats[h]+2,stats[h]+3,
	       stats[h]+4,stats[h]+5,stats[h]+6,stats[h]+7);

	if(o_ymedian==1){		/* need to compute the median */
	  double* vals = (double *) my_alloc(actualbinum*sizeof(double));
	  for (j=0;j<actualbinum;j++)
	    vals[j]=data[h][startindex+j];
	  qsort(vals,actualbinum,sizeof(double),sort_by_value);
	  stats[h][8] = (actualbinum%2 ==0 ? .5*(vals[actualbinum/2-1]+vals[actualbinum/2]) : vals[(actualbinum-1)/2]);
	  free(vals);
	}
      }


      /*

      output options

      0  "xmean"
      1  "xstd"
      2  "xmin"
      3  "xmax"
      4  "xmedian" 
      5  "ymean"
      6  "yadev"
      7  "ystd"
      8  "yvar"
      9  "yskew"
      10  "ykurt"
      11 "ymin"
      12 "ymax"
      13 "ymedian"
      14 "num"

      statistics in array 'stat'

      0 "mean"
      1 "adev"
      2 "std"
      3 "var"
      4 "skew"
      5 "kurt"
      6 "min"
      7 "max"
      8 "median"


      returned by moments

      0 ave
      1 adev
      2 sdev
      3 var
      4 skew
      5 curt
      6 min
      7 max

      returned by moments_short

      0 ave
      1 sdev
      2 min
      3 max
      
      */

      /* output */
      for(j=0;j<outnum-1;j++)
	if(outtype[h]==14)         /* number of observations */
	  printf(FLOAT_SEP, (double) actualbinum );
	else if(outtype[j]<4)      /* x values stat */
	  printf(FLOAT_SEP,stats[0][outtype[j]]);
	else if(outtype[j]==4)      /* x median */
	  printf(FLOAT_SEP,stats[0][8]);
	else
	  for(h=1;h<columns;h++)   /* y values stat */
	    printf(FLOAT_SEP,stats[h][outtype[j]-5]);
      if(outtype[outnum-1]==14)         /* number of observations */
	printf(FLOAT_NL, (double) actualbinum);
      else if(outtype[outnum-1]<4)      /* x values stat */
	printf(FLOAT_NL,stats[0][outtype[j]]);
      else if(outtype[j]==4)      /* x median */
	printf(FLOAT_SEP,stats[0][8]);
      else{
	for(h=1;h<columns-1;h++)   /* y values stat */
	  printf(FLOAT_SEP,stats[h][outtype[outnum-1]-5]);
	printf(FLOAT_NL,stats[columns-1][outtype[outnum-1]-5]);
      }

      startindex+=actualbinum;
    }

    /* deallocate stat */
    for(i=0;i<columns;i++) free(stats[i]);
    free(stats);

  }
  else{/* options -c, -x or -y */

    size_t i,j,h;

    /* check that the column spec are OK */
    for(h=0;h<printnum;h++)
      if(printcolumn[h]>columns){
	fprintf(stdout,
		"ERROR (%s): Col spec %zd larger than number of col\n",
		GB_PROGNAME,printcolumn[h]);
	exit(1);
      }

    for(j=0;j<binum;j++){
      for(i=0;i<steps;i++)
	for(h=0;h<printnum;h++)
	  printf(FLOAT_SEP,data[printcolumn[h]-1][j+i*binum]);
      printf("\n");
    }

    for(j=binum;j<binum+numpres/2;j++){
      for(h=0;h<printnum;h++)
	printf(FLOAT_SEP,data[printcolumn[h]-1][j]);

      for(i=0;i<steps-2;i++)
	for(h=0;h<printnum;h++)
	  printf(EMPTY_SEP,"nan");

      if(steps>1)
	for(h=0;h<printnum;h++)
	  printf(FLOAT_SEP,data[printcolumn[h]-1][j+(steps-1)*binum]);

      printf("\n");
    }

    for(j=binum+numpres/2;j<binum+numpres-numpres/2;j++){

      for(i=0;i<steps-1;i++)
	for(h=0;h<printnum;h++)
	  printf(EMPTY_SEP,"nan");

      for(h=0;h<printnum;h++)
	printf(FLOAT_SEP,data[printcolumn[h]-1][j+(steps-1)*binum]);
      
      printf("\n");
    }

  }
 
  return(0);
}