/*                                                +++++++++        */
# define FILENAME                                 "stats.c"
/*                                                +++++++++        */

# include "util.h"

/*          =========                                                */
Statistic * InitStats (Statistic * S)
/*          =========                                                */
{
  static i = 0;
  bzero(S,sizeof(Statistic));
  S->name = "NoName";
  S->is_initialized = true ;
  S->min =  INFINITY ;
  S->max = -INFINITY ;
  i++;
  return S ;
}

/*   =========                                                */
void DumpStats (Statistic * S, FILE * stream)
/*   =========                                                */
{
  char m[BUFSIZ], M[BUFSIZ];
  char * minP, *maxP;

  sprintf(m,"%f",S->min);
  sprintf(M,"%f",S->max);
  if (S->min == INFINITY) 
    minP = "INF" ;
  else 
    minP = m;
  if (S->max <= -INFINITY) 
    maxP = "-INF" ;
  else 
    maxP = M;
  
  fprintf(stream,
	  "name= %s num= %d mean= %f min= %s max= %s sum= %f dev= %f sumsq= %f AvgSumSq= %f FracDev= %f\n",
	  S->name,(int) S->num, S->mean, minP,maxP,S->sum,S->dev,S->sumsq,S->SqrtAvgSumSq,S->FracDev) ;
}

/*   ==========                                                */
void NDumpStats (Statistic * S, FILE * stream)
/*   ==========                                                */
{
# define HEAD(S) fprintf(stream,"\nNDumpStats on %s:  ",S->name)
# define FIELD(f) fprintf(stream,"%s= %f   ",MAKESTRING(f),S->f)
  
  HEAD(S);
  fprintf(stream,"error= %1d   ",S->error);  
  FIELD(num); FIELD(SumOfWt); FIELD(AvgWt);
  HEAD(S);
  FIELD(min); FIELD(max);
  HEAD(S);
  FIELD(sum); FIELD(mean); 
  HEAD(S);
  FIELD(Var); FIELD(dev); FIELD(FracDev);
  HEAD(S); 
  FIELD(WtSum); FIELD(WtAvg); 
  HEAD(S);
  FIELD(WtVar); FIELD(WtDev); FIELD(WtFracDev);
  fprintf(stream,"\n");
}

/*   ==============                                                */
void IncrementStats(Statistic * S,double v)
/*   ==============                                                */
{
  double var ;
  if (!S->is_initialized) InitStats(S);
  S->num++;
  S->sum += v;
  S->sumsq += v * v ;
  if (v > S->max) S->max = v ;
  if (v < S->min) S->min = v ;
  assert(S->num != 0);
  S->mean  = S->sum   / S->num ;
  var = S->sumsq/S->num - S->mean*S->mean ;
  if (var < 0) var = 0;
  assert(var >= 0);
  S->dev   = sqrt(var  ) ;
  S->SqrtAvgSumSq = sqrt(S->sumsq)/S->num ;
  S->rms = sqrt(S->sumsq/S->num);
  if (!FLOATEQ(S->mean,0))
    S->FracDev = S->dev / S->mean ;
}


/*   ===============                                                */
void IncStatsWithWts(Statistic * S,double v, double wt)
/*   ===============                                                */
{
  if (!S->is_initialized) InitStats(S);
  S->num ++;
  S->sum += v;
  S->WtSum += wt*v;
  S->sumsq += v * v ;
  S->WtSumSq += wt * v * v ;
  S->SumOfWt += wt;
  S->AvgWt = S->SumOfWt / S->num ;
  if (v > S->max) S->max = v ;
  if (v < S->min) S->min = v ;
  S->mean  = S->sum   / S->num ;
  S->WtAvg = S->WtSum / S->SumOfWt ;
  S->WtAvgSumSq = S->WtSumSq / S->SumOfWt ;
  {
    double d;

#   define MAX_ERROR 2

#   define SIGNAL_ERROR(S,d)  \
    { \
      S->error++ ; \
      if (S->error < MAX_ERROR) { \
	fprintf(stderr,"%s [%d] (IncStatsWithWts): Error %3d, Value= %f\n", \
		__FILE__,__LINE__,S->error,d); \
      } \
    } 

#   define CHECK_SET(S,cond,field,val) \
    { \
      if (cond) { \
	S->field = val ; \
      } else { \
        SIGNAL_ERROR(S,d); \
	S->field = MBG_NAN ; \
      } \
    } 

#   define FIX_LESS_THAN_0(d) \
    { \
      if (d<0) { \
        SIGNAL_ERROR(S,d); \
        d = fabs(d); \
      } \
    }

    S->WtVar = S->WtAvgSumSq - S->WtAvg * S->WtAvg ;
    FIX_LESS_THAN_0(S->WtVar);
    CHECK_SET(S,S->WtVar>=0,WtDev,sqrt(S->WtVar)) ;

    S->Var = S->sumsq/S->num - S->mean*S->mean ;
    FIX_LESS_THAN_0(S->Var);
    CHECK_SET(S,S->Var>=0,dev,sqrt(S->Var)) ;

    d = S->sumsq/S->num ;
    FIX_LESS_THAN_0(d);
    CHECK_SET(S,d>=0,SqrtAvgSumSq,sqrt(d)) ;
    
    CHECK_SET(S,!FLOATEQ(S->mean,0),FracDev,S->dev / S->mean);
    CHECK_SET(S,!FLOATEQ(S->WtAvg,0),WtFracDev,S->WtDev / S->WtAvg);
  }
}

/*   =====================                                                */
void ReComputeDerivedStats(Statistic * S)
/*   =====================                                                */
{
  S->mean  = S->sum   / S->num ;
  S->WtAvg = S->WtSum / S->SumOfWt ;
  S->WtAvgSumSq = S->WtSumSq / S->SumOfWt ;
  S->WtVar = S->WtAvgSumSq - S->WtAvg * S->WtAvg ;
  S->Var = S->sumsq/S->num - S->mean*S->mean ;
  ReComputeDerivedStats2(S) ;
}

/*   ======================                                                */
void ReComputeDerivedStats2 (Statistic * S)
/*   ======================                                                */
{
  double d;
  
#   define MAX_ERROR 2
  
#   define SIGNAL_ERROR(S,d)  \
    { \
      S->error++ ; \
      if (S->error < MAX_ERROR) { \
	fprintf(stderr,"%s [%d] (IncStatsWithWts): Error %3d, Value= %f\n", \
		__FILE__,__LINE__,S->error,d); \
      } \
    } 

#   define CHECK_SET(S,cond,field,val) \
    { \
      if (cond) { \
      S->field = val ; \
      } else { \
        SIGNAL_ERROR(S,d); \
	S->field = MBG_NAN ; \
      } \
    } 

#   define FIX_LESS_THAN_0(d) \
    { \
      if (d<0) { \
        SIGNAL_ERROR(S,d); \
        d = fabs(d); \
      } \
    }

    FIX_LESS_THAN_0(S->WtVar);
    CHECK_SET(S,S->WtVar>=0,WtDev,sqrt(S->WtVar)) ;

    FIX_LESS_THAN_0(S->Var);
    CHECK_SET(S,S->Var>=0,dev,sqrt(S->Var)) ;

    d = S->sumsq/S->num ;
    FIX_LESS_THAN_0(d);
    CHECK_SET(S,d>=0,SqrtAvgSumSq,sqrt(d)) ;
    
    CHECK_SET(S,!FLOATEQ(S->mean,0),FracDev,S->dev / S->mean);
    CHECK_SET(S,!FLOATEQ(S->WtAvg,0),WtFracDev,S->WtDev / S->WtAvg);
}