Changeset 30550 for trunk/psModules/src/detrend/pmFringeStats.c

Timestamp:

Feb 9, 2011, 4:28:33 PM (15 years ago)

Author:

watersc1

Message:

Updated fringe scale solver. This looks like it eliminates the serious issues Heather noticed, and does not seem to increase the error in current good fits.

File:

• trunk/psModules/src/detrend/pmFringeStats.c (modified) (2 diffs)

trunk/psModules/src/detrend/pmFringeStats.c

@@ -1008 +1008 @@
         }
     }
+    // Allocate array of vectors to hold data.
+    psArray *bins = psArrayAlloc(4000);
+    for (int j = 0; j < bins->n; j++) {
+      bins->data[j] = psVectorAllocEmpty(1,PS_TYPE_F32);
+      //      psVector *v = psVectorAllocEmpty(1,PS_TYPE_F32);
+      //      bins = psArrayAdd(bins,1,v);
+    }
+
+    // Fill vectors
+    pmFringeStats *fringe = fringes->data[0];
+    for (int j = 0; j < numRegions; j++) {
+      if (mask->data.PS_TYPE_VECTOR_MASK_DATA[j] == 0) {
+        int array_bin = (int) ((fringe->f->data.F32[j] - -0.1) / 5e-5);
+        psVector *bin = bins->data[array_bin];
+        psVectorAppend(bin,science->f->data.F32[j]);
+      }
+    }
     
+    psVector *fringe_positions = psVectorAllocEmpty(4000,PS_TYPE_F32);
+    psVector *science_values   = psVectorAllocEmpty(4000,PS_TYPE_F32);
+    psVector *science_errors   = psVectorAllocEmpty(4000,PS_TYPE_F32);
+
+    psStats *binStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    for (int i = 0; i < 4000; i++) {
+      psVector *bin = bins->data[i];
+      if (bin->n > 2) {
+        psStatsInit(binStats);
+
+        psVectorStats(binStats,bin,NULL,NULL,1);
+        
+        if (isfinite(binStats->robustStdev) &&
+            isfinite(binStats->robustMedian) &&
+            binStats->robustStdev > 0) {
+          psVectorAppend(fringe_positions,-0.1 + i * 5e-5);
+          psVectorAppend(science_values, binStats->robustMedian);
+          psVectorAppend(science_errors, binStats->robustStdev);
+        }
+      }
+      psFree(bins->data[i]);
+    }
+    psFree(bins);
+    psFree(binStats);
+
+    for (int i = 0; i < fringe_positions->n; i++) {
+      psTrace("psModules.detrend",7,"FITDATA: %f %f %f\n",
+              fringe_positions->data.F32[i],
+              science_values->data.F32[i],
+              science_errors->data.F32[i]);
+    }
 /*     // Begin switch from old outlier removal and fitting code. */
 
     psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
 
-    pmFringeStats *fringe = fringes->data[0];
-    psVector *errors = psVectorAlloc(science->df->n,PS_TYPE_F32);
-    for (int j = 0; j < errors->n; j++) {
-      errors->data.F32[j] = 1 / science->df->data.F32[j];
-    }
-    psVectorFitPolynomial1D(poly,mask,0xff,science->f,errors,fringe->f);
-
+    //    pmFringeStats *fringe = fringes->data[0];
+/*     psVector *errors = psVectorAlloc(science->df->n,PS_TYPE_F32); */
+/*     for (int j = 0; j < errors->n; j++) { */
+/*       errors->data.F32[j] = 1 / science->df->data.F32[j]; */
+/*     } */
+/*     psVectorFitPolynomial1D(poly,mask,0xff,science->f,errors,fringe->f); */
+    psVectorFitPolynomial1D(poly,NULL,0xff,science_values,science_errors,fringe_positions);
+    psFree(fringe_positions);
+    psFree(science_values);
+    psFree(science_errors);
+    
     for (int i = 0; i <= poly->nX; i++) {
       scale->coeff->data.F32[i] = poly->coeff[i];
@@ -1027 +1079 @@
     psFree(poly);
     //    psFree(fringe);
-    psFree(errors);
+    //    psFree(errors);
 
     psFree(median);