Changeset 36516

Timestamp:

Feb 15, 2014, 1:15:44 PM (12 years ago)

Author:

eugene

Message:

working on the parallel relphot

Location:

branches/eam_branches/ipp-20140206/Ohana/src

Files:

• libdvo/include/dvo.h (modified) (1 diff)
• libdvo/src/RegionHostTable.c (modified) (3 diffs)
• relphot/doc/parallel.txt (modified) (1 diff)
• relphot/include/relphot.h (modified) (1 diff)
• relphot/src/ImageOps.c (modified) (3 diffs)
• relphot/src/ImageTable.c (added)
• relphot/src/MeanMagIO.c (added)
• relphot/src/MosaicOps.c (modified) (2 diffs)
• relphot/src/StarOps.c (modified) (1 diff)
• relphot/src/args.c (modified) (2 diffs)
• relphot/src/assign_images.c (modified) (5 diffs)
• relphot/src/help.c (modified) (1 diff)
• relphot/src/indexCatalog.c (added)
• relphot/src/launch_region_hosts.c (added)
• relphot/src/load_catalogs.c (modified) (2 diffs)
• relphot/src/relphot_parallel_images.c (modified) (9 diffs)
• relphot/src/relphot_parallel_regions.c (modified) (1 diff)
• relphot/src/share_image_mags.c (added)
• relphot/src/share_mean_mags.c (modified) (2 diffs)
• uniphot/src/update_dvo_setphot.c (modified) (1 diff)

branches/eam_branches/ipp-20140206/Ohana/src/libdvo/include/dvo.h

@@ -295 +295 @@
   int pid;                    // remote process ID
   char *hostname;
-  ImageSubset *image;
+
   int line_numbers;
+
+  off_t Nimage;
+  off_t NIMAGE;
+  Image *image;
 } RegionHostInfo;
 

branches/eam_branches/ipp-20140206/Ohana/src/libdvo/src/RegionHostTable.c

@@ -7 +7 @@
   int i;
   for (i = Nhosts; i < NHOSTS; i++) {
+    hosts[i].Rmin = NAN;
+    hosts[i].Rmax = NAN;
+    hosts[i].Dmin = NAN;
+    hosts[i].Dmax = NAN;
+
     hosts[i].hostname = NULL;
     hosts[i].stdio[HOST_STDIN] = -1;
@@ -13 +18 @@
     hosts[i].pid = 0;
 
-    hosts[i].Rmin = NAN;
-    hosts[i].Rmax = NAN;
-    hosts[i].Dmin = NAN;
-    hosts[i].Dmax = NAN;
+    hosts[i].image = NULL;
+    hosts[i].Nimage = 0;
+    hosts[i].NIMAGE = 0;
   }
   return;
@@ -143 +147 @@
   return table;
 }
-

branches/eam_branches/ipp-20140206/Ohana/src/relphot/doc/parallel.txt

@@ +1 @@
+
+2014.02.14
+
+ more relphot notes:
+
+ * I load the catalogs, with a  
+
+2014.02.12
+
+I am making progess on the relphot -parallel-regions implementation.  some things I need to deal with:
+
+ * assign_images: select mosaics for all images saved for a given host (need to supply them as well)
+   ** I need to decide on the Image / Mosaic split.  If I am calibrating by image, then the image center defines ownership
+   ** if I am calibrating by exposure, then the mosaic center must define ownership
+ o determine SkyList covering the images for region host
+ o tag detections and objects which I own or do not own
+   o ImageOps.c / matchImage matches detections to images and needs to handle mine / not mine cases
+   o ditto for matchMosaics
+ o load_catalogs : I need to make my requested catalogs unique (just add uniquer to load_catalogs.c:130
+ o figure out what fields I'm saving in meanmags and how to construct it
+   -- mag, objID, catID, photcode, (dmag or other stats?)
+ * from which hosts do I slurp mean mags?
+ o function to merge a new meanmag array into the existing one
+ * how to go from (objID,catID) to a given catalog[i].average[j]
 
 2014.02.06

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

@@ -110 +110 @@
   StatType kronstats;
 } SetMrelInfo;
+
+typedef struct {
+  float M;
+  float dM;
+  float Xm;
+  int Nsec;
+  unsigned int objID;
+  unsigned int catID;
+} MeanMag;
 
 typedef struct {

branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/ImageOps.c

@@ -244 +244 @@
   
   int Nmatch = 0;
- for (i = 0; i < Ncatalog; i++) {
+  for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Nmeasure; j++) {
+      catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise
       ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode);
+
+      // skip measurements which do not match one of the requested photcodes (
+      // (do we not already exclude in bcatalog -- maybe needed for reload_objects?
       found = FALSE;
       for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
@@ -252 +256 @@
       }
       if (!found) continue; 
+
+      // if we match one of our images, myDet gets set to TRUE
       matchImage (catalog, j, i, doImageList);
       Nmatch ++;
@@ -322 +328 @@
     return;
   }
+  catalog[cat].measureT[meas].myDet = TRUE;
 
   if (USE_GRID) {

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

@@ -527 +527 @@
 
   if (!MOSAIC_ZEROPT) return (FALSE);
+  // if we are calibrating by mosaic, redefine myDet == on one of my mosaics
 
   int Nmatch = 0;
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Nmeasure; j++) {
+      catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise
+
       if (TimeSelect) {
         if (catalog[i].measureT[j].t < TSTART) continue;
@@ -582 +585 @@
   }
 
+  // this measurement is on one of my mosaics, mark it as mine.
+  catalog[cat].measureT[meas].myDet = TRUE;
   MeasureToMosaic[cat][meas] = mosID;
 

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

@@ -896 +896 @@
 }
 
+int markObjects (Catalog *catalog, int Ncatalog) {
+
+  // How strongly do I own this object?
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      int nOwn = 0;
+      int m = catalog[i].average[j].measureOffset;
+      for (n = 0; n < catalog[i].average[j].Nmeasure; n++) {
+        if (!catalog[i].measureT[m+n].myDet) continue;
+        nOwn ++;
+      }
+      catalog[i].averageT[j].nOwn = nOwn;
+    }
+  }
+  return TRUE;
+}
+
 void clean_stars (Catalog *catalog, int Ncatalog) {
 

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

@@ -343 +343 @@
   }
 
+  if ((N = get_argument (argc, argv, "-region-hosts"))) {
+    remove_argument (N, &argc, argv);
+    REGION_FILE = strcreate (argv[N])
+    remove_argument (N, &argc, argv);
+  }
+
+  REGION_HOST_ID = 0;
+  if ((N = get_argument (argc, argv, "-region-hostID"))) {
+    remove_argument (N, &argc, argv);
+    REGION_HOST_ID = atoi (argv[N])
+    remove_argument (N, &argc, argv);
+  }
+
   RelphotMode mode = MODE_ERROR;
   if ((N = get_argument (argc, argv, "-averages"))) {
@@ -355 +368 @@
     remove_argument (N, &argc, argv);
     mode = PARALLEL_IMAGES;
+    if (N >= argc) relphot_usage();
+    IMAGE_TABLE = strcreate (argv[N])
+    remove_argument (N, &argc, argv);
+    if (!REGION_FILE) relphot_usage();
   }
   if ((N = get_argument (argc, argv, "-parallel-regions"))) {
     remove_argument (N, &argc, argv);
     mode = PARALLEL_REGIONS;
-    if (N >= argc) relphot_usage();
-    REGION_FILE = strcreate (argv[N])
-    remove_argument (N, &argc, argv);
+    if (!REGION_FILE) relphot_usage();
   }
 

branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/assign_images.c

@@ -25 +25 @@
 
   // for each regionHost, select images which are contained by the region
-  // XXX this is probably too expensive : I need to loop over Nimage * Nregions
-  // inverting the order would probably be faster (need to track NIMAGE for each host)
   // even faster would be to use a tree to get to the real regions...
   select_images_hostregion (regionHosts, image, Nimage);
@@ -85 +83 @@
     i = find_host_for_coords (regionHosts, Rc, Dc);
 
+    // this is a bit memory expensive : I am making a complete copy of the image table here
     off_t Nsubset = regionHosts->host[i].Nimage;
-    regionHosts->host[i].image[Nsubset] = image[j]; 
-    regionHosts->host[i].line_number[Nsubset] = j;
+    regionHosts->host[i].image[Nsubset] = image[j];
+    // regionHosts->host[i].line_number[Nsubset] = j;
 
     regionHosts->host[i].Nimage ++;
@@ -93 +92 @@
       regionHosts->host[i].NIMAGE += D_NIMAGE;
       REALLOCATE (regionHosts->host[i].image, Image, regionHosts->host[i].NIMAGE);
-      REALLOCATE (regionHosts->host[i].line_number, off_t, regionHosts->host[i].NIMAGE);
+      // REALLOCATE (regionHosts->host[i].line_number, off_t, regionHosts->host[i].NIMAGE);
     }
   }
@@ -99 +98 @@
 }
 
+// XXX add a search tree to speed this up?
 int find_host_for_coords (RegionHostTable *regionHosts, double R, double d) {
 
@@ -112 +112 @@
   }
 }
+

branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/help.c

@@ -6 +6 @@
   fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
   fprintf (stderr, "       or:    relphot -apply-offsets -region RA RA DEC DEC\n");
-  fprintf (stderr, "       or:    relphot -parallel-regions (RegionFile)\n");
-  fprintf (stderr, "       or:    relphot -parallel-images ()\n");
+  fprintf (stderr, "       or:    relphot -parallel-regions -region-hosts (RegionFile)\n");
+  fprintf (stderr, "       or:    relphot -parallel-images (ImageTable) -region-hosts (RegionFile)\n");
   fprintf (stderr, "  use -h for more usage information\n");
   exit (2);

branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/load_catalogs.c

@@ -108 +108 @@
 Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) {
 
+  char uniquer[12];
+  int TIME = time(NULL);
+  int PID = getpid();
+  snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000);
+
   int Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
 
@@ -128 +133 @@
 
     ALLOCATE (table->hosts[i].results, char, 1024);
-    snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.subset.dat", table->hosts[i].pathname);
+    snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.%s.dat", table->hosts[i].pathname, uniquer);
 
     // options / arguments that can affect relphot_client -load:

branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/relphot_parallel_images.c

@@ -19 +19 @@
   INITTIME;
 
+  // load the RegionTable (UserRegion should not be used at this level)
+  RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE);
+
   // load the subset images belonging to this host
-  // XXX need to determine the file name for each host (here or in master)
   off_t Nimage;
-  ImageSubset *image = ImageSubsetLoad (IMAGES, &Nimage);
+  Image *image = ImageTableLoad (IMAGE_TABLE, &Nimage);
   if (!image) {
-    fprintf (stderr, "ERROR loading image subset %s\n", CATDIR);
+    fprintf (stderr, "ERROR loading image %s\n", IMAGE_TABLE);
     exit (2);
   }
-  client_logger_message ("loaded image subset data\n");
 
   // load the flat correction table (if defined)
@@ -34 +35 @@
   FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
 
-  // XXX need to define the skylist for this region
-  skylist = function(region);
+  SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (sky, CATDIR, "cpt");
+  SkyList *skylist = SkyListByBounds (sky, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
 
   /* load catalog data from region files (hostID is 0 since we are not a client */
@@ -41 +43 @@
   MARKTIME("-- load catalog data: %f sec\n", dtime);
   
+  // generate tables go from catID,objID -> catSeq,objSeq
+  indexCatalogs (catalog, Ncatalog);
+
   /* match measurements with images, mosaics */
   initImageBins  (catalog, Ncatalog, TRUE);
@@ -46 +51 @@
 
   initMosaicBins (catalog, Ncatalog, TRUE);
-  initGridBins   (catalog, Ncatalog);
   initMrel (catalog, Ncatalog);
 
@@ -54 +58 @@
   findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
   MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
+
+  markObjects (catalog, Ncatalog);
 
   SAVEPLOT = FALSE;
@@ -76 +82 @@
 
     // share mean mags for objects at the boundary (number of unowned meas > 0)
-    share_mean_mags (catalog, Ncatalog);
+    share_mean_mags (catalog, Ncatalog, regionHosts, i);
 
     // load mean mags from other region hosts
-    slurp_mean_mags (catalog, Ncatalog);
+    slurp_mean_mags (catalog, Ncatalog, regionHosts, i);
 
     // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements
@@ -87 +93 @@
     
     // share image mags for images with non-zero unowned detections
-    share_image_mags ();
+    share_image_mags (regionHosts, i);
 
-    slurp_image_mags ();
-
-    if (PLOTSTUFF) {
-      plot_scatter (catalog, Ncatalog, flatcorr); 
-      plot_mosaics ();
-      plot_images ();
-      plot_stars (catalog, Ncatalog);
-      plot_chisq (catalog, Ncatalog);
-    }
+    slurp_image_mags (regionHosts, i);
 
     if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); 
@@ -108 +106 @@
   }
 
-  if (PLOTSTUFF) {
-    plot_scatter (catalog, Ncatalog, flatcorr); 
-    plot_grid (catalog, flatcorr); 
-    plot_mosaics ();
-    plot_images ();
-    plot_stars (catalog, Ncatalog);
-    plot_chisq (catalog, Ncatalog);
-  }
-  
   /* set Mcal & Mmos for bad images */
   setMcal  (catalog, TRUE, flatcorr);
@@ -122 +111 @@
   MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 
-  setMcalFinal (); // copy per-mosaic calibrations to the images
+  share_image_mags (regionHosts, i);
 
-  /* at this point, we have correct cal coeffs in the image/mosaic structures */
-  for (i = 0; i < Ncatalog; i++) {
-    // these tiny values are set by BrightCatalogSplit from load_catalogs
-    free_tiny_values (&catalog[i]);
-    dvo_catalog_free (&catalog[i]);
-  }
-  freeImageBins (Ncatalog, TRUE);
-  freeMosaicBins (Ncatalog, TRUE);
-  freeGridBins (Ncatalog);
-
-  save_image_subset ();
+  // save_image_subset ();
 
   exit (0);
-
 }

branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/relphot_parallel_regions.c

@@ -42 +42 @@
 
   /* launch processing on the parallel region hosts */
-  launch_region_hosts (regionHosts, XXX);
+  launch_region_hosts (regionHosts);
 
-  /* retrieve updated image parameters from the remote hosts */
-  retrieve_images (regionHosts, XXX);
+  // retrieve updated image parameters from the remote hosts (also set Image.mcal)
+  slurp_image_mags (regionHosts, XXX);
 
   /* update catalogs (in parallel) */

branches/eam_branches/ipp-20140206/Ohana/src/relphot/src/share_mean_mags.c

@@ -3 +3 @@
 // we are sharing mean mags for all objects which (a) I own and (b) which have unowned detections
 
-int share_mean_mags () {
+# define D_NMEANMAGS 10000
+int share_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) {
 
-  SomeType *meanmags = NULL;
-  ALLOCATE (meanmags, SomeType, NMEANMAGS);
+  off_t Nmeanmags = 0;
+  off_t NMEANMAGS = D_NMEANMAGS;
 
-  // XXX skip some internal catalogs?
+  MeanMag *meanmags = NULL;
+  ALLOCATE (meanmags, MeanMag, NMEANMAGS);
+
+  int Ns;
+  int Nsecfilt = GetPhotcodeNsecfilt();
+
+  // XXX skip some catalogs based on UserPatch?
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
 
       // do I own this object? (in region range?)
-      if (catalog[i].average[j].R <  Rmin) continue;
-      if (catalog[i].average[j].R >= Rmax) continue;
-      if (catalog[i].average[j].D <  Dmin) continue;
-      if (catalog[i].average[j].D >= Dmax) continue;
+      if (catalog[i].average[j].R <  UserPatch.Rmin) continue;
+      if (catalog[i].average[j].R >= UserPatch.Rmax) continue;
+      if (catalog[i].average[j].D <  UserPatch.Dmin) continue;
+      if (catalog[i].average[j].D >= UserPatch.Dmax) continue;
 
       // does this object have missing detections (does someone else need it?)
-      if (catalog[i].averageExtra[j].Nmiss == 0) continue;
+      // XXX : sky objects without missing detections
+      // XXX watch out for detections which are not associated with an image (REF)
+      if (catalog[i].averageT[j].nOwn != catalog[i].averageT[j].Nmeasure) continue;
 
-      meanmags[Nmeanmags] = something(catalog[i].average[j]);
-      Nmeanmags ++;
-
-      // check realloc..
+      for (Ns = 0; Ns < Nphotcodes; Ns++) {
+        int thisCode = photcodes[Ns][0].code;
+        int Nsec = GetPhotcodeNsec(thisCode);
+        set_mean_mags (&meanmags[Nmeanmags], &catalog[i].average[j], &catalog[i].secfilt[Nsecfilt*j + Nsec], Nsec);
+        Nmeanmags ++;
+        CHECK_REALLOCATE (meanmags, MeanMag, NMEANMAGS, Nmeanmags, D_NMEANMAGS);
+      }
     }
   }
 
   // write out the meanmag fits table AND write state in some file
-  SaveMeanMags (meanmags, Nmeanmags, where);
+  int myHost = regionHosts->index[REGION_HOST_ID];
+  char *hostname = regionHosts->hosts[myHost].hostname;
+
+  char filename[1024];
+  snprintf (filename, 1024, "%s/%s.meanmags.fits", CATDIR, hostname);
+  MeanMagSave (filename, meanmags, Nmeanmags);
+
+  update_sync_file (hostname, nloop);
+
+  return TRUE;
 }
 
-int slurp_mean_mags () {
+int slurp_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) {
 
-  SomeType *meanmags = NULL;
-  ALLOCATE (meanmags, SomeType, NMEANMAGS);
+  int Nmeanmags = 0;
+  MeanMag *meanmags = NULL;
+  ALLOCATE (meanmags, MeanMag, 1);
 
   for (i = 0; i < Nhost; i++) {
     if (not_neighbor(host[i])) continue;
+    // XXX skip regions->hosts[i].hostID == REGION_HOST_ID
 
-    check_sync_file (host[i]);
+    check_sync_file (regionHosts->hosts[i].hostname, nloop);
     
-    meanmagsSubset = LoadMeanMags (host[i], &Nsubset);
+    char filename[1024];
+    snprintf (filename, 1024, "%s/%s.meanmags.fits", CATDIR, regionHosts->hosts[i].hostname);
 
-    merge_mean_mags ();
+    int Nsubset = 0;
+    MeanMag *meanmagsSubset = MeanMagLoad (host[i], &Nsubset);
+
+    // merge_mean_mags reallocs meanmags and frees the input meanmagsSubset
+    meanmags = merge_mean_mags (meanmags, &Nmeanmags, meanmagsSubset, Nsubset);
   }
 
@@ -51 +79 @@
     catID = meanmags[i].catID;
 
+    // set the mean mag
+    // XXX I need to figure out how to go from catID,objID to catSeq,objSeq
     catSeq = catalog_ID_to_seq (catID);
-    objSeq = object_ID_to_seq (objID);
+    objSeq = object_ID_to_seq (catID, objID);
 
-    // set the mean mag
     Nsecfit = photcode_to_secfilt (meanmags[i].photcode);
     catalog[catSeq].secfilt[objSeq*Nsec + Nsecfilt].mag = meanmags[i].mag;
   }
 }
+
+# define MSG_LENGTH 11
+
+int check_sync_file (char *hostname, int nloop) {
+
+  char filename[1024];
+  sprintf (filename, "%s/%s.meanmags.sync", CATDIR, hostname);
+
+  // XXX I need a rule to generate the filename for each host
+  FILE *f = NULL; 
+
+  while (retry) {
+
+    f = fopen (filename);
+    if (!f) {
+      usleep (2000000);
+      continue;
+    }
+
+    int Nread = fread (message, MSG_LENGTH, 1, f);
+    if (Nread < MSG_LENGTH) {
+      fclose (f);
+      usleep (2000000);
+      continue;
+    }
+    fclose (f);
+
+    // message is of the form: NLOOP: %03d
+    sscanf (message, "%s %d", &loop);
+    if (loop != nloop) {
+      usleep (2000000);
+      continue;
+    }
+    return TRUE;
+  }
+}
+
+int update_sync_file (char *hostname, int nloop) {
+
+  char filename[1024];
+  sprintf (filename, "%s/%s.meanmags.sync", CATDIR, hostname);
+
+  FILE *f = fopen (filename, "w");
+  snprintf (message, 11, "NLOOP: %03d\n", nloop);
+  
+  int Nwrite = fwrite (message, MSG_LENGTH, 1, f);
+  if (Nwrite != MSG_LENGTH) {
+    fprintf (stderr, "failure to write sync message\n");
+    exit (3);
+  }
+
+  fclose (f);
+  return TRUE;
+}
+
+int set_mean_mags (MeanMag *meanmags, Average *average, SecFilt *secfilt, int Nsec) {
+
+  meanmags->M  = secfilt->M;
+  meanmags->dM = secfilt->dM;
+  meanmags->Xm = secfilt->Xm;
+  meanmags->Nsec = Nsec; // key to secfilt entry
+
+  meanmags->objID = average->objID;
+  meanmags->catID = average->catID;
+
+  return TRUE;
+}
+
+MeanMag *merge_mean_mags (MeanMag *target, int *ntarget, MeanMag *source, int Nsource) {
+
+  REALLOCATE (target, MeanMag, *ntarget + Nsource);
+  for (i = 0; i < Nsource; i++) {
+    int n = i + *ntarget;
+    target[n] = source[i];
+  }
+  
+  free (source);
+
+  *ntarget += Nsource;
+  return (target);
+}
+

branches/eam_branches/ipp-20140206/Ohana/src/uniphot/src/update_dvo_setphot.c

@@ -2 +2 @@
 int update_dvo_setphot_parallel (SkyTable *sky, Image *image, off_t Nimage);
 
+// XXX convert this to use the ImageSubset to save space
 int update_dvo_setphot (Image *image, off_t Nimage, FlatCorrectionTable *flatcorr) {