Changeset 28405 for trunk/psphot/src/psphotApResid.c

Timestamp:

Jun 18, 2010, 2:25:38 PM (16 years ago)

Author:

Paul Price

Message:

Discovered a race condition in the threading process. The thread process was: psThreadJobAlloc, populate job with arguments, psThreadJobAddPending, psFree the job. This is not thread-safe, because we're decrementing the job's reference counter while the reference count is also being fiddled within the thread system (psThread.c; by pulling the job off the 'pending' queue and putting it on the 'done' queue). The reference count fiddling within the thread system was protected by a mutex, but the external free was not. Therefore, I'm pulling the free into the thread system so it can be protected (without the user having to worry about locks to do the simple case). This means that the user should no longer touch the job once he hands it over to the thread system (unless it's protected by calls to psThreadLock/psThreadUnlock, or it's on the other side of psThreadPoolWait from all the adds so that the threads aren't doing anything).

File:

• trunk/psphot/src/psphotApResid.c (modified) (14 diffs)

trunk/psphot/src/psphotApResid.c

@@ -22 +22 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (i == chisqNum) continue; // skip chisq image
-        if (!psphotApResidReadout (config, view, filerule, i, recipe)) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotApResidReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for %s entry %d", filerule, i);
-            return false;
-        }
+            return false;
+        }
     }
     return true;
@@ -56 +56 @@
 
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
+        return true;
     }
 
@@ -66 +66 @@
     int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     if (!status) {
-        nThreads = 0;
+        nThreads = 0;
     }
 
@@ -128 +128 @@
     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_APRESID_MAGS");
-
-            psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, psf);
-            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
-            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, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
-            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 false;
-            }
-            psFree(job);
-        }
-
-        // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-            return false;
-        }
-
-        // 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];
-                Nskip += scalar->data.S32;
-                scalar = job->args->data[6];
-                Nfail += scalar->data.S32;
-            }
-            psFree(job);
-        }
+        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_APRESID_MAGS");
+
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
+            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, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
+            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.");
+                return false;
+            }
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            return false;
+        }
+
+        // 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];
+                Nskip += scalar->data.S32;
+                scalar = job->args->data[6];
+                Nfail += scalar->data.S32;
+            }
+            psFree(job);
+        }
     }
 
@@ -184 +182 @@
     Npsf = 0;
 
-# ifdef DEBUG    
+# ifdef DEBUG
     FILE *f = fopen ("apresid.dat", "w");
     psAssert (f, "failed open");
@@ -199 +197 @@
         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
 
-        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
-        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
-        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
-            
+        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
+        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
+        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
+
         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
             continue;
         }
 
-        // XXX make this user-configurable?
-        if (source->errMag > 0.01) continue;
+        // XXX make this user-configurable?
+        if (source->errMag > 0.01) continue;
 
         // aperture residual for this source
@@ -221 +219 @@
 
 # ifdef DEBUG
-        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
-                 source->peak->xf, source->peak->yf, 
-                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], 
-                 source->psfMag, source->apMag, source->errMag,
-                 source->modelPSF->params->data.F32[PM_PAR_I0], 
-                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], 
-                 source->modelPSF->params->data.F32[PM_PAR_7]);
+        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
+                 source->peak->xf, source->peak->yf,
+                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
+                 source->psfMag, source->apMag, source->errMag,
+                 source->modelPSF->params->data.F32[PM_PAR_I0],
+                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
+                 source->modelPSF->params->data.F32[PM_PAR_7]);
 # endif
-        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
-        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
-        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
-        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
-        if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
+        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
+        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
+        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
+        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
+        if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
 
         psVectorAppend (mag, source->psfMag);
@@ -253 +251 @@
     if (Npsf < APTREND_NSTAR_MIN) {
         psWarning("Only %d valid aperture residual sources (need %d), giving up", Npsf, APTREND_NSTAR_MIN);
-        goto escape;
-    }
-
-    // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different 
-    // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a 
+        goto escape;
+    }
+
+    // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different
+    // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a
     // user-specified MAX order, which we should respect, regardless of the mode
 
@@ -270 +268 @@
     pmTrend2DMode mode = PM_TREND_MAP;
     if (mode == PM_TREND_MAP) {
-        MaxOrderForStars ++;
-    } 
+        MaxOrderForStars ++;
+    }
     APTREND_ORDER_MAX = PS_MIN (APTREND_ORDER_MAX, MaxOrderForStars);
 
@@ -283 +281 @@
     int NY = readout->image->numRows;
     for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
-        
-        int Nx, Ny;
-        if (NX > NY) {
-            Nx = i;
-            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
-        } else {
-            Ny = i;
-            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
-        }
-
-        float errorFloor;
+
+        int Nx, Ny;
+        if (NX > NY) {
+            Nx = i;
+            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
+        } else {
+            Ny = i;
+            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
+        }
+
+        float errorFloor;
         pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag);
-        if (!apTrend) {
-            continue;
-        }
-
-        // apply ApTrend results
-        // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
-        // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
-        // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
-        // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
+        if (!apTrend) {
+            continue;
+        }
+
+        // apply ApTrend results
+        // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
+        // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+        // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
+        // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
 
         // store the minimum errorFloor and best ApTrend to keep
         if (errorFloor < errorFloorMin) {
             errorFloorMin = errorFloor;
-            psFree (psf->ApTrend);
-            psf->ApTrend = psMemIncrRefCounter(apTrend);
-        }
-        psFree (apTrend);
+            psFree (psf->ApTrend);
+            psf->ApTrend = psMemIncrRefCounter(apTrend);
+        }
+        psFree (apTrend);
     }
     if (psf->ApTrend == NULL) {
         psWarning("Failed to find a valid aperture residual value");
-        goto escape;
+        goto escape;
     }
 
@@ -383 +381 @@
     if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) {
         psWarning("Failed to fit trend for %d x %d map", Nx, Ny);
-        psFree (apTrend);
-        return NULL;
+        psFree (apTrend);
+        return NULL;
     }
     if (apTrend->mode == PM_TREND_MAP) {
-        // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
-        psImageMapRepair (apTrend->map->map);
-        // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
+        // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
+        psImageMapRepair (apTrend->map->map);
+        // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
     }
 
@@ -400 +398 @@
     if (!isfinite(*apResidSysErr)) {
         psWarning("Failed to find systematic error for %d x %d map", Nx, Ny);
-        psFree (apTrend);
-        return NULL;
+        psFree (apTrend);
+        return NULL;
     }
 
@@ -408 +406 @@
 
     if (psTraceGetLevel("psphot") >= 4) {
-        char filename[64];
-        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
+        char filename[64];
+        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
         FILE *dumpFile = fopen (filename, "w");
         for (int i = 0; i < xPos->n; i++) {
@@ -457 +455 @@
         }
 
-        // clear the mask bit and set the circular mask pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
-
-        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
-
-        // clear the mask bit 
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        // clear the mask bit and set the circular mask pixels
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
+
+        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
+
+        // clear the mask bit
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
-        if (!status) {
-            Nskip ++;
-            psTrace ("psphot", 3, "skip : bad source mag");
-            continue;
-        }
-    
-        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
-            Nfail ++;
-            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
-            continue;
-        }
-        source->mode |= PM_SOURCE_MODE_AP_MAGS;
+        if (!status) {
+            Nskip ++;
+            psTrace ("psphot", 3, "skip : bad source mag");
+            continue;
+        }
+
+        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
+            Nfail ++;
+            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
+            continue;
+        }
+        source->mode |= PM_SOURCE_MODE_AP_MAGS;
     }