Index: /trunk/Ohana/src/delstar/Makefile =================================================================== --- /trunk/Ohana/src/delstar/Makefile (revision 35757) +++ /trunk/Ohana/src/delstar/Makefile (revision 35758) @@ -30,4 +30,8 @@ $(SRC)/delete_imagefile.$(ARCH).o \ $(SRC)/delete_duplicate_images.$(ARCH).o \ +$(SRC)/delete_duplicate_measures.$(ARCH).o \ +$(SRC)/delete_fix_LAP.$(ARCH).o \ +$(SRC)/delete_fix_LAP_edges.$(ARCH).o \ +$(SRC)/delete_fix_LAP_edges_delete.$(ARCH).o \ $(SRC)/delete_image_photcodes.$(ARCH).o \ $(SRC)/delete_photcodes.$(ARCH).o \ @@ -39,4 +43,7 @@ $(SRC)/parse_time.$(ARCH).o \ $(SRC)/ImageIDsOps.$(ARCH).o \ +$(SRC)/ImageOpsFixLAP.$(ARCH).o \ +$(SRC)/ImageSubsetFixLAP.$(ARCH).o \ +$(SRC)/MeasureEdgeOps.$(ARCH).o \ $(SRC)/check_permissions.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o @@ -48,7 +55,14 @@ $(SRC)/args.$(ARCH).o \ $(SRC)/delete_duplicate_images.$(ARCH).o \ +$(SRC)/delete_duplicate_measures.$(ARCH).o \ +$(SRC)/delete_fix_LAP.$(ARCH).o \ +$(SRC)/delete_fix_LAP_edges.$(ARCH).o \ +$(SRC)/delete_fix_LAP_edges_delete.$(ARCH).o \ $(SRC)/delete_photcodes.$(ARCH).o \ $(SRC)/delete_photcodes_single.$(ARCH).o \ $(SRC)/ImageIDsOps.$(ARCH).o \ +$(SRC)/ImageOpsFixLAP.$(ARCH).o \ +$(SRC)/ImageSubsetFixLAP.$(ARCH).o \ +$(SRC)/MeasureEdgeOps.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o Index: /trunk/Ohana/src/delstar/include/delstar.h =================================================================== --- /trunk/Ohana/src/delstar/include/delstar.h (revision 35757) +++ /trunk/Ohana/src/delstar/include/delstar.h (revision 35758) @@ -21,4 +21,27 @@ } IndexArray; +// to determine the image ID, we need time (tzero + trate) and photcode +typedef struct { + unsigned int imageID; + unsigned int externID; + 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) + e_time tmin; // start exp - 1 + e_time tmax; // stop exp + 1 + unsigned char trate; // scan rate (100 usec/pixel) + short photcode; +} ImageSubset; + +typedef struct { + double R; + double D; + int objID; + int catID; + int detID; + int imageID; +} MeasureEdge; + /* global variables set in parameter file */ char ImageCat[DVO_MAX_PATH]; @@ -33,4 +56,8 @@ char *IMAGENAME; +char *IMAGES; +char *IMSTATS_FILE; +char *MEASURE_EDGE_FILE; +char *EDGE_DELETIONS; double NSIGMA; @@ -52,4 +79,5 @@ int HOST_ID; char *HOSTDIR; +char *SINGLE_CPT; int SKIP_IMAGES; @@ -60,7 +88,8 @@ char *PHOTCODE_LIST; +char *UNIQUER; int MODE; -enum {MODE_NONE, MODE_IMAGENAME, MODE_IMAGEFILE, MODE_TIME, MODE_ORPHAN, MODE_MISSED, MODE_PHOTCODES, MODE_DUP_IMAGES}; +enum {MODE_NONE, MODE_IMAGENAME, MODE_IMAGEFILE, MODE_TIME, MODE_ORPHAN, MODE_MISSED, MODE_PHOTCODES, MODE_DUP_IMAGES, MODE_DUP_MEASURES, MODE_FIX_LAP, MODE_FIX_LAP_STATS, MODE_FIX_LAP_EDGES, MODE_FIX_LAP_EDGES_DELETE}; char DateKeyword[64], DateMode[64], UTKeyword[64], MJDKeyword[64], JDKeyword[64]; @@ -109,5 +138,5 @@ void SortAveMeasMatch (off_t *MEAS, off_t *AVE, off_t N); -int delete_photcodes (int hostID, char *hostpath); +int delete_photcodes (); int delete_photcodes_parallel (SkyList *sky); int delete_photcodes_catalog (Catalog *catalog, PhotCode **photcodes, int Nphotcodes); @@ -115,6 +144,6 @@ int delete_photcodes_single (char *cptname); -int delete_duplicate_images (int hostID, char *hostpath, FITS_DB *db); -int delete_duplicate_image_measures (int hostID, char *hostpath, IndexArray *imageID); +int delete_duplicate_images (FITS_DB *db); +int delete_duplicate_image_measures (IndexArray *imageID); int delete_duplicate_image_measures_parallel (SkyList *sky, IndexArray *imageID); int delete_duplicate_image_measures_catalog (Catalog *catalog, IndexArray *imageID); @@ -125,2 +154,43 @@ IndexArray *ImageIDLoad(char *filename); int ImageIDSave(char *filename, IndexArray *imageID); + +int delete_duplicate_measures (); +int delete_duplicate_measures_parallel (SkyList *sky); +int delete_duplicate_measures_catalog (Catalog *catalog); + +int delete_fix_LAP (ImageSubset *image, off_t Nimage); +int delete_fix_LAP_catalog (Catalog *catalog, ImageSubset *image, off_t Nimage); +int delete_fix_LAP_measures (off_t *measureDrop, Catalog *catalog, ImageSubset *image, off_t Nimage); +int delete_fix_LAP_parallel (SkyList *sky, ImageSubset *image, off_t Nimage); +int delete_fix_LAP_setstats (ImageSubset *image, off_t Nimage); + +void initImageIndex (ImageSubset *image, off_t Nimage_init); +int *getImageIndex (int *maxIndex); +int FindIDexp (off_t *ID, off_t *Seq, e_time time, short photcode); +int FindIDstk (off_t *ID, off_t *Seq, short *photcode, off_t extID); + +ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage); +int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage); +ImageSubset *ImagesToSubset (Image *image, off_t N); + +void sort_fullIDs (uint64_t *I, int *S, off_t N); + +int ImageValidSave(char *filename); +int ImageValidLoad(char *filename); +void SummaryImageStats (ImageSubset *image, off_t Nimage_comp); +void SummaryImageStats_Full (Image *image, off_t Nimage_comp); +void BumpInvalidImage (int Seq); +void BumpValidImage (int Seq); + +MeasureEdge *delete_fix_LAP_edges (off_t *nmeasure_edge); +MeasureEdge *delete_fix_LAP_edges_parallel (SkyList *sky, off_t *nmeasure_edge); +MeasureEdge *delete_fix_LAP_edges_get_measures (Catalog *catalog, MeasureEdge *measure_edge, off_t *Nmeasure_edge_in, off_t *NMEASURE_EDGE_IN); +int delete_fix_LAP_edges_find_dups (MeasureEdge *measure_edge, off_t Nmeasure_edge); + +MeasureEdge *MeasureEdgeMerge (MeasureEdge *measure_edge_all, off_t *Nmeasure_edge_in, off_t *NMEASURE_EDGE_IN, MeasureEdge *measure_edge, off_t Nmeasure_edge); +int MeasureEdgeSave(char *filename, MeasureEdge *measure_edge, off_t Nmeasure_edge); +MeasureEdge *MeasureEdgeLoad(char *filename, off_t *Nmeasure_edge); + +int delete_fix_LAP_edges_delete (); +int delete_fix_LAP_edges_drop_measures (Catalog *catalog, MeasureEdge *measure_edge, off_t Nmeasure_edge, int catIDcount); +int delete_fix_LAP_edges_delete_parallel (SkyList *sky); Index: /trunk/Ohana/src/delstar/src/ImageOpsFixLAP.c =================================================================== --- /trunk/Ohana/src/delstar/src/ImageOpsFixLAP.c (revision 35758) +++ /trunk/Ohana/src/delstar/src/ImageOpsFixLAP.c (revision 35758) @@ -0,0 +1,470 @@ +# include "delstar.h" +# define FT_BZERO_INT32 1.0*0x80000000 + +off_t Nimage = 0; +e_time *startImage_exp = NULL; +e_time *stopImage_exp = NULL; +short *photcodeImage_exp = NULL; +unsigned int *imageID_exp = NULL; +unsigned int *imageSeq_exp = NULL; + +unsigned int *externID_stk = NULL; +short *photcodeImage_stk = NULL; +unsigned int *imageID_stk = NULL; +unsigned int *imageSeq_stk = NULL; + +int maxID = 0; +int *imageIndex = NULL; + +int *Nimage_valid = NULL; +int *Nimage_invalid = NULL; + +void sort_image_externID (unsigned int *E, short *C, unsigned int *I, unsigned int *Q, off_t N); +void sort_image_times (e_time *S, e_time *E, short *C, unsigned int *I, unsigned int *Q, off_t N); + +// 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_exp, e_time, Nimage); + ALLOCATE (stopImage_exp, e_time, Nimage); + ALLOCATE (photcodeImage_exp, short, Nimage); + ALLOCATE (imageID_exp, unsigned int, Nimage); + ALLOCATE (imageSeq_exp, unsigned int, Nimage); + + ALLOCATE (externID_stk, unsigned int, Nimage); + ALLOCATE (photcodeImage_stk, short, Nimage); + ALLOCATE (imageID_stk, unsigned int, Nimage); + ALLOCATE (imageSeq_stk, unsigned int, Nimage); + + ALLOCATE (Nimage_valid, int, Nimage); + ALLOCATE (Nimage_invalid, int, Nimage); + memset (Nimage_valid, 0, sizeof(int)*Nimage); + memset (Nimage_invalid, 0, sizeof(int)*Nimage); + + // 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_exp[i] = start; + stopImage_exp[i] = stop; + photcodeImage_exp[i] = image[i].photcode; + imageID_exp[i] = image[i].imageID; + imageSeq_exp[i] = i; + + externID_stk[i] = image[i].externID; + photcodeImage_stk[i] = image[i].photcode; + imageID_stk[i] = image[i].imageID; + imageSeq_stk[i] = i; + + maxID = MAX(image[i].imageID, maxID); + } + + // sort the index, start, and stop by the start times: + sort_image_times (startImage_exp, stopImage_exp, photcodeImage_exp, imageID_exp, imageSeq_exp, Nimage); + + // sort the index, start, and stop by the start times: + sort_image_externID (externID_stk, photcodeImage_stk, imageID_stk, imageSeq_stk, Nimage); + + ALLOCATE (imageIndex, int, maxID + 1); + for (i = 0; i < maxID + 1; i++) { + imageIndex[i] = -1; + } + for (i = 0; i < Nimage; i++) { + if (image[i].imageID == 0) continue; + imageIndex[image[i].imageID] = i; + } + + return; +} + +int *getImageIndex (int *maxIDout) { + *maxIDout = maxID; + return imageIndex; +} + +// sort two times vectors and an index by first time vector +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_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]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +// sort two times vectors and an index by first time vector +void sort_image_externID (unsigned int *E, short *C, unsigned int *I, unsigned int *Q, off_t N) { + +# define SWAPFUNC(A,B){ unsigned int tmp_u; short tmp_c; \ + tmp_u = I[A]; I[A] = I[B]; I[B] = tmp_u; \ + tmp_u = E[A]; E[A] = E[B]; E[B] = tmp_u; \ + tmp_c = C[A]; C[A] = C[B]; C[B] = tmp_c; \ + tmp_u = Q[A]; Q[A] = Q[B]; Q[B] = tmp_u; \ +} +# define COMPARE(A,B)(E[A] < E[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +// 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 hav been sorted by start time +// 4) once we have (time < start[i]) we are past the desired exposure +int FindIDexp (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 (stopImage_exp[N] < time) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N + 1, Nimage); + } + } + + // check for the matched images starting from Nlo + for (N = Nlo; N < Nimage; N++) { + if (time > stopImage_exp[N]) continue; + if (time < startImage_exp[N]) return FALSE; + // a possible image match (start <= time <= stop) + for (Ni = N; (Ni < Nimage) && (time <= stopImage_exp[Ni]); Ni++) { + if (photcodeImage_exp[Ni] == photcode) { + *ID = imageID_exp[Ni]; + *Seq = imageSeq_exp[Ni]; + return TRUE; + } + } + return FALSE; + } + return FALSE; +} + +// find the image by extern ID +// 1) extID is unique among images +// 3) the image data have been sorted by extID +// 4) once we have (extID < externID[i]) we are past the desired exposure +int FindIDstk (off_t *ID, off_t *Seq, short *photcode, off_t extID) { + + off_t Nlo, Nhi, N, Ni; + + *ID = 0; + *Seq = -1; + *photcode = 0; + + // Find the last image with externID < extID + // Use bisection to find the overlapping images. + Nlo = 0; Nhi = Nimage; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (externID_stk[N] < extID) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N + 1, Nimage); + } + } + + // check for the matched images starting from Nlo + for (N = Nlo; N < Nimage; N++) { + if (extID > externID_stk[N]) continue; + if (extID < externID_stk[N]) return FALSE; + + // try to find other matches? + *photcode = photcodeImage_stk[N]; + *ID = imageID_stk[N]; + *Seq = imageSeq_stk[N]; + + for (Ni = N + 1; Ni < Nimage; Ni++) { + if (extID < externID_stk[Ni]) break; + if (externID_stk[Ni] == extID) { + fprintf (stderr, "WARNING: found duplicate extern IDs : ext ID %d, image ID 1: %d (%d), image ID 2: %d (%d)\n", externID_stk[N], imageID_stk[N], photcodeImage_stk[N], imageID_stk[Ni], photcodeImage_stk[Ni]); + } + } + return TRUE; + } + 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 Nphu = 0; + int Nexp = 0; + int Nstk = 0; + int Nskip = 0; + int Nbad_1 = 0; + int Nbad_2 = 0; + + for (i = 0; i < Nimage_comp; i++) { + if (image[i].photcode == 0) { + Nphu ++; + continue; // skip images with 0 photcode (eg, PHU) + } + + if ((image[i].photcode > 10000) && (image[i].photcode < 10600)) Nexp ++; + if ((image[i].photcode >= 11000) && (image[i].photcode <= 11600)) Nstk ++; + + if (Nimage_valid[i] + Nimage_invalid[i] == 0) Nskip ++; + + 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, "Nimage: "OFF_T_FMT"\n", Nimage_comp); + fprintf (stderr, "Nfixed: %d\n", Nfixed); + fprintf (stderr, "Nphu: %d\n", Nphu); + fprintf (stderr, "Nexp: %d\n", Nexp); + fprintf (stderr, "Nstk: %d\n", Nstk); + fprintf (stderr, "Nskip: %d\n", Nskip); + fprintf (stderr, "Ndrop (>1%%): %d\n", Nbad_1); + fprintf (stderr, "Ndrop (>5%%): %d\n", Nbad_2); +} + +void SummaryImageStats_Full (Image *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 Nphu = 0; + int Nexp = 0; + int Nstk = 0; + int Nskip = 0; + int Nbad_1 = 0; + int Nbad_2 = 0; + + BuildChipMatch (image, Nimage); + + for (i = 0; i < Nimage_comp; i++) { + if (image[i].photcode == 0) { + Nphu ++; + continue; // skip images with 0 photcode (eg, PHU) + } + + if ((image[i].photcode > 10000) && (image[i].photcode < 10600)) Nexp ++; + if ((image[i].photcode >= 11000) && (image[i].photcode <= 11600)) Nstk ++; + + float found_ratio = (Nimage_valid[i] + Nimage_invalid[i]) / (float) image[i].nstar; + + if ((Nimage_valid[i] + Nimage_invalid[i] == 0) || (found_ratio < 0.995)) { + double r, d; + if (!FindMosaicForImage (image, Nimage, i)) { + fprintf (stderr, "failed to get astrometry for %s\n", image[i].name); + } else { + if (!strcmp(&image[i].coords.ctype[4], "-DIS")) { + XY_to_RD (&r, &d, 0.0, 0.0, &image[i].coords); + } else { + XY_to_RD (&r, &d, 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords); + } + } + char *date = ohana_sec_to_date (image[i].tzero); + fprintf (stderr, "skip : %s %d %10.6f %10.6f %5.3f %s\n", date, image[i].photcode, r, d, found_ratio, image[i].name); + free (date); + Nskip ++; + } + + if (Nimage_invalid[i] > 0) Nfixed ++; + + + if (found_ratio < 0.99) Nbad_1 ++; + if (found_ratio < 0.95) Nbad_2 ++; + } + + fprintf (stderr, "Nimage: "OFF_T_FMT"\n", Nimage_comp); + fprintf (stderr, "Nfixed: %d\n", Nfixed); + fprintf (stderr, "Nphu: %d\n", Nphu); + fprintf (stderr, "Nexp: %d\n", Nexp); + fprintf (stderr, "Nstk: %d\n", Nstk); + fprintf (stderr, "Nskip: %d\n", Nskip); + fprintf (stderr, "Ndrop (>1%%): %d\n", Nbad_1); + fprintf (stderr, "Ndrop (>5%%): %d\n", Nbad_2); +} + +# 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: /trunk/Ohana/src/delstar/src/ImageSubsetFixLAP.c =================================================================== --- /trunk/Ohana/src/delstar/src/ImageSubsetFixLAP.c (revision 35758) +++ /trunk/Ohana/src/delstar/src/ImageSubsetFixLAP.c (revision 35758) @@ -0,0 +1,216 @@ +# include "delstar.h" +# define FT_BZERO_INT16 1.0*0x8000 +# define FT_BZERO_INT32 1.0*0x80000000 + +# 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"); + +ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage) { + + int Ncol; + off_t i; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nimage = 0; + ImageSubset *image = 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); + + char type[16]; + + GET_COLUMN (imageID, "IMAGE_ID", int); + GET_COLUMN (externID, "EXTERN_ID", int); + 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); + + ALLOCATE (image, ImageSubset, Nrow); + for (i = 0; i < Nrow; i++) { + image[i].imageID = imageID[i]; + image[i].externID = externID[i]; + 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]; + image[i].photcode = photcode[i]; + image[i].tmin = image[i].tzero - 1; + image[i].tmax = image[i].tzero + image[i].trate*image[i].NY + 1; + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (imageID); + free (externID); + free (NX); + free (NY); + free (nstar); + free (tzero); + free (trate); + free (photcode); + + *nimage = Nrow; + return image; +} + +int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage) { + + off_t 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_SUBSET"); + + // 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", "EXTERN_ID", "extern ID", "tmp", 1.0, FT_BZERO_INT32); + 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); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + unsigned int *imageID, *externID, *nstar; + short *NX, *NY; + e_time *tzero; + short *trate; + short *photcode; + + // create intermediate storage arrays + ALLOCATE (imageID, unsigned int, Nimage); + ALLOCATE (externID, unsigned 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 + for (i = 0; i < Nimage; i++) { + imageID[i] = image[i].imageID ; + externID[i] = image[i].externID; + 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 ; + photcode[i] = image[i].photcode; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, Nimage); + gfits_set_bintable_column (&theader, &ftable, "EXTERN_ID", externID, Nimage); + 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); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, Nimage); + + free (imageID ); + free (externID); + free (NX ); + free (NY ); + free (nstar ); + free (tzero ); + free (trate ); + free (photcode); + + 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; +} + +ImageSubset *ImagesToSubset (Image *image, 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 + ImageSubset *subset = NULL; + ALLOCATE (subset, ImageSubset, N); + + for (i = 0; i < N; i++) { + subset[i].imageID = image[i].imageID ; + subset[i].externID = image[i].externID ; + 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].tmin = image[i].tzero - 1; + subset[i].tmax = image[i].tzero + image[i].trate*image[i].NY + 1; + subset[i].photcode = image[i].photcode ; + } + return subset; +} + + Index: /trunk/Ohana/src/delstar/src/MeasureEdgeOps.c =================================================================== --- /trunk/Ohana/src/delstar/src/MeasureEdgeOps.c (revision 35758) +++ /trunk/Ohana/src/delstar/src/MeasureEdgeOps.c (revision 35758) @@ -0,0 +1,193 @@ +# include "delstar.h" +# define FT_BZERO_INT32 1.0*0x80000000 + +# 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"); + +MeasureEdge *MeasureEdgeLoad(char *filename, off_t *Nmeasure_edge) { + + int Ncol; + off_t i; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *Nmeasure_edge= 0; + MeasureEdge *measure_edge = NULL; + + 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 (R , "RA", double); + GET_COLUMN (D , "DEC", double); + GET_COLUMN (objID , "OBJ_ID", int); + GET_COLUMN (catID , "CAT_ID", int); + GET_COLUMN (detID , "DET_ID", int); + GET_COLUMN (imageID, "IMAGE_ID", int); + + ALLOCATE (measure_edge, MeasureEdge, Nrow); + + // merge the new arrays into the single array? + for (i = 0; i < Nrow; i++) { + measure_edge[i].R = R[i] ; + measure_edge[i].D = D[i] ; + measure_edge[i].objID = objID[i] ; + measure_edge[i].catID = catID[i] ; + measure_edge[i].detID = detID[i] ; + measure_edge[i].imageID = imageID[i]; + } + fprintf (stderr, "loaded data for %lld measures\n", (long long) Nrow); + + free (R ); + free (D ); + free (objID ); + free (catID ); + free (detID ); + free (imageID ); + + *Nmeasure_edge = Nrow; + return measure_edge; +} + +int MeasureEdgeSave(char *filename, MeasureEdge *measure_edge, off_t Nmeasure_edge) { + + 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", "MEASURE_EDGE"); + + // create the table layout + gfits_define_bintable_column (&theader, "D", "RA", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "DET_ID", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "tmp", "tmp", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + double *R, *D; + int *imageID, *objID, *catID, *detID; + + ALLOCATE (R , double, Nmeasure_edge); + ALLOCATE (D , double, Nmeasure_edge); + ALLOCATE (objID , int, Nmeasure_edge); + ALLOCATE (catID , int, Nmeasure_edge); + ALLOCATE (detID , int, Nmeasure_edge); + ALLOCATE (imageID, int, Nmeasure_edge); + + // assign the storage arrays + off_t i; + for (i = 0; i < Nmeasure_edge; i++) { + R[i] = measure_edge[i].R; + D[i] = measure_edge[i].D; + objID[i] = measure_edge[i].objID; + catID[i] = measure_edge[i].catID; + detID[i] = measure_edge[i].detID; + imageID[i] = measure_edge[i].imageID; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, Nmeasure_edge); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, Nmeasure_edge); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeasure_edge); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeasure_edge); + gfits_set_bintable_column (&theader, &ftable, "DET_ID", detID, Nmeasure_edge); + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, Nmeasure_edge); + + free (R ); + free (D ); + free (objID ); + free (catID ); + free (detID ); + free (imageID ); + + 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; +} + +MeasureEdge *MeasureEdgeMerge (MeasureEdge *measure_edge_all, off_t *Nmeasure_edge_in, off_t *NMEASURE_EDGE_IN, MeasureEdge *measure_edge, off_t Nmeasure_edge) { + + off_t Nmeasure_edge_all = *Nmeasure_edge_in; + off_t NMEASURE_EDGE_ALL = *NMEASURE_EDGE_IN; + + off_t Nstart = Nmeasure_edge_all; + + Nmeasure_edge_all += Nmeasure_edge; + NMEASURE_EDGE_ALL += Nmeasure_edge; + REALLOCATE (measure_edge_all, MeasureEdge, NMEASURE_EDGE_ALL); + + off_t i; + for (i = 0; i < Nmeasure_edge; i++) { + measure_edge_all[i + Nstart] = measure_edge[i]; + } + + *Nmeasure_edge_in = Nmeasure_edge_all; + *NMEASURE_EDGE_IN = NMEASURE_EDGE_ALL; + return measure_edge_all; +} + Index: /trunk/Ohana/src/delstar/src/Shutdown.c =================================================================== --- /trunk/Ohana/src/delstar/src/Shutdown.c (revision 35757) +++ /trunk/Ohana/src/delstar/src/Shutdown.c (revision 35758) @@ -22,5 +22,5 @@ SetProtect (TRUE); - gfits_db_close (db); + if (db) gfits_db_close (db); fprintf (stderr, "ERROR: delstar halted\n"); exit (1); Index: /trunk/Ohana/src/delstar/src/args.c =================================================================== --- /trunk/Ohana/src/delstar/src/args.c (revision 35757) +++ /trunk/Ohana/src/delstar/src/args.c (revision 35758) @@ -11,8 +11,9 @@ fprintf (stderr, " delstar -photcodes (list) : delete by photcode\n\n"); fprintf (stderr, " delstar -dup-images : delete duplicate images (by externID)\n\n"); + fprintf (stderr, " delstar -dup-measures : delete duplicate measures (by imageID + detID)\n\n"); fprintf (stderr, " optional flags:\n"); fprintf (stderr, " -v : verbose mode\n"); fprintf (stderr, " -photcode (code) : restrict by photcode\n"); - fprintf (stderr, " -update : apply changes (-dup-images only)\n"); + fprintf (stderr, " -update : apply changes (-dup-images or -dup-measures only)\n"); fprintf (stderr, " -image-details : list info about the deleted images (-dup-images only)\n"); fprintf (stderr, " -image-only : only examine the image table (changes are NOT saved; -dup-images only)\n"); @@ -70,10 +71,21 @@ } - // we should never delete the images before the measures (or we won't know what to delete) - // -image-only is for testing + // We should generally not delete the images before the measures (or we won't know what + // to delete). -image-only is for testing, -image-only-force -update WILL delete the images IMAGE_ONLY = FALSE; if ((N = get_argument (argc, argv, "-image-only"))) { IMAGE_ONLY = TRUE; UPDATE = FALSE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-image-only-force"))) { + IMAGE_ONLY = TRUE; + remove_argument (N, &argc, argv); + } + + SINGLE_CPT = NULL; + if ((N = get_argument (argc, argv, "-cpt"))) { + remove_argument (N, &argc, argv); + SINGLE_CPT = strcreate (argv[N]); remove_argument (N, &argc, argv); } @@ -104,4 +116,5 @@ // XXX for the moment, make this selection manual. it needs to be automatic // based on the state of the SkyTable + HOST_ID = 0; PARALLEL = FALSE; if ((N = get_argument (argc, argv, "-parallel"))) { @@ -160,4 +173,46 @@ MODE = MODE_DUP_IMAGES; remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-dup-measures"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_DUP_MEASURES; + remove_argument (N, &argc, argv); + SKIP_IMAGES = TRUE; // we do not need to load the images for -dup-measures + } + if ((N = get_argument (argc, argv, "-fix-LAP"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_FIX_LAP; + remove_argument (N, &argc, argv); + } + UNIQUER = NULL; + if ((N = get_argument (argc, argv, "-fix-LAP-imstats"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_FIX_LAP_STATS; + remove_argument (N, &argc, argv); + if (N == argc) { + fprintf (stderr, "USAGE: delstar -fix-LAP-imstats (uniquer)\n"); + exit (2); + } + UNIQUER = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-fix-LAP-edges"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_FIX_LAP_EDGES; + remove_argument (N, &argc, argv); + SKIP_IMAGES = TRUE; + } + EDGE_DELETIONS = NULL; + if ((N = get_argument (argc, argv, "-fix-LAP-edges-delete"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_FIX_LAP_EDGES_DELETE; + remove_argument (N, &argc, argv); + if (N == argc) { + fprintf (stderr, "USAGE: delstar -fix-LAP-edges-delete (deletions)\n"); + exit (2); + } + EDGE_DELETIONS = strcreate (argv[N]); + remove_argument (N, &argc, argv); + SKIP_IMAGES = TRUE; } if ((N = get_argument (argc, argv, "-photcodes"))) { @@ -198,4 +253,9 @@ } usage (); + } + + IMSTATS_FILE = NULL; + if ((MODE == MODE_FIX_LAP) && (!PARALLEL)) { + IMSTATS_FILE = strcreate ("delstar.fixLAP.stats.fits"); } @@ -236,4 +296,6 @@ }*/ + if (MODE == MODE_NONE) usage (); + if (argc != 1) usage (); return (TRUE); @@ -300,4 +362,11 @@ IMAGE_ONLY = TRUE; UPDATE = FALSE; + remove_argument (N, &argc, argv); + } + + SINGLE_CPT = NULL; + if ((N = get_argument (argc, argv, "-cpt"))) { + remove_argument (N, &argc, argv); + SINGLE_CPT = strcreate (argv[N]); remove_argument (N, &argc, argv); } @@ -320,4 +389,11 @@ } + IMAGES = NULL; + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + IMAGES = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + MODE = MODE_NONE; if ((N = get_argument (argc, argv, "-photcodes"))) { @@ -333,4 +409,61 @@ remove_argument (N, &argc, argv); } + if ((N = get_argument (argc, argv, "-dup-measures"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_DUP_MEASURES; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-fix-LAP"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_FIX_LAP; + remove_argument (N, &argc, argv); + if (!IMAGES) delstar_client_usage(); + } + if ((N = get_argument (argc, argv, "-fix-LAP-edges"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_FIX_LAP_EDGES; + remove_argument (N, &argc, argv); + } + EDGE_DELETIONS = NULL; + if ((N = get_argument (argc, argv, "-fix-LAP-edges-delete"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_FIX_LAP_EDGES_DELETE; + remove_argument (N, &argc, argv); + if (N == argc) { + fprintf (stderr, "USAGE: delstar -fix-LAP-edges-delete (deletions)\n"); + exit (2); + } + EDGE_DELETIONS = strcreate (argv[N]); + remove_argument (N, &argc, argv); + SKIP_IMAGES = TRUE; + } + + MEASURE_EDGE_FILE = NULL; + if (MODE == MODE_FIX_LAP_EDGES) { + N = get_argument (argc, argv, "-measures"); + if (N == 0) { + fprintf (stderr, "delstar_client -fix-LAP-edge needs -measures (FILE)\n"); + exit (2); + } + remove_argument (N, &argc, argv); + if (N == argc) { + fprintf (stderr, "USAGE: delstar -fix-LAP-edges -measures (FILE) [FILE missing]\n"); + exit (2); + } + MEASURE_EDGE_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + IMSTATS_FILE = NULL; + if (MODE == MODE_FIX_LAP) { + N = get_argument (argc, argv, "-imstats"); + if (N == 0) { + fprintf (stderr, "delstar_client -fix-LAP needs -imstats (FILE)\n"); + exit (2); + } + remove_argument (N, &argc, argv); + IMSTATS_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } if (MODE == MODE_NONE) delstar_client_usage (); Index: /trunk/Ohana/src/delstar/src/delete_duplicate_images.c =================================================================== --- /trunk/Ohana/src/delstar/src/delete_duplicate_images.c (revision 35757) +++ /trunk/Ohana/src/delstar/src/delete_duplicate_images.c (revision 35758) @@ -6,5 +6,5 @@ // find & delete duplicate images // duplicates based on externID (skip 0s) -int delete_duplicate_images (int hostID, char *hostpath, FITS_DB *db) { +int delete_duplicate_images (FITS_DB *db) { off_t i, Nimage; @@ -28,5 +28,5 @@ if (!IMAGE_ONLY) { - delete_duplicate_image_measures (hostID, hostpath, imageID); + delete_duplicate_image_measures (imageID); } @@ -63,5 +63,5 @@ } -int delete_duplicate_image_measures (int hostID, char *hostpath, IndexArray *imageID) { +int delete_duplicate_image_measures (IndexArray *imageID) { int i; @@ -83,5 +83,5 @@ // launch the remote jobs - if (PARALLEL && !hostID) { + if (PARALLEL && !HOST_ID) { int status = delete_duplicate_image_measures_parallel (skylist, imageID); return status; @@ -92,10 +92,10 @@ // does this host ID match the desired location for the table? - if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; // set up the basic catalog info char hostfile[1024]; - snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); - catalog.filename = hostID ? hostfile : skylist[0].filename[i]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; catalog.Nsecfilt = GetPhotcodeNsecfilt (); catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; @@ -241,4 +241,5 @@ // mark the measures to be dropped for (i = 0; i < Nmeasure; i++) { + measureDrop[i] = 0; if (measure[i].imageID == 0) continue; off_t Ni = measure[i].imageID - imageID->minID; Index: /trunk/Ohana/src/delstar/src/delete_duplicate_measures.c =================================================================== --- /trunk/Ohana/src/delstar/src/delete_duplicate_measures.c (revision 35758) +++ /trunk/Ohana/src/delstar/src/delete_duplicate_measures.c (revision 35758) @@ -0,0 +1,424 @@ +# include "delstar.h" + +// this function identifies detections to be deleted as being duplicates based on imageID + detID + +int delete_duplicate_measures () { + + int i; + Catalog catalog; + + // load the current sky table (layout of all SkyRegions) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); + if (!sky) { + fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); + exit (2); + } + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area (default depth) + SkyList *skylist = NULL; + if (SINGLE_CPT) { + skylist = SkyRegionByCPT (sky, SINGLE_CPT); + } else { + skylist = SkyListByPatch (sky, -1, &UserPatch); + } + if (!skylist) { + fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); + exit (2); + } + + // launch the remote jobs + if (PARALLEL && !HOST_ID) { + int status = delete_duplicate_measures_parallel (skylist); + return status; + } + + // delete detections from region + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + // set up the basic catalog info + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + if (VERBOSE) fprintf (stderr, "deleting from %s\n", catalog.filename); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "a")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + exit (2); + } + if (!catalog.Naves_disk) { + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + delete_duplicate_measures_catalog (&catalog); + + if (UPDATE) { + // XXX skip if nothing was deleted + dvo_catalog_save_complete (&catalog, VERBOSE2); + } + + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + return TRUE; +} + +// CATDIR is supplied globally +int delete_duplicate_measures_parallel (SkyList *sky) { + + // launch the delstar_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i, j; + for (i = 0; i < table->Nhosts; i++) { + + if (sky->Nregions < table->Nhosts) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname); + continue; + } + } + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + char command[1024]; + snprintf (command, 1024, "delstar_client -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -dup-measures", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); 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; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running delstar_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} + +int delete_duplicate_measures_catalog (Catalog *catalog) { + + off_t i, j, n, m, N, D, currentAve; + + Measure *measureOut = NULL; + Average *averageOut = NULL; + SecFilt *secfiltOut = NULL; + + off_t *measureDrop, *measureSeqRaw, *measureSeqOut, *measureRefOut, *measureAveRaw, *averageNmeas, *averageDmeas, *averageSeqOut, *averageRefOut, *measureAveOut; + uint64_t *fullID; + int *seq, *duplicates; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + /* internal counters */ + off_t Nmeasure = catalog[0].Nmeasure; + off_t Naverage = catalog[0].Naverage; + + Measure *measure = catalog[0].measure; + Average *average = catalog[0].average; + SecFilt *secfilt = catalog[0].secfilt; + + if (VERBOSE) fprintf (stderr, "starting with Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT"\n", catalog[0].Naverage, catalog[0].Nmeasure); + + // we have a table of average objects and an unsorted table of measurements. each measurement + // has a reference to the average object sequence (as well as an ID) + // measure[i].averef -> average[averef] + // measure[i].objID = average[averef].objID + // measure[i].catID = average[averef].catID + + // we want a sorted measure array with all averef entries in sequence, skipping the entries which match our criteria + + // arrays + ALLOCATE (measureDrop, off_t, Nmeasure); + ALLOCATE (measureSeqRaw, off_t, Nmeasure); + ALLOCATE (measureSeqOut, off_t, Nmeasure); + ALLOCATE (measureRefOut, off_t, Nmeasure); + ALLOCATE (measureAveRaw, off_t, Nmeasure); + ALLOCATE (averageNmeas, off_t, Naverage); + ALLOCATE (averageDmeas, off_t, Naverage); + ALLOCATE (averageSeqOut, off_t, Naverage); + ALLOCATE (averageRefOut, off_t, Naverage); + + ALLOCATE (fullID, uint64_t, Nmeasure); + ALLOCATE (seq, int, Nmeasure); + ALLOCATE (duplicates, int, Nmeasure); + + memset (measureDrop, 0, sizeof(off_t) * Nmeasure); + memset (duplicates, 0, sizeof(int) * Nmeasure); + + // make an array of the full IDs + for (i = 0; i < Nmeasure; i++) { + fullID[i] = ((uint64_t) (measure[i].imageID) << 32) | (measure[i].detID); + seq[i] = i; + } + + // sort the array and sequence number + sort_fullIDs (fullID, seq, catalog[0].Nmeasure); + + uint64_t current = fullID[0]; // current unique value + int count = 0; // number of duplicates for the current unique value + duplicates[0] = count; // duplicate sequence + + for (i = 1; i < catalog[0].Nmeasure; i++) { + off_t j = seq[i]; + int newID = (measure[j].imageID == 0) || (fullID[i] != current); + if (newID) { + count = 0; + current = fullID[i]; + } else { + count ++; + } + duplicates[i] = count; + } + + // mark the measures to be dropped + for (i = 0; i < Nmeasure; i++) { + if (!duplicates[i]) continue; + off_t j = seq[i]; + measureDrop[j] = TRUE; + off_t N = measure[j].averef; + fprintf (stderr, "0x%08x 0x%08x %8.4f %8.4f %5d\n", measure[j].imageID, measure[j].detID, average[N].R, average[N].D, measure[j].photcode); + } + + // set up the measure sequence lists + for (i = 0; i < Nmeasure; i++) { + measureSeqRaw[i] = i; + measureAveRaw[i] = measure[i].averef; + } + // sort measureSeqRaw and measureAveRaw in order of measureAveRaw + SortAveMeasMatch(measureSeqRaw, measureAveRaw, Nmeasure); + + // generate the measureSeqOut array + off_t NmeasOut = 0; + for (i = 0; i < Nmeasure; i++) { + j = measureSeqRaw[i]; + if (measureDrop[j]) { + continue; + } + measureSeqOut[NmeasOut] = j; + measureRefOut[j] = NmeasOut; + NmeasOut ++; + } + // n = measureSeqOut[i] : measureOut[i] = measure[n] (i = 0 -- NmeasOut, n = 0 -- Nmeasure) + // n = measureRefOut[i] : measureOut[n] = measure[i] + ALLOCATE (measureAveOut, off_t, NmeasOut); + + // count the number of measures for each averef + N = 0; + D = -1; + currentAve = measureAveRaw[0]; + for (i = 0; i < Nmeasure; i++) { + if (measureAveRaw[i] != currentAve) { + // we have hit the next entry in the list + averageNmeas[currentAve] = N; // number of measures + averageDmeas[currentAve] = D; // first measure + N = 0; + D = -1; + currentAve = measureAveRaw[i]; + } + j = measureSeqRaw[i]; + if (measureDrop[j]) continue; + N++; + if (D == -1) { + // first valid measureOut for this averef + D = measureRefOut[j]; + } + } + averageNmeas[currentAve] = N; // number of measures + averageDmeas[currentAve] = D; // first measure + + // generate the new average sequence, skipping entries with no measurements + off_t NaveOut = 0; + for (i = 0; i < Naverage; i++) { + if (averageNmeas[i] == 0) continue; + averageSeqOut[NaveOut] = i; + averageRefOut[i] = NaveOut; + NaveOut ++; + } + // n = averageSeqOut[i] : averageOut[i] = average[n] (i = 0 -- NaveOut, n = 0 -- Naverage) + // n = averageRefOut[i] : averageOut[n] = average[i] + + // generate the output averefs + for (i = 0; i < NmeasOut; i++) { + j = measureSeqOut[i]; // measure[j] = measureOut[i] + n = measureAveRaw[j]; // average[n] : measure[j] + N = averageRefOut[n]; // averageOut[N] = average[n]; + measureAveOut[i] = N; // measureOut[i].averef = measureAveOut[i] + } + + // copy the (kept) measurements in the sorted order + ALLOCATE (measureOut, Measure, NmeasOut); + for (i = 0; i < NmeasOut; i++) { + j = measureSeqOut[i]; + measureOut[i] = measure[j]; + measureOut[i].averef = measureAveOut[i]; + } + + // copy the (kept) average entries + ALLOCATE (averageOut, Average, NaveOut); + for (i = 0; i < NaveOut; i++) { + j = averageSeqOut[i]; + averageOut[i] = average[j]; + averageOut[i].Nmeasure = averageNmeas[j]; + averageOut[i].measureOffset = averageDmeas[j]; + } + + // copy the secfilt entries for the (kept) average entries + ALLOCATE (secfiltOut, SecFilt, NaveOut*Nsecfilt); + for (i = 0; i < NaveOut; i++) { + j = averageSeqOut[i]; + for (m = 0; m < Nsecfilt; m++) { + secfiltOut[i*Nsecfilt + m] = secfilt[j*Nsecfilt + m]; + } + } + + // update the values of average.measureOffset and average.Nmeasure + FREE(measure); + FREE(average); + FREE(secfilt); + catalog[0].measure = measureOut; + catalog[0].average = averageOut; + catalog[0].secfilt = secfiltOut; + catalog[0].Nmeasure = NmeasOut; + catalog[0].Naverage = NaveOut; + + // update catalog.Nmeasure and catalog.Naverage, double check + int NmeasureTotal = 0; + int measureOffsetOK = TRUE; + int averefOK = TRUE; + for (i = 0; i < NaveOut; i++) { + NmeasureTotal += catalog[0].average[i].Nmeasure; + if (VERBOSE2 && !(NmeasureTotal <= catalog[0].Nmeasure)) { + fprintf (stderr, "too few measurements: %d %d %d\n", (int) i, NmeasureTotal, (int) catalog[0].Nmeasure); + } + measureOffsetOK &= (catalog[0].average[i].measureOffset < catalog[0].Nmeasure); + if (VERBOSE2 && !(catalog[0].average[i].measureOffset < catalog[0].Nmeasure)) { + fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure); + } + measureOffsetOK &= (catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure); + if (VERBOSE2 && !(catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure)) { + fprintf (stderr, "orrset + Nmeasure too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].measureOffset, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure); + } + m = catalog[0].average[i].measureOffset; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { + averefOK &= (catalog[0].measure[m+j].averef == i); + if (VERBOSE2 && !(catalog[0].measure[m+j].averef == i)) { + fprintf (stderr, "averef broken: %d vs %d (measure %d)\n", (int) i, catalog[0].measure[m+j].averef, (int) (m+j)); + } + } + } + + if (!measureOffsetOK) { + fprintf (stderr, "ERROR: catalog %s has an invalid measureOffset\n", catalog[0].filename); + } + + if (!averefOK) { + fprintf (stderr, "ERROR: catalog %s has invalid averefs\n", catalog[0].filename); + } + + if (NmeasureTotal != catalog[0].Nmeasure) { + fprintf (stderr, "ERROR: catalog %s has an invalid Nmeasure\n", catalog[0].filename); + } + + // MARKTIME(" match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeasure, (long long) Naverage); + + catalog[0].sorted = TRUE; + + if (VERBOSE) fprintf (stderr, "ending with Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT"\n", catalog[0].Naverage, catalog[0].Nmeasure); + + off_t NdelAves = Naverage - catalog[0].Naverage; + off_t NdelMeas = Nmeasure - catalog[0].Nmeasure; + if (NdelAves || NdelMeas) { + fprintf (stderr, "deleting "OFF_T_FMT" measures and "OFF_T_FMT" averages : %s\n", NdelMeas, NdelAves, catalog[0].filename); + } + + FREE (fullID); + FREE (seq); + FREE (duplicates); + + FREE (measureDrop); + FREE (measureSeqRaw); + FREE (measureSeqOut); + FREE (measureRefOut); + FREE (measureAveRaw); + FREE (averageNmeas); + FREE (averageDmeas); + FREE (averageSeqOut); + FREE (averageRefOut); + FREE (measureAveOut); + + return TRUE; +} + +// sort fullId vector and an index vector +void sort_fullIDs (uint64_t *I, int *S, off_t N) { + +# define SWAPFUNC(A,B){ uint64_t tmp_i; int tmp_s; \ + tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \ + tmp_s = S[A]; S[A] = S[B]; S[B] = tmp_s; \ +} +# define COMPARE(A,B)(I[A] < I[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + Index: /trunk/Ohana/src/delstar/src/delete_fix_LAP.c =================================================================== --- /trunk/Ohana/src/delstar/src/delete_fix_LAP.c (revision 35758) +++ /trunk/Ohana/src/delstar/src/delete_fix_LAP.c (revision 35758) @@ -0,0 +1,667 @@ +# include "delstar.h" + +// this function identifies detections to be deleted as being duplicates based on imageID + detID + +int delete_fix_LAP (ImageSubset *image, off_t Nimage) { + + int i; + Catalog catalog; + + // load the current sky table (layout of all SkyRegions) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); + if (!sky) { + fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); + exit (2); + } + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area (default depth) + SkyList *skylist = SINGLE_CPT ? SkyRegionByCPT (sky, SINGLE_CPT) : SkyListByPatch (sky, -1, &UserPatch); + if (!skylist) { + fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); + exit (2); + } + + // this is needed in the master and client jobs (in master, just so Nimage is set in ImageOpsFixLAP.c) + initImageIndex (image, Nimage); + + // launch the remote jobs + if (PARALLEL && !HOST_ID) { + int status = delete_fix_LAP_parallel (skylist, image, Nimage); + return status; + } + + // delete detections from region + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + // set up the basic catalog info + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + if (VERBOSE) fprintf (stderr, "deleting from %s\n", catalog.filename); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "a")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + exit (2); + } + if (!catalog.Naves_disk) { + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + off_t NaverageStart = catalog.Naverage; + off_t NmeasureStart = catalog.Nmeasure; + delete_fix_LAP_catalog (&catalog, image, Nimage); + + // skip if nothing was deleted + int noChange = (NaverageStart == catalog.Naverage) && (NmeasureStart == catalog.Nmeasure); + if (UPDATE && !noChange) { + dvo_catalog_save_complete (&catalog, VERBOSE2); + } + + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + + ImageValidSave (IMSTATS_FILE); + + return TRUE; +} + +// CATDIR is supplied globally +int delete_fix_LAP_parallel (SkyList *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/delstar.fixLAP.%s.dat", CATDIR, uniquer); + + if (!ImageSubsetSave (imageFile, image, Nimage)) { + fprintf (stderr, "failed to write image subset\n"); + exit (1); + } + + // launch the delstar_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i, j; + for (i = 0; i < table->Nhosts; i++) { + + if (sky->Nregions < table->Nhosts) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname); + continue; + } + } + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/delstar.fixLAP.imstats.%s.dat", table->hosts[i].pathname, uniquer); + + char command[1024]; + snprintf (command, 1024, "delstar_client -hostID %d -D CATDIR %s -hostdir %s -images %s -imstats %s -region %f %f %f %f -fix-LAP", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, imageFile, table->hosts[i].results, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); 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; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running delstar_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].results) continue; + if (!ImageValidLoad (table->hosts[i].results)) { + fprintf (stderr, "failed to read data from %s\n", table->hosts[i].hostname); + } + } + + return TRUE; +} + +// CATDIR is supplied globally +int delete_fix_LAP_setstats (ImageSubset *image, off_t Nimage) { + + // load the current sky table (layout of all SkyRegions) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); + if (!sky) { + fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); + exit (2); + } + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area (default depth) + SkyList *skylist = SINGLE_CPT ? SkyRegionByCPT (sky, SINGLE_CPT) : SkyListByPatch (sky, -1, &UserPatch); + if (!skylist) { + fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); + exit (2); + } + fprintf (stderr, "setup sky\n"); + + // this is needed in the master and client jobs (in master, just so Nimage is set in ImageOpsFixLAP.c) + initImageIndex (image, Nimage); + fprintf (stderr, "make index\n"); + + // write out the subset table of image information + char imageFile[512]; + snprintf (imageFile, 512, "%s/delstar.fixLAP.%s.dat", CATDIR, UNIQUER); + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/delstar.fixLAP.imstats.%s.dat", table->hosts[i].pathname, UNIQUER); + + fprintf (stderr, "read %s\n", table->hosts[i].results); + if (!ImageValidLoad (table->hosts[i].results)) { + fprintf (stderr, "failed to read data from %s\n", table->hosts[i].hostname); + } + } + + return TRUE; +} + +int UnpackPSPSStackDetectionID(int *sourceID, int *imageID, int *detID, uint64_t pspsStackID); + +int Nvalid = 0; +int Ninvalid = 0; + +int delete_fix_LAP_catalog (Catalog *catalog, ImageSubset *image, off_t Nimage) { + + off_t i, j, n, m, N, D, currentAve; + + Measure *measureOut = NULL; + Average *averageOut = NULL; + SecFilt *secfiltOut = NULL; + + off_t *measureDrop, *measureSeqRaw, *measureSeqOut, *measureRefOut, *measureAveRaw, *averageNmeas, *averageDmeas, *averageSeqOut, *averageRefOut, *measureAveOut; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + /* internal counters */ + off_t Nmeasure = catalog[0].Nmeasure; + off_t Naverage = catalog[0].Naverage; + + Measure *measure = catalog[0].measure; + Average *average = catalog[0].average; + SecFilt *secfilt = catalog[0].secfilt; + + Nvalid = 0; + Ninvalid = 0; + if (VERBOSE) fprintf (stderr, "starting with Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT"\n", catalog[0].Naverage, catalog[0].Nmeasure); + + // we have a table of average objects and an unsorted table of measurements. each measurement + // has a reference to the average object sequence (as well as an ID) + // measure[i].averef -> average[averef] + // measure[i].objID = average[averef].objID + // measure[i].catID = average[averef].catID + + // we want a sorted measure array with all averef entries in sequence, skipping the entries which match our criteria + + // arrays + ALLOCATE (measureDrop, off_t, Nmeasure); + delete_fix_LAP_measures (measureDrop, catalog, image, Nimage); + + ALLOCATE (measureSeqRaw, off_t, Nmeasure); + ALLOCATE (measureSeqOut, off_t, Nmeasure); + ALLOCATE (measureRefOut, off_t, Nmeasure); + ALLOCATE (measureAveRaw, off_t, Nmeasure); + ALLOCATE (averageNmeas, off_t, Naverage); + ALLOCATE (averageDmeas, off_t, Naverage); + ALLOCATE (averageSeqOut, off_t, Naverage); + ALLOCATE (averageRefOut, off_t, Naverage); + + // set up the measure sequence lists + for (i = 0; i < Nmeasure; i++) { + measureSeqRaw[i] = i; + measureAveRaw[i] = measure[i].averef; + } + // sort measureSeqRaw and measureAveRaw in order of measureAveRaw + SortAveMeasMatch(measureSeqRaw, measureAveRaw, Nmeasure); + + // generate the measureSeqOut array + off_t NmeasOut = 0; + for (i = 0; i < Nmeasure; i++) { + j = measureSeqRaw[i]; + if (measureDrop[j]) { + continue; + } + measureSeqOut[NmeasOut] = j; + measureRefOut[j] = NmeasOut; + NmeasOut ++; + } + // n = measureSeqOut[i] : measureOut[i] = measure[n] (i = 0 -- NmeasOut, n = 0 -- Nmeasure) + // n = measureRefOut[i] : measureOut[n] = measure[i] + ALLOCATE (measureAveOut, off_t, NmeasOut); + + // count the number of measures for each averef + N = 0; + D = -1; + currentAve = measureAveRaw[0]; + for (i = 0; i < Nmeasure; i++) { + if (measureAveRaw[i] != currentAve) { + // we have hit the next entry in the list + averageNmeas[currentAve] = N; // number of measures + averageDmeas[currentAve] = D; // first measure + N = 0; + D = -1; + currentAve = measureAveRaw[i]; + } + j = measureSeqRaw[i]; + if (measureDrop[j]) continue; + N++; + if (D == -1) { + // first valid measureOut for this averef + D = measureRefOut[j]; + } + } + averageNmeas[currentAve] = N; // number of measures + averageDmeas[currentAve] = D; // first measure + + // generate the new average sequence, skipping entries with no measurements + off_t NaveOut = 0; + for (i = 0; i < Naverage; i++) { + if (averageNmeas[i] == 0) continue; + averageSeqOut[NaveOut] = i; + averageRefOut[i] = NaveOut; + NaveOut ++; + } + // n = averageSeqOut[i] : averageOut[i] = average[n] (i = 0 -- NaveOut, n = 0 -- Naverage) + // n = averageRefOut[i] : averageOut[n] = average[i] + + // generate the output averefs + for (i = 0; i < NmeasOut; i++) { + j = measureSeqOut[i]; // measure[j] = measureOut[i] + n = measureAveRaw[j]; // average[n] : measure[j] + N = averageRefOut[n]; // averageOut[N] = average[n]; + measureAveOut[i] = N; // measureOut[i].averef = measureAveOut[i] + } + + // copy the (kept) measurements in the sorted order + ALLOCATE (measureOut, Measure, NmeasOut); + for (i = 0; i < NmeasOut; i++) { + j = measureSeqOut[i]; + measureOut[i] = measure[j]; + measureOut[i].averef = measureAveOut[i]; + } + + // copy the (kept) average entries + ALLOCATE (averageOut, Average, NaveOut); + for (i = 0; i < NaveOut; i++) { + j = averageSeqOut[i]; + averageOut[i] = average[j]; + averageOut[i].Nmeasure = averageNmeas[j]; + averageOut[i].measureOffset = averageDmeas[j]; + } + + // copy the secfilt entries for the (kept) average entries + ALLOCATE (secfiltOut, SecFilt, NaveOut*Nsecfilt); + for (i = 0; i < NaveOut; i++) { + j = averageSeqOut[i]; + for (m = 0; m < Nsecfilt; m++) { + secfiltOut[i*Nsecfilt + m] = secfilt[j*Nsecfilt + m]; + } + } + + // update the values of average.measureOffset and average.Nmeasure + FREE(measure); + FREE(average); + FREE(secfilt); + catalog[0].measure = measureOut; + catalog[0].average = averageOut; + catalog[0].secfilt = secfiltOut; + catalog[0].Nmeasure = NmeasOut; + catalog[0].Naverage = NaveOut; + + // update catalog.Nmeasure and catalog.Naverage, double check + int NmeasureTotal = 0; + int measureOffsetOK = TRUE; + int averefOK = TRUE; + for (i = 0; i < NaveOut; i++) { + NmeasureTotal += catalog[0].average[i].Nmeasure; + if (VERBOSE2 && !(NmeasureTotal <= catalog[0].Nmeasure)) { + fprintf (stderr, "too few measurements: %d %d %d\n", (int) i, NmeasureTotal, (int) catalog[0].Nmeasure); + } + measureOffsetOK &= (catalog[0].average[i].measureOffset < catalog[0].Nmeasure); + if (VERBOSE2 && !(catalog[0].average[i].measureOffset < catalog[0].Nmeasure)) { + fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure); + } + measureOffsetOK &= (catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure); + if (VERBOSE2 && !(catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure)) { + fprintf (stderr, "offset + Nmeasure too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].measureOffset, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure); + } + m = catalog[0].average[i].measureOffset; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { + averefOK &= (catalog[0].measure[m+j].averef == i); + if (VERBOSE2 && !(catalog[0].measure[m+j].averef == i)) { + fprintf (stderr, "averef broken: %d vs %d (measure %d)\n", (int) i, catalog[0].measure[m+j].averef, (int) (m+j)); + } + } + } + + if (!measureOffsetOK) { + fprintf (stderr, "ERROR: catalog %s has an invalid measureOffset\n", catalog[0].filename); + } + + if (!averefOK) { + fprintf (stderr, "ERROR: catalog %s has invalid averefs\n", catalog[0].filename); + } + + if (NmeasureTotal != catalog[0].Nmeasure) { + fprintf (stderr, "ERROR: catalog %s has an invalid Nmeasure\n", catalog[0].filename); + } + + // MARKTIME(" match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeasure, (long long) Naverage); + + catalog[0].sorted = TRUE; + + if (VERBOSE) fprintf (stderr, "ending with Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT"\n", catalog[0].Naverage, catalog[0].Nmeasure); + + off_t NdelAves = Naverage - catalog[0].Naverage; + off_t NdelMeas = Nmeasure - catalog[0].Nmeasure; + if (NdelAves || NdelMeas) { + fprintf (stderr, "deleting "OFF_T_FMT" measures and "OFF_T_FMT" averages (%d valid, %d invalid): %s\n", NdelMeas, NdelAves, Nvalid, Ninvalid, catalog[0].filename); + } + + FREE (measureDrop); + FREE (measureSeqRaw); + FREE (measureSeqOut); + FREE (measureRefOut); + FREE (measureAveRaw); + FREE (averageNmeas); + FREE (averageDmeas); + FREE (averageSeqOut); + FREE (averageRefOut); + FREE (measureAveOut); + + return TRUE; +} + +// this function makes the needed modifications to the measure table (to imageIDs) and updates +// the array measureDrop to mark the entries we want to remove + +// needed modifications: +// * delete all photcodes corresponding to 2mass, synth, or wise (not gpc1 or stack) +// * delete all duplicates +// * check the imageID +// * re-assign the imageID if needed + +int delete_fix_LAP_measures (off_t *measureDrop, Catalog *catalog, ImageSubset *image, off_t Nimage) { + + /* internal counters */ + off_t Nmeasure = catalog[0].Nmeasure; + + Measure *measure = catalog[0].measure; + // Average *average = catalog[0].average; + + uint64_t *fullID; + int *seq, *duplicates; + + ALLOCATE (fullID, uint64_t, Nmeasure); + ALLOCATE (seq, int, Nmeasure); + ALLOCATE (duplicates, int, Nmeasure); + + memset (measureDrop, 0, sizeof(off_t) * Nmeasure); + memset (duplicates, 0, sizeof(int) * Nmeasure); + + // make an array of the full IDs + off_t i; + for (i = 0; i < Nmeasure; i++) { + fullID[i] = ((uint64_t) (measure[i].imageID) << 32) | (measure[i].detID); + seq[i] = i; + } + + // sort the array and sequence number + sort_fullIDs (fullID, seq, catalog[0].Nmeasure); + + uint64_t current = fullID[0]; // current unique value + int count = 0; // number of duplicates for the current unique value + duplicates[0] = count; // duplicate sequence + + for (i = 1; i < catalog[0].Nmeasure; i++) { + off_t j = seq[i]; + int newID = (measure[j].imageID == 0) || (fullID[i] != current); + if (newID) { + count = 0; + current = fullID[i]; + } else { + count ++; + } + duplicates[i] = count; + } + + // mark the measures to be dropped + for (i = 0; i < Nmeasure; i++) { + off_t j = seq[i]; + + int keep = TRUE; + if (duplicates[i]) keep = FALSE; + if (measure[j].photcode < 10000) keep = FALSE; // Drop 2MASS, WISE, SYNTH here + + if (!keep) { + measureDrop[j] = TRUE; + // off_t N = measure[j].averef; + // fprintf (stderr, "0x%08x 0x%08x %8.4f %8.4f %5d\n", measure[j].imageID, measure[j].detID, average[N].R, average[N].D, measure[j].photcode); + } + } + + int maxID = 0; + int *imageIndex = getImageIndex(&maxID); + + // check the image IDs for the remaining measures + for (i = 0; i < Nmeasure; i++) { + off_t j = seq[i]; + if (measureDrop[j]) continue; + + // for gpc1 exposures, 10000 < photcode < 10600 + if ((10000 < measure[j].photcode) && (measure[j].photcode < 10600)) { + int resetID = FALSE; + off_t imageID = measure[j].imageID; + if (imageID == 0) { + resetID = TRUE; + goto resetID_exp; + } + if (imageID > maxID) { + fprintf (stderr, "invalid image ID "OFF_T_FMT" for catalog %s\n", imageID, catalog[0].filename); + exit (2); + } + off_t imageN = imageIndex[imageID]; + if (imageN < 0) { + // image ID for this measure is not in the image table + resetID = TRUE; + } + if ((measure[j].t < image[imageN].tmin) || (measure[j].t > image[imageN].tmax)) { + // measure.imageID does not match correct image + resetID = TRUE; + } + if (measure[j].photcode != image[imageN].photcode) { + // measure.imageID does not match correct image + resetID = TRUE; + } + + if (!resetID) { + // off_t N = measure[j].averef; + // fprintf (stderr, "no repair 0x%08x 0x%08x %8.4f %8.4f %5d\n", measure[j].imageID, measure[j].detID, average[N].R, average[N].D, measure[j].photcode); + BumpValidImage (imageN); + Nvalid ++; + continue; + } + + resetID_exp: + if (resetID) { + // off_t N = measure[j].averef; + // fprintf (stderr, "repair 0x%08x 0x%08x %8.4f %8.4f %5d\n", measure[j].imageID, measure[j].detID, average[N].R, average[N].D, measure[j].photcode); + // update the image ID + short photcode = measure[j].photcode; + e_time time = measure[j].t; + + off_t imageID, Seq; + if (!FindIDexp (&imageID, &Seq, time, photcode)) { + fprintf (stderr, "error cannot find matching image ID\n"); + exit (5); + } + if (imageID <= 0) { + fprintf (stderr, "error cannot find matching image ID\n"); + exit (5); + } + measure[j].imageID = imageID; + BumpInvalidImage (Seq); + Ninvalid ++; + continue; + } + } + + // for gpc1 exposures, 10000 < photcode < 10600 + if ((11000 <= measure[j].photcode) && (measure[j].photcode < 11600)) { + int resetID = FALSE; + off_t imageID = measure[j].imageID; + if (imageID == 0) { + resetID = TRUE; + goto resetID_stk; + } + if (imageID > maxID) { + fprintf (stderr, "invalid image ID "OFF_T_FMT" for catalog %s\n", imageID, catalog[0].filename); + exit (2); + } + + off_t imageN = imageIndex[imageID]; + if (imageN < 0) { + // image ID for this measure is not in the image table + resetID = TRUE; + } + int sourceID, externID, detID; + UnpackPSPSStackDetectionID (&sourceID, &externID, &detID, measure[j].extID); + if (externID != image[imageN].externID) { + // measure.imageID does not match correct image + resetID = TRUE; + } + if (measure[j].photcode != image[imageN].photcode) { + // measure.imageID does not match correct image + resetID = TRUE; + } + + if (!resetID) { + // off_t N = measure[j].averef; + // fprintf (stderr, "no repair 0x%08x 0x%08x %8.4f %8.4f %5d\n", measure[j].imageID, measure[j].detID, average[N].R, average[N].D, measure[j].photcode); + BumpValidImage (imageN); + Nvalid ++; + continue; + } + + resetID_stk: + if (resetID) { + // off_t N = measure[j].averef; + // fprintf (stderr, "repair 0x%08x 0x%08x %8.4f %8.4f %5d\n", measure[j].imageID, measure[j].detID, average[N].R, average[N].D, measure[j].photcode); + // update the image ID + off_t imageID, Seq; + short Photcode; + if (!FindIDstk (&imageID, &Seq, &Photcode, externID)) { + fprintf (stderr, "error : cannot find matching stack image ID\n"); + exit (6); + } + measure[j].imageID = imageID; + BumpInvalidImage (Seq); + Ninvalid ++; + continue; + } + } + } + + FREE (fullID); + FREE (seq); + FREE (duplicates); + + return TRUE; +} + +int UnpackPSPSStackDetectionID(int *sourceID, int *imageID, int *detID, uint64_t pspsStackID) +{ + // sourceID : ID of database + table that tracked the image (< 0x7f = 127) + // imageID : external ID of the image which provided the detections (< 0x1000.0000 ~ 2.7e8) + // detID : detection sequence in image (< 0x1000.0000 ~ 2.7e8) + + // 0x0000.0000.0000.0000 + + // bits 0 - 27 : 0x0000.0000.0fff.ffff + // bits 28 - 55 : 0x00ff.ffff.f000.0000 + // bits 56 - 63 : 0xff00.0000.0000.0000 + + *sourceID = (pspsStackID & (uint64_t) 0xff00000000000000) >> 56; + *imageID = (pspsStackID & (uint64_t) 0x00fffffff0000000) >> 28; + *detID = (pspsStackID & (uint64_t) 0x000000000fffffff); + + return TRUE; +} Index: /trunk/Ohana/src/delstar/src/delete_fix_LAP_edges.c =================================================================== --- /trunk/Ohana/src/delstar/src/delete_fix_LAP_edges.c (revision 35758) +++ /trunk/Ohana/src/delstar/src/delete_fix_LAP_edges.c (revision 35758) @@ -0,0 +1,315 @@ +# include "delstar.h" + +// this function identifies detections to be deleted as being duplicates based on imageID +// + detID *** at the edges of catalogs *** + +// unlike many other versions of dvo operations, this function needs to gather detections +// back to the master program and then do the analysis (more like relphot). the first +// block of the code slurps the data from the catalogs, the second block looks for duplicates + +MeasureEdge *delete_fix_LAP_edges (off_t *nmeasure_edge) { + + int i; + Catalog catalog; + + // load the current sky table (layout of all SkyRegions) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); + if (!sky) { + fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); + exit (2); + } + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area (default depth) + SkyList *skylist = SINGLE_CPT ? SkyRegionByCPT (sky, SINGLE_CPT) : SkyListByPatch (sky, -1, &UserPatch); + if (!skylist) { + fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); + exit (2); + } + + // launch the remote jobs + if (PARALLEL && !HOST_ID) { + MeasureEdge *measure_edge = delete_fix_LAP_edges_parallel (skylist, nmeasure_edge); + return measure_edge; + } + + // XXX accumulate detections in a local structure: + off_t Nmeasure_edge = 0; + off_t NMEASURE_EDGE = 1000; + MeasureEdge *measure_edge = NULL; + ALLOCATE (measure_edge, MeasureEdge, NMEASURE_EDGE); + + // delete detections from region + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + // set up the basic catalog info + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + if (VERBOSE) fprintf (stderr, "deleting from %s\n", catalog.filename); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "a")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + exit (2); + } + if (!catalog.Naves_disk) { + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + // off_t NaverageStart = catalog.Naverage; + // off_t NmeasureStart = catalog.Nmeasure; + measure_edge = delete_fix_LAP_edges_get_measures (&catalog, measure_edge, &Nmeasure_edge, &NMEASURE_EDGE); + + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + + *nmeasure_edge = Nmeasure_edge; + return measure_edge; +} + +// CATDIR is supplied globally +MeasureEdge *delete_fix_LAP_edges_parallel (SkyList *sky, off_t *nmeasure_edge) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + // launch the delstar_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i, j; + for (i = 0; i < table->Nhosts; i++) { + + if (sky->Nregions < table->Nhosts) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname); + continue; + } + } + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/delstar.fixLAPedges.measures.%s.dat", table->hosts[i].pathname, uniquer); + + char command[1024]; + snprintf (command, 1024, "delstar_client -hostID %d -D CATDIR %s -hostdir %s -measures %s -region %f %f %f %f -fix-LAP-edges", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, table->hosts[i].results, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); 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; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running delstar_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + off_t Nmeasure_edge_all = 0; + off_t NMEASURE_EDGE_ALL = 1000; + MeasureEdge *measure_edge_all = NULL; + ALLOCATE (measure_edge_all, MeasureEdge, NMEASURE_EDGE_ALL); + + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].results) continue; + off_t Nmeasure_edge = 0; + MeasureEdge *measure_edge = MeasureEdgeLoad (table->hosts[i].results, &Nmeasure_edge); + if (!measure_edge) { + fprintf (stderr, "failed to read data from %s\n", table->hosts[i].hostname); + continue; + } + // merge new measures into a single array + measure_edge_all = MeasureEdgeMerge (measure_edge_all, &Nmeasure_edge_all, &NMEASURE_EDGE_ALL, measure_edge, Nmeasure_edge); + } + + *nmeasure_edge = Nmeasure_edge_all; + return measure_edge_all; +} + +MeasureEdge *delete_fix_LAP_edges_get_measures (Catalog *catalog, MeasureEdge *measure_edge, off_t *Nmeasure_edge_in, off_t *NMEASURE_EDGE_IN) { + + off_t i, j, m; + + off_t Nmeasure_edge = *Nmeasure_edge_in; + off_t NMEASURE_EDGE = *NMEASURE_EDGE_IN; + + // boundaries of the catalog file + double Rmin, Rmax, Dmin, Dmax; + gfits_scan (&catalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&catalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax); + +# define GAP 2.0 + double dD = GAP / 3600.0; + double dR = dD / cos(RAD_DEG * 0.5 * (Dmin + Dmax)); + + double Rgapmin = Rmin + dR; + double Rgapmax = Rmax - dR; + double Dgapmin = Dmin + dD; + double Dgapmax = Dmax - dD; + + // XXX should deal with pole, but not yet... + for (i = 0; i < catalog[0].Naverage; i++) { + if (catalog[0].average[i].R < Rgapmin) goto save_measures; + if (catalog[0].average[i].R > Rgapmax) goto save_measures; + if (catalog[0].average[i].D < Dgapmin) goto save_measures; + if (catalog[0].average[i].D > Dgapmax) goto save_measures; + continue; + + save_measures: + m = catalog[0].average[i].measureOffset; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { + measure_edge[Nmeasure_edge].R = catalog[0].average[i].R; + measure_edge[Nmeasure_edge].D = catalog[0].average[i].D; + measure_edge[Nmeasure_edge].objID = catalog[0].measure[m+j].objID; + measure_edge[Nmeasure_edge].catID = catalog[0].measure[m+j].catID; + measure_edge[Nmeasure_edge].detID = catalog[0].measure[m+j].detID; + measure_edge[Nmeasure_edge].imageID = catalog[0].measure[m+j].imageID; + } + Nmeasure_edge ++; + CHECK_REALLOCATE (measure_edge, MeasureEdge, NMEASURE_EDGE, Nmeasure_edge, 1000); + } + + *Nmeasure_edge_in = Nmeasure_edge; + *NMEASURE_EDGE_IN = NMEASURE_EDGE; + + return measure_edge; +} + +// this function makes the needed modifications to the measure table (to imageIDs) and updates +// the array measureDrop to mark the entries we want to remove + +// needed modifications: +// * delete all photcodes corresponding to 2mass, synth, or wise (not gpc1 or stack) +// * delete all duplicates +// * check the imageID +// * re-assign the imageID if needed + +int delete_fix_LAP_edges_find_dups (MeasureEdge *measure_edge, off_t Nmeasure_edge) { + + off_t i; + + // XXX TEST +# if (0) + FILE *f = fopen ("delstar.fixLAPedge.dat", "w"); + for (i = 0; i < Nmeasure_edge; i++) { + fprintf (f, "0x%08x 0x%08x | 0x%08x 0x%08x | %8.4f %8.4f\n", measure_edge[i].imageID, measure_edge[i].detID, measure_edge[i].objID, measure_edge[i].catID, measure_edge[i].R, measure_edge[i].D); + } + fclose (f); +# endif + + uint64_t *fullID; + int *seq, *duplicates; + + ALLOCATE (fullID, uint64_t, Nmeasure_edge); + ALLOCATE (seq, int, Nmeasure_edge); + ALLOCATE (duplicates, int, Nmeasure_edge); + + // memset (measureDrop, 0, sizeof(off_t) * Nmeasure_edge); + memset (duplicates, 0, sizeof(int) * Nmeasure_edge); + + // make an array of the full IDs + for (i = 0; i < Nmeasure_edge; i++) { + fullID[i] = ((uint64_t) (measure_edge[i].imageID) << 32) | (measure_edge[i].detID); + seq[i] = i; + } + + // sort the array and sequence number + sort_fullIDs (fullID, seq, Nmeasure_edge); + + uint64_t current = fullID[0]; // current unique value + int count = 0; // number of duplicates for the current unique value + duplicates[0] = count; // duplicate sequence + + for (i = 1; i < Nmeasure_edge; i++) { + off_t j = seq[i]; + int newID = (measure_edge[j].imageID == 0) || (fullID[i] != current); + if (newID) { + count = 0; + current = fullID[i]; + } else { + count ++; + } + duplicates[i] = count; + } + + FILE *f = fopen ("delstar.fixLAPedge.dat", "w"); + + // mark the measures to be dropped + for (i = 0; i < Nmeasure_edge; i++) { + off_t j = seq[i]; + + int keep = TRUE; + if (duplicates[i]) keep = FALSE; + + if (!keep) { + // measureDrop[j] = TRUE; + fprintf (f, "0x%08x 0x%08x | 0x%08x 0x%08x | %8.4f %8.4f\n", measure_edge[j].imageID, measure_edge[j].detID, measure_edge[j].objID, measure_edge[j].catID, measure_edge[j].R, measure_edge[j].D); + } + } + fclose (f); + + FREE (fullID); + FREE (seq); + FREE (duplicates); + + return TRUE; +} Index: /trunk/Ohana/src/delstar/src/delete_fix_LAP_edges_delete.c =================================================================== --- /trunk/Ohana/src/delstar/src/delete_fix_LAP_edges_delete.c (revision 35758) +++ /trunk/Ohana/src/delstar/src/delete_fix_LAP_edges_delete.c (revision 35758) @@ -0,0 +1,484 @@ +# include "delstar.h" + +// this function deletes identified detections + +int delete_fix_LAP_edges_delete () { + + int i; + Catalog catalog; + + // load the current sky table (layout of all SkyRegions) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); + if (!sky) { + fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); + exit (2); + } + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area (default depth) + SkyList *skylist = SINGLE_CPT ? SkyRegionByCPT (sky, SINGLE_CPT) : SkyListByPatch (sky, -1, &UserPatch); + if (!skylist) { + fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); + exit (2); + } + + // launch the remote jobs + if (PARALLEL && !HOST_ID) { + int status = delete_fix_LAP_edges_delete_parallel (skylist); + return status; + } + + // load the detections to delete + FILE *f = fopen (EDGE_DELETIONS, "r"); + if (!f) { + fprintf (stderr, "ERROR: failed to read the list of detections to be deleted\n"); + exit (1); + } + + // XXX accumulate detections in a local structure: + off_t Nmeasure_edge = 0; + off_t NMEASURE_EDGE = 1000; + MeasureEdge *measure_edge = NULL; + ALLOCATE (measure_edge, MeasureEdge, NMEASURE_EDGE); + char line[1024]; + while (scan_line(f, line) != EOF) { + double ra, dec; + int imageID, detID, objID, catID; + iparse (&imageID, 1, line); + iparse (&detID, 2, line); + iparse (&objID, 4, line); + iparse (&catID, 5, line); + dparse (&ra, 7, line); + dparse (&dec, 8, line); + measure_edge[Nmeasure_edge].imageID = imageID; + measure_edge[Nmeasure_edge].detID = detID; + measure_edge[Nmeasure_edge].objID = objID; + measure_edge[Nmeasure_edge].catID = catID; + measure_edge[Nmeasure_edge].R = ra; + measure_edge[Nmeasure_edge].D = dec; + Nmeasure_edge ++; + CHECK_REALLOCATE (measure_edge, MeasureEdge, NMEASURE_EDGE, Nmeasure_edge, 1000); + } + // sort_measure_edge_by_catID (measure_edge, Nmeasure_edge); + + // create a catID count table + int maxCatID = 0; + for (i = 0; i < Nmeasure_edge; i++) { + maxCatID = MAX(maxCatID, measure_edge[i].catID); + } + + // catIDcount[N] = count(catID == N) + int *catIDcount = NULL; + ALLOCATE (catIDcount, int, maxCatID + 1); + memset (catIDcount, 0, (maxCatID + 1) * sizeof(int)); + for (i = 0; i < Nmeasure_edge; i++) { + catIDcount[measure_edge[i].catID] ++; + } + + // delete detections from region + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + int myCatID = skylist[0].regions[i][0].index; + if (myCatID > maxCatID) continue; + if (!catIDcount[myCatID]) continue; + + // set up the basic catalog info + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + if (VERBOSE) fprintf (stderr, "deleting from %s\n", catalog.filename); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "a")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + exit (2); + } + if (!catalog.Naves_disk) { + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + off_t NaverageStart = catalog.Naverage; + off_t NmeasureStart = catalog.Nmeasure; + delete_fix_LAP_edges_drop_measures (&catalog, measure_edge, Nmeasure_edge, catIDcount[myCatID]); + + // skip if nothing was deleted + int noChange = (NaverageStart == catalog.Naverage) && (NmeasureStart == catalog.Nmeasure); + if (UPDATE && !noChange) { + // XXX save a backup copy first + dvo_catalog_backup (&catalog, TRUE); + dvo_catalog_save_complete (&catalog, VERBOSE2); + } + + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + free (measure_edge); + free (catIDcount); + + return TRUE; +} + +// CATDIR is supplied globally +int delete_fix_LAP_edges_delete_parallel (SkyList *sky) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + // launch the delstar_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + { + // ensure that the paths are absolute path names + char *tmppath = abspath (EDGE_DELETIONS, DVO_MAX_PATH); + free (EDGE_DELETIONS); + EDGE_DELETIONS = tmppath; + } + + int i, j; + for (i = 0; i < table->Nhosts; i++) { + + if (sky->Nregions < table->Nhosts) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname); + continue; + } + } + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + // ALLOCATE (table->hosts[i].results, char, 1024); + // snprintf (table->hosts[i].results, 1024, "%s/delstar.fixLAPedges.measures.%s.dat", table->hosts[i].pathname, uniquer); + + char command[1024]; + snprintf (command, 1024, "delstar_client -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -fix-LAP-edges-delete %s", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, + EDGE_DELETIONS); + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); 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; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running delstar_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} + +int delete_fix_LAP_edges_drop_measures (Catalog *catalog, MeasureEdge *measure_edge, off_t Nmeasure_edge, int catIDcount) { + + off_t i, j, n, m, N, D, currentAve; + + myAssert (catIDcount > 0, "oops : nothing to delete?"); + + // get a list of detections values corresponding to this catID. + // there should be catIDcount of them + + int myCatID = catalog[0].catID; + + int NbadMeasures = 0; + int *badMeasures = NULL; + ALLOCATE (badMeasures, int, catIDcount); + + for (i = 0; i < Nmeasure_edge; i++) { + if (measure_edge[i].catID != myCatID) continue; + myAssert (NbadMeasures < catIDcount, "too many bad measures?"); + badMeasures[NbadMeasures] = i; + NbadMeasures ++; + } + + Measure *measureOut = NULL; + Average *averageOut = NULL; + SecFilt *secfiltOut = NULL; + + off_t *measureDrop, *measureSeqRaw, *measureSeqOut, *measureRefOut, *measureAveRaw, *averageNmeas, *averageDmeas, *averageSeqOut, *averageRefOut, *measureAveOut; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + /* internal counters */ + off_t Nmeasure = catalog[0].Nmeasure; + off_t Naverage = catalog[0].Naverage; + + Measure *measure = catalog[0].measure; + Average *average = catalog[0].average; + SecFilt *secfilt = catalog[0].secfilt; + + if (VERBOSE) fprintf (stderr, "starting with Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT"\n", catalog[0].Naverage, catalog[0].Nmeasure); + + // we have a table of average objects and an unsorted table of measurements. each measurement + // has a reference to the average object sequence (as well as an ID) + // measure[i].averef -> average[averef] + // measure[i].objID = average[averef].objID + // measure[i].catID = average[averef].catID + + // we want a sorted measure array with all averef entries in sequence, skipping the entries which match our criteria + + // arrays + ALLOCATE (measureDrop, off_t, Nmeasure); + memset (measureDrop, 0, sizeof(off_t) * Nmeasure); + + int invalid = FALSE; + int Ndrop = 0; + for (i = 0; i < Nmeasure; i++) { + for (j = 0; j < NbadMeasures; j++) { + int seq = badMeasures[j]; + if (measure[i].imageID != measure_edge[seq].imageID) continue; + if (measure[i].detID != measure_edge[seq].detID) continue; + // this should be a measure to delete + if (measure[i].objID != measure_edge[seq].objID) { + fprintf (stderr, "ERROR: this is NOT the valid match!\n"); + invalid = TRUE; + continue; + } + if (measure[i].catID != measure_edge[seq].catID) { + fprintf (stderr, "ERROR: this is NOT the valid match!\n"); + invalid = TRUE; + continue; + } + measureDrop[i] = TRUE; + Ndrop ++; + } + } + + if (Ndrop != NbadMeasures) { + fprintf (stderr, "ERROR: wrong number of bad detections!\n"); + invalid = TRUE; + } + + if (invalid) { + fprintf (stderr, "ERROR: bad detection selection seems wrong, aborting\n"); + exit (2); + } + + ALLOCATE (measureSeqRaw, off_t, Nmeasure); + ALLOCATE (measureSeqOut, off_t, Nmeasure); + ALLOCATE (measureRefOut, off_t, Nmeasure); + ALLOCATE (measureAveRaw, off_t, Nmeasure); + ALLOCATE (averageNmeas, off_t, Naverage); + ALLOCATE (averageDmeas, off_t, Naverage); + ALLOCATE (averageSeqOut, off_t, Naverage); + ALLOCATE (averageRefOut, off_t, Naverage); + + // set up the measure sequence lists + for (i = 0; i < Nmeasure; i++) { + measureSeqRaw[i] = i; + measureAveRaw[i] = measure[i].averef; + } + // sort measureSeqRaw and measureAveRaw in order of measureAveRaw + SortAveMeasMatch(measureSeqRaw, measureAveRaw, Nmeasure); + + // generate the measureSeqOut array + off_t NmeasOut = 0; + for (i = 0; i < Nmeasure; i++) { + j = measureSeqRaw[i]; + if (measureDrop[j]) { + continue; + } + measureSeqOut[NmeasOut] = j; + measureRefOut[j] = NmeasOut; + NmeasOut ++; + } + // n = measureSeqOut[i] : measureOut[i] = measure[n] (i = 0 -- NmeasOut, n = 0 -- Nmeasure) + // n = measureRefOut[i] : measureOut[n] = measure[i] + ALLOCATE (measureAveOut, off_t, NmeasOut); + + // count the number of measures for each averef + N = 0; + D = -1; + currentAve = measureAveRaw[0]; + for (i = 0; i < Nmeasure; i++) { + if (measureAveRaw[i] != currentAve) { + // we have hit the next entry in the list + averageNmeas[currentAve] = N; // number of measures + averageDmeas[currentAve] = D; // first measure + N = 0; + D = -1; + currentAve = measureAveRaw[i]; + } + j = measureSeqRaw[i]; + if (measureDrop[j]) continue; + N++; + if (D == -1) { + // first valid measureOut for this averef + D = measureRefOut[j]; + } + } + averageNmeas[currentAve] = N; // number of measures + averageDmeas[currentAve] = D; // first measure + + // generate the new average sequence, skipping entries with no measurements + off_t NaveOut = 0; + for (i = 0; i < Naverage; i++) { + if (averageNmeas[i] == 0) continue; + averageSeqOut[NaveOut] = i; + averageRefOut[i] = NaveOut; + NaveOut ++; + } + // n = averageSeqOut[i] : averageOut[i] = average[n] (i = 0 -- NaveOut, n = 0 -- Naverage) + // n = averageRefOut[i] : averageOut[n] = average[i] + + // generate the output averefs + for (i = 0; i < NmeasOut; i++) { + j = measureSeqOut[i]; // measure[j] = measureOut[i] + n = measureAveRaw[j]; // average[n] : measure[j] + N = averageRefOut[n]; // averageOut[N] = average[n]; + measureAveOut[i] = N; // measureOut[i].averef = measureAveOut[i] + } + + // copy the (kept) measurements in the sorted order + ALLOCATE (measureOut, Measure, NmeasOut); + for (i = 0; i < NmeasOut; i++) { + j = measureSeqOut[i]; + measureOut[i] = measure[j]; + measureOut[i].averef = measureAveOut[i]; + } + + // copy the (kept) average entries + ALLOCATE (averageOut, Average, NaveOut); + for (i = 0; i < NaveOut; i++) { + j = averageSeqOut[i]; + averageOut[i] = average[j]; + averageOut[i].Nmeasure = averageNmeas[j]; + averageOut[i].measureOffset = averageDmeas[j]; + } + + // copy the secfilt entries for the (kept) average entries + ALLOCATE (secfiltOut, SecFilt, NaveOut*Nsecfilt); + for (i = 0; i < NaveOut; i++) { + j = averageSeqOut[i]; + for (m = 0; m < Nsecfilt; m++) { + secfiltOut[i*Nsecfilt + m] = secfilt[j*Nsecfilt + m]; + } + } + + // update the values of average.measureOffset and average.Nmeasure + FREE(measure); + FREE(average); + FREE(secfilt); + catalog[0].measure = measureOut; + catalog[0].average = averageOut; + catalog[0].secfilt = secfiltOut; + catalog[0].Nmeasure = NmeasOut; + catalog[0].Naverage = NaveOut; + + // update catalog.Nmeasure and catalog.Naverage, double check + int NmeasureTotal = 0; + int measureOffsetOK = TRUE; + int averefOK = TRUE; + for (i = 0; i < NaveOut; i++) { + NmeasureTotal += catalog[0].average[i].Nmeasure; + if (VERBOSE2 && !(NmeasureTotal <= catalog[0].Nmeasure)) { + fprintf (stderr, "too few measurements: %d %d %d\n", (int) i, NmeasureTotal, (int) catalog[0].Nmeasure); + } + measureOffsetOK &= (catalog[0].average[i].measureOffset < catalog[0].Nmeasure); + if (VERBOSE2 && !(catalog[0].average[i].measureOffset < catalog[0].Nmeasure)) { + fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure); + } + measureOffsetOK &= (catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure); + if (VERBOSE2 && !(catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure)) { + fprintf (stderr, "offset + Nmeasure too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].measureOffset, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure); + } + m = catalog[0].average[i].measureOffset; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { + averefOK &= (catalog[0].measure[m+j].averef == i); + if (VERBOSE2 && !(catalog[0].measure[m+j].averef == i)) { + fprintf (stderr, "averef broken: %d vs %d (measure %d)\n", (int) i, catalog[0].measure[m+j].averef, (int) (m+j)); + } + } + } + + if (!measureOffsetOK) { + fprintf (stderr, "ERROR: catalog %s has an invalid measureOffset\n", catalog[0].filename); + } + + if (!averefOK) { + fprintf (stderr, "ERROR: catalog %s has invalid averefs\n", catalog[0].filename); + } + + if (NmeasureTotal != catalog[0].Nmeasure) { + fprintf (stderr, "ERROR: catalog %s has an invalid Nmeasure\n", catalog[0].filename); + } + + // MARKTIME(" match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeasure, (long long) Naverage); + + catalog[0].sorted = TRUE; + + if (VERBOSE) fprintf (stderr, "ending with Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT"\n", catalog[0].Naverage, catalog[0].Nmeasure); + + off_t NdelAves = Naverage - catalog[0].Naverage; + off_t NdelMeas = Nmeasure - catalog[0].Nmeasure; + if (NdelAves || NdelMeas) { + fprintf (stderr, "deleting "OFF_T_FMT" measures and "OFF_T_FMT" averages: %s\n", NdelMeas, NdelAves, catalog[0].filename); + } + + FREE (measureDrop); + FREE (measureSeqRaw); + FREE (measureSeqOut); + FREE (measureRefOut); + FREE (measureAveRaw); + FREE (averageNmeas); + FREE (averageDmeas); + FREE (averageSeqOut); + FREE (averageRefOut); + FREE (measureAveOut); + FREE (badMeasures); + + return TRUE; +} Index: /trunk/Ohana/src/delstar/src/delete_photcodes.c =================================================================== --- /trunk/Ohana/src/delstar/src/delete_photcodes.c (revision 35757) +++ /trunk/Ohana/src/delstar/src/delete_photcodes.c (revision 35758) @@ -1,5 +1,5 @@ # include "delstar.h" -int delete_photcodes (int hostID, char *hostpath) { +int delete_photcodes () { int i; @@ -21,5 +21,5 @@ // launch the remote jobs - if (PARALLEL && !hostID) { + if (PARALLEL && !HOST_ID) { int status = delete_photcodes_parallel (skylist); return status; @@ -34,10 +34,10 @@ // does this host ID match the desired location for the table? - if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; // set up the basic catalog info char hostfile[1024]; - snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); - catalog.filename = hostID ? hostfile : skylist[0].filename[i]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; catalog.Nsecfilt = GetPhotcodeNsecfilt (); catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; Index: /trunk/Ohana/src/delstar/src/delstar.c =================================================================== --- /trunk/Ohana/src/delstar/src/delstar.c (revision 35757) +++ /trunk/Ohana/src/delstar/src/delstar.c (revision 35758) @@ -21,7 +21,68 @@ switch (MODE) { case MODE_DUP_IMAGES: - if (!delete_duplicate_images (0, NULL, &db)) exit (1); + if (!delete_duplicate_images (&db)) exit (1); exit (0); break; + case MODE_DUP_MEASURES: + if (!delete_duplicate_measures ()) exit (1); + exit (0); + break; + case MODE_FIX_LAP: { + off_t Nimage; + Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + dvo_image_unlock (&db); + + // the imageSubset here is a reduced set of fields, not a reduced set of images + ImageSubset *subset = ImagesToSubset (image, Nimage); + + if (!delete_fix_LAP (subset, Nimage)) exit (1); + SummaryImageStats (subset, Nimage); + + exit (0); + break; + } + case MODE_FIX_LAP_EDGES: { + off_t Nmeasure_edge = 0; + MeasureEdge *measure_edge = delete_fix_LAP_edges(&Nmeasure_edge); + if (!measure_edge) { + fprintf (stderr, "problem loading edge measures\n"); + exit (1); + } + delete_fix_LAP_edges_find_dups (measure_edge, Nmeasure_edge); + exit (0); + break; + } + case MODE_FIX_LAP_EDGES_DELETE: { + if (!delete_fix_LAP_edges_delete()) { + fprintf (stderr, "problem loading edge measures\n"); + exit (1); + } + exit (0); + break; + } + case MODE_FIX_LAP_STATS: { + off_t Nimage; + Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + // dvo_image_unlock (&db); + fprintf (stderr, "read images\n"); + + // the imageSubset here is a reduced set of fields, not a reduced set of images + ImageSubset *subset = ImagesToSubset (image, Nimage); + fprintf (stderr, "make subset\n"); + + if (!delete_fix_LAP_setstats (subset, Nimage)) exit (1); + SummaryImageStats_Full (image, Nimage); + + exit (0); + break; + } case MODE_IMAGEFILE: delete_imagefile (&db); Index: /trunk/Ohana/src/delstar/src/delstar_client.c =================================================================== --- /trunk/Ohana/src/delstar/src/delstar_client.c (revision 35757) +++ /trunk/Ohana/src/delstar/src/delstar_client.c (revision 35758) @@ -22,7 +22,49 @@ exit (1); } - if (!delete_duplicate_image_measures (HOST_ID, HOSTDIR, imageID)) exit (1); + if (!delete_duplicate_image_measures (imageID)) exit (1); exit (0); break; + case MODE_DUP_MEASURES: + if (!delete_duplicate_measures ()) exit (1); + exit (0); + break; + + case MODE_FIX_LAP: { + off_t Nimage; + ImageSubset *image = ImageSubsetLoad (IMAGES, &Nimage); + if (!image) { + fprintf (stderr, "failed to read image table\n"); + exit (1); + } + + if (!delete_fix_LAP (image, Nimage)) exit (1); + exit (0); + break; + } + + case MODE_FIX_LAP_EDGES: { + off_t Nmeasure_edge = 0; + MeasureEdge *measure_edge = delete_fix_LAP_edges(&Nmeasure_edge); + if (!measure_edge) { + fprintf (stderr, "problem loading edge measures\n"); + exit (1); + } + if (!MeasureEdgeSave (MEASURE_EDGE_FILE, measure_edge, Nmeasure_edge)) { + fprintf (stderr, "problem saving edge measures\n"); + exit (1); + } + exit (0); + break; + } + + case MODE_FIX_LAP_EDGES_DELETE: { + if (!delete_fix_LAP_edges_delete()) { + fprintf (stderr, "problem loading edge measures\n"); + exit (1); + } + exit (0); + break; + } + case MODE_IMAGEFILE: break; @@ -32,5 +74,5 @@ break; case MODE_PHOTCODES: - if (!delete_photcodes (HOST_ID, HOSTDIR)) { + if (!delete_photcodes ()) { fprintf (stderr, "failure deleting measurements from %s\n", HOSTDIR); exit (1);