Changeset 21166 for trunk/psphot/src/psphotSourceStats.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/psphotSourceStats.c (modified) (7 diffs)

trunk/psphot/src/psphotSourceStats.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, pmDetections *detections)
-{
-    bool     status  = false;
+psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
+
+    bool status = false;
     psArray *sources = NULL;
-    float BIG_RADIUS;
 
     psTimerStart ("psphot.stats");
 
-    // bit-masks to test for good/bad pixels
-    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
-
-    // bit-mask to mark pixels not used in analysis
-    psMaskType markVal = psMetadataLookupU8(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
-
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
+    // 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
 
     // determine properties (sky, moments) of initial sources
-    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    if (!status) return NULL;
     float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
-    if (!status) return NULL;
-    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
-    if (!status) return NULL;
-    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
     if (!status) return NULL;
     char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
@@ -38 +31 @@
     }
 
+    // generate the array of sources, define the associated pixel
     sources = psArrayAllocEmpty (peaks->n);
 
+    for (int i = 0; i < peaks->n; i++) {
+
+        pmPeak *peak = peaks->data[i];
+        if (peak->assigned) continue;
+
+        // create a new source
+        pmSource *source = pmSourceAlloc();
+
+        // add the peak
+        source->peak = psMemIncrRefCounter(peak);
+
+        // allocate space for moments
+        source->moments = pmMomentsAlloc();
+
+        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
+        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
+
+        peak->assigned = true;
+        psArrayAdd (sources, 100, source);
+        psFree (source);
+    }
+
+    if (!strcasecmp (breakPt, "PEAKS")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
+        psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
+        psphotVisualShowMoments (sources);
+        return sources;
+    }
+
+    // threaded measurement of the source magnitudes
     int Nfail = 0;
     int Nmoments = 0;
-    for (int i = 0; i < peaks->n; i++) {
-
-        pmPeak *peak = peaks->data[i];
-        if (peak->assigned) continue;
-
-        // create a new source, add peak
-        pmSource *source = pmSourceAlloc();
-        source->peak = psMemIncrRefCounter(peak);
-
-        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
-        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
-        if (!strcasecmp (breakPt, "PEAKS")) {
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
-            continue;
-        }
+
+    // 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_SOURCE_STATS");
+
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, recipe);
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
+            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[2];
+                Nmoments += scalar->data.S32;
+                scalar = job->args->data[3];
+                Nfail += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+
+    psFree (cellGroups);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot.stats"));
+
+    psphotVisualShowMoments (sources);
+
+    return (sources);
+}
+
+bool psphotSourceStats_Threaded (psThreadJob *job) {
+
+    bool status = false;
+    float BIG_RADIUS;
+    psScalar *scalar = NULL;
+
+    psArray *sources                = job->args->data[0];
+    psMetadata *recipe              = job->args->data[1];
+
+    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
+    if (!status) return false;
+    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
+    if (!status) return false;
+    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    if (!status) return false;
+
+    // bit-masks to test for good/bad pixels
+    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+
+    // bit-mask to mark pixels not used in analysis
+    psMaskType markVal = psMetadataLookupU8(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // threaded measurement of the sources moments
+    int Nfail = 0;
+    int Nmoments = 0;
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
 
         // skip faint sources for moments measurement
         if (source->peak->SN < MIN_SN) {
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
             continue;
         }
@@ -72 +169 @@
         status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
         if (!status) {
-          psFree (source);
-          Nfail ++;
-          psErrorClear();
-          continue;
+            psErrorClear(); // XXX re-consider the errors raised here
+            Nfail ++;
+            continue;
         }
 
@@ -82 +178 @@
         status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
         if (!status) {
-          psFree (source);
-          Nfail ++;
-          psErrorClear();
-          continue;
+            Nfail ++;
+            psErrorClear();
+            continue;
         }
 
@@ -91 +186 @@
         status = pmSourceMoments (source, RADIUS);
         if (status) {
-            // add to the source array
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
             Nmoments ++;
             continue;
@@ -105 +196 @@
         status = pmSourceMoments (source, BIG_RADIUS);
         if (status) {
-            // add to the source array
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
             Nmoments ++;
             continue;
         }
 
-        psFree (source);
         Nfail ++;
         psErrorClear();
@@ -119 +205 @@
     }
 
-    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot.stats"));
-
-    psphotVisualShowMoments (sources);
-
-    return (sources);
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[2];
+    scalar->data.S32 = Nmoments;
+
+    scalar = job->args->data[3];
+    scalar->data.S32 = Nfail;
+    
+    return true;
 }
-
-// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky)