Changeset 41914

Timestamp:

Nov 11, 2021, 11:10:31 AM (5 years ago)

Author:

eugene

Message:

if too few sample cells are available, calculate the stdev of 3 subsamples

File:

• tags/ipp-ps1-20211104/psModules/src/astrom/pmAstrometryObjects.c (modified) (5 diffs)

tags/ipp-ps1-20211104/psModules/src/astrom/pmAstrometryObjects.c

@@ -41 +41 @@
 float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
 float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value);
+float pmAstromSubsetSystematics (psVector *value);
 
 #define PM_ASTROMETRYOBJECTS_DEBUG 1
@@ -442 +443 @@
     // VAL_COUNT (9) is the (min) goal density, but a single cell may have only MIN_COUNT (7)
     // we have N points.  require a min of 9 pts per cell (configurable?).  grid is square.
-    // Ncell*Ncell*9 = Npts, Ncell = MIN(sqrt(Npts/9), 5)
-
-    int Ncell = PS_MIN(sqrt((float)(xPos->n / VAL_COUNT)), 2);  // have to at least have a 2x2 grid
+    // Ncell*Ncell*9 = Npts, Ncell = MAX(sqrt(Npts/9), 5)
+
+    // XXX : ERROR -- use MAX not MIN
+    int Ncell = PS_MAX(sqrt((float)(xPos->n / VAL_COUNT)), 2);  // have to at least have a 2x2 grid
 
     // find the range of x,y values
-    float xMin = 1e9, xMax = -1e9, yMin = 1e9, yMax = 1e9;
+    float xMin = 1e9, xMax = -1e9, yMin = 1e9, yMax = -1e9;
     for (int i = 0; i < xPos->n; i++) {
         xMin = PS_MIN(xPos->data.F32[i],xMin);
@@ -479 +481 @@
         int ix = PS_MIN(PS_MAX(0, (xPos->data.F32[i] - xMin) * xStep), Ncell - 1);
         int iy = PS_MIN(PS_MAX(0, (yPos->data.F32[i] - yMin) * yStep), Ncell - 1);
-    
+
         psArray *yBin = xBin->data[ix];
         psVector *Bin = yBin->data[iy];
@@ -504 +506 @@
     }    
   
+    // The MIN_COUNT value means that sometimes fewer than 3 bins are any good.  
+    // For so few good bins, the stdev is not very meaningful.  In that case, let's
+    // just make 3 subgroups and ask about the stdev of those subgroups.
+    if (sample->n < 3) {
+      psFree (xBin);
+      psFree (stats);
+      psFree (sample);
+
+      float result = pmAstromSubsetSystematics (value);
+      return result;
+    }
+
     stats->options = PS_STAT_SAMPLE_STDEV;
     psVectorStats (stats, sample, NULL, NULL, 0); 
@@ -517 +531 @@
     psFree (sample);
   
+    return result;
+}
+
+# define SUBSET_NSAMPLE 3
+// last-ditch attempt to assess the systematic error in the data.  Just split into 3 bins,
+// calculate median in each, and calculate stdev of the subset.  The minimum number of
+// values needed to make this measurement in any sensible way : 3*2 = 6.
+float pmAstromSubsetSystematics (psVector *value) {
+
+    // give up if too few stars:
+    if (value->n < 2*SUBSET_NSAMPLE) return NAN;
+
+    psVector *sample = psVectorAlloc (SUBSET_NSAMPLE, PS_TYPE_F32);
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    // note that this drops the last 1 or 2 values if value->n % 3 = 1 or 2
+    int nSubset = value->n / SUBSET_NSAMPLE;
+    psVector *subset = psVectorAlloc (nSubset, PS_TYPE_F32);
+    for (int iter = 0; iter < SUBSET_NSAMPLE; iter++) {
+        for (int i = 0; i < nSubset; i++) {
+            subset->data.F32[i] = value->data.F32[i + nSubset*iter];
+        }  
+
+        // psVectorStats resets stats so we can call this repeatedly
+        psVectorStats (stats, subset, NULL, NULL, 0); 
+        if (isfinite(stats->sampleMedian)) {
+            sample->data.F32[iter] = stats->sampleMedian;
+        }
+    }
+
+    stats->options = PS_STAT_SAMPLE_STDEV;
+    psVectorStats (stats, sample, NULL, NULL, 0); 
+    float result = stats->sampleStdev;
+
+    if (!isfinite(stats->sampleStdev)) {
+      fprintf (stderr, "*** bad solution (2) ***\n");
+    }
+
+    psFree (stats);
+    psFree (sample);
+    psFree (subset); 
+
     return result;
 }