Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 34316) +++ trunk/psphot/src/Makefile.am (revision 34317) @@ -117,4 +117,5 @@ psphotStackObjects.c \ psphotStackPSF.c \ + psphotStackAllocateOutput.c \ psphotCleanup.c @@ -219,5 +220,6 @@ psphotPetrosianVisual.c \ psphotEfficiency.c \ - psphotSetNFrames.c + psphotSetNFrames.c \ + psphotSourceMemory.c # not currently used Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 34316) +++ trunk/psphot/src/psphot.h (revision 34317) @@ -360,4 +360,5 @@ pmReadout **chiReadout, int index); +bool psphotStackAllocateOutput( const pmConfig *config, pmFPAview *view, psMetadata *recipe); bool psphotStackRemoveChisqFromInputs (pmConfig *config, const char *filerule); @@ -483,6 +484,6 @@ 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,3 +510,6 @@ ); +bool psphotSourceMemory(pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index); + #endif Index: trunk/psphot/src/psphotKronIterate.c =================================================================== --- trunk/psphot/src/psphotKronIterate.c (revision 34316) +++ trunk/psphot/src/psphotKronIterate.c (revision 34317) @@ -5,5 +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,5 +42,5 @@ // 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,7 +54,19 @@ 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,5 +76,4 @@ psTimerStart ("psphot.kron"); - // determine the number of allowed threads @@ -236,4 +247,5 @@ 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,4 +310,8 @@ 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,4 +325,9 @@ 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,20 +343,20 @@ } - // 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,4 +366,11 @@ } 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,4 +399,22 @@ } +#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,7 +536,9 @@ 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,5 +594,5 @@ } - source->moments->Mrf = Mrf; + source->moments->Mrf = Mrf; source->moments->KronFlux = Sum; source->moments->KronFluxErr = sqrt(Var); @@ -564,7 +612,7 @@ 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; Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 34316) +++ trunk/psphot/src/psphotReadout.c (revision 34317) @@ -194,5 +194,5 @@ // 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,5 +312,5 @@ // 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 Index: trunk/psphot/src/psphotReadoutMinimal.c =================================================================== --- trunk/psphot/src/psphotReadoutMinimal.c (revision 34316) +++ trunk/psphot/src/psphotReadoutMinimal.c (revision 34317) @@ -88,5 +88,5 @@ // 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 Index: trunk/psphot/src/psphotSetThreads.c =================================================================== --- trunk/psphot/src/psphotSetThreads.c (revision 34316) +++ trunk/psphot/src/psphotSetThreads.c (revision 34317) @@ -30,5 +30,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 13); + task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 14); task->function = &psphotKronIterate_Threaded; psThreadTaskAdd(task); Index: trunk/psphot/src/psphotSourceMemory.c =================================================================== --- trunk/psphot/src/psphotSourceMemory.c (revision 34317) +++ trunk/psphot/src/psphotSourceMemory.c (revision 34317) @@ -0,0 +1,54 @@ +# include "psphotInternal.h" + + +bool psphotSourceMemory (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotSourceMemoryReadout (config, view, filerule, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on memory usage measurement for %s entry %d", filerule, i); + return false; + } + } + return true; +} +bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index) +{ + bool status; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (!detections) { + return true; + } + + psArray *sources = detections->allSources; + if (!sources) { + return true; + } + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping memory measurement"); + return true; + } + + psU64 bytes = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + bytes += pmSourceMemoryUse(source); + } + + psLogMsg ("psphot", PS_LOG_INFO, "input %s %d: %ld sources %.1f MB. %.1f bytes per source.\n", + filerule, index, sources->n, (psF32)bytes/1024/1024, (psF32) bytes / sources->n); + + return true; +} Index: trunk/psphot/src/psphotStackAllocateOutput.c =================================================================== --- trunk/psphot/src/psphotStackAllocateOutput.c (revision 34317) +++ trunk/psphot/src/psphotStackAllocateOutput.c (revision 34317) @@ -0,0 +1,31 @@ +# include "psphotInternal.h" + +// allocate readout and images for PSPHOT.STACK.INPUT.CNV +// This is called only if we are skiping the psf matching +bool psphotStackAllocateOutput (const pmConfig *config, pmFPAview *view, psMetadata *recipe) { + bool status = false; + + bool useRaw = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.USE.RAW"); + char *fileruleSrc = useRaw ? "PSPHOT.STACK.INPUT.RAW" : "PSPHOT.STACK.INPUT.CNV"; + + int num = psphotFileruleCount(config, fileruleSrc); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + // find hte source readout + pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, fileruleSrc, i); + psAssert (fileSrc, "missing file?"); + pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa); + psAssert (readoutSrc, "missing readout?"); + + pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); + psAssert (fileOut, "missing output file?"); + pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa); + if (readoutOut == NULL) { + readoutOut = pmFPAGenerateReadout(config, view, "PSPHOT.STACK.OUTPUT.IMAGE", fileSrc->fpa, NULL, i); + psAssert (readoutOut, "missing readout?"); + } + } + + return true; +} Index: trunk/psphot/src/psphotStackImageLoop.c =================================================================== --- trunk/psphot/src/psphotStackImageLoop.c (revision 34316) +++ trunk/psphot/src/psphotStackImageLoop.c (revision 34317) @@ -21,4 +21,5 @@ 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,8 +38,12 @@ } - 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,10 +68,18 @@ 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"); Index: trunk/psphot/src/psphotStackMatchPSFsNext.c =================================================================== --- trunk/psphot/src/psphotStackMatchPSFsNext.c (revision 34316) +++ trunk/psphot/src/psphotStackMatchPSFsNext.c (revision 34317) @@ -81,4 +81,7 @@ 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,4 +126,27 @@ 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"); Index: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 34316) +++ trunk/psphot/src/psphotStackReadout.c (revision 34317) @@ -2,4 +2,5 @@ static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule); +static void logMemStats(const char *heading); // we have 3 possible real filesets: @@ -51,4 +52,6 @@ // 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,6 +85,10 @@ } - // 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,4 +158,5 @@ // psphotDumpTest (config, view, STACK_SRC); psMemDump("sourcestats"); + logMemStats("sourcestats"); // classify sources based on moments, brightness @@ -196,5 +204,7 @@ // 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,4 +218,5 @@ psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources) + logMemStats("pass1"); // if we only do one pass, skip to extended source analysis @@ -299,4 +310,6 @@ } + 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,5 +344,7 @@ // 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,5 +368,4 @@ return psphotReadoutCleanup (config, view, STACK_SRC); } - bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES"); @@ -419,5 +433,5 @@ // 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,4 +443,10 @@ } 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,4 +494,36 @@ 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: