/* See {pz_histogram.h} */
/* Last edited on 2015-01-20 16:46:38 by stolfilocal */

#include <math.h>
#include <pz_types.h>

#include <pz_histogram.h>

pz_histogram_t *pz_histogram_new(double lo, double hi, unsigned m)
  {
    affirm(m >= 2);
    affirm(hi > lo);

    { /* with */ 
      ??? r = NEW(REF pz_histogram_t);
      ??? h = r^;
      /* do */
      h = pz_histogram_t{
        lo = lo, hi = hi,
        m = m,
        c = double_vec_new(m),
        wSum = 0.0e0,
        wxSum = 0.0e0,
        wxxSum = 0.0e0
      };
      for (i = 0; i < (h.c^.ne ) ; i++){ h.c^[i] = 0.0e0 ;};
      return r;
    }
  } /* pz_histogram_new */


void pz_histogram_add(pz_histogram_t *h, double x, double  w)
  {

    { /* with */
      ??? xx = MIN(h.hi, MAX(h.lo, x));
      ??? k = ROUND(((double)h.m-1)*(xx - h.lo)/(h.hi - h.lo));
      ??? bin = h.c^[k];
      /* do */
      bin = bin + w;
      h.wSum = h.wSum + w;
      h.wxSum = h.wxSum + w*x;
      h.wxxSum = h.wxxSum + w*x*x;
    }
  } /* pz_histogram_add */

void pz_histogram_output(FILE *wr, pz_histogram_t *h)
  {

    { /* with */
      ??? m = h.m;
      ??? c = h.c^;
      ??? step = (h.hi - h.lo)/((double)m-1);
      /* do */
      for (i = 0; i <= m-1 ; i++){

        { /* with */
          ??? fi = ((double)i);
          ??? xLo = h.lo + step * (fi - 0.5e0);
          ??? xMd = h.lo + step * fi;
          ??? xHi = h.lo + step * (fi + 0.5e0);
          /* do */
          fprintf(wr, " ");
          fprintf(wr, FLR(xLo, 13,6));
          fprintf(wr, " ");
          fprintf(wr, FLR(xMd, 13,6));
          fprintf(wr, " ");
          fprintf(wr, FLR(xHi, 13,6));
          fprintf(wr, " ");
          fprintf(wr, FLR(c[i], 13,6));
          fprintf(wr, "\n");;
        };
      };;
    }
  } /* pz_histogram_output */


double pz_histogram_avg(pz_histogram_t *h)
  {
    if (( h.wSum == 0.0e0 )){
      return 0.0e0
    }else{ 
      return h.wxSum/h.wSum;
    }
  } /* pz_histogram_avg */
  

double pz_histogram_var(pz_histogram_t *h)
  {
    if (( h.wSum == 0.0e0 )){
      return 0.0e0
    }else{
      return (h.wxxSum - h.wxSum*h.wxSum/h.wSum)/h.wSum;
    }
  } /* pz_histogram_var */
  

double pz_histogram_dev( pz_histogram_t *h)
  {
    return Math.sqrt(pz_histogram_var(h))
  } /* pz_histogram_dev */


char *FLR( double x, unsigned w, unsigned p)
  {
    return Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, prec = p), w)
  } /* FLR */