# include "relphot.h"

void liststats_setmode (StatType *stats, char *strmode) {

  stats->statmode = -1;
  if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
  if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
  if (!strcmp (strmode, "WT_MEAN"))          { stats->statmode = STATS_WT_MEAN; return; }
  if (!strcmp (strmode, "INNER_MEAN"))       { stats->statmode = STATS_INNER_MEAN; return; }
  if (!strcmp (strmode, "INNER_WTMEAN"))     { stats->statmode = STATS_INNER_WTMEAN; return; }
  if (!strcmp (strmode, "CHI_INNER_MEAN"))   { stats->statmode = STATS_CHI_INNER_MEAN; return; }
  if (!strcmp (strmode, "CHI_INNER_WTMEAN")) { stats->statmode = STATS_CHI_INNER_WTMEAN; return; }

  fprintf (stderr, "ERROR: invalid stats mode: %s\n", strmode);
  exit (1);
}

int liststats_init (StatType *stats) {
  stats->median  = NAN;
  stats->mean    = NAN;
  stats->sigma   = NAN;
  stats->error   = NAN;
  stats->chisq   = NAN;
  stats->min     = NAN;
  stats->max     = NAN;
  stats->Upper80 = NAN;
  stats->Lower20 = NAN;
  stats->total   = NAN;
  stats->Nmeas = 0;
  return TRUE;
}

int liststats (double *value, double *dvalue, double *weight, int N, StatType *stats) {
  
  int i, ks, ke;
  double Mo, dMo, M, dM, Nm, X2, dS, R, W, *chi;

  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");

  liststats_init (stats);

  ke = ks = dMo = 0;

  stats[0].Nmeas = N;
  stats[0].mean  = 0;
  stats[0].sigma = 0;
  stats[0].error = 0;
  stats[0].chisq = 0;
  if (N < 1) return (FALSE);

  int N80 = MIN (N-1, 0.8*N);
  int N20 = MAX (0, 0.2*N);

  if (weight) {
    dsortthree (value, dvalue, weight, N);
  } else {
    if (dvalue) {
      dsortpair (value, dvalue, N);
    } else {
      dsort (value, N);
    }
  }

  // these values do not depend on the errors or weighting scheme
  stats[0].median  = value[(int)(0.5*N)];
  stats[0].min     = value[0];
  stats[0].max     = value[N-1];
  stats[0].Upper80 = value[N80];
  stats[0].Lower20 = value[N20];

  switch (stats->statmode) {
    case STATS_MEDIAN:
      ks = 0;
      ke = N;
      Mo = stats[0].median;
      Nm = N;
      goto chisq;
      break;
    case STATS_MEAN:
    case STATS_WT_MEAN:
      ks = 0;
      ke = N;
      break;
    case STATS_INNER_MEAN:
    case STATS_INNER_WTMEAN:
    case STATS_CHI_INNER_MEAN:
    case STATS_CHI_INNER_WTMEAN:
      ks = 0.25*N + 0.50;
      ke = 0.75*N + 0.25;
      if (N <= 3) {
	ks = 0;
	ke = N;
      }
      break;
    case STATS_NONE:
      myAbort ("undefined stats");
  }    

  // for these two modes, I need a vector of the chi-square contribution
  // I'm actually just using chisq to get the correct sorting order
  if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) {
    ALLOCATE (chi, double, N);
    for (i = 0; i < N; i++) {
      chi[i] = (value[i] - stats[0].median) / dvalue[i];
    }
    if (weight) {
      dsortfour (chi, value, dvalue, weight, N);
    } else {
      dsortthree (chi, value, dvalue, N);
    }
    free (chi);
  }

  int WeightedMean = FALSE;
  WeightedMean |= (stats->statmode == STATS_WT_MEAN);
  WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
  WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);

  /* calculating the per-star offset based on the desired weighting scheme */
  M = dM = Nm = W = R = 0;
  if (weight) {
    // the weight value is multiplied by whichever nominal weighting scheme is provided 
    // thus the user should set weight to 1.0 for nominal weight, and 100 for heavy weight (or so)
    // and 0.01 for under-weight
    if (WeightedMean) {
      for (i = ks; i < ke; i++) {
	M  += value[i] * weight[i] / SQ(dvalue[i]);
	W  +=            weight[i] / SQ(dvalue[i]);
	dM += SQ (weight[i] / dvalue[i]);
	R  += weight[i] / SQ(dvalue[i]);
	Nm += 1.0;  
      }	
      Mo  = M / W;
      dMo = sqrt (dM) / R;
    } else {
      for (i = ks; i < ke; i++) {
	M  += value[i] * weight[i];
	W  +=            weight[i];
	dM += SQ (weight[i] * dvalue[i]);
	R  += weight[i];
	Nm += 1.0;  
      }	
      Mo  = M / W;
      dMo = sqrt (dM) / R;
    }
  } else {
    // NULL weight vector is supplied, revert to standard form (above reverts if weight[i] == 1)
    if (WeightedMean) {
      // weighted by inverse-variance
      for (i = ks; i < ke; i++) {
	M   += value[i] / SQ (dvalue[i]);
	dM  += 1.0 / SQ (dvalue[i]);
	Nm  += 1.0;  
      }	
      Mo = M / dM;
      dMo = sqrt (1.0 / dM);
    } else {
      // pure un-weighted
      for (i = ks; i < ke; i++) {
	M   += value[i];
	dM  += SQ (dvalue[i]);
	Nm  += 1.0;  
      }	
      Mo = M / Nm;
      dMo = sqrt (dM) / Nm;
    }
  }

 chisq:
  /* find sigma and chisq */
  X2 = dS = 0;
  for (i = ks; i < ke; i++) {
    M  = SQ (value[i] - Mo);
    dM = SQ (dvalue[i]);
    X2 += M / dM;
    dS += M;
  }
  X2 = X2 / (Nm - 1);
  dS = sqrt (dS / (Nm - 1));

  stats[0].mean  = Mo;
  stats[0].Nmeas = Nm;
  stats[0].chisq = X2;
  stats[0].sigma = dS;
  stats[0].error = dMo;

  return (TRUE);
}

// we could define the weight to be the only scale factor:
// \mu      = \sum (value_i * weight_i) / \sum (weight_i)
// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
// R = \sum (weight_i^2)

// or, we could define the weight to be a scale factor times the inverse error:
// \mu      = \sum (value_i * weight_i / sigma_i) / \sum (weight_i)
// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
// R = \sum (weight_i^2)
