IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotGuessModels.c

    r20453 r21166  
    11# include "psphotInternal.h"
    22
    3 bool psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads);
    4 psphotGuessModelForRegionArgs *psphotGuessModelForRegionArgsAlloc();
    5 bool psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args);
     3// XXX : the threading here is not great.  this may be due to blocks between elements, but
     4// the selection of the objects in a cell is not optimal.  To fix:
     5// 1) define the boundaries of the cells up front
     6// 2) loop over the sources once and associate them with their cell
     7// 3) define the threaded function to work with sources for a given cell
    68
    79// A guess for when the moments aren't available
    … …  
    2022}
    2123
    22 # define NFILL 4
    23 
    2424// construct an initial PSF model for each object
    2525bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
    … …  
    3333    assert (recipe);
    3434
    35     // int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    36     // if (!status) {
    37     // nThreads = 0;
    38     // }
    39     int nThreads = 0;
     35    // determine the number of allowed threads
     36    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     37    if (!status) {
     38        nThreads = 0;
     39    }
     40    // nThreads = 0; // XXX until testing is complete, do not thread this function
    4041
    4142    // bit-masks to test for good/bad pixels
    … …  
    5354    psphotInitRadiusPSF (recipe, psf->type);
    5455
    55     // the strategy here is to divide the image into 2x2 blocks of cells and cycle through
    56     // the four discontiguous sets of cells, threading all within a set and blocking between
    57     // sets
    58 
    59     // we divide the image region into 2*2 blocks of size Nx*Ny, the image will have
    60     // Cx*Cy blocks so that (2Nx)Cx = numCols, (2Ny)Cy = numRows.  We want to choose Cx and
    61     // Cy so that (2Nx)Cx * (2Ny)Cy = 4 * NFILL * nThreads -- each of the four sets of cells
    62     // has enough cells to allow NFILL cells for each thread (to better distribute heavy and
    63     // light load cells
    64    
    65     // choose Cx, Cy:
     56    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    6657    int Cx = 1, Cy = 1;
    6758    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
    6859
    69     // the array runs from readout->image->col0 to readout->image->col0 + readout->image->numCols
    70     int Xo = readout->image->col0;
    71     int Yo = readout->image->row0;
    72     int Nx = readout->image->numCols / (2*Cx);
    73     int Ny = readout->image->numRows / (2*Cy);
    74 
    75     // we can thread all of the cells within a set, but need to block between sets
    76     for (int jy = 0; jy < 2; jy++) {
    77         for (int jx = 0; jx < 2; jx++) {
    78 
    79             // generate jobs for all of the Cx*Cy cells in this set
    80             for (int ix = 0; ix < Cx; ix++) {
    81                 for (int iy = 0; iy < Cy; iy++) {
    82                
    83                     int x0 = (2*ix + jx)*Nx + Xo;
    84                     int x1 = x0 + Nx;
    85 
    86                     int y0 = (2*iy + jy)*Ny + Yo;
    87                     int y1 = y0 + Ny;
    88 
    89                     if (readout->image->numCols + Xo - x1 < Nx) {
    90                         x1 = readout->image->numCols + Xo;
    91                     }
    92                     if (readout->image->numRows + Yo - y1 < Ny) {
    93                         y1 = readout->image->numRows + Yo;
    94                     }
    95 
    96                     // fprintf (stderr, "launch %d,%d - %d,%d\n", x0, y0, x1, y1);
    97                    
    98                     psRegion *cellRegion = psRegionAlloc (x0, x1, y0, y1);
    99                     *cellRegion = psRegionForImage (readout->image, *cellRegion);
    100                     // if we got the math above right, this should be a NOP
    101                                                                              
    102                     psphotGuessModelForRegionArgs *args = psphotGuessModelForRegionArgsAlloc();
    103                     args->readout = psMemIncrRefCounter(readout);
    104                     args->sources = psMemIncrRefCounter(sources);
    105                     args->psf     = psMemIncrRefCounter(psf);
    106                     args->region  = psMemIncrRefCounter(cellRegion);
    107 
    108                     args->maskVal = maskVal;
    109                     args->markVal = markVal;
    110 
    111                     // allocate a job -- if threads are not defined, this just runs the job
    112                     psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
    113                     psArrayAdd(job->args, 1, args);
    114                     if (!psThreadJobAddPending(job)) {
    115                         psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    116                         return false;
    117                     }
    118                     psFree(cellRegion);
    119                     psFree(job);
    120                     psFree(args);
    121                 }
    122             }
    123 
    124             // wait for the threads to finish and manage results
    125             // wait here for the threaded jobs to finish
    126             // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
    127             if (!psThreadPoolWait (false)) {
     60    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     61
     62    for (int i = 0; i < cellGroups->n; i++) {
     63
     64        psArray *cells = cellGroups->data[i];
     65
     66        for (int j = 0; j < cells->n; j++) {
     67
     68            // allocate a job -- if threads are not defined, this just runs the job
     69            psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
     70            psArrayAdd(job->args, 1, readout);
     71            psArrayAdd(job->args, 1, cells->data[j]); // sources
     72            psArrayAdd(job->args, 1, psf);
     73
     74            // XXX change these to use abstract mask type info
     75            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8);
     76            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_U8);
     77
     78            if (!psThreadJobAddPending(job)) {
    12879                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     80                psFree (job);
    12981                return false;
    13082            }
    131 
    132             // we have only supplied one type of job, so we can assume the types here
    133             psThreadJob *job = NULL;
    134             while ((job = psThreadJobGetDone()) != NULL) {
    135                 // we have no returned data from this operation
    136                 if (job->args->n < 1) {
    137                     fprintf (stderr, "error with job\n");
    138                 }
    139                 psFree(job);
    140             }
    141         }
    142     }
    143    
     83            psFree(job);
     84        }
     85
     86        // wait for the threads to finish and manage results
     87        // wait here for the threaded jobs to finish
     88        // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
     89        if (!psThreadPoolWait (false)) {
     90            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     91            return false;
     92        }
     93
     94        // we have only supplied one type of job, so we can assume the types here
     95        psThreadJob *job = NULL;
     96        while ((job = psThreadJobGetDone()) != NULL) {
     97            // we have no returned data from this operation
     98            if (job->args->n < 1) {
     99                fprintf (stderr, "error with job\n");
     100            }
     101            psFree(job);
     102        }
     103    }
     104
     105    // XXX check that all sources that should have models
     106    int nMiss = 0;
     107    for (int i = 0; i < sources->n; i++) {
     108
     109        pmSource *source = sources->data[i];
     110        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     111            source->mode &= ~PM_SOURCE_MODE_EXT_LIMIT;
     112            continue;
     113        }
     114
     115        nMiss ++;
     116    }
     117    psLogMsg ("psphot.models", 4, "failed to build models for %d objects\n", nMiss);
     118
     119    psFree (cellGroups);
     120
    144121    psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models"));
    145122    return true;
    146123}
    147124
    148 static void psphotGuessModelForRegionArgsFree(psphotGuessModelForRegionArgs *args)
    149 {
    150     psFree(args->readout);
    151     psFree(args->sources);
    152     psFree(args->psf);
    153     psFree(args->region);
    154     return;
    155 }
    156 
    157 psphotGuessModelForRegionArgs *psphotGuessModelForRegionArgsAlloc()
    158 {
    159     psphotGuessModelForRegionArgs *args = psAlloc(sizeof(psphotGuessModelForRegionArgs));
    160     psMemSetDeallocator(args, (psFreeFunc)psphotGuessModelForRegionArgsFree);
    161 
    162     args->readout = NULL;
    163     args->sources = NULL;
    164     args->psf = NULL;
    165     args->region = NULL;
    166 
    167     args->maskVal = 0;
    168     args->markVal = 0;
    169     return args;
    170 }
    171 
    172125// construct models only for sources in the specified region
    173 bool psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args) {
    174 
    175     pmReadout *readout = args->readout;
    176     psArray *sources = args->sources;
    177     pmPSF *psf = args->psf;
    178     psRegion *region = args->region;
    179     psMaskType maskVal = args->maskVal;
    180     psMaskType markVal = args->markVal;
     126bool psphotGuessModel_Threaded (psThreadJob *job) {
     127
     128    pmReadout *readout = job->args->data[0];
     129    psArray *sources   = job->args->data[1];
     130    pmPSF *psf         = job->args->data[2];
     131   
     132    psMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],U8);
     133    psMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],U8);
    181134
    182135    int nSrc = 0;
    … …  
    184137    for (int i = 0; i < sources->n; i++) {
    185138        pmSource *source = sources->data[i];
     139
     140        // XXXX this is just for a test: use this to mark sources for which the model is measured
     141        // check later that all are used.
     142        source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    186143
    187144        // skip non-astronomical objects (very likely defects)
    … …  
    189146        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    190147        if (!source->peak) continue;
    191 
    192         if (source->peak->xf <  region->x0) continue;
    193         if (source->peak->yf <  region->y0) continue;
    194         if (source->peak->xf >= region->x1) continue;
    195         if (source->peak->yf >= region->y1) continue;
    196148
    197149        nSrc ++;
    … …  
    256208    }
    257209
    258     // fprintf (stderr, "%d for region %lf,%lf - %lf,%lf\n", nSrc, region->x0, region->y0, region->x1, region->y1);
    259210    return true;
    260211}
    261212
    262 bool psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads) {
    263 
    264     int nCells = nThreads * NFILL; // number of cells in a single set
    265     int C = sqrt(nCells) + 0.5;
    266    
    267     // we need to assign Cx and Cy based on the dimensionality of the image
    268     // crude way to find most evenly balanced factors of nCells:
    269     for (int i = C; i >= 1; i--) {
    270         int C1 = nCells / C;
    271         int C2 = nCells / C1;
    272         if (C1*C2 != nCells) continue;
    273 
    274         if (readout->image->numRows > readout->image->numCols) {
    275             *Cx = PS_MAX (C1, C2);
    276             *Cy = PS_MIN (C1, C2);
    277         } else {
    278             *Cx = PS_MAX (C1, C2);
    279             *Cy = PS_MIN (C1, C2);
    280         }
    281         return true;
    282     }
    283     *Cx = 1;
    284     *Cy = 1;
    285 
    286     return true;
    287 }
Note: See TracChangeset for help on using the changeset viewer.