Changeset 20399

Timestamp:

Oct 26, 2008, 2:21:25 PM (18 years ago)

Author:

eugene

Message:

adding threading for psphotGuessModels (threads still clashing...)

Location:

branches/eam_branch_20081026/psphot

Files:

• doc/notes.txt (modified) (1 diff)
• src/Makefile.am (modified) (1 diff)
• src/psphot.c (modified) (1 diff)
• src/psphot.h (modified) (3 diffs)
• src/psphotArguments.c (modified) (1 diff)
• src/psphotCleanup.c (modified) (1 diff)
• src/psphotGuessModels.c (modified) (8 diffs)
• src/psphotReadout.c (modified) (2 diffs)
• src/psphotSetThreads.c (added)

branches/eam_branch_20081026/psphot/doc/notes.txt

@@ +1 @@
+
+2008.10.26
+
+  (note: this is unoptimized)
+  0 threads: built models for 3921 objects: 15.223455 sec
+             built models for 2868 objects: 10.698732 sec
+
+  2 threads: built models for 3921 objects: 37.529655 sec
+             built models for 2862 objects: 28.142142 sec
+
+  (clearly, the threading is doing bad things, at least without optimization..)
+
+  old code (no threading):  
+             built models for 3923 objects: 15.265046 sec
+             built models for 2842 objects: 10.840331 sec
+
+  (so, the threaded version without threading is OK)
 
 2008.06.25

branches/eam_branch_20081026/psphot/src/Makefile.am

@@ -18 +18 @@
 
 libpsphot_la_SOURCES = \
+        psphotSetThreads.c             \
         psphotErrorCodes.c             \
         psphotImageQuality.c           \

branches/eam_branch_20081026/psphot/src/psphot.c

@@ -12 +12 @@
     psphotErrorRegister();              // register our error codes/messages
     psphotModelClassInit ();            // load implementation-specific models
+    psphotSetThreads ();
+
 
     // load command-line arguments, options, and system config data

branches/eam_branch_20081026/psphot/src/psphot.h

@@ -13 +13 @@
 #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs
 
+typedef struct {
+    pmReadout *readout;
+    psArray *sources;
+    pmPSF *psf;
+    psRegion *region;
+    psMaskType maskVal;
+    psMaskType markVal;
+} psphotGuessModelForRegionArgs;
 
 // top-level psphot functions
@@ -41 +49 @@
 bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
 bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
-bool            psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
+bool            psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
+bool            psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args);
 bool            psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
 bool            psphotReplaceUnfitSources (psArray *sources, psMaskType maskVal);
@@ -51 +60 @@
 bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
+bool            psphotSetThreads ();
 
 // used by psphotFindDetections

branches/eam_branch_20081026/psphot/src/psphotArguments.c

@@ -29 +29 @@
     // these options override the PSPHOT recipe values loaded from recipe files
     psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE);
+
+    // Number of threads
+    if ((N = psArgumentGet(argc, argv, "-threads"))) {
+        psArgumentRemove(N, &argc, argv);
+        int nThreads = atoi(argv[N]);
+        psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NTHREADS", 0, "number of psphot threads", nThreads);
+        psArgumentRemove(N, &argc, argv);
+
+        // create the thread pool with number of desired threads, supplying our thread launcher function
+        // XXX need to determine the number of threads from the config data
+        psThreadPoolInit (nThreads);
+    }
+    psphotSetThreads();
 
     // run the test model (requires X,Y coordinate)

branches/eam_branch_20081026/psphot/src/psphotCleanup.c

@@ -5 +5 @@
     psFree (config);
 
+    psThreadPoolFinalize ();
     psTimerStop ();
     psMemCheckCorruption (stderr, true);

branches/eam_branch_20081026/psphot/src/psphotGuessModels.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+bool psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads);
+psphotGuessModelForRegionArgs *psphotGuessModelForRegionArgsAlloc();
+bool psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args);
 
 // A guess for when the moments aren't available
@@ -16 +20 @@
 }
 
+# define NFILL 4
+
 // construct an initial PSF model for each object
-bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
+bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
     bool status;
 
     psTimerStart ("psphot");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
 
     // bit-masks to test for good/bad pixels
@@ -37 +52 @@
     psphotInitRadiusPSF (recipe, psf->type);
 
+    // the strategy here is to divide the image into 2x2 blocks of cells and cycle through
+    // the four discontiguous sets of cells, threading all within a set and blocking between
+    // sets 
+
+    // we divide the image region into 2*2 blocks of size Nx*Ny, the image will have 
+    // Cx*Cy blocks so that (2Nx)Cx = numCols, (2Ny)Cy = numRows.  We want to choose Cx and
+    // Cy so that (2Nx)Cx * (2Ny)Cy = 4 * NFILL * nThreads -- each of the four sets of cells
+    // has enough cells to allow NFILL cells for each thread (to better distribute heavy and
+    // light load cells
+    
+    // choose Cx, Cy:
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    // the array runs from readout->image->col0 to readout->image->col0 + readout->image->numCols 
+    int Xo = readout->image->col0;
+    int Yo = readout->image->row0;
+    int Nx = readout->image->numCols / (2*Cx);
+    int Ny = readout->image->numRows / (2*Cy);
+
+    // we can thread all of the cells within a set, but need to block between sets
+    for (int jy = 0; jy < 2; jy++) {
+        for (int jx = 0; jx < 2; jx++) {
+
+            // generate jobs for all of the Cx*Cy cells in this set
+            for (int ix = 0; ix < Cx; ix++) {
+                for (int iy = 0; iy < Cy; iy++) {
+                
+                    int x0 = (2*ix + jx)*Nx + Xo;
+                    int x1 = x0 + Nx;
+
+                    int y0 = (2*iy + jy)*Ny + Yo;
+                    int y1 = y0 + Ny;
+
+                    if (readout->image->numCols + Xo - x1 < Nx) {
+                        x1 = readout->image->numCols + Xo;
+                    }
+                    if (readout->image->numRows + Yo - y1 < Ny) {
+                        y1 = readout->image->numRows + Yo;
+                    }
+
+                    // fprintf (stderr, "launch %d,%d - %d,%d\n", x0, y0, x1, y1);
+                    
+                    psRegion *cellRegion = psRegionAlloc (x0, x1, y0, y1);
+                    *cellRegion = psRegionForImage (readout->image, *cellRegion);
+                    // if we got the math above right, this should be a NOP
+                                                                             
+                    psphotGuessModelForRegionArgs *args = psphotGuessModelForRegionArgsAlloc();
+                    args->readout = psMemIncrRefCounter(readout);
+                    args->sources = psMemIncrRefCounter(sources);
+                    args->psf     = psMemIncrRefCounter(psf);
+                    args->region  = psMemIncrRefCounter(cellRegion);
+
+                    args->maskVal = maskVal;
+                    args->markVal = markVal;
+
+                    // allocate a job -- if threads are not defined, this just runs the job
+                    psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
+                    psArrayAdd(job->args, 1, args);
+                    if (!psThreadJobAddPending(job)) {
+                        psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                        return false;
+                    }
+                    psFree(cellRegion);
+                    psFree(job);
+                    psFree(args);
+                }
+            }
+
+            // wait for the threads to finish and manage results
+            // wait here for the threaded jobs to finish
+            // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
+            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) {
+                // we have no returned data from this operation
+                if (job->args->n < 1) {
+                    fprintf (stderr, "error with job\n");
+                }
+                psFree(job);
+            }
+        }
+    }
+    
+    psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
+    return true;
+}
+
+static void psphotGuessModelForRegionArgsFree(psphotGuessModelForRegionArgs *args)
+{
+    psFree(args->readout);
+    psFree(args->sources);
+    psFree(args->psf);
+    psFree(args->region);
+    return;
+}
+
+psphotGuessModelForRegionArgs *psphotGuessModelForRegionArgsAlloc()
+{
+    psphotGuessModelForRegionArgs *args = psAlloc(sizeof(psphotGuessModelForRegionArgs));
+    psMemSetDeallocator(args, (psFreeFunc)psphotGuessModelForRegionArgsFree);
+
+    args->readout = NULL;
+    args->sources = NULL;
+    args->psf = NULL;
+    args->region = NULL;
+
+    args->maskVal = 0;
+    args->markVal = 0;
+    return args;
+}
+
+// construct models only for sources in the specified region
+bool psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args) {
+
+    pmReadout *readout = args->readout;
+    psArray *sources = args->sources;
+    pmPSF *psf = args->psf;
+    psRegion *region = args->region;
+    psMaskType maskVal = args->maskVal;
+    psMaskType markVal = args->markVal;
+
+    int nSrc = 0;
+
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
@@ -43 +187 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-
+        if (!source->peak) continue;
+
+        if (source->peak->xf <  region->x0) continue;
+        if (source->peak->yf <  region->y0) continue;
+        if (source->peak->xf >= region->x1) continue;
+        if (source->peak->yf >= region->y1) continue;
+
+        nSrc ++;
+        
         // XXX if a source is faint, it will not have moments measured.
         // it must be modelled as a PSF.  In this case, we need to use
@@ -53 +205 @@
         } else {
             // use the source moments, etc to guess basic model parameters
-            modelEXT = pmSourceModelGuess (source, psf->type);
+            modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC
             if (!modelEXT) {
                 modelEXT = wildGuess(source, psf);
@@ -68 +220 @@
 
         // set PSF parameters for this model (apply 2D shape model)
-        pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
+        pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC
         if (modelPSF == NULL) {
             psError(PSPHOT_ERR_PSF, false,
@@ -92 +244 @@
         // XXX need to define the guess flux?
         // set the fit radius based on the object flux limit and the model
+        // this function affects the mask pixels
         psphotCheckRadiusPSF (readout, source, modelPSF, markVal);
 
@@ -99 +252 @@
 
         pmSourceCacheModel (source, maskVal);
-    }
-    psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
+
+    }
+
+    // fprintf (stderr, "%d for region %lf,%lf - %lf,%lf\n", nSrc, region->x0, region->y0, region->x1, region->y1);
     return true;
 }
 
-// XXX do we always know which model is supposed to be used?
+bool psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads) {
+
+    int nCells = nThreads * NFILL; // number of cells in a single set
+    int C = sqrt(nCells) + 0.5;
+    
+    // we need to assign Cx and Cy based on the dimensionality of the image
+    // crude way to find most evenly balanced factors of nCells:
+    for (int i = C; i >= 1; i--) {
+        int C1 = nCells / C;
+        int C2 = nCells / C1;
+        if (C1*C2 != nCells) continue;
+
+        if (readout->image->numRows > readout->image->numCols) {
+            *Cx = PS_MAX (C1, C2);
+            *Cy = PS_MIN (C1, C2);
+        } else {
+            *Cx = PS_MAX (C1, C2);
+            *Cy = PS_MIN (C1, C2);
+        }
+        return true;
+    }
+    *Cx = 1;
+    *Cy = 1; 
+
+    return true;
+}

branches/eam_branch_20081026/psphot/src/psphotReadout.c

@@ -149 +149 @@
 
     // construct an initial model for each object
-    psphotGuessModels (readout, sources, recipe, psf);
+    psphotGuessModels (config, readout, sources, psf);
 
     // XXX test output of models
@@ -207 +207 @@
 
     // create full input models
-    psphotGuessModels (readout, newSources, recipe, psf);
+    psphotGuessModels (config, readout, newSources, psf);
 
     // replace all sources so fit below applies to all at once