#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "hbook.h"

/*====================================================*/

#define DEBUG 0
#define MAXX 64
#define MAXY  25

char* hprint(SHM_HIST  *hp, int id) {  
  char *cplot;
  cplot=mprint(hp[id].hist, hp[id].nbin, hp[id].title);
  sprintf(&cplot[strlen(cplot)],"\t\tID=%d, Stat: Entries=%d, Underflow=%d, Overflow=%d \n"
	  ,hp[id].id,hp[id].entries,hp[id].under,hp[id].over);
  return cplot;
}
char* mprint(unsigned int  *hist, int nbin, const char *title) {
  int i,j,k,v,s,s0,s1,s2,s3;    
  int ybin,y5,miny=15,maxy=MAXY,mod,entry;
  float sum,sum2,mean,mean2,ey2,rms,x;
  char st[MAXX+1];
  static char chist[MAXX*MAXY*2];
  /*----------------------------------------------------*/
  int sum_bin,max,max1;
  unsigned int pd[MAXX];
  printf("\n");

  //printf("hprint:: hist=%p len=%d \n",hist,nbin);
  memset(chist, 0, sizeof(chist));
  sprintf(chist,"----------------------------------------------------------------------------\n");
   /*---  X scale (rebin) ---*/

  for (i=1;i<32768;i++)  if ((nbin/i)<=MAXX) break;
  sum_bin=i;
  if (DEBUG) sprintf(&chist[strlen(chist)],"  Xscale: nbin=%d sum_bin=%d  \n",nbin,sum_bin);

  /*---  Y scale ---*/

  i=0; j=0; max=0; max1=0; entry=0;
  while(1) {
    pd[i]=0;
    for (k=0;k<sum_bin;k++) { 
	if (j<nbin) { pd[i]+=hist[j]; if ((int)hist[j]>max1) max1=(int)hist[j]; } else pd[i]+=0; 
      j++; 
    }
    if ((int)pd[i]>max) max=pd[i];  entry+=pd[i];
    if (DEBUG>1) 
      printf(" pd[%d]=%d max=%d \n",i,pd[i],max);
    i++; if (i==MAXX) break;
  }

  for (j=10;j>0;j--) {
    mod=j%10;
    for (i=maxy;i>miny;i--) {
      ybin=max/i+1; 
      if (DEBUG>1) 
        printf(" Yscale: maxY=%d, max=%d, ybin=%d mod=%d\n",i,max,ybin,mod);
      if ((ybin%10)==mod) goto a;
    }
  }
a:
  if      (max<=25) { maxy=max; ybin=1; }
  else if (i!=miny) { maxy=i; 
                      if (mod) {   if (ybin>=7)  ybin+=10-mod; 
			maxy=max/ybin+1;
			if (maxy<16) maxy=16;
		      } 
                    } 
  else              { ybin=max/maxy+1; }
  if (ybin==0)  ybin=1;

  /*-------------  Statistics   ------------------------*/
  mean=0; rms=0;
  if (entry>0) {
    sum=0; sum2=0;
    for (i=0;i<nbin;i++) {     x=i+0.5;
      sum +=(float)hist[i]*x;
      sum2+=(float)hist[i]*x*x;
    } 
    mean=sum/entry; mean2=mean*mean; ey2=sum2/entry;
    if (DEBUG>1) 
      printf(" sum=%f, sum2=%f, \n mean=%f, mean2=%f, E(y**2)=%f\n"
	     ,sum,sum2,mean,mean2,ey2);
    rms=(float)sqrt((double)((float)nbin/(float)(nbin-1)*(ey2-mean2)));
  }
  sprintf(&chist[strlen(chist)],"  Title: %s \n",title);
  sprintf(&chist[strlen(chist)],"  Entry=%d, Max=%d(%d), Mean=%.2f, RMS=%.2f  \n",entry,max,max1,mean,rms);
  /*----------------------------------------------------*/
  if (DEBUG>1) printf(" Yscale: maxY=%d  max=%d  ybin=%d\n",maxy,max,ybin);
  max=maxy*ybin; y5=ybin/5;
  if (DEBUG) 
      sprintf(&chist[strlen(chist)],"  Yscale: maxY=%d  max=%d  ybin=%d  y5=%d \n",maxy,max,ybin,y5);

  /*----- boba -----------------------------------------*/

  for (j=1;j<=maxy;j++) {
    s=max-ybin*j;   st[MAXX]=0;
    s0=s+y5; s1=s0+y5;  s2=s1+y5;  s3=s2+y5;
    for (i=0;i<MAXX;i++) {   /*  fill line  */
      v=*(pd+i);      st[i]=' ';
      if (0<v && v<=ybin && j==maxy && ybin!=1)   {st[i]='.';  }
      if (v>s)       {st[i]='.';     }
      if (v>s0)      {st[i]=',';     }
      if (v>s1)      {st[i]='+';     }
      if (v>s2)      {st[i]=';';     }
      if (v>s3)      {st[i]='X';     }
      if (DEBUG>2) 
	printf(" v=%d, s0123=%d %d %d %d %d ybin=%d, y5=%d, i=%d st[i]=%c \n"
	      ,v,s,s0,s1,s2,s3,ybin,y5,i,st[i]);
    } 
    st[MAXX]=0;
    sprintf(&chist[strlen(chist)],"%8d   I%s\n",s+ybin,st);
    //printf("%8d   I%s\n",s+ybin,st);
  }
  for (i=0;i<MAXX;i++) { st[i]='0'+i%10; } st[MAXX]=0;
  sprintf(&chist[strlen(chist)],"            %s\n",st);
  sprintf(&chist[strlen(chist)],"   ");
  //printf("            %s\n",st);
  //printf("   ");
  for (i=0;i<MAXX/9;i++) { sprintf(&chist[strlen(chist)],"%10d",10*i*sum_bin); /*printf("%10d",10*i*sum_bin);*/ }
  sprintf(&chist[strlen(chist)],"\n");
  //printf("\n");
  //puts(chist);
  //for (int k=0;k<strlen(chist);k++) printf(" %d: %d %x %c \n",k,chist[k],chist[k],chist[k]);
  printf("\n strlen(chist)=%d  size(chist)=%d\n",(int)strlen(chist),(int)sizeof(chist));
  return chist;
}

/*=============================================================*/
/*           fill                                              */
/*=============================================================*/
void mfill(unsigned int *hist, int n, float x, float w) {
  int ich;
  ich=(int)x; 
  if (ich>=n) { ich=n-1;/* printf(" mfill: Over:  n=%d x=%f \n",n,x);*/}
  if (ich<0)  { ich=0; /*  printf(" mfill: Under: n=%d x=%f \n",n,x);*/}
  hist[ich]+=(int)w;
  if (DEBUG>2) printf(" MFILL: x=%f, ich=%d w=%f \n",x,ich,w);
  /* printf(" MFILL: hist=%p Len=%d \n",hist,n); */
}
void mf1(unsigned int *hist, int n, int ich) {
    if (ich>=n) { ich=n-1; /*  printf(" mfill: Over:  n=%d x=%f \n",n,x);*/}
    if (ich<0)  { ich=0;   /*  printf(" mfill: Under: n=%d x=%f \n",n,x);*/}
    hist[ich]++;
    //if (DEBUG>2) printf(" MFILL_1: x=%f, ich=%d w=%d \n",x,ich,w);
    /* printf(" MFILL_1: hist=%p Len=%d \n",hist,n); */
}
inline void hf1(SHM_HIST *hp, int id, unsigned int ich) {
  hp[id].entries++;
  if (ich>=hp[id].nbin) { hp[id].over++; /*  printf(" mfill: Over:  n=%d x=%f \n",n,x);*/  return;}
  //if (ich<0)  { hp[id].under++;   /*  printf(" mfill: Under: n=%d x=%f \n",n,x);*/     return;}
  hp[id].hist[ich]++;
  //if (DEBUG>2) printf(" MFILL_1: x=%f, ich=%d w=%d \n",x,ich,w);
  /* printf(" MFILL_1: hist=%p Len=%d \n",hist,n); */
}
/*=============================================================*/
/*           book                                              */
/*=============================================================*/
void hbook1(SHM_HIST *hp, int id, char *title, int nbin, float  xmin,  float xmax) {
  
  strncpy(hp[id].title,title,shm_hist_LT);  
  if (strlen(title)>=shm_hist_LT) hp[id].title[shm_hist_LT-1]=0;
  if (nbin>shm_hist_NB) nbin=shm_hist_NB;
  hp[id].id=id;
  hp[id].nbin=nbin;  hp[id].xmin=xmin;  hp[id].xmax=xmax;  hp[id].xbin=(xmax-xmin)/nbin;
  hp[id].under=0; hp[id].over=0;  hp[id].integral=0;  hp[id].entries=0;
  mreset(hp[id].hist, hp[id].nbin);
}
/*=============================================================*/
/*           reset                                             */
/*=============================================================*/
void hreset0(SHM_HIST *hp, int nh) {
  int i;
  for (i=0; i<nh; i++) {
    mreset(hp[i].hist, hp[i].nbin);
    hp[i].under=0; hp[i].over=0; hp[i].integral=0;  hp[i].entries=0;
  }
}
void hreset(SHM_HIST *hp, int id) {
  mreset(hp[id].hist, hp[id].nbin);
  hp[id].under=0; hp[id].over=0; hp[id].integral=0;  hp[id].entries=0;
}
void mreset(unsigned int *hist, int n) {
  int i;
  //printf("mreset:: hist=%p len=%d \n",hist,n);
  for (i=0;i<n;i++)  hist[i]=0;
}
/*=============================================================*/
/*           copy                                              */
/*=============================================================*/
void mcopy(unsigned int *hist1, unsigned int *hist2, int n) {
  int i;
  for (i=0;i<n;i++)  hist2[i]=hist1[i];
}
/*=============================================================*/
/*           statistics                                        */
/*=============================================================*/
float hstati(int sw, SHM_HIST  *hp, int id) {
  return mstati(sw, hp[id].hist, hp[id].nbin);
}
int hnoent(SHM_HIST  *hp,  int id) {
  return hp[id].entries;
}
float mstati(int sw, unsigned int  *hist, int nbin) {
  int i,entry;
  float sum,sum2,mean,mean2,ey2,rms,x;
  /*-------------  Statistics   ------------------------*/
  entry=0; sum=0; sum2=0;
  for (i=0;i<nbin;i++) {     x=i+0.5;
    entry+=hist[i];
    sum +=(float)hist[i]*x;
    sum2+=(float)hist[i]*x*x;
  } 
  if (entry==0) return 0;
  mean=sum/entry; mean2=mean*mean; ey2=sum2/entry;
  if (DEBUG>1) 
    printf(" sum=%f, sum2=%f, \n mean=%f, mean2=%f, E(y**2)=%f\n"
	   ,sum,sum2,mean,mean2,ey2);
  rms=(float)sqrt((double)((float)nbin/(float)(nbin-1)*(ey2-mean2)));
  /*printf(" Entry=%d, Mean=%.2f, RMS=%.2f \n",entry,mean,rms);*/  
  /*----------------------------------------------------*/
  if     (sw==1)  return (float)entry; 
  else if (sw==3) return rms;   
  else            return mean;  

}
/*=============================================================*/
/*                                                             */
/*=============================================================*/
float gauss()
{
  int i;
  float s = -6.0;
  for(i=0;i<12;i++) s += ((float)rand())/((float)RAND_MAX);
  return s;
}
/*=============================================================*/
/*                                                             */
/*=============================================================*/