# include "data.h"

int vstat (int argc, char **argv) {
  
  int i, N;
  double max, min, sum, var, dvar, mean, stdev;
  float *X, IgnoreValue;
  int Ignore, Quiet;

  int *Nval, bin, Nmode, Nmed;
  double dx, mode, median;
  Vector *vec;

  Ignore = FALSE;
  if (N = get_argument (argc, argv, "-ignore")) {
    Ignore = TRUE;
    remove_argument (N, &argc, argv);
    IgnoreValue = atof(argv[N]);
    remove_argument (N, &argc, argv);
  }

  Quiet = FALSE;
  if (N = get_argument (argc, argv, "-q")) {
    Quiet = TRUE;
    remove_argument (N, &argc, argv);
  }
  if (N = get_argument (argc, argv, "-quiet")) {
    Quiet = TRUE;
    remove_argument (N, &argc, argv);
  }

  if (argc != 2) {
    fprintf (stderr, "USAGE: vstat (vector)\n");
    return (FALSE);
  }

  if ((vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);

  /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */

  /* calculate max, min, mean, sum, npix */
  X = vec[0].elements;
  max = min = *X;
  sum = N = 0;
  for (i = 0; i < vec[0].Nelements; i++, X++) {
    if (!finite (*X)) continue;
    if (Ignore && (*X == IgnoreValue)) continue;
    max = MAX (*X, max);
    min = MIN (*X, min);
    sum += *X;
    N++;
  }      
  mean = sum / N;

  /* calculate median and mode with resolution of (max - min) / 1000 */ 
  dx = (max - min) / 1000;
  if (dx == 0) {
    median = mode = min;
    stdev = 0.0;
    goto skip;
  }

  ALLOCATE (Nval, int, 1002);
  bzero (Nval, 1000*sizeof(int));
  X = vec[0].elements;
  var = 0;
  for (i = 0; i < vec[0].Nelements; i++, X++) {
    if (!finite (*X)) continue;
    if (Ignore && (*X == IgnoreValue)) continue;
    bin = MAX (0, MIN (1000, (*X - min) / dx));
    Nval[bin] ++;
    dvar = (*X - mean);
    var += dvar*dvar;
  }      
  stdev = sqrt (var / N);

  Nmode = 0;
  mode = Nval[Nmode];
  median = 0;
  Nmed = -1;
  for (i = 0; i < 1001; i++) {
    if (Nmed == -1) {
      median += Nval[i];
      if (median >= N / 2.0) {
	Nmed = i;
	median = i * dx + min;
      }
    }
    if (mode < Nval[i]) {
      Nmode = i;
      mode = Nval[Nmode];
    }
  }
  mode = Nmode * dx + min;
  free (Nval);

skip:
  if (!Quiet) {
    fprintf (stderr, "mean: %f, stdev: %f, min: %f, max: %f, median: %f, mode: %f, Npts: %d\n", 
	     mean, stdev, min, max, median, mode, N);
  }

  set_variable ("MIN",      min);
  set_variable ("MAX",      max);
  set_variable ("MEDIAN",   median);
  set_variable ("MEAN",     mean);
  set_variable ("MODE",     mode);
  set_variable ("TOTAL",    sum);
  set_int_variable ("NPIX", N);
  set_variable ("SIGMA",    stdev);
  return (TRUE);
}

