Changeset 39621

Timestamp:

Jul 1, 2016, 4:38:39 PM (10 years ago)

Author:

eugene

Message:

add clipped version of the mean stats for warp photometry

Location:

trunk/Ohana/src/relphot/src

Files:

• relphot_objects.c (modified) (1 diff)
• setMrelCatalog.c (modified) (9 diffs)

trunk/Ohana/src/relphot/src/relphot_objects.c

@@ -129 +129 @@
     initMrel (&catalog, 1);
     setMrelFinal (&catalog, NULL, TRUE);
+
+    // uncomment for extra verbosity
+    // MARKTIME("setMrelFinal for "OFF_T_FMT" average "OFF_T_FMT" measure : %f sec\n", catalog.Naverage, catalog.Nmeasure, dtime);
+    // gettimeofday (&startTimer, NULL); // reset timer
+
     // XXX if we want to have options for setting warp, chip, stack independently, we need to init only the desired ones
     // NOTE flatcorr == NULL, but it should have been applied already by setphot

trunk/Ohana/src/relphot/src/setMrelCatalog.c

@@ -1 +1 @@
 # include "relphot.h"
 
-/*
-# define TEST_OBJ_ID 0x3ae2
-# define TEST_CAT_ID 0x16c5f
-*/
+# if (0)
+# define TEST_OBJ_ID 0x0000000e
+# define TEST_CAT_ID 0x000076ee
+# else
 # define TEST_OBJ_ID 0
 # define TEST_CAT_ID 0
-
+# endif
+
+int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats);
 int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
 void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag);
+void sort_entry_by_offset (double *offset, int *entry, int N);
+void sort_StatDataSet (StatDataSet *dataset);
 
 # define UBERCAL_WEIGHT 100.0
-void sort_StatDataSet (StatDataSet *dataset);
 
 # define SKIP_THIS_MEAS(REASON) {                       \
@@ -903 +906 @@
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *psfstats = &results->psfstats;
-    int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats);
+    int Nranking = magStatsByRankingClipped (&results->psfData[Nsec], psfstats);
     if (Nranking) {
       secfilt[Nsec].FpsfWrp  = psfstats->mean;
@@ -916 +919 @@
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *apstats = &results->apstats;
-    Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
+    Nranking = magStatsByRankingClipped (&results->aperData[Nsec], apstats);
     if (Nranking) {
       secfilt[Nsec].FapWrp     = apstats->mean;
@@ -927 +930 @@
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *kronstats = &results->kronstats;
-    Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
+    Nranking = magStatsByRankingClipped (&results->kronData[Nsec], kronstats);
     if (Nranking) {
       secfilt[Nsec].FkronWrp     = kronstats->mean;
@@ -963 +966 @@
 }
 
-# if (0)
 // outlier warp measurements are driving bad mean values.  
 int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats) {
@@ -985 +987 @@
   // since we are only allowed to clip at most 25% of the points, we need to have at least 4 points to do this clipping
   if (Nranking > 3) {
+
+    double *fluxes = NULL;
+    double *offset = NULL;
+    int    *entry  = NULL;
+
+    // now find the distances
+    ALLOCATE (entry,  int,    Nranking);
+    ALLOCATE (offset, double, Nranking);
+    ALLOCATE (fluxes, double, Nranking);
+
     // first, measure the median of the flxlist:
-    ALLOCATE (fluxes, double, Nranking);
     for (i = 0; i < Nranking; i++) {
       fluxes[i] = dataset->flxlist[i];
@@ -993 +1004 @@
     int Nmidpoint = (int)(0.5*Nranking); // if Nranking is odd (e.g., 5), Nmidpoint is central bin (e.g., 2); if Nranking is even (e.g., 6), Nmidpoint is central bin + 1 (e.g., 3);
     double medianFlux = (Nranking % 2) ? fluxes[Nmidpoint] : 0.5*(fluxes[Nmidpoint] + fluxes[Nmidpoint - 1]);
+
+    for (i = 0; i < Nranking; i++) {
+      offset[i] = fabs((dataset->flxlist[i] - medianFlux) / dataset->errlist[i]);
+      entry[i] = i;
+    }
+
+    sort_entry_by_offset (offset, entry, Nranking);
+
+    // we are only going to clip at most 25% of the points.  Start at the 75% point and check for outliers
+    int Nreject = 0;
+    for (i = 0.75*Nranking; i < Nranking; i++) {
+      if (offset[i] > 4.0) {
+        // this is an outlier -- add 10 to the ranking
+        dataset->ranking[entry[i]] += 10;
+        Nreject ++;
+      }
+    }
+    Nranking -= Nreject; // the excluded points are now removed from the count
+
+    // re-sort by the rank so liststats below only operates on the unrejected points
+    sort_StatDataSet (dataset);
+
+    free (fluxes);
+    free (offset);
+    free (entry);
   }
 
@@ -1000 +1036 @@
   return (Nranking);
 }
-# endif
 
 int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags) {
@@ -1035 +1070 @@
 }
 
+void sort_entry_by_offset (double *offset, int *entry, int N) {
+
+# define SWAPFUNC(A,B){ double dtmp; int itmp; \
+    dtmp = offset[A]; offset[A] = offset[B]; offset[B] = dtmp;     \
+    itmp = entry[A];  entry[A]  = entry[B];  entry[B]  = itmp;     \
+  }
+# define COMPARE(A,B)(offset[A] < offset[B])
+  
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+  
+# undef SWAPFUNC
+# undef COMPARE
+}
+
 void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag) {