Changeset 33610

Timestamp:

Mar 24, 2012, 10:27:28 PM (14 years ago)

Author:

eugene

Message:

setMrelFinal was marking all ref photometry to be skipped, so objects with only SYNTH phot were not getting average photometry; bcatalog sort by Nmeasure

Location:

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

Files:

• StarOps.c (modified) (6 diffs)
• bcatalog.c (modified) (6 diffs)
• relphot.c (modified) (1 diff)
• setMrelFinal.c (modified) (3 diffs)

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

@@ -153 +153 @@
 }
 
+int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
+
+  off_t k;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  off_t m = average[0].measureOffset;
+
+  for (k = 0; k < average[0].Nmeasure; k++, m++) {
+    fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
+  }
+
+  int Ns;
+  for (Ns = 0; Ns < Nsecfilt; Ns++) {
+    fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
+  }
+  return 1;
+}
+
+# define MARK_SKIP_MEAS \
+  catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \
+  catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
+
 // setMrel and setMrelOutput are extremely similar, but have slightly different implications:
 // * setMrel uses the internal Tiny structures only
@@ -190 +213 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
 
+      if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) {
+        fprintf (stderr, "test obj\n");
+        print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure);
+      }
+
       int minUbercalDist = 1000;
 
@@ -205 +233 @@
         int haveSynth = FALSE;
         int haveUbercal = FALSE;
+        int isRefPhot = FALSE;
 
         N = 0;
@@ -221 +250 @@
             // these detetions have no image (eg, ref values such as 2MASS)
             Mcal = Mmos = Mgrid = 0;
+            isRefPhot = TRUE;
           } else {
             Mcal  = getMcal  (m, i, flatcorr, catalog);
-            if (isnan(Mcal)) continue;
+            if (isnan(Mcal)) {
+              MARK_SKIP_MEAS;
+              continue;
+            }
             Mmos  = getMmos  (m, i);
-            if (isnan(Mmos)) continue;
+            if (isnan(Mmos)) {
+              MARK_SKIP_MEAS;
+              continue;
+            }
             Mgrid = getMgrid (m, i);
-            if (isnan(Mgrid)) continue;
+            if (isnan(Mgrid)) {
+              MARK_SKIP_MEAS;
+              continue;
+            }
           }
 
           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
-          if (isnan(Msys)) { continue; }
+          if (isnan(Msys)) { 
+            if (!isRefPhot) {
+              MARK_SKIP_MEAS;
+            }
+            continue; 
+          }
           list[N] = Msys - Mcal - Mmos - Mgrid;
 
@@ -255 +299 @@
             if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) {
               // detections without one of these bits should only be used in PASS_1
+              MARK_SKIP_MEAS;
               continue;
             }
@@ -260 +305 @@
           // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars)
           if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) {
-            if (pass < 4) continue;
+            if (pass < 4) {
+              MARK_SKIP_MEAS;
+              continue;
+            }
             haveSynth = TRUE;
           }

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

@@ -1 +1 @@
 # include "relphot.h"
+
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
@@ -212 +214 @@
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
-    LimitDensityCatalog (subcatalog, catalog);
+    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
   }
 
@@ -224 +226 @@
 }
 
-int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+// sort by decreasing Nmeasure (X)
+void sort_by_Nmeasure (int *X, off_t *Y, off_t N) {
+
+# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;   \
+  tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
+  tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
+}
+# define COMPARE(A,B)(X[A] > X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) {
 
   Catalog tmpcatalog;
@@ -232 +250 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  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);
+  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&oldcatalog[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);
@@ -255 +273 @@
 
   // we are going to select Nmax entries by generating a random-sorted index list
-  off_t *index, tmp, i, j, ave;
+  int *value;
+  off_t *index, i, j, ave;
   ALLOCATE (index, off_t, Naverage);
+  ALLOCATE (value, int, 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;
-  }
+    value[i] = subcatalog[0].averageT[i].Nmeasure;
+  }
+  sort_by_Nmeasure (value, index, Naverage);
 
   // count the number of measurements this selection will yield
@@ -298 +314 @@
 
   if (VERBOSE) {
+    char *basename = filebasename (oldcatalog[0].filename);
+    fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", 
+             Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
+    free (basename);
+  }
+
+  free (index);
+  free (value);
+  free (subcatalog[0].averageT);
+  free (subcatalog[0].measureT);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].averageT = tmpcatalog.averageT;
+  subcatalog[0].measureT = tmpcatalog.measureT;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nmax;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
+
+  return (TRUE);
+}
+
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+
+  Catalog tmpcatalog;
+
+  double Rmin, Rmax, Dmin, Dmax;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  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].averageT[ave].Nmeasure;
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.averageT, AverageTiny, Nmax);
+  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt,  SecFilt, Nmax * Nsecfilt);
+
+  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.averageT[i] = subcatalog[0].averageT[ave];
+    tmpcatalog.averageT[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.averageT[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].averageT[ave].measureOffset + j;
+      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
+      tmpcatalog.measureT[Nmeasure].averef = i;
+      Nmeasure ++;
+    }
+    for (j = 0; j < Nsecfilt; j++) {
+      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
+    }
+  }
+
+  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);

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

@@ -230 +230 @@
       MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
     }
+
+    /* at this point, we have correct cal coeffs in the image/mosaic structures */
+    for (i = 0; i < Ncatalog; i++) {
+      free_tiny_values (&catalog[i]);
+      dvo_catalog_free (&catalog[i]);
+    }
+    freeImageBins (Ncatalog, TRUE);
+    freeMosaicBins (Ncatalog, TRUE);
+    freeGridBins (Ncatalog);
+
   } // (NLOOP > 0) : this loop determines the offsets per chip 
 
-  // only change the real database if -update is requested
+  reload_images (&db);
+
+  // only change the real database files if -update is requested
   if (!UPDATE) exit (0);
   
-  reload_images (&db);
-
-  /* at this point, we have correct cal coeffs in the image/mosaic structures */
-  for (i = 0; i < Ncatalog; i++) {
-    free_tiny_values (&catalog[i]);
-    dvo_catalog_free (&catalog[i]);
-  }
-  freeImageBins (Ncatalog, TRUE);
-  freeMosaicBins (Ncatalog, TRUE);
-  freeGridBins (Ncatalog);
-
   /* load catalog data from region files, update Mrel include all data */
   reload_catalogs (skylist, flatcorr, 0, NULL);

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c

@@ -123 +123 @@
 // non-outliers with Imag range
 
+// int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
+// 
+//   off_t k;
+// 
+//   int Nsecfilt = GetPhotcodeNsecfilt ();
+// 
+//   off_t m = average[0].measureOffset;
+// 
+//   for (k = 0; k < average[0].Nmeasure; k++, m++) {
+//     fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
+//   }
+// 
+//   int Ns;
+//   for (Ns = 0; Ns < Nsecfilt; Ns++) {
+//     fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
+//   }
+//   return 1;
+// }
+
+int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
+
 void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) {
 
@@ -137 +158 @@
   
   /* allow measures marked as outliers (POOR) and off image region (AREA) */
-  if (pass >= 2) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
+  if (pass >= 2) {
+    MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
+  } else {
+    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
+  }
 
   /* mark measures which should be ignored on this pass */
   for (i = 0; i < catalog[0].Naverage; i++) {
     Ntot += catalog[0].average[i].Nmeasure;
+
+    if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
+      fprintf (stderr, "test obj\n");
+      print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
+    }
 
     // mark measurements for each secfilt separately
@@ -173 +203 @@
         // clear SKIP for all measures at first 
         catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-
-        // skip measurements from BAD images and mosaics 
-        // do NOT skip measurements without a matching image 
-        if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
-        if (isnan(getMmos (m, 0))) goto skip;
-
-        // PASS 4 : skip measurements by inst mag limit 
+        catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
+
+        // skip measurements from BAD images and mosaics (not REF mags)
+        // do NOT skip measurements without a matching image (REF mags)
+        off_t Nim = getImageEntry (m, 0);
+        if (Nim > -1) {
+          if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
+          if (isnan(getMmos (m, 0))) goto skip;
+        }
+
+        // PASS 4 : skip measurements by inst mag limit (not REF mags)
         if ((pass < 4) && ImagSelect) {
-          mag = PhotInst (&catalog[0].measure[m]);
-          if (mag < ImagMin) goto skip;
-          if (mag > ImagMax) goto skip;
+          if (Nim > -1) {
+            mag = PhotInst (&catalog[0].measure[m]);
+            if (mag < ImagMin) goto skip;
+            if (mag > ImagMax) goto skip;
+          }
         }