Changeset 39641 for trunk

Timestamp:

Jul 29, 2016, 8:07:10 AM (10 years ago)

Author:

eugene

Message:

add 2 test image options; add manual uniquer for easy restart; set to use IRLS for image and object averages

Location:

trunk/Ohana/src/relphot

Files:

• include/relphot.h (modified) (3 diffs)
• src/ImageOps.c (modified) (10 diffs)
• src/args.c (modified) (4 diffs)
• src/bcatalog.c (modified) (1 diff)
• src/liststats.c (modified) (5 diffs)
• src/load_catalogs.c (modified) (1 diff)
• src/relphot_images.c (modified) (4 diffs)

trunk/Ohana/src/relphot/include/relphot.h

@@ -276 +276 @@
 int    PARALLEL_MANUAL;
 int    PARALLEL_SERIAL;
+char  *MANUAL_UNIQUER;
 
 int    PARALLEL_REGIONS_MANUAL;
@@ -286 +287 @@
 int    MOSAIC_ZEROPT;
 
-int    TEST_IMAGE;
+int    TEST_IMAGE1;
+int    TEST_IMAGE2;
 
 int    NLOOP;
@@ -309 +311 @@
 int    FREEZE_IMAGES;
 int    FREEZE_MOSAICS;
+int    CALIBRATE_STACKS_AND_WARPS;
 int    USE_GRID;
 int    KEEP_UBERCAL;

trunk/Ohana/src/relphot/src/ImageOps.c

@@ -1 +1 @@
 # include "relphot.h"
+# define BASIC_STATS 0
+
 void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
 
@@ -614 +616 @@
 
   off_t i, j, m, c, n, N, Nmax;
-  int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys, Nbright;
+  int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys;
   float Msys, Mrel, Mmos, Mgrid, Mflat;
-  double *list, *dlist, *Mlist, *dMlist;
+  // double *list, *dlist, *Mlist, *dMlist;
 
   StatType stats;
@@ -637 +639 @@
     Nmax = MAX (Nmax, N_onImage[i]);
   }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-  ALLOCATE (Mlist, double, Nmax);
-  ALLOCATE (dMlist, double, Nmax);
+  StatDataSet *refStars = StatDataSetAlloc (1, Nmax);
+  StatDataSet *brightStars = StatDataSetAlloc (1, Nmax);
 
   Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
@@ -647 +647 @@
   for (i = 0; i < Nimage; i++) {
     
+    if (image[i].imageID == TEST_IMAGE1) {
+      fprintf (stderr, "test image 1\n");
+    }
+    if (image[i].imageID == TEST_IMAGE2) {
+      fprintf (stderr, "test image 1\n");
+    }
+
     /* on PoorImages run, skip good images */
     if (PoorImages) {
@@ -664 +671 @@
 
     int minUbercalDist = 1000;
-
-    // number of stars to measure the bright-end scatter
-    Nbright = 0;
-
-    N = 0;
+    
+    off_t Nref = 0;  // total number of reference stars on the image
+    int Nbright = 0; // number of stars to measure the bright-end scatter
+
     for (j = 0; j < N_onImage[i]; j++) {
       
@@ -716 +722 @@
 
     skip:
-      list[N] = Msys - Mrel - Mmos - Mgrid + Mflat;
-      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-      if (fabs(list[N]) > 0.03) {
+      refStars->flxlist[Nref] = Msys - Mrel - Mmos - Mgrid + Mflat;
+      refStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+      refStars->wgtlist[Nref] = 1;
+      refStars->msklist[Nref] = 0;
+      if (fabs(refStars->flxlist[Nref]) > 0.03) {
         // fprintf (stderr, "deviant\n");
       }
 
-      if (image[i].imageID == TEST_IMAGE) {
-        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat, list[N]);
+      if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
+        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat, refStars->flxlist[Nref]);
       }
 
       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-        Mlist[Nbright] = list[N];
-        dMlist[Nbright] = dlist[N];
+        brightStars->flxlist[Nbright] = refStars->flxlist[Nref];
+        brightStars->errlist[Nbright] = refStars->errlist[Nref];
+        brightStars->wgtlist[Nbright] = 1;
+        brightStars->msklist[Nbright] = 0;
         Nbright ++;
       }
-      N++;
+      Nref++;
     }
     /* N_onImage[i] is all measurements, N is good measurements */
@@ -740 +750 @@
     mark = FALSE;
     if (!PoorImages) {
-      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onImage[i]);
+      mark = (Nref < IMAGE_TOOFEW) || (Nref < IMAGE_GOOD_FRACTION*N_onImage[i]);
       if (mark) {
         image[i].flags |= ID_IMAGE_PHOTOM_FEW;
@@ -750 +760 @@
 
     // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
-    liststats (list, dlist, NULL, N, &stats);
-    image[i].Mcal  = stats.mean;
-    image[i].dMcal = stats.error;
-    image[i].nFitPhotom = N;
-    image[i].Xm    = 100.0*log10(stats.chisq);
+# if (BASIC_STATS) 
+    liststats (refStars->flxlist, refStars->errlist, NULL, Nref, &stats);
+# else
+    liststats_irls (refStars, Nref, &stats);
+# endif
+    image[i].Mcal       = stats.mean;
+    image[i].dMcal      = stats.error;
+    image[i].nFitPhotom = Nref;
+    image[i].Xm         = 100.0*log10(stats.chisq);
     Ncalibrated ++;
 
-    if (image[i].imageID == TEST_IMAGE) {
+    if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
+      for (j = 0; j < Nref; j++) {
+        fprintf (stderr, "%1d, %8d : %6.3f  %6.3f  %6.3f  %d\n", (int) i, (int) image[i].imageID, refStars->flxlist[j], refStars->errlist[j], refStars->wgtlist[j], refStars->msklist[j]);
+      }
+    }
+
+    if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
       fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
     }
@@ -767 +787 @@
     if (PLOTSTUFF) {
       fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
-      plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
+      plot_setMcal (refStars->flxlist, N, &stats, CLOUD_TOLERANCE);
     }
 
     // bright end scatter
-    liststats (Mlist, dMlist, NULL, Nbright, &stats);
+    liststats (brightStars->flxlist, brightStars->errlist, NULL, Nbright, &stats);
     image[i].dMagSys = stats.sigma;
 
@@ -782 +802 @@
     image[i].ubercalDist = minUbercalDist + 1;
   }
-  free (list);
-  free (dlist);
-  free (Mlist);
-  free (dMlist);
+
+  StatDataSetFree (brightStars, 1);
+  StatDataSetFree (refStars, 1);
 
   fprintf (stderr, "%d images calibrated\n", Ncalibrated);

trunk/Ohana/src/relphot/src/args.c

@@ -128 +128 @@
   }
 
+  // elements needed for parallel regions / parallel images
+  MANUAL_UNIQUER = NULL;
+  if ((N = get_argument (argc, argv, "-manual-uniquer"))) {
+    remove_argument (N, &argc, argv);
+    MANUAL_UNIQUER = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   PLOTSTUFF = FALSE;
   if ((N = get_argument (argc, argv, "-plot"))) {
@@ -176 +184 @@
   }
 
-  TEST_IMAGE = -1;
-  if ((N = get_argument (argc, argv, "-test-image"))) {
-    remove_argument (N, &argc, argv);
-    TEST_IMAGE = atoi(argv[N]);
+  TEST_IMAGE1 = -1;
+  if ((N = get_argument (argc, argv, "-test-image1"))) {
+    char *endptr;
+    remove_argument (N, &argc, argv);
+    TEST_IMAGE1 = strtol(argv[N], &endptr, 0);
+    if (*endptr) relphot_usage (); 
+    remove_argument (N, &argc, argv);
+  }
+  TEST_IMAGE2 = -1;
+  if ((N = get_argument (argc, argv, "-test-image2"))) {
+    char *endptr;
+    remove_argument (N, &argc, argv);
+    TEST_IMAGE2 = strtol(argv[N], &endptr, 0);
+    if (*endptr) relphot_usage (); 
     remove_argument (N, &argc, argv);
   }
@@ -283 +301 @@
     FREEZE_IMAGES = TRUE;
   }
-
   FREEZE_MOSAICS = FALSE;
   if ((N = get_argument (argc, argv, "-mosfreeze"))) {
     remove_argument (N, &argc, argv);
     FREEZE_MOSAICS = TRUE;
+  }
+  CALIBRATE_STACKS_AND_WARPS = FALSE;
+  if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) {
+    remove_argument (N, &argc, argv);
+    CALIBRATE_STACKS_AND_WARPS = TRUE;
   }
 
@@ -620 +642 @@
   }
 
-  TEST_IMAGE = -1;
+  TEST_IMAGE1 = -1;
+  TEST_IMAGE2 = -1;
 
   RESET = FALSE;

trunk/Ohana/src/relphot/src/bcatalog.c

@@ -106 +106 @@
 
       // skip garbage measurements
-      if (catalog[0].measure[offset].psfQF < 0.85) { Npsfqf ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF < 0.85))     { Npsfqf ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.85)) { Npsfqf ++; continue; }
       if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
 

trunk/Ohana/src/relphot/src/liststats.c

@@ -210 +210 @@
   
   // XXX add to dataset elements?
-  ALLOCATE_PTR (wt, double, Npoints);
+  ALLOCATE_PTR (wtvals, double, Npoints);
+  ALLOCATE_PTR (wtlist, double, Npoints);
 
   int converged = FALSE;
@@ -217 +218 @@
     for (int i = 0; i < Npoints; i++) {
       // we are only including the formal error, not the weight in the definition of wt[]
-      wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+      wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
     }
     
     double oldValue = value;
-    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, wt, Npoints)) {
+    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, wtvals, Npoints)) {
       value = oldValue;
       break;
@@ -235 +236 @@
       
   // calculate the weight thresholds to mask the bad points:
-  double Sum_W = 0.0;
+  // double Sum_W = 0.0;
   for (int i = 0; i < Npoints; i++) {
-    wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
-    Sum_W += wt[i];
-  }
-  double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
+    wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+    wtlist[i] = wtvals[i];
+    // Sum_W += wtvals[i];
+  }
+  dsort (wtlist, Npoints);
+  int midpt = 0.5 * Npoints;
+  double WtMedian = (Npoints % 2) ? wtlist[midpt] : 0.5*(wtlist[midpt] + wtlist[midpt-1]);
+  // double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
+  double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
 
   // generate the unmasked subset
@@ -255 +261 @@
   double dSig = 0.0;
   for (int i = 0; i < Npoints; i++) {
-    if ((wt[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
+    if ((wtvals[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
       dataset->msklist[i] = TRUE; // mark the masked points
       continue;
@@ -313 +319 @@
   free (dykeep);
   free (wtkeep);
-  free (wt);
+  free (wtvals);
+  free (wtlist);
 
   return TRUE;

trunk/Ohana/src/relphot/src/load_catalogs.c

@@ -115 +115 @@
   int TIME = time(NULL);
   int PID = getpid();
-  snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000);
+
+  if (PARALLEL_MANUAL && MANUAL_UNIQUER) {
+    snprintf (uniquer, 12, "%11s", MANUAL_UNIQUER);
+  } else {
+    snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000);
+  }
 
   int Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog

trunk/Ohana/src/relphot/src/relphot_images.c

@@ -53 +53 @@
   FlatCorrectionTable *flatcorr = NULL;
 
-  if (NLOOP > 0) {
+  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
     /* load catalog data from region files (hostID is 0 since we are not a client */
     catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL);
@@ -154 +154 @@
 
     /* set Mcal & Mmos for bad images */
-    setMcal  (catalog, TRUE, flatcorr);
-    setMmos  (catalog, TRUE, flatcorr);
+    int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
+    setMcal  (catalog, onlyPoorImages, flatcorr);
+    setMmos  (catalog, onlyPoorImages, flatcorr);
     MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 
@@ -248 +249 @@
      context, this function writes the image parameters as a subset table for the remote
      clients */
-  reload_catalogs (skylist, flatcorr, 0, NULL);
-  MARKTIME("-- updated all catalogs: %f sec\n", dtime);
+  if (!CALIBRATE_STACKS_AND_WARPS) {
+    reload_catalogs (skylist, flatcorr, 0, NULL);
+    MARKTIME("-- updated all catalogs: %f sec\n", dtime);
+  }
 
   if (UPDATE_CATFORMAT) {
@@ -255 +258 @@
     db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
   }
-  if (NLOOP > 0) {
+  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
     // do not save changes if we did not make changes. 
     dvo_image_update (&db, VERBOSE);