Changeset 38753

Timestamp:

Sep 17, 2015, 8:27:18 PM (11 years ago)

Author:

eugene

Message:

substantial rework DeleteImagesByExternID to handle lensing

Location:

branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src

Files:

• dvomergeUpdate_catalogs.c (modified) (2 diffs)
• dvorepairDeleteImagesByExternID.c (modified) (7 diffs)

branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c

@@ -207 +207 @@
       }
 
-      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 +218 @@
       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]);
       }

branches/eam_branches/ipp-20150625/Ohana/src/dvomerge/src/dvorepairDeleteImagesByExternID.c

@@ -16 +16 @@
 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 +24 @@
   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 +187 @@
   }
   
-  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 +216 @@
   
     // 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 +281 @@
   }
 
-  fprintf (stderr, "Deleted %d detections\n", NdeleteTotal);
+  fprintf (stderr, "Deleted %d measure, %d lensing detections\n", NmeasureDelTotal, NlensingDelTotal);
 
   off_t NnewImage;
@@ -408 +295 @@
   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 +671 @@
 }
 
+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;
+}
+