Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 33651) +++ trunk/Ohana/src/relphot/Makefile (revision 33963) @@ -46,4 +46,5 @@ $(SRC)/relphot_objects.$(ARCH).o \ $(SRC)/select_images.$(ARCH).o \ +$(SRC)/client_logger.$(ARCH).o \ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ @@ -73,4 +74,5 @@ $(SRC)/relphot_objects.$(ARCH).o \ $(SRC)/relphot_client.$(ARCH).o \ +$(SRC)/client_logger.$(ARCH).o \ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ Index: trunk/Ohana/src/relphot/doc/mosaic.txt =================================================================== --- trunk/Ohana/src/relphot/doc/mosaic.txt (revision 33963) +++ trunk/Ohana/src/relphot/doc/mosaic.txt (revision 33963) @@ -0,0 +1,30 @@ + +For 'mosaic' zero points, I define the mosaics by grouping the images +by observation time. I do this by generating a list of the image +start times. I sort this list and then find the unique entries. Each +time defines a new mosaic entry. I then assign images to the mosaics +by matching their start times. + +The above currently happens on the image subset. I need to identify +the images which are NOT in the subset but for which a mosaic has been +constructed -- these are partially accepted images. I then need to be +able to do something with those images and mosaics: either reject all +partial mosaics or accept all images from partial mosaics. The +advantage of the latter is that the choice of accepted images depends +on the catalog / regions loaded -- I need to accept images only if all +overlapping region files are loaded, right? + +So, a few things are needed: + +1) identify the images which are NOT selected by have matching start +times. +2) exclude them from the subset somehow. + +Some questions: + + * the image subset is currently a *copy* of the Image array, not a + view. does it help to change this? (beyond the memory savings?) + + * I'm still a bit unclear on what 'load_images' vs 'reload_images' + does. + Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 33651) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 33963) @@ -3,4 +3,5 @@ # include # include +# include /* # define GRID_V1 */ @@ -26,11 +27,23 @@ float dMcal; float dMsys; - short nFitPhotom; + unsigned short nFitPhotom; short Xm; float secz; float ubercalDist; unsigned int flags; + char skipCal; // if TRUE, this mosaic is incomplete and should not be calibrated Coords coords; } Mosaic; + +typedef enum { + STATS_NONE, + STATS_MEAN, + STATS_MEDIAN, + STATS_WT_MEAN, + STATS_INNER_MEAN, + STATS_INNER_WTMEAN, + STATS_CHI_INNER_MEAN, + STATS_CHI_INNER_WTMEAN +} ListStatsMode; typedef struct { @@ -46,4 +59,5 @@ double total; int Nmeas; + ListStatsMode statmode; } StatType; @@ -79,15 +93,15 @@ /* global variables set in parameter file */ -# define MAX_PATH_LENGTH 1024 -char ImageCat[MAX_PATH_LENGTH]; -char ImageTemplate[MAX_PATH_LENGTH]; -char CatTemplate[MAX_PATH_LENGTH]; -char GSCFILE[MAX_PATH_LENGTH]; +# define DVO_MAX_PATH 1024 +char ImageCat[DVO_MAX_PATH]; +char ImageTemplate[DVO_MAX_PATH]; +char CatTemplate[DVO_MAX_PATH]; +char GSCFILE[DVO_MAX_PATH]; char *CATDIR; char CATMODE[16]; /* raw, mef, split, mysql */ char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ -char CameraConfig[MAX_PATH_LENGTH]; +char CameraConfig[DVO_MAX_PATH]; char CAMERA[64]; /* eg, gpc1 */ -char SKY_TABLE[MAX_PATH_LENGTH]; +char SKY_TABLE[DVO_MAX_PATH]; int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ @@ -111,4 +125,6 @@ int PARALLEL_MANUAL; int PARALLEL_SERIAL; + +int NTHREADS; int VERBOSE; @@ -141,4 +157,5 @@ int PLOTDELAY; int UpdateAverages; +int ApplyOffsets; char *PhotcodeList; @@ -194,5 +211,5 @@ int args PROTO((int argc, char **argv)); int args_client PROTO((int argc, char **argv)); -int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); +int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog, int Ncat)); void clean_images PROTO((void)); void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr)); @@ -236,15 +253,15 @@ 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)); void initialize_client PROTO((int argc, char **argv)); -void initstats PROTO((char *mode)); -int liststats PROTO((double *value, double *dvalue, int N, StatType *stats)); +void liststats_setmode PROTO((StatType *stats, char *strmode)); +int liststats PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats)); Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath)); Catalog *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog)); 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, SkyRegion *region)); int main PROTO((int argc, char **argv)); @@ -273,5 +290,5 @@ int setMmos PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); -void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage)); int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr)); int setMave PROTO((Catalog *catalog, int Ncatalog)); @@ -319,2 +336,6 @@ int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage); ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage); + +int client_logger_init (); +int client_logger_message (char *format,...); + Index: trunk/Ohana/src/relphot/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relphot/src/ConfigInit.c (revision 33651) +++ trunk/Ohana/src/relphot/src/ConfigInit.c (revision 33963) @@ -35,7 +35,7 @@ // force CATDIR to be absolute (so parallel mode will work) char *tmpcatdir = NULL; - ALLOCATE (tmpcatdir, char, MAX_PATH_LENGTH); + ALLOCATE (tmpcatdir, char, DVO_MAX_PATH); GetConfig (config, "CATDIR", "%s", 0, tmpcatdir); - CATDIR = abspath (tmpcatdir, MAX_PATH_LENGTH); + CATDIR = abspath (tmpcatdir, DVO_MAX_PATH); free (tmpcatdir); Index: trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- trunk/Ohana/src/relphot/src/GridOps.c (revision 33651) +++ trunk/Ohana/src/relphot/src/GridOps.c (revision 33963) @@ -527,5 +527,7 @@ double *list, *dlist; float Msys, Mrel, Mcal, Mmos; + StatType stats; + liststats_setmode (&stats, STATMODE); if (!USE_GRID) return; @@ -595,5 +597,5 @@ } - liststats (list, dlist, N, &stats); + liststats (list, dlist, NULL, N, &stats); gridM[i] = stats.mean; gridS[i] = stats.sigma; Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 33651) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 33963) @@ -41,7 +41,7 @@ // relationships between the measure,catalog set and the images: -static IDX_T **MeasureToImage; // image index from measure,catalog : MeasureToImage[cat][meas] = ImageIndex -static IDX_T **ImageToCatalog; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) -static IDX_T **ImageToMeasure; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) +static IDX_T **MeasureToImage = NULL; // image index from measure,catalog : MeasureToImage[cat][meas] = ImageIndex +static IDX_T **ImageToCatalog = NULL; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) +static IDX_T **ImageToMeasure = NULL; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) // MeasureToImage was 'bin' @@ -322,4 +322,6 @@ off_t i; + if (!MeasureToImage) return -1; + i = MeasureToImage[cat][meas]; return (i); @@ -355,4 +357,6 @@ short distance; + if (!MeasureToImage) return -1; + i = MeasureToImage[cat][meas]; if (i == -1) return (1000); @@ -395,5 +399,7 @@ float Msys, Mrel, Mmos, Mgrid, Mflat; double *list, *dlist, *Mlist, *dMlist; + StatType stats; + liststats_setmode (&stats, STATMODE); if (FREEZE_IMAGES) return; @@ -514,5 +520,6 @@ } - liststats (list, dlist, N, &stats); + // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not) + liststats (list, dlist, NULL, N, &stats); image[i].Mcal = stats.mean; image[i].dMcal = stats.error; @@ -523,5 +530,5 @@ // bright end scatter - liststats (Mlist, dMlist, Nbright, &stats); + liststats (Mlist, dMlist, NULL, Nbright, &stats); image[i].dMagSys = stats.sigma; @@ -556,5 +563,4 @@ double *mlist, *slist, *dlist; double MaxOffset, MaxScatter, MedOffset; - StatType stats; if (FREEZE_IMAGES) return; @@ -574,9 +580,14 @@ N++; } - initstats ("MEAN"); - liststats (mlist, dlist, N, &stats); + + // use a straight mean to find the global image statistics (no weighting) + StatType stats; + liststats_setmode (&stats, "MEAN"); + + liststats (mlist, dlist, NULL, N, &stats); MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma); MedOffset = stats.median; - liststats (slist, dlist, N, &stats); + + liststats (slist, dlist, NULL, N, &stats); MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median); fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); @@ -599,5 +610,4 @@ fprintf (stderr, "%d images marked poor\n", Nmark); - initstats (STATMODE); free (mlist); free (slist); @@ -720,6 +730,6 @@ double *list, *dlist; float Mcal, Mmos, Mgrid; + StatType stats; - bzero (&stats, sizeof (StatType)); if (FREEZE_IMAGES) return (stats); @@ -751,5 +761,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -779,5 +791,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -807,5 +821,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -835,5 +851,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); Index: trunk/Ohana/src/relphot/src/ImageSubset.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageSubset.c (revision 33651) +++ trunk/Ohana/src/relphot/src/ImageSubset.c (revision 33963) @@ -92,4 +92,11 @@ } +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage) { @@ -167,8 +174,16 @@ } - gfits_fwrite_header (f, &header); - gfits_fwrite_matrix (f, &matrix); - gfits_fwrite_Theader (f, &theader); - gfits_fwrite_table (f, &ftable); + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for image subset %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for image subset %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for image subset %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for image subset %s\n", filename); gfits_free_header (&header); @@ -177,6 +192,14 @@ gfits_free_table (&ftable); - fclose (f); - fflush (f); + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file image subset %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file image subset %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing image subset file %s\n", filename); return TRUE; Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 33651) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 33963) @@ -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,14 @@ } -/* find mosaic frames (unique time periods & photcode name matches mosaic) */ -void initMosaics (Image *image, off_t Nimage) { +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stopTimer, (void *) NULL); \ + dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +/* 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; @@ -159,5 +167,9 @@ char *pname; + struct timeval startTimer, stopTimer; + if (!MOSAIC_ZEROPT) return; + + gettimeofday (&startTimer, (void *) NULL); /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no @@ -166,10 +178,11 @@ */ - // 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); + MARKTIME("create array of all image obstimes: %f sec\n", dtime); Nmosaic = 0; @@ -178,5 +191,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 +206,5 @@ } Nmosaic ++; + MARKTIME("create subset array of mosaic obstimes: %f sec\n", dtime); // now I have a list of uniq start times, and they are in order @@ -212,4 +227,7 @@ mosaic[i].secz = NAN; mosaic[i].photcode = 0; + mosaic[i].skipCal = FALSE; + + memset (&mosaic[i].coords, 0, sizeof(Coords)); MosaicN_IMAGE[i] = 10; @@ -219,25 +237,51 @@ } - ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs - - // assign each image to a mosaic + int Nskip, Nmark; + Nskip = Nmark = 0; + // find any mosaics (startTimesMosaic) which match unselected images for (i = 0; i < Nimage; i++) { - ImageToMosaic[i] = -1; - + 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 = 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); + /* find a matching mosaic */ + j = findMosaic(startTimesMosaic, Nmosaic, start); + if (j != -1) { + // mark this mosaic as bad + mosaic[j].skipCal = TRUE; + Nmark ++; + } + } + fprintf (stderr, "%d total images, %d overlap skyregion, %d do not overlap, but match overlapping mosaics\n", (int) Nimage, (int) Nskip, (int) Nmark); + MARKTIME("find unselected images matching selected mosaics: %f sec\n", dtime); + + ALLOCATE (ImageToMosaic, off_t, Nsubset); // mosaic to which image belongs + + // assign each image to a mosaic + for (i = 0; i < Nsubset; i++) { + ImageToMosaic[i] = -1; + + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (subset[i].photcode); + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + + start = subset[i].tzero; + stop = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1); + j = findMosaic(startTimesMosaic, Nmosaic, start); if (j == -1) { - fprintf (stderr, "error?\n"); + fprintf (stderr, "programming error? all subset images should belong to a mosaic\n"); abort(); } @@ -268,7 +312,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,11 +322,13 @@ 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; } -// use bisection to find the overlapping mosaic +// use bisection to find the overlapping mosaic (returns exact match) +// startTimes is a sorted, unique list of times +// start might not be in the list off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start) { @@ -290,7 +336,12 @@ // find the last mosaic before start - Nlo = 0; // startTimes[Nlo] guaranteed to be <= start - Nhi = Nmosaic - 1; // startTimes[Nhi] guaranteed to be >= start - while (Nhi - Nlo > 10) { + Nlo = 0; // first valid startTimes value + Nhi = Nmosaic - 1; // last valid startTimes value + + // if start is not in this range, return -1 + if (start < startTimes[Nlo]) return (-1); + if (start > startTimes[Nhi]) return (-1); + + while (Nhi - Nlo > 4) { N = 0.5*(Nlo + Nhi); if (startTimes[N] < start) { @@ -300,9 +351,9 @@ } } - - // check for the matched mosaic starting from Nlo - for (N = Nlo; N < Nmosaic; N++) { - if (start > startTimes[N]) continue; - return (N); + // we now have : startTimes[Nlo] < start <= starTimes[Nhi] + + // find a matched mosaic starting from Nlo, or return -1 + for (N = Nlo; N <= Nhi; N++) { + if (startTimes[N] == start) return N; } return (-1); @@ -394,4 +445,5 @@ } +// XXX : what about mosaics with skipCal == TRUE? void setMcalFinal () { @@ -569,11 +621,97 @@ } +typedef struct { + int Nfew; + int Nbad; + int Ncal; + int Ngrid; + int Nrel; + int Nsys; + int Nskip; + 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->Nskip = 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; + summary->Nskip += results->Nskip; +} + +static int npass_output = 0; + +// 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 (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { - off_t i, j, m, c, n, N, Nmax; - int mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nbright; - float Msys, Mrel, Mcal, Mgrid, Mflat; - double *list, *dlist, *Mlist, *dMlist; - StatType stats; + off_t i, N, Nmax; Image *image; @@ -581,9 +719,12 @@ 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); - - int Nsecfilt = GetPhotcodeNsecfilt (); if (PoorImages) { @@ -597,142 +738,16 @@ Nmax = MAX (Nmax, N_onMosaic[i]); } - ALLOCATE (list, double, Nmax); - ALLOCATE (dlist, double, Nmax); - ALLOCATE (Mlist, double, Nmax); - ALLOCATE (dMlist, double, Nmax); - - Nfew = Nbad = Ncal = Nrel = Ngrid = Nsys = 0; - - for (i = 0; i < Nmosaic; i++) { - - /* on PoorImages run, skip good images */ - if (PoorImages) { - bad = mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); - if (!bad) continue; - } - - // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0: - // we optionally do not recalibrate images with UBERCAL zero points - if (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL) { - mosaic[i].ubercalDist = 0; - // propagate ubercalDist to the images - for (j = 0; j < MosaicN_Image[i]; j++) { - off_t im = MosaicToImage[i][j]; - image[im].ubercalDist = mosaic[i].ubercalDist; - // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist); - } - if (KEEP_UBERCAL) continue; - } - - int minUbercalDist = 1000; - - // number of stars to measure the bright-end scatter - Nbright = 0; - - N = 0; - for (j = 0; j < N_onMosaic[i]; j++) { - - m = MosaicToMeasure[i][j]; - c = MosaicToCatalog[i][j]; - - if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { - Nbad ++; - continue; - } - Mcal = getMcal (m, c, flatcorr, catalog); - if (isnan(Mcal)) { - Ncal++; - continue; - } - Mgrid = getMgrid (m, c); - if (isnan(Mgrid)) { - Ngrid ++; - continue; - } - Mrel = getMrel (catalog, m, c); - if (isnan(Mrel)) { - Nrel ++; - continue; - } - - // 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); - - n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); - if (isnan(Msys)) { - Nsys++; - continue; - } - - 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[i] is all measurements, N is good measurements */ - - /* too few good measurements or too many bad measurements (skip in PoorImages run) */ - if (!PoorImages) { - mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[i]); - if (mark) { - if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[MosaicToImage[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, N_onMosaic[i]); } - mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW; - Nfew ++; - } else { - mosaic[i].flags &= ~ID_IMAGE_PHOTOM_FEW; - } - } - liststats (list, dlist, N, &stats); - if (VERBOSE2 && PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas, N); - - mosaic[i].Mcal = stats.mean; - mosaic[i].dMcal = stats.error; - mosaic[i].nFitPhotom = N; - mosaic[i].Xm = 100.0*log10(stats.chisq); - - plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); - - // bright end scatter - liststats (Mlist, dMlist, Nbright, &stats); - mosaic[i].dMsys = stats.sigma; - - if (mosaic[i].Mcal < -CLOUD_TOLERANCE) { - mosaic[i].Mcal = 0.0; - } - - // 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[i].ubercalDist = minUbercalDist + 1; - - // propagate ubercalDist to the images - for (j = 0; j < MosaicN_Image[i]; j++) { - off_t im = MosaicToImage[i][j]; - image[im].ubercalDist = mosaic[i].ubercalDist; - // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist); - } - } - free (list); - free (dlist); - free (Mlist); - free (dMlist); - - fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Ngrid, Nrel, Nsys); + + SetMmosInfo info; + SetMmosInfoInit (&info, Nmax, TRUE, PoorImages); + + for (i = 0; i < Nmosaic; i++) { + setMmos_mosaic (&mosaic[i], 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) { @@ -744,4 +759,341 @@ } +// 'mosaic' is a pointer to the current mosaic of interest (Nmos) +int setMmos_mosaic (Mosaic *myMosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) { + + Image *imageReal; + off_t j, NimageReal; + + StatType stats; + liststats_setmode (&stats, STATMODE); + + double *list = info->list; + double *dlist = info->dlist; + double *Mlist = info->Mlist; + double *dMlist = info->dMlist; + + assert (Nmos >= 0); + assert (Nmos < Nmosaic); + + imageReal = getimages (&NimageReal, NULL); + + /* on PoorImages run, skip good images */ + if (info->PoorImages) { + int bad = myMosaic[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 (myMosaic[0].flags & ID_IMAGE_PHOTOM_UBERCAL) { + myMosaic[0].ubercalDist = 0; + // propagate ubercalDist to the images + for (j = 0; j < MosaicN_Image[Nmos]; j++) { + off_t im = MosaicToImage[Nmos][j]; + assert (im < NimageReal); + assert (im >= 0); + image[im].ubercalDist = myMosaic[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; + + // do not modify the calibration for mosaics with partial images loaded (skipCal TRUE) + if (myMosaic[0].skipCal) { + info->Nskip ++; + return TRUE; + } + + int testImage = FALSE; + // testImage |= (abs(myMosaic[0].start - 1324104046) < 10); + // testImage |= (abs(myMosaic[0].start - 1324103823) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003245) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003069) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003125) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003300) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003365) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003191) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003014) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003484) < 10); + // testImage |= (abs(myMosaic[0].start - 1323003419) < 10); + // testImage |= (abs(myMosaic[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 ++; + continue; + } + Mcal = getMcal (m, c, flatcorr, catalog); + if (isnan(Mcal)) { + info->Ncal++; + continue; + } + Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) { + info->Ngrid ++; + continue; + } + Mrel = getMrel (catalog, m, c); + if (isnan(Mrel)) { + info->Nrel ++; + continue; + } + + // 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++; + continue; + } + + 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: + assert (N < info->Nmax); + assert (N >= 0); + assert (Nbright < info->Nmax); + assert (Nbright >= 0); + + 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]); } + myMosaic[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 { + myMosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW; + } + } + + liststats (list, dlist, NULL, N, &stats); + if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N); + + myMosaic[0].Mcal = stats.mean; + myMosaic[0].dMcal = stats.error; + myMosaic[0].nFitPhotom = N; + myMosaic[0].Xm = 100.0*log10(stats.chisq); + + if (testImage) { + fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic[0].start, stats.mean, stats.error, myMosaic[0].nFitPhotom); + } + + plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); + + // bright end scatter + liststats (Mlist, dMlist, NULL, Nbright, &stats); + myMosaic[0].dMsys = stats.sigma; + + if (myMosaic[0].Mcal < -CLOUD_TOLERANCE) { + myMosaic[0].Mcal = 0.0; + } + + if (testImage) { + fprintf (stderr, "%f %f : %d %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Xm, pow(10.0, 0.01*myMosaic[0].Xm)); + } + + // 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 + myMosaic[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 = myMosaic[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); + + // each time this function is called, we cycle through the available mosaics + // make sure we start at 0 + nextMosaic = 0;; + + // 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), %d partials skipped\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, + threadinfo[i].info.Nskip); + 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), %d partials skipped\n", + summary.Nfew, + summary.Nbad, + summary.Ncal, + summary.Ngrid, + summary.Nrel, + summary.Nsys, + summary.Nskip); + free (threadinfo); + + npass_output ++; + + 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; +} + +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], 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, @@ -955,4 +1307,5 @@ for (i = 0; i < Nmosaic; i++) { if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].skipCal) continue; list[n] = mosaic[i].Mcal; dlist[n] = 1; @@ -960,5 +1313,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -981,6 +1336,6 @@ n = 0; for (i = 0; i < Nmosaic; i++) { - if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].skipCal) continue; list[n] = mosaic[i].dMcal; dlist[n] = 1; @@ -988,5 +1343,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -1011,4 +1368,5 @@ for (i = 0; i < Nmosaic; i++) { if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].skipCal) continue; N = 0; @@ -1032,5 +1390,7 @@ // fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n", Nmosaic, n); - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -1053,6 +1413,6 @@ n = 0; for (i = 0; i < Nmosaic; i++) { - if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].skipCal) continue; list[n] = pow(10.0, 0.01*mosaic[i].Xm); dlist[n] = 1; @@ -1060,5 +1420,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -1072,5 +1434,4 @@ double *mlist, *slist, *dlist; double MaxOffset, MedOffset, MaxScatter; - StatType stats; if (!MOSAIC_ZEROPT) return; @@ -1085,4 +1446,5 @@ for (i = N = 0; i < Nmosaic; i++) { if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].skipCal) continue; mlist[N] = mosaic[i].Mcal; slist[N] = mosaic[i].dMcal; @@ -1090,9 +1452,13 @@ N++; } - initstats ("MEAN"); - liststats (mlist, dlist, N, &stats); + + StatType stats; + liststats_setmode (&stats, "MEAN"); + + liststats (mlist, dlist, NULL, N, &stats); MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma); MedOffset = stats.median; - liststats (slist, dlist, N, &stats); + + liststats (slist, dlist, NULL, N, &stats); MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median); fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); @@ -1102,4 +1468,7 @@ // if we are keeping ubercal sacrosanct, then we should not be allowed to break them... if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + + if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue; + if (mosaic[i].skipCal) continue; mark = FALSE; @@ -1121,5 +1490,5 @@ fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Nscatter, Noffset); - initstats (STATMODE); + free (mlist); free (slist); Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 33651) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 33963) @@ -2,13 +2,39 @@ 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 Ncode; + int Nsys; + int Nbad; + int Ncal; + int Nmos; + int Ngrid; + double *list; + double *dlist; + double *wlist; + double *aplist; + double *daplist; +} SetMrelInfo; + +enum {THREAD_RUN, THREAD_DONE}; + +typedef struct { + int entry; + int state; + Catalog *catalog; + int Ncatalog; + FlatCorrectionTable *flatcorr; + SetMrelInfo summary; +} ThreadInfo; + +// ** internal functions: +void *setMrel_worker (void *data); +int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr); +int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); +int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure); + +// 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 +47,4 @@ } } - - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); - ALLOCATE (Moffset, double, MAX (1, Nmax)); } @@ -47,132 +69,54 @@ } -int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { - - off_t j, k, m; - int i, N, Nfew, Nsys, Nbad, Ncal, Nmos, Ngrid; - float Msys, Mcal, Mmos, Mgrid; - StatType stats; - - 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); -} - -int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) { - - off_t k; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - off_t m = average[0].measureOffset; - - for (k = 0; k < average[0].Nmeasure; k++, m++) { - fprintf (stderr, "meas: %08x\n", measure[m].dbFlags); - } - - int Ns; - for (Ns = 0; Ns < Nsecfilt; Ns++) { - fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags); - } - return 1; -} - -# define MARK_SKIP_MEAS \ - catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; +void SetMrelInfoInit (SetMrelInfo *results, int allocLists) { + results->Nfew = 0; + results->Ncode = 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); + ALLOCATE (results->wlist, double, Nmax); + } +} + +void SetMrelInfoFree (SetMrelInfo *results) { + free (results->list); + free (results->dlist); + free (results->wlist); +} + +void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) { + summary->Nfew += results->Nfew ; + summary->Ncode += results->Ncode ; + 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); +} // setMrel and setMrelOutput are extremely similar, but have slightly different implications: @@ -180,411 +124,495 @@ // * setMrelOutput skips stars for which there are too few good measurements // * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass - -// setMave is also similar to the above. but: -// * setMave is called by relphot -update-objects -// * setMave excludes all detections with (PSF_QF < 0.85), setMrelOutput allows these for PASS > 2 -// * setMave updates 2MASS average flags -// * setMave updates average EXT flags (PS1 and 2MASS) +// * setMrelOutput updates 2MASS average flags +// * setMrelOutput updates average EXT flags (PS1 and 2MASS) + +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++) { + // pass == -1 for anything other than the final pass + setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); + SetMrelInfoAccum (&summary, &results); + } + if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + + SetMrelInfoFree (&results); + + return (TRUE); +} + int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) { - int i, N; - off_t j, k, m, Nmax; - float Msys, Mcal, Mmos, Mgrid, Map; - double *list, *dlist, *aplist, *daplist; + int i; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + SetMrelInfo summary, results; + SetMrelInfoInit (&summary, FALSE); + SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist + ALLOCATE (results.aplist, double, Nmax); + ALLOCATE (results.daplist, double, Nmax); + + for (i = 0; i < Ncatalog; i++) { + setMrel_catalog (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options + SetMrelInfoAccum (&summary, &results); + } + if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + + SetMrelInfoFree (&results); + free (results.aplist); + free (results.daplist); + 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); + + // each time this function is called, we cycle through the available catalogs. + // make sure we start at 0 + nextCatalog = 0;; + + // 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++) { + if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", + i, + threadinfo[i].summary.Ncode, + 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); + } + if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, 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 (); + + 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; + + // pass == -1 for anything other than the final pass + setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt); + SetMrelInfoAccum (&threadinfo->summary, &results); + } + + SetMrelInfoFree (&results); + return NULL; +} + +# define SKIP_THIS_MEAS(REASON) { \ + catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + if (catalog[Nc].measure) { \ + catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + } \ + results->REASON ++; \ + continue; } + + +// setMrel_catalog is used in 3 different contexts: + +// * during the main relphot iterations: +// ** operations only apply to measureT / averageT +// ** we are applying the analysis to the bright subset catalog +// ** we skip any stars found to be bad (STAR_BAD) + +// * during the final pass +// ** operations are applied to all objects +// ** we have special tests for PS1 extended data, 2MASS extended & good data, & synthetic photometry + +// * during the stand-alone '-averages' mode (relphot_objects) +// ** no image data is loaded +// ** we only use the provided calibration (measure.Mcal) +// ** we skip outlier rejection of measurements + +// 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, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { + + off_t j, k, m; + int N; + float Msys, Mcal, Mmos, Mgrid; + StatType stats, apstats; - int Nsec; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - /* Nmeasure is now different, need to reallocate */ - Nmax = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure); - } - } - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); - - ALLOCATE (aplist, double, MAX (1, Nmax)); - ALLOCATE (daplist, double, MAX (1, Nmax)); - - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - - if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) { - fprintf (stderr, "test obj\n"); - print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure); - } + liststats_setmode (&stats, STATMODE); + liststats_setmode (&apstats, STATMODE); + + double *list = results->list; + double *dlist = results->dlist; + double *wlist = results->wlist; + double *aplist = results->aplist; + double *daplist = results->daplist; + + SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays + + int isSetMrelFinal = (pass >= 0); + + for (j = 0; j < catalog[Nc].Naverage; j++) { + // XXX accumulate all secfilt values in a single pass? + + // option for a test print + if (FALSE && (catalog[Nc].average[j].objID == 0x46a4) && (catalog[Nc].average[j].catID == 0xf40e)) { + fprintf (stderr, "test obj\n"); + print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure); + } + + int GoodPS1 = FALSE; + int Good2MASS = FALSE; + int Galaxy2MASS = FALSE; + + int NextPS1 = 0; + int NpsfPS1 = 0; + + 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 */ + + /* star/photcodes already calibrated */ + if ( isSetMrelFinal && catalog[Nc].found[Nsecfilt*j+Nsec]) continue; + + // skip bad stars + if (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue; + + int Ncode = 0; + int Next = 0; + int haveSynth = FALSE; + + int forceSynth = FALSE; + int forceSynthEntry = -1; + + int haveUbercal = FALSE; int minUbercalDist = 1000; - - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - Nsec = GetPhotcodeNsec(thisCode); - - /* star/photcodes already calibrated */ - if (catalog[i].found[Nsecfilt*j+Nsec]) continue; - - int Ncode = 0; - int Next = 0; - int haveSynth = FALSE; - int haveUbercal = FALSE; - int isRefPhot = FALSE; - - 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; } - Ncode ++; - - if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - - if (getImageEntry (m, i) < 0) { - // these detetions have no image (eg, ref values such as 2MASS) - Mcal = Mmos = Mgrid = 0; - isRefPhot = TRUE; - } else { - Mcal = getMcal (m, i, flatcorr, catalog); - if (isnan(Mcal)) { - MARK_SKIP_MEAS; - continue; - } - Mmos = getMmos (m, i); - if (isnan(Mmos)) { - MARK_SKIP_MEAS; - continue; - } - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) { - MARK_SKIP_MEAS; - continue; - } - } - - Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); - if (isnan(Msys)) { - if (!isRefPhot) { - MARK_SKIP_MEAS; - } - continue; - } - list[N] = Msys - Mcal - Mmos - Mgrid; - - Map = PhotAper (&catalog[i].measure[m]); + + 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; } + Ncode ++; + + if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); + + if (getImageEntry (m, Nc) < 0) { + // measurements without an image are either external reference photometry or + // data for which the associated image has not been loaded (probably because of + // overlaps). Msys + measure.Mcal is our best guess of the true magnitude + Mmos = Mgrid = 0; + Mcal = catalog[Nc].measureT[m].Mcal; // check that this is zero for loaded REF value + } else { + Mcal = getMcal (m, Nc, flatcorr, catalog); + if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); + Mmos = getMmos (m, Nc); + if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); + Mgrid = getMgrid (m, Nc); + if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); + } + + Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]); + if (isnan(Msys)) SKIP_THIS_MEAS(Nsys); + + list[N] = Msys - Mcal - Mmos - Mgrid; + + int myUbercalDist = getUbercalDist(m,Nc); + minUbercalDist = MIN(minUbercalDist, myUbercalDist); + + if (isSetMrelFinal) { + float Map = PhotAper (&catalog[Nc].measure[m]); aplist[N] = Map - Mcal - Mmos - Mgrid; - int myUbercalDist = getUbercalDist(m,i); - minUbercalDist = MIN(minUbercalDist, myUbercalDist); - - // count the extended detections - if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) { - if (!isnan(catalog[i].measure[m].Map)) { - if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { + // special options for PS1 data + if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) { + // count the extended detections + if (!isnan(catalog[Nc].measure[m].Map)) { + if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) { Next ++; + NextPS1 ++; + } else { + NpsfPS1 ++; } } } // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013) - if ((catalog[i].measure[m].photcode >= 2011) && (catalog[i].measure[m].photcode <= 2013)) { - if (catalog[i].measure[m].photFlags & 0x00c00000) { + if ((catalog[Nc].measure[m].photcode >= 2011) && (catalog[Nc].measure[m].photcode <= 2013)) { + if (catalog[Nc].measure[m].photFlags & 0x00c00000) { Next ++; + Galaxy2MASS = TRUE; } - if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) { - // detections without one of these bits should only be used in PASS_1 - MARK_SKIP_MEAS; - continue; - } - } - // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars) - if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) { - if (pass < 4) { - MARK_SKIP_MEAS; - continue; - } - haveSynth = TRUE; - } - - // 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...) (XXX make this optional?) - if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { - dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR); - haveUbercal = TRUE; - } - - // 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++; - } - if (N < 1) continue; - - for (m = 0; m < N; m++) { - daplist[m] = dlist[m]; - } - - // XXX force WT_MEAN or MEAN here? - liststats (list, dlist, N, &stats); - liststats (aplist, daplist, N, &apstats); - catalog[i].found[Nsecfilt*j+Nsec] = TRUE; - switch (pass) { - case 0: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0; - break; - case 1: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1; - break; - case 2: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2; - break; - case 3: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3; - break; - case 4: - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4; - break; - } - - if (haveSynth) { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH; - } - if (haveUbercal) { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL; - } - - /* use sigma or error in dM for output? */ - catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; - catalog[i].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean; - catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma); - catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; - catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode; - catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas; - - catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80; - catalog[i].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20; - catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; - - if ((Next > 0) && (Next > 0.5*N)) { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT; - } - } - } - } - - free (list); - free (dlist); - - free (aplist); - free (daplist); - return (TRUE); -} - -/* grab Nsec for named photcode */ -# define NAMED_PHOTCODE_NSEC(MY_NSEC,NAME) \ - short MY_NSEC = -1; \ - { \ - PhotCode *code = GetPhotcodebyName (NAME); \ - if (code) { \ - MY_NSEC = GetPhotcodeNsec (code->equiv); \ - } } - -// For each average object, set the average mags based on existing equiv photometry. -// NOTE: this function operates on the real Measure & Average structures, not the -// MeasureTiny & AverageTiny structures -// NOTE: this function is called on the remote machine -- make sure any (global) options -// are passed to the relphot_client program -int setMave (Catalog *catalog, int Ncatalog) { - - off_t j, k, m, Nmax; - int i, Ns, Nsecfilt, N, Nc; - float Msys; - double *list, *dlist; - StatType stats; - DVOAverageFlags flagBits; - - flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT; - - // pre-allocate a list for stats purposes - Nmax = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure); - } - } - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); - - Nsecfilt = GetPhotcodeNsecfilt (); - - // we want to raise some bits on the 2MASS (JHK) secfilt flags, if we have 2MASS data - NAMED_PHOTCODE_NSEC (Nsec_J, "2MASS_J"); - NAMED_PHOTCODE_NSEC (Nsec_H, "2MASS_H"); - NAMED_PHOTCODE_NSEC (Nsec_K, "2MASS_K"); - - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - - // update average photometry for each of the average filters - - // XXX Note that this would be faster if we had an array of results and accumulated - // them in a single pass - - int Next = 0; - - for (Ns = 0; Ns < Nsecfilt; Ns++) { - - PhotCode *code = GetPhotcodebyNsec (Ns); - Nc = code[0].code; - - N = 0; - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; - if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue; - - Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - if (isnan(Msys)) continue; - - // reject POOR detections (PSF_QF < 0.85) and count extended detections - // XXX only apply this filter for psphot data from GPC1 for now... - if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) { - if (catalog[i].measure[m].psfQual < 0.85) continue; - if (!isnan(catalog[i].measure[m].Map)) { - if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { - Next ++; + if (pass == 0) { + if (catalog[Nc].measure[m].photFlags & 0x00000007) { + Good2MASS = TRUE; + } else { + // detections without one of these bits should only be used in PASS_1 + SKIP_THIS_MEAS(Nbad); + continue; } } } - // XXX make it optional to apply Mcal (do I need an extra field for advisory Mcal?) - list[N] = Msys - catalog[i].measure[m].Mcal; - dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); - N++; - } - if (N < 1) continue; - - liststats (list, dlist, N, &stats); - - /* use sigma or error in dM for output? */ - catalog[i].secfilt[Nsecfilt*j+Ns].M = stats.mean; - catalog[i].secfilt[Nsecfilt*j+Ns].dM = MAX (stats.error, stats.sigma); - catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; - catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N; - catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas; - - catalog[i].secfilt[Nsecfilt*j+Ns].M_80 = 1000 * stats.Upper80; - catalog[i].secfilt[Nsecfilt*j+Ns].M_20 = 1000 * stats.Lower20; - - if ((Next > 0) && (Next > 0.5*N)) { - catalog[i].secfilt[Nsecfilt*j+Ns].flags |= ID_SECF_OBJ_EXT; - } - } - - // update average flags based on the detection stats. - - int Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000]) - int goodPS1 = FALSE; // true if any PS1 measurements have psfQual > 0.85 - int good2MASS = FALSE; // true if 2MASS J measurements have significant detections - int nEXT = 0; - int nPSF = 0; // number of PS1 PSF vs EXT measurements - int have2MASS = FALSE; - - // count, flag good and extended detections - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - - // PS1 data : - if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) { - if (catalog[i].measure[m].psfQual > 0.85) { - goodPS1 = TRUE; - if (!isnan(catalog[i].measure[m].Map)) { - if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { - nEXT ++; - } else { - nPSF ++; + // Blindly accepth the SYNTH mags if we are above saturation, otherwise, + // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars) + if ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) { + // something of a hack: force object to use synth values if synth mags >> + // saturation (3pi instrumental mags < -15) + float MaxMagForceSynth = NAN; + switch (catalog[Nc].measure[m].photcode) { + case 3001: + MaxMagForceSynth = 13.64; + break; + case 3002: + MaxMagForceSynth = 13.76; + break; + case 3003: + MaxMagForceSynth = 13.74; + break; + case 3004: + MaxMagForceSynth = 12.94; + break; + case 3005: + MaxMagForceSynth = 12.01; + break; + } + if (catalog[Nc].measureT[m].M < MaxMagForceSynth) { + forceSynth = TRUE; + forceSynthEntry = N; + } else { + if (pass < 4) { + SKIP_THIS_MEAS(Nbad); + continue; } + haveSynth = TRUE; } } } - - // 2MASS data J-band flags - if (catalog[i].measure[m].photcode == 2011) { - // only need to do this once (always have JHK triplet; galaxy flag is same for all 3) - have2MASS = TRUE; - if (catalog[i].measure[m].photFlags & 0x00c00000) { - Galaxy2MASS = TRUE; + + // dlist gives the error per measurement, wlist gives the weight + // we can modify the error and 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); + wlist[N] = 1.0; + + // up-weight the ubercal values (or convergence can take a long time...) + if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + wlist[N] = 10.0; + } + + // 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) { + wlist[N] = 100.0; } - if (catalog[i].measure[m].photFlags & 0x00000007) { - good2MASS = TRUE; - if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_0; - } else { - if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_1; - } - } - // 2MASS data H-band flags - if (catalog[i].measure[m].photcode == 2012) { - if (catalog[i].measure[m].photFlags & 0x00000007) { - good2MASS = TRUE; - if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_0; - } else { - if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_1; - } - } - // 2MASS data K-band flags - if (catalog[i].measure[m].photcode == 2013) { - if (catalog[i].measure[m].photFlags & 0x00000007) { - good2MASS = TRUE; - if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_0; - } else { - if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_1; - } - } - } + } + N++; + } + + int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1; + + // when performing the grid analysis, STAR_TOOFEW should be set to 1; + if (N < Nminmeas) { /* 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; + if (Ncode == 0) { + results->Ncode ++; + } else { + results->Nfew ++; + } + continue; + } else { + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW; + } + + if (forceSynth) { + // use the single SYNTH value instead of the other mags here + myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error"); + list[0] = list[forceSynthEntry]; + dlist[0] = dlist[forceSynthEntry]; + wlist[0] = wlist[forceSynthEntry]; + N = 1; + } + liststats (list, dlist, wlist, N, &stats); + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].M = stats.mean; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; + + // when running -averages, we have no information about the images, so we cannot set this + if (minUbercalDist > -1) { + catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; + } + + if (isSetMrelFinal) { + catalog[Nc].found[Nsecfilt*j+Nsec] = TRUE; + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision) + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas; + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20; + + // NOTE : use the modified weight for apmags as well as psf mags + liststats (aplist, daplist, wlist, N, &apstats); + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean; + + // NOTE: for 2MASS measurements, Next should be 1, as should N + if ((Next > 0) && (Next > 0.5*N)) { + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT; + } + + switch (pass) { + case 0: + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0; + GoodPS1 = TRUE; + break; + case 1: + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1; + GoodPS1 = TRUE; + break; + case 2: + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2; + GoodPS1 = TRUE; + break; + case 3: + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3; + break; + case 4: + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4; + break; + } + if (haveSynth) { + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH; + } + if (haveUbercal) { + catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL; + } + } + } + + if (isSetMrelFinal) { + DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT; // we attempt to set a few flags here; reset those bits before trying: - catalog[i].average[j].flags &= ~flagBits; - - // XXX set the secfilt bits? - if (nEXT && (nEXT > nPSF)) { - catalog[i].average[j].flags |= ID_OBJ_EXT; - } - if (goodPS1) { - catalog[i].average[j].flags |= ID_OBJ_GOOD; + catalog[Nc].average[j].flags &= ~flagBits; + + if (NextPS1 && (NextPS1 > NpsfPS1)) { + catalog[Nc].average[j].flags |= ID_OBJ_EXT; + } + if (GoodPS1) { + catalog[Nc].average[j].flags |= ID_OBJ_GOOD; } if (Galaxy2MASS) { - catalog[i].average[j].flags |= ID_OBJ_EXT_ALT; - if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_SECF_OBJ_EXT; - if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_SECF_OBJ_EXT; - if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_SECF_OBJ_EXT; - } - if (good2MASS) { - catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT; - } - } - } - - free (list); - free (dlist); + catalog[Nc].average[j].flags |= ID_OBJ_EXT_ALT; + } + if (Good2MASS) { + catalog[Nc].average[j].flags |= ID_OBJ_GOOD_ALT; + } + } + } return (TRUE); +} + +int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) { + + off_t k; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + off_t m = average[0].measureOffset; + + for (k = 0; k < average[0].Nmeasure; k++, m++) { + fprintf (stderr, "meas: %08x\n", measure[m].dbFlags); + } + + int Ns; + for (Ns = 0; Ns < Nsecfilt; Ns++) { + fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags); + } + return 1; } @@ -612,4 +640,7 @@ if (isnan(Mgrid)) continue; + // note that measurements for which the image is not selected will not be modified + // (that is a good thing! -- we keep a prior calibration) + // set the output calibration catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid; @@ -630,5 +661,7 @@ double Chisq, MaxScatter, MaxChisq; double *xlist, *slist, *dlist; + StatType stats; + liststats_setmode (&stats, "MEAN"); if (VERBOSE) fprintf (stderr, "marking poor stars\n"); @@ -644,11 +677,5 @@ int Nsecfilt = GetPhotcodeNsecfilt (); - // XX int oldPLOTSTUFF = PLOTSTUFF; - // XX PLOTSTUFF = TRUE; - // XX plot_chisq (catalog, Ncatalog); - // XX PLOTSTUFF = oldPLOTSTUFF; - // eliminate bad stars using the stats for a single secfilt at a time - // XXX DEP replace average.flags with secfilt flags for (Ns = 0; Ns < Nphotcodes; Ns ++) { @@ -669,8 +696,8 @@ } - initstats ("MEAN"); - liststats (xlist, dlist, Ntot, &stats); + liststats (xlist, dlist, NULL, Ntot, &stats); MaxChisq = MAX (STAR_CHISQ, 2*stats.median); - liststats (slist, dlist, Ntot, &stats); + + liststats (slist, dlist, NULL, Ntot, &stats); MaxScatter = MAX (STAR_SCATTER, 2*stats.median); fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq); @@ -696,5 +723,4 @@ } fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave); - initstats (STATMODE); } free (xlist); @@ -702,4 +728,9 @@ free (dlist); } + +// clean_measures examines the stats for a single star. It measures the INNER 50% mean +// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3) +// sigma of the INNER 50% mean. it then flags any measurements which are more than +// NSIGMA_REJECT (5) sigma of the mean # define NSIGMA_CLIP 3.0 @@ -712,5 +743,4 @@ double *tlist, *list, *dlist; float Msys, Mcal, Mmos, Mgrid; - StatType stats; int Ncal, Nmos, Ngrid, Nfew; @@ -732,4 +762,9 @@ /* it makes no sense to mark 3-sigma outliers with <5 measurements */ TOOFEW = MAX (5, STAR_TOOFEW); + + // stats structures for inner and full stats + StatType instats, stats; + liststats_setmode (&instats, "INNER_MEAN"); + liststats_setmode (&stats, "MEAN"); Ndel = Nave = 0; @@ -779,11 +814,10 @@ // calculated mean of inner 50% - initstats ("INNER_MEAN"); - liststats (list, dlist, N, &stats); - stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */ + liststats (list, dlist, NULL, N, &instats); + instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */ // ignore entries > 3sigma from inner mean for (k = m = 0; k < N; k++) { - if (fabs (list[k] - stats.median) < NSIGMA_CLIP*stats.sigma) { + if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) { list[m] = list[k]; m++; @@ -791,6 +825,5 @@ } // recalculate the mean & sigma of the accepted measurements - initstats ("MEAN"); - liststats (list, dlist, m, &stats); + liststats (list, dlist, NULL, m, &stats); stats.sigma = MAX (MIN_ERROR, stats.sigma); @@ -841,12 +874,10 @@ } } - initstats (STATMODE); if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave); + free (list); free (dlist); free (ilist); free (tlist); - - } @@ -898,6 +929,7 @@ } - // fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0); - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -938,5 +970,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); @@ -977,5 +1011,7 @@ } - liststats (list, dlist, n, &stats); + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); free (list); free (dlist); Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 33651) +++ trunk/Ohana/src/relphot/src/args.c (revision 33963) @@ -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 @@ -300,4 +307,10 @@ } + ApplyOffsets = FALSE; + if ((N = get_argument (argc, argv, "-apply-offsets"))) { + remove_argument (N, &argc, argv); + ApplyOffsets = TRUE; + } + if (UpdateAverages && (argc == 1)) return TRUE; if (argc != 2) relphot_usage (); @@ -331,12 +344,4 @@ } if (!HOSTDIR) relphot_client_usage(); - - // CATDIR[0] = 0; - // if ((N = get_argument (argc, argv, "-catdir"))) { - // remove_argument (N, &argc, argv); - // strcpy (CATDIR, argv[N]); - // remove_argument (N, &argc, argv); - // } - // if (!CATDIR[0]) relphot_client_usage(); IMAGES = NULL; // used in -update mode Index: trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/bcatalog.c (revision 33651) +++ trunk/Ohana/src/relphot/src/bcatalog.c (revision 33963) @@ -3,5 +3,5 @@ int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog); -int bcatalog (Catalog *subcatalog, Catalog *catalog) { +int bcatalog (Catalog *subcatalog, Catalog *catalog, int Ncat) { off_t i, j, offset; @@ -74,4 +74,5 @@ subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map = NAN; subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT; subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm = NAN; subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20 = NAN_S_SHORT; @@ -158,7 +159,15 @@ subcatalog[0].measureT[Nmeasure].averef = Naverage; if (RESET) { - if (!KEEP_UBERCAL) { - subcatalog[0].measureT[Nmeasure].Mcal = 0; - subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + // only reset Mcal for measures with a matching image + // do not reset Mcal for ubercal images unless explicitly requested + if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + if (!KEEP_UBERCAL) { + subcatalog[0].measureT[Nmeasure].Mcal = 0.0; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + } + } else { + if (getImageEntry (Nmeasure, Ncat) >= 0) { + subcatalog[0].measureT[Nmeasure].Mcal = 0.0; + } } subcatalog[0].measureT[Nmeasure].dbFlags &= ~PHOTOM_FLAGS; Index: trunk/Ohana/src/relphot/src/client_logger.c =================================================================== --- trunk/Ohana/src/relphot/src/client_logger.c (revision 33963) +++ trunk/Ohana/src/relphot/src/client_logger.c (revision 33963) @@ -0,0 +1,35 @@ +# include "relphot.h" + +// I'm getting unlogged errors and failures. I need a log ouput that the clients can +// write independent of the master + +static FILE *logfile = NULL; +int client_logger_init () { + + char filename[DVO_MAX_PATH]; + + snprintf (filename, DVO_MAX_PATH, "%s/log.rlpc.XXXXXX", HOSTDIR); + + int fd = mkstemp (filename); + if (fd == -1) { + fprintf (stderr, "failed to open client logger, exiting\n"); + exit (50); + } + + logfile = fdopen (fd, "w"); + if (!logfile) { + fprintf (stderr, "failed to fdopen client logger, exiting\n"); + exit (51); + } + return TRUE; +} + +int client_logger_message (char *format,...) { + + va_list argp; + + va_start (argp, format); + vfprintf (logfile, format, argp); + va_end (argp); + return TRUE; +} Index: trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- trunk/Ohana/src/relphot/src/global_stats.c (revision 33651) +++ trunk/Ohana/src/relphot/src/global_stats.c (revision 33963) @@ -1,3 +1,9 @@ # include "relphot.h" + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stopTimer, (void *) NULL); \ + dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { @@ -5,8 +11,10 @@ StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD; - initstats ("MEAN"); + // struct timeval startTimer, stopTimer; fprintf (stderr, "\n"); fprintf (stderr, "STATS median mean sigma min max Nmeas\n"); + + // gettimeofday (&startTimer, (void *) NULL); int Ns; @@ -17,7 +25,11 @@ int seccode = photcodes[Ns][0].code; + // MARKTIME("start star stats: %f sec\n", dtime); stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr); + // MARKTIME("done star N: %f sec\n", dtime); stX = statsStarX (catalog, Ncatalog, Nsec); + // MARKTIME("done star X: %f sec\n", dtime); stS = statsStarS (catalog, Ncatalog, Nsec); + // MARKTIME("done star S: %f sec\n", dtime); fprintf (stderr, " --- stats for %s ---\n", photcodes[Ns][0].name); @@ -46,8 +58,4 @@ fprintf (stderr, "dMcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas); fprintf (stderr, "chisq mosaic: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas); - - - initstats (STATMODE); - } Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 33651) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 33963) @@ -12,11 +12,25 @@ PhotcodeList = strcreate (argv[1]); photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes - } + } else { + char tmpline1[256]; + int Ns; + Nphotcodes = GetPhotcodeNsecfilt (); + ALLOCATE (photcodes, PhotCode *, Nphotcodes); + ALLOCATE (PhotcodeList, char, 256); + for (Ns = 0; Ns < Nphotcodes; Ns++) { + photcodes[Ns] = GetPhotcodebyNsec (Ns); + if (Ns > 0) { + snprintf (tmpline1, 256, "%s,%s", PhotcodeList, photcodes[Ns][0].name); + } else { + snprintf (tmpline1, 256, "%s", photcodes[Ns][0].name); + } + strcpy (PhotcodeList, tmpline1); + } + } + if (USE_GRID && (Nphotcodes > 1)) { fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n"); exit (1); } - - initstats (STATMODE); IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; Index: trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- trunk/Ohana/src/relphot/src/liststats.c (revision 33651) +++ trunk/Ohana/src/relphot/src/liststats.c (revision 33963) @@ -1,30 +1,25 @@ # include "relphot.h" -enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, - M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN}; +void liststats_setmode (StatType *stats, char *strmode) { -static int statmode; + stats->statmode = -1; + if (!strcmp (strmode, "MEAN")) { stats->statmode = STATS_MEAN; return; } + if (!strcmp (strmode, "MEDIAN")) { stats->statmode = STATS_MEDIAN; return; } + if (!strcmp (strmode, "WT_MEAN")) { stats->statmode = STATS_WT_MEAN; return; } + if (!strcmp (strmode, "INNER_MEAN")) { stats->statmode = STATS_INNER_MEAN; return; } + if (!strcmp (strmode, "INNER_WTMEAN")) { stats->statmode = STATS_INNER_WTMEAN; return; } + if (!strcmp (strmode, "CHI_INNER_MEAN")) { stats->statmode = STATS_CHI_INNER_MEAN; return; } + if (!strcmp (strmode, "CHI_INNER_WTMEAN")) { stats->statmode = STATS_CHI_INNER_WTMEAN; return; } -void initstats (char *mode) { - - statmode = -1; - if (!strcmp (mode, "MEAN")) statmode = M_MEAN; - if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN; - if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN; - if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN; - if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN; - if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN; - if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN; - - if (statmode == -1) { - fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode); - exit (1); - } + fprintf (stderr, "ERROR: invalid stats mode: %s\n", strmode); + exit (1); } -int liststats (double *value, double *dvalue, int N, StatType *stats) { +int liststats (double *value, double *dvalue, double *weight, int N, StatType *stats) { - int i, ks, ke, Nm; - double Mo, dMo, M, dM, X2, dS, *chi; + int i, ks, ke; + double Mo, dMo, M, dM, Nm, X2, dS, R, W, *chi; + + myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set"); ke = ks = dMo = 0; @@ -37,66 +32,116 @@ if (N < 1) return (FALSE); - dsortpair (value, dvalue, N); - stats[0].median = value[(int)(0.5*N)]; - stats[0].min = value[0]; - stats[0].max = value[N-1]; int N80 = MIN (N-1, 0.8*N); int N20 = MAX (0, 0.2*N); + + if (weight) { + dsortthree (value, dvalue, weight, N); + } else { + dsortpair (value, dvalue, N); + } + + // these values do not depend on the errors or weighting scheme + stats[0].median = value[(int)(0.5*N)]; + stats[0].min = value[0]; + stats[0].max = value[N-1]; stats[0].Upper80 = value[N80]; stats[0].Lower20 = value[N20]; - switch (statmode) { - case M_MEDIAN: - ks = 0; - ke = N; - Mo = stats[0].median; - Nm = N; - goto chisq; - break; - case M_MEAN: - case M_WT_MEAN: - ks = 0; - ke = N; - break; - case M_INNER_MEAN: - case M_INNER_WTMEAN: - case M_CHI_INNER_MEAN: - case M_CHI_INNER_WTMEAN: - ks = 0.25*N + 0.50; - ke = 0.75*N + 0.25; - if (N <= 3) { + switch (stats->statmode) { + case STATS_MEDIAN: ks = 0; ke = N; - } - break; + Mo = stats[0].median; + Nm = N; + goto chisq; + break; + case STATS_MEAN: + case STATS_WT_MEAN: + ks = 0; + ke = N; + break; + case STATS_INNER_MEAN: + case STATS_INNER_WTMEAN: + case STATS_CHI_INNER_MEAN: + case STATS_CHI_INNER_WTMEAN: + ks = 0.25*N + 0.50; + ke = 0.75*N + 0.25; + if (N <= 3) { + ks = 0; + ke = N; + } + break; + case STATS_NONE: + myAbort ("undefined stats"); } - if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) { + // for these two modes, I need a vector of the chi-square contribution + // I'm actually just using chisq to get the correct sorting order + if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) { ALLOCATE (chi, double, N); for (i = 0; i < N; i++) { chi[i] = (value[i] - stats[0].median) / dvalue[i]; } - dsortthree (chi, value, dvalue, N); + if (weight) { + dsortfour (chi, value, dvalue, weight, N); + } else { + dsortthree (chi, value, dvalue, N); + } free (chi); } - /* calculating the per-star offset based on the weighted average */ - M = dM = Nm = 0; - if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) { - for (i = ks; i < ke; i++) { - M += value[i] / SQ (dvalue[i]); - dM += 1.0 / SQ (dvalue[i]); - Nm ++; - } - Mo = M / dM; - dMo = sqrt (1.0 / dM); + int WeightedMean = FALSE; + WeightedMean |= (stats->statmode == STATS_WT_MEAN); + WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN); + WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN); + + /* calculating the per-star offset based on the desired weighting scheme */ + M = dM = Nm = W = R = 0; + if (weight) { + // the weight value is multiplied by whichever nominal weighting scheme is provided + // thus the user should set weight to 1.0 for nominal weight, and 100 for heavy weight (or so) + // and 0.01 for under-weight + if (WeightedMean) { + for (i = ks; i < ke; i++) { + M += value[i] * weight[i] / SQ(dvalue[i]); + W += weight[i] / SQ(dvalue[i]); + dM += SQ (weight[i] / dvalue[i]); + R += weight[i] / SQ(dvalue[i]); + Nm += 1.0; + } + Mo = M / W; + dMo = sqrt (dM) / R; + } else { + for (i = ks; i < ke; i++) { + M += value[i] * weight[i]; + W += weight[i]; + dM += SQ (weight[i] * dvalue[i]); + R += weight[i]; + Nm += 1.0; + } + Mo = M / W; + dMo = sqrt (dM) / R; + } } else { - for (i = ks; i < ke; i++) { - M += value[i]; - dM += SQ (dvalue[i]); - Nm ++; - } - Mo = M / (double) Nm; - dMo = sqrt (dM / (double) Nm); + // NULL weight vector is supplied, revert to standard form (above reverts if weight[i] == 1) + if (WeightedMean) { + // weighted by inverse-variance + for (i = ks; i < ke; i++) { + M += value[i] / SQ (dvalue[i]); + dM += 1.0 / SQ (dvalue[i]); + Nm += 1.0; + } + Mo = M / dM; + dMo = sqrt (1.0 / dM); + } else { + // pure un-weighted + for (i = ks; i < ke; i++) { + M += value[i]; + dM += SQ (dvalue[i]); + Nm += 1.0; + } + Mo = M / Nm; + dMo = sqrt (dM) / Nm; + } } @@ -111,5 +156,5 @@ } X2 = X2 / (Nm - 1); - dS = sqrt (dS / Nm); + dS = sqrt (dS / (Nm - 1)); stats[0].mean = Mo; @@ -122,2 +167,11 @@ } +// we could define the weight to be the only scale factor: +// \mu = \sum (value_i * weight_i) / \sum (weight_i) +// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) +// R = \sum (weight_i^2) + +// or, we could define the weight to be a scale factor times the inverse error: +// \mu = \sum (value_i * weight_i / sigma_i) / \sum (weight_i) +// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) +// R = \sum (weight_i^2) Index: trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/load_catalogs.c (revision 33651) +++ trunk/Ohana/src/relphot/src/load_catalogs.c (revision 33963) @@ -68,5 +68,5 @@ // select only the brighter stars - bcatalog (&catalog[i], &tcatalog); + bcatalog (&catalog[i], &tcatalog, i); dvo_catalog_unlock (&tcatalog); dvo_catalog_free (&tcatalog); @@ -123,5 +123,5 @@ // ensure that the paths are absolute path names - char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH); + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); free (table->hosts[i].pathname); table->hosts[i].pathname = tmppath; @@ -190,5 +190,10 @@ } if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { - HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + // if any one of the remote jobs fails, we should fail and exit + int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + if (!status) { + fprintf (stderr, "at least one remote client job failed to load data, exiting\n"); + exit (3); + } } Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 33651) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 33963) @@ -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, NULL); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); } else { @@ -54,10 +59,7 @@ skylist = SkyListByPatch (sky, -1, &subregion); - subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0; - SkyList *extraList = SkyListByPatch (sky, -1, &subregion); - // 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, &subregion); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); @@ -66,5 +68,8 @@ off_t *LineNumberExtra, i; - subsetExtra = select_images (extraList, image, Nimage, &LineNumberExtra, &NsubsetExtra); + subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0; + SkyList *extraList = SkyListByPatch (sky, -1, &subregion); + + subsetExtra = select_images (extraList, image, Nimage, inSubset, &LineNumberExtra, &NsubsetExtra, &subregion); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); @@ -83,17 +88,26 @@ // 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, region); MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); } } +# if (1) + FILE *ftest = fopen ("skydump.dat", "w"); + int i; + for (i = 0; i < skylist[0].Nregions; i++) { + fprintf (ftest, "%s : %f %f : %f %f\n", skylist[0].regions[i][0].name, skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, skylist[0].regions[i][0].Dmax); + } + fclose (ftest); +# endif + + // 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: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 33651) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 33963) @@ -86,5 +86,6 @@ TIMESTAMP(time4); - setMrelFinal (&catalog, flatcorr); + initMrel (&catalog, 1); + setMrelFinal (&catalog, flatcorr, FALSE); TIMESTAMP(time5); @@ -104,4 +105,8 @@ freeGridBins (1); TIMESTAMP(time7); + + if (hostID) { + client_logger_message ("updated catalog file %s\n", catalog.filename); + } } @@ -119,15 +124,24 @@ int reload_catalogs_parallel (SkyList *sky) { - off_t Nimage; - ImageSubset *image = getimages_subset (&Nimage); - - // write out the subset table of image information + // name of image subset file: char imageFile[512]; snprintf (imageFile, 512, "%s/Images.subset.dat", CATDIR); - if (!ImageSubsetSave (imageFile, image, Nimage)) { - fprintf (stderr, "failed to write image subset\n"); - exit (1); - } - free (image); + + // the ApplyOffsets option re-uses an existing Images.subset.dat file + if (!ApplyOffsets) { + off_t Nimage; + ImageSubset *image = getimages_subset (&Nimage); + if (!ImageSubsetSave (imageFile, image, Nimage)) { + fprintf (stderr, "failed to write image subset\n"); + exit (1); + } + free (image); + } else { + // load the current sky table (layout of all SkyRegions) + SkyTable *skytable = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + // XXX kind of hackish... + ALLOCATE (sky, SkyList, 1); + strcpy (sky->hosts, skytable->hosts); + } // now launch the relphot_client jobs to the parallel hosts @@ -140,5 +154,5 @@ // ensure that the paths are absolute path names - char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH); + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); free (table->hosts[i].pathname); table->hosts[i].pathname = tmppath; @@ -201,5 +215,9 @@ } if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { - HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + if (!status) { + fprintf (stderr, "at least one remote client job failed to load data, exiting\n"); + exit (3); + } } return (TRUE); Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 33651) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 33963) @@ -23,4 +23,14 @@ if (UpdateAverages) { relphot_objects (0, NULL); + exit (0); + } + + if (ApplyOffsets) { + // re-run this step from a previous attempt (assumes an existing Images.subset.dat file) + if (!PARALLEL) { + fprintf (stderr, "-apply-offsets only makes sense in an parallel context\n"); + exit (2); + } + reload_catalogs (NULL, NULL, 0, NULL); exit (0); } @@ -121,5 +131,5 @@ plot_scatter (catalog, Ncatalog, flatcorr); } - setMrel (catalog, Ncatalog, flatcorr); + setMrel (catalog, Ncatalog, flatcorr); // threaded if (PLOTSTUFF) { plot_scatter (catalog, Ncatalog, flatcorr); @@ -128,4 +138,5 @@ setMmos (catalog, FALSE, flatcorr); setMgrid (catalog, flatcorr); + MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); if (PLOTSTUFF) { @@ -139,8 +150,14 @@ // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); // if (i % 6 == 1) rationalize_images (); - if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE, flatcorr); - if (i % 6 == 3) clean_stars (catalog, Ncatalog); - if (i % 6 == 5) clean_mosaics (); - if (i % 6 == 5) clean_images (); + + // NOTE : in the past, I was not iterating enough before cleaning. make sure we do + // at least 8 loops first -- that should get the systematic errors down to the ~1% + // level, even in cases where we have an even split between photometric data and + // data with 1 mag of extinction. + + if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); + if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); + if ((i > 8) && (i % 8 == 5)) clean_mosaics (); + if ((i > 8) && (i % 8 == 5)) clean_images (); // if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); @@ -148,5 +165,5 @@ // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); - global_stats (catalog, Ncatalog, flatcorr); + if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr); MARKTIME("-- finished loop %d: %f sec\n", i, dtime); } Index: trunk/Ohana/src/relphot/src/relphot_client.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_client.c (revision 33651) +++ trunk/Ohana/src/relphot/src/relphot_client.c (revision 33963) @@ -19,4 +19,5 @@ // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc) initialize_client (argc, argv); + client_logger_init (); // load the current sky table (layout of all SkyRegions) @@ -33,4 +34,5 @@ exit (2); } + client_logger_message ("loaded sky table and defined patch\n"); switch (MODE) { @@ -42,9 +44,14 @@ exit (2); } + client_logger_message ("loaded catalogs\n"); + BrightCatalog *bcatalog = BrightCatalogMerge (catalog, Ncatalog); + client_logger_message ("generated subset data\n"); + if (!BrightCatalogSave (BCATALOG, bcatalog)) { fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); exit (2); } + client_logger_message ("generated subset table\n"); break; } @@ -58,4 +65,5 @@ exit (2); } + client_logger_message ("loaded image subset data\n"); // save the available image information in the static array in ImageOps.c @@ -66,6 +74,8 @@ snprintf (flatcorrFile, 1024, "%s/flatcorr.fits", CATDIR); FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrFile, VERBOSE); + client_logger_message ("loaded flat-field correction data\n"); reload_catalogs (skylist, flatcorr, HOST_ID, HOSTDIR); + client_logger_message ("updated catalogs\n"); break; } @@ -73,4 +83,5 @@ case MODE_UPDATE_OBJECTS: { relphot_objects (HOST_ID, HOSTDIR); + client_logger_message ("updated objects\n"); break; } @@ -80,5 +91,7 @@ abort(); } + client_logger_message ("done with relphot_client\n"); exit (0); } + Index: trunk/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_objects.c (revision 33651) +++ trunk/Ohana/src/relphot/src/relphot_objects.c (revision 33963) @@ -78,4 +78,5 @@ catalog.secfilt[j*Nsecfilt + k].Map = NAN; catalog.secfilt[j*Nsecfilt + k].dM = NAN; + catalog.secfilt[j*Nsecfilt + k].Mstdev = NAN_S_SHORT; catalog.secfilt[j*Nsecfilt + k].Xm = NAN; catalog.secfilt[j*Nsecfilt + k].M_20 = NAN_S_SHORT; @@ -89,8 +90,12 @@ } - setMave (&catalog, 1); + populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE); + + initMrel (&catalog, 1); + setMrelFinal (&catalog, NULL, TRUE); if (!UPDATE) { dvo_catalog_unlock (&catalog); + free_tiny_values(&catalog); dvo_catalog_free (&catalog); continue; @@ -104,4 +109,5 @@ dvo_catalog_save (&catalog, VERBOSE); dvo_catalog_unlock (&catalog); + free_tiny_values(&catalog); dvo_catalog_free (&catalog); } @@ -127,5 +133,5 @@ // ensure that the paths are absolute path names - char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH); + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); free (table->hosts[i].pathname); table->hosts[i].pathname = tmppath; Index: trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- trunk/Ohana/src/relphot/src/select_images.c (revision 33651) +++ trunk/Ohana/src/relphot/src/select_images.c (revision 33963) @@ -13,4 +13,5 @@ } SkyRegionCoords; +void MakeSkyCoordIndex (SkyRegionCoords **skycoords_out, double **RmaxSky_out, off_t **index_out, SkyList *skylist); void dsortindex (double *X, off_t *Y, int N); off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); @@ -22,5 +23,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, SkyRegion *region) { Image *image; @@ -28,14 +29,8 @@ off_t *line_number; int InRange, ecode, found; - double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; + double Ri[5], Di[5], Xi[5], Yi[5]; Coords tcoords; - SkyRegionCoords *skycoords; struct timeval start, stop; - double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion; - - double *RmaxSky; - off_t *index; - if (skylist[0].Nregions < 1) { *Nimage = 0; @@ -54,60 +49,59 @@ strcpy (tcoords.ctype, "RA---TAN"); - ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); - - ALLOCATE (RmaxSky, double, skylist[0].Nregions); - ALLOCATE (index, off_t, skylist[0].Nregions); - - RminSkyRegion = +360.0; - RmaxSkyRegion = -360.0; - DminSkyRegion = +90.0; - DmaxSkyRegion = -90.0; - - D_NIMAGE = 6000; - - // FILE *ftest = fopen ("relphot.dump.dat", "w"); - - /* compare with each region file */ + double RminSkyRegion = region[0].Rmin; + double RmaxSkyRegion = region[0].Rmax; + double DminSkyRegion = region[0].Dmin; + double DmaxSkyRegion = region[0].Dmax; + + double RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); + + // generate a set of Rmin/Rmax values using this RmidSkyRegion for DEC bands of width + // dDecBand (0.1 deg). DecIndex = (Dec + 90) / dDecBand + double *RminBand, *RmaxBand; + float dDecBand = 0.1; + int NDecBands = 180 / dDecBand + 1; + ALLOCATE (RminBand, double, NDecBands); + ALLOCATE (RmaxBand, double, NDecBands); + for (i = 0; i < NDecBands; i++) { + RminBand[i] = +360.0; + RmaxBand[i] = -360.0; + } + + /* compare with each region file */ for (i = 0; i < skylist[0].Nregions; i++) { - - /* we make positional comparisons in the projection of catalog */ - skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin); - skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin); - tcoords.crval1 = skycoords[i].Rc; - tcoords.crval2 = skycoords[i].Dc; - - /* define catalog corners */ - RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords); - RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords); - RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords); - RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords); - skycoords[i].Xc[4] = skycoords[i].Xc[0]; - skycoords[i].Yc[4] = skycoords[i].Yc[0]; - - RmaxSky[i] = skylist[0].regions[i][0].Rmax; - index[i] = i; - - dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]); - dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]); - skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy; - skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy; - skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy; - skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; - skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; - - RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin); - RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax); - DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin); - DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax); - } - RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); - MARKTIME("create sky region coords: %f sec\n", dtime); - - dsortindex (RmaxSky, index, skylist[0].Nregions); - MARKTIME("sort sky coords: %f sec\n", dtime); + int iDecBandMin = (skylist[0].regions[i][0].Dmin + 90.0) / dDecBand; + int iDecBandMax = (skylist[0].regions[i][0].Dmax + 90.0) / dDecBand; + + double RminAlt = ohana_normalize_angle_to_midpoint (skylist[0].regions[i][0].Rmin, RmidSkyRegion); + double RmaxAlt = ohana_normalize_angle_to_midpoint (skylist[0].regions[i][0].Rmax, RmidSkyRegion); + + for (j = iDecBandMin; j <= iDecBandMax; j++) { + RminBand[j] = MIN(RminBand[j], RminAlt); + RmaxBand[j] = MAX(RmaxBand[j], RmaxAlt); + } + } + +# if (1) + // XXX quick test of the dec bands: + FILE *fout = fopen ("dec.bands.dat", "w"); + for (i = 0; i < NDecBands; i++) { + fprintf (fout, "%d %f : %f - %f\n", (int) i, dDecBand * i - 90.0, RminBand[i], RmaxBand[i]); + } + fclose (fout); +# endif + + SkyRegionCoords *skycoords = NULL; + double *RmaxSky = NULL; + off_t *index = NULL; + + if (!USE_BASIC_CHECK) { + MakeSkyCoordIndex (&skycoords, &RmaxSky, &index, skylist); + } if (VERBOSE) fprintf (stderr, "finding images\n"); BuildChipMatch (timage, Ntimage); - MARKTIME("build chip match: %f sec\n", dtime); + MARKTIME("build chip match for %d images: %f sec\n", (int) Ntimage, dtime); + + D_NIMAGE = 6000; nimage = 0; @@ -185,4 +179,29 @@ if (RminImage > RmaxSkyRegion) continue; if (RmaxImage < RminSkyRegion) continue; + + // the above checks are only valid for the outermost region. however, at a given + // declination, the range of RA is limited by the actual catalog boundaries + // check if this image is in range for the Dmax and Dmin locations + + int iDecBandMin = (DminImage + 90.0) / dDecBand; + int iDecBandMax = (DmaxImage + 90.0) / dDecBand; + + // the sky region RA is defined to be 0 - 360.0 + if (RminImage < RminBand[iDecBandMin]) { + if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage)); + continue; + } + if (RminImage < RminBand[iDecBandMax]) { + if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage)); + continue; + } + if (RmaxImage > RmaxBand[iDecBandMin]) { + if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage)); + continue; + } + if (RmaxImage > RmaxBand[iDecBandMax]) { + if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage)); + continue; + } // image overlaps region, keep it @@ -230,4 +249,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) */ @@ -255,13 +275,14 @@ MARKTIME("finish image selection: %f sec\n", dtime); - // fclose (ftest); - if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage); REALLOCATE (image, Image, MAX (nimage, 1)); REALLOCATE (line_number, off_t, MAX (nimage, 1)); - free (skycoords); - free (RmaxSky); - free (index); + + if (!USE_BASIC_CHECK) { + free (skycoords); + free (RmaxSky); + free (index); + } *Nimage = nimage; @@ -379,2 +400,61 @@ return (Nlo); } + +void MakeSkyCoordIndex (SkyRegionCoords **skycoords_out, double **RmaxSky_out, off_t **index_out, SkyList *skylist) { + + off_t i; + SkyRegionCoords *skycoords; + Coords tcoords; + double *RmaxSky; + off_t *index; + double dx, dy; + + // the comparison is made in the catalog local projection. below we set crval1,2 + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + + /* compare with each region file */ + ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); + + ALLOCATE (RmaxSky, double, skylist[0].Nregions); + ALLOCATE (index, off_t, skylist[0].Nregions); + + for (i = 0; i < skylist[0].Nregions; i++) { + + /* we make positional comparisons in the projection of catalog */ + skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin); + skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin); + tcoords.crval1 = skycoords[i].Rc; + tcoords.crval2 = skycoords[i].Dc; + + /* define catalog corners */ + RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords); + RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords); + skycoords[i].Xc[4] = skycoords[i].Xc[0]; + skycoords[i].Yc[4] = skycoords[i].Yc[0]; + + RmaxSky[i] = skylist[0].regions[i][0].Rmax; + index[i] = i; + + dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]); + dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]); + skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy; + skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy; + skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy; + skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; + skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; + } + + dsortindex (RmaxSky, index, skylist[0].Nregions); + + *RmaxSky_out = RmaxSky; + *index_out = index; + *skycoords_out = skycoords; + + return; +} Index: trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 33651) +++ trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 33963) @@ -5,5 +5,5 @@ // output dbFlags values -void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) { +void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) { off_t i; @@ -43,4 +43,5 @@ catalog[0].secfilt[N].Map = NAN; catalog[0].secfilt[N].dM = NAN; + catalog[0].secfilt[N].Mstdev = NAN_S_SHORT; catalog[0].secfilt[N].Xm = NAN; catalog[0].secfilt[N].M_20 = NAN_S_SHORT; @@ -65,7 +66,16 @@ } - if (!KEEP_UBERCAL) { - catalog[0].measure[m].Mcal = 0; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + // only reset Mcal for measures with a matching image + // do not reset Mcal for ubercal images unless explicitly requested + + if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + if (!KEEP_UBERCAL) { + catalog[0].measure[m].Mcal = 0.0; + catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + } + } else { + if (getImageEntry (m, 0) >= 0) { + catalog[0].measure[m].Mcal = 0.0; + } } catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS; @@ -84,10 +94,11 @@ } - clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ + // XXX make this optional? (do not clean for -averages?) + if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ for (i = 0; i < 5; i++) { skip_measurements (catalog, i, flatcorr); /* set ID_MEAS_SKIP for measures to be skipped */ setMrelOutput (catalog, 1, i, flatcorr); /* set Mrel using remaining measures */ } - setMcalOutput (catalog, 1, flatcorr); + if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr); /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */ @@ -248,5 +259,5 @@ } } - if (VERBOSE2) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n", + if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n", pass, Ntot, Ntry, Nskip, Nkeep); }