Changeset 34317

Timestamp:

Aug 16, 2012, 2:38:37 PM (14 years ago)

Author:

bills

Message:

Several changes to psphot and psphotStack

1. Add function to compute the memory used by all of the sources
2. Add function to log memory stats in psphotStack
3. Don't load or psf match the convolved images into psphotStack unless they are needed.

This saves a lot of time if RADIAL_APERTURES are not enabled and we are measuring on the raw images

4. Several improvements to the error handling in the Kron measurements. In particular the Mrf value

was being treated as positive when both the numerator and denomiator were negative. This is the cause
of many of the very large kron radius measurements. Set Mrf to NAN in these cases.

5. Skip kron measurments for source if the previous Mrf measurement was NAN earlier in the process.

This saves memory by not expanding the source images when the kron measurements were actually going
to be skipped

6. Remove the hack cut on fwhmMaj. Managing the kron radius measurements better fixes the memory explosion

that that helped to reduce.

Location:

trunk/psphot/src

Files:

• Makefile.am (modified) (2 diffs)
• psphot.h (modified) (3 diffs)
• psphotKronIterate.c (modified) (13 diffs)
• psphotReadout.c (modified) (2 diffs)
• psphotReadoutMinimal.c (modified) (1 diff)
• psphotSetThreads.c (modified) (1 diff)
• psphotSourceMemory.c (added)
• psphotStackAllocateOutput.c (added)
• psphotStackImageLoop.c (modified) (3 diffs)
• psphotStackMatchPSFsNext.c (modified) (2 diffs)
• psphotStackReadout.c (modified) (12 diffs)

trunk/psphot/src/Makefile.am

@@ -117 +117 @@
         psphotStackObjects.c          \
         psphotStackPSF.c              \
+        psphotStackAllocateOutput.c   \
         psphotCleanup.c
 
@@ -219 +220 @@
         psphotPetrosianVisual.c        \
         psphotEfficiency.c             \
-        psphotSetNFrames.c
+        psphotSetNFrames.c             \
+        psphotSourceMemory.c
 
 # not currently used

trunk/psphot/src/psphot.h

@@ -360 +360 @@
                                      pmReadout **chiReadout,
                                      int index);
+bool psphotStackAllocateOutput( const pmConfig *config, pmFPAview *view, psMetadata *recipe);
 
 bool psphotStackRemoveChisqFromInputs (pmConfig *config, const char *filerule);
@@ -483 +484 @@
 bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal);
 
-bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index);
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass);
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass);
 bool psphotKronIterate_Threaded (psThreadJob *job);
 
@@ -509 +510 @@
     );
 
+bool psphotSourceMemory(pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index);
+
 #endif

trunk/psphot/src/psphotKronIterate.c

@@ -5 +5 @@
 
 
-bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
 {
     bool status = true;
@@ -42 +42 @@
         // psAssert (psf, "missing psf?");
 
-        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) {
+        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
             return false;
@@ -54 +54 @@
 bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
 
-bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
+#ifdef DUMP_KRS
+FILE *dumpFile = NULL;
+#endif
+
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
 
     bool status = false;
+
+#ifdef DUMP_KRS
+    if (!dumpFile) {
+        dumpFile = fopen("kr.txt", "w");
+        psAssert (dumpFile, "failed to open kr.txt");
+    }
+    fprintf(dumpFile, "\n\n Input %d\n", index);
+#endif
 
     if (!sources->n) {
@@ -64 +76 @@
 
     psTimerStart ("psphot.kron");
-
 
     // determine the number of allowed threads
@@ -236 +247 @@
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
 
 // set this to 0 to run without threading
@@ -298 +310 @@
     float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
     float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
+    int pass                        = PS_SCALAR_VALUE(job->args->data[13],S32);
+#ifndef REVERT_ON_BAD_MEASUREMENT
+    (void) pass;
+#endif
 
     for (int j = 0; j < KRON_ITERATIONS; j++) {
@@ -309 +325 @@
             if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
             if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+            if (!isfinite(source->moments->Mrf)) {
+                // Once we save a bad Mrf measurement we give up on this source
+                // checking here allows us to avoid adding and subtracting the model
+                continue;
+            }
 
             // replace object in image
@@ -322 +343 @@
             }
 
-            // On first iteration set window radius to sky radius (if valid) on second iteration
-            // use a factor times the previous radial moment value up to a maximum value that
-            // depends on the surface brightness of the source
-            float maxWindow;
-            if (j == 0) {
-                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
-            } else {
+            // On first iteration set window radius to sky radius (if valid). We also use this on subsequent
+            // iterations if we cannot find a better limit
+            float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            if (j > 0) {
+                // on subsequent iterations we use a factor times the previous radial moment value
+                // limited to a maximum value that depends on the surface brightness of the source
                 if (KRON_SB_MIN_DIVISOR) {
                     // Limit window radius based on surface brightness if we have a good measurement of kron flux
-                    if (isfinite(source->moments->Mrf) && source->moments->Mrf > 0 && 
-                        isfinite(source->moments->KronFlux)  && (source->moments->KronFlux > 0)) {
+                    if (isfinite(source->moments->KronFlux)  && (source->moments->KronFlux > 0)) {
                         float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR;
 
-                        maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
-                    } else {
-                        maxWindow = RADIUS;
+                        if (isfinite(source->moments->Mrf) && source->moments->Mrf > 0) {
+                            maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
+                        } else {
+                            maxWindow = PS_MIN(Rmax, maxWindow);
+                        }
                     }
                 } else {
@@ -345 +366 @@
             }
             float windowRadius = PS_MAX(RADIUS, maxWindow);
+
+#ifdef REVERT_ON_BAD_MEASURMENT
+            // save previous measurements. We might revert back to them if this round fails
+            float MrfPrior = source->moments->Mrf;
+            float KronFluxPrior = source->moments->KronFlux;
+            float KronFluxErrPrior = source->moments->KronFluxErr;
+#endif
 
             // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
@@ -371 +399 @@
             }
 
+#ifdef REVERT_ON_BAD_MEASUREMENT
+            // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously
+            // in pass 1 keep that measurement
+            bool reverted = false;
+            if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) {
+                source->moments->Mrf = MrfPrior;    // This is finite otherwise we wouldn't have gotten here
+                source->moments->KronFlux = KronFluxPrior;
+                source->moments->KronFluxErr = KronFluxErrPrior;
+                reverted = true;
+            }
+#endif
+#ifdef DUMP_KRS
+#ifndef REVERT_ON_BAD_MEASUREMENT
+            bool reverted = false;
+#endif
+            fprintf(dumpFile, "%7d %1d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, reverted, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
+#endif
+
             // if we subtracted it above, re-subtract the object, leave local sky
             if (reSubtract) {
@@ -490 +536 @@
 
     float MrfTry = RF/RS;
-    if (!isfinite(MrfTry)) {
-        // We did not get a successul measurement of the kron radius. 
-        // Leave the current values unchanged.
+    if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
+        // We did not get a good measurement
+        source->moments->Mrf = NAN;
+        source->moments->KronFlux  = NAN;
+        source->moments->KronFluxErr  = NAN;
         return false;
     }
@@ -546 +594 @@
     }
 
-    source->moments->Mrf = Mrf;
+    source->moments->Mrf         = Mrf;
     source->moments->KronFlux    = Sum;
     source->moments->KronFluxErr = sqrt(Var);
@@ -564 +612 @@
     psAssert(kronWindow, "need a window");
 
-    // If we are not applying the window then we don't need to check for valid Mrf here.
+    // XXX: If we are not applying the window then we don't need to check for valid Mrf here.
     // We should give the this module a chance to measure a good value. 
-    // Not yet though. We need to test the effect of this
+    // However experiments show that it hardly ever succeeds in getting a better value
     if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
 

trunk/psphot/src/psphotReadout.c

@@ -194 +194 @@
     // but this is chosen above to be appropriate for the PSF objects (not galaxies)
     // psphotKronMasked(config, view, filerule);
-    psphotKronIterate(config, view, filerule);
+    psphotKronIterate(config, view, filerule, 1);
 
     // identify CRs and extended sources (only unmeasured sources are measured)
@@ -312 +312 @@
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, filerule);
-    psphotKronIterate(config, view, filerule);
+    psphotKronIterate(config, view, filerule, 2);
 
     // measure source size for the remaining sources

trunk/psphot/src/psphotReadoutMinimal.c

@@ -88 +88 @@
 
     // re-measure the kron mags with models subtracted and more appropriate windows
-    psphotKronIterate(config, view, filerule);
+    psphotKronIterate(config, view, filerule, 1);
 
     // measure source size for the remaining sources

trunk/psphot/src/psphotSetThreads.c

@@ -30 +30 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 13);
+    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 14);
     task->function = &psphotKronIterate_Threaded;
     psThreadTaskAdd(task);

trunk/psphot/src/psphotStackImageLoop.c

@@ -21 +21 @@
     psMemDump("startloop");
 
+    pmFPAview *view = pmFPAviewAlloc (0);
     pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
     pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
@@ -37 +38 @@
     }
 
-    pmFPAview *view = pmFPAviewAlloc (0);
-
-
-    // XXX for now, just load the full set of images up front except for EXPNUM which we defer
+    bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
+    bool match_psfs = radial_apertures;
+    bool needConvolved = radial_apertures || !useRaw;
+    if (!needConvolved) {
+        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.INPUT.CNV");
+    }
+
+    // just load the full set of images up front except for EXPNUM which we defer
     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
@@ -63 +68 @@
                 psMemDump("load");
 
-                // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
-                if (!psphotStackMatchPSFs (config, view)) {
-                    psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                    psFree (view);
-                    return false;
-                }
+                if (match_psfs) {
+                    // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
+                    if (!psphotStackMatchPSFs (config, view)) {
+                        psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                        psFree (view);
+                        return false;
+                    }
+                } else {
+                    if (!psphotStackAllocateOutput (config, view, recipe)) {
+                        psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                        psFree (view);
+                        return false;
+                    }
+                }
                 psMemDump("stackmatch");
 

trunk/psphot/src/psphotStackMatchPSFsNext.c

@@ -81 +81 @@
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     psAssert (maskVal, "missing mask value?");
+    psImageMaskType maskSat = psMetadataLookupImageMask(&status, recipe, "MASK.SAT"); // Mask value for bad pixels
+    psAssert (maskSat, "missing mask value?");
+
 
     float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
@@ -123 +126 @@
     psImageSmoothMask_Threaded(readout->variance, readout->variance, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss);
     psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth"));
+
+    // Insure that invalid pixels are masked
+    // XXX: the smoothing seems to generate nan pixels in the variance image
+    // XXX: We may need to loop over the cached sources and redefine the maskObj images...
+    pmReadoutMaskInvalid(readout, maskVal, maskSat);
+    {
+        // Now go rebuild the sources' copies of the mask
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->allSources;
+
+        if (sources && sources->n) {
+            for (int i = 0 ; i < sources->n; i++) {
+                // XXX: move this to a function in pmSource.c
+                pmSource *source = sources->data[i];
+                if (source->maskObj && source->maskView) {
+                    psFree(source->maskObj);
+                    source->maskObj = psImageCopy (source->maskObj, source->maskView, PS_TYPE_IMAGE_MASK);
+                }
+            }
+        }
+    }
 
     psLogMsg("psphot", PS_LOG_INFO, "smoothed");

trunk/psphot/src/psphotStackReadout.c

@@ -2 +2 @@
 
 static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
+static void logMemStats(const char *heading);
 
 // we have 3 possible real filesets:
@@ -51 +52 @@
     // by the multiple threads, not the total time used by all threads.
     psTimerStart ("psphotReadout");
+
+    logMemStats("Start");
 
     pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); // allow models to have ugly fits (eg, central cusp)
@@ -82 +85 @@
     }
 
-    // Generate the mask and weight images (if not supplied) and set mask bits
+    // Generate the mask and weight images (if not supplied) and set mask bits. 
+    // This also insures that all invalid pixels are masked (this is done for STACK_CNV in psphotStackMatchPSFs)
     if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!psphotSetMaskAndVariance (config, view, STACK_OUT)) {
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
@@ -151 +158 @@
     // psphotDumpTest (config, view, STACK_SRC);
     psMemDump("sourcestats");
+    logMemStats("sourcestats");
 
     // classify sources based on moments, brightness
@@ -196 +204 @@
     // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     // but iterates to an appropriately larger size
-    psphotKronIterate(config, view, STACK_SRC);
+    logMemStats("before.kron.1");
+    psphotKronIterate(config, view, STACK_SRC, 1);
+    logMemStats("after.kron.1");
 
     // identify CRs and extended sources
@@ -208 +218 @@
     psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
 
+    logMemStats("pass1");
 
     // if we only do one pass, skip to extended source analysis
@@ -299 +310 @@
     }
 
+    logMemStats("prematch");
+
     // generate the objects (objects unify the sources from the different images) NOTE: could
     // this just match the detections for the chisq image, and not bother measuring the source
@@ -331 +344 @@
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, STACK_SRC);
-    psphotKronIterate(config, view, STACK_SRC);
+    logMemStats("before.kron.2");
+    psphotKronIterate(config, view, STACK_SRC, 2);
+    logMemStats("after.kron.2");
 
     // measure source size for the remaining sources
@@ -353 +368 @@
         return psphotReadoutCleanup (config, view, STACK_SRC);
     }
-
 
     bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
@@ -419 +433 @@
         // psphotEfficiency wants to have the PSF of the image, but since we are measuring on
         // the convolved images we need to generate PSFs for the DET images
-        if (!psphotChoosePSF (config, view, STACK_DET, false)) { // pass 1
+        if (!psphotChoosePSF (config, view, STACK_DET, false)) {
             psLogMsg ("psphot", 3, "failure to construct a psf model for raw input");
             return psphotReadoutCleanup (config, view, STACK_DET);
@@ -429 +443 @@
     }
     psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
+
+    logMemStats("final");
+#if (1)
+    psphotSourceMemory(config, view, STACK_SRC);
+    psphotSourceMemory(config, view, STACK_OUT);
+#endif
 
     // replace failed sources?
@@ -474 +494 @@
     return true;
 }
+
+// read the memory usage data from /proc and log them out
+// This will only work on a system that has /proc (not MacOS for instance) but since this function just
+// tries to open and read from a file it is safe
+static void logMemStats(const char *heading) {
+
+    // file containing memory statistics for this process proc. See proc(5)
+    const char* statm_path = "/proc/self/statm";
+
+    FILE *f = fopen(statm_path,"r");
+    if (!f) {
+        psLogMsg ("psphot", PS_LOG_WARN, "failed to open %s", statm_path);
+        return;
+    }
+
+    unsigned long vmSize, resident, share, text, lib, data;
+
+    int nread;
+    nread = fscanf(f,"%ld %ld %ld %ld %ld %ld", &vmSize, &resident, &share, &text, &lib, &data);
+    fclose(f);
+    if (nread != 6) {
+        psLogMsg ("psphot", PS_LOG_WARN, "failed to read 6 items from %s", statm_path);
+        return;
+    }
+  
+    // assuming 4 KB page size here
+#   define PAGES_TO_MB(_v) (_v * 4.096 / 1024.)
+    psLogMsg ("psphot", PS_LOG_INFO, "Memory usage at %20s: Total VmSize: %8.2f MB   Resident %8.2f MB   Data: %8.2f MB\n",
+        heading, PAGES_TO_MB(vmSize), PAGES_TO_MB(resident), PAGES_TO_MB(data));
+}
+
+
 
 /* here is the process: