Changeset 42788 for trunk/Ohana/src/relphot

Timestamp:

Feb 7, 2025, 3:07:30 PM (17 months ago)

Author:

eugene

Message:

modify grid correction to use median of the data vectors

Location:

trunk/Ohana/src/relphot

Files:

• Makefile (modified) (1 diff)
• include/relphot.h (modified) (3 diffs)
• src/DumpAllMags.c (added)
• src/GridIO.c (modified) (1 diff)
• src/GridOps.c (modified) (9 diffs)
• src/relphot_images.c (modified) (1 diff)

trunk/Ohana/src/relphot/Makefile

@@ -82 +82 @@
 $(SRC)/setExclusions.$(ARCH).o   \
 $(SRC)/setMrelFinal.$(ARCH).o    \
-$(SRC)/setMrelCatalog.$(ARCH).o          \
-$(SRC)/BoundaryTreeOps.$(ARCH).o         \
+$(SRC)/setMrelCatalog.$(ARCH).o   \
+$(SRC)/BoundaryTreeOps.$(ARCH).o  \
+$(SRC)/DumpAllMags.$(ARCH).o      \
 $(SRC)/write_coords.$(ARCH).o
 

trunk/Ohana/src/relphot/include/relphot.h

@@ -160 +160 @@
   float dY; // bin = int(Ychip * dY), dY = Ny / NyChip
   // NxChip, NyChip = 4900,4900 for now
+
+    int **nAlloc; // allocated vector length
+  double ***dMval; // values used to calculate corrections
 } GridCorrectionType;
 
@@ -607 +610 @@
 void          initGridBins        PROTO((void));
 GridCorrectionType *getGridCorrByCode PROTO((int code));
+GridCorrectionType *newGridCorrByCode PROTO((int code));
 GridCorrectionType *getGridCorrNext   PROTO((int *Nlast));
 int           GridCorrectionSave  PROTO((void));
@@ -894 +898 @@
 double get_median_zpt_tgroups (short photcode);
 void   set_median_zpt_tgroups (short photcode, double zpt);
+
+int DumpAllMags(char *filename, Catalog *catalog, int Ncatalog);

trunk/Ohana/src/relphot/src/GridIO.c

@@ -151 +151 @@
 
     if (!gfits_scan (&header, "PHOTCODE", "%d", 1, &photcode))    myAbort ("failed to get PHOTCODE");
-    GridCorrectionType *GridCorr = getGridCorrByCode (photcode);
+    GridCorrectionType *GridCorr = newGridCorrByCode (photcode);
     GridCorr->photcode = photcode;
 

trunk/Ohana/src/relphot/src/GridOps.c

@@ -30 +30 @@
 # define NX_CHIP_DEFAULT 1000
 # define NY_CHIP_DEFAULT 1000
-# define NX_BIN_DEFAULT 2
-# define NY_BIN_DEFAULT 2
+# define NX_BIN_DEFAULT 4
+# define NY_BIN_DEFAULT 4
 
 void initGridBins (void) {
@@ -76 +76 @@
       if (GridCorr[code]->dMgrid) FREE (GridCorr[code]->dMgrid[ix]);
       if (GridCorr[code]->nMgrid) FREE (GridCorr[code]->nMgrid[ix]);
+      if (GridCorr[code]->nAlloc) FREE (GridCorr[code]->nAlloc[ix]);
+
+      // dMval is a 2D array of vectors
+      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
+        FREE (GridCorr[code]->dMval[ix][iy]);
+      }
+      FREE (GridCorr[code]->dMval[ix]);
     }      
 
@@ -81 +88 @@
     FREE (GridCorr[code]->dMgrid);
     FREE (GridCorr[code]->nMgrid);
+    FREE (GridCorr[code]->nAlloc);
+    FREE (GridCorr[code]-> dMval);
 
     FREE(GridCorr[code]);
@@ -109 +118 @@
 }
 
-GridCorrectionType *getGridCorrByCode (int code) {
+// generate a new GridCorr structure in the array
+GridCorrectionType *newGridCorrByCode (int code) {
 
   if (GridCorr == NULL) return NULL;
@@ -116 +126 @@
   if (code <          0) return NULL;
 
+  myAssert (GridCorr[code] == NULL, "oops, should not already be allocated!");
   ALLOCATE(GridCorr[code], GridCorrectionType, 1)
+  GridCorrectionType *result = GridCorr[code];
+  return result;
+}
+
+// return the GridCode structure corresponding to 'code'
+GridCorrectionType *getGridCorrByCode (int code) {
+
+  if (GridCorr == NULL) return NULL;
+
+  if (code >= NGridCorr) return NULL;
+  if (code <          0) return NULL;
+
   GridCorrectionType *result = GridCorr[code];
   return result;
@@ -163 +186 @@
     // we are normally accessing this array randomly, so there is no advantage to
     // doing Mgrid[y][x] vs Mgrid[x][y]
+    ALLOCATE (GridCorr[code]-> dMval, double **, NX_BIN);
     ALLOCATE (GridCorr[code]-> Mgrid, float *, NX_BIN);
     ALLOCATE (GridCorr[code]->dMgrid, float *, NX_BIN);
     ALLOCATE (GridCorr[code]->nMgrid,   int *, NX_BIN);
+    ALLOCATE (GridCorr[code]->nAlloc,   int *, NX_BIN);
     for (int ix = 0; ix < NX_BIN; ix++) {
       ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NY_BIN);
       ALLOCATE (GridCorr[code]->dMgrid[ix], float, NY_BIN);
       ALLOCATE (GridCorr[code]->nMgrid[ix],   int, NY_BIN);
+      ALLOCATE (GridCorr[code]->nAlloc[ix],   int, NY_BIN);
+
+      // when we first allocate this array, set the elements to 0 (NULL)
+      // these arrays will be allocated or reallocated in resetMgrid
+      ALLOCATE_ZERO (GridCorr[code]-> dMval[ix], double *, NX_BIN);
     }      
   }
@@ -190 +220 @@
         GridCorr[code]->dMgrid[ix][iy] = 0.0;
         GridCorr[code]->nMgrid[ix][iy] =   0;
+        GridCorr[code]->nAlloc[ix][iy] = 100;
+
+        // allocate or reset vector length to default
+        if (GridCorr[code]->dMval[ix][iy]) {
+          REALLOCATE (GridCorr[code]->dMval[ix][iy], double, GridCorr[code]->nAlloc[ix][iy]);
+        } else {
+          ALLOCATE (GridCorr[code]->dMval[ix][iy], double, GridCorr[code]->nAlloc[ix][iy]);
+        }
       }
     }      
@@ -263 +301 @@
         if (fabs(dM) > 0.5) continue;
 
-        grid-> Mgrid[ix][iy] += dM;
-        grid->dMgrid[ix][iy] += dM*dM;
+        int ival = grid->nMgrid[ix][iy];
+        grid-> dMval[ix][iy][ival] = dM;
         grid->nMgrid[ix][iy] ++;
+        CHECK_REALLOCATE (grid-> dMval[ix][iy], double, grid->nAlloc[ix][iy], grid->nMgrid[ix][iy], 100);
+
+        /* XXX old code to calculate simple means:
+           grid-> Mgrid[ix][iy] += dM;
+           grid->dMgrid[ix][iy] += dM*dM;
+           grid->nMgrid[ix][iy] ++;
+        */
       }
     }
   }
+
+  StatType stats;
+  liststats_setmode (&stats, "MEDIAN");
+
+  // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt
+  for (int code = 0; code < NGridCorr; code++) {
+    if (!GridCorr[code]) continue;
+
+    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
+      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
+
+        // cells without sufficient coverage stay at 0.0
+        if (GridCorr[code]->nMgrid[ix][iy] < 5) {
+          GridCorr[code]-> Mgrid[ix][iy] = 0.0;
+          GridCorr[code]->dMgrid[ix][iy] = 0.0;
+          continue;
+        }
+
+        liststats_init (&stats);
+
+        liststats (GridCorr[code]->dMval[ix][iy], NULL, NULL, GridCorr[code]->nMgrid[ix][iy], &stats);
+        double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
+
+        GridCorr[code]-> Mgrid[ix][iy] = stats.median;
+        GridCorr[code]->dMgrid[ix][iy] = altSigma; // sample stdev
+      }
+    }
+  }
+  return;
+}
+
+# if (0)
 
   // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt
@@ -300 +377 @@
     // float GridAve = GridSum / GridCnt;
     // fprintf (stderr, "grid average: %f\n", GridAve);
-  }
-  return;
-}
+# endif
 
 float getMgrid (Measure *measure) {

trunk/Ohana/src/relphot/src/relphot_images.c

@@ -147 +147 @@
     MagResidSave ("mag.resid.fits", catalog);
 
+    // DumpAllMags ("all.mags.fits", catalog, Ncatalog);
+
     save_images_updates (&db);