Index: /trunk/Ohana/src/opihi/dvo/mmatch.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/mmatch.c	(revision 34462)
+++ /trunk/Ohana/src/opihi/dvo/mmatch.c	(revision 34462)
@@ -0,0 +1,314 @@
+# include "dvoshell.h"
+
+/* This function uses the 'find_match' algorithm to select the objects of interest.
+   All measurements for matched objects result in an element in the output vectors.
+   if -index (index) is supplied, the index of the input vectors are stored to match
+   the output vectors.
+ */
+
+int mmatch (int argc, char **argv) {
+  
+  off_t i, j, k, n, m, *index;
+  int N, Ncat, Npts, NPTS, last, Nfields, Nsecfilt, Ninvec;
+  int VERBOSE;
+  char name[1024];
+  void *Signal;
+  float RADIUS;
+
+  Catalog catalog;
+
+  Vector **vec, **invec, *RAvec, *DECvec, *IDXvec;
+  dbField *fields;
+  dbValue *values;
+  SkyList *skylist;
+
+  /* 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;
+
+  VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
+
+  IDXvec = NULL;
+  if ((N = get_argument (argc, argv, "-index"))) {
+    remove_argument (N, &argc, argv);
+    IDXvec = SelectVector (argv[N], ANYVECTOR, TRUE);
+    if (IDXvec == NULL) goto help;
+    remove_argument (N, &argc, argv);
+  }
+
+  int PARALLEL = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel"))) {
+    remove_argument (N, &argc, argv);
+    PARALLEL = TRUE;
+  }
+
+  // read RA,DEC coords from a 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);
+
+  /* load photcode information */
+  if (!InitPhotcodes ()) goto escape;
+  Nsecfilt = GetPhotcodeNsecfilt ();
+
+  // init locally static variables (time refs)
+  dbExtractMeasuresInit();
+
+  // 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, 0, 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 (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);
+
+  // parse the fields to be extracted and returned
+  fields = dbCmdlineFields (argc, argv, DVO_TABLE_MEASURE, &last, &Nfields);
+  if (fields == NULL) goto help;
+  if ((Nfields == 0) || (last != argc)) {
+    dbFreeFields (fields, Nfields);
+    dvo_catalog_free (&catalog);
+    goto help;
+  }
+
+  // load image data if needed (for fields listed below)
+  int loadImages = FALSE;
+  int mosaicMode = FALSE;
+  for (i = 0; !loadImages && (i < Nfields); i++) {
+    if (!MEASURE_HAS_XCCD) {
+      // I'm keeping this code because it gives a way of handling dvo dbs that don't have
+      // measure.xccd if we need it
+      if (fields[i].ID == MEAS_XCCD) loadImages = TRUE;
+      if (fields[i].ID == MEAS_YCCD) loadImages = TRUE;
+    }
+    if (fields[i].ID == MEAS_XMOSAIC) 	    loadImages = mosaicMode = TRUE;
+    if (fields[i].ID == MEAS_YMOSAIC) 	    loadImages = mosaicMode = TRUE;
+    if (fields[i].ID == MEAS_EXTERN_ID)     loadImages = mosaicMode = TRUE;
+    if (fields[i].ID == MEAS_FLAT)          loadImages = mosaicMode = TRUE;
+    if (fields[i].ID == MEAS_CENTER_OFFSET) loadImages = mosaicMode = TRUE;
+  }
+
+  // use the whole sky (since we select random points around the sky)
+  SkyRegionSelection selection;
+  selection.useDisplay = FALSE;
+  selection.useSkyregion = FALSE;
+  if (loadImages && !SetImageSelection (mosaicMode, &selection)) goto escape;
+
+  /* load regions which contain all supplied RA,DEC coordinates */
+  if ((skylist = SelectRegionsByCoordVectors (RAvec, DECvec)) == NULL) goto escape;
+
+  /* create output storage vectors */
+  Npts = 0;
+  NPTS = 1000;
+  ALLOCATE (values, dbValue, Nfields);
+  ALLOCATE (vec, Vector *, Nfields);
+  for (i = 0; i < Nfields; i++) {
+    if (ISNUM(fields[i].name[0])) {
+      sprintf (name, "v_%s", fields[i].name);
+    } else {
+      sprintf (name, "%s", fields[i].name);
+    }
+    if ((vec[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) goto escape;
+    ResetVector (vec[i], fields[i].type, NPTS);
+  }
+  if (IDXvec) {
+    ResetVector (IDXvec, OPIHI_INT, NPTS);
+  }
+
+  // vectors to track recovered values
+  int Nelem = RAvec->Nelements;
+  ALLOCATE (index, off_t, Nelem);
+
+  // grab data from all selected sky regions
+  Signal = signal (SIGINT, handle_interrupt);
+  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 */
+    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 | LOAD_MEAS;
+    catalog.Nsecfilt = Nsecfilt;
+
+    if (VERBOSE) gprint (GP_ERR, "trying %s ("OFF_T_FMT" of "OFF_T_FMT")\n", catalog.filename,  i,  skylist[0].Nregions);
+      
+    // an error exit status here is a significant error
+    if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) {
+      gprint (GP_ERR, "ERROR: failure to open catalog file %s\n", catalog.filename);
+      exit (2);
+    }
+    dvo_catalog_unlock (&catalog);  /// we can unlock here since this is read-only (do not block other access)
+    if (catalog.Naverage == 0) {
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+
+    // returns the matches to AVERAGE objects; what do we do for each MEASURE?
+    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
+
+    for (j = 0; (j < Nelem) && !interrupt; j++) {
+      Ncat = index[j];
+
+      if (Ncat == -1) continue;
+      if (Ncat == -2) continue;
+
+      m = catalog.average[Ncat].measureOffset;
+
+      dbExtractMeasuresInitAve (); // reset counters for saved fields (costs very little)
+
+      for (k = 0; (k < catalog.average[Ncat].Nmeasure); k++, m++) {
+	if (catalog.measure[m].averef != Ncat) {
+	  gprint (GP_ERR, "ERROR: inconsistent measure->average link.  Unsorted database?\n");
+	  goto escape;
+	}
+
+	// extract the relevant values for this measurement
+	dbExtractMeasuresInitMeas (); // reset counters for saved fields  (costs very little
+	for (n = 0; n < Nfields; n++) {
+	  values[n] = dbExtractMeasures (&catalog.average[Ncat], &catalog.secfilt[Ncat*Nsecfilt], &catalog.measure[m], &fields[n]);
+	}
+
+	// XXX if we are allowed to return more rows than the supplied RA,DEC we will need to create an output RA,DEC
+	for (n = 0; n < Nfields; n++) {
+	  if (vec[n][0].type == OPIHI_FLT) {
+	    vec[n][0].elements.Flt[Npts] = values[n].Flt;
+	  } else {
+	    vec[n][0].elements.Int[Npts] = values[n].Int;
+	  }
+	}
+	if (IDXvec) {
+	  IDXvec[0].elements.Int[Npts] = j;
+	}
+
+	// extend length of output vectors
+	Npts++;
+	if (Npts >= NPTS) {
+	  NPTS += 2000;
+	  for (n = 0; n < Nfields; n++) {
+	    REALLOCATE (vec[n][0].elements.Flt, opihi_flt, NPTS);
+	  }
+	  if (IDXvec) {
+	    REALLOCATE (IDXvec[0].elements.Int, opihi_int, NPTS);
+	  }
+	}
+      }
+    }
+
+    dvo_catalog_free (&catalog);
+  }
+  signal (SIGINT, Signal);
+  interrupt = FALSE;
+
+  for (n = 0; n < Nfields; n++) {
+    ResetVector (vec[n], fields[n].type, Npts);
+  }
+  if (IDXvec) {
+    ResetVector (IDXvec, IDXvec->type, Npts);
+  }
+
+  // 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] = IDXvec;
+
+    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields + 1, FALSE, NULL);
+    if (!status) goto escape;
+  }
+
+  if (vec) free (vec);
+  if (values) free (values);
+  if (invec) FreeVectorArray (invec, Ninvec);
+  dbFreeFields (fields, Nfields);
+  SkyListFree (skylist);
+  return (TRUE);
+
+ escape:
+  if (vec) free (vec);
+  if (values) free (values);
+  if (invec) FreeVectorArray (invec, Ninvec);
+  dbFreeFields (fields, Nfields);
+  SkyListFree (skylist);
+  return (FALSE);
+
+ help:
+  gprint (GP_ERR, "USAGE: mmatch (RA) (DEC) (RADIUS) field[,field,field...] [-index index]\n");
+  gprint (GP_ERR, "   OR: mmatch -coords (filename.fits) (RADIUS) field[,field,field...] [-index index]\n");
+
+  if ((argc > N + 1) && !strcasecmp (argv[N+1], "fields")) {
+    gprint (GP_ERR, "  RA : right ascension (J2000) [degrees]\n");
+    gprint (GP_ERR, "  DEC : declination [degrees]\n");
+    return (FALSE);
+  }
+  gprint (GP_ERR, " mmatch --help fields : for a complete listing of allowed fields\n");
+  return (FALSE);
+}