Changeset 35640 for branches/eam_branches/ipp-20130509/psphot/src

Timestamp:

Jun 9, 2013, 3:31:46 PM (13 years ago)

Author:

eugene

Message:

thread add/sub noise, thread model background, pass chisq conv options to LMM fitting

Location:

branches/eam_branches/ipp-20130509/psphot/src

Files:

• psphot.h (modified) (3 diffs)
• psphotAddNoise.c (modified) (6 diffs)
• psphotBlendFit.c (modified) (6 diffs)
• psphotChoosePSF.c (modified) (2 diffs)
• psphotExtendedSourceFits.c (modified) (6 diffs)
• psphotModelBackground.c (modified) (9 diffs)
• psphotSetThreads.c (modified) (2 diffs)
• psphotSourceFits.c (modified) (6 diffs)

branches/eam_branches/ipp-20130509/psphot/src/psphot.h

@@ -67 +67 @@
 bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
+bool            psphotModelBackground_Threaded (psThreadJob *job);
 
 bool            psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -119 +120 @@
 bool            psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add);
 bool            psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add);
+bool            psphotAddOrSubNoise_Threaded (psThreadJob *job);
 
 bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -243 +245 @@
 bool            psphotFitInit (int nThreads);
 bool            psphotFitSummary (void);
+bool            psphotFitSummaryExtended (void);
+bool            psphotFitInitExtended (void);
 
 bool            psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule);

branches/eam_branches/ipp-20130509/psphot/src/psphotAddNoise.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+bool psphotAddOrSubNoise_Threaded (psThreadJob *job);
 bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal);
 
@@ -32 +33 @@
 }
 
-static int Nmasked = 0;
-
 // the return state indicates if any sources were actually replaced
 bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) {
 
     bool status = false;
+
+    psTimerStart ("psphot.noise");
 
     // find the currently selected readout
@@ -55 +56 @@
     if (!sources) return false;
 
-    psTimerStart ("psphot.noise");
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -84 +89 @@
 
     psphotVisualShowImage (readout);
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_ADD_NOISE");
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, FACTOR,   PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, SIZE,     PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, add,      PS_TYPE_U8);
+
+# if (1)
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise.");
+                return false;
+            }
+# else
+            if (!psphotAddOrSubNoise_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise.");
+                return false;
+            }
+            psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psFree(cellGroups);
+            psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise.");
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            // we have no returned data from this operation
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            }
+            psFree(job);
+        }
+    }
+
+    if (add) {
+        psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
+    } else {
+        psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
+    }
+
+    psFree (cellGroups);
+
+    psphotVisualShowImage (readout);
+
+    return true;
+}
+
+bool psphotAddOrSubNoise_Threaded (psThreadJob *job) {
+
+    psArray *sources = job->args->data[0];
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[1], PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2], PS_TYPE_IMAGE_MASK_DATA);
+    psF32 FACTOR = PS_SCALAR_VALUE(job->args->data[3], F32);
+    psF32 SIZE = PS_SCALAR_VALUE(job->args->data[4], F32);
+    psBool add  = PS_SCALAR_VALUE(job->args->data[5], U8);
 
     // loop over all source
@@ -104 +184 @@
         psphotMaskSource (source, add, markVal);
     }
-    if (add) {
-        psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
-    } else {
-        psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
-    }
-    fprintf (stderr, "masked %d objects\n", Nmasked);
-
-    psphotVisualShowImage (readout);
-
     return true;
 }
@@ -142 +213 @@
         }
     }
-    Nmasked ++;
-
-    return true;
-}
-
+    return true;
+}
+

branches/eam_branches/ipp-20130509/psphot/src/psphotBlendFit.c

@@ -86 +86 @@
     if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
         fitMaxTol = 1.0;
+    }
+
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        gainFactorMode = 0;
     }
 
@@ -119 +131 @@
     fitOptions->mode          = PM_SOURCE_FIT_PSF;
     fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+
+    fitOptions->gainFactorMode = gainFactorMode;
+    fitOptions->chisqConvergence = chisqConvergence;
 
     psphotInitLimitsPSF (recipe, readout);
@@ -211 +226 @@
 
     psLogMsg ("psphot.psphotBlendFit", PS_LOG_WARN, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
+    psphotFitSummary ();
 
     psphotVisualShowResidualImage (readout, false);
@@ -271 +287 @@
 
         // skip non-astronomical objects (very likely defects)
-        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
-        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (source->mode &  PM_SOURCE_MODE_BLEND) goto skip_blend;
+        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) goto skip_cr;
+        if (source->type == PM_SOURCE_TYPE_DEFECT) goto skip_defect;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) goto skip_sat;
 
         // skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) goto skip_sat;
 
         // skip DBL second sources (ie, added by psphotFitBlob
-        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
+        if (source->mode &  PM_SOURCE_MODE_PAIR) goto skip_blend;
 
         // do not include MOMENTS_FAILURES in the fit
-        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) goto skip_generic;
 
         // limit selection to some SN limit
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) continue;
+            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) goto skip_generic;
         } else {
-            if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
+            if (sqrt(source->peak->detValue) < FIT_SN_LIM) goto skip_generic;
         }
         // exclude sources outside optional analysis region
-        if (source->peak->xf < AnalysisRegion.x0) continue;
-        if (source->peak->yf < AnalysisRegion.y0) continue;
-        if (source->peak->xf > AnalysisRegion.x1) continue;
-        if (source->peak->yf > AnalysisRegion.y1) continue;
+        if (source->peak->xf < AnalysisRegion.x0) goto skip_generic;
+        if (source->peak->yf < AnalysisRegion.y0) goto skip_generic;
+        if (source->peak->xf > AnalysisRegion.x1) goto skip_generic;
+        if (source->peak->yf > AnalysisRegion.y1) goto skip_generic;
 
         // if model is NULL, we don't have a starting guess
-        if (source->modelPSF == NULL) continue;
+        if (source->modelPSF == NULL) goto skip_generic;
 
         // skip sources which are insignificant flux?
@@ -307 +323 @@
                      source->modelPSF->params->data.F32[2],
                      source->modelPSF->params->data.F32[3]);
-            continue;
+            goto skip_generic;
         }
 
@@ -357 +373 @@
         pmSourceCacheModel (source, maskVal);
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        continue;
+
+    skip_blend:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_cr:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_defect:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_sat:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_generic:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
     }
 

branches/eam_branches/ipp-20130509/psphot/src/psphotChoosePSF.c

@@ -158 +158 @@
     }
     float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance
+
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        gainFactorMode = 0;
+    }
 
     // options which modify the behavior of the model fitting
@@ -169 +181 @@
     options->fitOptions->mode          = PM_SOURCE_FIT_PSF;
     options->fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
-
     
+    options->fitOptions->gainFactorMode   = gainFactorMode;
+    options->fitOptions->chisqConvergence = chisqConvergence;
+
     psArray *stars = psArrayAllocEmpty (sources->n);
 

branches/eam_branches/ipp-20130509/psphot/src/psphotExtendedSourceFits.c

@@ -53 +53 @@
     psTimerStart ("psphot.extended");
 
+    psphotFitInitExtended();
+
     // find the currently selected readout
     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
@@ -107 +109 @@
         fitMaxTol = 1.0;
     }
+
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        gainFactorMode = 0;
+    }
+
+    // Define source fitting parameters for extended source fits
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode           = PM_SOURCE_FIT_EXT;
+    fitOptions->saveCovariance = true;  // XXX make this a user option?
+    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    fitOptions->nIter          = fitIter;
+    fitOptions->minTol         = fitMinTol;
+    fitOptions->maxTol         = fitMaxTol;
+
+    fitOptions->gainFactorMode   = gainFactorMode;
+    fitOptions->chisqConvergence = chisqConvergence;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -196 +222 @@
             psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
             psArrayAdd(job->args, 1, iter);
-            psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            psArrayAdd(job->args, 1, fitOptions);
 
             PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitIter, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitMinTol, PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitMaxTol, PS_TYPE_F32);
-
             PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
@@ -291 +314 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit, %d failed\n", Nfaint, Nfail);
+
+    psphotFitSummaryExtended();
+
     return true;
 }
@@ -314 +340 @@
     psMetadataIterator *iter = job->args->data[3];
     psRegion *region        = job->args->data[4];
-    int psfSize             = PS_SCALAR_VALUE(job->args->data[5],S32);
-    int fitIter             = PS_SCALAR_VALUE(job->args->data[6],S32);
-    float fitMinTol         = PS_SCALAR_VALUE(job->args->data[7],F32);
-    float fitMaxTol         = PS_SCALAR_VALUE(job->args->data[8],F32);
-    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[9],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[10],PS_TYPE_IMAGE_MASK_DATA);
-
-    // Define source fitting parameters for extended source fits
-    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode           = PM_SOURCE_FIT_EXT;
-    fitOptions->saveCovariance = true;  // XXX make this a user option?
-    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
-    fitOptions->nIter          = fitIter;
-    fitOptions->minTol         = fitMinTol;
-    fitOptions->maxTol         = fitMaxTol;
+    pmSourceFitOptions *fitOptions = job->args->data[5];
+
+    int psfSize             = PS_SCALAR_VALUE(job->args->data[6],S32);
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[8],PS_TYPE_IMAGE_MASK_DATA);
 
     // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
@@ -540 +556 @@
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     }
-    psFree (fitOptions);
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[9];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[10];
+    scalar->data.S32 = Nconvolve;
+
     scalar = job->args->data[11];
-    scalar->data.S32 = Next;
+    scalar->data.S32 = NconvolvePass;
 
     scalar = job->args->data[12];
-    scalar->data.S32 = Nconvolve;
+    scalar->data.S32 = Nplain;
 
     scalar = job->args->data[13];
-    scalar->data.S32 = NconvolvePass;
+    scalar->data.S32 = NplainPass;
 
     scalar = job->args->data[14];
-    scalar->data.S32 = Nplain;
+    scalar->data.S32 = Nfaint;
 
     scalar = job->args->data[15];
-    scalar->data.S32 = NplainPass;
-
-    scalar = job->args->data[16];
-    scalar->data.S32 = Nfaint;
-
-    scalar = job->args->data[17];
     scalar->data.S32 = Nfail;
 

branches/eam_branches/ipp-20130509/psphot/src/psphotModelBackground.c

@@ -28 +28 @@
 }
 
+// track background cell failures for log reporting purposes
+static int nFailures = 0;
 
 // Generate the background model
@@ -33 +35 @@
 // linear interpolation to generate full-scale model
 //
-// NOTE that the 'analysis' metedata pass in here is used to store the binning information.
+// NOTE that the 'analysis' metedata passed in here is used to store the binning information.
 // This may be the analysis for this readout, but it may be the analysis for the pmFPAfile
 // corresponding to the model.  Other information about the background model is saved on the
@@ -109 +111 @@
 
     const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth;
-    psStats *stats = psStatsAlloc (statsOption);
+    psStats *statsDefaults = psStatsAlloc (statsOption);
 
     // set range for old-version of sky statistic
     if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) {
-        stats->min = 0.25;
-        stats->max = 0.75;
+        statsDefaults->min = 0.25;
+        statsDefaults->max = 0.75;
     }
 
     // set user-option for number of pixels per region
-    stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
+    statsDefaults->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
     if (!status) {
-        stats->nSubsample = 1000;
+        statsDefaults->nSubsample = 1000;
     }
 
     // optionally set the binsize
-    stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
+    statsDefaults->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
     if (status) {
-        stats->options |= PS_STAT_USE_BINSIZE;
+        statsDefaults->options |= PS_STAT_USE_BINSIZE;
     }
 
     // optionally set the sigma clipping
-    stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
+    statsDefaults->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
     if (!status) {
-        if (stats->options & PS_STAT_FITTED_MEAN) {
-            stats->clipSigma = 1.0;
+        if (statsDefaults->options & PS_STAT_FITTED_MEAN) {
+            statsDefaults->clipSigma = 1.0;
         } else {
-            stats->clipSigma = 3.0;
+            statsDefaults->clipSigma = 3.0;
         }
     }
-
-    // stats is not initialized by psStats???  use this to save the input options
-    // XXX re-work this to use the new psStatsInit function
-    psStats *statsDefaults = psStatsAlloc (statsOption);
-    *statsDefaults = *stats;
 
     // we save the binning structure for use in psphotMagnitudes
@@ -154 +151 @@
 
     // XXXX we can thread this here by running blocks in parallel
-
-    int nFailures = 0;
     psImageBackgroundInit();
+    nFailures = 0; // reset for this pass
 
     // we have Nx * Ny model points, but we can use a window which is larger (or smaller) than
@@ -163 +159 @@
 
     // measure clipped median for subimages
-    psRegion ruffRegion = psRegionSet (0,0,0,0);
-    psRegion fineRegion = psRegionSet (0,0,0,0);
     for (int iy = 0; iy < model->numRows; iy++) {
         for (int ix = 0; ix < model->numCols; ix++) {
-
-            // convert the ruff grid cell to the equivalent fine grid cell
-            ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, 
-                                      iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample);
-            fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
-            fineRegion = psRegionForImage (image, fineRegion);
-
-            psImage *subset  = psImageSubset (image, fineRegion);
-            if (!subset->numCols || !subset->numRows) {
-                psFree (subset);
-                continue;
+            
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_MODEL_BACKGROUND");
+
+            psArrayAdd(job->args, 1, image);
+            psArrayAdd(job->args, 1, mask);
+            psArrayAdd(job->args, 1, binning);
+            psArrayAdd(job->args, 1, rng);
+            psArrayAdd(job->args, 1, statsDefaults);
+
+            psArrayAdd(job->args, 1, modelData);
+            psArrayAdd(job->args, 1, modelStdevData);
+
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+
+            PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, dXsample, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, dYsample, PS_TYPE_F32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, statsOptionLocation, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, statsOptionWidth, PS_TYPE_S32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, NAN, PS_TYPE_F32); // this is used as a return value for dQvalue
+
+# if (1)
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                return NULL;
             }
-            psImage *submask = psImageSubset (mask, fineRegion);
-
-            // reset the default values
-            *stats = *statsDefaults;
-
-            // Use the selected background statistic for the first pass
-            // If it fails, fall back on the "ROBUST_MEDIAN" version
-            // If both fail, set the pixel to NAN and (below) interpolate
-            // XXX psImageBackground will probably be renamed psImageStats
-            // XXX don't bother trying if there are no valid pixels...
-
-            psVector *sample = NULL;
-
-            // turn on stats tracing in desired cells
-            # if (0)
-            psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS");
-            assert (plots);
-
-            int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X");
-            assert (status);
-            int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y");
-            assert (status);
-
-            bool gotX = (xPlot < 0) || (xPlot == ix);
-            bool gotY = (yPlot < 0) || (yPlot == iy);
-
-            if (gotX && gotY) {
-                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6);
-                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
-            } else {
-                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0);
-                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
+# else
+            if (!psphotModelBackground_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Failure to model background.");
+                return NULL;
             }
-            # endif
-
-            if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
-                if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-                    modelData[iy][ix] = stats->robustMedian;
-                } else {
-                    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
-                }
-                modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
-                // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
-                psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
-
-                // supply sample to plotting routing
-                psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
-            } else {
-                psStatsOptions currentOptions = stats->options;
-                stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
-                if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
-                    if ((nFailures < 3) || (nFailures % 100 == 0)) {
-                        psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
-                                  "(%dx%d, (row0,col0) = (%d,%d)",
-                                  subset->numRows, subset->numCols, subset->row0, subset->col0);
-                    }
-                    nFailures ++;
-                    modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
-                } else {
-                    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
-                    modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
-
-                    // supply sample to plotting routing
-                    psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
-                }
-                // drop errors caused by psImageBackground failures
-                // XXX we probably should trap and exit on serious failures
-                psErrorClear();
-                stats->options = currentOptions;
-            }
-            psFree(sample);
-            modelData[iy][ix] += SKY_BIAS;
-            psFree (subset);
-            psFree (submask);
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = job->args->data[14];
+                float dQvalue = scalar->data.F32;
+                psVectorAppend (dQ, dQvalue);
+            }
+            psFree(job);
+# endif
         }
+    }
+
+    // wait for the threads to finish and manage results
+    if (!psThreadPoolWait (false, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+        return NULL;
+    }
+
+    // we have only supplied one type of job, so we can assume the types here
+    psThreadJob *job = NULL;
+    while ((job = psThreadJobGetDone()) != NULL) {
+        if (job->args->n < 1) {
+            fprintf (stderr, "error with job\n");
+        } else {
+            psScalar *scalar = job->args->data[14];
+            float dQvalue = scalar->data.F32;
+            psVectorAppend (dQ, dQvalue);
+        }
+        psFree(job);
     }
 
@@ -302 +274 @@
     if (Count == 0) {
         psError (PSPHOT_ERR_DATA, true, "failed to build background image");
-        psFree(stats);
         psFree(statsDefaults);
         psFree(binning);
@@ -318 +289 @@
             modelData[iy][ix] = Value;
             modelStdevData[iy][ix] = ValueStdev;
+        }
+    }
+
+    // apply artificial offset if desired
+    for (int iy = 0; iy < model->numRows; iy++) {
+        for (int ix = 0; ix < model->numCols; ix++) {
+            modelData[iy][ix] += SKY_BIAS;
         }
     }
@@ -364 +342 @@
     psFree(dQ);
 
-    psFree(stats);
     psFree(statsDefaults);
     psFree(binning);
@@ -433 +410 @@
     return true;
 }
+
+bool psphotModelBackground_Threaded (psThreadJob *job) {
+
+    psImage *image          = job->args->data[0];
+    psImage *mask           = job->args->data[1];
+
+    psImageBinning *binning = job->args->data[2];
+
+    psRandom *rng           = job->args->data[3];
+    psStats *statsDefaults  = job->args->data[4];
+
+    psF32 **modelData       = job->args->data[5];
+    psF32 **modelStdevData  = job->args->data[6];
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
+
+    int ix                  = PS_SCALAR_VALUE(job->args->data[8],S32);
+    int iy                  = PS_SCALAR_VALUE(job->args->data[9],S32);
+
+    float dXsample          = PS_SCALAR_VALUE(job->args->data[10],F32);
+    float dYsample          = PS_SCALAR_VALUE(job->args->data[11],F32);
+
+    psStatsOptions statsOptionLocation = PS_SCALAR_VALUE(job->args->data[12],S32);
+    psStatsOptions statsOptionWidth    = PS_SCALAR_VALUE(job->args->data[13],S32);
+
+    // convert the ruff grid cell to the equivalent fine grid cell
+    psRegion ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample);
+    psRegion fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
+    fineRegion = psRegionForImage (image, fineRegion);
+
+    psImage *subset  = psImageSubset (image, fineRegion);
+    if (!subset->numCols || !subset->numRows) {
+        psFree (subset);
+        return false; // XXX do we / should we fail on this?
+    }
+    psImage *submask = psImageSubset (mask, fineRegion);
+
+    psStats *stats = psStatsAlloc (PS_STAT_NONE);
+    *stats = *statsDefaults;
+
+    // Use the selected background statistic for the first pass
+    // If it fails, fall back on the "ROBUST_MEDIAN" version
+    // If both fail, set the pixel to NAN and (later) interpolate
+
+    psVector *sample = NULL;
+    float dQvalue = NAN;
+
+    if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
+        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+            modelData[iy][ix] = stats->robustMedian;
+        } else {
+            modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
+        }
+        modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
+
+        // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
+        // XXX this operation is not thread safe -- move out somehow...
+        // psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
+        dQvalue = stats->robustUQ + stats->robustLQ - 2*stats->robustMedian;
+        // return dQvalue to main thread
+
+        // supply sample to plotting routing
+        // only allow in a non-threaded context -- can plot more than one cell
+        // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
+    } else {
+        // psStatsOptions currentOptions = stats->options;
+        stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
+        if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
+                          "(%dx%d, (row0,col0) = (%d,%d)",
+                          subset->numRows, subset->numCols, subset->row0, subset->col0);
+            }
+            nFailures ++; // static
+            modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
+        } else {
+            modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
+            modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
+
+            // supply sample to plotting routing
+            // only allow in a non-threaded context -- can plot more than one cell
+            // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
+        }
+        // drop errors caused by psImageBackground failures
+        // NOTE : psStats raises errors when it cannot find a solution; this 
+        // probably should not be errors but exit stats info only, but c'est la code
+        // XXX we probably should trap and exit on serious failures
+        psErrorClear();
+        // stats->options = currentOptions; // this is not needed (set by init above)
+    }
+    psFree (stats);
+    psFree (sample);
+    psFree (subset);
+    psFree (submask);
+
+    // return the dQvalue to the calling thread
+    psScalar *scalar = job->args->data[14];
+    scalar->data.F32 = dQvalue;
+
+    return true;
+}
+
+// code for in-line plotting from above
+// turn on stats tracing in desired cells
+// XXX this code may need to be re-worked for a threaded context
+# if (0)
+    psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS");
+    assert (plots);
+
+    int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X");
+    assert (status);
+    int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y");
+    assert (status);
+
+    bool gotX = (xPlot < 0) || (xPlot == ix);
+    bool gotY = (yPlot < 0) || (yPlot == iy);
+
+    if (gotX && gotY) {
+        (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6);
+        (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
+    } else {
+        (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0);
+        (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
+    }
+# endif

branches/eam_branches/ipp-20130509/psphot/src/psphotSetThreads.c

@@ -5 +5 @@
     psThreadTask *task = NULL;
 
+    pmPSFThreads ();
+
+    task = psThreadTaskAlloc("PSPHOT_MODEL_BACKGROUND", 15);
+    task->function = &psphotModelBackground_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
     task = psThreadTaskAlloc("PSPHOT_GUESS_MODEL", 5);
     task->function = &psphotGuessModel_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_ADD_NOISE", 6);
+    task->function = &psphotAddOrSubNoise_Threaded;
     psThreadTaskAdd(task);
     psFree(task);
@@ -40 +52 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 18);
+    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 16);
     task->function = &psphotExtendedSourceFits_Threaded;
     psThreadTaskAdd(task);

branches/eam_branches/ipp-20130509/psphot/src/psphotSourceFits.c

@@ -3 +3 @@
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
 
+static int NfitBlend = 0;
+
 static int NfitPSF = 0;
-static int NfitBlend = 0;
+static int NfitIterPSF = 0;
+
 static int NfitDBL = 0;
+static int NfitIterDBL = 0;
+static int NfitPixDBL = 0;
+
 static int NfitEXT = 0;
+static int NfitIterEXT = 0;
+static int NfitPixEXT = 0;
+
 static int NfitPCM = 0;
+static int NfitIterPCM = 0;
+static int NfitPixPCM = 0;
 
 bool psphotFitInit (int nThreads) {
@@ -17 +28 @@
 bool psphotFitSummary () {
 
-    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
-              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf (%d iter), %5d blend: %5d ext (%d iter, %d pix), %5d dbl (%d iter, %d pix) : %6.2f sec\n",
+              NfitPSF, NfitIterPSF, NfitBlend, NfitEXT, NfitIterEXT, NfitPixEXT, NfitDBL, NfitIterDBL, NfitPixDBL, psTimerMark ("psphot.fits"));
+    return true;
+}
+
+bool psphotFitSummaryExtended () {
+
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d ext (%d iter, %d pix), %5d pcm (%d iter, %d pix)\n",
+              NfitEXT, NfitIterEXT, NfitPixEXT, NfitPCM, NfitIterPCM, NfitPixPCM);
+    return true;
+}
+
+bool psphotFitInitExtended () {
+    NfitEXT = 0;
+    NfitIterEXT = 0;
+    NfitPixEXT = 0;
+    NfitPCM = 0;
+    NfitIterPCM = 0;
+    NfitPixPCM = 0;
     return true;
 }
@@ -169 +197 @@
     options.mode = PM_SOURCE_FIT_PSF;
     pmSourceFitModel (source, PSF, &options, maskVal);
+    NfitIterPSF += PSF->nIter;
 
     if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
@@ -440 +469 @@
     options.mode = PM_SOURCE_FIT_PSF;
     pmSourceFitSet (source, modelSet, &options, maskVal);
+    NfitIterDBL += DBL->nIter;
+    NfitPixDBL += DBL->nDOF;
+
     return (modelSet);
 }
@@ -507 +539 @@
 
     pmSourceFitModel (source, model, &options, maskVal);
+    NfitIterEXT += model->nIter;
+    NfitPixEXT += model->nDOF;
 
 # if (PS_TRACE_ON) 
@@ -587 +621 @@
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);  // NOTE : 1687 allocs in here
+    NfitIterPCM += pcm->modelConv->nIter;
+    NfitPixPCM += pcm->modelConv->nDOF;
     if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }