Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 38752) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 38753) @@ -207,5 +207,5 @@ } - if (!dvo_catalog_backup (&outcatalog, TRUE)) { + if (!dvo_catalog_backup (&outcatalog, "~", TRUE)) { fprintf (stderr, "ERROR: failed to make backup for catalog %s\n", outlist[0].filename[j]); exit (1); @@ -218,5 +218,5 @@ SetProtect (FALSE); - if (!dvo_catalog_unlink_backup (&outcatalog, TRUE)) { + if (!dvo_catalog_unlink_backup (&outcatalog, "~", TRUE)) { fprintf (stderr, "WARNING: failed to remove backup for catalog %s\n", outlist[0].filename[j]); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c (revision 38752) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c (revision 38753) @@ -16,4 +16,7 @@ int SaveImageTable (Image *image, off_t Nimage, char *catdir, FITS_DB *oldDB); Image *DeleteSelectedImages (Image *image, off_t Nimage, int *deleteImage, myIndexType *imageIDindex, off_t *NnewImage); +int DeleteLensing (Catalog *catalog, myIndexType *imageIDindex, int *deleteImage, int *nDelete); +int DeleteMeasure (Catalog *catalog, myIndexType *imageIDindex, int *deleteImage, int *nDelete); +int RepairAverage (Catalog *catalog); int dvorepairDeleteImagesByExternID (int argc, char **argv) { @@ -21,22 +24,10 @@ FITS_DB db; // database handle pointing to input image table - int i, j, N, NmeasureNew, Ndelete, Nvalid; - off_t Nimage, Nmeasure; + int i, N; + off_t Nimage; int nPass, raPass; Image *image; - Measure *measure; - Measure *measureNew; - - char *cpmFilenameSrc = NULL; - char *cptFilenameSrc = NULL; - char *cpsFilenameSrc = NULL; - char *cpmFilenameTgt = NULL; - char *cptFilenameTgt = NULL; - char *cpsFilenameTgt = NULL; - char *imageFilename = NULL; - - DVOCatFormat catformat; N = get_argument (argc, argv, "-delete-images-by-extern-id"); @@ -196,12 +187,6 @@ } - ALLOCATE(cpmFilenameSrc, char, strlen(catdir) + 64); - ALLOCATE(cptFilenameSrc, char, strlen(catdir) + 64); - ALLOCATE(cpsFilenameSrc, char, strlen(catdir) + 64); - ALLOCATE(cpmFilenameTgt, char, strlen(catdir) + 64); - ALLOCATE(cptFilenameTgt, char, strlen(catdir) + 64); - ALLOCATE(cpsFilenameTgt, char, strlen(catdir) + 64); - - int NdeleteTotal = 0; + int NmeasureDelTotal = 0; + int NlensingDelTotal = 0; for (raPass = 0; raPass < nPass; raPass++) { @@ -231,162 +216,64 @@ // loop over the populated input regions - int Ncheck = 0; for (i = 0; i < inlist[0].Nregions; i++) { if (!inlist[0].regions[i][0].table) continue; - sprintf (cpmFilenameSrc, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name); - sprintf (cptFilenameSrc, "%s/%s.cpt", catdir, inlist[0].regions[i][0].name); - sprintf (cpsFilenameSrc, "%s/%s.cps", catdir, inlist[0].regions[i][0].name); - sprintf (cpmFilenameTgt, "%s/%s.cpm.broken", catdir, inlist[0].regions[i][0].name); - sprintf (cptFilenameTgt, "%s/%s.cpt.broken", catdir, inlist[0].regions[i][0].name); - sprintf (cpsFilenameTgt, "%s/%s.cps.broken", catdir, inlist[0].regions[i][0].name); - - Header cpmHeaderPHU; - Header cpmHeaderTBL; - FTable cpmFtable; - - /*** read and examine the CPM file ***/ - { - if (VERBOSE) fprintf (stderr, "check %s\n", cpmFilenameSrc); - - cpmFtable.header = &cpmHeaderTBL; - - // open cpm file - FILE *cpmFile = fopen(cpmFilenameSrc, "r"); - if (!cpmFile) continue; - // myAssert(cpmFile, "failed to open cpm file"); - - // load the cpm header - if (!gfits_fread_header (cpmFile, &cpmHeaderPHU)) { - myAbort("failure to cpm header"); - } - - // move to TBL header - int Nbytes = cpmHeaderPHU.datasize + gfits_data_size (&cpmHeaderPHU); - fseeko (cpmFile, Nbytes, SEEK_SET); - - // read cpm TBL header - if (!gfits_fread_header (cpmFile, &cpmHeaderTBL)) { - myAbort("can't read header for cpm table"); - } - - char nativeOrder; - Catalog catalog; - dvo_catalog_init (&catalog, TRUE); - catalog.f = cpmFile; - if (!gfits_fread_uncompressed (&catalog, &cpmFtable, &nativeOrder, TRUE)) { - myAbort("can't read data for cpm table"); - } - - int NbytesPerRow; - gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow); - - char compressMode[256]; - if (gfits_scan (&cpmHeaderTBL, "DVO_CMP", "%s", 1, compressMode)) { - if (strcmp (compressMode, "NONE")) { - myAbort ("fix compression"); - } - } - - // this function can only handle PS1_V5 and later formats - // NOTE: FtableToMeasure (and equivalent) free the data associated with cpmFtable. - // We are left with the data in measure. - measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat, nativeOrder); - myAssert(measure, "failed to convert ftable to measure data"); - - Nvalid = (int)(cpmFtable.validsize / NbytesPerRow); - Nvalid = MIN(Nmeasure, Nvalid); - - // close the input cpm file - fclose(cpmFile); - - // allocate an output array of measures (to replace, if needed) - ALLOCATE (measureNew, Measure, Nvalid); - - NmeasureNew = 0; - Ndelete = 0; - - // examine all measurements: find ones that need to be deleted - for (j = 0; j < Nvalid; j++) { - int imageID = measure[j].imageID; - myAssert(imageID, "measure is missing an image ID"); - // XXX this case is valid if we have REF detections (no associated image) - - int N = myIndexGetEntry(imageIDindex, imageID); - if (N < 0) { - // this detection comes from a non-existant image, delete - Ndelete ++; - continue; - } - - if (deleteImage[N]) { - Ndelete ++; - continue; - } - measureNew[NmeasureNew] = measure[j]; - NmeasureNew ++; - } + char filename[DVO_MAX_PATH]; + snprintf (filename, DVO_MAX_PATH, "%s/%s.cpt", catdir, inlist[0].regions[i][0].name); + + Catalog catalog; + + // set up the basic catalog info + dvo_catalog_init (&catalog, TRUE); + catalog.filename = filename; // XXX only allow non-parallel dbs for now + catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_SECFILT | DVO_LOAD_LENSING; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + if (!dvo_catalog_open (&catalog, inlist[0].regions[i], VERBOSE, "r")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); + exit (1); } - - NdeleteTotal += Ndelete; - - // if we actually want to delete any measurements, write out a new cpm file - if (Ndelete > 0) { - - fprintf (stderr, "deleting %d of %d (keep %d) detections from %s -> %s\n", (int) Ndelete, (int) Nvalid, (int) NmeasureNew, cpmFilenameSrc, cpmFilenameTgt); + if (!catalog.Naverage_disk) { + if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } - // the CPT and CPS tables need to be regenerated. This must happen first because, in the process, we also update measure->averef - // what about other tables? - RepairTableCPT_V1(cptFilenameSrc, cptFilenameTgt, cpsFilenameSrc, cpsFilenameTgt, measureNew, NmeasureNew, image, Nimage, imageIDindex, catformat); - - // convert internal to external format - if (!MeasureToFtable (&cpmFtable, NULL, measureNew, NmeasureNew, catformat, TRUE)) { - myAbort("trouble converting format"); - } - - // rename the old cpm file: - if (rename (cpmFilenameSrc, cpmFilenameTgt)) { - perror ("tried to rename file"); - exit (2); - } - - // create and write the output file - FILE *cpmFile = fopen(cpmFilenameSrc, "w"); - myAssert(cpmFile, "failed to open cpt file"); - - // write PHU header - if (!gfits_fwrite_header (cpmFile, &cpmHeaderPHU)) { - myAbort("can't write primary header"); - } - - Matrix matrix; - - // write the PHU matrix; this is probably a NOP, do I have to keep it in? - gfits_create_matrix (&cpmHeaderPHU, &matrix); - if (!gfits_fwrite_matrix (cpmFile, &matrix)) { - myAbort("can't write primary matrix"); - } - gfits_free_matrix (&matrix); - - // write the table data - if (!gfits_fwrite_ftable_range (cpmFile, &cpmFtable, 0, NmeasureNew, 0, NmeasureNew)) { - myAbort("can't write table data"); - } - fclose (cpmFile); - - gfits_free_table (&cpmFtable); - } else { - if (VERBOSE) fprintf (stderr, "nothing to delete in %s\n", cpmFilenameSrc); + /*** delete the measure and lensing entries matching the list of imageIDs ***/ + int NmeasureDel = 0; + DeleteMeasure (&catalog, imageIDindex, deleteImage, &NmeasureDel); + + int NlensingDel = 0; + DeleteLensing (&catalog, imageIDindex, deleteImage, &NlensingDel); + + fprintf (stderr, "deleting %d measure, %d lensing (keep %d, %d)\n", NmeasureDel, NlensingDel, (int) catalog.Nmeasure, (int) catalog.Nlensing); + + if (!NmeasureDel && !NlensingDel) { + if (VERBOSE) fprintf (stderr, "nothing to delete in %s\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; } - gfits_free_header (&cpmHeaderPHU); - gfits_free_header (&cpmHeaderTBL); - free (measure); - free (measureNew); - - Ncheck ++; - if (Ncheck % 1000 == 0) { - fprintf (stderr, "%s...", inlist[0].regions[i][0].name); + NmeasureDelTotal += NmeasureDel; + NlensingDelTotal += NlensingDel; + + // the CPT and CPS tables need to be regenerated. This must happen first because, in the process, we also update measure->averef + // what about other tables? + RepairAverage (&catalog); + + if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename); + + if (!dvo_catalog_backup (&catalog, ".undel", TRUE)) { + fprintf (stderr, "ERROR: failed to make backup for catalog %s\n", catalog.filename); + exit (1); } + + SetProtect (TRUE); + if (!dvo_catalog_save (&catalog, VERBOSE)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog.filename); exit (1); } + if (!dvo_catalog_unlock (&catalog)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog.filename); exit (1); } + SetProtect (FALSE); + dvo_catalog_free (&catalog); } fprintf (stderr, "\n"); @@ -394,5 +281,5 @@ } - fprintf (stderr, "Deleted %d detections\n", NdeleteTotal); + fprintf (stderr, "Deleted %d measure, %d lensing detections\n", NmeasureDelTotal, NlensingDelTotal); off_t NnewImage; @@ -408,14 +295,9 @@ free (deleteIDs); - free(cpmFilenameSrc); - free(cptFilenameSrc); - free(cpsFilenameSrc); - free(cpmFilenameTgt); - free(cptFilenameTgt); - free(cpsFilenameTgt); - SkyTableFree(insky); gfits_db_free (&db); + + ohana_memcheck (TRUE); exit (0); @@ -789,2 +671,223 @@ } +int RepairAverage (Catalog *catalog) { + + off_t i, j, Nave; + + // average is sorted so averef is valid + // secfilt is also sorted so sequence is valid (but this is kind of moot since it will be re-calculated) + // average.measureOffset is NOT valid + // average.lensingOffset is NOT valid + + Average *average = catalog->average; + Measure *measure = catalog->measure; + Lensing *lensing = catalog->lensing; + + // reset the values of average.Nmeasure, average.Nlensing + for (i = 0; i < catalog->Naverage; i++) { + average[i].Nmeasure = 0; + average[i].Nlensing = 0; + } + + // re-calculate average.Nmeasure + for (i = 0; i < catalog->Nmeasure; i++) { + Nave = measure[i].averef; + myAssert (measure[i].objID == average[Nave].objID, "invalid measure:average match"); + average[Nave].Nmeasure ++; + } + + // re-calculate average.Nlensing + for (i = 0; i < catalog->Nlensing; i++) { + Nave = lensing[i].averef; + myAssert (lensing[i].objID == average[Nave].objID, "invalid lensing:average match"); + average[Nave].Nlensing ++; + } + + // create a copy of the average table, keeping only the entries with Nmeasure & Nlensing > 0 + ALLOCATE_PTR (averageNew, Average, catalog->Naverage); + ALLOCATE_PTR (averefNew, int, catalog->Naverage); + + Nave = 0; + for (i = 0; i < catalog->Naverage; i++) { + averefNew[i] = -1; + if (!average[i].Nmeasure && !average[i].Nlensing) continue; + + averageNew[Nave] = average[i]; + averefNew[i] = Nave; + Nave ++; + } + off_t NaverageNew = Nave; + + // update measure.averef values + for (i = 0; i < catalog->Nmeasure; i++) { + Nave = averefNew[measure[i].averef]; + myAssert (Nave >= 0, "oops"); + measure[i].averef = Nave; + } + + // update lensing.averef values + for (i = 0; i < catalog->Nlensing; i++) { + Nave = averefNew[lensing[i].averef]; + myAssert (Nave >= 0, "oops"); + lensing[i].averef = Nave; + } + + // XXX need to update measureOffset and lensingOffset + + // measure[] should be blocked and sequential: + // measure[i+1].averef >= measure[i].averef + + // update average.measureOffset values + Nave = -1; + for (i = 0; i < catalog->Nmeasure; i++) { + if (measure[i].averef == Nave) continue; + Nave = measure[i].averef; + averageNew[Nave].measureOffset = i; + } + + // update average.lensingOffset values + Nave = -1; + for (i = 0; i < catalog->Nlensing; i++) { + if (lensing[i].averef == Nave) continue; + Nave = lensing[i].averef; + averageNew[Nave].lensingOffset = i; + } + + // check the result (measure -> average) + for (i = 0; i < catalog->Nmeasure; i++) { + Nave = measure[i].averef; + myAssert(averageNew[Nave].objID == measure[i].objID, "objIDs do not match"); + myAssert(averageNew[Nave].catID == measure[i].catID, "catIDs do not match"); + } + // check the result (average -> measure) + for (i = 0; i < NaverageNew; i++) { + int m = averageNew[i].measureOffset; + for (j = 0; j < averageNew[i].Nmeasure; j++) { + myAssert(averageNew[i].objID == measure[j+m].objID, "objIDs do not match"); + myAssert(averageNew[i].catID == measure[j+m].catID, "catIDs do not match"); + myAssert(measure[j+m].averef == i, "averef broken"); + } + } + + // check the result (lensing -> average) + for (i = 0; i < catalog->Nlensing; i++) { + Nave = lensing[i].averef; + myAssert(averageNew[Nave].objID == lensing[i].objID, "objIDs do not match"); + myAssert(averageNew[Nave].catID == lensing[i].catID, "catIDs do not match"); + } + // check the result (average -> lensing) + for (i = 0; i < NaverageNew; i++) { + int m = averageNew[i].lensingOffset; + for (j = 0; j < averageNew[i].Nlensing; j++) { + myAssert(averageNew[i].objID == lensing[j+m].objID, "objIDs do not match"); + myAssert(averageNew[i].catID == lensing[j+m].catID, "catIDs do not match"); + myAssert(lensing[j+m].averef == i, "averef broken"); + } + } + + free (catalog->secfilt); + ALLOCATE (catalog->secfilt, SecFilt, NaverageNew*catalog->Nsecfilt); + for (i = 0; i < NaverageNew*catalog->Nsecfilt; i++) { + dvo_secfilt_init (&catalog->secfilt[i], SECFILT_RESET_ALL); + } + + free (averefNew); + free (catalog->average); + catalog->average = averageNew; + catalog->Naverage = NaverageNew; + + return (TRUE); +} + +// delete measure entries +int DeleteMeasure (Catalog *catalog, myIndexType *imageIDindex, int *deleteImage, int *nDelete) { + + int j; + + Measure *measure = catalog->measure; + + // allocate an output array of measures (to replace, if needed) + ALLOCATE_PTR (measureNew, Measure, catalog->Nmeasure); + + int NmeasureNew = 0; + int NmeasureDel = 0; + + // examine all measurements: find ones that need to be deleted + for (j = 0; j < catalog->Nmeasure; j++) { + int imageID = measure[j].imageID; + myAssert(imageID, "measure is missing an image ID"); + // XXX this case is valid if we have REF detections (no associated image) + + int N = myIndexGetEntry(imageIDindex, imageID); + if (N < 0) { + // this detection comes from a non-existant image; delete + NmeasureDel ++; + continue; + } + + // measure matches a bad image; delete + if (deleteImage[N]) { + NmeasureDel ++; + continue; + } + + // keep this measure + measureNew[NmeasureNew] = measure[j]; + NmeasureNew ++; + } + + free (catalog->measure); + catalog->measure = measureNew; + catalog->Nmeasure = NmeasureNew; + + *nDelete = NmeasureDel; + + return TRUE; +} + +// delete lensing entries +int DeleteLensing (Catalog *catalog, myIndexType *imageIDindex, int *deleteImage, int *nDelete) { + + int j; + + Lensing *lensing = catalog->lensing; + + // allocate an output array of measures (to replace, if needed) + ALLOCATE_PTR (lensingNew, Lensing, catalog->Nlensing); + + int NlensingNew = 0; + int NlensingDel = 0; + + // examine all lensingments: find ones that need to be deleted + for (j = 0; j < catalog->Nlensing; j++) { + int imageID = lensing[j].imageID; + myAssert(imageID, "lensing is missing an image ID"); + // XXX this case is valid if we have REF detections (no associated image) + + int N = myIndexGetEntry(imageIDindex, imageID); + if (N < 0) { + // this detection comes from a non-existant image; delete + NlensingDel ++; + continue; + } + + // lensing matches a bad image; delete + if (deleteImage[N]) { + NlensingDel ++; + continue; + } + + // keep this lensing + lensingNew[NlensingNew] = lensing[j]; + NlensingNew ++; + } + + free (catalog->lensing); + catalog->lensing = lensingNew; + catalog->Nlensing = NlensingNew; + + *nDelete = NlensingDel; + + return TRUE; +} +