Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/Makefile =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/Makefile (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/Makefile (revision 35593) @@ -1,3 +1,3 @@ -default: uniphot setphot setphot_client setfwhm +default: uniphot setphot setphot_client setfwhm setposangle setposangle_client setastrom setastrom_client fiximids fiximids_client help: @echo "make options: uniphot setphot setfwhm default help install default (uniphot setphot setfwhm)" @@ -25,6 +25,8 @@ setastrom: $(BIN)/setastrom.$(ARCH) setastrom_client: $(BIN)/setastrom_client.$(ARCH) +fiximids: $(BIN)/fiximids.$(ARCH) +fiximids_client: $(BIN)/fiximids_client.$(ARCH) -install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setphot $(DESTBIN)/setphot_client $(DESTBIN)/setposangle $(DESTBIN)/setposangle_client $(DESTBIN)/setastrom $(DESTBIN)/setastrom_client +install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setphot $(DESTBIN)/setphot_client $(DESTBIN)/setposangle $(DESTBIN)/setposangle_client $(DESTBIN)/setastrom $(DESTBIN)/setastrom_client $(DESTBIN)/fiximids $(DESTBIN)/fiximids_client UNIPHOT = \ @@ -137,2 +139,29 @@ $(SETFWHM): $(INC)/uniphot.h $(BIN)/setfwhm.$(ARCH): $(SETFWHM) + +FIXIMIDS = \ +$(SRC)/fiximids.$(ARCH).o \ +$(SRC)/initialize_fiximids.$(ARCH).o \ +$(SRC)/load_images_fiximids.$(ARCH).o \ +$(SRC)/ImageSubsetFixImageIDs.$(ARCH).o \ +$(SRC)/update_dvo_fiximids.$(ARCH).o \ +$(SRC)/update_catalog_fiximids.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o + +$(FIXIMIDS): $(INC)/fiximids.h +$(BIN)/fiximids.$(ARCH): $(FIXIMIDS) + +FIXIMIDS_CLIENT = \ +$(SRC)/fiximids_client.$(ARCH).o \ +$(SRC)/initialize_fiximids.$(ARCH).o \ +$(SRC)/load_images_fiximids.$(ARCH).o \ +$(SRC)/ImageSubsetFixImageIDs.$(ARCH).o \ +$(SRC)/update_dvo_fiximids.$(ARCH).o \ +$(SRC)/update_catalog_fiximids.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o + +$(FIXIMIDS_CLIENT): $(INC)/fiximids.h +$(BIN)/fiximids_client.$(ARCH): $(FIXIMIDS_CLIENT) + Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/include/fiximids.h =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/include/fiximids.h (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/include/fiximids.h (revision 35593) @@ -2,10 +2,4 @@ # include # include - -typedef struct { - unsigned int start; - unsigned int stop; - Coords coords; -} Mosaic; // to determine the image ID, we need time (tzero + trate) and photcode @@ -14,4 +8,5 @@ unsigned short NX; // image width unsigned short NY; // image height + unsigned int nstar; // number of stars on the image e_time tzero; // readout time (row 0) unsigned char trate; // scan rate (100 usec/pixel) @@ -28,5 +23,7 @@ char *SINGLE_CPT; int VERBOSE; -int RESET; +int VERBOSE_IMSTATS; +int SUMMARY_IMSTATS; +char *IMSTATS_FILE; // location of partial imstats file for remote client int UPDATE; int PARALLEL; @@ -40,22 +37,16 @@ void ConfigInit PROTO((int *argc, char **argv)); -void usage_fiximids PROTO(()); -void initialize_fiximids PROTO((int argc, char **argv)); -int args_fiximids PROTO((int argc, char **argv)); - -void usage_fiximids_client PROTO(()); -void initialize_fiximids_client PROTO((int argc, char **argv)); -int args_fiximids_client PROTO((int argc, char **argv)); - -// int update_dvo_fiximids_client PROTO((ImageSubset *image, off_t Nimage, FlatCorrectionTable *flatcorr)); +void usage_fiximids PROTO(()); +void initialize_fiximids PROTO((int argc, char **argv)); +int args_fiximids PROTO((int argc, char **argv)); + +void usage_fiximids_client PROTO(()); +void initialize_fiximids_client PROTO((int argc, char **argv)); +int args_fiximids_client PROTO((int argc, char **argv)); ImageSubset *ImageSubsetLoad PROTO((char *filename, off_t *nimage)); int ImageSubsetSave PROTO((char *filename, ImageSubset *image, off_t Nimage)); -Image *ImagesFromSubset PROTO((ImageSubset *subset, off_t N)); ImageSubset *ImagesToSubset PROTO((Image *image, off_t N)); -void lock_image_db PROTO((FITS_DB *db, char *filename)); -void unlock_image_db PROTO((FITS_DB *db)); -void create_image_db PROTO((FITS_DB *db)); void set_db PROTO((FITS_DB *in)); int Shutdown PROTO((char *format, ...) OHANA_FORMAT(printf, 1, 2) ); @@ -64,13 +55,19 @@ int SetSignals PROTO((void)); -Image *load_images_fiximids PROTO((FITS_DB *db, off_t *Nimage)); -int update_dvo_fiximids PROTO((ImageSubset *image, off_t Nimage)); -int update_dvo_fiximids_parallel PROTO((SkyTable *sky, ImageSubset *image, off_t Nimage)); +Image *load_images_fiximids PROTO((FITS_DB *db, off_t *Nimage)); +int update_dvo_fiximids PROTO((ImageSubset *image, off_t Nimage)); +int update_dvo_fiximids_parallel PROTO((SkyTable *sky, ImageSubset *image, off_t Nimage)); + +void update_catalog_fiximids PROTO((Catalog *catalog)); -void update_catalog_fiximids PROTO((Catalog *catalog, ImageSubset *image, off_t *index, off_t Nimage)); -int fiximids_local_astrometry PROTO((float *posAngle, float *pltScale, double x, double y, Coords *mosaic, Coords *coords)); +void sort_image_times PROTO((e_time *S, e_time *E, short *C, unsigned int *I, unsigned int *Q, off_t N)); +int FindImageID PROTO((off_t *ID, off_t *Seq, e_time time, short photcode)); +void initImageIndex PROTO((ImageSubset *image, off_t Nimage_init)); -off_t getMosaicByTimes PROTO((unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos)); -void sort_mosaic_times PROTO((unsigned int *S, unsigned int *E, off_t *I, off_t N)); -void initMosaics PROTO((ImageSubset *image, off_t Nimage)); -Mosaic *getMosaicForImage PROTO((off_t im)); +void BumpValidImage PROTO((int Seq)); +void BumpInvalidImage PROTO((int Seq)); +void CompareImageCounts PROTO((ImageSubset *image, off_t Nimage_comp)); +void SummaryImageStats PROTO((ImageSubset *image, off_t Nimage_comp)); + +int ImageValidSave PROTO((char *filename)); +int ImageValidLoad PROTO((char *filename)); Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/ImageSubsetFixImageIDs.c =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/ImageSubsetFixImageIDs.c (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/ImageSubsetFixImageIDs.c (revision 35593) @@ -1,3 +1,5 @@ # include "fiximids.h" +# define FT_BZERO_INT16 1.0*0x8000 +# define FT_BZERO_INT32 1.0*0x80000000 # define GET_COLUMN(OUT,NAME,TYPE) \ @@ -8,5 +10,5 @@ int Ncol; - off_t i, j, k; + off_t i; off_t Nrow; Header header; @@ -55,8 +57,9 @@ GET_COLUMN (imageID, "IMAGE_ID", int); - GET_COLUMN (NX, "NX", int); - GET_COLUMN (NY, "NY", int); - GET_COLUMN (tzero, "TZERO", float); - GET_COLUMN (trate, "TRATE", float); + GET_COLUMN (NX, "NX", short); + GET_COLUMN (NY, "NY", short); + GET_COLUMN (nstar, "nstar", int); + GET_COLUMN (tzero, "TZERO", int); + GET_COLUMN (trate, "TRATE", short); GET_COLUMN (photcode, "PHOTCODE" , short); @@ -66,4 +69,5 @@ image[i].NX = NX[i]; image[i].NY = NY[i]; + image[i].nstar = nstar[i]; image[i].tzero = tzero[i]; image[i].trate = trate[i]; @@ -75,4 +79,5 @@ free (NX); free (NY); + free (nstar); free (tzero); free (trate); @@ -85,5 +90,5 @@ int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage) { - off_t i, j, k; + off_t i; Header header; Header theader; @@ -100,8 +105,9 @@ // an unsigned int needs to have bzero of 0x8000 gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", "tmp", 1.0, FT_BZERO_INT32); - gfits_define_bintable_column (&theader, "J", "NX", "tmp", "tmp", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NY", "tmp", "tmp", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "TZERO", "tmp", "tmp", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "TRATE", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NX", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NY", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NSTAR", "tmp", "tmp", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "TZERO", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "TRATE", "tmp", "tmp", 1.0, 0.0); gfits_define_bintable_column (&theader, "I", "PHOTCODE", "tmp", "tmp", 1.0, 0.0); @@ -109,16 +115,18 @@ gfits_create_table (&theader, &ftable); - unsigned int *imageID; - int *NX, *NY; - float *tzero, *trate; + unsigned int *imageID, *nstar; + short *NX, *NY; + e_time *tzero; + short *trate; short *photcode; // create intermediate storage arrays ALLOCATE (imageID, unsigned int, Nimage); - ALLOCATE (NX, int, Nimage); - ALLOCATE (NY, int, Nimage); - ALLOCATE (tzero, float, Nimage); - ALLOCATE (trate, float, Nimage); - ALLOCATE (photcode, int, Nimage); + ALLOCATE (NX, short, Nimage); + ALLOCATE (NY, short, Nimage); + ALLOCATE (nstar, unsigned int, Nimage); + ALLOCATE (tzero, e_time, Nimage); + ALLOCATE (trate, short, Nimage); + ALLOCATE (photcode, short, Nimage); // assign the storage arrays @@ -127,4 +135,5 @@ NX[i] = image[i].NX ; NY[i] = image[i].NY ; + nstar[i] = image[i].nstar ; tzero[i] = image[i].tzero ; trate[i] = image[i].trate ; @@ -136,4 +145,5 @@ gfits_set_bintable_column (&theader, &ftable, "NX", NX, Nimage); gfits_set_bintable_column (&theader, &ftable, "NY", NY, Nimage); + gfits_set_bintable_column (&theader, &ftable, "NSTAR", nstar, Nimage); gfits_set_bintable_column (&theader, &ftable, "TZERO", tzero, Nimage); gfits_set_bintable_column (&theader, &ftable, "TRATE", trate, Nimage); @@ -143,4 +153,5 @@ free (NX ); free (NY ); + free (nstar ); free (tzero ); free (trate ); @@ -169,26 +180,4 @@ } -// XXX drop these two? - -# if (0) -Image *ImagesFromSubset (ImageSubset *subset, off_t N) { - - off_t i; - - // we have been given an ImageSubset array, containing a reduced set of image fields - // create full a Image array and save the needed values - Image *image = NULL; - ALLOCATE (image, Image, N); - for (i = 0; i < N; i++) { - image[i].imageID = subset[i].imageID; - image[i].NX = subset[i].NX ; - image[i].NY = subset[i].NY ; - image[i].tzero = subset[i].tzero ; - image[i].trate = subset[i].trate ; - image[i].coords = subset[i].coords ; - } - return image; -} - ImageSubset *ImagesToSubset (Image *image, off_t N) { @@ -201,13 +190,14 @@ for (i = 0; i < N; i++) { - subset[i].imageID = image[i].imageID; - subset[i].NX = image[i].NX ; - subset[i].NY = image[i].NY ; - subset[i].tzero = image[i].tzero ; - subset[i].trate = image[i].trate ; - subset[i].coords = image[i].coords ; + subset[i].imageID = image[i].imageID ; + subset[i].NX = image[i].NX ; + subset[i].NY = image[i].NY ; + subset[i].nstar = image[i].nstar ; + subset[i].tzero = image[i].tzero ; + subset[i].trate = image[i].trate ; + subset[i].photcode = image[i].photcode ; } return subset; } -# endif (0) + Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids.c =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids.c (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids.c (revision 35593) @@ -24,4 +24,12 @@ status = update_dvo_fiximids (subset, Nimage); + if (VERBOSE_IMSTATS) { + CompareImageCounts (subset, Nimage); + } + + if (SUMMARY_IMSTATS) { + SummaryImageStats (subset, Nimage); + } + if (!status) exit (1); exit (0); Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids_client.c =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids_client.c (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids_client.c (revision 35593) @@ -1,8 +1,8 @@ -# include "setposangle.h" +# include "fiximids.h" -// setposangle_client is run on a remote host and is responsible for updating the catalogs +// fiximids_client is run on a remote host and is responsible for updating the catalogs // owned by that host. -// setposangle_client -hostname (hostname) -catdir (catdir) -hostdir (hostdir) -images (images) +// fiximids_client -hostname (hostname) -catdir (catdir) -hostdir (hostdir) -images (images) // load the SkyTable and HostTable from (catdir) [contains the host responsibilities] // load the image data from (images) [contains the zp infomation to apply] @@ -13,5 +13,5 @@ // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, IMAGES) - initialize_setposangle_client (argc, argv); + initialize_fiximids_client (argc, argv); // load the image subset table from the specified location @@ -19,5 +19,5 @@ ImageSubset *subset = ImageSubsetLoad (IMAGES, &Nimage); - update_dvo_setposangle (subset, Nimage); + update_dvo_fiximids (subset, Nimage); exit (0); Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/initialize_fiximids.c =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/initialize_fiximids.c (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/initialize_fiximids.c (revision 35593) @@ -54,4 +54,16 @@ if ((N = get_argument (argc, argv, "-v"))) { VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + SUMMARY_IMSTATS = FALSE; + if ((N = get_argument (argc, argv, "-imstats"))) { + SUMMARY_IMSTATS = TRUE; + remove_argument (N, &argc, argv); + } + VERBOSE_IMSTATS = FALSE; + if ((N = get_argument (argc, argv, "-vimstats"))) { + SUMMARY_IMSTATS = TRUE; + VERBOSE_IMSTATS = TRUE; remove_argument (N, &argc, argv); } @@ -168,4 +180,12 @@ } + VERBOSE_IMSTATS = FALSE; + IMSTATS_FILE = NULL; + if ((N = get_argument (argc, argv, "-imstats"))) { + remove_argument (N, &argc, argv); + IMSTATS_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + UPDATE = FALSE; if ((N = get_argument (argc, argv, "-update"))) { Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/load_images_fiximids.c =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/load_images_fiximids.c (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/load_images_fiximids.c (revision 35593) @@ -1,10 +1,15 @@ # include "fiximids.h" - -off_t Nimage; -e_time *startImage; -e_time *stopImage; -short *photcodeImage; - -Image *load_images_fiximids (FITS_DB *db, off_t *N) { +# define FT_BZERO_INT32 1.0*0x80000000 + +off_t Nimage = 0; +e_time *startImage = NULL; +e_time *stopImage = NULL; +short *photcodeImage = NULL; +unsigned int *imageID = NULL; +unsigned int *imageSeq = NULL; +int *Nimage_valid = NULL; +int *Nimage_invalid = NULL; + +Image *load_images_fiximids (FITS_DB *db, off_t *Nimage_load) { Image *image; @@ -16,5 +21,5 @@ /* use a vtable to keep the images to be calibrated */ - image = gfits_table_get_Image (&db[0].ftable, Nimage, &db[0].swapped); + image = gfits_table_get_Image (&db[0].ftable, Nimage_load, &db[0].swapped); if (!image) { fprintf (stderr, "ERROR: failed to read images\n"); @@ -22,9 +27,45 @@ } - // XXX init lookup tables - - fprintf (stderr, "loaded "OFF_T_FMT" images\n", *Nimage); + fprintf (stderr, "loaded "OFF_T_FMT" images\n", *Nimage_load); return (image); +} + +// find mosaic frames (unique time periods & photcode name matches mosaic) +void initImageIndex (ImageSubset *image, off_t Nimage_init) { + + off_t i; + + Nimage = Nimage_init; + ALLOCATE (startImage, e_time, Nimage); + ALLOCATE (stopImage, e_time, Nimage); + ALLOCATE (photcodeImage, short, Nimage); + ALLOCATE (imageID, unsigned int, Nimage); + ALLOCATE (imageSeq, unsigned int, Nimage); + + ALLOCATE (Nimage_valid, int, Nimage); + ALLOCATE (Nimage_invalid, int, Nimage); + + memset (Nimage_valid, 0, Nimage*sizeof(int)); + memset (Nimage_invalid, 0, Nimage*sizeof(int)); + + // save the image data in the local static arrays + for (i = 0; i < Nimage; i++) { + + /* set image time range (small boundary buffer) */ + e_time start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); + e_time stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + + startImage[i] = start; + stopImage[i] = stop; + photcodeImage[i] = image[i].photcode; + imageID[i] = image[i].imageID; + imageSeq[i] = i; + } + + // sort the index, start, and stop by the start times: + sort_image_times (startImage, stopImage, photcodeImage, imageID, imageSeq, Nimage); + + return; } @@ -32,41 +73,113 @@ // This code is copied from relastro/src/MosaicOps.c -off_t FindImageID (e_time time, short photcode) { - - // use bisection to find the overlapping imagea - - off_t Nlo, Nhi, N; - - // find the last mosaic before start - Nlo = 0; Nhi = Nmosaic; +// find the image by time. this search makes the following assumptions: +// 1) there is only a single EXPOSURE from a single camera with the given exposure time +// 2) there may be multiple IMAGES from that exposure, distinguished by photcode +// 3) the image data has been sorted by start time +// 4) once we have (time > start[i]) we are past the images of interest +int FindImageID (off_t *ID, off_t *Seq, e_time time, short photcode) { + + off_t Nlo, Nhi, N, Ni; + + *ID = 0; + *Seq = -1; + + // Find the last image ending BEFORE 'time' (time > start[i]). + // Use bisection to find the overlapping images. + Nlo = 0; Nhi = Nimage; while (Nhi - Nlo > 10) { N = 0.5*(Nlo + Nhi); - if (startMos[N] < start) { + if (stopImage[N] < time) { Nlo = MAX(N, 0); } else { - Nhi = MIN(N + 1, Nmosaic); + Nhi = MIN(N + 1, Nimage); } } - // check for the matched mosaic starting from Nlo - // we may have to go much beyond Nlo since stop is not sorted - // can we use a sorted version of stop to check when we are beyond the valid range?? - for (N = Nlo; N < Nmosaic; N++) { - if (stop < stopMos[N]) continue; - if (start > startMos[N]) continue; - if (stop < startMos[N]) return (-1); - return (indexMos[N]); - } - - return (-1); + // check for the matched images starting from Nlo + for (N = Nlo; N < Nimage; N++) { + if (time > stopImage[N]) continue; + if (time < startImage[N]) return FALSE; + // a possible image match (start <= time <= stop) + for (Ni = N; (Ni < Nimage) && (time <= stopImage[Ni]); Ni++) { + if (photcodeImage[Ni] == photcode) { + *ID = imageID[Ni]; + *Seq = imageSeq[Ni]; + return TRUE; + } + } + return FALSE; + } + return FALSE; +} + +// increment the count of valid detections for this image +void BumpValidImage (int Seq) { + myAssert (Seq >= 0, "impossible seq"); + myAssert (Seq < Nimage, "impossible seq"); + Nimage_valid[Seq] ++; +} + +// increment the count of invalid detections for this image +void BumpInvalidImage (int Seq) { + myAssert (Seq >= 0, "impossible seq"); + myAssert (Seq < Nimage, "impossible seq"); + Nimage_invalid[Seq] ++; +} + +void CompareImageCounts (ImageSubset *image, off_t Nimage_comp) { + + off_t i; + + for (i = 0; i < Nimage_comp; i++) { + // XXX can I limit the images to those in range of the user region? + + if (Nimage_valid[i] + Nimage_invalid[i] == 0) continue; // XXX for now, skip (optionally?) images with no matched detections + if (image[i].photcode == 0) continue; // skip images with 0 photcode (eg, PHU) + + char *date = ohana_sec_to_date (image[i].tzero); + char *code = GetPhotcodeNamebyCode (image[i].photcode); + fprintf (stderr, "image %s %s : %d %d %d : %f\n", date, code, image[i].nstar, Nimage_valid[i], Nimage_invalid[i], + (Nimage_valid[i] + Nimage_invalid[i]) / (float) image[i].nstar); + } +} + +void SummaryImageStats (ImageSubset *image, off_t Nimage_comp) { + + off_t i; + + // let's report the following: + // Nimages with corrected image IDs + // Nimages with Nsum / Nstar < 0.95 + // Nimages with Nsum / Nstar < 0.99 + + int Nfixed = 0; + int Nbad_1 = 0; + int Nbad_2 = 0; + + for (i = 0; i < Nimage_comp; i++) { + if (Nimage_valid[i] + Nimage_invalid[i] == 0) continue; // XXX for now, skip (optionally?) images with no matched detections + if (image[i].photcode == 0) continue; // skip images with 0 photcode (eg, PHU) + + if (Nimage_invalid[i] > 0) Nfixed ++; + + float found_ratio = (Nimage_valid[i] + Nimage_invalid[i]) / (float) image[i].nstar; + + if (found_ratio < 0.99) Nbad_1 ++; + if (found_ratio < 0.95) Nbad_2 ++; + } + + fprintf (stderr, "Nfixed: %d, Ndrop (>1%%) %d, Ndrop (>5%%) %d\n", Nfixed, Nbad_1, Nbad_2); } // sort two times vectors and an index by first time vector -void sort_mosaic_times (unsigned int *S, unsigned int *E, off_t *I, off_t N) { - -# define SWAPFUNC(A,B){ unsigned int tmp_t; off_t tmp_i; \ +void sort_image_times (e_time *S, e_time *E, short *C, unsigned int *I, unsigned int *Q, off_t N) { + +# define SWAPFUNC(A,B){ e_time tmp_t; unsigned int tmp_u; short tmp_c; \ tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \ tmp_t = E[A]; E[A] = E[B]; E[B] = tmp_t; \ - tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \ + tmp_c = C[A]; C[A] = C[B]; C[B] = tmp_c; \ + tmp_u = I[A]; I[A] = I[B]; I[B] = tmp_u; \ + tmp_u = Q[A]; Q[A] = Q[B]; Q[B] = tmp_u; \ } # define COMPARE(A,B)(S[A] < S[B]) @@ -79,95 +192,118 @@ } -// use the time to make the match, but use bracketing to make it faster: - -// find mosaic frames (unique time periods & photcode name matches mosaic) -void initMosaics (ImageSubset *image, off_t Nimage) { - - off_t i, Nmos, NMOSAIC; - unsigned int start, stop; - - unsigned int *startMos, *stopMos; - off_t *indexMos; - - Nmosaic = 0; - NMOSAIC = 10; - ALLOCATE (mosaic, Mosaic, NMOSAIC); - - ALLOCATE (startMos, unsigned int, NMOSAIC); - ALLOCATE (stopMos, unsigned int, NMOSAIC); - ALLOCATE (indexMos, off_t, NMOSAIC); - - /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */ - for (i = 0; i < Nimage; i++) { - if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue; - - /* set image time range */ - start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); - stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); - - /* a new mosaic, define ranges */ - mosaic[Nmosaic].start = start; - mosaic[Nmosaic].stop = stop; - mosaic[Nmosaic].coords = image[i].coords; - - startMos[Nmosaic] = start; - stopMos[Nmosaic] = stop; - indexMos[Nmosaic] = Nmosaic; - - Nmosaic ++; - if (Nmosaic == NMOSAIC) { - NMOSAIC += 10; - REALLOCATE (mosaic, Mosaic, NMOSAIC); - REALLOCATE (startMos, unsigned int, NMOSAIC); - REALLOCATE (stopMos, unsigned int, NMOSAIC); - REALLOCATE (indexMos, off_t, NMOSAIC); - } - } - - // sort the index, start, and stop by the start times: - sort_mosaic_times (startMos, stopMos, indexMos, Nmosaic); - - // array to store image->mosaic index - Nmosaic_for_images = Nimage; - ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images); - - /* now assign the WRP images to these mosaics */ - for (i = 0; i < Nimage; i++) { - mosaic_for_images[i] = -1; // default value for no mosaic found - - if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue; - - /* set image time range */ - start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); - stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); - - Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos); - if (Nmos == -1) { - fprintf (stderr, "cannot match mosaic for %d (ID = %d)\n", (int) i, (int) image[i].imageID); - continue; - } - - // mosaic corresponding to this image - mosaic_for_images[i] = Nmos; - } - - free (startMos); - free (stopMos); - free (indexMos); - return; -} - -Mosaic *getMosaicForImage (off_t im) { - - off_t mos; - - if (im < 0) abort(); - if (im >= Nmosaic_for_images) abort(); - - // search for the mosaic that - mos = mosaic_for_images[im]; - if (mos < 0) return NULL; - - return &mosaic[mos]; -} - +# 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"); + +int ImageValidLoad(char *filename) { + + int Ncol; + off_t i; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return FALSE; + } + + /* 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 FALSE; + } + 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 FALSE; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + if (VERBOSE) fprintf (stderr, "can't read table header\n"); + fclose (f); + return FALSE; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + if (VERBOSE) fprintf (stderr, "can't read table data\n"); + fclose (f); + return FALSE; + } + fclose (f); + + char type[16]; + + GET_COLUMN (Nvalid, "NVALID", int); + GET_COLUMN (Ninvalid, "NINVALID", int); + + myAssert (Nrow == Nimage, "wrong size imstats result file\n"); + + // merge the new arrays into the single array? + for (i = 0; i < Nrow; i++) { + Nimage_valid[i] += Nvalid[i]; + Nimage_invalid[i] += Ninvalid[i]; + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (Nvalid); + free (Ninvalid); + + return TRUE; +} + +int ImageValidSave(char *filename) { + + 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_VALID"); + + // create the table layout + gfits_define_bintable_column (&theader, "J", "NVALID", "tmp", "tmp", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "NINVALID", "tmp", "tmp", 1.0, FT_BZERO_INT32); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "NVALID", Nimage_valid, Nimage); + gfits_set_bintable_column (&theader, &ftable, "NINVALID", Nimage_invalid, Nimage); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename); + return FALSE; + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + fflush (f); + + return TRUE; +} Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_catalog_fiximids.c =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_catalog_fiximids.c (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_catalog_fiximids.c (revision 35593) @@ -1,10 +1,12 @@ # include "fiximids.h" -void update_catalog_fiximids (Catalog *catalog, ImageSubset *image, off_t *index, off_t Nimage) { +void update_catalog_fiximids (Catalog *catalog) { - float posAngle, pltScale; - off_t i, j, found; + off_t i; - found = 0; + off_t found = 0; + off_t Nvalid = 0; + off_t Ninvalid = 0; + for (i = 0; i < catalog[0].Nmeasure; i++) { @@ -17,19 +19,20 @@ e_time time = catalog[0].measure[i].t; - off_t imageID = FindImageID (time, photcode); - if (imageID <= 0) continue; // detections with imageID == 0 do not have a valid image (eg, ref photcode) + // skip detections with no valid imageID (eg, ref photcode) + + off_t imageID, Seq; + if (!FindImageID (&imageID, &Seq, time, photcode)) continue; + if (imageID <= 0) continue; if (imageID == catalog[0].measure[i].imageID) { // existing image ID is valid Nvalid ++; - onImage[imageID] ++; - continue; + BumpValidImage (Seq); + } else { + // existing image ID is NOT valid + Ninvalid ++; + catalog[0].measure[i].imageID = imageID; + BumpInvalidImage (Seq); } - - // existing image ID is NOT valid - catalog[0].measure[i].imageID = imageID; - Ninvalid ++; - onImage[imageID] ++; - continue; found ++; @@ -37,5 +40,5 @@ if (found) { - fprintf (stderr, "found "OFF_T_FMT" matches\n", found); + fprintf (stderr, "found "OFF_T_FMT" matches: "OFF_T_FMT" valid, "OFF_T_FMT" invalid\n", found, Nvalid, Ninvalid); } } Index: branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_dvo_fiximids.c =================================================================== --- branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_dvo_fiximids.c (revision 35592) +++ branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_dvo_fiximids.c (revision 35593) @@ -6,5 +6,5 @@ SkyList *skylist = NULL; Catalog catalog; - off_t i, maxID, *index; + off_t i; // load the current sky table (layout of all SkyRegions) @@ -12,25 +12,11 @@ SkyTableSetFilenames (sky, CATDIR, "cpt"); + // generate index tables + initImageIndex (image, Nimage); + if (PARALLEL && !HOST_ID) { update_dvo_fiximids_parallel (sky, image, Nimage); return TRUE; } - - // create an index for the image IDs - maxID = 0; - for (i = 0; i < Nimage; i++) { - maxID = MAX(maxID, image[i].imageID); - } - ALLOCATE (index, off_t, maxID + 1); - for (i = 0; i < maxID + 1; i++) { - index[i] = -1; - } - for (i = 0; i < Nimage; i++) { - if (!image[i].imageID) continue; // images with ID == 0 are virtual - index[image[i].imageID] = i; - } - - // set up the links between images and mosaics - initMosaics (image, Nimage); // determine the populated SkyRegions overlapping the requested area (default depth) @@ -69,5 +55,5 @@ } - update_catalog_fiximids (&catalog, image, index, Nimage); + update_catalog_fiximids (&catalog); if (!UPDATE) { @@ -83,4 +69,10 @@ dvo_catalog_free (&catalog); } + + // for the remote clients, we need to dump out the results table (Nvalid, Ninvalid) + if (HOST_ID && IMSTATS_FILE) { + ImageValidSave (IMSTATS_FILE); + } + return (TRUE); } @@ -90,7 +82,12 @@ int update_dvo_fiximids_parallel (SkyTable *sky, ImageSubset *image, off_t Nimage) { + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + // write out the subset table of image information char imageFile[512]; - snprintf (imageFile, 512, "%s/Images.subset.dat", CATDIR); + snprintf (imageFile, 512, "%s/fiximids.images.%s.dat", CATDIR, uniquer); if (!ImageSubsetSave (imageFile, image, Nimage)) { @@ -112,15 +109,21 @@ table->hosts[i].pathname = tmppath; + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/fiximids.results.%s.dat", table->hosts[i].pathname, uniquer); + char command[1024]; snprintf (command, 1024, "fiximids_client -hostID %d -catdir %s -hostdir %s -images %s -region %f %f %f %f", - table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, imageFile, + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, imageFile, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); - fprintf (stderr, "command: %s\n", command); + int saveImstats = VERBOSE_IMSTATS || SUMMARY_IMSTATS; char tmpline[1024]; - if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } - if (SINGLE_CPT) { snprintf (tmpline, 1024, "%s -cpt %s", SINGLE_CPT, command); strcpy (command, tmpline); } - if (UPDATE) { snprintf (tmpline, 1024, "%s -update", command); strcpy (command, tmpline); } + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (saveImstats) { snprintf (tmpline, 1024, "%s -imstats %s", command, table->hosts[i].results); strcpy (command, tmpline); } + if (SINGLE_CPT) { snprintf (tmpline, 1024, "%s -cpt %s", command, SINGLE_CPT); strcpy (command, tmpline); } + if (UPDATE) { snprintf (tmpline, 1024, "%s -update", command); strcpy (command, tmpline); } + + fprintf (stderr, "command: %s\n", command); if (PARALLEL_MANUAL) continue; @@ -154,4 +157,17 @@ } + if (VERBOSE_IMSTATS || SUMMARY_IMSTATS) { + for (i = 0; i < table->Nhosts; i++) { + while (!ImageValidLoad (table->hosts[i].results)) { + // failed to get the data from this host. This can happen for various reasons. Give the user a chance to try again... + fprintf (stderr, "failed to read data from %s, stopping operations until this can be fixed\n", table->hosts[i].hostname); + fprintf (stderr, "you may run the command manually and send this process the CONT signal\n"); + int pid = getpid(); + kill (pid, SIGSTOP); + fprintf (stderr, "retrying %s\n", table->hosts[i].results); + } + } + } + return (TRUE); }