Changeset 21166 for trunk/psphot/src/psphotBlendFit.c

Timestamp:

Jan 24, 2009, 10:54:29 AM (17 years ago)

Author:

eugene

Message:

clean up threading model for psphotGuessModel (defined standard way to split sources into regions); extending the threading to psphotMagnitudes, psphotApReset, psphotBlendFit, psphotSourceStats

File:

• trunk/psphot/src/psphotBlendFit.c (modified) (8 diffs)

trunk/psphot/src/psphotBlendFit.c

@@ -2 +2 @@
 
 // XXX I don't like this name
-bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
+bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
     int Nfit = 0;
@@ -12 +12 @@
     psTimerStart ("psphot.fit.nonlinear");
 
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+    // nThreads = 0; // XXX until testing is complete, do not thread this function
+
+    psphotInitLimitsPSF (recipe, readout);
+    psphotInitLimitsEXT (recipe);
+    psphotInitRadiusPSF (recipe, psf->type);
+
+    // starts the timer, sets up the array of fitSets 
+    psphotFitInit (nThreads);
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // 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);
+
+    fprintf (stderr, "starting with %ld sources\n", sources->n);
+
+    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_BLEND_FIT");
+            psArray *newSources = psArrayAllocEmpty(16);
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, recipe);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+            psArrayAdd(job->args, 1, newSources); // return for new sources
+            psFree (newSources);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (job);
+                return NULL;
+            }
+            psFree(job);
+
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            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 = NULL;
+                scalar = job->args->data[5];
+                Nfit += scalar->data.S32;
+                scalar = job->args->data[6];
+                Npsf += scalar->data.S32;
+                scalar = job->args->data[7];
+                Next += scalar->data.S32;
+                scalar = job->args->data[8];
+                Nfail += scalar->data.S32;
+
+                // add these back onto sources
+                psArray *newSources = job->args->data[4];
+                for (int j = 0; j < newSources->n; j++) {
+                    psArrayAdd (sources, 16, newSources->data[j]);
+                }
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+
+    if (psTraceGetLevel("psphot") >= 6) {
+      psphotSaveImage (NULL, readout->image,  "image.v2.fits");
+    }
+
+    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "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);
+
+    psphotVisualShowResidualImage (readout);
+    psphotVisualShowFlags (sources);
+
+    return true;
+}
+
+bool psphotBlendFit_Threaded (psThreadJob *job) {
+
+    bool status = false;
+    int Nfit = 0;
+    int Npsf = 0;
+    int Next = 0;
+    int Nfail = 0;
+    psScalar *scalar = NULL;
+
+    pmReadout *readout  = job->args->data[0];
+    psMetadata *recipe  = job->args->data[1];
+    psArray *sources    = job->args->data[2];
+    pmPSF *psf          = job->args->data[3];
+    psArray *newSources = job->args->data[4];
+
     // bit-masks to test for good/bad pixels
     psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
@@ -23 +143 @@
     maskVal |= markVal;
 
-    // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortBySN);
-
     // S/N limit to perform full non-linear fits
     float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
-
-    psphotInitLimitsPSF (recipe, readout);
-    psphotInitLimitsEXT (recipe);
-    psphotInitRadiusPSF (recipe, psf->type);
-
-    psphotFitInit ();
 
     // option to limit analysis to a specific region
@@ -41 +152 @@
 
     for (int i = 0; i < sources->n; i++) {
-        // if (i%100 == 0) psphotFitSummary ();
-
         pmSource *source = sources->data[i];
 
@@ -57 +166 @@
         if (source->peak->SN < FIT_SN_LIM) continue;
 
-        // XXX this should use peak?
+        // exclude sources outside optional analysis region
         if (source->peak->xf < AnalysisRegion.x0) continue;
         if (source->peak->yf < AnalysisRegion.y0) continue;
@@ -66 +175 @@
         if (source->modelPSF == NULL) continue;
 
-        // skip sources which are insignificant flux?
+        // skip sources which are insignificant flux? 
+        // XXX this is somewhat ad-hoc
         if (source->modelPSF->params->data.F32[1] < 0.1) {
             psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
@@ -81 +191 @@
         Nfit ++;
 
-        // XXX TEST
-        // if ((fabs(source->peak->x - 1202) < 2) && (fabs(source->peak->y - 1065) < 2)) {
-        // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-        // }
-
-        // try fitting PSFs, then try extended sources
-        // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
-        // XXX re-consider conditions under which the source has EXT fit:
-        // I should try EXT if the source size measurement says it is large
+        // try fitting PSFs or extended sources depending on source->mode
+        // these functions subtract the resulting fitted source
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            if (psphotFitBlob (readout, source, sources, psf, maskVal, markVal)) {
+            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_EXTENDED;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
@@ -115 +218 @@
     }
 
-    if (psTraceGetLevel("psphot") >= 6) {
-      psphotSaveImage (NULL, readout->image,  "image.v2.fits");
-    }
-
-    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "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);
-
-    psphotVisualShowResidualImage (readout);
-    psphotVisualShowFlags (sources);
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[5];
+    scalar->data.S32 = Nfit;
+
+    scalar = job->args->data[6];
+    scalar->data.S32 = Npsf;
+
+    scalar = job->args->data[7];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[8];
+    scalar->data.S32 = Nfail;
 
     return true;