Changeset 31258

Timestamp:

Apr 10, 2011, 10:20:54 AM (15 years ago)

Author:

eugene

Message:

add a max-density option to limit the operating data volume (randomly selected from reduced subset); add a step to rationalize the zero points: negative clouds get pushed to zero and the impact redistrubted to the rest of the images weighted by overlaps; fixed are error matching images for the 0,360 boundary

Location:

branches/eam_branches/ipp-20110404/Ohana/src/relphot

Files:

• include/relphot.h (modified) (2 diffs)
• src/MosaicOps.c (modified) (6 diffs)
• src/StarOps.c (modified) (2 diffs)
• src/args.c (modified) (1 diff)
• src/bcatalog.c (modified) (4 diffs)
• src/initialize.c (modified) (1 diff)
• src/load_images.c (modified) (1 diff)
• src/relphot.c (modified) (1 diff)
• src/select_images.c (modified) (6 diffs)

branches/eam_branches/ipp-20110404/Ohana/src/relphot/include/relphot.h

@@ -90 +90 @@
 
 PhotCode      *refPhotcode;
+
+int MaxDensityUse;
+double MaxDensityValue;
 
 int AreaSelect;
@@ -223 +226 @@
 void relphot_usage (void);
 void relphot_help (int argc, char **argv);
+
+int rationalize_mosaics ();
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog);

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/MosaicOps.c

@@ -302 +302 @@
     return;
   }
-  fprintf (stderr, "missed measurement\n");
+  // if we miss a measurement, it might be a REF type of value 
+  // fprintf (stderr, "missed measurement\n");
   return;
 }
@@ -403 +404 @@
       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
       if (mark) {
-        if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]); }
+        if (1 || VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]); }
         mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW;
         Nfew ++;
@@ -433 +434 @@
 }
   
+// When we rationalize the images/mosaics, we are driving the negative cloud images back
+// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
+// driven by the coupling of common stars. 
+int rationalize_mosaics (Catalog *catalog, int Ncatalog) {
+
+  double *imageOffset, *starOffset;
+  int *starNcount;
+  int **Slist, *NSlist, *NSLIST;
+  int i, j, k, m, nMos;
+
+  off_t Nimage;
+  Image *image;
+
+  // set a test value for now
+  float CLOUD_TOLERANCE = 0.025;
+
+  if (!MOSAIC_ZEROPT) return (FALSE);
+  if (FREEZE_MOSAICS) return (FALSE);
+
+  image = getimages (&Nimage);
+
+  ALLOCATE (imageOffset, double, Nmosaic);
+
+  ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic
+  ALLOCATE (NSlist, int,   Nmosaic); // number of stars on mosaic
+  ALLOCATE (NSLIST, int,   Nmosaic); // number of Slist entries allocated
+
+  // find the images / mosaics with negative clouds and save their offset
+  for (i = 0; i < Nmosaic; i++) {
+  
+    NSlist[i] =   0;
+    NSLIST[i] = 100;
+    ALLOCATE (Slist[i], int, NSLIST[i]);
+
+    imageOffset[i] = 0.0;
+
+    if (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE) {
+      fprintf (stderr, "cloud-free: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
+    } 
+    if (mosaic[i].Mcal < -CLOUD_TOLERANCE) {
+      imageOffset[i] = -mosaic[i].Mcal; 
+      // NOTE the negative sign: down below, we are going to add in the negative of Mcal
+      // to this image, and the propagated mean values for other images
+      fprintf (stderr, "anti-clouds: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
+    } 
+    if (mosaic[i].Mcal > CLOUD_TOLERANCE) {
+      fprintf (stderr, "cloudy    : %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
+    } 
+  }
+
+  // allocate an array for star offset
+  int Nstars = 0;
+  int NSTARS = 1000;
+  ALLOCATE (starOffset, double, NSTARS);
+  ALLOCATE (starNcount, int,    NSTARS);
+
+  // find the mean offset for each star
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      starOffset[Nstars] = 0.0;
+      starNcount[Nstars] = 0;
+
+      // skip bad stars
+      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      m = catalog[i].average[j].measureOffset;
+
+      // determine the mosaic for each measurement
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+
+        // skip unused measurements
+        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+
+        // skip REF measurements (not tied to an image)
+        if (getImageEntry (m, i) < 0) continue;
+
+        // find the source of this measurement (skip unassigned measurements)
+        nMos = bin[i][m];
+        if (nMos == -1) continue;
+
+        // accumulate the offsets from the negative cloud images (others have 0.0 value)
+        starOffset[Nstars] += imageOffset[nMos];
+        starNcount[Nstars] ++;
+        
+        // record the mosaic->star reference
+        Slist[nMos][NSlist[nMos]] = Nstars;
+        NSlist[nMos] ++;
+        if (NSlist[nMos] == NSLIST[nMos]) {
+          NSLIST[nMos] += 100;
+          REALLOCATE (Slist[nMos], int, NSLIST[nMos]);
+        }         
+      }
+      Nstars ++;
+      if (Nstars == NSTARS) {
+        NSTARS += 1000;
+        REALLOCATE (starOffset, double, NSTARS);
+        REALLOCATE (starNcount, int,    NSTARS);
+      }
+    }
+  }
+
+  // find the mean offset of the images without negative clouds
+  for (i = 0; i < Nmosaic; i++) {
+
+    // a negative cloud image (cloud: Mcal > 0; anti-clouds: Mcal < 0; imageOffset = -Mcal)
+    if (imageOffset[i] > 0.0) continue;
+
+    // we need to actually have cross-references to count
+    if (NSlist[i] < 2) continue;
+
+    float dM = 0.0;
+    for (j = 0; j < NSlist[i]; j++) {
+      Nstars = Slist[i][j];
+      if (starNcount[Nstars] > 1) {
+        dM += (starOffset[Nstars] / starNcount[Nstars]);
+      }
+    }
+    imageOffset[i] = dM / NSlist[i];
+    fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[imlist[i][0]].name, dM, NSlist[i], imageOffset[i]);
+  }
+
+  for (i = 0; i < Nmosaic; i++) {
+    fprintf (stderr, "correction: %s : %f\n", image[imlist[i][0]].name, imageOffset[i]);
+  } 
+
+  // apply all offset values to the mosaics
+  // find the images / mosaics with negative clouds and save their offset
+  for (i = 0; i < Nmosaic; i++) {
+    mosaic[i].Mcal += imageOffset[i];
+  }
+
+  for (i = 0; i < Nmosaic; i++){
+    free (Slist[i]);
+  }
+  free (NSlist);
+  free (NSLIST);
+  free (imageOffset);
+
+  return (TRUE);
+}
+
 StatType statsMosaicM (Catalog *catalog) {
 
@@ -563 +704 @@
 void clean_mosaics () {
 
-  off_t i, N, mark, Nmark;
+  off_t i, N, mark, Nmark, Nscatter, Noffset;
   double *mlist, *slist, *dlist;
   double MaxOffset, MedOffset, MaxScatter;
@@ -592 +733 @@
   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
   
-  Nmark = 0;
+  Nmark = Nscatter = Noffset = 0;
   for (i = 0; i < Nmosaic; i++) {
     mark = FALSE;
-    mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset);
+    if (mosaic[i].dMcal > MaxScatter) {
+      mark = TRUE;
+      Nscatter ++;
+    }
+    if (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset) {
+      mark = TRUE;
+      Noffset ++;
+    }
     if (mark) { 
       Nmark ++;
@@ -604 +752 @@
   }
 
-  fprintf (stderr, OFF_T_FMT" mosaics marked poor\n",  Nmark);
+  fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);
   initstats (STATMODE);
   free (mlist);

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c

@@ -5 +5 @@
 static double *dlist;
 
+// When we rationalize the images/mosaics, we are driving the negative cloud images back
+// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
+// driven by the coupling of common stars.  This array carries the impact of those offsets
+// on each star
+static double *Moffset;
+
 void initMrel (Catalog *catalog, int Ncatalog) {
 
@@ -16 +22 @@
   }
 
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
+  ALLOCATE (list,    double, MAX (1, Nmax));
+  ALLOCATE (dlist,   double, MAX (1, Nmax));
+  ALLOCATE (Moffset, double, MAX (1, Nmax));
 }  
 

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/args.c

@@ -128 +128 @@
   }
 
+  MaxDensityUse = FALSE;
+  if ((N = get_argument (argc, argv, "-max-density"))) {
+    remove_argument (N, &argc, argv);
+    MaxDensityValue = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+    MaxDensityUse = TRUE;
+  }
+
   SHOW_PARAMS = FALSE;
   if ((N = get_argument (argc, argv, "-params"))) {

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/bcatalog.c

@@ -1 +1 @@
 # include "relphot.h"
+
+extern double drand48();
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
@@ -70 +72 @@
       // check for galaxies
       if (!isnan(catalog[0].measure[offset].Map)) {
-          if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) {
-              nEXT ++;
-          } else {
-              nPSF ++;
-          }
+        if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) {
+          nEXT ++;
+        } else {
+          nPSF ++;
+        }
       }
 
@@ -139 +141 @@
   assert (PhotNsec == catalog[0].Nsecfilt);
 
+  // limit the total number of stars in the catalog
+  if (MaxDensityUse) {
+    LimitDensityCatalog (subcatalog, catalog);
+  }
+
   if (VERBOSE) {
     fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", 
@@ -147 +154 @@
   return (TRUE);
 }
+
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+
+  Catalog tmpcatalog;
+
+  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);
+
+  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
+
+  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
+  assert (AREA > 0);
+
+  off_t Nmax = MaxDensityValue * AREA;
+  if (subcatalog[0].Naverage <= Nmax) {
+    if (VERBOSE) {
+      fprintf (stderr, "subcatalog has less than the max density\n");
+    }
+    return (TRUE);
+  }
+
+  off_t Naverage = subcatalog[0].Naverage;
+
+  // select a random subset of Nmax stars from subcatalog using Fisher-Yates
+
+  // we are going to select Nmax entries by generating a random-sorted index list
+  off_t *index, tmp, i, j, ave;
+  ALLOCATE (index, off_t, Naverage);
+  for (i = 0; i < Naverage; i++) {
+    index[i] = i;
+  }
+  for (i = 0; i < Naverage; i++) {
+    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
+    tmp = index[j];
+    index[j] = index[i];
+    index[i] = tmp;
+  }
+
+  // count the number of measurements this selection will yield
+  off_t NMEASURE = 0;
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    NMEASURE += subcatalog[0].average[ave].Nmeasure;
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.average, Average, Nmax);
+  ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * PhotNsec);
+  ALLOCATE (tmpcatalog.measure, Measure, NMEASURE);
+
+  off_t Nmeasure = 0;
+
+  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    tmpcatalog.average[i] = subcatalog[0].average[ave];
+    tmpcatalog.average[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].average[ave].measureOffset + j;
+      tmpcatalog.measure[Nmeasure] = subcatalog[0].measure[offset];
+      tmpcatalog.measure[Nmeasure].averef = i;
+      Nmeasure ++;
+    }
+  }
+
+  if (VERBOSE) {
+    fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", 
+             Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  catalog[0].filename);
+  }
+
+  free (subcatalog[0].average);
+  free (subcatalog[0].measure);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].average = tmpcatalog.average;
+  subcatalog[0].measure = tmpcatalog.measure;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nmax;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
+
+  return (TRUE);
+}

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/initialize.c

@@ -55 +55 @@
     exit (0);
   }
+
+  // init the random seed
+  long A, B;
+  A = time(NULL);
+  for (B = 0; A == time(NULL); B++);
+  srand48(B);
 }
-

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/load_images.c

@@ -44 +44 @@
   // select the images which overlap the selected sky regions
   subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
-  MARKTIME("selected images: %f sec\n", dtime);
+  MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
 
   // generate db->vtable from db->ftable based on the selection

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/relphot.c

@@ -140 +140 @@
       plot_chisq (catalog, Ncatalog);
     }
-    if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); 
-    if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
-    if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
-    if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
+    if (i % 6 == 1) rationalize_mosaics (catalog, Ncatalog);
+    // if (i % 6 == 1) rationalize_images ();
+    if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE); 
+    if (i % 6 == 3) clean_stars (catalog, Ncatalog);
+    if (i % 6 == 5) clean_mosaics ();
+    if (i % 6 == 5) clean_images ();
+
+    // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); 
+    // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
+    // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
+    // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
     global_stats (catalog, Ncatalog);
   }
 
-  SAVEPLOT = TRUE;
-  plot_scatter (catalog, Ncatalog); 
-  plot_grid (catalog); 
-  plot_mosaics ();
-  plot_images ();
-  plot_stars (catalog, Ncatalog);
-  plot_chisq (catalog, Ncatalog);
-
+  if (PLOTSTUFF) {
+    plot_scatter (catalog, Ncatalog); 
+    plot_grid (catalog); 
+    plot_mosaics ();
+    plot_images ();
+    plot_stars (catalog, Ncatalog);
+    plot_chisq (catalog, Ncatalog);
+  }
+  
   if (USE_GRID) dump_grid ();
   if (!UPDATE) exit (0);

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/select_images.c

@@ -64 +64 @@
   DmaxSkyRegion = -90.0;
 
+  // FILE *ftest = fopen ("relphot.dump.dat", "w");
+
   /* compare with each region file */
   for (i = 0; i < skylist[0].Nregions; i++) { 
@@ -131 +133 @@
     }
 
+    // XXX temporary test : record center coords for each accepted and rejected chip/mosaic
+    // double RAo, DECo;
+
     /* define image corners - note the DIS images (mosaic phu) are special */
     if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) {
@@ -138 +143 @@
       Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY;
       Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY;
+      // XY_to_RD(&RAo, &DECo, 0.0, 0.0, &timage[i].coords);
     } else {
       Xi[0] = 0;            Yi[0] = 0;
@@ -144 +150 @@
       Xi[3] = 0;            Yi[3] = timage[i].NY;
       Xi[4] = 0;            Yi[4] = 0;
+      // XY_to_RD(&RAo, &DECo, 0.5*timage[i].NX, 0.5*timage[i].NY, &timage[i].coords);
     }
     found = FALSE;
 
     /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */
-    double RminImage = 360.0;
-    double RmaxImage =   0.0;
+    double RminImage = RmidSkyRegion + 180.0;
+    double RmaxImage = RmidSkyRegion - 180.0;
     double DminImage = +90.0;
     double DmaxImage = -90.0;
@@ -218 +225 @@
 
   found_it:
+    // XXX We claim this is a good image: write to a test file
+    // fprintf (ftest, "%s : %lf %lf  : %f %f %d %x\n", timage[i].name, RAo, DECo, timage[i].Mcal, timage[i].dMcal, timage[i].nstar, timage[i].flags);
+
     image[nimage] = timage[i]; 
     /* always allow 'few' images to succeed, if possible */
@@ -239 +249 @@
   }
   MARKTIME("finish image selection: %f sec\n", dtime);
+
+  // fclose (ftest);
 
   if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);