Changeset 42090

Timestamp:

Feb 28, 2022, 2:45:39 PM (4 years ago)

Author:

eugene

Message:

add a function to measure the systematic errors in cases with small numbers of stars

File:

• trunk/psModules/src/astrom/pmAstrometryObjects.c (modified) (3 diffs)

trunk/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
@@ -433 +434 @@
 
     // pre-filter the values to ensure no NANs, other invalid
-    if (xPos->n < VAL_COUNT*2*2) {
-      return NAN;
+
+    // if we do not have enough measurements (< 25), use pmAstromSubsetSystematics instead (OK to 9 values)
+    if (xPos->n < MIN_COUNT*2*2) {
+      float result = pmAstromSubsetSystematics (value);
+      return result;
     }
 
@@ -517 +521 @@
     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*3 = 6.
+float pmAstromSubsetSystematics (psVector *value) {
+
+    // give up if too few stars:
+    if (value->n < 3*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;
 }