Changeset 41557

Timestamp:

Apr 27, 2021, 4:25:16 PM (5 years ago)

Author:

eugene

Message:

common function to soften fit errors; use IRLS for MagResidSave; add Median and altSigma to MagResidSave

Location:

trunk/Ohana/src/relphot

Files:

• include/relphot.h (modified) (1 diff)
• src/ImageOps.c (modified) (1 diff)
• src/MagResidSave.c (modified) (14 diffs)
• src/MosaicOps.c (modified) (2 diffs)
• src/TGroupOps.c (modified) (1 diff)
• src/bcatalog.c (modified) (1 diff)
• src/fit1d_irls.c (modified) (1 diff)
• src/setMrelCatalog.c (modified) (1 diff)

trunk/Ohana/src/relphot/include/relphot.h

@@ -598 +598 @@
 void          FitDataSetAlloc     PROTO((FitDataSet *dataset, int Nmax, int order, int Nbootstrap));
 void          FitDataSetAddPriors PROTO((FitDataSet *dataset));
+StatType      FitDataSetSoften    PROTO((FitDataSet *dataset, int Nvalues));
 
 unsigned int *ReadTGroupFile      PROTO((FILE *f, int *nelem));

trunk/Ohana/src/relphot/src/ImageOps.c

@@ -784 +784 @@
     }
 
-    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
-    // XXX: apply airmass from above and fit only zp?
-
-    // XXX soften the stats?
+    // soften the errors based on the scatter
+    FitDataSetSoften (&psfStars, Nref);
 
     fit1d_irls (&psfStars, Nref);

trunk/Ohana/src/relphot/src/MagResidSave.c

@@ -12 +12 @@
   float dMresid;
   float  Mchisq;
+  float  Mmedian;
+  float  Maltsig;
   int    Nfit;
 } MagResidResult;
@@ -39 +41 @@
   }
 
-  // we are making a 0-order fit and not doing bootstrap analysis, using OLS
+  // we are making a 0-order fit and not doing bootstrap analysis
   FitDataSet psfStars;
   FitDataSetAlloc (&psfStars, Nmax, 0, 0);
-  psfStars.MaxIterations = 0;
+  // By default, we use IRLS; if MaxIterations = 0, then we use OLS
+  // psfStars.MaxIterations = 0;
 
   gfits_init_header (&header);
@@ -56 +59 @@
   gfits_define_bintable_column (&theader, "E", "MAG_ERROR", "magnitude error", "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "E", "MAG_CHISQ", "magnitude resid chisq", "unitless", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MAG_MEDIAN", "median magnitude residual", "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MAG_ALTSIG", "68 percentile range as sigma", "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "E", "AIRMASS",   "", "unitless", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "D", "MJD",       "", "unitless", 1.0, 0.0);
@@ -62 +67 @@
   gfits_define_bintable_column (&theader, "J", "FLAGS",     "image flags", "unitless",   1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "NFIT",      "number of stars fitted", NULL, 1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "J", "NMEAS",     "number of stars on image", NULL, 1.0, FT_BZERO_INT32);
 
   // generate the output array that carries the data
@@ -73 +79 @@
   ALLOCATE_PTR (dMresid,  float,  Nimage);
   ALLOCATE_PTR (Mchisq,   float,  Nimage);
+  ALLOCATE_PTR (Mmedian,  float,  Nimage);
+  ALLOCATE_PTR (Maltsig,  float,  Nimage);
   ALLOCATE_PTR (airmass,  float,  Nimage);
   ALLOCATE_PTR (mjd,      double, Nimage);
@@ -79 +87 @@
   ALLOCATE_PTR (flags,    int,    Nimage);
   ALLOCATE_PTR (Nfit,     int,    Nimage);
+  ALLOCATE_PTR (Nmeas,    int,    Nimage);
 
   // XXXX fill in the vectors here
@@ -94 +103 @@
     Mchisq[i]   = NAN;
     dMresid[i]  = NAN;
+    Mmedian[i]  = NAN;
+    Maltsig[i]  = NAN;
     Nfit[i]     = 0;
     Nmeas[i]    = N_onImage[i];
@@ -110 +121 @@
     Mchisq[i]   = result.Mchisq;
     dMresid[i]  = result.dMresid;
+    Mmedian[i]  = result.Mmedian;
+    Maltsig[i]  = result.Maltsig;
     Nfit[i]     = result.Nfit;
   }
 
   // add the columns to the output array
-  gfits_set_bintable_column (&theader, &ftable, "RA",        Rs,   Nimage);
-  gfits_set_bintable_column (&theader, &ftable, "DEC",       Ds,   Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "RA",        Rs,       Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "DEC",       Ds,       Nimage);
   gfits_set_bintable_column (&theader, &ftable, "MAG_RESID", Mresid,   Nimage);
   gfits_set_bintable_column (&theader, &ftable, "MAG_ERROR", dMresid,  Nimage);
   gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Mchisq,   Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "MAG_MEDIAN", Mmedian, Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "MAG_ALTSIG", Maltsig, Nimage);
   gfits_set_bintable_column (&theader, &ftable, "AIRMASS",   airmass,  Nimage);
   gfits_set_bintable_column (&theader, &ftable, "MJD",       mjd,      Nimage);
@@ -125 +140 @@
   gfits_set_bintable_column (&theader, &ftable, "FLAGS",     flags,    Nimage);
   gfits_set_bintable_column (&theader, &ftable, "NFIT",      Nfit,     Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "NMEAS",     Nmeas,    Nimage);
 
   free (Rs);
@@ -131 +147 @@
   free (dMresid );
   free (Mchisq  );
+  free (Mmedian );
+  free (Maltsig );
   free (airmass );
   free (mjd     );
@@ -137 +155 @@
   free (flags   );
   free (Nfit    );
+  free (Nmeas   );
 
   FILE *f = fopen (filename, "w");
@@ -185 +204 @@
   result.dMresid = NAN;
   result.Mchisq  = NAN;
+  result.Mmedian = NAN;
+  result.Maltsig = NAN;
   result.Nfit    = 0;
 
@@ -239 +260 @@
     result.dMresid = NAN;
     result.Mchisq  = NAN;
+    result.Mmedian = NAN;
+    result.Maltsig = NAN;
     result.Nfit    = Nref;
     return result;
   }
+
+  StatType stats = FitDataSetSoften (psfStars, Nref);
+  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;
 
   // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
@@ -249 +275 @@
   result.Nfit    = psfStars->Nmeas;
   result.Mchisq  = psfStars->chisq;
+  result.Mmedian = stats.median;
+  result.Maltsig = altSigma;
 
   return result;

trunk/Ohana/src/relphot/src/MosaicOps.c

@@ -1161 +1161 @@
 
   // in GOOD_MOSAIC mode, we fit good mosaics ignoring night state
+  // This case should only be used when not fitting nights / tgroups 
   if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_GOOD_MOSAIC)) {
     if (badMosaic) useMmos = FALSE; // do not fit bad mosaics
@@ -1312 +1313 @@
   }
 
-  StatType stats;
-  liststats_setmode (&stats, "MEDIAN");
-  liststats (psfStars->alldata->yVector, NULL, NULL, N, &stats);
-  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
-
-  // soften the individual errors with 10% of the scatter above
-  for (j = 0; j < N; j++) {
-    double newSigma = hypot(psfStars->alldata->dyVector[j], 0.1*altSigma);
-    psfStars->alldata->dyVector[j] = newSigma;
-  }
-
+  // soften the errors based on the scatter
+  FitDataSetSoften (psfStars, N);
+  
   fit1d_irls (psfStars, N);
   fit1d_irls (brightStars, Nbright);

trunk/Ohana/src/relphot/src/TGroupOps.c

@@ -927 +927 @@
   */
 
-  // use liststats to find the 20-pct, median, 80-pct points
-  StatType stats;
-  liststats_setmode (&stats, "MEDIAN");
-  liststats (psfStars->alldata->yVector, NULL, NULL, N, &stats);
-  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
+  StatType stats = FitDataSetSoften (psfStars, N);
+  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;
 
   TGTimes *mygroup = (TGTimes *) myTGroup->parent;
   fprintf (stderr, "TGroup Stats %f,%d : %d %d %6.3f %6.3f\n", ohana_sec_to_mjd(mygroup->start), myTGroup->photcode, (int) N, (int) myTGroup->Nmeasure, stats.median, altSigma);
-
-  // soften the individual errors with 10% of the scatter above
-  for (j = 0; j < N; j++) {
-    double newSigma = hypot(psfStars->alldata->dyVector[j], 0.1*altSigma);
-    psfStars->alldata->dyVector[j] = newSigma;
-  }
 
   // do anything special with identified good nights?

trunk/Ohana/src/relphot/src/bcatalog.c

@@ -106 +106 @@
 
       // skip garbage measurements
-      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF < 0.85))     { Npsfqf ++; continue; }
-      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.85)) { Npsfqf ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF < 0.95))     { Npsfqf ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.95)) { Npsfqf ++; continue; }
       if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
 

trunk/Ohana/src/relphot/src/fit1d_irls.c

@@ -368 +368 @@
   }
 }
+
+StatType FitDataSetSoften (FitDataSet *dataset, int Nvalues) {
+
+  // use liststats to find the 20-pct, median, 80-pct points
+  StatType stats;
+  liststats_setmode (&stats, "MEDIAN");
+  liststats (dataset->alldata->yVector, NULL, NULL, Nvalues, &stats);
+
+  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
+
+  // soften the individual errors with 10% of the scatter above
+  for (int j = 0; j < Nvalues; j++) {
+    double newSigma = hypot(dataset->alldata->dyVector[j], 0.1*altSigma);
+    dataset->alldata->dyVector[j] = newSigma;
+  }
+  
+  return stats;
+}

trunk/Ohana/src/relphot/src/setMrelCatalog.c

@@ -267 +267 @@
     if (USE_REFERENCE_WEIGHT && (code->type == PHOT_REF)) modifiedWeight = code->photomErrSys; // we are overloading this field for now
 
+    // XXX I should probably do something more clever with the measured scatter in the nights, mosaics, images
+
     // combine the various errors in quadrature here:
     float dMsys = hypot (code->photomErrSys, dMgrp);