/*
  gbhisto (ver. 5.6.1) -- Produce histogram from 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"


void printhist(double const *data,int const size,
	       int const steps,int const o_output,
	       int const o_bintype,double ratio,
	       double const refpoint,double xmin, double xmax) {

  int i=0;
  double stepsize;
  int *histo;
  double *midpoint,*width;
  int index=0;

  histo = (int *) malloc(steps*sizeof(int));
  midpoint = (double *) malloc(steps*sizeof(double));
  width = (double *) malloc(steps*sizeof(double));

  for(i=0;i<steps;i++){
    histo[i]=0;
  }

  if(o_bintype == 0){		/* linear binning */    
    stepsize=(xmax-xmin)/steps;
    
    for(i=0;i<steps;i++)
      midpoint[i] = xmin+stepsize*(i+refpoint);

    for(i=0;i<steps;i++)
      width[i] = stepsize;
    
    for(i=0;i<size;i++){
      const double val = data[i];
      if(val<xmax){
	index = (int) floor((val-xmin)/stepsize);
	index = (index<0?0:index);
	histo[index]++;
      }
      else if(val==xmax){
	histo[steps-1]++;
      }
    }
  }
  else if(o_bintype==1){	/* log binning */
    int ignored=0;
    double lxmin;

    if(xmax <=0) {
      fprintf(stderr,"ERROR (%s): no positive values! Exiting...\n",GB_PROGNAME);
      exit(-1);
    }
    else if(xmin <=0){      
      /* determine the smallest positive value */
      xmin=xmax;
      for(i=0;i<size;i++){
	if( data[i]>0 && xmin > data[i]) xmin = data[i]; 
      }
    }

    lxmin=log(xmin);

    stepsize=(log(xmax)-log(xmin))/steps;
    
    for(i=0;i<steps;i++)
      midpoint[i] = xmin*exp(stepsize*(i+refpoint));

    for(i=0;i<steps;i++)
      width[i] = xmin*( exp(stepsize*(i+1))-exp(stepsize*i) );
    
    for(i=0;i<size;i++){
      const double val = data[i];
      if(val>=xmin && val<xmax){
	index = (int) floor((log(val)-lxmin)/stepsize);
	index = (index<0?0:index);
	histo[index]++;
      }
      else if(val==xmax){
	histo[steps-1]++;
      }
      else if (val<xmin){
	ignored++;
      }
    }

    if(ignored>0)
      fprintf(stderr,"WARNING (%s): %d non positive values ignored\n",GB_PROGNAME,ignored);

  }
  else if(o_bintype==2){	/* geometric binning */

    double lratio = log(ratio);

    stepsize=(xmax-xmin)/(exp(lratio*steps)-1.);

    for(i=0;i<steps;i++){
      midpoint[i] = xmin+(exp(lratio*(refpoint+(double) i))-1.)*stepsize;
      width[i] = exp(lratio*i)*(ratio-1.)*stepsize;
    }
    
    for(i=0;i<size;i++){
      const double val = data[i];
      if(val<xmax){
	index = (int) floor( log(1.+(val-xmin)/stepsize)/lratio );
	index = (index<0?0:index);
	histo[index]++;
      }
      else if(val==xmax){
	histo[steps-1]++;
      }
    }
  }
  
  if(o_output==2){
    for(i=0;i<steps;i++){
      printf(FLOAT_SEP,midpoint[i]);
      printf(FLOAT_NL,(double) histo[i]/(size*width[i]));
    }
  }
  else if(o_output==1){
    for(i=0;i<steps;i++){
      printf(FLOAT_SEP,midpoint[i]);
      printf(FLOAT_NL,(double) histo[i]/size);
    }
  }
  else if(o_output==5){
    double cumprob=0;
    for(i=0;i<steps;i++){
      cumprob += (double) histo[i]/size;
      printf(FLOAT_SEP,midpoint[i]);
      printf(FLOAT_NL,cumprob);
    }
  }
  else if(o_output==6){
    double cumprob=0;
    for(i=0;i<steps;i++){
      cumprob += (double) histo[i]/size;
      printf(FLOAT_SEP,midpoint[i]);
      printf(FLOAT_NL,1-cumprob);
    }
  }
  else {
    for(i=0;i<steps;i++){
      printf(FLOAT_SEP,midpoint[i]);
      printf(INT_NL,histo[i]);
    }
  }
  
}


void printhist_discrete(double *data,int size,int const o_output) {


  qsort(data,size,sizeof(double),sort_by_value);

  int i=0,count=1;

  if(o_output==3){
  
    while(i<size){
      if(i<size-1 && data[i+1]==data[i]) count++;
      else{
	if(floor(data[i])==data[i])
	  printf(INT_SEP,(int) data[i]);
	else 
	  printf(FLOAT_SEP,data[i]);
	printf(INT_NL,count);
	count=1;
      }
      i++;
    }

  }
  else if(o_output==4){

    while(i<size){
      if(i<size-1 && data[i+1]==data[i]) count++;
      else{
	if(floor(data[i])==data[i])
	  printf(INT_SEP,(int) data[i]);
	else 
	  printf(FLOAT_SEP,data[i]);
	printf(FLOAT_NL,(double) count/size);
	count=1;
      }
      i++;
    }

  }

}




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

  unsigned steps=10;    
  int o_output=0;/*output type*/
  int o_bintype=0;/*binning type*/
  double refpoint=0.5;/* bin reference point  */
  int o_table=0;/*if input is of tabular type*/
  int o_window=0; /* if the windows to bin is provided */

  double wmin,wmax;

  double ratio=1.618; /*golden ratio by default... meaningless but nice*/
  
  int o_verbose=0;

  char *splitstring = strdup(" \t");
  
  /* COMMAND LINE PROCESSING */
  
  /* variables for reading command line options */
  /* ------------------------------------------ */
  int opt;
  /* ------------------------------------------ */
  
  /* read the command line */    
  while((opt=getopt_long(argc,argv,"n:M:tvhB:r:F:w:S:",gb_long_options, &gb_option_index))!=EOF){
    if(opt==0){
      gbutils_header(argv[0],stdout);
      exit(0);
    }
    else if(opt=='?'){
      fprintf(stderr,"ERROR (%s): option %c not recognized\n",GB_PROGNAME,optopt);
      return(-1);
    }
    else if(opt=='F'){
      /*set the fields separator string*/
      free(splitstring);
      splitstring = strdup(optarg);
    }
    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=='n'){
      /*set the number of steps*/
      steps = (unsigned) atoi(optarg);
    }
    else if(opt=='M'){
	/*set the type of output*/
      o_output = atoi(optarg);
      if(o_output<0 || o_output>6){
	fprintf(stderr,"ERROR (%s): output option '%d' not recognized. Try option -h.\n",
		GB_PROGNAME,o_output);
	exit(-1);
      } 
    }
    else if(opt=='B'){
      /*set the type of binning*/
      o_bintype = atoi(optarg);
      if(o_bintype <0 || o_bintype>2 ) {
	fprintf(stderr,"ERROR (%s): binning option '%d' not recognized. Try option -h.\n",
		GB_PROGNAME,o_bintype);
	exit(-1);
      }
    }
    else if(opt=='S'){
      /*set the type of binning*/
      refpoint = atof(optarg);
      if(refpoint <0 || refpoint >1 ) {
	fprintf(stderr,"ERROR (%s): unsupported reference point value '%f'. Try option -h.\n",
		GB_PROGNAME,refpoint);
	exit(-1);
      }
    }
    else if(opt=='r'){
	/*set the ration in the geometric binning*/
      ratio = atof(optarg);
      if(ratio <=1) {
	fprintf(stderr,"ERROR (%s): option -r requires a value larger than one. Try option -h.\n",
		GB_PROGNAME);
	exit(-1);
      }
    }
    else if(opt=='t'){
      /*set the table form*/
      o_table=1;
    }
    else if(opt=='v'){
      /*set verbose output*/
      o_verbose=1;
    }
    else if(opt=='h'){
      /*print short help*/
      fprintf(stdout,"Compute histogram counting the occurrences of each value (discrete data) or\n");
      fprintf(stdout,"binning the values on a regular grid (continuous data). Data are read from \n");
      fprintf(stdout,"standard input. For discrete data, values and their occurrences are printed;\n");
      fprintf(stdout,"For binned data, the bin reference point and the relative statistics are\n");
      fprintf(stdout,"printed. The type of variable and the statistics can be set with option -M.\n");
      fprintf(stdout,"The binning grid can be linear, logarithmic or geometric. The bin reference\n");
      fprintf(stdout,"point is the algebraic midpoint for linear binning and the geometric midpoint\n");
      fprintf(stdout,"for log or geometric binning. The reference point can be moved using option\n");
      fprintf(stdout,"-S. Set it to 0 for the lower bound or to 1 for the upper bound.\n");
      fprintf(stdout,"\nUsage: %s [options]\n\n",argv[0]);
      fprintf(stdout,"Options:\n");
      fprintf(stdout," -n  number of equispaced bins where the histogram is computed (default 10)\n");
      fprintf(stdout," -w  min,max set manually the binning window                               \n");
      fprintf(stdout," -M  set the type of output (default 0)                                    \n");
      fprintf(stdout,"      0  occurrences number (continuous binned variable)\n");
      fprintf(stdout,"      1  relative frequency (continuous binned variable)\n");
      fprintf(stdout,"      2  empirical density  (continuous binned variable)\n");
      fprintf(stdout,"      3  occurrences number (discrete variable)\n");
      fprintf(stdout,"      4  relative frequency (discrete variable)\n");
      fprintf(stdout,"      5  CDF (left cumulated probability; continuous binned variable)\n");
      fprintf(stdout,"      6  inverse CDF (1-CDF; continuous binned variable)\n");
      fprintf(stdout," -B  set the type of binning (default 0)                                   \n");
      fprintf(stdout,"      0  linear\n");
      fprintf(stdout,"      1  logarithmic\n");
      fprintf(stdout,"      2  geometric\n");
      fprintf(stdout," -S  set the bin reference point (default 0.5)                              \n");
      fprintf(stdout," -r  set the ratio for geometrical binning  (default 1.618)               \n");
      fprintf(stdout," -t  print the histogram for each imput column\n");
      fprintf(stdout," -F  specify the input fields separators (default \" \\t\")                \n");
      fprintf(stdout," -v  verbose mode\n");
      exit(0);
    }
  }
  /* END OF COMMAND LINE PROCESSING */

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

  /* ++++++++++++++++++++++++++++ */
  if(o_verbose){
    fprintf(stdout,"#variable type:\t");
    switch(o_output){
    case 0:
    case 1:
    case 2:
      fprintf(stdout,"continuous");
      break;
    case 3:
      fprintf(stdout,"discrete");
      break;
    }
    fprintf(stdout,"\n");

    fprintf(stdout,"#output type:\t");
    switch(o_output){
    case 0:
      fprintf(stdout,"occurencies number");
      break;
    case 1:
      fprintf(stdout,"relative frequency");
      break;
    case 2:
      fprintf(stdout,"empirical density");
      break;
    case 3:
      fprintf(stdout,"occurencies number");
      break;
    case 4:
      fprintf(stdout,"relative frequency");
      break;
    case 5:
      fprintf(stdout,"left cumulated probability");
      break;
    case 6:
      fprintf(stdout,"right cumulated probability");
      break;
    }
    fprintf(stdout,"\n");

    fprintf(stdout,"#binning:\t");
    switch(o_bintype){
    case 0:
      fprintf(stdout,"linear");
    case 1:
      fprintf(stdout,"logarithmic");
    case 2:
      fprintf(stdout,"geometric with ratio %f",ratio);
      break;
    }
    fprintf(stdout,"\n");

    fprintf(stdout,"#reference point: %f\n",refpoint);

    fprintf(stdout,"#bins number:\t");
    switch(o_output){
    case 0:
    case 1:
    case 2:
      fprintf(stdout,"%d",steps);
      break;
    }
    fprintf(stdout,"\n");

    fprintf(stdout,"#tabular data:\t");
    switch(o_table){
    case 0:
      fprintf(stdout,"off");
      break;
    case 1:
      fprintf(stdout,"on");
      break;
    }
    fprintf(stdout,"\n");

    fprintf(stdout,"#data statistics:\n");
    fprintf(stdout,
	    "# mean      stdev     skewness  kurtosis  ave.dev.  min       max\n");
  }  
  /* ++++++++++++++++++++++++++++ */


  
  if(o_table){
    size_t rows=0,columns=0,i;
    double **data=NULL;
    
   loadtable(&data,&rows,&columns,0,splitstring);
    

    /* ++++++++++++++++++++++++++++ */
    if(o_verbose){
      for(i=0;i<columns;i++){
	double xmean,xvar,xsd,xskew,xkurt,xadev,xmin,xmax;
	moment(data[i],rows,&xmean,&xadev,&xsd,&xvar,&xskew,&xkurt,&xmin,&xmax);
	fprintf(stdout,"# %+.2e %+.2e %+.2e %+.2e %+.2e %+.2e %+.2e %zd\n",
		xmean,xsd,xskew,xkurt,xadev,xmin,xmax,rows);
      }
    }
    /* ++++++++++++++++++++++++++++ */

    if(o_output>=3 && o_output<=4)
      for(i=0;i<columns;i++)
	printhist_discrete(data[i],rows,o_output);
    else      
      for(i=0;i<columns;i++){
	size_t j;
	double xmin = DBL_MAX;  /*set maximum double */
	double xmax = -DBL_MAX;  /*set to minimum double */
	
	/* determination of xmax and xmin */
	if(o_window){
	  xmin=wmin;
	  xmax=wmax;
	}
	else
	  for(j=0;j<rows;j++){
	    if( xmin > data[i][j]) xmin = data[i][j]; 
	    if( xmax < data[i][j]) xmax = data[i][j];
	  }
	/* ------------------------------ */

	printhist(data[i],rows,steps,o_output,o_bintype,ratio,refpoint,xmin,xmax);
	if(i<columns-1) printf("\n\n");
      }
   
  }
  else{
    double *data=NULL;
    size_t size;
    
   load(&data,&size,0,splitstring);
    
    /* ++++++++++++++++++++++++++++ */
    if(o_verbose){
      double xmean,xvar,xsd,xskew,xkurt,xadev,xmin,xmax;
      moment(data,size,&xmean,&xadev,&xsd,&xvar,&xskew,&xkurt,&xmin,&xmax);
      fprintf(stdout,"# %+.2e %+.2e %+.2e %+.2e %+.2e %+.2e %+.2e %zd\n",
	      xmean,xsd,xskew,xkurt,xadev,xmin,xmax,size);
    }
    /* ++++++++++++++++++++++++++++ */
    
    if(o_output>=3 && o_output<=4) printhist_discrete(data,size,o_output);
    else {
      
      size_t i;
      double xmin = DBL_MAX;  /*set maximum double */
      double xmax = -DBL_MAX;  /*set to minimum double */

      
      /* determination of xmax and xmin */
      if(o_window){
	xmin=wmin;
	xmax=wmax;
      }
      else
	for(i=0;i<size;i++){
	  if( xmin > data[i]) xmin = data[i]; 
	  if( xmax < data[i]) xmax = data[i];
	}
      /* ------------------------------ */

      printhist(data,size,steps,o_output,o_bintype,ratio,refpoint,xmin,xmax);
    }

  }

  
  return(0);
}