Changeset 34041 for branches/meh_branches/ppstack_test/Ohana/src/opihi/dvo/avmatch.c

Timestamp:

Jun 19, 2012, 5:24:19 PM (14 years ago)

Author:

mhuber

Message:

merging latest r34040 trunk changes to branch

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• Ohana/src/opihi/dvo/avmatch.c (modified) (13 diffs)

branches/meh_branches/ppstack_test/Ohana/src/opihi/dvo/avmatch.c

@@ -14 +14 @@
   
   off_t i, j, n, m, *index;
-  int N, Ncat, Npts, NPTS, last, Nfields, Nsecfilt;
+  int N, Ncat, Npts, NPTS, last, Nfields, Nsecfilt, Ninvec;
   int VERBOSE;
-  char name[1024];
+  char name[1024], *found;
   void *Signal;
   float RADIUS;
@@ -22 +22 @@
   Catalog catalog;
 
-  Vector **vec, *RAvec, *DECvec;
+  Vector **vec, **invec, *RAvec, *DECvec;
   dbField *fields;
   dbValue *values;
@@ -29 +29 @@
   /* defaults */
   vec = NULL;
+  invec = NULL;
   fields = NULL;
   values = NULL;
   skylist = NULL;
+  Ninvec = 0;
 
   if ((N = get_argument (argc, argv, "-h"))) goto help;
   if ((N = get_argument (argc, argv, "--help"))) goto help;
-  if (argc < 5) goto help;
 
   VERBOSE = FALSE;
@@ -43 +44 @@
   }
 
+  int PARALLEL = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel"))) {
+    remove_argument (N, &argc, argv);
+    PARALLEL = TRUE;
+  }
+
+  // dump results directly to fits file (esp for parallel dvo)
+  char *CoordsFile = NULL;
+  if ((N = get_argument (argc, argv, "-coords"))) {
+    remove_argument (N, &argc, argv);
+    CoordsFile = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (!CoordsFile && (argc < 5)) goto help;
+  if ( CoordsFile && (argc < 3)) goto help;
+
   dvo_catalog_init (&catalog, TRUE);
 
@@ -49 +67 @@
   Nsecfilt = GetPhotcodeNsecfilt ();
 
+  // init locally static variables (time refs)
+  dbExtractAveragesInit (); 
+
+  // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results
+  RAvec  = NULL;
+  DECvec = NULL;
+  if (PARALLEL && !HOST_ID) {
+    if (!CoordsFile) {
+      // get vectors corresponding to coordinates of interest
+      if ((RAvec  = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help;
+      if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help;
+      
+      ALLOCATE (vec, Vector *, 2);
+      vec[0] = RAvec;
+      vec[1] = DECvec;
+
+      CoordsFile = abspath("coords.fits", 1024);
+      int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, NULL);
+      if (!status) goto escape;
+    }
+
+    char *targv1 = argv[1];
+    char *targv2 = argv[2];
+    argv[1] = strcreate ("-coords");
+    argv[2] = strcreate (CoordsFile);
+    free (CoordsFile);
+
+    // I need to pass the RA & DEC vectors to the remote clients...
+    int status = HostTableParallelOps (argc, argv, RESULT_FILE, RAvec->Nelements, VERBOSE);
+    if (vec) free (vec);
+    
+    free (argv[1]);
+    free (argv[2]);
+    argv[1] = targv1;
+    argv[2] = targv2;
+
+    return status;
+  }
+
   // get vectors corresponding to coordinates of interest
-  if ((RAvec  = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help;
-  if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help;
-  RADIUS = atof (argv[3]);
-
-  // strip off RA, DEC, RADIUS arguments
-  remove_argument (1, &argc, argv);
-  remove_argument (1, &argc, argv);
+  if (CoordsFile) {
+    // read RAvec, DECvec from coords file (1st 2 fields?)
+    Ninvec = 0;
+    invec = ReadVectorTableFITS (CoordsFile, "COORDS", &Ninvec);
+    RAvec = invec[0];
+    DECvec = invec[1];
+  } else {
+    if ((RAvec  = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help;
+    if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help;
+    // strip off RA & DEC
+    remove_argument (1, &argc, argv);
+    remove_argument (1, &argc, argv);
+  }
+  RADIUS = atof (argv[1]);
   remove_argument (1, &argc, argv);
 
@@ -66 +130 @@
     dvo_catalog_free (&catalog);
     goto help;
+  }
+
+  // check the requested fields : are all average/secfilt entries, or do we need measures?
+  int needMeasures = FALSE;
+  for (i = 0; !needMeasures && (i < Nfields); i++) {
+    if (fields[i].magMode == MAG_NONE) continue;
+    if (fields[i].photcode == NULL) continue; // assert this?
+    if (fields[i].photcode[0].type == PHOT_REF) needMeasures = TRUE;
+    if (fields[i].photcode[0].type == PHOT_DEP) needMeasures = TRUE;
   }
 
@@ -92 +165 @@
   }
   ALLOCATE (index, off_t, NPTS);
+  ALLOCATE (found, char, NPTS);
+  memset (found, 0, NPTS*sizeof(char));
 
   // grab data from all selected sky regions
@@ -97 +172 @@
   interrupt = FALSE;
   for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) {
+
+    // does this host ID match the desired location for the table?
+    if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue;
+
     /* lock, load, unlock catalog */
-    catalog.filename = skylist[0].filename[i];
-    catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
+    char hostfile[1024];
+    snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name);
+    catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
+    catalog.catflags = LOAD_AVES | LOAD_SECF;
+    catalog.catflags |= needMeasures ? LOAD_MEAS : SKIP_MEAS;
     catalog.Nsecfilt = 0;
 
@@ -124 +206 @@
       if (Ncat == -2) continue;
 
-      // XXX set a 'found' vector to double check we catch everything?
       m = catalog.average[Ncat].measureOffset;
 
       // reset counters for saved fields, extract fields
-      dbExtractAveragesInit (); 
+      dbExtractAveragesInitAve (); 
       for (n = 0; n < Nfields; n++) {
         values[n] = dbExtractAverages (&catalog.average[Ncat], &catalog.secfilt[Ncat*Nsecfilt], &catalog.measure[m], &fields[n]);
@@ -141 +222 @@
         }
       }
+      found[Npts] = TRUE;
     }
     dvo_catalog_free (&catalog);
@@ -147 +229 @@
   interrupt = FALSE;
 
+  // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
+  // only write the fields which were in a valid catalog
+  if (RESULT_FILE) {
+    // extend the array by one to hold index array
+    REALLOCATE (vec, Vector *, Nfields + 1);
+    vec[Nfields] = InitVector();
+    strcpy (vec[Nfields]->name, "index");
+    ResetVector (vec[Nfields], OPIHI_INT, NPTS);
+    Vector *idxVec = vec[Nfields];
+
+    // only write out the rows which were found
+    Npts = 0;
+    for (i = 0; i < NPTS; i++) {
+      if (!found[i]) continue;
+      idxVec->elements.Int[Npts] = i;
+      Npts ++;
+    }
+    int Nfound = Npts;
+    idxVec->Nelements = Nfound;
+
+    fprintf (stderr, "found %d of %d pts\n", Nfound, NPTS);
+
+    for (i = 0; i < Nfields; i++) {
+      if (vec[i][0].type == OPIHI_FLT) {
+        opihi_flt *tmp = NULL;
+        ALLOCATE (tmp, opihi_flt, Nfound);
+        Npts = 0;
+        for (j = 0; j < NPTS; j++) {
+          if (!found[j]) continue;
+          tmp[Npts] = vec[i][0].elements.Flt[j];
+          Npts++;
+        }
+        free (vec[i][0].elements.Flt);
+        vec[i][0].elements.Flt = tmp;
+      } else {
+        opihi_int *tmp = NULL;
+        ALLOCATE (tmp, opihi_int, Nfound);
+        Npts = 0;
+        for (j = 0; j < NPTS; j++) {
+          if (!found[j]) continue;
+          tmp[Npts] = vec[i][0].elements.Int[j];
+          Npts++;
+        }
+        free (vec[i][0].elements.Int);
+        vec[i][0].elements.Int = tmp;
+      }
+      vec[i][0].Nelements = Nfound;
+    }
+    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields + 1, FALSE, NULL);
+    free (vec[Nfields]->elements.Int);
+    free (vec[Nfields]);
+    if (!status) goto escape;
+  }
+
   if (vec) free (vec);
   if (values) free (values);
+  if (invec) FreeVectorArray (invec, Ninvec);
   dbFreeFields (fields, Nfields);
   SkyListFree (skylist);
@@ -156 +293 @@
   if (vec) free (vec);
   if (values) free (values);
+  if (invec) FreeVectorArray (invec, Ninvec);
   dbFreeFields (fields, Nfields);
   SkyListFree (skylist);
@@ -162 +300 @@
  help:
   gprint (GP_ERR, "USAGE: avmatch (RA) (DEC) (RADIUS) field[,field,field...]\n");
+  gprint (GP_ERR, "   OR: avmatch -coords (filename.fits) (RADIUS) field[,field,field...]\n");
 
   if ((argc > N + 1) && !strcasecmp (argv[N+1], "fields")) {