# include "data.h"

int vgroup (int argc, char **argv) {
  
  int i, j, N;
  Vector *xin, *yin, *xout, *yout;
  opihi_flt xmin, xmax, sum;
  double *values;

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

  if (argc != 5) {
    gprint (GP_ERR, "USAGE: vbin <xin> <yin> <xout> <yout>\n");
    gprint (GP_ERR, " group x,y values in bins defined by <xout>\n");
    return (FALSE);
  }

  if ((xin  = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((yin  = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((xout = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((yout = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);

  // re-binning creates a float vector
  ResetVector (yout, OPIHI_FLT, xout[0].Nelements);

  // storage vector
  ALLOCATE (values, double, xin[0].Nelements);

  for (i = 0; i < xout[0].Nelements - 1; i++) {
    xmin = xout[0].elements.Flt[i];
    xmax = xout[0].elements.Flt[i+1];

    N = 0;
    for (j = 0; j < xin[0].Nelements; j++) {
      if (xin[0].elements.Flt[j] < xmin) continue;
      if (xin[0].elements.Flt[j] > xmax) continue;
      values[N] = yin[0].elements.Flt[j];
      N++;
    }
    
    dsort (values, N);
    if (N > 1) {
      sum = (N % 2) ? values[(int)(0.5*N)] : 0.5*(values[N/2] + values[N/2 + 1]);
    } else {
      sum = values[0];
    }

    // measure the stat for this bin
    // sum = 0.0;
    // for (j = 0; j < N; j++) {
    //   sum += values[j];
    // }
    // sum /= N;
    yout[0].elements.Flt[i] = sum;
  }
  free (values);
  return (TRUE);
}
