Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 36528) +++ trunk/Ohana/src/relastro/Makefile (revision 36630) @@ -69,4 +69,16 @@ $(SRC)/resort_catalog.$(ARCH).o \ $(SRC)/BrightCatalog.$(ARCH).o \ +$(SRC)/assign_images.$(ARCH).o \ +$(SRC)/launch_region_hosts.$(ARCH).o \ +$(SRC)/relastro_parallel_images.$(ARCH).o \ +$(SRC)/relastro_parallel_regions.$(ARCH).o \ +$(SRC)/MeanPosIO.$(ARCH).o \ +$(SRC)/share_mean_pos.$(ARCH).o \ +$(SRC)/share_images_pos.$(ARCH).o \ +$(SRC)/ImagePosIO.$(ARCH).o \ +$(SRC)/ImageTable.$(ARCH).o \ +$(SRC)/markObjects.$(ARCH).o \ +$(SRC)/indexCatalogs.$(ARCH).o \ +$(SRC)/syncfile.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 36528) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 36630) @@ -6,4 +6,13 @@ # include +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); + // choose off_t or int depending on full-scale relphot analysis resources // # define IDX_T off_t @@ -20,5 +29,7 @@ typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode; -typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS, TARGET_HPM} FitTarget; +typedef enum {OP_NONE, OP_IMAGES, OP_HIGH_SPEED, OP_MERGE_SOURCE, OP_UPDATE_OBJECTS, OP_UPDATE_OFFSETS, OP_LOAD_OBJECTS, OP_HPM, OP_PARALLEL_REGIONS, OP_PARALLEL_IMAGES} RelastroOp; + +typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget; typedef enum { @@ -30,4 +41,20 @@ MARK_BIG_OFFSET = 0x0010, } MeasurementMask; + +typedef struct { + float R; + float D; + unsigned int objID; + unsigned int catID; +} MeanPos; + +typedef struct { + Coords coords; + float dXpixSys; + float dYpixSys; + unsigned int imageID; + int nFitAstrom; + int flags; +} ImagePos; typedef struct { @@ -93,4 +120,5 @@ unsigned int start; unsigned int stop; + off_t myImage; float Mcal; float dMcal; @@ -113,10 +141,4 @@ } StatType; -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } - /* global variables set in parameter file */ # define DVO_MAX_PATH 1024 @@ -129,4 +151,10 @@ char SKY_TABLE[DVO_MAX_PATH]; 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 PARALLEL_REGIONS_MANUAL; int HOST_ID; @@ -208,4 +236,5 @@ FitMode FIT_MODE; +RelastroOp RELASTRO_OP; FitTarget FIT_TARGET; @@ -450,2 +479,44 @@ int hpm_catalogs_parallel (SkyList *skylist); int hpm_objects (SkyRegion *region, Catalog *catalog); + +int strextend (char *input, char *format,...); +int launch_region_hosts (RegionHostTable *regionHosts); + +int assign_images (FITS_DB *db, RegionHostTable *regionHosts); +int select_images_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage); +int calculate_image_bounds (Image *image, double *rmin, double *rmax, double *dmin, double *dmax, double Rmid); +int calculate_host_image_bounds (RegionHostTable *regionHosts); +int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc); + +int relastro_parallel_regions (); +int relastro_parallel_images (); + +char *make_filename (char *dirname, char *hostname, int hostID, char *tailname); +int check_sync_file (char *filename, int nloop); +int clear_sync_file (char *filename); +int update_sync_file (char *filename, int nloop); + +int share_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int slurp_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int set_mean_pos (MeanPos *meanpos, Average *average); +MeanPos *merge_mean_pos (MeanPos *target, int *ntarget, MeanPos *source, int Nsource); + +int MeanPosSave(char *filename, MeanPos *meanpos, off_t Nmeanpos); +MeanPos *MeanPosLoad(char *filename, off_t *nmeanpos); + +int indexCatalogs (Catalog *catalog, int Ncatalog); +int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq); + +int markObjects (Catalog *catalog, int Ncatalog); + +int ImagePosSave(char *filename, ImagePos *image_pos, off_t Nimage_pos); +ImagePos *ImagePosLoad(char *filename, off_t *nimage_pos); + +int share_image_pos (RegionHostTable *regionHosts, int nloop); +int slurp_image_pos (RegionHostTable *regionHosts, int nloop); +ImagePos *merge_image_pos (ImagePos *target, int *ntarget, ImagePos *source, int Nsource); +int set_image_pos (ImagePos *image_pos, Image *image); + +Image *ImageTableLoad(char *filename, off_t *nimage); +int ImageTableSave (char *filename, Image *images, off_t Nimages); +int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage); Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 36528) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 36630) @@ -18,6 +18,7 @@ if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); - GetConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM); // exclude measurements on this basis - GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW); + // set defaults for all of these if they are not used by parallel / remote clients + if (!ScanConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM)) SIGMA_LIM = 0.01; + if (!ScanConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW)) SRC_MEAS_TOOFEW = 3; if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NITER", "%d", 0, &IMFIT_CLIP_NITER)) IMFIT_CLIP_NITER = 3; @@ -25,18 +26,14 @@ if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01; - GetConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN); - GetConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW); - GetConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW); + if (!ScanConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN)) PM_DT_MIN = 0.25; + if (!ScanConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW)) PM_TOOFEW = 4; + if (!ScanConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW)) POS_TOOFEW = 1; + if (!ScanConfig (config, "PAR_FACTOR_MIN", "%lf", 0, &PAR_FACTOR_MIN)) PAR_FACTOR_MIN = 0.2; + if (!ScanConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP)) NX_MAP = 5; + if (!ScanConfig (config, "RELASTRO_MAP_NY", "%d", 0, &NY_MAP)) NY_MAP = 5; + if (!ScanConfig (config, "RELASTRO_DPOS_MAX", "%lf", 0, &DPOS_MAX)) DPOS_MAX = 6.0; + if (!ScanConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS)) ADDSTAR_RADIUS = 1.0; - GetConfig (config, "PAR_FACTOR_MIN", "%lf", 0, &PAR_FACTOR_MIN); - - GetConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP); - GetConfig (config, "RELASTRO_MAP_NY", "%d", 0, &NY_MAP); - GetConfig (config, "RELASTRO_DPOS_MAX", "%lf", 0, &DPOS_MAX); - GetConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS); - - if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS)) { - USE_FIXED_PIXCOORDS = FALSE; - } + if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS)) USE_FIXED_PIXCOORDS = FALSE; // force CATDIR to be absolute (so parallel mode will work) @@ -56,10 +53,6 @@ sprintf (ImageCat, "%s/Images.dat", CATDIR); - if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) { - SKY_DEPTH = 2; - } - if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) { - SKY_TABLE[0] = 0; - } + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) SKY_DEPTH = 2; + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) SKY_TABLE[0] = 0; if (*CATMODE == 0) strcpy (CATMODE, "RAW"); Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 36528) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 36630) @@ -226,6 +226,6 @@ Measure *measureBig = &catalog[cat].measure[meas]; int TESTPT = FALSE; - TESTPT |= (measureBig->imageID == CAT_ID_SRC) && (measureBig->detID == OBJ_ID_SRC); - TESTPT |= (measureBig->imageID == CAT_ID_DST) && (measureBig->detID == OBJ_ID_DST); + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measureBig->imageID == CAT_ID_SRC) && (measureBig->detID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measureBig->imageID == CAT_ID_DST) && (measureBig->detID == OBJ_ID_DST); if (TESTPT) { fprintf (stderr, "got test det\n"); @@ -462,6 +462,6 @@ int TESTPT = FALSE; - TESTPT |= (measure->imageID == CAT_ID_SRC) && (measure->detID == OBJ_ID_SRC); - TESTPT |= (measure->imageID == CAT_ID_DST) && (measure->detID == OBJ_ID_DST); + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measure->imageID == CAT_ID_SRC) && (measure->detID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measure->imageID == CAT_ID_DST) && (measure->detID == OBJ_ID_DST); if (TESTPT) { fprintf (stderr, "got test det\n"); Index: trunk/Ohana/src/relastro/src/ImagePosIO.c =================================================================== --- trunk/Ohana/src/relastro/src/ImagePosIO.c (revision 36630) +++ trunk/Ohana/src/relastro/src/ImagePosIO.c (revision 36630) @@ -0,0 +1,315 @@ +# include "relastro.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) +ImagePos *ImagePosLoad(char *filename, off_t *nimage_pos) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nimage_pos = 0; + ImagePos *image_pos = 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 (crval1, "CRVAL1", double); + GET_COLUMN (crval2, "CRVAL2", double); + GET_COLUMN (crpix1, "CRPIX1", float); + GET_COLUMN (crpix2, "CRPIX2", float); + GET_COLUMN (cdelt1, "CDELT1", float); + GET_COLUMN (cdelt2, "CDELT2", float); + GET_COLUMN (pc1_1, "PC1_1", float); + GET_COLUMN (pc1_2, "PC1_2", float); + GET_COLUMN (pc2_1, "PC2_1", float); + GET_COLUMN (pc2_2, "PC2_2", float); + GET_COLUMN (polyterms, "POLYTERMS", float); // verify that we got 14 columns? + GET_COLUMN (ctype, "CTYPE", char); // verify that we got 15 columns? + GET_COLUMN (Npolyterms, "NPOLYTERMS", char); + GET_COLUMN (dXpixSys, "XPIX_SYS_ERR", float); + GET_COLUMN (dYpixSys, "YPIX_SYS_ERR", float); + GET_COLUMN (imageID, "ID", int); + GET_COLUMN (nFitAstrom, "NFIT", int); + GET_COLUMN (flags, "FLAGS", 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_pos, ImagePos, Nrow); + for (i = 0; i < Nrow; i++) { + image_pos[i].coords.crval1 = crval1 [i]; + image_pos[i].coords.crval2 = crval2 [i]; + image_pos[i].coords.crpix1 = crpix1 [i]; + image_pos[i].coords.crpix2 = crpix2 [i]; + image_pos[i].coords.cdelt1 = cdelt1 [i]; + image_pos[i].coords.cdelt2 = cdelt2 [i]; + image_pos[i].coords.pc1_1 = pc1_1 [i]; + image_pos[i].coords.pc1_2 = pc1_2 [i]; + image_pos[i].coords.pc2_1 = pc2_1 [i]; + image_pos[i].coords.pc2_2 = pc2_2 [i]; + image_pos[i].coords.Npolyterms = Npolyterms[i]; + image_pos[i].dXpixSys = dXpixSys [i]; + image_pos[i].dYpixSys = dYpixSys [i]; + image_pos[i].imageID = imageID [i]; + image_pos[i].nFitAstrom = nFitAstrom[i]; + image_pos[i].flags = flags [i]; + + // polyterms and ctype are a bit different + memcpy (&image_pos[i].coords.polyterms, &polyterms[i*14], 14*sizeof(float)); + memcpy (&image_pos[i].coords.ctype , &ctype [i*15], 15*sizeof(char)); + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (crval1 ); + free (crval2 ); + free (crpix1 ); + free (crpix2 ); + free (cdelt1 ); + free (cdelt2 ); + free (pc1_1 ); + free (pc1_2 ); + free (pc2_1 ); + free (pc2_2 ); + free (polyterms ); + free (ctype ); + free (Npolyterms); + free (dXpixSys ); + free (dYpixSys ); + free (imageID ); + free (nFitAstrom); + free (flags ); + + *nimage_pos = Nrow; + return image_pos; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int ImagePosSave(char *filename, ImagePos *image_pos, off_t Nimage_pos) { + + 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_POS"); + + gfits_define_bintable_column (&theader, "D", "CRVAL1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "CRVAL2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CRPIX1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CRPIX2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CDELT1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CDELT2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC1_1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC1_2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC2_1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC2_2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "14E", "POLYTERMS", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "15A", "CTYPE", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "L", "NPOLYTERMS", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "XPIX_SYS_ERR", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "YPIX_SYS_ERR", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ID", "image ID", "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); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + double *crval1 ; + double *crval2 ; + float *crpix1 ; + float *crpix2 ; + float *cdelt1 ; + float *cdelt2 ; + float *pc1_1 ; + float *pc1_2 ; + float *pc2_1 ; + float *pc2_2 ; + float *polyterms ; + char *ctype ; + char *Npolyterms; + float *dXpixSys ; + float *dYpixSys ; + int *imageID ; + int *nFitAstrom; + int *flags ; + + // create intermediate storage arrays + ALLOCATE (crval1 , double, Nimage_pos); + ALLOCATE (crval2 , double, Nimage_pos); + ALLOCATE (crpix1 , float , Nimage_pos); + ALLOCATE (crpix2 , float , Nimage_pos); + ALLOCATE (cdelt1 , float , Nimage_pos); + ALLOCATE (cdelt2 , float , Nimage_pos); + ALLOCATE (pc1_1 , float , Nimage_pos); + ALLOCATE (pc1_2 , float , Nimage_pos); + ALLOCATE (pc2_1 , float , Nimage_pos); + ALLOCATE (pc2_2 , float , Nimage_pos); + ALLOCATE (polyterms , float , 14*Nimage_pos); + ALLOCATE (ctype , char , 15*Nimage_pos); + ALLOCATE (Npolyterms, char , Nimage_pos); + ALLOCATE (dXpixSys , float , Nimage_pos); + ALLOCATE (dYpixSys , float , Nimage_pos); + ALLOCATE (imageID , int , Nimage_pos); + ALLOCATE (nFitAstrom, int , Nimage_pos); + ALLOCATE (flags , int , Nimage_pos); + + // assign the storage arrays + for (i = 0; i < Nimage_pos; i++) { + crval1 [i] = image_pos[i].coords.crval1 ; + crval2 [i] = image_pos[i].coords.crval2 ; + crpix1 [i] = image_pos[i].coords.crpix1 ; + crpix2 [i] = image_pos[i].coords.crpix2 ; + cdelt1 [i] = image_pos[i].coords.cdelt1 ; + cdelt2 [i] = image_pos[i].coords.cdelt2 ; + pc1_1 [i] = image_pos[i].coords.pc1_1 ; + pc1_2 [i] = image_pos[i].coords.pc1_2 ; + pc2_1 [i] = image_pos[i].coords.pc2_1 ; + pc2_2 [i] = image_pos[i].coords.pc2_2 ; + Npolyterms[i] = image_pos[i].coords.Npolyterms; + dXpixSys [i] = image_pos[i].dXpixSys ; + dYpixSys [i] = image_pos[i].dYpixSys ; + imageID [i] = image_pos[i].imageID ; + nFitAstrom[i] = image_pos[i].nFitAstrom ; + flags [i] = image_pos[i].flags ; + + // polyterms and ctype are a bit different + memcpy (&polyterms[i*14], &image_pos[i].coords.polyterms, 14*sizeof(float)); + memcpy (&ctype [i*15], &image_pos[i].coords.ctype , 15*sizeof(char)); + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "CRVAL1", crval1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRVAL2", crval2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRPIX1", crpix1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRPIX2", crpix2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CDELT1", cdelt1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CDELT2", cdelt2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC1_1", pc1_1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC1_2", pc1_2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC2_1", pc2_1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC2_2", pc2_2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "POLYTERMS", polyterms , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CTYPE", ctype , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "NPOLYTERMS", Npolyterms, Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "XPIX_SYS_ERR",dXpixSys , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "YPIX_SYS_ERR",dYpixSys , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "ID", imageID , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "NFIT", nFitAstrom, Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nimage_pos); + + free (crval1 ); + free (crval2 ); + free (crpix1 ); + free (crpix2 ); + free (cdelt1 ); + free (cdelt2 ); + free (pc1_1 ); + free (pc1_2 ); + free (pc2_1 ); + free (pc2_2 ); + free (polyterms ); + free (ctype ); + free (Npolyterms); + free (dXpixSys ); + free (dYpixSys ); + free (imageID ); + free (nFitAstrom); + free (flags ); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image_pos 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_pos %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for image_pos %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for image_pos %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for image_pos %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_pos %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file image_pos %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing image_pos file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/ImageTable.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageTable.c (revision 36630) +++ trunk/Ohana/src/relastro/src/ImageTable.c (revision 36630) @@ -0,0 +1,62 @@ +# include "relastro.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/relastro/src/MeanPosIO.c =================================================================== --- trunk/Ohana/src/relastro/src/MeanPosIO.c (revision 36630) +++ trunk/Ohana/src/relastro/src/MeanPosIO.c (revision 36630) @@ -0,0 +1,183 @@ +# include "relastro.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) +MeanPos *MeanPosLoad(char *filename, off_t *nmeanpos) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nmeanpos = 0; + MeanPos *meanpos = 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 (R, "RA", double); + GET_COLUMN (D, "DEC", double); + GET_COLUMN (objID, "OBJ_ID", int); + GET_COLUMN (catID, "CAT_ID", int); + + ALLOCATE (meanpos, MeanPos, Nrow); + for (i = 0; i < Nrow; i++) { + meanpos[i].R = R [i]; + meanpos[i].D = D [i]; + meanpos[i].objID = objID[i]; + meanpos[i].catID = catID[i]; + } + fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); + + free (R ); + free (D ); + free (objID); + free (catID); + + // free FITS table pieces... + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + *nmeanpos = Nrow; + return meanpos; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int MeanPosSave(char *filename, MeanPos *meanpos, off_t Nmeanpos) { + + 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", "MEANPOS"); + + gfits_define_bintable_column (&theader, "D", "RA", "mean position, ra", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "mean position, dec", "degrees", 1.0, 0.0); + 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); + + double *R, *D; + unsigned int *objID, *catID; + + // create intermediate storage arrays + ALLOCATE (R, double, Nmeanpos); + ALLOCATE (D, double, Nmeanpos); + ALLOCATE (objID, unsigned int, Nmeanpos); + ALLOCATE (catID, unsigned int, Nmeanpos); + + // assign the storage arrays + for (i = 0; i < Nmeanpos; i++) { + R [i] = meanpos[i].R ; + D [i] = meanpos[i].D ; + objID[i] = meanpos[i].objID; + catID[i] = meanpos[i].catID; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, Nmeanpos); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, Nmeanpos); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeanpos); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeanpos); + + free (R ); + free (D ); + 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 meanpos %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for meanpos %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for meanpos %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 meanpos %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file meanpos %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing meanpos file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relastro/src/MosaicOps.c (revision 36528) +++ trunk/Ohana/src/relastro/src/MosaicOps.c (revision 36630) @@ -11,5 +11,5 @@ // list of mosaic associated with each image -static off_t Nmosaic_for_images; // number of images (for off_ternal checks) +static off_t Nmosaic_for_images; // number of images (for internal checks) static off_t *mosaic_for_images; // array of: image -> mosaic @@ -105,4 +105,5 @@ mosaic[Nmosaic].secz = image[i].secz; mosaic[Nmosaic].coords = image[i].coords; + mosaic[Nmosaic].myImage = i; // init the mosaic_own_images array data @@ -135,4 +136,8 @@ ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images); + // emit an error if we miss mosaics, but stop if we miss too many + int NmissMosaic = 0; + int NtestMosaic = 0; + /* now assign the WRP images to these mosaics */ for (i = 0; i < Nimage; i++) { @@ -140,4 +145,5 @@ if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue; + NtestMosaic ++; /* set image time range */ @@ -147,5 +153,8 @@ Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos); if (Nmos == -1) { - fprintf (stderr, "cannot match mosaic for %s\n", image[i].name); + if (NmissMosaic < 1000) { + fprintf (stderr, "cannot match mosaic for %s\n", image[i].name); + } + NmissMosaic ++; continue; } @@ -164,4 +173,10 @@ } + fprintf (stderr, "mosaic matching : %d of possible %d failed to match\n", NmissMosaic, NtestMosaic); + if (NmissMosaic > 0.5*NtestMosaic) { + fprintf (stderr, "serious problem with mosaic matching\n"); + exit (5); + } + free (startMos); free (stopMos); @@ -250,5 +265,5 @@ if (im >= Nmosaic_for_images) abort(); - // search for the mosaic that + // search for the mosaic that matches this image mos = mosaic_for_images[im]; if (mos < 0) return NULL; @@ -256,2 +271,60 @@ return &mosaic[mos]; } + +// extend each host image table to include the mosaic 'images' needed by the host +int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage) { + + int i; + off_t j; + char *mosaicUsed; + + ALLOCATE (mosaicUsed, char, Nmosaic); + + // we need to add the mosaics to each of the region hosts lists of images + for (i = 0; i < regionHosts->Nhosts; i++) { + + int Nadd = 0; + int NADD = 100; + off_t *addMosaic = NULL; + ALLOCATE (addMosaic, off_t, NADD); + + // reset the mosaicUsed flags (valid only for this host) + memset (mosaicUsed, 0, Nmosaic * sizeof(char)); + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // find the mosaics associated with a given + for (j = 0; j < host->Nimage; j++) { + + int im = host->imseq[j]; + + if (im < 0) abort(); + if (im >= Nmosaic_for_images) abort(); + + // search for the mosaic that matches this image (skip unmatched images) + off_t mos = mosaic_for_images[im]; + if (mos < 0) continue; + + if (mosaicUsed[mos]) continue; + + mosaicUsed[mos] = TRUE; + addMosaic[Nadd] = mos; + Nadd ++; + + CHECK_REALLOCATE (addMosaic, off_t, NADD, Nadd, 100); + } + + REALLOCATE (host->image, Image, host->Nimage + Nadd); + + for (j = 0; j < Nadd; j++) { + off_t mos = addMosaic[j]; + off_t mos_im = mosaic[mos].myImage; + + host->image[host->Nimage + j] = image[mos_im]; + } + + host->Nimage += Nadd; + } + return TRUE; +} + Index: trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/StarMaps.c (revision 36528) +++ trunk/Ohana/src/relastro/src/StarMaps.c (revision 36630) @@ -59,7 +59,6 @@ off_t i, N, Nimages; int xbin, ybin; - struct timeval start, stop; - gettimeofday (&start, (void *) NULL); + INITTIME; // Images *images = getimages(&Nimages, NULL); return value ignored Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 36528) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 36630) @@ -1,5 +1,5 @@ # include "relastro.h" -static off_t Nmax; +static off_t Nmax; static double *X, *dX; static double *Y, *dY; @@ -36,4 +36,22 @@ ALLOCATE (pX, double, MAX (1, Nmax)); ALLOCATE (pY, double, MAX (1, Nmax)); +} + +void freeObjectData () { + + free (R); + free (D); + free (T); + free (X); + free (Y); + + free (dR); + free (dD); + free (dT); + free (dX); + free (dY); + + free (pX); + free (pY); } @@ -205,5 +223,5 @@ if (((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) && (N <= PM_TOOFEW)) mode = FIT_AVERAGE; - if (FIT_TARGET == TARGET_HIGH_SPEED) { + if (RELASTRO_OP == OP_HIGH_SPEED) { Tmean = 0.5*(Tmax - Tmin); } else { @@ -412,4 +430,6 @@ } + freeObjectData (); + if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (NaveSum + NpmSum + NparSum), NaveSum, NpmSum, NparSum, NskipSum, NoffSum); return (TRUE); Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 36528) +++ trunk/Ohana/src/relastro/src/args.c (revision 36630) @@ -13,4 +13,5 @@ /* possible operations */ FIT_TARGET = TARGET_NONE; + RELASTRO_OP = OP_NONE; FIT_MODE = FIT_AVERAGE; @@ -21,5 +22,5 @@ if (N > argc - 6) usage_merge_source(); if (strcmp(argv[N+3], "into")) usage_merge_source(); - FIT_TARGET = TARGET_MERGE_SOURCE; + RELASTRO_OP = OP_MERGE_SOURCE; remove_argument (N, &argc, argv); OBJ_ID_SRC = strtol(argv[N], &endptr, 0); @@ -43,10 +44,51 @@ if ((N = get_argument (argc, argv, "-update-objects"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OBJECTS; + RELASTRO_OP = OP_UPDATE_OBJECTS; } if ((N = get_argument (argc, argv, "-update-offsets"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OFFSETS; + RELASTRO_OP = OP_UPDATE_OFFSETS; + } + + // elements needed for parallel regions / parallel images + 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); + } + + IMAGE_TABLE = NULL; + if ((N = get_argument (argc, argv, "-parallel-images"))) { + remove_argument (N, &argc, argv); + RELASTRO_OP = OP_PARALLEL_IMAGES; + if (N >= argc) usage(); + IMAGE_TABLE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + if (!REGION_FILE) usage(); + } + + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + RELASTRO_OP = OP_IMAGES; + } + + PARALLEL_REGIONS_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-regions"))) { + remove_argument (N, &argc, argv); + RELASTRO_OP = OP_PARALLEL_REGIONS; + if (!REGION_FILE) usage(); + if ((N = get_argument (argc, argv, "-parallel-regions-manual"))) { + remove_argument (N, &argc, argv); + PARALLEL_REGIONS_MANUAL = TRUE; + } } @@ -90,5 +132,5 @@ // XXX include a parallax / no-parallax option if (N >= argc - 4) usage(); - FIT_TARGET = TARGET_HIGH_SPEED; + RELASTRO_OP = OP_HIGH_SPEED; remove_argument (N, &argc, argv); PHOTCODE_A_LIST = strcreate(argv[N]); @@ -104,5 +146,5 @@ if ((N = get_argument (argc, argv, "-hpm"))) { if (N >= argc - 2) usage(); - FIT_TARGET = TARGET_HPM; + RELASTRO_OP = OP_HPM; remove_argument (N, &argc, argv); RADIUS = atof(argv[N]); @@ -116,5 +158,5 @@ remove_argument (N, &argc, argv); PARALLEL_OUTPUT = TRUE; - if ((FIT_TARGET != TARGET_HIGH_SPEED) && (FIT_TARGET != TARGET_HPM)) { + if ((RELASTRO_OP != OP_HIGH_SPEED) && (RELASTRO_OP != OP_HPM)) { fprintf (stderr, "-parallel-output only valid for -high-speed or -hpm modes\n"); exit (1); @@ -143,5 +185,7 @@ } - if (FIT_TARGET == TARGET_NONE) usage(); + if (RELASTRO_OP == OP_NONE) usage(); + + if (((RELASTRO_OP == OP_IMAGES) || (RELASTRO_OP == OP_PARALLEL_REGIONS) || (RELASTRO_OP == OP_PARALLEL_IMAGES)) && (FIT_TARGET == TARGET_NONE)) usage(); /* specify portion of the sky : allow default of all sky? */ @@ -427,5 +471,5 @@ /* possible operations */ - FIT_TARGET = TARGET_NONE; + RELASTRO_OP = TARGET_NONE; FIT_MODE = FIT_AVERAGE; @@ -437,4 +481,9 @@ BCATALOG = NULL; + REGION_FILE = NULL; + REGION_HOST_ID = 0; + IMAGE_TABLE = NULL; + PARALLEL_REGIONS_MANUAL = FALSE; + HOST_ID = 0; if ((N = get_argument (argc, argv, "-hostID"))) { @@ -457,15 +506,15 @@ BCATALOG = strcreate(argv[N]); remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_LOAD_OBJECTS; + RELASTRO_OP = OP_LOAD_OBJECTS; } if ((N = get_argument (argc, argv, "-update-objects"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OBJECTS; + RELASTRO_OP = OP_UPDATE_OBJECTS; } if ((N = get_argument (argc, argv, "-update-offsets"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OFFSETS; + RELASTRO_OP = OP_UPDATE_OFFSETS; } @@ -487,5 +536,5 @@ // XXX include a parallax / no-parallax option if (N >= argc - 5) usage_client(); - FIT_TARGET = TARGET_HIGH_SPEED; + RELASTRO_OP = OP_HIGH_SPEED; remove_argument (N, &argc, argv); PHOTCODE_A_LIST = strcreate(argv[N]); @@ -501,5 +550,5 @@ if ((N = get_argument (argc, argv, "-hpm"))) { if (N >= argc - 3) usage(); - FIT_TARGET = TARGET_HPM; + RELASTRO_OP = OP_HPM; remove_argument (N, &argc, argv); RADIUS = atof(argv[N]); @@ -539,5 +588,5 @@ } - if (FIT_TARGET == TARGET_NONE) usage_client(); + if (RELASTRO_OP == OP_NONE) usage_client(); /* specify portion of the sky : allow default of all sky? */ @@ -718,7 +767,13 @@ void usage () { - fprintf (stderr, "ERROR: USAGE: relastro -update-simple [options]\n"); - fprintf (stderr, " OR: relastro -update-chips [options]\n"); - fprintf (stderr, " OR: relastro -update-mosaic [options]\n"); + fprintf (stderr, "ERROR: USAGE: relastro -images -update-simple [options]\n"); + fprintf (stderr, " OR: relastro -images -update-chips [options]\n"); + fprintf (stderr, " OR: relastro -images -update-mosaic [options]\n"); + fprintf (stderr, " OR: relastro -parallel-regions -update-simple [options]\n"); + fprintf (stderr, " OR: relastro -parallel-regions -update-chips [options]\n"); + fprintf (stderr, " OR: relastro -parallel-regions -update-mosaic [options]\n"); + fprintf (stderr, " OR: relastro -parallel-images -update-simple [options]\n"); + fprintf (stderr, " OR: relastro -parallel-images -update-chips [options]\n"); + fprintf (stderr, " OR: relastro -parallel-images -update-mosaic [options]\n"); fprintf (stderr, " OR: relastro -update-objects [options]\n"); fprintf (stderr, " OR: relastro -high-speed [options]\n"); Index: trunk/Ohana/src/relastro/src/assign_images.c =================================================================== --- trunk/Ohana/src/relastro/src/assign_images.c (revision 36630) +++ trunk/Ohana/src/relastro/src/assign_images.c (revision 36630) @@ -0,0 +1,280 @@ +# include "relastro.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); + + initMosaics (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); + + // supply the mosaics to the image table for the regionHosts : we already have the image + // <-> mosaic relationship, we just need to select these mosaics (once per regionHost) + select_mosaics_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) { + + 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); + ALLOCATE (regionHosts->hosts[i].imseq, off_t, regionHosts->hosts[i].NIMAGE); + } + + for (j = 0; j < Nimage; j++) { + + /* select images by photcode, or equiv photcode, if specified */ + if (NphotcodesKeep > 0) { + int found = FALSE; + // XXX this bit of code excludes DIS mosaics and should be fixed + for (i = 0; (i < NphotcodesKeep) && !found; i++) { + if (photcodesKeep[i][0].code == image[j].photcode) found = TRUE; + if (photcodesKeep[i][0].code == GetPhotcodeEquivCodebyCode(image[j].photcode)) found = TRUE; + } + if (!found) continue; + } + if (NphotcodesSkip > 0) { + int found = FALSE; + for (i = 0; (i < NphotcodesSkip) && !found; i++) { + if (photcodesSkip[i][0].code == image[j].photcode) found = TRUE; + if (photcodesSkip[i][0].code == GetPhotcodeEquivCodebyCode(image[j].photcode)) 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 + 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; + + 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 image center thus, an image + // with center 0.5 will have chips bounds ranging from ~ -1.5 to +2.5 or so, while an + // image 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].imseq[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].imseq, 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) { + + 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 = 360.0; + double Rmax = 0.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, 180.0); + + Rmin = MIN (Rmin, Rc); + Dmin = MIN (Dmin, Rc); + + Rmax = MAX (Rmax, Rc); + Dmax = MAX (Dmax, Rc); + } + + } + + 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; +} + + +// XXX add a search tree to speed this up? +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/relastro/src/indexCatalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/indexCatalogs.c (revision 36630) +++ trunk/Ohana/src/relastro/src/indexCatalogs.c (revision 36630) @@ -0,0 +1,66 @@ +# include "relastro.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; + + if (!Ncatalog) return TRUE; + + // 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].average[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].average[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/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 36528) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 36630) @@ -6,5 +6,5 @@ args (argc, argv); - if (FIT_TARGET == TARGET_MERGE_SOURCE) return; + if (RELASTRO_OP == OP_MERGE_SOURCE) return; fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); Index: trunk/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 36630) +++ trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 36630) @@ -0,0 +1,144 @@ +# include "relastro.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, "relastro -parallel-images %s", filename); + strextend (command, "-region-hosts %s", REGION_FILE); + strextend (command, "-region-hostID %d", host->hostID); + strextend (command, "-D CATDIR %s", CATDIR); + strextend (command, "-region %f %f %f %f", host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat); + strextend (command, "-statmode %s", STATMODE); + strextend (command, "-minerror %f", MIN_ERROR); + strextend (command, "-nloop %d", NLOOP); + strextend (command, "-threads %d", NTHREADS); + + switch (FIT_TARGET) { + case TARGET_SIMPLE: + strextend (command, "-update-simple"); + break; + case TARGET_CHIPS: + strextend (command, "-update-chips"); + break; + case TARGET_MOSAICS: + strextend (command, "-update-mosaics"); + break; + case TARGET_NONE: + abort(); + } + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (UPDATE) strextend (command, "-update"); + if (PARALLEL) strextend (command, "-parallel"); + if (PARALLEL_MANUAL) strextend (command, "-parallel-manual"); + if (PARALLEL_SERIAL) strextend (command, "-parallel-serial"); + if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST); + + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (PhotFlagSelect) strextend (command, "+photflags"); + if (PhotFlagBad) strextend (command, "+photflagbad %d", PhotFlagBad); + if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); + + if (MinBadQF > 0.0) strextend (command, "-min-bad-psfqf %f", MinBadQF); + if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset %f", MaxMeanOffset); + + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + strextend (command, "-time %s %s", tstart, tstop); + free (tstart); + free (tstop); + } + + 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 + + // remove client is done, go ahead with next client + check_sync_file (syncfile, 1); + free (syncfile); + } + + if (PARALLEL_REGIONS_MANUAL) { + fprintf (stderr, "run the relastro_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/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 36528) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 36630) @@ -115,4 +115,9 @@ 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); + // load the list of hosts HostTable *table = HostTableLoad (CATDIR, sky->hosts); @@ -131,5 +136,5 @@ ALLOCATE (table->hosts[i].results, char, 1024); - snprintf (table->hosts[i].results, 1024, "%s/relastro.catalog.subset.dat", table->hosts[i].pathname); + snprintf (table->hosts[i].results, 1024, "%s/relastro.catalog.%s.dat", table->hosts[i].pathname, uniquer); // options / arguments that can affect relastro_client -load: Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 36528) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 36630) @@ -1,9 +1,3 @@ # include "relastro.h" - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap) { @@ -12,7 +6,6 @@ off_t Nimage, Nsubset; off_t *LineNumber; - struct timeval start, stop; - gettimeofday (&start, (void *) NULL); + INITTIME; // convert database table to internal structure @@ -35,5 +28,5 @@ /* unlock, if we can (else, unlocked below) */ - int unlockImages = !UPDATE || (FIT_TARGET == TARGET_UPDATE_OFFSETS); + int unlockImages = !UPDATE || (RELASTRO_OP == OP_UPDATE_OFFSETS); if (unlockImages) dvo_image_unlock (db); Index: trunk/Ohana/src/relastro/src/markObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/markObjects.c (revision 36630) +++ trunk/Ohana/src/relastro/src/markObjects.c (revision 36630) @@ -0,0 +1,24 @@ +# include "relastro.h" + +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++) { + ALLOCATE (catalog[i].nOwn, int, catalog[i].Naverage); + memset (catalog[i].nOwn, 0, catalog[i].Naverage*sizeof(int)); + for (j = 0; j < catalog[i].Naverage; j++) { + int nOwn = 0; + int m = catalog[i].average[j].measureOffset; + for (n = 0; n < catalog[i].average[j].Nmeasure; n++) { + if (!catalog[i].measureT[m+n].myDet) continue; + nOwn ++; + } + catalog[i].nOwn[j] = nOwn; + } + } + return TRUE; +} + Index: trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro.c (revision 36528) +++ trunk/Ohana/src/relastro/src/relastro.c (revision 36630) @@ -10,33 +10,41 @@ SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - switch (FIT_TARGET) { - case TARGET_UPDATE_OBJECTS: + switch (RELASTRO_OP) { + case OP_UPDATE_OBJECTS: /* the object analysis is a separate process iterating over catalogs */ relastro_objects (skylist, 0, NULL); exit (0); - case TARGET_HIGH_SPEED: + case OP_HIGH_SPEED: /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */ high_speed_catalogs (sky, skylist, 0, NULL); exit (0); - case TARGET_HPM: + case OP_HPM: hpm_catalogs (sky, skylist, 0, NULL); exit (0); - case TARGET_MERGE_SOURCE: + case OP_MERGE_SOURCE: /* a special method to manually merge unlinked detections of sources togther (not parallel) */ relastro_merge_source (sky); exit (0); - case TARGET_SIMPLE: - case TARGET_CHIPS: - case TARGET_MOSAICS: + case OP_IMAGES: relastro_images (skylist); exit (0); - case TARGET_UPDATE_OFFSETS: + case OP_UPDATE_OFFSETS: // iterate over catalogs to make detection coordinates consistant UpdateObjectOffsets (skylist, 0, NULL); + exit (0); + + case OP_PARALLEL_REGIONS: + // run image updates in parallel across multiple remote machines + relastro_parallel_regions (); + exit (0); + + case OP_PARALLEL_IMAGES: + // operation on the remote machines in the PARALLEL_REGION mode + relastro_parallel_images (); exit (0); Index: trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_client.c (revision 36528) +++ trunk/Ohana/src/relastro/src/relastro_client.c (revision 36630) @@ -28,23 +28,23 @@ SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - switch (FIT_TARGET) { + switch (RELASTRO_OP) { - case TARGET_LOAD_OBJECTS: { + case OP_LOAD_OBJECTS: { // USAGE: relastro_client -load-objects int Ncatalog; Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR); if (!catalog) { - fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); - exit (2); + fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); + exit (2); } BrightCatalog *bcatalog = BrightCatalogMerge (catalog, Ncatalog); if (!BrightCatalogSave (BCATALOG, bcatalog)) { - fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); - exit (2); + fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); + exit (2); } break; } - case TARGET_UPDATE_OBJECTS: { + case OP_UPDATE_OBJECTS: { // USAGE: relastro_client -update-objects relastro_objects (skylist, HOST_ID, HOSTDIR); @@ -52,5 +52,5 @@ } - case TARGET_HIGH_SPEED: { + case OP_HIGH_SPEED: { // USAGE: relastro_client -high-speed high_speed_catalogs (sky, skylist, HOST_ID, HOSTDIR); @@ -58,5 +58,5 @@ } - case TARGET_HPM: { + case OP_HPM: { // USAGE: relastro_client -high-speed hpm_catalogs (sky, skylist, HOST_ID, HOSTDIR); @@ -64,6 +64,6 @@ } - // XXX loading the images is fairly costly -- see if we can do an image subset - case TARGET_UPDATE_OFFSETS: { + // XXX loading the images is fairly costly -- see if we can do an image subset? + case OP_UPDATE_OFFSETS: { FITS_DB db; Index: trunk/Ohana/src/relastro/src/relastro_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_images.c (revision 36528) +++ trunk/Ohana/src/relastro/src/relastro_images.c (revision 36630) @@ -6,7 +6,6 @@ Catalog *catalog; FITS_DB db; - struct timeval start, stop; - gettimeofday (&start, (void *) NULL); + INITTIME; /* register database handle with shutdown procedure */ Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 36528) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 36630) @@ -144,5 +144,4 @@ // PM_TOOFEW // SRC_MEAS_TOOFEW - // FIT_TARGET char command[1024]; Index: trunk/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 36630) +++ trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 36630) @@ -0,0 +1,131 @@ +# include "relastro.h" + +/* This function is essentially identical to relastro_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 relastro_parallel_images () { + + int i, Ncatalog; + Catalog *catalog = NULL; + + INITTIME; + + // 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); + + // 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... + initMosaics (image, Nimage); + + initImages (image, NULL, Nimage); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByBounds (sky, -1, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + initCoords(); + + /* load catalog data from region files (hostID is 0 since we are not a client */ + catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL); + MARKTIME("-- load catalog data: %f sec\n", dtime); + + // NOTE: if I let all hosts load blindly, I saturate the data clients with too many + // relastro_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); + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("-- make image bins: %f sec\n", dtime); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("-- set up image indexes: %f sec\n", dtime); + + // set test points based on the starmap + createStarMap (catalog, Ncatalog); + + markObjects (catalog, Ncatalog); + + SAVEPLOT = FALSE; + + /* major modes */ + switch (FIT_TARGET) { + case TARGET_SIMPLE: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + share_mean_pos (catalog, Ncatalog, regionHosts, i); + slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + UpdateSimple (catalog, Ncatalog); + share_image_pos (regionHosts, i); + slurp_image_pos (regionHosts, i); + } + break; + + case TARGET_CHIPS: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + share_mean_pos (catalog, Ncatalog, regionHosts, i); + slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + UpdateChips (catalog, Ncatalog); + share_image_pos (regionHosts, i); + slurp_image_pos (regionHosts, i); + MARKTIME("update chips: %f sec\n", dtime); + } + // create summary plots of the process + // relastroVisualSummaryChips(); + break; + + case TARGET_MOSAICS: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + share_mean_pos (catalog, Ncatalog, regionHosts, i); + slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + UpdateMosaic (catalog, Ncatalog); + share_image_pos (regionHosts, i); + slurp_image_pos (regionHosts, i); + } + break; + + default: + fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__); + exit (2); + } + + // 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); + } + + share_image_pos (regionHosts, -1); + + exit (0); +} Index: trunk/Ohana/src/relastro/src/relastro_parallel_regions.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 36630) +++ trunk/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 36630) @@ -0,0 +1,79 @@ +# include "relastro.h" + +int relastro_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); + } + + // 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_pos (regionHosts, -1)) Shutdown ("error loading image updates"); + + if (!UPDATE) { + dvo_image_unlock (&db); + MARKTIME ("finished relastro -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); + + if (PARALLEL) { + // save the updated image parameters + dvo_image_save (&db, VERBOSE); + dvo_image_unlock (&db); + } + + // iterate over catalogs to make detection coordinates consistant + UpdateObjectOffsets (skylist, 0, NULL); + + if (!PARALLEL) { + // save the changes to the image parameters + dvo_image_save (&db, VERBOSE); + dvo_image_unlock (&db); + } + + MARKTIME ("finished relastro -parallel-regions: %f sec total\n", dtime); + + exit (0); +} + Index: trunk/Ohana/src/relastro/src/resort_catalog.c =================================================================== --- trunk/Ohana/src/relastro/src/resort_catalog.c (revision 36528) +++ trunk/Ohana/src/relastro/src/resort_catalog.c (revision 36630) @@ -30,6 +30,5 @@ if (catalog[0].sorted == TRUE) return; - // struct timeval start, stop; - // gettimeofday (&start, NULL); + // INITTIME; /* internal counters */ Index: trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- trunk/Ohana/src/relastro/src/select_images.c (revision 36528) +++ trunk/Ohana/src/relastro/src/select_images.c (revision 36630) @@ -15,10 +15,4 @@ void dsortindex (double *X, off_t *Y, int N); off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap) { @@ -31,5 +25,4 @@ Coords tcoords; SkyRegionCoords *skycoords; - struct timeval start, stop; double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion; @@ -50,5 +43,5 @@ } - gettimeofday (&start, (void *) NULL); + INITTIME; // the comparison is made in the catalog local projection. below we set crval1,2 @@ -132,5 +125,6 @@ if (NphotcodesKeep > 0) { found = FALSE; - // XXX this bit of code excludes DIS mosaics and should be fixed + // we have to keep DIS mosaics explicitly (photcode = 0) + if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) found = TRUE; for (k = 0; (k < NphotcodesKeep) && !found; k++) { if (photcodesKeep[k][0].code == timage[i].photcode) found = TRUE; Index: trunk/Ohana/src/relastro/src/share_images_pos.c =================================================================== --- trunk/Ohana/src/relastro/src/share_images_pos.c (revision 36630) +++ trunk/Ohana/src/relastro/src/share_images_pos.c (revision 36630) @@ -0,0 +1,117 @@ +# include "relastro.h" + +// we are sharing image astrometry calibrations for all images which (a) I own and (b) which have unowned detections + +# define D_NIMAGEPOS 1000 +int share_image_pos (RegionHostTable *regionHosts, int nloop) { + + off_t i, Nimages; + Image *images = getimages (&Nimages, NULL); + + off_t Nimage_pos = 0; + off_t NIMAGE_POS = D_NIMAGEPOS; + + ImagePos *image_pos = NULL; + ALLOCATE (image_pos, ImagePos, NIMAGE_POS); + + 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_pos (&image_pos[Nimage_pos], &images[i]); + Nimage_pos ++; + + CHECK_REALLOCATE (image_pos, ImagePos, NIMAGE_POS, Nimage_pos, D_NIMAGEPOS); + } + + // 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 *iposfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagepos.fits"); + ImagePosSave (iposfile, image_pos, Nimage_pos); + free (image_pos); + free (iposfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagepos.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_image_pos (RegionHostTable *regionHosts, int nloop) { + + off_t Nimage, i; + Image *images = getimages (&Nimage, NULL); + + int Nimage_pos = 0; + ImagePos *image_pos = NULL; + ALLOCATE (image_pos, ImagePos, 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, "imagepos.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset; + char *iposfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagepos.fits"); + ImagePos *image_pos_subset = ImagePosLoad (iposfile, &Nsubset); + free (iposfile); + + image_pos = merge_image_pos (image_pos, &Nimage_pos, image_pos_subset, Nsubset); + } + + for (i = 0; i < Nimage_pos; i++) { + off_t seq = getImageByID (image_pos[i].imageID); + if (seq < 0) { + // XXX is this a problem? (no, other hosts don't know which images I own) + continue; + } + images[seq].coords = image_pos[i].coords ; + images[seq].dXpixSys = image_pos[i].dXpixSys ; + images[seq].dYpixSys = image_pos[i].dYpixSys ; + images[seq].imageID = image_pos[i].imageID ; + images[seq].nFitAstrom = image_pos[i].nFitAstrom; + images[seq].flags = image_pos[i].flags ; + } + free (image_pos); + + fprintf (stderr, "DONE grabbing image mags from other hosts\n"); + + return TRUE; +} + +int set_image_pos (ImagePos *image_pos, Image *image) { + + image_pos->coords = image->coords; + image_pos->dXpixSys = image->dXpixSys; + image_pos->dYpixSys = image->dYpixSys; + image_pos->imageID = image->imageID; + image_pos->nFitAstrom = image->nFitAstrom; + image_pos->flags = image->flags; + + return TRUE; +} + +ImagePos *merge_image_pos (ImagePos *target, int *ntarget, ImagePos *source, int Nsource) { + + off_t i; + + REALLOCATE (target, ImagePos, *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/relastro/src/share_mean_pos.c =================================================================== --- trunk/Ohana/src/relastro/src/share_mean_pos.c (revision 36630) +++ trunk/Ohana/src/relastro/src/share_mean_pos.c (revision 36630) @@ -0,0 +1,132 @@ +# include "relastro.h" +// we are sharing mean positions for all objects which (a) I own and (b) which have unowned detections + +# define D_NMEANPOS 10000 +int share_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + int i; + off_t j; + + off_t Nmeanpos = 0; + off_t NMEANPOS = D_NMEANPOS; + + MeanPos *meanpos = NULL; + ALLOCATE (meanpos, MeanPos, NMEANPOS); + + 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?) --- CAREFUL HERE!! + if (catalog[i].average[j].R < Rmin) continue; + if (catalog[i].average[j].R >= Rmax) continue; + if (catalog[i].average[j].D < Dmin) continue; + if (catalog[i].average[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].nOwn[j] == catalog[i].average[j].Nmeasure) continue; + + set_mean_pos (&meanpos[Nmeanpos], &catalog[i].average[j]); + Nmeanpos ++; + CHECK_REALLOCATE (meanpos, MeanPos, NMEANPOS, Nmeanpos, D_NMEANPOS); + } + } + + // write out the meanmag fits table AND write state in some file + char *hostname = regionHosts->hosts[myHost].hostname; + + char *posfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanpos.fits"); + MeanPosSave (posfile, meanpos, Nmeanpos); + free (meanpos); + free (posfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanpos.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + off_t i; + + int Nmeanpos = 0; + MeanPos *meanpos = NULL; + ALLOCATE (meanpos, MeanPos, 1); + + fprintf (stderr, "grabbing mean object pos 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, "meanpos.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset = 0; + char *posfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanpos.fits"); + MeanPos *meanposSubset = MeanPosLoad (posfile, &Nsubset); + free (posfile); + + // merge_mean_pos reallocs meanpos and frees the input meanposSubset + meanpos = merge_mean_pos (meanpos, &Nmeanpos, meanposSubset, Nsubset); + } + + for (i = 0; i < Nmeanpos; i++) { + int objID = meanpos[i].objID; + int catID = meanpos[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; + } + + catalog[catSeq].average[objSeq].R = meanpos[i].R; + catalog[catSeq].average[objSeq].D = meanpos[i].D; + } + free (meanpos); + + fprintf (stderr, "DONE grabbing mean object pos from other hosts...\n"); + + return TRUE; +} + +int set_mean_pos (MeanPos *meanpos, Average *average) { + + meanpos->R = average->R; + meanpos->D = average->D; + meanpos->objID = average->objID; + meanpos->catID = average->catID; + + return TRUE; +} + +MeanPos *merge_mean_pos (MeanPos *target, int *ntarget, MeanPos *source, int Nsource) { + + off_t i; + + REALLOCATE (target, MeanPos, *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/relastro/src/syncfile.c =================================================================== --- trunk/Ohana/src/relastro/src/syncfile.c (revision 36630) +++ trunk/Ohana/src/relastro/src/syncfile.c (revision 36630) @@ -0,0 +1,77 @@ +# include "relastro.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; +}