Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 36528) +++ trunk/Ohana/src/relphot/Makefile (revision 36630) @@ -47,5 +47,16 @@ $(SRC)/relphot_images.$(ARCH).o \ $(SRC)/relphot_objects.$(ARCH).o \ +$(SRC)/relphot_parallel_regions.$(ARCH).o \ +$(SRC)/relphot_parallel_images.$(ARCH).o \ $(SRC)/select_images.$(ARCH).o \ +$(SRC)/assign_images.$(ARCH).o \ +$(SRC)/launch_region_hosts.$(ARCH).o \ +$(SRC)/ImageTable.$(ARCH).o \ +$(SRC)/ImageMagIO.$(ARCH).o \ +$(SRC)/syncfile.$(ARCH).o \ +$(SRC)/share_image_mags.$(ARCH).o \ +$(SRC)/share_mean_mags.$(ARCH).o \ +$(SRC)/MeanMagIO.$(ARCH).o \ +$(SRC)/indexCatalog.$(ARCH).o \ $(SRC)/client_logger.$(ARCH).o \ $(SRC)/setExclusions.$(ARCH).o \ Index: trunk/Ohana/src/relphot/doc/parallel.txt =================================================================== --- trunk/Ohana/src/relphot/doc/parallel.txt (revision 36528) +++ trunk/Ohana/src/relphot/doc/parallel.txt (revision 36630) @@ -1,2 +1,71 @@ + +2014.02.15 + +Nearly done with the relphot mods. some outstanding questions; + + * image vs mosaic? + * how do decide which images need to be shared? + +2014.02.14 + + more relphot notes: + + * I load the catalogs, with a + +2014.02.12 + +I am making progess on the relphot -parallel-regions implementation. some things I need to deal with: + + * assign_images: select mosaics for all images saved for a given host (need to supply them as well) + ** I need to decide on the Image / Mosaic split. If I am calibrating by image, then the image center defines ownership + ** if I am calibrating by exposure, then the mosaic center must define ownership + o determine SkyList covering the images for region host + o tag detections and objects which I own or do not own + o ImageOps.c / matchImage matches detections to images and needs to handle mine / not mine cases + o ditto for matchMosaics + o load_catalogs : I need to make my requested catalogs unique (just add uniquer to load_catalogs.c:130 + o figure out what fields I'm saving in meanmags and how to construct it + -- mag, objID, catID, photcode, (dmag or other stats?) + * from which hosts do I slurp mean mags? + o function to merge a new meanmag array into the existing one + o how to go from (objID,catID) to a given catalog[i].average[j] + +2014.02.06 + +** a single region-level host owns images for which the center lands + in its region. + +** it owns detections which come from images which it owns + +** it owns objects which land in its region + +Extending parallel relphot processing to split the sky (and images) +into regions, each of which runs in parallel at the same time: + +Top Level (relphot -parallel-images): + * define regions of the sky -> hosts + * load images, assign to hosts + * launch region-level jobs on remote hosts + + Region Level (relphot -parallel-images-region) + * load my image subset table + * request objects and detections for my skyregion + * match images & objects, etc + + * update my image parameters + * update my detections + * write out detections + * load detections from my border hosts + * match to my objects + * update my objects + * write out objects + * load objects from my border hosts + * update images + (iterate N times) + * write out image parameters + +Top Level + * read image parameters, update + * update objects 2012.02.13 Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 36528) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 36630) @@ -23,4 +23,13 @@ // # define IDX_T off_t # define IDX_T int + +typedef enum { + MODE_ERROR = 0, + UPDATE_IMAGES, + UPDATE_AVERAGES, + PARALLEL_REGIONS, + PARALLEL_IMAGES, + APPLY_OFFSETS, +} RelphotMode; typedef enum { @@ -108,4 +117,25 @@ typedef struct { + float M; + float dM; + float Xm; + int Nsec; + unsigned int objID; + unsigned int catID; + int photcode; +} MeanMag; + +typedef struct { + float Mcal; + float dMcal; + float dMagSys; + float Xm; + int nFitPhotom; + int flags; + unsigned int imageID; + short ubercalDist; +} ImageMag; + +typedef struct { AverageTiny *average; // array of (minimal) average data MeasureTiny *measure; // array of (minimal) measure data @@ -155,4 +185,9 @@ int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ +// globals for parallel region operations +char *REGION_FILE; +char *IMAGE_TABLE; +int REGION_HOST_ID; + int HOST_ID; char *HOSTDIR; @@ -179,4 +214,6 @@ int PARALLEL_MANUAL; int PARALLEL_SERIAL; + +int PARALLEL_REGIONS_MANUAL; int NTHREADS; @@ -266,5 +303,5 @@ void InterpolateGrid PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords)); off_t *SelectRefMosaic PROTO((Mosaic **refmosaic, off_t *Nimage)); -int args PROTO((int argc, char **argv)); +RelphotMode args PROTO((int argc, char **argv)); int args_client PROTO((int argc, char **argv)); int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog, int Ncat)); @@ -279,4 +316,8 @@ void findImages PROTO((Catalog *catalog, int Ncatalog, int doImageList)); int findMosaics PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); + +void makeMosaics (Image *image, off_t Nimage); +Mosaic *getMosaicForImage (off_t im); +void setMosaicCenters (Image *image, off_t Nimage); void set_db (FITS_DB *in); @@ -315,5 +356,5 @@ 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)); +RelphotMode initialize PROTO((int argc, char **argv)); void initialize_client PROTO((int argc, char **argv)); void liststats_setmode PROTO((StatType *stats, char *strmode)); @@ -426,2 +467,44 @@ int init_synthetic_mags (); int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm); + +int relphot_parallel_regions (); + +int assign_images (FITS_DB *db, RegionHostTable *regionHosts); +int select_images_hostregion (RegionHostTable *hosts, Image *image, off_t Nimage); +int find_host_for_coords (RegionHostTable *regionHosts, double R, double d); +int calculate_image_bounds (Image *image, double *rmin, double *rmax, double *dmin, double *dmax, double Rmid); + +int launch_region_hosts (RegionHostTable *regionHosts); +int strextend (char *input, char *format,...); + +Image *ImageTableLoad(char *filename, off_t *nimage); +int ImageTableSave (char *filename, Image *images, off_t Nimages); + +int indexCatalogs (Catalog *catalog, int Ncatalog); +int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq); + +int check_sync_file (char *filename, int nloop); +int clear_sync_file (char *filename); +int update_sync_file (char *filename, int nloop); +char *make_filename (char *dirname, char *hostname, int hostID, char *tailname); + +int share_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int slurp_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int set_mean_mags (MeanMag *meanmags, AverageTiny *average, SecFilt *secfilt, int Nsec); +MeanMag *merge_mean_mags (MeanMag *target, int *ntarget, MeanMag *source, int Nsource); + +MeanMag *MeanMagLoad(char *filename, off_t *nmeanmags); +int MeanMagSave(char *filename, MeanMag *meanmags, off_t Nmeanmags); + +int share_image_mags (RegionHostTable *regionHosts, int nloop); +int slurp_image_mags (RegionHostTable *regionHosts, int nloop); +int set_image_mags (ImageMag *image_mags, Image *image); +ImageMag *merge_image_mags (ImageMag *target, int *ntarget, ImageMag *source, int Nsource); + +ImageMag *ImageMagLoad(char *filename, off_t *nimage_mags); +int ImageMagSave(char *filename, ImageMag *image_mags, off_t Nimage_mags); + +int markObjects (Catalog *catalog, int Ncatalog); + +int relphot_parallel_images (); +int relphot_parallel_regions (); Index: trunk/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 36528) +++ trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 36630) @@ -101,4 +101,5 @@ measure[i].catID = catID[i]; measure[i].photcode = photcode[i]; + measure[i].myDet = FALSE; } fprintf (stderr, "loaded data for %lld measure\n", (long long) Nrow); @@ -141,4 +142,5 @@ GET_COLUMN(flags, "FLAGS", int); GET_COLUMN(catID, "CAT_ID", int); + GET_COLUMN(objID, "OBJ_ID", int); gfits_free_header (&theader); gfits_free_table (&ftable); @@ -153,4 +155,6 @@ average[i].flags = flags[i]; average[i].catID = catID[i]; + average[i].objID = objID[i]; + average[i].nOwn = 0; } fprintf (stderr, "loaded data for %lld average\n", (long long) Nrow); @@ -162,4 +166,5 @@ free (flags ); free (catID ); + free (objID ); catalog->average = average; @@ -370,4 +375,5 @@ gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ref", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ref", NULL, 1.0, 0.0); // generate the output array that carries the data @@ -381,4 +387,5 @@ int *flags ; ALLOCATE (flags, int, catalog->Naverage); int *catID ; ALLOCATE (catID, int, catalog->Naverage); + int *objID ; ALLOCATE (objID, int, catalog->Naverage); // assign the storage arrays @@ -391,4 +398,5 @@ flags[i] = average[i].flags; catID[i] = average[i].catID; + objID[i] = average[i].objID; } @@ -400,4 +408,5 @@ gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, catalog->Naverage); free (R ); @@ -407,4 +416,5 @@ free (flags ); free (catID ); + free (objID ); gfits_fwrite_Theader (f, &theader); Index: trunk/Ohana/src/relphot/src/ImageMagIO.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 36630) +++ trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 36630) @@ -0,0 +1,221 @@ +# include "relphot.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +ImageMag *ImageMagLoad(char *filename, off_t *nimage_mags) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nimage_mags = 0; + ImageMag *image_mags = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // we read the entire block of data, then extract the columns, then set the image structure values. + // I free the FITS table data after extracting the colums to avoid having 3 copies in memory. + + char type[16]; + + GET_COLUMN (Mcal, "MCAL", float); + GET_COLUMN (dMcal, "MCAL_ERR", float); + GET_COLUMN (dMagSys, "MCAL_SYSERR", float); + GET_COLUMN (Xm, "CHISQ", float); + GET_COLUMN (nFitPhotom, "NFIT", int); + GET_COLUMN (flags, "FLAGS", int); + GET_COLUMN (ubercalDist, "UDIST", int); + GET_COLUMN (imageID, "ID", int); + + // free the memory associated with the FITS files + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + ALLOCATE (image_mags, ImageMag, Nrow); + for (i = 0; i < Nrow; i++) { + image_mags[i].Mcal = Mcal [i]; + image_mags[i].dMcal = dMcal [i]; + image_mags[i].dMagSys = dMagSys [i]; + image_mags[i].Xm = Xm [i]; + image_mags[i].nFitPhotom = nFitPhotom [i]; + image_mags[i].flags = flags [i]; + image_mags[i].ubercalDist = ubercalDist[i]; + image_mags[i].imageID = imageID [i]; + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (Mcal ); + free (dMcal ); + free (dMagSys ); + free (Xm ); + free (nFitPhotom ); + free (flags ); + free (ubercalDist); + free (imageID ); + + *nimage_mags = Nrow; + return image_mags; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int ImageMagSave(char *filename, ImageMag *image_mags, off_t Nimage_mags) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "IMAGE_MAGS"); + + gfits_define_bintable_column (&theader, "E", "MCAL", "cal offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "cal error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR", "systematic error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ", "cal chisq", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "UDIST", "distance to ubercal images", "images", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + float *Mcal ; + float *dMcal ; + float *dMagSys ; + float *Xm ; + int *nFitPhotom ; + int *flags ; + int *ubercalDist ; + int *imageID ; + + // create intermediate storage arrays + ALLOCATE (Mcal , float, Nimage_mags); + ALLOCATE (dMcal , float, Nimage_mags); + ALLOCATE (dMagSys , float, Nimage_mags); + ALLOCATE (Xm , float, Nimage_mags); + ALLOCATE (nFitPhotom , int, Nimage_mags); + ALLOCATE (flags , int, Nimage_mags); + ALLOCATE (ubercalDist , int, Nimage_mags); + ALLOCATE (imageID , int, Nimage_mags); + + // assign the storage arrays + for (i = 0; i < Nimage_mags; i++) { + Mcal [i] = image_mags[i].Mcal ; + dMcal [i] = image_mags[i].dMcal ; + dMagSys [i] = image_mags[i].dMagSys ; + Xm [i] = image_mags[i].Xm ; + nFitPhotom [i] = image_mags[i].nFitPhotom ; + flags [i] = image_mags[i].flags ; + ubercalDist[i] = image_mags[i].ubercalDist; + imageID [i] = image_mags[i].imageID ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MCAL", Mcal , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "MCAL_ERR", dMcal , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "MCAL_SYSERR", dMagSys , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "CHISQ", Xm , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "NFIT", nFitPhotom , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "UDIST", ubercalDist, Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "ID", imageID , Nimage_mags); + + free (Mcal ); + free (dMcal ); + free (dMagSys ); + free (Xm ); + free (nFitPhotom ); + free (flags ); + free (ubercalDist); + free (imageID ); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image_mags file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for image_mags %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for image_mags %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for image_mags %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for image_mags %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file image_mags %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file image_mags %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing image_mags file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 36528) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 36630) @@ -244,7 +244,11 @@ int Nmatch = 0; - for (i = 0; i < Ncatalog; i++) { + for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Nmeasure; j++) { + catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode); + + // skip measurements which do not match one of the requested photcodes ( + // (do we not already exclude in bcatalog -- maybe needed for reload_objects? found = FALSE; for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { @@ -252,4 +256,6 @@ } if (!found) continue; + + // if we match one of our images, myDet gets set to TRUE matchImage (catalog, j, i, doImageList); Nmatch ++; @@ -322,4 +328,5 @@ return; } + catalog[cat].measureT[meas].myDet = TRUE; if (USE_GRID) { Index: trunk/Ohana/src/relphot/src/ImageTable.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageTable.c (revision 36630) +++ trunk/Ohana/src/relphot/src/ImageTable.c (revision 36630) @@ -0,0 +1,62 @@ +# include "relphot.h" + +Image *ImageTableLoad(char *filename, off_t *nimage) { + + int status; + off_t Nimage; + FITS_DB db; + + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, filename, 60.0, LCK_SOFT); + if (!status) { + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + // if the file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + Shutdown ("ERROR: failed to read images"); + } + + *nimage = Nimage; + return image; +} + +int ImageTableSave (char *filename, Image *images, off_t Nimages) { + + int status; + FITS_DB db; + + // setup image table format and lock + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&db, filename, 60.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + // load or create the image table + if (db.dbstate != LCK_EMPTY) { + Shutdown ("image catalog already exists %s", db.filename); + } + + dvo_image_create (&db, GetZeroPoint()); + + /* add the new images and save */ + dvo_image_addrows (&db, images, Nimages); + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/MeanMagIO.c =================================================================== --- trunk/Ohana/src/relphot/src/MeanMagIO.c (revision 36630) +++ trunk/Ohana/src/relphot/src/MeanMagIO.c (revision 36630) @@ -0,0 +1,200 @@ +# include "relphot.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +MeanMag *MeanMagLoad(char *filename, off_t *nmeanmags) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nmeanmags = 0; + MeanMag *meanmags = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // a bit annoying : we read the entire block of data, then extract the columns, then set the image structure values. + // this means I need 3 copies in memory at some point. ugh. + + char type[16]; + + GET_COLUMN (M, "MAG", float); + GET_COLUMN (dM, "MAG_ERR", float); + GET_COLUMN (Xm, "MAG_CHISQ", float); + GET_COLUMN (Nsec, "NSEC", int); + GET_COLUMN (objID, "OBJ_ID", int); + GET_COLUMN (catID, "CAT_ID", int); + + // free the memory associated with the FITS files + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + ALLOCATE (meanmags, MeanMag, Nrow); + for (i = 0; i < Nrow; i++) { + meanmags[i].M = M [i]; + meanmags[i].dM = dM [i]; + meanmags[i].Xm = Xm [i]; + meanmags[i].Nsec = Nsec [i]; + meanmags[i].objID = objID[i]; + meanmags[i].catID = catID[i]; + } + fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); + + free (M ); + free (dM ); + free (Xm ); + free (Nsec ); + free (objID); + free (catID); + + *nmeanmags = Nrow; + return meanmags; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int MeanMagSave(char *filename, MeanMag *meanmags, off_t Nmeanmags) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "MEANMAGS"); + + gfits_define_bintable_column (&theader, "E", "MAG", "mean magnitude", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "mean magnitude error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CHISQ", "mean magnitude chisq", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NSEC", "secfilt sequence", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ID", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ID", NULL, 1.0, 1.0*0x8000); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + float *M, *dM, *Xm; + int *Nsec; + unsigned int *objID, *catID; + + // create intermediate storage arrays + ALLOCATE (M, float, Nmeanmags); + ALLOCATE (dM, float, Nmeanmags); + ALLOCATE (Xm, float, Nmeanmags); + ALLOCATE (Nsec, int, Nmeanmags); + ALLOCATE (objID, unsigned int, Nmeanmags); + ALLOCATE (catID, unsigned int, Nmeanmags); + + // assign the storage arrays + for (i = 0; i < Nmeanmags; i++) { + M [i] = meanmags[i].M ; + dM [i] = meanmags[i].dM ; + Xm [i] = meanmags[i].Xm ; + Nsec [i] = meanmags[i].Nsec ; + objID[i] = meanmags[i].objID; + catID[i] = meanmags[i].catID; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MAG", M, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Xm, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "NSEC", Nsec, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeanmags); + + free (M ); + free (dM ); + free (Xm ); + free (Nsec ); + free (objID); + free (catID); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for meanmags %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanmags %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for meanmags %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for meanmags %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file meanmags %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file meanmags %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing meanmags file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 36528) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 36630) @@ -320,4 +320,157 @@ } +/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes but we + match images by MOSAICNAME. this last point is weak: it forces a single camera at a time. + we can extend the logic to multiple cameras if we make list of MOSAICNAMES (better to assign a camera ID) + */ +void makeMosaics (Image *image, off_t Nimage) { + + off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; + unsigned int start, stop, *startTimes, *startTimesMosaic; + char *pname; + + if (!MOSAIC_ZEROPT) return; + + INITTIME; + + /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no + * entry in the image table that represents this mosaic. Instead, it is an + * internal construct that defines a group of related images + */ + + // generate a list of all image start times + ALLOCATE (startTimes, unsigned int, Nimage); + for (i = 0; i < Nimage; i++) { + startTimes[i] = image[i].tzero; + } + sort_times (startTimes, Nimage); + MARKTIME("create array of all image obstimes: %f sec\n", dtime); + + Nmosaic = 0; + NMOSAIC = 1000; + ALLOCATE (startTimesMosaic, unsigned int, NMOSAIC); + startTimesMosaic[0] = startTimes[0]; + + // generate a list of the unique start times (these define the mosaics) + for (i = 0; i < Nimage; i++) { + if (startTimes[i] < startTimesMosaic[Nmosaic]) { + fprintf (stderr, "error?\n"); + abort(); + } + if (startTimes[i] == startTimesMosaic[Nmosaic]) continue; + Nmosaic ++; + if (Nmosaic >= NMOSAIC) { + NMOSAIC += 1000; + REALLOCATE (startTimesMosaic, unsigned int, NMOSAIC); + } + startTimesMosaic[Nmosaic] = startTimes[i]; + } + Nmosaic ++; + MARKTIME("create array of mosaic obstimes: %f sec\n", dtime); + + // now I have a list of uniq start times, and they are in order + // create the mosaic arrays for these times + ALLOCATE (mosaic, Mosaic, Nmosaic); + + ALLOCATE (MosaicToImage, off_t *, Nmosaic); + ALLOCATE (MosaicN_Image, off_t, Nmosaic); + ALLOCATE (MosaicN_IMAGE, off_t, Nmosaic); + + // init the mosaic array values + for (i = 0; i < Nmosaic; i++) { + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].Mcal = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].dMsys = 0.0; + mosaic[i].Xm = 0.0; + mosaic[i].flags = 0; + mosaic[i].secz = NAN; + mosaic[i].photcode = 0; + mosaic[i].skipCal = FALSE; + + memset (&mosaic[i].coords, 0, sizeof(Coords)); + + MosaicN_IMAGE[i] = 10; + MosaicN_Image[i] = 0; + ALLOCATE (MosaicToImage[i], off_t, MosaicN_IMAGE[i]); + MosaicToImage[i][0] = -1; + } + + ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs + + // assign each image to a mosaic + for (i = 0; i < Nimage; i++) { + ImageToMosaic[i] = -1; + + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (image[i].photcode); + if (!pname) continue; + + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + + start = image[i].tzero; + stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + + j = findMosaic(startTimesMosaic, Nmosaic, start); + if (j == -1) { + fprintf (stderr, "programming error? all image images should belong to a mosaic\n"); + abort(); + } + + // add reference from image to mosaic + ImageToMosaic[i] = j; + + // have we already found this mosaic? + found = (MosaicN_Image[j] > 0); + + /* add image to mosaic image list */ + MosaicToImage[j][MosaicN_Image[j]] = i; + MosaicN_Image[j] ++; + if (MosaicN_Image[j] == MosaicN_IMAGE[j]) { + MosaicN_IMAGE[j] += 10; + REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]); + } + if (found) continue; + + /* a new mosaic, define ranges */ + if (mosaic[j].start != start) { + fprintf (stderr, "error?\n"); + abort(); + } + mosaic[j].stop = stop; + mosaic[j].Mcal = 0.0; + mosaic[j].dMcal = 0.0; + mosaic[j].Xm = 0.0; + mosaic[j].dMsys = image[i].flags; + mosaic[j].flags = image[i].flags; + mosaic[j].secz = image[i].secz; + mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + } + MARKTIME("assign images to mosaic: %f sec\n", dtime); + + // free this or not? + free (MosaicN_IMAGE); + free (startTimes); + free (startTimesMosaic); + + fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic); + + return; +} + +Mosaic *getMosaicForImage (off_t im) { + + if (im < 0) return NULL; + if (!ImageToMosaic) return NULL; + + off_t m = ImageToMosaic[im]; + if (m < 0) return NULL; + if (m >= Nmosaic) return NULL; + + return (&mosaic[m]); +} + // use bisection to find the overlapping mosaic (returns exact match) // startTimes is a sorted, unique list of times @@ -350,4 +503,123 @@ } return (-1); +} + +void setMosaicCenters (Image *image, off_t Nimage) { + + /* find max dR, dD range for all mosaics */ + /* define mosaic.coords to cover dR, dD */ + /* send results to initGridBins */ + + off_t i, j, m, NX, NY, NC, Nc; + double R, D, Rmid, Dmid; + double Mcal, dMcal, Xm; + double *Rc, *Dc; + + NC = 100; + ALLOCATE (Rc, double, NC); + ALLOCATE (Dc, double, NC); + + for (i = 0; i < Nmosaic; i++) { + Nc = 0; + Rmid = Dmid = NAN; + Mcal = dMcal = Xm = 0; + for (j = 0; j < MosaicN_Image[i]; j++) { + m = MosaicToImage[i][j]; + + if (!FindMosaicForImage (image, Nimage, m)) { + if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + + NX = image[m].NX; + NY = image[m].NY; + XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &image[m].coords); + R = ohana_normalize_angle_to_midpoint (R, 180.0); + + // Exclude images with crazy astrometry + // XXX NOTE : this is gpc1-specific + { + double dP1 = hypot(image[m].coords.pc1_1, image[m].coords.pc1_2); + double dP2 = hypot(image[m].coords.pc2_1, image[m].coords.pc2_2); + if (fabs(dP1 - 1.0) > 0.02) continue; + if (fabs(dP2 - 1.0) > 0.02) continue; + + double X00, Y00, X10, Y10, X01, Y01; + XY_to_LM (&X00, &Y00, 0.0, 0.0, &image[m].coords); + XY_to_LM (&X10, &Y10, image[m].NX, 0.0, &image[m].coords); + XY_to_LM (&X01, &Y01, 0.0, image[m].NY, &image[m].coords); + double dS0 = hypot ((X00 - X10), (Y00 - Y10)); + double dS1 = hypot ((X00 - X01), (Y00 - Y01)); + if (dS0 > 6000) continue; + if (dS1 > 6500) continue; + } + + Rc[Nc] = R; + Dc[Nc] = D; + Nc ++; + if (Nc >= NC) { + NC += 100; + REALLOCATE (Rc, double, NC); + REALLOCATE (Dc, double, NC); + } + + Mcal += image[m].Mcal; + dMcal += image[m].dMcal; + Xm += image[m].Xm; + + // for ubercal images, we (elsewhere) keep Mcal frozen + + /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ + image[m].Mcal = 0.0; + image[m].dMcal = NAN; + image[m].Xm = NAN_S_SHORT; + } + + if (Nc > 0) { + dsort (Rc, Nc); + if (Rc[Nc-1] - Rc[0] > 180.0) { + // in our list, Rc is in the range 0.0 to 360.0. + // any mosaic which is close to the 0.0, 360.0 boundary may have some on + // one side or the other. count how many have values more than Rc[0] + 180. + // if more than half are at the large end, re-normalize to that range + int Nbig = 0; + for (j = 1; j < Nc; j++) { + if (Rc[j] - Rc[0] > 180.0) Nbig ++; + } + if (Nbig > 0.5*Nc) { + for (j = 0; j < Nc; j++) { + Rc[j] = ohana_normalize_angle_to_midpoint (Rc[j], 360.0); + } + dsort (Rc, Nc); + } else if (Nbig > 0) { + for (j = 0; j < Nc; j++) { + Rc[j] = ohana_normalize_angle_to_midpoint (Rc[j], 0.0); + } + dsort (Rc, Nc); + } + } + dsort (Dc, Nc); + + Rmid = Rc[(int)(0.5*Nc)]; + Dmid = Dc[(int)(0.5*Nc)]; + } + + strcpy (mosaic[i].coords.ctype, "DEC--TAN"); + mosaic[i].coords.crval1 = Rmid; + mosaic[i].coords.crval2 = Dmid; + mosaic[i].coords.crpix1 = 0.0; + mosaic[i].coords.crpix2 = 0.0; + mosaic[i].coords.cdelt1 = 1.0 / 3600.0; + mosaic[i].coords.cdelt2 = 1.0 / 3600.0; + mosaic[i].coords.pc1_1 = 1.0; + mosaic[i].coords.pc2_2 = 1.0; + mosaic[i].coords.pc1_2 = 0.0; + mosaic[i].coords.pc2_1 = 0.0; + + mosaic[i].Mcal = Mcal / MosaicN_Image[i]; + mosaic[i].dMcal = dMcal / MosaicN_Image[i]; + mosaic[i].Xm = Xm / MosaicN_Image[i]; + } + return; } @@ -527,8 +799,11 @@ if (!MOSAIC_ZEROPT) return (FALSE); + // if we are calibrating by mosaic, redefine myDet == on one of my mosaics int Nmatch = 0; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Nmeasure; j++) { + catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise + if (TimeSelect) { if (catalog[i].measureT[j].t < TSTART) continue; @@ -582,4 +857,6 @@ } + // this measurement is on one of my mosaics, mark it as mine. + catalog[cat].measureT[meas].myDet = TRUE; MeasureToMosaic[cat][meas] = mosID; Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 36528) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 36630) @@ -908,4 +908,43 @@ } +int markObjects (Catalog *catalog, int Ncatalog) { + + int i, n; + off_t j; + + // How strongly do I own this object? + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + int nOwn = 0; + int m = catalog[i].averageT[j].measureOffset; + for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) { + if (!catalog[i].measureT[m+n].myDet) continue; + nOwn ++; + } + catalog[i].averageT[j].nOwn = nOwn; + } + } + return TRUE; +} + +int dumpObjects (char *filename, Catalog *catalog, int Ncatalog) { + + int i, n; + off_t j; + + FILE *ftest = fopen (filename, "w"); + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + int m = catalog[i].averageT[j].measureOffset; + for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) { + fprintf (ftest, "%08x %08x %10.6f %10.6f %3d %1d\n", catalog[i].averageT[j].catID, catalog[i].averageT[j].objID, catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].averageT[j].Nmeasure, catalog[i].measureT[m+n].myDet); + } + } + } + fclose (ftest); + return TRUE; +} + void clean_stars (Catalog *catalog, int Ncatalog) { Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 36528) +++ trunk/Ohana/src/relphot/src/args.c (revision 36630) @@ -1,8 +1,538 @@ # include "relphot.h" -int args (int argc, char **argv) { +RelphotMode args (int argc, char **argv) { int N; double trange; + + /* define time */ + TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &TSTART)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (MODE_ERROR); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_dtime (argv[N], &trange)) { + if (!ohana_str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (MODE_ERROR); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + } + + /* specify portion of the sky */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + /* specify region file by name (eg n0000/0000.00) */ + UserCatalog = NULL; + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserCatalog = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + USE_BASIC_CHECK = FALSE; + if ((N = get_argument (argc, argv, "-basic-image-search"))) { + remove_argument (N, &argc, argv); + USE_BASIC_CHECK = TRUE; + } + + USE_FULL_OVERLAP = TRUE; + if ((N = get_argument (argc, argv, "-sloppy-image-overlap"))) { + remove_argument (N, &argc, argv); + USE_FULL_OVERLAP = FALSE; + } + + SET_MREL_VERSION = 1; + if ((N = get_argument (argc, argv, "-set-mrel-version"))) { + remove_argument (N, &argc, argv); + SET_MREL_VERSION = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE2 = VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + 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 + HOST_ID = 0; + PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + PARALLEL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the remote jobs and waiting for completion, + // relphot will simply list the remote command and wait for the user to signal completion + PARALLEL_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-manual"))) { + PARALLEL = TRUE; // -parallel-manual implies -parallel + PARALLEL_MANUAL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the relphot_client jobs remotely, they are + // run in serial via 'system' + PARALLEL_SERIAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-serial"))) { + if (PARALLEL_MANUAL) { + fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); + exit (1); + } + PARALLEL = TRUE; // -parallel-serial implies -parallel + PARALLEL_SERIAL = TRUE; + remove_argument (N, &argc, argv); + } + + PLOTSTUFF = FALSE; + if ((N = get_argument (argc, argv, "-plot"))) { + PLOTSTUFF = TRUE; + remove_argument (N, &argc, argv); + } + + PLOTDELAY = 500000; + if ((N = get_argument (argc, argv, "-plotdelay"))) { + remove_argument (N, &argc, argv); + PLOTDELAY = 1e6*atof(argv[N]); + PLOTSTUFF = TRUE; // always turn on plotting if i request a plot delay + remove_argument (N, &argc, argv); + } + + if ((N = get_argument (argc, argv, "-outroot"))) { + remove_argument (N, &argc, argv); + OUTROOT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } else { + OUTROOT = strcreate ("relphot"); + } + + strcpy (STATMODE, "WT_MEAN"); + if ((N = get_argument (argc, argv, "-statmode"))) { + remove_argument (N, &argc, argv); + strcpy (STATMODE, argv[N]); + remove_argument (N, &argc, argv); + } + + NLOOP = 8; + if ((N = get_argument (argc, argv, "-n"))) { + remove_argument (N, &argc, argv); + NLOOP = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-nloop"))) { + remove_argument (N, &argc, argv); + NLOOP = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + NGRID = 8; + if ((N = get_argument (argc, argv, "-ngrid"))) { + remove_argument (N, &argc, argv); + NGRID = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + RESET = FALSE; + if ((N = get_argument (argc, argv, "-reset"))) { + remove_argument (N, &argc, argv); + RESET = TRUE; + } + + RESET_ZEROPTS = FALSE; + if ((N = get_argument (argc, argv, "-reset-zpts"))) { + remove_argument (N, &argc, argv); + RESET_ZEROPTS = TRUE; + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + UPDATE_CATFORMAT = NULL; + if ((N = get_argument (argc, argv, "-update-catformat"))) { + remove_argument (N, &argc, argv); + UPDATE_CATFORMAT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + SAVE_IMAGE_UPDATES = TRUE; + if ((N = get_argument (argc, argv, "-skip-image-updates"))) { + remove_argument (N, &argc, argv); + SAVE_IMAGE_UPDATES = FALSE; + } + + MaxDensityUse = FALSE; + if ((N = get_argument (argc, argv, "-max-density"))) { + remove_argument (N, &argc, argv); + MaxDensityValue = atof(argv[N]); + remove_argument (N, &argc, argv); + MaxDensityUse = TRUE; + } + + CLOUD_TOLERANCE = 0.02; + if ((N = get_argument (argc, argv, "-cloud-limit"))) { + remove_argument (N, &argc, argv); + CLOUD_TOLERANCE = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + // XXX should we load a tree from CATDIR by default? + // NOTE: a given catdir needs an appropriate boundary tree + BOUNDARY_TREE = NULL; + if ((N = get_argument (argc, argv, "-boundary-tree"))) { + remove_argument (N, &argc, argv); + BOUNDARY_TREE = strcreate(argv[N]); + load_tess (BOUNDARY_TREE); + remove_argument (N, &argc, argv); + } + + SHOW_PARAMS = FALSE; + if ((N = get_argument (argc, argv, "-params"))) { + remove_argument (N, &argc, argv); + SHOW_PARAMS = TRUE; + } + + PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0; + if ((N = get_argument (argc, argv, "-plrange"))) { + remove_argument (N, &argc, argv); + PlotMmin = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotMmax = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotdMmin = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotdMmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + /* XXX this argument used to do two things: specify the camera name and tell the analysis to + calculate a common zero point for a single mosaic. I've moved the MOSAICNAME concept into the + config system, but need to use this argument to specify that the mosaic zeropoints should be + calculated. */ + MOSAIC_ZEROPT = FALSE; + if ((N = get_argument (argc, argv, "-mosaic"))) { + remove_argument (N, &argc, argv); + MOSAIC_ZEROPT = TRUE; + if (!strcasecmp (MOSAICNAME, "none")) { + fprintf (stderr, "mosaic astrometry selected by MOSAICNAME not defined\n"); + exit (2); + } + } + + FREEZE_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-imfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_IMAGES = TRUE; + } + + FREEZE_MOSAICS = FALSE; + if ((N = get_argument (argc, argv, "-mosfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_MOSAICS = TRUE; + } + + // USE_GRID is not valid for all cases, probably should be its own mode... + USE_GRID = FALSE; + if ((N = get_argument (argc, argv, "-grid"))) { + remove_argument (N, &argc, argv); + USE_GRID = TRUE; + } + + KEEP_UBERCAL = TRUE; + if ((N = get_argument (argc, argv, "-reset-ubercal"))) { + remove_argument (N, &argc, argv); + KEEP_UBERCAL = FALSE; + } + + MIN_ERROR = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = atof (argv[N]); + remove_argument (N, &argc, argv); + /* require MIN_ERROR > 0 */ + } + + AreaSelect = FALSE; + if ((N = get_argument (argc, argv, "-area"))) { + remove_argument (N, &argc, argv); + AreaXmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaXmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaSelect = TRUE; + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; + } + + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = TRUE; + init_synthetic_mags(); + } + + refPhotcode = NULL; + if ((N = get_argument (argc, argv, "-refcode"))) { + remove_argument (N, &argc, argv); + refPhotcode = GetPhotcodebyName (argv[N]); + if (!refPhotcode) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); + exit (1); + } + remove_argument (N, &argc, argv); + } + + REGION_FILE = NULL; + if ((N = get_argument (argc, argv, "-region-hosts"))) { + remove_argument (N, &argc, argv); + REGION_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + REGION_HOST_ID = 0; + if ((N = get_argument (argc, argv, "-region-hostID"))) { + remove_argument (N, &argc, argv); + REGION_HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + RelphotMode mode = MODE_ERROR; + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + mode = UPDATE_IMAGES; + } + if ((N = get_argument (argc, argv, "-averages"))) { + remove_argument (N, &argc, argv); + mode = UPDATE_AVERAGES; + } + if ((N = get_argument (argc, argv, "-apply-offsets"))) { + remove_argument (N, &argc, argv); + mode = APPLY_OFFSETS; + } + IMAGE_TABLE = NULL; + if ((N = get_argument (argc, argv, "-parallel-images"))) { + remove_argument (N, &argc, argv); + mode = PARALLEL_IMAGES; + if (N >= argc) relphot_usage(); + IMAGE_TABLE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + if (!REGION_FILE) relphot_usage(); + } + + PARALLEL_REGIONS_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-regions"))) { + remove_argument (N, &argc, argv); + mode = PARALLEL_REGIONS; + if (!REGION_FILE) relphot_usage(); + if ((N = get_argument (argc, argv, "-parallel-regions-manual"))) { + remove_argument (N, &argc, argv); + PARALLEL_REGIONS_MANUAL = TRUE; + } + } + + switch (mode) { + case UPDATE_AVERAGES: + if (argc != 1) relphot_usage(); + break; + + case UPDATE_IMAGES: + case PARALLEL_IMAGES: + case PARALLEL_REGIONS: + PhotcodeList = strcreate (argv[1]); + photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes + remove_argument (1, &argc, argv); + break; + + default: + fprintf (stderr, "no valid mode selected\n"); + relphot_usage(); + break; + } + if (argc != 1) relphot_usage (); + + return mode; +} + +int args_client (int argc, char **argv) { + + int N; + double trange; + + // by definition, the client is not parallel + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + HOST_ID = 0; + if ((N = get_argument (argc, argv, "-hostID"))) { + remove_argument (N, &argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOST_ID) relphot_client_usage(); + + HOSTDIR = NULL; + if ((N = get_argument (argc, argv, "-hostdir"))) { + remove_argument (N, &argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOSTDIR) relphot_client_usage(); + + IMAGES = NULL; // used in -update mode + BCATALOG = NULL; // used in -load mode + MODE = MODE_NONE; + if ((N = get_argument (argc, argv, "-load"))) { + MODE = MODE_LOAD; + remove_argument (N, &argc, argv); + BCATALOG = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-update-catalogs"))) { + if (MODE) { + fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); + relphot_client_usage(); + } + MODE = MODE_UPDATE; + remove_argument (N, &argc, argv); + IMAGES = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-update-objects"))) { + if (MODE) { + fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); + relphot_client_usage(); + } + MODE = MODE_UPDATE_OBJECTS; + remove_argument (N, &argc, argv); + } + if (!MODE) relphot_client_usage(); + + strcpy (STATMODE, "WT_MEAN"); + if ((N = get_argument (argc, argv, "-statmode"))) { + remove_argument (N, &argc, argv); + strcpy (STATMODE, argv[N]); + remove_argument (N, &argc, argv); + } + + BOUNDARY_TREE = NULL; + if ((N = get_argument (argc, argv, "-boundary-tree"))) { + remove_argument (N, &argc, argv); + BOUNDARY_TREE = strcreate(argv[N]); + load_tess (BOUNDARY_TREE); + remove_argument (N, &argc, argv); + } + + /* specify portion of the sky */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + SET_MREL_VERSION = 1; + if ((N = get_argument (argc, argv, "-set-mrel-version"))) { + remove_argument (N, &argc, argv); + SET_MREL_VERSION = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE2 = VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + RESET = FALSE; + if ((N = get_argument (argc, argv, "-reset"))) { + remove_argument (N, &argc, argv); + RESET = TRUE; + } + + RESET_ZEROPTS = FALSE; + if ((N = get_argument (argc, argv, "-reset-zpts"))) { + remove_argument (N, &argc, argv); + RESET_ZEROPTS = TRUE; + } + + KEEP_UBERCAL = TRUE; + if ((N = get_argument (argc, argv, "-reset-ubercal"))) { + remove_argument (N, &argc, argv); + KEEP_UBERCAL = FALSE; + } /* define time */ @@ -33,153 +563,12 @@ } - /* specify portion of the sky */ - UserPatch.Rmin = 0; - UserPatch.Rmax = 360; - UserPatch.Dmin = -90; - UserPatch.Dmax = +90; - if ((N = get_argument (argc, argv, "-region"))) { - remove_argument (N, &argc, argv); - UserPatch.Rmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Rmax = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmax = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - /* specify region file by name (eg n0000/0000.00) */ - UserCatalog = NULL; - if ((N = get_argument (argc, argv, "-catalog"))) { - remove_argument (N, &argc, argv); - UserCatalog = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - - USE_BASIC_CHECK = FALSE; - if ((N = get_argument (argc, argv, "-basic-image-search"))) { - remove_argument (N, &argc, argv); - USE_BASIC_CHECK = TRUE; - } - - USE_FULL_OVERLAP = TRUE; - if ((N = get_argument (argc, argv, "-sloppy-image-overlap"))) { - remove_argument (N, &argc, argv); - USE_FULL_OVERLAP = FALSE; - } - - SET_MREL_VERSION = 1; - if ((N = get_argument (argc, argv, "-set-mrel-version"))) { - remove_argument (N, &argc, argv); - SET_MREL_VERSION = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - VERBOSE = VERBOSE2 = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { - VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-vv"))) { - VERBOSE2 = VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - - 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 - HOST_ID = 0; - PARALLEL = FALSE; - if ((N = get_argument (argc, argv, "-parallel"))) { - PARALLEL = TRUE; - remove_argument (N, &argc, argv); - } - // this is a test mode : rather than launching the remote jobs and waiting for completion, - // relphot will simply list the remote command and wait for the user to signal completion - PARALLEL_MANUAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-manual"))) { - PARALLEL = TRUE; // -parallel-manual implies -parallel - PARALLEL_MANUAL = TRUE; - remove_argument (N, &argc, argv); - } - // this is a test mode : rather than launching the relphot_client jobs remotely, they are - // run in serial via 'system' - PARALLEL_SERIAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-serial"))) { - if (PARALLEL_MANUAL) { - fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); - exit (1); - } - PARALLEL = TRUE; // -parallel-serial implies -parallel - PARALLEL_SERIAL = TRUE; - remove_argument (N, &argc, argv); - } - - PLOTSTUFF = FALSE; - if ((N = get_argument (argc, argv, "-plot"))) { - PLOTSTUFF = TRUE; - remove_argument (N, &argc, argv); - } - - PLOTDELAY = 500000; - if ((N = get_argument (argc, argv, "-plotdelay"))) { - remove_argument (N, &argc, argv); - PLOTDELAY = 1e6*atof(argv[N]); - PLOTSTUFF = TRUE; // always turn on plotting if i request a plot delay - remove_argument (N, &argc, argv); - } - - if ((N = get_argument (argc, argv, "-outroot"))) { - remove_argument (N, &argc, argv); - OUTROOT = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } else { - OUTROOT = strcreate ("relphot"); - } - - strcpy (STATMODE, "WT_MEAN"); - if ((N = get_argument (argc, argv, "-statmode"))) { - remove_argument (N, &argc, argv); - strcpy (STATMODE, argv[N]); - remove_argument (N, &argc, argv); - } - - NLOOP = 8; - if ((N = get_argument (argc, argv, "-n"))) { - remove_argument (N, &argc, argv); - NLOOP = atof (argv[N]); - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-nloop"))) { - remove_argument (N, &argc, argv); - NLOOP = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - NGRID = 8; - if ((N = get_argument (argc, argv, "-ngrid"))) { - remove_argument (N, &argc, argv); - NGRID = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - RESET = FALSE; - if ((N = get_argument (argc, argv, "-reset"))) { - remove_argument (N, &argc, argv); - RESET = TRUE; - } - - RESET_ZEROPTS = FALSE; - if ((N = get_argument (argc, argv, "-reset-zpts"))) { - remove_argument (N, &argc, argv); - RESET_ZEROPTS = TRUE; - } + + MIN_ERROR = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = atof (argv[N]); + remove_argument (N, &argc, argv); + /* require MIN_ERROR > 0 */ + } UPDATE = FALSE; @@ -196,8 +585,41 @@ } - SAVE_IMAGE_UPDATES = TRUE; - if ((N = get_argument (argc, argv, "-skip-image-updates"))) { - remove_argument (N, &argc, argv); - SAVE_IMAGE_UPDATES = FALSE; + AreaSelect = FALSE; + if ((N = get_argument (argc, argv, "-area"))) { + remove_argument (N, &argc, argv); + AreaXmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaXmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaSelect = TRUE; + } + + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = FALSE; + init_synthetic_mags(); + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; } @@ -210,376 +632,8 @@ } - CLOUD_TOLERANCE = 0.02; - if ((N = get_argument (argc, argv, "-cloud-limit"))) { - remove_argument (N, &argc, argv); - CLOUD_TOLERANCE = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - // XXX should we load a tree from CATDIR by default? - // NOTE: a given catdir needs an appropriate boundary tree - BOUNDARY_TREE = NULL; - if ((N = get_argument (argc, argv, "-boundary-tree"))) { - remove_argument (N, &argc, argv); - BOUNDARY_TREE = strcreate(argv[N]); - load_tess (BOUNDARY_TREE); - remove_argument (N, &argc, argv); - } - - SHOW_PARAMS = FALSE; - if ((N = get_argument (argc, argv, "-params"))) { - remove_argument (N, &argc, argv); - SHOW_PARAMS = TRUE; - } - - PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0; - if ((N = get_argument (argc, argv, "-plrange"))) { - remove_argument (N, &argc, argv); - PlotMmin = atof (argv[N]); - remove_argument (N, &argc, argv); - PlotMmax = atof (argv[N]); - remove_argument (N, &argc, argv); - PlotdMmin = atof (argv[N]); - remove_argument (N, &argc, argv); - PlotdMmax = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - /* XXX this argument used to do two things: specify the camera name and tell the analysis to - calculate a common zero point for a single mosaic. I've moved the MOSAICNAME concept into the - config system, but need to use this argument to specify that the mosaic zeropoints should be - calculated. */ - MOSAIC_ZEROPT = FALSE; - if ((N = get_argument (argc, argv, "-mosaic"))) { - remove_argument (N, &argc, argv); - MOSAIC_ZEROPT = TRUE; - if (!strcasecmp (MOSAICNAME, "none")) { - fprintf (stderr, "mosaic astrometry selected by MOSAICNAME not defined\n"); - exit (2); - } - } - - FREEZE_IMAGES = FALSE; - if ((N = get_argument (argc, argv, "-imfreeze"))) { - remove_argument (N, &argc, argv); - FREEZE_IMAGES = TRUE; - } - - FREEZE_MOSAICS = FALSE; - if ((N = get_argument (argc, argv, "-mosfreeze"))) { - remove_argument (N, &argc, argv); - FREEZE_MOSAICS = TRUE; - } - - USE_GRID = FALSE; - if ((N = get_argument (argc, argv, "-grid"))) { - remove_argument (N, &argc, argv); - USE_GRID = TRUE; - } - - KEEP_UBERCAL = TRUE; - if ((N = get_argument (argc, argv, "-reset-ubercal"))) { - remove_argument (N, &argc, argv); - KEEP_UBERCAL = FALSE; - } - - MIN_ERROR = 0.001; - if ((N = get_argument (argc, argv, "-minerror"))) { - remove_argument (N, &argc, argv); - MIN_ERROR = atof (argv[N]); - remove_argument (N, &argc, argv); - /* require MIN_ERROR > 0 */ - } - - AreaSelect = FALSE; - if ((N = get_argument (argc, argv, "-area"))) { - remove_argument (N, &argc, argv); - AreaXmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaXmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaSelect = TRUE; - } - - ImagSelect = FALSE; - if ((N = get_argument (argc, argv, "-instmag"))) { - remove_argument (N, &argc, argv); - ImagMin = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagMax = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagSelect = TRUE; - } - - DophotSelect = FALSE; - if ((N = get_argument (argc, argv, "-dophot"))) { - remove_argument (N, &argc, argv); - DophotValue = atof (argv[N]); - remove_argument (N, &argc, argv); - DophotSelect = TRUE; - } - - SyntheticPhotometry = FALSE; - if ((N = get_argument (argc, argv, "-synthphot"))) { - remove_argument (N, &argc, argv); - SyntheticPhotometry = TRUE; - init_synthetic_mags(); - } - - refPhotcode = NULL; - if ((N = get_argument (argc, argv, "-refcode"))) { - remove_argument (N, &argc, argv); - refPhotcode = GetPhotcodebyName (argv[N]); - if (!refPhotcode) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); - exit (1); - } - remove_argument (N, &argc, argv); - } - - UpdateAverages = FALSE; - if ((N = get_argument (argc, argv, "-averages"))) { - remove_argument (N, &argc, argv); - UpdateAverages = TRUE; - } - - 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 (); + if ((MODE == MODE_UPDATE_OBJECTS) && (argc == 1)) return TRUE; + if (argc != 2) relphot_client_usage (); return TRUE; } -int args_client (int argc, char **argv) { - - int N; - double trange; - - // by definition, the client is not parallel - PARALLEL = FALSE; - PARALLEL_MANUAL = FALSE; - PARALLEL_SERIAL = FALSE; - - HOST_ID = 0; - if ((N = get_argument (argc, argv, "-hostID"))) { - remove_argument (N, &argc, argv); - HOST_ID = atoi (argv[N]); - remove_argument (N, &argc, argv); - } - if (!HOST_ID) relphot_client_usage(); - - HOSTDIR = NULL; - if ((N = get_argument (argc, argv, "-hostdir"))) { - remove_argument (N, &argc, argv); - HOSTDIR = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if (!HOSTDIR) relphot_client_usage(); - - IMAGES = NULL; // used in -update mode - BCATALOG = NULL; // used in -load mode - MODE = MODE_NONE; - if ((N = get_argument (argc, argv, "-load"))) { - MODE = MODE_LOAD; - remove_argument (N, &argc, argv); - BCATALOG = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-update-catalogs"))) { - if (MODE) { - fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); - relphot_client_usage(); - } - MODE = MODE_UPDATE; - remove_argument (N, &argc, argv); - IMAGES = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-update-objects"))) { - if (MODE) { - fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); - relphot_client_usage(); - } - MODE = MODE_UPDATE_OBJECTS; - remove_argument (N, &argc, argv); - } - if (!MODE) relphot_client_usage(); - - strcpy (STATMODE, "WT_MEAN"); - if ((N = get_argument (argc, argv, "-statmode"))) { - remove_argument (N, &argc, argv); - strcpy (STATMODE, argv[N]); - remove_argument (N, &argc, argv); - } - - BOUNDARY_TREE = NULL; - if ((N = get_argument (argc, argv, "-boundary-tree"))) { - remove_argument (N, &argc, argv); - BOUNDARY_TREE = strcreate(argv[N]); - load_tess (BOUNDARY_TREE); - remove_argument (N, &argc, argv); - } - - /* specify portion of the sky */ - UserPatch.Rmin = 0; - UserPatch.Rmax = 360; - UserPatch.Dmin = -90; - UserPatch.Dmax = +90; - if ((N = get_argument (argc, argv, "-region"))) { - remove_argument (N, &argc, argv); - UserPatch.Rmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Rmax = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmax = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - SET_MREL_VERSION = 1; - if ((N = get_argument (argc, argv, "-set-mrel-version"))) { - remove_argument (N, &argc, argv); - SET_MREL_VERSION = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - VERBOSE = VERBOSE2 = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { - VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-vv"))) { - VERBOSE2 = VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - - RESET = FALSE; - if ((N = get_argument (argc, argv, "-reset"))) { - remove_argument (N, &argc, argv); - RESET = TRUE; - } - - RESET_ZEROPTS = FALSE; - if ((N = get_argument (argc, argv, "-reset-zpts"))) { - remove_argument (N, &argc, argv); - RESET_ZEROPTS = TRUE; - } - - KEEP_UBERCAL = TRUE; - if ((N = get_argument (argc, argv, "-reset-ubercal"))) { - remove_argument (N, &argc, argv); - KEEP_UBERCAL = FALSE; - } - - /* define time */ - TimeSelect = FALSE; - if ((N = get_argument (argc, argv, "-time"))) { - remove_argument (N, &argc, argv); - if (!ohana_str_to_time (argv[N], &TSTART)) { - fprintf (stderr, "ERROR: syntax error\n"); - return (FALSE); - } - remove_argument (N, &argc, argv); - if (!ohana_str_to_dtime (argv[N], &trange)) { - if (!ohana_str_to_time (argv[N], &TSTOP)) { - fprintf (stderr, "ERROR: syntax error\n"); - return (FALSE); - } - } else { - if (trange < 0) { - trange = fabs (trange); - TSTOP = TSTART; - TSTART -= trange; - } else { - TSTOP = TSTART + trange; - } - } - remove_argument (N, &argc, argv); - TimeSelect = TRUE; - } - - - MIN_ERROR = 0.001; - if ((N = get_argument (argc, argv, "-minerror"))) { - remove_argument (N, &argc, argv); - MIN_ERROR = atof (argv[N]); - remove_argument (N, &argc, argv); - /* require MIN_ERROR > 0 */ - } - - UPDATE = FALSE; - if ((N = get_argument (argc, argv, "-update"))) { - remove_argument (N, &argc, argv); - UPDATE = TRUE; - } - - UPDATE_CATFORMAT = NULL; - if ((N = get_argument (argc, argv, "-update-catformat"))) { - remove_argument (N, &argc, argv); - UPDATE_CATFORMAT = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - - AreaSelect = FALSE; - if ((N = get_argument (argc, argv, "-area"))) { - remove_argument (N, &argc, argv); - AreaXmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaXmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaSelect = TRUE; - } - - SyntheticPhotometry = FALSE; - if ((N = get_argument (argc, argv, "-synthphot"))) { - remove_argument (N, &argc, argv); - SyntheticPhotometry = FALSE; - init_synthetic_mags(); - } - - ImagSelect = FALSE; - if ((N = get_argument (argc, argv, "-instmag"))) { - remove_argument (N, &argc, argv); - ImagMin = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagMax = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagSelect = TRUE; - } - - DophotSelect = FALSE; - if ((N = get_argument (argc, argv, "-dophot"))) { - remove_argument (N, &argc, argv); - DophotValue = atof (argv[N]); - remove_argument (N, &argc, argv); - DophotSelect = TRUE; - } - - MaxDensityUse = FALSE; - if ((N = get_argument (argc, argv, "-max-density"))) { - remove_argument (N, &argc, argv); - MaxDensityValue = atof(argv[N]); - remove_argument (N, &argc, argv); - MaxDensityUse = TRUE; - } - - if ((MODE == MODE_UPDATE_OBJECTS) && (argc == 1)) return TRUE; - if (argc != 2) relphot_client_usage (); - - return TRUE; -} - Index: trunk/Ohana/src/relphot/src/assign_images.c =================================================================== --- trunk/Ohana/src/relphot/src/assign_images.c (revision 36630) +++ trunk/Ohana/src/relphot/src/assign_images.c (revision 36630) @@ -0,0 +1,289 @@ +# include "relphot.h" + +// This function generates a subset of the images based on selections. Input db has already +// been loaded with the raw fits table data +int assign_images (FITS_DB *db, RegionHostTable *regionHosts) { + + off_t Nimage; + + INITTIME; + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + Image *image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + MARKTIME("convert image table to internal structure: %f sec\n", dtime); + + // *** NOTE : for the moment, regions must be in the range 0 - 360, -90 - +90 + + // generate the chip match here so we can define the mosaic centers (if needed) + BuildChipMatch (image, Nimage); + MARKTIME("build chip match for %d images: %f sec\n", (int) Nimage, dtime); + + if (MOSAIC_ZEROPT) { + makeMosaics (image, Nimage); + + // center coords and Mcal, dMcal, Xm for the mosaics + setMosaicCenters (image, Nimage); + MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime); + } + + // register the image array with ImageOps.c for later getimageByID calls + initImages (image, NULL, Nimage); + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + // for each regionHost, select images which are contained by the region + // even faster would be to use a tree to get to the real regions... + select_images_hostregion (regionHosts, image, Nimage); + + return TRUE; +} + +# define D_NIMAGE 1000 + +// assign images to the region hosts; at the end, each host will have its list of images +int select_images_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage) { + + int ecode, found; + off_t i, j; + + // INITTIME; + + for (i = 0; i < regionHosts->Nhosts; i++) { + regionHosts->hosts[i].Nimage = 0; + regionHosts->hosts[i].NIMAGE = D_NIMAGE; + ALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE); + } + + for (j = 0; j < Nimage; j++) { + + /* exclude images by photcode */ + ecode = GetPhotcodeEquivCodebyCode (image[j].photcode); + found = FALSE; + int Ns; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) continue; + + /* exclude images by time */ + if (TimeSelect) { + if (image[j].tzero < TSTART) continue; + if (image[j].tzero > TSTOP) continue; + } + + // do not include DIS (PHU-level mosaics) in the output list + if (!strcmp(&image[j].coords.ctype[4], "-DIS")) continue; + + // match the image to its corresponding ASTROMETRIC mosaic (not the same as the Mosaic above) + if (!FindMosaicForImage (image, Nimage, j)) { + if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + + // Exclude images with crazy astrometry + // XXX NOTE : this is gpc1-specific + { + double dP1 = hypot(image[j].coords.pc1_1, image[j].coords.pc1_2); + double dP2 = hypot(image[j].coords.pc2_1, image[j].coords.pc2_2); + if (fabs(dP1 - 1.0) > 0.02) continue; + if (fabs(dP2 - 1.0) > 0.02) continue; + + double X00, Y00, X10, Y10, X01, Y01; + XY_to_LM (&X00, &Y00, 0.0, 0.0, &image[j].coords); + XY_to_LM (&X10, &Y10, image[j].NX, 0.0, &image[j].coords); + XY_to_LM (&X01, &Y01, 0.0, image[j].NY, &image[j].coords); + double dS0 = hypot ((X00 - X10), (Y00 - Y10)); + double dS1 = hypot ((X00 - X01), (Y00 - Y01)); + if (dS0 > 6000) continue; + if (dS1 > 6500) continue; + } + + // use a reference coordinate for each image to assign to hosts + // define image center - note the DIS images (mosaic phu) are special, but we have + // already excluded them above. we also save the image centers for reference + double Xc, Yc; + double Rc, Dc; + + Xc = 0.5*image[j].NX; + Yc = 0.5*image[j].NY; + + XY_to_RD (&Rc, &Dc, Xc, Yc, &image[j].coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); + image[j].RAo = Rc; + image[j].DECo = Dc; + + // NOTE: if we are NOT using mosaic centers, we don't need to keep chips with centers + // on opposite sides of 0,360 together. in which case, it is OK for the range of Rc + // to be 0.0 to 360.0. also NOTE: RAo,DECo are only used for debugging reference. + + if (MOSAIC_ZEROPT) { + // use the coords of the associated mosaic to select + Mosaic *mosaic = getMosaicForImage (j); + Rc = mosaic->coords.crval1; + Dc = mosaic->coords.crval2; + // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the + // chips are located. but, for host assignment, we rationalize to 0.0 - 360.0 + Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); + } + + i = find_host_for_coords (regionHosts, Rc, Dc); + + if (i == -1) continue; + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // image bounds are defined for a range centered on the mosaic center thus, a mosaic + // with center 0.5 will have chips bounds ranging from ~ -1.5 to +2.5 or so, while a mosaic + // with center 359.5 will have chips bounds ranging from ~ 357.5 to 361.5 or so + double Rmin, Rmax, Dmin, Dmax; + calculate_image_bounds (&image[j], &Rmin, &Rmax, &Dmin, &Dmax, Rc); + + host->RminCat = MIN(Rmin, host->RminCat); + host->RmaxCat = MAX(Rmax, host->RmaxCat); + host->DminCat = MIN(Dmin, host->DminCat); + host->DmaxCat = MAX(Dmax, host->DmaxCat); + + // regionHosts needs to have the full outer boundary + // (so reload_catalogs covers the correct region) + regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin); + regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax); + regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin); + regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax); + + // this is a bit memory expensive : I am making a complete copy of the image table here + off_t Nsubset = regionHosts->hosts[i].Nimage; + regionHosts->hosts[i].image[Nsubset] = image[j]; + // regionHosts->hosts[i].line_number[Nsubset] = j; + + regionHosts->hosts[i].Nimage ++; + if (regionHosts->hosts[i].Nimage == regionHosts->hosts[i].NIMAGE) { + regionHosts->hosts[i].NIMAGE += D_NIMAGE; + REALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE); + // REALLOCATE (regionHosts->hosts[i].line_number, off_t, regionHosts->hosts[i].NIMAGE); + } + } + + return TRUE; +} + +double Xf[] = {0.0, 1.0, 0.0, 1.0}; +double Yf[] = {0.0, 0.0, 1.0, 1.0}; + +int calculate_image_bounds (Image *image, double *rmin, double *rmax, double *dmin, double *dmax, double Rmid) { + + int n; + + double Rmin = 360.0; + double Rmax = 0.0; + double Dmin = +90.0; + double Dmax = -90.0; + + // define image corners + for (n = 0; n < 4; n++) { + double Xc, Yc, Rc, Dc; + Xc = Xf[n]*image->NX; + Yc = Yf[n]*image->NY; + XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, Rmid); + + Rmin = MIN (Rmin, Rc); + Rmax = MAX (Rmax, Rc); + Dmin = MIN (Dmin, Dc); + Dmax = MAX (Dmax, Dc); + } + + *rmin = Rmin; + *rmax = Rmax; + *dmin = Dmin; + *dmax = Dmax; + + return TRUE; +} + +int calculate_host_image_bounds (RegionHostTable *regionHosts, double Rmid) { + + int i, n; + off_t j; + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // XXX clear the chip match? + + BuildChipMatch (host->image, host->Nimage); + + double Rmin = 720.0; + double Rmax = -360.0; + double Dmin = +90.0; + double Dmax = -90.0; + + for (j = 0; j < host->Nimage; j++) { + + Image *image = &host->image[j]; + + if (!FindMosaicForImage (host->image, host->Nimage, j)) { + fprintf (stderr, "missing astrometry? programming error?\n"); + abort (); + } + + // define image corners + for (n = 0; n < 4; n++) { + double Xc, Yc, Rc, Dc; + Xc = Xf[n]*image->NX; + Yc = Yf[n]*image->NY; + XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, Rmid); + + Rmin = MIN (Rmin, Rc); + Dmin = MIN (Dmin, Rc); + + Rmax = MAX (Rmax, Rc); + Dmax = MAX (Dmax, Rc); + } + + } + + // RminCat, RmaxCat may extended beyond 0.0 - 360.0 + host->RminCat = Rmin; + host->DminCat = Dmin; + host->RmaxCat = Rmax; + host->DmaxCat = Dmax; + + // regionHosts needs to have the full outer boundary + // (so reload_catalogs covers the correct region) + regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin); + regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax); + regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin); + regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax); + } + return TRUE; +} + + +// Rc is in range 0.0 - 360.0, hosts Rmin,Rmax in range 0.0 - 360.0 +int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc) { + + int i; + + if (isnan(Rc)) return -1; + if (isnan(Dc)) return -1; + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + if (Rc < host->Rmin) continue; + if (Rc >= host->Rmax) continue; + if (Dc < host->Dmin) continue; + if (Dc >= host->Dmax) continue; + + return i; + } + return -1; +} + Index: trunk/Ohana/src/relphot/src/client_logger.c =================================================================== --- trunk/Ohana/src/relphot/src/client_logger.c (revision 36528) +++ trunk/Ohana/src/relphot/src/client_logger.c (revision 36630) @@ -5,13 +5,13 @@ static FILE *logfile = NULL; -int client_logger_init () { +int client_logger_init (char *dirname) { char filename[DVO_MAX_PATH]; - snprintf (filename, DVO_MAX_PATH, "%s/log.rlpc.XXXXXX", HOSTDIR); + snprintf (filename, DVO_MAX_PATH, "%s/log.rlpc.XXXXXX", dirname); int fd = mkstemp (filename); if (fd == -1) { - fprintf (stderr, "failed to open client logger, exiting\n"); + fprintf (stderr, "failed to open client logger %s, exiting\n", filename); exit (50); } @@ -32,4 +32,6 @@ vfprintf (logfile, format, argp); va_end (argp); + + fflush (logfile); return TRUE; } Index: trunk/Ohana/src/relphot/src/help.c =================================================================== --- trunk/Ohana/src/relphot/src/help.c (revision 36528) +++ trunk/Ohana/src/relphot/src/help.c (revision 36630) @@ -2,7 +2,11 @@ void relphot_usage (void) { - fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -region RA RA DEC DEC\n"); - fprintf (stderr, " or: relphot (photcodes) -catalog (name)\n"); - fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); + fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -images\n"); + fprintf (stderr, " or: relphot -averages\n"); + fprintf (stderr, " or: relphot -apply-offsets\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-regions -region-hosts (RegionFile)\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-images (ImageTable) -region-hosts (RegionFile)\n\n"); + fprintf (stderr, " regions: -region RA RA DEC DEC)\n"); + fprintf (stderr, " or: -catalog (name)\n"); fprintf (stderr, " use -h for more usage information\n"); exit (2); @@ -18,7 +22,11 @@ show_help: - fprintf (stderr, "ERROR: USAGE: relphot (photcode) -region RA RA DEC DEC\n"); - fprintf (stderr, " or: relphot (photcode) -catalog (name)\n"); - fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); + fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -images\n"); + fprintf (stderr, " or: relphot -averages\n"); + fprintf (stderr, " or: relphot -apply-offsets\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-regions -region-hosts (RegionFile)\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-images (ImageTable) -region-hosts (RegionFile)\n\n"); + fprintf (stderr, " regions: -region RA RA DEC DEC)\n"); + fprintf (stderr, " or: -catalog (name)\n"); fprintf (stderr, " options: \n"); fprintf (stderr, " -time (start) (stop)\n"); Index: trunk/Ohana/src/relphot/src/indexCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/indexCatalog.c (revision 36630) +++ trunk/Ohana/src/relphot/src/indexCatalog.c (revision 36630) @@ -0,0 +1,64 @@ +# include "relphot.h" + +static int catIDmax = 0; +static int *catIDseq = NULL; +static int *objIDmax = NULL; +static int **objIDseq = NULL; + +int indexCatalogs (Catalog *catalog, int Ncatalog) { + + int i; + off_t j; + + // find the max value of catID + for (i = 0; i < Ncatalog; i++) { + catIDmax = MAX (catalog[i].catID, catIDmax); + } + + ALLOCATE (catIDseq, int, catIDmax + 1); + for (i = 0; i < catIDmax + 1; i++) { + catIDseq[i] = -1; + } + + for (i = 0; i < Ncatalog; i++) { + int catID = catalog[i].catID; + catIDseq[catID] = i; + } + + ALLOCATE (objIDmax, int, Ncatalog); + ALLOCATE (objIDseq, int *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + objIDmax[i] = 0; + for (j = 0; j < catalog[i].Naverage; j++) { + objIDmax[i] = MAX (catalog[i].averageT[j].objID, objIDmax[i]); + } + + ALLOCATE (objIDseq[i], int, objIDmax[i] + 1); + for (j = 0; j < objIDmax[i] + 1; j++) { + objIDseq[i][j] = -1; + } + + for (j = 0; j < catalog[i].Naverage; j++) { + int objID = catalog[i].averageT[j].objID; + objIDseq[i][objID] = j; + } + } + return TRUE; +} + +int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq) { + + if (catID > catIDmax) return FALSE; + + int cat = catIDseq[catID]; + if (cat < 1) return FALSE; + + if (objID > objIDmax[cat]) return FALSE; + + int obj = objIDseq[cat][objID]; + if (obj < 1) return FALSE; + + *catSeq = cat; + *objSeq = obj; + return TRUE; +} Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 36528) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 36630) @@ -1,16 +1,13 @@ # include "relphot.h" -void initialize (int argc, char **argv) { +RelphotMode initialize (int argc, char **argv) { relphot_help (argc, argv); ConfigInit (&argc, argv); - args (argc, argv); + RelphotMode mode = args (argc, argv); + if (!mode) exit (2); - if (!UpdateAverages) { - // load the list of photcodes into the globals (photcodes, Nphotcodes) - // only a single remaining argument in this mode (the list of photcodes, eg g,r,i) - PhotcodeList = strcreate (argv[1]); - photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes - } else { + // UPDATE_AVERAGES always operates on all photcodes? + if (mode == UPDATE_AVERAGES) { char tmpline1[256]; int Ns; @@ -67,4 +64,6 @@ for (B = 0; A == time(NULL); B++); srand48(B); + + return mode; } Index: trunk/Ohana/src/relphot/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 36630) +++ trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 36630) @@ -0,0 +1,117 @@ +# include "relphot.h" +# define DEBUG 0 + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} + +int launch_region_hosts (RegionHostTable *regionHosts) { + + int i; + + // clear the I/O files + for (i = 0; i < regionHosts->Nhosts; i++) { + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.sync"); + truncate (syncfile, 0); + free (syncfile); + + char *fitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.fits"); + truncate (fitsfile, 0); + free (fitsfile); + + char *imsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.sync"); + truncate (imsyncfile, 0); + free (imsyncfile); + + char *imfitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.fits"); + truncate (imfitsfile, 0); + free (imfitsfile); + + char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); + truncate (loopsyncfile, 0); + free (loopsyncfile); + } + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // communication files: + // subset images per host : CATDIR/Image.HOSTNAME.fits + char filename[1024]; + snprintf (filename, 1024, "%s/Image.%d.fits", CATDIR, host->hostID); + + // write the image subset for this host + ImageTableSave (filename, host->image, host->Nimage); + + char command[1024]; + snprintf (command, 1024, "relphot %s -parallel-images %s -region-hosts %s -region-hostID %d -D CATDIR %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f", + PhotcodeList, filename, REGION_FILE, host->hostID, CATDIR, host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR); + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (RESET_ZEROPTS) strextend (command, "-reset-zpts"); + if (!KEEP_UBERCAL) strextend (command, "-reset-ubercal"); + if (DophotSelect) strextend (command, "-dophot %d", DophotValue); + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); + if (SyntheticPhotometry) strextend (command, "-synthphot"); + + if (UPDATE) strextend (command, "-update"); + if (MOSAIC_ZEROPT) strextend (command, "-mosaic"); + if (FREEZE_IMAGES) strextend (command, "-imfreeze"); + if (FREEZE_MOSAICS) strextend (command, "-mosfreeze"); + if (PARALLEL) strextend (command, "-parallel"); + if (PARALLEL_MANUAL) strextend (command, "-parallel-manual"); + if (PARALLEL_SERIAL) strextend (command, "-parallel-serial"); + + // XXX deprecate this if we are happy with the new version + if (SET_MREL_VERSION != 1) strextend (command, "-set-mrel-version %d", SET_MREL_VERSION); + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_REGIONS_MANUAL) continue; + + // launch the job, then wait for it to be done loading catalogs. force the remote + // client to generate the file + char *syncfile = make_filename (CATDIR, host->hostname, host->hostID, "loadcat.sync"); + clear_sync_file (syncfile); + + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", host->hostname, command, host->stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", host->hostname, errorInfo); + exit (1); + } + host->pid = pid; // save for future reference + + // remote client is done, go ahead with next client + check_sync_file (syncfile, 1); + free (syncfile); + } + + if (PARALLEL_REGIONS_MANUAL) { + fprintf (stderr, "run the relphot_client commands above. when these are done, hit return\n"); + getchar(); + } else { + RegionHostTableWaitJobsGetIO (regionHosts, __FILE__, __LINE__, VERBOSE); + } + + int status = TRUE; + for (i = 0; i < regionHosts->Nhosts; i++) { + status = status && (regionHosts->hosts[i].status == 0); + } + + return status; +} Index: trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/load_catalogs.c (revision 36528) +++ trunk/Ohana/src/relphot/src/load_catalogs.c (revision 36630) @@ -89,7 +89,6 @@ fprintf (stderr, "using "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measurements)\n", Nstar, Nstar_total, Nmeas, Nmeas_total); - if (!hostID && (Nstar < 1)) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n"); + if (!hostID && !REGION_HOST_ID && (Nstar < 1)) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n"); // in regular relphot, we shutdown here; in relphot_client, we generate and return an empty table (for consistency) - // XXX consider only returning the populated catalogs @@ -107,4 +106,9 @@ # define DEBUG 1 Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); int Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog @@ -128,5 +132,5 @@ ALLOCATE (table->hosts[i].results, char, 1024); - snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.subset.dat", table->hosts[i].pathname); + snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.%s.dat", table->hosts[i].pathname, uniquer); // options / arguments that can affect relphot_client -load: Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 36528) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 36630) @@ -156,4 +156,8 @@ // load the list of hosts HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: problem with parallel host table\n"); + exit (2); + } if (BOUNDARY_TREE) { Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 36528) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 36630) @@ -3,26 +3,42 @@ int main (int argc, char **argv) { - /* get configuration info, args */ - initialize (argc, argv); + // get configuration info, args + RelphotMode mode = initialize (argc, argv); + if (!mode) exit (2); - /* the object analysis is a separate process iterating over catalogs */ - if (UpdateAverages) { - relphot_objects (0, NULL); - exit (0); + switch (mode) { + case UPDATE_IMAGES: + relphot_images (); + exit (0); + + case UPDATE_AVERAGES: + // take the current set of detections and set the mean magnitudes + relphot_objects (0, NULL); + exit (0); + + case PARALLEL_REGIONS: + // run image updates in parallel across multiple remote machines + relphot_parallel_regions (); + exit (0); + + case PARALLEL_IMAGES: + // operation on the remote machines in the PARALLEL_REGION mode + relphot_parallel_images (); + exit (0); + + case APPLY_OFFSETS: + // 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); + + default: + fprintf (stderr, "ERROR: no valid relphot mode chosen\n"); + exit (2); } - - 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); - } - - relphot_images (); - - exit (0); + fprintf (stderr, "IMPOSSIBLE: skipped out of switch?\n"); + exit (1); } - Index: trunk/Ohana/src/relphot/src/relphot_client.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_client.c (revision 36528) +++ trunk/Ohana/src/relphot/src/relphot_client.c (revision 36630) @@ -19,5 +19,5 @@ // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc) initialize_client (argc, argv); - client_logger_init (); + client_logger_init (HOSTDIR); // load the current sky table (layout of all SkyRegions) Index: trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_images.c (revision 36528) +++ trunk/Ohana/src/relphot/src/relphot_images.c (revision 36630) @@ -244,5 +244,7 @@ if (!UPDATE) exit (0); - /* load catalog data from region files, update Mrel include all data */ + /* Load catalog data from region files, update Mrel include all data. In a parallel + context, this function writes the image parameters as a subset table for the remote + clients */ reload_catalogs (skylist, flatcorr, 0, NULL); MARKTIME("-- updated all catalogs: %f sec\n", dtime); Index: trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 36630) +++ trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 36630) @@ -0,0 +1,165 @@ +# include "relphot.h" + +/* This function is essentially identical to relphot_images, except: + + * load the subset images saved by the master node + * distinguish detections we own (touch our images) and those we don't + * distinguish objects we own (in region) and those we don't + * update the unowned detections for owned objects from neighbor regions + * update the unowned objects for owned detections + */ + +int dumpObjects (char *filename, Catalog *catalog, int Ncatalog); + +int relphot_parallel_images () { + + int i, Ncatalog; + Catalog *catalog = NULL; + + INITTIME; + + client_logger_init (CATDIR); + + // load the RegionTable (UserRegion should not be used at this level) + RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE); + int myHost = regionHosts->index[REGION_HOST_ID]; + RegionHostFindNeighbors (regionHosts, myHost); + client_logger_message ("started parallel images on %s\n", regionHosts->hosts[myHost].hostname); + + // load the subset images belonging to this host + off_t Nimage; + Image *image = ImageTableLoad (IMAGE_TABLE, &Nimage); + if (!image) { + fprintf (stderr, "ERROR loading image %s\n", IMAGE_TABLE); + exit (2); + } + + // once we have read this table, we should remove it for repeat runs + // unlink (IMAGE_TABLE); // XXX a bit risky, add some protection? + + // XXX need to deal with mosaic vs image... + makeMosaics (image, Nimage); + + initImages (image, NULL, Nimage); + + // load the flat correction table (if defined) + char flatcorrfile[256]; + sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // UserPatch.Rmin,Rmax may have range from a few degrees < 0.0 to few degrees > 360.0. + // the following function correctly chooses the sky regions on the 0,360 boundary + SkyList *skylist = SkyListByBounds (sky, -1, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + /* load catalog data from region files (hostID is 0 since we are not a client */ + catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); + MARKTIME("-- load catalog data: %f sec\n", dtime); + client_logger_message ("loaded catalog data\n"); + + // NOTE: if I let all hosts load blindly, I saturate the data clients with too many + // relphot_client requests. I need to have the master mediate this. the master + // will not launch the next remote job until this one says it is done + char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); + update_sync_file (syncfile, 1); + + // generate tables go from catID,objID -> catSeq,objSeq + indexCatalogs (catalog, Ncatalog); + client_logger_message ("indexed Catalogs\n"); + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("-- make image bins: %f sec\n", dtime); + + initMosaicBins (catalog, Ncatalog, TRUE); + initMrel (catalog, Ncatalog); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("-- set up image indexes: %f sec\n", dtime); + + findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); + + client_logger_message ("done setting up indexes\n"); + + // dumpObjects ("test.obj.dat", catalog, Ncatalog); + markObjects (catalog, Ncatalog); + + SAVEPLOT = FALSE; + + setExclusions (catalog, Ncatalog, TRUE); + + global_stats (catalog, Ncatalog, flatcorr); + + if (PLOTSTUFF) { + plot_star_coords (catalog, Ncatalog); + // plot_mosaic_fields (catalog); + } + + // XXX : add this in the loop at various points + // if (PLOTSTUFF) plot_scatter (catalog, Ncatalog, flatcorr); + + /* determine fit values */ + client_logger_message ("starting the loops : %d \n", NLOOP); + for (i = 0; i < NLOOP; i++) { + + // set the mean stellar mags given the measurements and the image calibrations + setMrel (catalog, Ncatalog, flatcorr); + + // share mean mags for objects at the boundary (number of unowned meas > 0) + share_mean_mags (catalog, Ncatalog, regionHosts, i); + client_logger_message ("shared mean mag data : loop %d \n", i); + + // load mean mags from other region hosts + slurp_mean_mags (catalog, Ncatalog, regionHosts, i); + client_logger_message ("slurped mean mag data : loop %d \n", i); + + // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements + setMcal (catalog, FALSE, flatcorr); + setMmos (catalog, FALSE, flatcorr); + MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); + + // share image mags for images with non-zero unowned detections + share_image_mags (regionHosts, i); + client_logger_message ("shared image data : loop %d \n", i); + + slurp_image_mags (regionHosts, i); + client_logger_message ("slurped image data : loop %d \n", i); + + 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 % 3 == 2) global_stats (catalog, Ncatalog, flatcorr); + MARKTIME("-- finished loop %d: %f sec\n", i, dtime); + } + client_logger_message ("done with loops\n"); + + // this is a checkpoint to make sure all hosts have finished the loop above + char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loop.sync"); + update_sync_file (loopsyncfile, 0); + free (loopsyncfile); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + char *loopsync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); + check_sync_file (loopsync, 0); + free (loopsync); + } + client_logger_message ("all hosts are done the loops\n"); + + /* set Mcal & Mmos for bad images */ + setMcal (catalog, TRUE, flatcorr); + setMmos (catalog, TRUE, flatcorr); + MARKTIME("-- finalize Mcal values: %f sec\n", dtime); + + setMcalFinal (); // copy per-mosaic calibrations to the images + + share_image_mags (regionHosts, -1); + client_logger_message ("done with parallel images\n"); + + exit (0); +} Index: trunk/Ohana/src/relphot/src/relphot_parallel_regions.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 36630) +++ trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 36630) @@ -0,0 +1,75 @@ +# include "relphot.h" + +int relphot_parallel_regions () { + + int status; + FITS_DB db; + + INITTIME; + + // load the RegionTable (UserRegion should not be used at this level) + RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE); + if (!regionHosts) { + fprintf (stderr, "ERROR: problem with region host table\n"); + exit (2); + } + + // load the flat correction table (if defined) + char flatcorrfile[256]; + sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); + + // register database handle with shutdown procedure + set_db (&db); + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + // lock and load the image db table + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status && UPDATE) { + fprintf (stderr, "error\n"); + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + // if the Image.dat file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + MARKTIME("-- load image data: %f sec\n", dtime); + + /* assign the images to the different region hosts */ + if (!assign_images (&db, regionHosts)) Shutdown ("error assigning images to region hosts"); + MARKTIME("-- assign images: %f sec\n", dtime); + + /* launch processing on the parallel region hosts */ + if (!launch_region_hosts (regionHosts)) Shutdown ("error launching region hosts"); + + // retrieve updated image parameters from the remote hosts (also set Image.mcal) + if (!slurp_image_mags (regionHosts, -1)) Shutdown ("error loading image updates"); + + if (!UPDATE) { + dvo_image_unlock (&db); + MARKTIME ("finished relphot -parallel-regions: %f sec total\n", dtime); + fprintf (stderr, "NOTE: UPDATE is OFF (results are not saved)\n"); + exit (0); + } + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByBounds (sky, -1, regionHosts->Rmin, regionHosts->Rmax, regionHosts->Dmin, regionHosts->Dmax); + + /* update catalogs (in parallel) */ + reload_catalogs (skylist, flatcorr, 0, NULL); + + // save the changes to the image parameters + dvo_image_save (&db, VERBOSE); + dvo_image_unlock (&db); + MARKTIME ("finished relphot -parallel-regions: %f sec total\n", dtime); + + exit (0); +} + Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 36528) +++ trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 36630) @@ -255,5 +255,5 @@ // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter - float Map = PhotAper (&measure[k]); + float Map = PhotAperCat (&measure[k]); if (!isnan(Map)) { Maplist[Nap] = Map - Mcal - Mmos - Mgrid; @@ -263,5 +263,5 @@ } - float Mkron = PhotKron (&measure[k]); + float Mkron = PhotKronCat (&measure[k]); if (!isnan(Mkron)) { Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; @@ -282,4 +282,8 @@ wpsflist[Npsf] = isUbercal ? 10.0 : 1.0; + // NOTE: + // Msys is measure[i].M + zp corrections + // Mcal is image[j].Mcal + // Mmos and Mgrid are offsets for mosaic and grid // tie down reference photometry if the -refcode (code) option is selected Index: trunk/Ohana/src/relphot/src/share_image_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/share_image_mags.c (revision 36630) +++ trunk/Ohana/src/relphot/src/share_image_mags.c (revision 36630) @@ -0,0 +1,120 @@ +# include "relphot.h" + +// we are sharing image calibrations mags for all images which (a) I own and (b) which have unowned detections + +# define D_NIMAGEMAGS 1000 +int share_image_mags (RegionHostTable *regionHosts, int nloop) { + + off_t i, Nimages; + Image *images = getimages (&Nimages, NULL); + + off_t Nimage_mags = 0; + off_t NIMAGE_MAGS = D_NIMAGEMAGS; + + ImageMag *image_mags = NULL; + ALLOCATE (image_mags, ImageMag, NIMAGE_MAGS); + + for (i = 0; i < Nimages; i++) { + // XXX does this image have missing detections (does someone else need it?) + // XXX : NOTE NEED TO FIX THIS: if (imageExtra[i].Nmiss == 0) continue; + + set_image_mags (&image_mags[Nimage_mags], &images[i]); + Nimage_mags ++; + + CHECK_REALLOCATE (image_mags, ImageMag, NIMAGE_MAGS, Nimage_mags, D_NIMAGEMAGS); + } + + // write out the image_mag fits table AND write state in some file + int myHost = regionHosts->index[REGION_HOST_ID]; + char *hostname = regionHosts->hosts[myHost].hostname; + + char *imagfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagemags.fits"); + ImageMagSave (imagfile, image_mags, Nimage_mags); + free (image_mags); + free (imagfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagemags.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_image_mags (RegionHostTable *regionHosts, int nloop) { + + off_t Nimage, i; + Image *images = getimages (&Nimage, NULL); + + int Nimage_mags = 0; + ImageMag *image_mags = NULL; + ALLOCATE (image_mags, ImageMag, 1); + + fprintf (stderr, "grabbing image mags from other hosts...\n"); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset; + char *imagfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.fits"); + ImageMag *image_mags_subset = ImageMagLoad (imagfile, &Nsubset); + free (imagfile); + + image_mags = merge_image_mags (image_mags, &Nimage_mags, image_mags_subset, Nsubset); + } + + for (i = 0; i < Nimage_mags; i++) { + off_t seq = getImageByID (image_mags[i].imageID); + if (seq < 0) { + // XXX is this a problem? (no, other hosts don't know which images I own) + continue; + } + images[seq].Mcal = image_mags[i].Mcal; + images[seq].dMcal = image_mags[i].dMcal; + images[seq].dMagSys = image_mags[i].dMagSys; + images[seq].Xm = image_mags[i].Xm; + images[seq].nFitPhotom = image_mags[i].nFitPhotom; + images[seq].flags = image_mags[i].flags; + images[seq].ubercalDist = image_mags[i].ubercalDist; + } + free (image_mags); + + fprintf (stderr, "DONE grabbing image mags from other hosts\n"); + + return TRUE; +} + +int set_image_mags (ImageMag *image_mags, Image *image) { + + image_mags->Mcal = image->Mcal; + image_mags->dMcal = image->dMcal; + image_mags->dMagSys = image->dMagSys; + image_mags->Xm = image->Xm; + image_mags->nFitPhotom = image->nFitPhotom; + image_mags->flags = image->flags; + image_mags->ubercalDist = image->ubercalDist; + image_mags->imageID = image->imageID; + + return TRUE; +} + +ImageMag *merge_image_mags (ImageMag *target, int *ntarget, ImageMag *source, int Nsource) { + + off_t i; + + REALLOCATE (target, ImageMag, *ntarget + Nsource); + for (i = 0; i < Nsource; i++) { + off_t n = i + *ntarget; + target[n] = source[i]; + } + + free (source); + + *ntarget += Nsource; + return (target); +} + Index: trunk/Ohana/src/relphot/src/share_mean_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 36630) +++ trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 36630) @@ -0,0 +1,149 @@ +# include "relphot.h" +// we are sharing mean mags for all objects which (a) I own and (b) which have unowned detections + +# define D_NMEANMAGS 10000 +int share_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + int i; + off_t j; + + off_t Nmeanmags = 0; + off_t NMEANMAGS = D_NMEANMAGS; + + MeanMag *meanmags = NULL; + ALLOCATE (meanmags, MeanMag, NMEANMAGS); + + int Ns; + int Nsecfilt = GetPhotcodeNsecfilt(); + + int myHost = regionHosts->index[REGION_HOST_ID]; + double Rmin = regionHosts->hosts[myHost].Rmin; + double Rmax = regionHosts->hosts[myHost].Rmax; + double Dmin = regionHosts->hosts[myHost].Dmin; + double Dmax = regionHosts->hosts[myHost].Dmax; + + // XXX skip some catalogs based on UserPatch? + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + // do I own this object? (in region range?) + if (catalog[i].averageT[j].R < Rmin) continue; + if (catalog[i].averageT[j].R >= Rmax) continue; + if (catalog[i].averageT[j].D < Dmin) continue; + if (catalog[i].averageT[j].D >= Dmax) continue; + + // does this object have missing detections (does someone else need it?) + // XXX : sky objects without missing detections + // XXX watch out for detections which are not associated with an image (REF) + if (catalog[i].averageT[j].nOwn == catalog[i].averageT[j].Nmeasure) continue; + + for (Ns = 0; Ns < Nphotcodes; Ns++) { + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + set_mean_mags (&meanmags[Nmeanmags], &catalog[i].averageT[j], &catalog[i].secfilt[Nsecfilt*j + Nsec], Nsec); + Nmeanmags ++; + CHECK_REALLOCATE (meanmags, MeanMag, NMEANMAGS, Nmeanmags, D_NMEANMAGS); + } + } + } + + // write out the meanmag fits table AND write state in some file + char *hostname = regionHosts->hosts[myHost].hostname; + + char *magsfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanmags.fits"); + MeanMagSave (magsfile, meanmags, Nmeanmags); + free (meanmags); + free (magsfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanmags.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + off_t i; + + int Nmeanmags = 0; + MeanMag *meanmags = NULL; + ALLOCATE (meanmags, MeanMag, 1); + + fprintf (stderr, "grabbing mean object mags from other hosts...\n"); + + int Nsecfilt = GetPhotcodeNsecfilt(); + + for (i = 0; i < regionHosts->Nhosts; i++) { + // if (not_neighbor(host[i])) continue; + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset = 0; + char *magsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.fits"); + MeanMag *meanmagsSubset = MeanMagLoad (magsfile, &Nsubset); + free (magsfile); + + // merge_mean_mags reallocs meanmags and frees the input meanmagsSubset + meanmags = merge_mean_mags (meanmags, &Nmeanmags, meanmagsSubset, Nsubset); + } + + for (i = 0; i < Nmeanmags; i++) { + int objID = meanmags[i].objID; + int catID = meanmags[i].catID; + + // set the mean mag + int catSeq; + off_t objSeq; + if (!catID_and_objID_to_seq (catID, objID, &catSeq, &objSeq)) { + // XXX what should I do if this does not match? + continue; + } + + int ecode = GetPhotcodeEquivCodebyCode (meanmags[i].photcode); + if (ecode < 0) continue; + int Nsec = GetPhotcodeNsec (ecode); + if (Nsec < 0) continue; + + catalog[catSeq].secfilt[objSeq*Nsecfilt + Nsec].M = meanmags[i].M; + } + free (meanmags); + + fprintf (stderr, "DONE grabbing mean object mags from other hosts...\n"); + + return TRUE; +} + +int set_mean_mags (MeanMag *meanmags, AverageTiny *average, SecFilt *secfilt, int Nsec) { + + meanmags->M = secfilt->M; + meanmags->dM = secfilt->dM; + meanmags->Xm = secfilt->Xm; + meanmags->Nsec = Nsec; // key to secfilt entry + + meanmags->objID = average->objID; + meanmags->catID = average->catID; + + return TRUE; +} + +MeanMag *merge_mean_mags (MeanMag *target, int *ntarget, MeanMag *source, int Nsource) { + + off_t i; + + REALLOCATE (target, MeanMag, *ntarget + Nsource); + for (i = 0; i < Nsource; i++) { + off_t n = i + *ntarget; + target[n] = source[i]; + } + + free (source); + + *ntarget += Nsource; + return (target); +} + Index: trunk/Ohana/src/relphot/src/syncfile.c =================================================================== --- trunk/Ohana/src/relphot/src/syncfile.c (revision 36630) +++ trunk/Ohana/src/relphot/src/syncfile.c (revision 36630) @@ -0,0 +1,77 @@ +# include "relphot.h" +# define MSG_LENGTH 12 + +char *make_filename (char *dirname, char *hostname, int hostID, char *tailname) { + + char tmp[10], *line; + int Nchar = snprintf (tmp, 0, "%s/%s.%03d.%s", dirname, hostname, hostID, tailname); + + ALLOCATE (line, char, Nchar + 1); + snprintf (line, Nchar + 1, "%s/%s.%03d.%s", dirname, hostname, hostID, tailname); + + return line; +} + +int check_sync_file (char *filename, int nloop) { + + char message[MSG_LENGTH]; + + FILE *f = NULL; + + while (TRUE) { + + f = fopen (filename, "r"); + if (!f) { + usleep (2000000); + continue; + } + + // XXX MSG_LENGTH : 0 EOL byte? + int Nread = fread (message, 1, MSG_LENGTH, f); + if (Nread < MSG_LENGTH) { + fclose (f); + usleep (2000000); + continue; + } + fclose (f); + + // message is of the form: NLOOP: %03d + int loop; + sscanf (message, "%*s %d", &loop); + if (loop != nloop) { + usleep (2000000); + continue; + } + return TRUE; + } + return FALSE; +} + +int clear_sync_file (char *filename) { + // delete file contents + truncate (filename, 0); + + return TRUE; +} + +int update_sync_file (char *filename, int nloop) { + + char message[MSG_LENGTH]; + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "failure to open sync file for write\n"); + exit (4); + } + + snprintf (message, MSG_LENGTH, "NLOOP: %03d\n", nloop); + + int Nwrite = fwrite (message, 1, MSG_LENGTH, f); + if (Nwrite != MSG_LENGTH) { + fprintf (stderr, "failure to write sync message\n"); + exit (3); + } + + fclose (f); + return TRUE; +}