Changeset 41813 for branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.c

Timestamp:

Sep 27, 2021, 6:31:19 AM (5 years ago)

Author:

eugene

Message:

add function to calculate stdev of median residuals per bin

File:

• branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.c (modified) (3 diffs)

branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.c

@@ -40 +40 @@
 // XXX this is defined in pmPSFtry.h, which makes no sense
 float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
+float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value);
 
 #define PM_ASTROMETRYOBJECTS_DEBUG 1
@@ -375 +376 @@
       results->dXrange = pmAstromVectorRange (xResGood, 0.1, 0.9, results->xStats->clippedStdev);
       results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev);
+
+      results->dXstdev = pmAstrom2DSystematics (xFit, yFit, xRes);
+      results->dYstdev = pmAstrom2DSystematics (xFit, yFit, yRes);
     }
 
@@ -395 +399 @@
 
     return (results);
+}
+
+# define VAL_COUNT 9.0
+# define MIN_COUNT 7.0
+
+float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value) {
+
+    // pre-filter the values to ensure no NANs, other invalid
+    if (xPos->n < VAL_COUNT*2*2) { return NAN; }
+
+    // generate a grid covering the full range of x,y and measure the median value in each
+    // grid cell.  calculate the stdev of those median values.
+
+    // 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
+
+    // find the range of x,y values
+    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);
+        xMax = PS_MAX(xPos->data.F32[i],xMax);
+        yMin = PS_MIN(yPos->data.F32[i],yMin);
+        yMax = PS_MAX(yPos->data.F32[i],yMax);
+    }
+
+    float xStep = Ncell / (xMax - xMin);
+    float yStep = Ncell / (yMax - yMin);
+  
+    if (isnan(xStep)) { return NAN; }
+    if (isnan(yStep)) { return NAN; }
+
+    psArray *xBin = psArrayAlloc (Ncell);
+    for (int ix = 0; ix < Ncell; ix++) {
+        psArray *yBin = psArrayAlloc (Ncell); 
+        xBin->data[ix] = yBin;
+        for (int iy = 0; iy < Ncell; iy++) {
+            yBin->data[iy] = psVectorAllocEmpty(128, PS_TYPE_F32);
+        }
+    }    
+
+    // xValue = ix/xStep + xMin -> ix = (xValue - xMin) * xStep;
+
+    for (int i = 0; i < xPos->n; i++) {
+        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];
+
+        psVectorAppend(Bin, value->data.F32[i]);
+    }
+
+    psVector *sample = psVectorAllocEmpty (Ncell*Ncell, PS_TYPE_F32);
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    // calculate the median for each vector and save on a vector
+    for (int ix = 0; ix < Ncell; ix++) {
+        psArray *yBin = xBin->data[ix];
+        for (int iy = 0; iy < Ncell; iy++) {
+            psVector *Bin = yBin->data[iy];
+      
+            // psVectorStats resets stats so we can call this repeatedly
+            psVectorStats (stats, Bin, NULL, NULL, 0); 
+            if (isfinite(stats->sampleMedian)) {
+                psVectorAppend (sample, stats->sampleMedian);
+            }
+        }
+    }    
+  
+    stats->options = PS_STAT_SAMPLE_STDEV;
+    psVectorStats (stats, sample, NULL, NULL, 0); 
+    float result = stats->sampleStdev;
+
+    // NOTE: the elements of xBin are freed automatically, which extends down to the vectors
+    psFree (xBin); 
+    psFree (stats);
+    psFree (sample);
+  
+    return result;
 }