# include "relastro.h"

static int Nskip1 = 0;
static int Nskip2 = 0;
static unsigned int Tref = 1323353985;
static short Cref = 10355;

int bcatalog (Catalog *subcatalog, Catalog *catalog) {

  off_t i, j, offset;
  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
  int Nsecfilt;

  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
  Nsecfilt = GetPhotcodeNsecfilt ();
  assert (catalog[0].Nsecfilt == Nsecfilt);

  /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
  memset(subcatalog, 0, sizeof(Catalog));
  NAVERAGE = 50;
  NMEASURE = 1000;
  ALLOCATE (subcatalog[0].average,  Average,     NAVERAGE);
  ALLOCATE (subcatalog[0].secfilt,  SecFilt,     NAVERAGE*Nsecfilt);
  ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
  Nmeasure = Naverage = 0;

  DVOAverageFlags averageBits = 
    ID_STAR_FEW     |
    ID_STAR_POOR    |
    ID_STAR_FIT_AVE |
    ID_STAR_FIT_PM  |
    ID_STAR_FIT_PAR |
    ID_STAR_USE_AVE |
    ID_STAR_USE_PM  |
    ID_STAR_USE_PAR; 

  /* exclude stars not in range or with too few measurements */
  for (i = 0; i < catalog[0].Naverage; i++) {
    if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
      continue;
      // XXX can't (easily) count this for a given chip
    }
    
    /* start with all stars good */
    subcatalog[0].average[Naverage] = catalog[0].average[i];
    subcatalog[0].average[Naverage].measureOffset = Nmeasure;
    for (j = 0; j < Nsecfilt; j++) {
      subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j];
    }

    if (RESET) {
      subcatalog[0].average[Naverage].flags &= ~averageBits;
    }

    Nm = 0;
    for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {

      offset = catalog[0].average[i].measureOffset + j;
      
      // filter objects based on user supplied criteria, including SIGMA_LIM
      if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) {
	catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
	if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
	  Nskip1 ++;
	}
	continue;
      }

      // filter out outliers
      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
	catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
	if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
	  Nskip2 ++;
	}
	continue;
      }
      catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP;

      // re-assess on each run of relastro if a measurement should be used

      // NOCAL : this flag is used to mark measurements which must be ignored for this
      // analysis.  an example would be time ranges or photcodes which are ignored.  NOCAL
      // is internal only; it is not written out.

      // SKIP : this flag is used to mark measurements currently ignored, but may be
      // allowed before the analysis is complete.  for example, a star with measurements
      // only near the edge of images will have this flag set for most of the analysis,
      // but before the final average properties are calculated, these measurements may be
      // allowed.

      CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
      // subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset];
      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
      subcatalog[0].measureT[Nmeasure].averef   = Naverage;
      if (RESET) {
        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
      }

      Nmeasure ++;
      Nm ++;
      if (Nmeasure == NMEASURE) {
        NMEASURE += 1000;
        REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
      }
    }
    subcatalog[0].average[Naverage].Nmeasure = Nm;
    Naverage ++;
    if (Naverage == NAVERAGE) {
      NAVERAGE += 50;
      REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
      REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
    }
  }
  REALLOCATE (subcatalog[0].average,  Average,     MAX (Naverage, 1));
  REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1));
  REALLOCATE (subcatalog[0].secfilt,  SecFilt,     MAX (Naverage, 1)*Nsecfilt);
  subcatalog[0].Naverage = Naverage;
  subcatalog[0].Nmeasure = Nmeasure;
  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
  subcatalog[0].catID    = catalog[0].catID;
  assert (Nsecfilt == catalog[0].Nsecfilt);

  // limit the total number of stars in the catalog
  if (MaxDensityUse) {
    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
  } else {
    if (VERBOSE2) {
      char *basename = filebasename (catalog[0].filename);
      fprintf (stderr, "subset of "OFF_T_FMT" ("OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" total) measures for catalog %s\n", 
	       subcatalog[0].Naverage, catalog[0].Naverage, subcatalog[0].Nmeasure, catalog[0].Nmeasure, basename);
      free (basename);
    }
  }
  return (TRUE);
}

void bcatalog_show_skips () {
  fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
}

// sort by decreasing Nmeasure (X)
void sort_by_Nmeasure (int *X, off_t *Y, off_t N) {

# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;	\
  tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
  tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
}
# define COMPARE(A,B)(X[A] > X[B])

  OHANA_SORT (N, COMPARE, SWAPFUNC);

# undef SWAPFUNC
# undef COMPARE

}

/* this version does NOT use AverageTiny, MeasureTiny */ 
int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) {

  Catalog tmpcatalog;

  double Rmin, Rmax, Dmin, Dmax;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
  gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax);

  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);

  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
  assert (AREA > 0);

  off_t Nmax = MaxDensityValue * AREA;
  if (subcatalog[0].Naverage <= Nmax) {
    if (VERBOSE) {
      fprintf (stderr, "subcatalog has less than the max density\n");
    }
    return (TRUE);
  }

  off_t Naverage = subcatalog[0].Naverage;

  // select a random subset of Nmax stars from subcatalog using Fisher-Yates

  // we are going to select Nmax entries by choosing the brightest objects
  int *value;
  off_t *index, i, j, ave;
  ALLOCATE (index, off_t, Naverage);
  ALLOCATE (value, int, Naverage);
  for (i = 0; i < Naverage; i++) {
    index[i] = i;
    value[i] = subcatalog[0].average[i].Nmeasure;
  }
  sort_by_Nmeasure (value, index, Naverage);

  // count the number of measurements this selection will yield
  off_t NMEASURE = 0;
  for (i = 0; i < Nmax; i++) {
    ave = index[i];
    NMEASURE += subcatalog[0].average[ave].Nmeasure;
  }

  // allocate the output data 
  ALLOCATE (tmpcatalog.average,  Average,     Nmax);
  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
  ALLOCATE (tmpcatalog.secfilt,  SecFilt,     Nmax * Nsecfilt);

  off_t Nmeasure = 0;

  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
  for (i = 0; i < Nmax; i++) {
    ave = index[i];
    tmpcatalog.average[i] = subcatalog[0].average[ave];
    tmpcatalog.average[i].measureOffset = Nmeasure;
    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
      off_t offset = subcatalog[0].average[ave].measureOffset + j;
      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
      tmpcatalog.measureT[Nmeasure].averef = i;
      Nmeasure ++;
    }
    for (j = 0; j < Nsecfilt; j++) {
      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
    }
  }

  if (VERBOSE2) {
    char *basename = filebasename (oldcatalog[0].filename);
    fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", 
	     Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
    free (basename);
  }

  free (index);
  free (value);
  free (subcatalog[0].average);
  free (subcatalog[0].measureT);
  free (subcatalog[0].secfilt);

  subcatalog[0].average = tmpcatalog.average;
  subcatalog[0].measureT = tmpcatalog.measureT;
  subcatalog[0].secfilt = tmpcatalog.secfilt;
  subcatalog[0].Naverage = Nmax;
  subcatalog[0].Nmeasure = Nmeasure;
  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;

  return (TRUE);
}

/* this version does NOT use AverageTiny, MeasureTiny */ 
int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog) {

  Catalog tmpcatalog;

  double Rmin, Rmax, Dmin, Dmax;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);

  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);

  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
  assert (AREA > 0);

  off_t Nmax = MaxDensityValue * AREA;
  if (subcatalog[0].Naverage <= Nmax) {
    if (VERBOSE2) {
      fprintf (stderr, "subcatalog has less than the max density\n");
    }
    return (TRUE);
  }

  off_t Naverage = subcatalog[0].Naverage;

  // select a random subset of Nmax stars from subcatalog using Fisher-Yates

  // we are going to select Nmax entries by generating a random-sorted index list
  off_t *index, tmp, i, j, ave;
  ALLOCATE (index, off_t, Naverage);
  for (i = 0; i < Naverage; i++) {
    index[i] = i;
  }
  for (i = 0; i < Naverage; i++) {
    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
    tmp = index[j];
    index[j] = index[i];
    index[i] = tmp;
  }

  // count the number of measurements this selection will yield
  off_t NMEASURE = 0;
  for (i = 0; i < Nmax; i++) {
    ave = index[i];
    NMEASURE += subcatalog[0].average[ave].Nmeasure;
  }

  // allocate the output data 
  ALLOCATE (tmpcatalog.average,  Average,     Nmax);
  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
  ALLOCATE (tmpcatalog.secfilt,  SecFilt,     Nmax * Nsecfilt);

  off_t Nmeasure = 0;

  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
  for (i = 0; i < Nmax; i++) {
    ave = index[i];
    tmpcatalog.average[i] = subcatalog[0].average[ave];
    tmpcatalog.average[i].measureOffset = Nmeasure;
    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
      off_t offset = subcatalog[0].average[ave].measureOffset + j;
      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
      tmpcatalog.measureT[Nmeasure].averef = i;
      Nmeasure ++;
    }
    for (j = 0; j < Nsecfilt; j++) {
      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
    }
  }

  if (VERBOSE2) {
    fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", 
	     Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  catalog[0].filename);
  }

  free (subcatalog[0].average);
  free (subcatalog[0].measureT);
  free (subcatalog[0].secfilt);

  subcatalog[0].average = tmpcatalog.average;
  subcatalog[0].measureT = tmpcatalog.measureT;
  subcatalog[0].secfilt = tmpcatalog.secfilt;
  subcatalog[0].Naverage = Nmax;
  subcatalog[0].Nmeasure = Nmeasure;
  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
  subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;

  return (TRUE);
}

