Changeset 39648

Timestamp:

Aug 5, 2016, 1:08:51 PM (10 years ago)

Author:

eugene

Message:

make the initial guess for IRLS stats use the median

Location:

trunk/Ohana/src/relphot

Files:

• include/relphot.h (modified) (2 diffs)
• src/StatDataSetOps.c (modified) (2 diffs)
• src/liststats.c (modified) (8 diffs)
• src/select_images.c (modified) (1 diff)

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

@@ -10 +10 @@
 
 # define ID_SECF_STACK_PRIMARY 0x00004000
+
+# define NBOOTSTRAP 100
 
 # ifndef MAX_INT
@@ -121 +123 @@
   int    *msklist;            // mask modifications
   int     Nlist;
+
+  double *values;
+  double *wtvals;
+  double *wtlist;
+  double *ykeep;
+  double *dykeep;
+  double *wtkeep;
+  double *ysample;
+  double *dysample;
+  double *wtsample;
+  double *bvalue;
+
 } StatDataSet;
 

trunk/Ohana/src/relphot/src/StatDataSetOps.c

@@ -14 +14 @@
     ALLOCATE (dataset[i].measSeq,    int, Nmax);
     ALLOCATE (dataset[i].msklist,    int, Nmax);
+
+    ALLOCATE (dataset[i].values,   double, Nmax);
+    ALLOCATE (dataset[i].wtvals,   double, Nmax);
+    ALLOCATE (dataset[i].wtlist,   double, Nmax);
+    ALLOCATE (dataset[i].ykeep,    double, Nmax);
+    ALLOCATE (dataset[i].dykeep,   double, Nmax);
+    ALLOCATE (dataset[i].wtkeep,   double, Nmax);
+    ALLOCATE (dataset[i].ysample,  double, Nmax);
+    ALLOCATE (dataset[i].dysample, double, Nmax);
+    ALLOCATE (dataset[i].wtsample, double, Nmax);
+    ALLOCATE (dataset[i].bvalue,   double, NBOOTSTRAP);
   }  
+
   return dataset;
 }
@@ -29 +41 @@
     FREE (dataset[i].measSeq);
     FREE (dataset[i].msklist);
+
+    FREE (dataset[i].values);
+    FREE (dataset[i].wtvals);
+    FREE (dataset[i].wtlist);
+    FREE (dataset[i].ykeep);
+    FREE (dataset[i].dykeep);
+    FREE (dataset[i].wtkeep);
+    FREE (dataset[i].ysample);
+    FREE (dataset[i].dysample);
+    FREE (dataset[i].wtsample);
+    FREE (dataset[i].bvalue);
   }  
   FREE (dataset);

trunk/Ohana/src/relphot/src/liststats.c

@@ -193 +193 @@
 # define FLT_TOLERANCE 1e-6
 # define WEIGHT_THRESHOLD 0.3
-# define NBOOTSTRAP 100
 
 int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts);
@@ -202 +201 @@
 int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats) {
 
-  double value;
-
   liststats_init (stats);
 
-  // OLS
-  if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, NULL, Npoints)) return FALSE;
+  if (Npoints == 0) {
+    double value = NAN;
+    stats->mean = value;
+    stats->min  = value;
+    stats->max  = value;
+    stats->Nmeas = Npoints;
+    stats->chisq = NAN;
+    stats->sigma = NAN;
+    stats->error = NAN;
+    return TRUE;
+  }
+
+  if (Npoints == 1) {
+    double value = dataset->flxlist[0];
+    stats->mean = value;
+    stats->min  = value;
+    stats->max  = value;
+    stats->Nmeas = Npoints;
+    stats->chisq = NAN;
+    stats->sigma = NAN;
+    stats->error = dataset->errlist[0];
+    return TRUE;
+  }
+
+  int midpt = 0.5 * Npoints;
+
+  // make the initial guess based on the median (not weighted mean)
+  // ALLOCATE_PTR (values, double, Npoints);
+  for (int i = 0; i < Npoints; i++) {
+    dataset->values[i] = dataset->flxlist[i];
+  }
+  dsort (dataset->values, Npoints);
+  double value = (Npoints % 2) ? dataset->values[midpt] : 0.5*(dataset->values[midpt] + dataset->values[midpt-1]);
+
+  // OLS (replace by median above)
+  // if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, NULL, Npoints)) return FALSE;
   
   // XXX add to dataset elements?
-  ALLOCATE_PTR (wtvals, double, Npoints);
-  ALLOCATE_PTR (wtlist, double, Npoints);
+  // ALLOCATE_PTR (wtvals, double, Npoints);
+  // ALLOCATE_PTR (wtlist, double, Npoints);
 
   int converged = FALSE;
@@ -218 +249 @@
     for (int i = 0; i < Npoints; i++) {
       // we are only including the formal error, not the weight in the definition of wt[]
-      wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+      dataset->wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
     }
     
     double oldValue = value;
-    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, wtvals, Npoints)) {
+    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->wtvals, Npoints)) {
       value = oldValue;
       break;
@@ -238 +269 @@
   // double Sum_W = 0.0;
   for (int i = 0; i < Npoints; i++) {
-    wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
-    wtlist[i] = wtvals[i];
-    // Sum_W += wtvals[i];
-  }
-  dsort (wtlist, Npoints);
-  int midpt = 0.5 * Npoints;
-  double WtMedian = (Npoints % 2) ? wtlist[midpt] : 0.5*(wtlist[midpt] + wtlist[midpt-1]);
+    dataset->wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+    dataset->wtlist[i] = dataset->wtvals[i];
+    // Sum_W += dataset->wtvals[i];
+  }
+  dsort (dataset->wtlist, Npoints);
+  double WtMedian = (Npoints % 2) ? dataset->wtlist[midpt] : 0.5*(dataset->wtlist[midpt] + dataset->wtlist[midpt-1]);
   // double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
   double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
@@ -250 +280 @@
   // generate the unmasked subset
   // XXX add these to the dataset elements?
-  ALLOCATE_PTR ( ykeep, double, Npoints);
-  ALLOCATE_PTR (dykeep, double, Npoints);
-  ALLOCATE_PTR (wtkeep, double, Npoints);
+  // ALLOCATE_PTR ( ykeep, double, Npoints);
+  // ALLOCATE_PTR (dykeep, double, Npoints);
+  // ALLOCATE_PTR (wtkeep, double, Npoints);
 
   // save unmasked points
@@ -261 +291 @@
   double dSig = 0.0;
   for (int i = 0; i < Npoints; i++) {
-    if ((wtvals[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
+    if ((dataset->wtvals[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
       dataset->msklist[i] = TRUE; // mark the masked points
       continue;
     }
-     ykeep[Nkeep] = dataset->flxlist[i];
-    dykeep[Nkeep] = dataset->errlist[i];
-    wtkeep[Nkeep] = dataset->wgtlist[i]; // externally-supplied weight
+     dataset->ykeep[Nkeep] = dataset->flxlist[i];
+    dataset->dykeep[Nkeep] = dataset->errlist[i];
+    dataset->wtkeep[Nkeep] = dataset->wgtlist[i]; // externally-supplied weight
     Nkeep ++;
     
@@ -283 +313 @@
   // bootstrap resampling to generate the errorbars
   // XXX add these to the dataset elements?
-  ALLOCATE_PTR (ysample,  double, Nkeep);
-  ALLOCATE_PTR (dysample, double, Nkeep);
-  ALLOCATE_PTR (wtsample, double, Nkeep);
-  ALLOCATE_PTR (bvalue,   double, NBOOTSTRAP); // vector to save the bootstrap values
+  // ALLOCATE_PTR (ysample,  double, Nkeep);
+  // ALLOCATE_PTR (dysample, double, Nkeep);
+  // ALLOCATE_PTR (wtsample, double, Nkeep);
+  // ALLOCATE_PTR (bvalue,   double, NBOOTSTRAP); // vector to save the bootstrap values
 
   int Nboot = 0;
@@ -295 +325 @@
       // I need to draw Npoints random entries from 'points' with replacement:
       int N = Nkeep * drand48();
-       ysample[i] =  ykeep[N];
-      dysample[i] = dykeep[N];
-      wtsample[i] = wtkeep[N];
-    }
-
-    if (!fit_least_squares (&value, ysample, dysample, wtsample, NULL, Nkeep)) continue;
-
-    bvalue[Nboot] = value;
+      dataset->ysample[i]  = dataset->ykeep[N];
+      dataset->dysample[i] = dataset->dykeep[N];
+      dataset->wtsample[i] = dataset->wtkeep[N];
+    }
+
+    if (!fit_least_squares (&value, dataset->ysample, dataset->dysample, dataset->wtsample, NULL, Nkeep)) continue;
+
+    dataset->bvalue[Nboot] = value;
     Nboot ++;
   }
 
-  dsort (bvalue, Nboot);
+  dsort (dataset->bvalue, Nboot);
   
-  double Slo = VectorFractionInterpolate (bvalue, 0.158655, Nboot);
-  double Shi = VectorFractionInterpolate (bvalue, 0.841345, Nboot);
+  double Slo = VectorFractionInterpolate (dataset->bvalue, 0.158655, Nboot);
+  double Shi = VectorFractionInterpolate (dataset->bvalue, 0.841345, Nboot);
   stats->error = (Shi - Slo) / 2.0;
-
-  free (bvalue);
-  free ( ysample);
-  free (dysample);
-  free (wtsample);
-  free ( ykeep);
-  free (dykeep);
-  free (wtkeep);
-  free (wtvals);
-  free (wtlist);
 
   return TRUE;

trunk/Ohana/src/relphot/src/select_images.c

@@ -46 +46 @@
   double DmaxSkyRegion = region[0].Dmax;
 
-  double dD = CALIBRATE_STACKS_AND_WARPS ? 0.5 : 5.0;
+  double dD = CALIBRATE_STACKS_AND_WARPS ? 0.0 : 5.0;
   double dR = dD / cos(RAD_DEG*DminSkyRegion);
   RminSkyRegion -= dR;