Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h (revision 33796) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h (revision 33797) @@ -3,4 +3,5 @@ # include # include +# include /* # define GRID_V1 */ @@ -111,4 +112,6 @@ int PARALLEL_MANUAL; int PARALLEL_SERIAL; + +int NTHREADS; int VERBOSE; @@ -236,5 +239,5 @@ void initMosaicBins PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); void initMosaicGrid PROTO((Image *image, off_t Nimage)); -void initMosaics PROTO((Image *image, off_t Nimage)); +void initMosaics PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage)); void initMrel PROTO((Catalog *catalog, int Ncatalog)); void initialize PROTO((int argc, char **argv)); @@ -246,5 +249,5 @@ SkyList *load_images PROTO((FITS_DB *db, char *regionName, SkyRegion *region)); -Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage)); +Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage)); int main PROTO((int argc, char **argv)); Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c (revision 33796) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c (revision 33797) @@ -50,5 +50,5 @@ ALLOCATE (mosaic, Mosaic, NMOSAIC); - ALLOCATE (MosaicToImage, off_t *, NMOSAIC); + ALLOCATE (MosaicToImage, off_t *, NMOSAIC); ALLOCATE (MosaicN_Image, off_t, NMOSAIC); ALLOCATE (MosaicN_IMAGE, off_t, NMOSAIC); @@ -95,9 +95,9 @@ if (found) continue; - /* a new mosaic, define ranges */ + /* a new mosaic, define ranges -- preserve the original values incase this image is not used */ mosaic[Nmosaic].start = start; mosaic[Nmosaic].stop = stop; - mosaic[Nmosaic].Mcal = 0.0; - mosaic[Nmosaic].dMcal = 0.0; + mosaic[Nmosaic].Mcal = 0.0; // note : mosaic stores only offsets relative to the original image values + mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images mosaic[Nmosaic].dMsys = 0.0; mosaic[Nmosaic].Xm = 0.0; @@ -152,6 +152,8 @@ } -/* find mosaic frames (unique time periods & photcode name matches mosaic) */ -void initMosaics (Image *image, off_t Nimage) { +/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...) + this function will also identify the images NOT in the subset which belong to a selected mosaic + */ +void initMosaics (Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage) { off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; @@ -166,10 +168,10 @@ */ - // generate a list of unique start times (these define the mosaics) - ALLOCATE (startTimes, unsigned int, Nimage); - for (i = 0; i < Nimage; i++) { - startTimes[i] = image[i].tzero; - } - sort_times (startTimes, Nimage); + // generate a list of all subset image start times + ALLOCATE (startTimes, unsigned int, Nsubset); + for (i = 0; i < Nsubset; i++) { + startTimes[i] = subset[i].tzero; + } + sort_times (startTimes, Nsubset); Nmosaic = 0; @@ -178,5 +180,6 @@ startTimesMosaic[0] = startTimes[0]; - for (i = 0; i < Nimage; i++) { + // generate a list of the unique start times (these define the mosaics) + for (i = 0; i < Nsubset; i++) { if (startTimes[i] < startTimesMosaic[Nmosaic]) { fprintf (stderr, "error?\n"); @@ -192,4 +195,35 @@ } Nmosaic ++; + + int Nskip, Nmiss, Nmark; + Nskip = Nmiss = Nmark = 0; + // find any start times which match unselected images + for (i = 0; i < Nimage; i++) { + if (inSubset[i]) { + Nskip ++; + continue; + } + + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (image[i].photcode); + if (!pname) continue; + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + + /* set image time range */ + start = subset[i].tzero; + stop = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1); + + /* find a matching mosaic */ + j = findMosaic(startTimesMosaic, Nmosaic, start); + if (j == -1) { + // no matching mosaic + Nmiss ++; + } else { + // mark this mosaic as bad + Nmark ++; + } + } + fprintf (stderr, "%d images, %d skip, %d miss, %d mark\n", (int) Nimage, (int) Nskip, (int) Nmiss, (int) Nmark); // now I have a list of uniq start times, and they are in order @@ -219,21 +253,21 @@ } - ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs + ALLOCATE (ImageToMosaic, off_t, Nsubset); // mosaic to which image belongs // assign each image to a mosaic - for (i = 0; i < Nimage; i++) { + for (i = 0; i < Nsubset; i++) { ImageToMosaic[i] = -1; /* select valid mosaic images by photcode */ - pname = GetPhotcodeNamebyCode (image[i].photcode); + pname = GetPhotcodeNamebyCode (subset[i].photcode); status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); if (status) continue; /* set image time range */ - // start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); - // stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); - - start = image[i].tzero; - stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + // start = subset[i].tzero - MAX(0.01*subset[i].trate*subset[i].NY, 1); + // stop = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1); + + start = subset[i].tzero; + stop = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1); j = findMosaic(startTimesMosaic, Nmosaic, start); @@ -268,7 +302,7 @@ mosaic[j].dMsys = 0.0; mosaic[j].Xm = 0.0; - mosaic[j].flags = image[i].flags; - mosaic[j].secz = image[i].secz; - mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + mosaic[j].flags = subset[i].flags; + mosaic[j].secz = subset[i].secz; + mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode); } @@ -278,7 +312,7 @@ free (startTimesMosaic); - initMosaicGrid (image, Nimage); - - fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic); + initMosaicGrid (subset, Nsubset); + + fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic); return; } @@ -569,7 +603,92 @@ } +typedef struct { + int Nfew; + int Nbad; + int Ncal; + int Ngrid; + int Nrel; + int Nsys; + off_t Nmax; + int PoorImages; + double *list; + double *dlist; + double *Mlist; + double *dMlist; +} SetMmosInfo; + +enum {THREAD_RUN, THREAD_DONE}; + +typedef struct { + int entry; + int state; + Catalog *catalog; + Image *image; + FlatCorrectionTable *flatcorr; + SetMmosInfo info; +} ThreadInfo; + +int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr); +void *setMmos_worker (void *data); +int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr); + +void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) { + info->Nfew = 0; + info->Nbad = 0; + info->Ncal = 0; + info->Nrel = 0; + info->Ngrid = 0; + info->Nsys = 0; + + info->Nmax = Nmax; + info->PoorImages = PoorImages; + + if (allocLists) { + ALLOCATE (info->list, double, Nmax); + ALLOCATE (info->dlist, double, Nmax); + ALLOCATE (info->Mlist, double, Nmax); + ALLOCATE (info->dMlist, double, Nmax); + } +} + +void SetMmosInfoFree (SetMmosInfo *info) { + free (info->list); + free (info->dlist); + free (info->Mlist); + free (info->dMlist); +} + +void SetMmosInfoAccum (SetMmosInfo *summary, SetMmosInfo *results) { + summary->Nfew += results->Nfew ; + summary->Nsys += results->Nsys ; + summary->Nbad += results->Nbad ; + summary->Ncal += results->Ncal ; + summary->Nrel += results->Nrel ; + summary->Ngrid += results->Ngrid; +} + static int npass_output = 0; -int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { +// mutex to lock setMmos_worker operations +static pthread_mutex_t setMmos_mutex = PTHREAD_MUTEX_INITIALIZER; +static int nextMosaic = 0; + +// we have an array of mosaics (mosaic, Nmosaic). we need to hand out mosaics one at a time to +// the worker threads as they need +off_t getNextMosaicForThread () { + + pthread_mutex_lock (&setMmos_mutex); + if (nextMosaic >= Nmosaic) { + pthread_mutex_unlock (&setMmos_mutex); + return (-1); + } + int thisMosaic = nextMosaic; + nextMosaic ++; + + pthread_mutex_unlock (&setMmos_mutex); + return (thisMosaic); +} + +int setMmos_old (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { off_t i, j, m, c, n, N, Nmax; @@ -798,4 +917,351 @@ } +int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + off_t i, N, Nmax; + Image *image; + + if (!MOSAIC_ZEROPT) return (FALSE); + if (FREEZE_MOSAICS) return (FALSE); + + if (NTHREADS) { + int status = setMmos_threaded (catalog, PoorImages, flatcorr); + return status; + } + + image = getimages (&N, NULL); + + fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); + + if (PoorImages) { + // XXX use bad stars and measurements for PoorImages? or not? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + IMAGE_BAD = 0; + } + + Nmax = 0; + for (i = 0; i < Nmosaic; i++) { + Nmax = MAX (Nmax, N_onMosaic[i]); + } + + SetMmosInfo info; + SetMmosInfoInit (&info, Nmax, TRUE, PoorImages); + + for (i = 0; i < Nmosaic; i++) { + setMmos_mosaic (mosaic, i, image, catalog, &info, flatcorr); + } + SetMmosInfoFree (&info); + + npass_output ++; + + fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys); + + if (PoorImages) { + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; + STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; + } + return (TRUE); +} + +int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) { + + off_t j; + StatType stats; + + double *list = info->list; + double *dlist = info->dlist; + double *Mlist = info->Mlist; + double *dMlist = info->dMlist; + + /* on PoorImages run, skip good images */ + if (info->PoorImages) { + int bad = mosaic[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); + if (!bad) return TRUE; + } + + // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0: + // we optionally do not recalibrate images with UBERCAL zero points + if (mosaic[0].flags & ID_IMAGE_PHOTOM_UBERCAL) { + mosaic[0].ubercalDist = 0; + // propagate ubercalDist to the images + for (j = 0; j < MosaicN_Image[Nmos]; j++) { + off_t im = MosaicToImage[Nmos][j]; + image[im].ubercalDist = mosaic[0].ubercalDist; + // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist); + } + if (KEEP_UBERCAL) return TRUE; + } + + int Nsecfilt = GetPhotcodeNsecfilt (); + + int minUbercalDist = 1000; + + int testImage = FALSE; + // testImage |= (abs(mosaic[0].start - 1323003245) < 10); + // testImage |= (abs(mosaic[0].start - 1323003069) < 10); + // testImage |= (abs(mosaic[0].start - 1323003125) < 10); + // testImage |= (abs(mosaic[0].start - 1323003300) < 10); + // testImage |= (abs(mosaic[0].start - 1323003365) < 10); + // testImage |= (abs(mosaic[0].start - 1323003191) < 10); + // testImage |= (abs(mosaic[0].start - 1323003014) < 10); + // testImage |= (abs(mosaic[0].start - 1323003484) < 10); + // testImage |= (abs(mosaic[0].start - 1323003419) < 10); + // testImage |= (abs(mosaic[0].start - 1323002949) < 10); + + FILE *fout = NULL; + if (testImage) { + char filename[64]; + snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output); + fout = fopen (filename, "w"); + } + + // number of stars to measure the bright-end scatter + int Nbright = 0; + + int N = 0; + for (j = 0; j < N_onMosaic[Nmos]; j++) { + float Msys, Mrel, Mcal, Mgrid, Mflat; + + off_t m = MosaicToMeasure[Nmos][j]; + off_t c = MosaicToCatalog[Nmos][j]; + + if (fout) { + Mcal = getMcal (m, c, flatcorr, catalog); + Mgrid = getMgrid (m, c); + Mrel = getMrel (catalog, m, c); + Mflat = getMflat (m, c, flatcorr, catalog); + + off_t n = catalog[c].measureT[m].averef; + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); + + float delta = Msys - Mrel - Mcal - Mgrid + Mflat; + + int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); + + fprintf (fout, "%f %f : %f %f %f %f %f : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mcal, Mgrid, Mflat, delta, isBad); + } + + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { + info->Nbad ++; + return TRUE; + } + Mcal = getMcal (m, c, flatcorr, catalog); + if (isnan(Mcal)) { + info->Ncal++; + return TRUE; + } + Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) { + info->Ngrid ++; + return TRUE; + } + Mrel = getMrel (catalog, m, c); + if (isnan(Mrel)) { + info->Nrel ++; + return TRUE; + } + + // image.Mcal is not supposed to include the flat-field correction, so we need to + // apply that offset as well here for this image (in other words, each detection is + // being compared to the model, excluding the zero point, Mcal. The model includes + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) + + Mflat = getMflat (m, c, flatcorr, catalog); + + off_t n = catalog[c].measureT[m].averef; + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); + if (isnan(Msys)) { + info->Nsys++; + return TRUE; + } + + PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); + if (!code) goto skip; + if (code->equiv < 1) goto skip; + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec == -1) goto skip; + minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist); + + skip: + list[N] = Msys - Mrel - Mcal - Mgrid + Mflat; + dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { + Mlist[Nbright] = list[N]; + dMlist[Nbright] = dlist[N]; + Nbright ++; + } + N++; + } + /* N_onMosaic[Nmos] is all measurements, N is good measurements */ + + if (fout) { + fclose (fout); + } + + /* too few good measurements or too many bad measurements (skip in PoorImages run) */ + if (!info->PoorImages) { + int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]); + if (mark) { + if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[MosaicToImage[Nmos][0]].name, Nmos, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, N_onMosaic[Nmos]); } + mosaic[0].flags |= ID_IMAGE_PHOTOM_FEW; + info->Nfew ++; + if (testImage) { + fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos])); + } + } else { + mosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW; + } + } + liststats (list, dlist, N, &stats); + if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N); + + mosaic[0].Mcal = stats.mean; + mosaic[0].dMcal = stats.error; + mosaic[0].nFitPhotom = N; + mosaic[0].Xm = 100.0*log10(stats.chisq); + + if (testImage) { + fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, mosaic[0].start, stats.mean, stats.error, mosaic[0].nFitPhotom); + } + + plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); + + // bright end scatter + liststats (Mlist, dMlist, Nbright, &stats); + mosaic[0].dMsys = stats.sigma; + + if (mosaic[0].Mcal < -CLOUD_TOLERANCE) { + mosaic[0].Mcal = 0.0; + } + + if (testImage) { + fprintf (stderr, "%f %f\n", mosaic[0].Mcal, mosaic[0].dMsys); + } + + // minUbercalDist calculated here is the min value for any star owned by this image + // since this particular image is tied to that star, bump its distance by 1 + mosaic[0].ubercalDist = minUbercalDist + 1; + + // propagate ubercalDist to the images + for (j = 0; j < MosaicN_Image[Nmos]; j++) { + off_t im = MosaicToImage[Nmos][j]; + image[im].ubercalDist = mosaic[0].ubercalDist; + // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist); + } + return TRUE; +} + +int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + int i; + off_t N; + + Image *image = getimages (&N, NULL); + + fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); + + if (PoorImages) { + // XXX use bad stars and measurements for PoorImages? or not? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + IMAGE_BAD = 0; + } + + off_t Nmax = 0; + for (i = 0; i < Nmosaic; i++) { + Nmax = MAX (Nmax, N_onMosaic[i]); + } + + SetMmosInfo summary; + SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages); + + pthread_attr_t attr; + pthread_attr_init (&attr); + pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED); + + pthread_t *threads; + ALLOCATE (threads, pthread_t, NTHREADS); + + ThreadInfo *threadinfo; + ALLOCATE (threadinfo, ThreadInfo, NTHREADS); + + // launch N worker threads + for (i = 0; i < NTHREADS; i++) { + threadinfo[i].entry = i; + threadinfo[i].state = THREAD_RUN; + threadinfo[i].catalog = catalog; + threadinfo[i].image = image; + threadinfo[i].flatcorr = flatcorr; + SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages); + pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]); + } + pthread_attr_destroy (&attr); + + // wait until all threads have finished + while (1) { + int allDone = TRUE; + for (i = 0; i < NTHREADS; i++) { + if (threadinfo[i].state == THREAD_RUN) allDone = FALSE; + } + if (allDone) { + break; + } + usleep (500000); + } + + // all threads are done, free the threads array and grab the info + free (threads); + + // report stats & summary from the threads + for (i = 0; i < NTHREADS; i++) { + fprintf (stderr, "setMmos thread %d : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", + i, + threadinfo[i].info.Nfew, + threadinfo[i].info.Nbad, + threadinfo[i].info.Ncal, + threadinfo[i].info.Ngrid, + threadinfo[i].info.Nrel, + threadinfo[i].info.Nsys); + SetMmosInfoAccum (&summary, &threadinfo[i].info); + } + fprintf (stderr, "total : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", + summary.Nfew, + summary.Nbad, + summary.Ncal, + summary.Ngrid, + summary.Nrel, + summary.Nsys); + free (threadinfo); + + return TRUE; +} + +void *setMmos_worker (void *data) { + + ThreadInfo *threadinfo = data; + + SetMmosInfo results; + SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here + + while (1) { + + off_t i = getNextMosaicForThread(); + if (i == -1) { + threadinfo->state = THREAD_DONE; + return NULL; + } + + Catalog *catalog = threadinfo->catalog; + FlatCorrectionTable *flatcorr = threadinfo->flatcorr; + Image *image = threadinfo->image; + + setMmos_mosaic (mosaic, i, image, catalog, &results, flatcorr); + SetMmosInfoAccum (&threadinfo->info, &results); + } + + SetMmosInfoFree (&results); + return NULL; +} + // When we rationalize the images/mosaics, we are driving the negative cloud images back // to 0.0. At the same time, we make a guess to the effective impact on all other images, Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c (revision 33796) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c (revision 33797) @@ -2,13 +2,33 @@ static int Nmax; -static double *list; -static double *dlist; - -// When we rationalize the images/mosaics, we are driving the negative cloud images back -// to 0.0. At the same time, we make a guess to the effective impact on all other images, -// driven by the coupling of common stars. This array carries the impact of those offsets -// on each star -static double *Moffset; - + +typedef struct { + int Nfew; + int Nsys; + int Nbad; + int Ncal; + int Nmos; + int Ngrid; + double *list; + double *dlist; +} SetMrelInfo; + +enum {THREAD_RUN, THREAD_DONE}; + +typedef struct { + int entry; + int state; + Catalog *catalog; + int Ncatalog; + FlatCorrectionTable *flatcorr; + SetMrelInfo summary; +} ThreadInfo; + +void *setMrel_worker (void *data); +int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr); +int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); + +// we want to allocate the stats list,dlist arrays only once (or once per thread). +// this function finds the largest array so we can allocate that max size when needed void initMrel (Catalog *catalog, int Ncatalog) { @@ -21,8 +41,4 @@ } } - - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); - ALLOCATE (Moffset, double, MAX (1, Nmax)); } @@ -47,108 +63,273 @@ } +void SetMrelInfoInit (SetMrelInfo *results, int allocLists) { + results->Nfew = 0; + results->Nsys = 0; + results->Nbad = 0; + results->Ncal = 0; + results->Nmos = 0; + results->Ngrid = 0; + if (allocLists) { + ALLOCATE (results->list, double, Nmax); + ALLOCATE (results->dlist, double, Nmax); + } +} + +void SetMrelInfoFree (SetMrelInfo *results) { + free (results->list); + free (results->dlist); +} + +void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) { + summary->Nfew += results->Nfew ; + summary->Nsys += results->Nsys ; + summary->Nbad += results->Nbad ; + summary->Ncal += results->Ncal ; + summary->Nmos += results->Nmos ; + summary->Ngrid += results->Ngrid; +} + +// mutex to lock setMrel_worker operations +static pthread_mutex_t setMrel_mutex = PTHREAD_MUTEX_INITIALIZER; +static int nextCatalog = 0; + +// we have an array of catalogs (catalog, Ncatalog). we need to hand out catalogs one at a time to +// the worker threads as they need +off_t getNextCatalogForThread (int Ncatalog) { + + pthread_mutex_lock (&setMrel_mutex); + if (nextCatalog >= Ncatalog) { + pthread_mutex_unlock (&setMrel_mutex); + return (-1); + } + int thisCatalog = nextCatalog; + nextCatalog ++; + + pthread_mutex_unlock (&setMrel_mutex); + return (thisCatalog); +} + int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { + int i; + + if (NTHREADS) { + int status = setMrel_threaded (catalog, Ncatalog, flatcorr); + return status; + } + + int Nsecfilt = GetPhotcodeNsecfilt (); + + SetMrelInfo summary, results; + SetMrelInfoInit (&summary, FALSE); + SetMrelInfoInit (&results, TRUE); + + for (i = 0; i < Ncatalog; i++) { + setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt); + SetMrelInfoAccum (&summary, &results); + } + fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + + SetMrelInfoFree (&results); + + return (TRUE); +} + +// set the Mrel values for the specified catalog +// NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array +int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { + off_t j, k, m; - int i, N, Nfew, Nsys, Nbad, Ncal, Nmos, Ngrid; + int N; float Msys, Mcal, Mmos, Mgrid; StatType stats; + double *list = results->list; + double *dlist = results->dlist; + + SetMrelInfoInit (results, FALSE); // do not allocate list,dlist arrays + + for (j = 0; j < catalog[Nc].Naverage; j++) { + // XXX accumulate all secfilt values in a single pass? + + int minUbercalDist = 1000; + + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + /* calculate the average mag in this SEC photcode for a single star */ + + // skip bad stars + if (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + + N = 0; + m = catalog[Nc].averageT[j].measureOffset; + for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) { + + // skip measurements that do not match the current photcode + PhotCode *code = GetPhotcodebyCode (catalog[Nc].measureT[m].photcode); + if (!code) continue; + if (code->equiv != thisCode) { continue; } + + if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) { results->Nbad ++; continue; } + + if (getImageEntry (m, Nc) < 0) { + Mcal = Mmos = Mgrid = 0; + } else { + Mcal = getMcal (m, Nc, flatcorr, catalog); + if (isnan(Mcal)) { results->Ncal ++; continue; } + Mmos = getMmos (m, Nc); + if (isnan(Mmos)) { results->Nmos ++; continue; } + Mgrid = getMgrid (m, Nc); + if (isnan(Mgrid)) { results->Ngrid++; continue; } + } + + Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]); + if (isnan(Msys)) { results->Nsys++; continue; } + list[N] = Msys - Mcal - Mmos - Mgrid; + + int myUbercalDist = getUbercalDist(m,Nc); + minUbercalDist = MIN(minUbercalDist, myUbercalDist); + + // dlist gives the error, which is used as the weight in WT_MEAN. + // we can modify the resulting weight in a few ways: + // 1) MIN_ERROR guarantees a floor + // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode + // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages + // 4) some reference photcode of some kind can be specified as fixed and have a high weight + dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR); + + // up-weight the ubercal values (or convergence can take a long time...) + if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + dlist[N] = MAX (0.1*catalog[Nc].measureT[m].dM, MIN_ERROR); + } + + // tie down reference photometry if the -refcode (code) option is selected + // eg, -refcode g_SDSS + // this probably makes no sense in the context of multifilter analysis + if (refPhotcode) { + if (code->code == refPhotcode->code) { + // tiny error -> large weight + // dlist[N] = MAX (0.01*catalog[Nc].measureT[m].dM, MIN_ERROR); + dlist[N] = 0.0001; + } + } + N++; + } + + // when performing the grid analysis, STAR_TOOFEW will be set to 1; + if (N <= STAR_TOOFEW) { /* too few measurements */ + // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D, catalog[Nc].measureT[catalog[Nc].averageT[j].measureOffset].imageID, N, STAR_TOOFEW); + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW; + results->Nfew ++; + } else { + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW; + } + + liststats (list, dlist, N, &stats); + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].M = stats.mean; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; + } + } + + // values which need to be passed out: Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys + return (TRUE); +} + +int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { + + int i; + + SetMrelInfo summary; + SetMrelInfoInit (&summary, FALSE); + + pthread_attr_t attr; + pthread_attr_init (&attr); + pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED); + + pthread_t *threads; + ALLOCATE (threads, pthread_t, NTHREADS); + + ThreadInfo *threadinfo; + ALLOCATE (threadinfo, ThreadInfo, NTHREADS); + + // launch N worker threads + for (i = 0; i < NTHREADS; i++) { + threadinfo[i].entry = i; + threadinfo[i].state = THREAD_RUN; + threadinfo[i].catalog = catalog; + threadinfo[i].Ncatalog = Ncatalog; + threadinfo[i].flatcorr = flatcorr; + SetMrelInfoInit (&threadinfo[i].summary, FALSE); + pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]); + } + pthread_attr_destroy (&attr); + + // wait until all threads have finished + while (1) { + int allDone = TRUE; + for (i = 0; i < NTHREADS; i++) { + if (threadinfo[i].state == THREAD_RUN) allDone = FALSE; + } + if (allDone) { + break; + } + usleep (500000); + } + + // all threads are done, free the threads array and grab the info + free (threads); + + // report stats & summary from the threads + for (i = 0; i < NTHREADS; i++) { + fprintf (stderr, "setMrel thread %d : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", + i, + threadinfo[i].summary.Nfew, + threadinfo[i].summary.Nbad, + threadinfo[i].summary.Ncal, + threadinfo[i].summary.Nmos, + threadinfo[i].summary.Ngrid, + threadinfo[i].summary.Nsys); + SetMrelInfoAccum (&summary, &threadinfo[i].summary); + } + fprintf (stderr, "total : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + free (threadinfo); + + return TRUE; +} + +void *setMrel_worker (void *data) { + + ThreadInfo *threadinfo = data; + int Nsecfilt = GetPhotcodeNsecfilt (); - Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0; - - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - // XXX accumulate all secfilt values in a single pass? - - int minUbercalDist = 1000; - - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - int Nsec = GetPhotcodeNsec(thisCode); - - /* calculate the average mag in this SEC photcode for a single star */ - - // skip bad stars - if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - - N = 0; - m = catalog[i].averageT[j].measureOffset; - for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { - - // skip measurements that do not match the current photcode - PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode); - if (!code) continue; - if (code->equiv != thisCode) { continue; } - - if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; } - - if (getImageEntry (m, i) < 0) { - Mcal = Mmos = Mgrid = 0; - } else { - Mcal = getMcal (m, i, flatcorr, catalog); - if (isnan(Mcal)) { Ncal ++; continue; } - Mmos = getMmos (m, i); - if (isnan(Mmos)) { Nmos ++; continue; } - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) { Ngrid++; continue; } - } - - Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); - if (isnan(Msys)) { Nsys++; continue; } - list[N] = Msys - Mcal - Mmos - Mgrid; - - int myUbercalDist = getUbercalDist(m,i); - minUbercalDist = MIN(minUbercalDist, myUbercalDist); - - // dlist gives the error, which is used as the weight in WT_MEAN. - // we can modify the resulting weight in a few ways: - // 1) MIN_ERROR guarantees a floor - // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode - // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages - // 4) some reference photcode of some kind can be specified as fixed and have a high weight - dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR); - - // up-weight the ubercal values (or convergence can take a long time...) - if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { - dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR); - } - - // tie down reference photometry if the -refcode (code) option is selected - // eg, -refcode g_SDSS - // this probably makes no sense in the context of multifilter analysis - if (refPhotcode) { - if (code->code == refPhotcode->code) { - // tiny error -> large weight - // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR); - dlist[N] = 0.0001; - } - } - N++; - } - - // when performing the grid analysis, STAR_TOOFEW will be set to 1; - if (N <= STAR_TOOFEW) { /* too few measurements */ - // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].measureT[catalog[i].averageT[j].measureOffset].imageID, N, STAR_TOOFEW); - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW; - Nfew ++; - } else { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW; - } - - liststats (list, dlist, N, &stats); - - catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; - catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma; - catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; - - catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; - } - } - } - fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys); - - return (TRUE); + + SetMrelInfo results; + SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here + + while (1) { + + off_t i = getNextCatalogForThread(threadinfo->Ncatalog); + if (i == -1) { + threadinfo->state = THREAD_DONE; + return NULL; + } + + Catalog *catalog = threadinfo->catalog; + FlatCorrectionTable *flatcorr = threadinfo->flatcorr; + + setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt); + SetMrelInfoAccum (&threadinfo->summary, &results); + } + + SetMrelInfoFree (&results); + return NULL; } Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/args.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/args.c (revision 33796) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/args.c (revision 33797) @@ -74,4 +74,11 @@ } + NTHREADS = 0; + if ((N = get_argument (argc, argv, "-threads"))) { + remove_argument (N, &argc, argv); + NTHREADS = atof (argv[N]); + remove_argument (N, &argc, argv); + } + // XXX for the moment, make this selection manual. it needs to be automatic // based on the state of the SkyTable Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/load_images.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/load_images.c (revision 33796) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/load_images.c (revision 33797) @@ -15,4 +15,5 @@ off_t *LineNumber; struct timeval start, stop; + char *inSubset; SkyTable *sky = NULL; @@ -34,4 +35,8 @@ MARKTIME("read image table: %f sec\n", dtime); + // allocate and init an array to identify the images included in the subset + ALLOCATE (inSubset, char, Nimage); + memset (inSubset, 0, Nimage); + // determine the populated SkyRegions overlapping the requested area @@ -45,5 +50,5 @@ // select the images which overlap the selected sky regions // 'subset' points to a new copy of the data (different from 'image') - subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); } else { @@ -59,5 +64,5 @@ // select the images which overlap the selected sky regions // 'subset' points to a new copy of the data (different from 'image') - subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); @@ -66,5 +71,5 @@ off_t *LineNumberExtra, i; - subsetExtra = select_images (extraList, image, Nimage, &LineNumberExtra, &NsubsetExtra); + subsetExtra = select_images (extraList, image, Nimage, inSubset, &LineNumberExtra, &NsubsetExtra); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); @@ -83,17 +88,17 @@ // select the images which overlap the selected sky regions // 'subset' points to a new copy of the data (different from 'image') - subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); } } + // match chips to mosaics (if applicable) + initMosaics (subset, Nsubset, image, inSubset, Nimage); + MARKTIME("init mosaics: %f sec\n", dtime); + // save the subset of images in the static reference in ImageOps, set up indexes initImages (subset, LineNumber, Nsubset); MARKTIME("init images: %f sec\n", dtime); - // match chips to mosaics (if applicable) - initMosaics (subset, Nsubset); - MARKTIME("init mosaics: %f sec\n", dtime); - return (skylist); } Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot.c (revision 33796) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot.c (revision 33797) @@ -121,5 +121,5 @@ plot_scatter (catalog, Ncatalog, flatcorr); } - setMrel (catalog, Ncatalog, flatcorr); + setMrel (catalog, Ncatalog, flatcorr); // threaded if (PLOTSTUFF) { plot_scatter (catalog, Ncatalog, flatcorr); Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/select_images.c =================================================================== --- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/select_images.c (revision 33796) +++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/select_images.c (revision 33797) @@ -22,5 +22,5 @@ fprintf (stderr, MSG, __VA_ARGS__); } -Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) { +Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage) { Image *image; @@ -125,4 +125,6 @@ double RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); MARKTIME("create sky region coords: %f sec\n", dtime); + + // if we overlap 0,360, RminSkyRegion will be 180 (or so) dsortindex (RmaxSky, index, skylist[0].Nregions); @@ -252,4 +254,5 @@ found_it: image[nimage] = timage[i]; + inSubset[i] = TRUE; /* always allow 'few' images to succeed, if possible (new images / detections may have * been added) */