# include "relphot.h"

// we've just reloaded the data from disk; we now need to apply the Image/Mosaic/Grid
// calibrations determined by the rest of the program.  We also need to set the final
// output dbFlags values

void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) {

  off_t i;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  /*** RESET photometry flags and Mcal values ***/
  if (RESET) {
    // flags used by the photometry analysis (excluding UBERCAL)
    unsigned int PHOTOM_FLAGS = 
      ID_MEAS_NOCAL |          // detection ignored for this analysis (photcode, time range)
      ID_MEAS_POOR_PHOTOM |    // detection is photometry outlier 
      ID_MEAS_SKIP_PHOTOM |    // detection was ignored for photometry measurement 
      ID_MEAS_AREA |	       // detetion was outside acceptable area of device
      ID_MEAS_SYNTH_MAG |      // magnitude is synthetic
      ID_MEAS_STACK_PRIMARY |  // this stack measurement is in the primary skycell
      ID_MEAS_STACK_PHOT_SRC;  // this measurement supplied the stack photometry

    // ID_MEAS_PHOTOM_UBERCAL -- externally-supplied zero point from ubercal analysis
    // this is set by 'setphot', do not reset here

    // flags used by the photometry analysis (excluding UBERCAL)
    // unsigned int secfiltFlags = 
    //   ID_PHOTOM_PASS_0 | // average measured at pass 0
    //   ID_PHOTOM_PASS_1 | // average measured at pass 1
    //   ID_PHOTOM_PASS_2 | // average measured at pass 2
    //   ID_PHOTOM_PASS_3 | // average measured at pass 3
    //   ID_PHOTOM_PASS_4 | // average measured at pass 3
    //   ID_SECF_USE_SYNTH | // average measured at pass 3
    //   ID_SECF_USE_UBERCAL | // average measured at pass 3
    //   ID_SECF_OBJ_EXT; // average measured at pass 3
    
    // XXX I should really deprecate the concept of applying the average 
    // calculation to a limited set of photcodes.  
    // for now, just do all photcodes here (
    // as it stands, only stacks are limited by photcode; mean exp and forced warp 
    // are applied to all Nsecfilt

    for (i = 0; i < catalog[0].Naverage; i++) {
      int Ns;
      for (Ns = 0; Ns < Nsecfilt; Ns++) {
	
	off_t N = Nsecfilt*i+Ns;
	dvo_secfilt_init (&catalog[0].secfilt[N], SECFILT_RESET_ALL);

	off_t m = catalog[0].average[i].measureOffset;
	off_t j;
	for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
	
	  /* select measurements by time */
	  if (TimeSelect) {
	    if (catalog[0].measure[m].t < TSTART) continue;
	    if (catalog[0].measure[m].t > TSTOP) continue;
	  }
	
	  // only reset Mcal for measures with a matching image
	  // do not reset Mcal for ubercal images unless explicitly requested

	  if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
	    if (!KEEP_UBERCAL) {
	      catalog[0].measure[m].Mcal = 0.0;
	      catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
	    } 
	  } else {
	    if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
	      catalog[0].measure[m].Mcal = 0.0;
	    }
	  }
	  catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS;
	}
      }
    }
  }

  // this sets flags in the measureT element, not the measure element
  setExclusions (catalog, 1, VERBOSE);  /* mark by area */

  /* set catalog[0].found[i] = FALSE */
  ALLOCATE (catalog[0].found_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
  ALLOCATE (catalog[0].foundWarp_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
  for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) {
    catalog[0].found_t[i] = FALSE;
    catalog[0].foundWarp_t[i] = FALSE;
  }

  for (i = 0; i < catalog[0].Naverage; i++) {
    catalog[0].average[i].psfQF     = NAN;	// force recalculation below
    catalog[0].average[i].psfQFperf = NAN;	// force recalculation below
    catalog[0].average[i].stargal   = NAN;	// force recalculation below
    catalog[0].average[i].photFlagsUpper = 0;	// reset (will be re-calculated here)
    catalog[0].average[i].photFlagsLower = 0;	// reset (will be re-calculated here)
    catalog[0].average[i].NwarpOK        = 0;	// reset (will be re-calculated here)
  }

  // XXX make this optional? (do not clean for -averages?)
  // XXX : for the moment, disable clean_measures
  if (!simpleAverage && FALSE) clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */

  ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure);
  setMeasureRank (catalog);
  setMrelOutput (catalog, 1, flatcorr); // sets the values secfilt.M = <measure.M + measure.Mflat - image.Mcal>
  setMcalOutput (catalog, 1, flatcorr); // sets measure.Mcal = image.Mcal

  /* clear ID_OBJ_POOR, ID_OBJ_FEW values before writing ??? */
  /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */
  // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found)
  for (i = 0; i < catalog[0].Naverage; i++) {
    catalog[0].average[i].flags &= ~ID_OBJ_FEW;
    catalog[0].average[i].flags &= ~ID_OBJ_POOR;
    off_t j;
    for (j = 0; j < Nsecfilt; j++) {
	catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_FEW; 
	catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_POOR; 
    }
    off_t m = catalog[0].average[i].measureOffset;
    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
      catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
    }
  }
}

/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
// For each star & photcode, we are selecting the measurements to use in the average
// photometry value.  This function is called several times, lowering the bar to
// acceptance on each pass. The rejections at the different passes are:

// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD)
// 1 : accept measurements thought to be POOR based on photflags
// 2 : accept the measurements marked as outliers; accept images marked as outliers
// 3 : accept measurements thought to be BAD based on photflags
// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT)

// XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all
// non-outliers with Imag range

// int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
// 
//   off_t k;
// 
//   int Nsecfilt = GetPhotcodeNsecfilt ();
// 
//   off_t m = average[0].measureOffset;
// 
//   for (k = 0; k < average[0].Nmeasure; k++, m++) {
//     fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
//   }
// 
//   int Ns;
//   for (Ns = 0; Ns < Nsecfilt; Ns++) {
//     fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
//   }
//   return 1;
// }

int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);

void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) {

  off_t i, k, m;
  off_t Ntot, Ntry, Nkeep, Nskip;
  float mag;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  Ntot = Ntry = Nskip = Nkeep = 0;

  /* allow measures from images marked POOR and FEW */
  if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
  
  /* allow measures marked as outliers (POOR) and off image region (AREA) */
  if (pass >= 3) {
    MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
  } else {
    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
  }

  /* mark measures which should be ignored on this pass */
  for (i = 0; i < catalog[0].Naverage; i++) {
    Ntot += catalog[0].average[i].Nmeasure;

    if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
      fprintf (stderr, "test obj\n");
      print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
    }

    // mark measurements for each secfilt separately
    int Ns;
    for (Ns = 0; Ns < Nphotcodes; Ns++) {
      
      int thisCode = photcodes[Ns][0].code;
      int Nsec = GetPhotcodeNsec(thisCode);
      
      /* star/photcodes already calibrated */
      if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue;  

      m = catalog[0].average[i].measureOffset;
      for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
	Ntry++;

	// skip measurements not related to this photcode
	PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
	if (!code) continue;
	if (code->equiv != thisCode) continue;

	// skip measurements by time range (mark as skipped for this secfilt)
	// XXX note that this is a bit dangerous : some objects may never get calibrated 
	if (TimeSelect) {
	  if (catalog[0].measure[m].t < TSTART) goto skip;
	  if (catalog[0].measure[m].t > TSTOP) goto skip;
	}
	
	// we now have a measurement relavant to this photcode (and time range)

	// clear SKIP for all measures at first 
	catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
	catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;

	// skip measurements from BAD images and mosaics (not REF mags)
	// do NOT skip measurements without a matching image (REF mags)
	off_t Nim = getImageEntry (m, 0);
	if (Nim > -1) {
	  if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
	  if (isnan(getMmos (m, 0))) goto skip;
	}

	// PASS 4 : skip measurements by inst mag limit (not REF mags)
	if ((pass < 4) && ImagSelect) {
	  if (Nim > -1) {
	    mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF);
	    if (mag < ImagMin) goto skip;
	    if (mag > ImagMax) goto skip;
	  }
	}

	// PASS 3 : accept bad measurements (eg, SAT, CR), internal outliers
	if (pass < 3) {
	  if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip;
	  if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) {
	    if (catalog[0].measure[m].psfQF < 0.85) goto skip;
	  }
	}
	
	// PASS 2 : (no additional cut)

	// PASS 1 : accept poor measurements as well (eg, POOR FIT, etc)
	if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip;

	// PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc)

	Nkeep ++;
	
	continue;

      skip:
	catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
	catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
	Nskip ++;
      }
    }
  }
  if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
			pass, Ntot, Ntry, Nskip, Nkeep);
}

void setMeasureRank (Catalog *catalog) {

  int i;

  // make these global
  int IMAGE_BAD  = ID_IMAGE_PHOTOM_NOCAL;
  int IMAGE_POOR = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
  int MEAS_BAD   = ID_MEAS_NOCAL        | ID_MEAS_SKIP_PHOTOM;
  int MEAS_POOR  = ID_MEAS_POOR_PHOTOM  | ID_MEAS_AREA;

  Measure     *measure     = catalog[0].measure;
  MeasureTiny *measureT    = catalog[0].measureT;
  char        *measureRank = catalog[0].measureRank;

  /* set measureRank[] based on various quality measurements */
  for (i = 0; i < catalog[0].Nmeasure; i++) {
    measureRank[i] = 11; // start at a low rank

    // clear SKIP for all measures at first 
    measureT[i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    measure [i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;

    // skip measurements without a valid photcode
    PhotCode *code = GetPhotcodebyCode (measure[i].photcode);
    if (!code) continue;

    // measurements outside time range have poor rank
    if (TimeSelect) {
      if (measure[i].t < TSTART) 		     { measureRank[i] = 10; continue; }
      if (measure[i].t > TSTOP)  		     { measureRank[i] = 10; continue; }
    }
    
    // measurements from BAD images and mosaics (not REF mags)
    off_t Nim = getImageEntry (i, 0);
    int imageFlags = getImageFlags (i, 0);
    int mosaicFlags = getMosaicFlags (i, 0);

    if (Nim > -1) {
      if (imageFlags  & IMAGE_BAD)                   { measureRank[i] = 9; continue; }

      if (mosaicFlags & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
    
      // measurements ranked by inst mag limit (not REF mags)
      if (ImagSelect) {
	float mag = PhotInst (&measure[i], MAG_CLASS_PSF);
	if (mag < ImagMin) 			     { measureRank[i] = 8; continue; }
	if (mag > ImagMax) 			     { measureRank[i] = 8; continue; }
      }
    }

    // RANK 1 : BAD photFlags
    if (measure[i].dbFlags & MEAS_BAD)               { measureRank[i] = 7; continue; }

    // RANK 3 : BAD photFlags (eg, SAT, CR), internal outliers
    if (measure[i].photFlags & code->photomBadMask)  { measureRank[i] = 6; continue; }

    if (Nim > -1) {
      if (imageFlags  & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
      if (mosaicFlags & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
    }

    // RANK 2 : psfQF value
    if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; }
	
    // RANK 1 : POOR photFlags
    if (measure[i].dbFlags & MEAS_POOR)              { measureRank[i] = 3; continue; }

    // RANK 1 : POOR photFlags
    if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; }
    
    // RANK 1 : psfQFperfect value
    if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; }
	
    // RANK 0 : perfect measurement:
    measureRank[i] = 0;
  }
}
