Index: /branches/eam_branches/ipp-20130307/Ohana/src/relastro/src/hpm_objects.c
===================================================================
--- /branches/eam_branches/ipp-20130307/Ohana/src/relastro/src/hpm_objects.c	(revision 35286)
+++ /branches/eam_branches/ipp-20130307/Ohana/src/relastro/src/hpm_objects.c	(revision 35286)
@@ -0,0 +1,317 @@
+# include "relastro.h"
+
+# define NEXT_I { if (Ngroup < 2) slowMoving[ni] = TRUE; newI = TRUE; i++; continue; }
+# define NEXT_J { j++; continue; }
+
+int hpm_objects (SkyRegion *region, Catalog *catalog) {
+
+  off_t i, j, m, J, ni, nj, *N1;
+  off_t Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad;
+  off_t l, i1, NAVERAGE, NMEASURE;
+  int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas, Nepoch, nv[2], Nmatchmeasobj;
+  int foundA, foundB, XVERB;
+  double *X1, *Y1;
+  double dX, dY, dR, RADIUS2;
+  Coords tcoords;
+  Catalog catalogOut;
+  Catalog testcat;
+
+  int Nsecfilt;
+  char filename[1024];
+
+  // XXX are we saving these in an hpm dvodb?
+  snprintf (filename, 1024, "%s/%s.cpt", HIGH_SPEED_DIR, region[0].name);
+  fprintf (stderr, "%s\n",filename);
+
+  dvo_catalog_init (&catalogOut, TRUE); /*initialise new catalogue*/
+  catalogOut.filename = strcreate(filename);
+
+  SIGMA_LIM = 0.0;
+
+  Nsecfilt = GetPhotcodeNsecfilt();
+  // off_t Naverage = catalog[0].Naverage;
+  // off_t Nmeasure = catalog[0].Nmeasure;
+  catalogOut.filename = filename; // based on the input name, need to keep everything below the catdir portion
+  catalogOut.Nsecfilt = Nsecfilt;
+  catalogOut.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; // load all data
+  
+  catalogOut.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
+  catalogOut.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
+
+  if (!dvo_catalog_open (&catalogOut, region, VERBOSE2,"w")) {
+    fprintf (stderr, "ERROR: failure to open catalog file %s\n",
+	     filename);
+    exit (2);
+  }
+
+  NAVERAGE = 1000;
+  NMEASURE = 10000;
+  REALLOCATE (catalogOut.average, Average, NAVERAGE);
+  REALLOCATE (catalogOut.measure, Measure, NMEASURE);
+  REALLOCATE (catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
+
+  // testcat is used to determine the fit for a single object group
+  // objects which do not have a high-quality testcat fit are not kept
+  dvo_catalog_init (&testcat, TRUE);
+  testcat.Naverage = 1; // this is fixed -- only one obj in testcat
+  NMEASURE_TEST = 1000;
+  REALLOCATE (catalogOut.average, Average, 1);
+  REALLOCATE (catalogOut.measure, Measure, NMEASURE_TEST);
+  REALLOCATE (catalogOut.secfilt, SecFilt, Nsecfilt);
+
+  // we need at least 2 objects if we are going to match anything...
+  if (catalog[0].Naverage < 2) return (TRUE);
+
+  // mask with which to mark objects to be ignored
+  ALLOCATE (slowMoving, int, catalog[0].Naverage);
+  memset (slowMoving, 0, catalog[0].Naverage*sizeof(int));
+
+  if (VERBOSE) fprintf (stderr, "checking "OFF_T_FMT" objects\n",  catalog[0].Naverage);
+  Nslow = 0;
+  Ninvalid = 0;
+  NgroupAbad = 0;
+  NgroupBbad = 0;
+
+  // mark (exclude) objects with both sets of target photcodes
+  for (i = 0; i < catalog[0].Naverage; i++) {
+
+    XVERB  = (catalog[0].average[i].objID == OBJ_ID_SRC) && (catalog[0].average[i].catID == CAT_ID_SRC);
+    XVERB |= (catalog[0].average[i].objID == OBJ_ID_DST) && (catalog[0].average[i].catID == CAT_ID_DST);
+    if (XVERB) {
+      fprintf (stderr, "test object\n");
+    }
+
+    // selection criteria:
+    // (Nps1 > 2) && (Trange < 180)
+    if (catalog[0].average[i].Trange > MAX_TRANGE) {
+      slowMoving[i] = TRUE;
+      Nslow ++;
+      continue;
+    }
+    // count the PS1 detections via explicit photcode ranges?
+    // XXX this is a total hard-wired hack...
+    int Nps1 = 0;
+    for (j = 0; i < 5; j++) {
+      Nps1 += catalog[0].secfilt[Nsecfilt*i+j].Ncode;
+    }
+    if (Nps1 < MIN_PS1_DET) {
+      slowMoving[i] = TRUE;
+      Nslow ++;
+      continue;
+    }      
+  }
+
+  fprintf (stderr, OFF_T_FMT" slow, "OFF_T_FMT" total objects; "OFF_T_FMT" possible fast\n",  Nslow,  catalog[0].Naverage,  catalog[0].Naverage - Nslow);
+    
+  // double loop over unmarked objects (sorted in RA / X)
+  // record all pairs within the desired match distance
+
+  // create tmp local positional index
+  ALLOCATE (X1, double, catalog[0].Naverage);
+  ALLOCATE (Y1, double, catalog[0].Naverage);
+  ALLOCATE (N1, off_t,  catalog[0].Naverage);
+
+  // define a local projection
+  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
+  if (region[0].Dmax < 90) {
+    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
+  } else {
+    tcoords.crval2 = 90.0;
+  }
+  tcoords.crpix1 = 0;
+  tcoords.crpix2 = 0;
+  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+  tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
+  tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+  tcoords.Npolyterms = 1;
+  strcpy (tcoords.ctype, "RA---ARC");
+
+  /* build spatial index (RA sort) referencing input array sequence */
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    status = RD_to_XY (&X1[i], &Y1[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
+    N1[i] = i;
+    assert (status);
+  }
+  sort_coords_index (X1, Y1, N1, catalog[0].Naverage);
+
+  Nmatch = 0;
+  Nmatchmeas = 0;
+  RADIUS2 = SQ(RADIUS);
+
+  // group is a list of objects that are within a clump.  this set will be tested
+  // via a clipped fit to the measurements after the group is identified
+  Ngroup = 0;
+  NGROUP = 100;
+  ALLOCATE (group, int, NGROUP);
+
+  // in the loop below, we need to do a bunch of things when we go to the next main object
+  newI = TRUE;
+
+  // mark (skip) objects with both sets of target photcodes
+  // the loop below is attempting to find associations of multiple objects which have 
+  // passed the cuts above.  the index i is following the primary object of interest
+  // the index j is used to explore possible near neighbors.  
+  // When we go to the next object 'i', Nmatch is reset 
+  for (i = j = 0; (i < catalog[0].Naverage) && (j < catalog[0].Naverage);) {
+
+    ni = N1[i];
+    nj = N1[j];
+
+    if (newI) {
+      Nmatch = 0; 
+      Ngroup = 1;
+      group[0] = ni;
+      newI = FALSE;
+    }
+
+    XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
+    XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST);
+    if (XVERB) {
+      fprintf (stderr, "test object %d or %d\n", (int) ni, (int) nj);
+    }
+
+    if (slowMoving[ni]) NEXT_I;
+    if (slowMoving[nj]) NEXT_J;
+
+    if (!finite(X1[i]) || !finite(Y1[i])) NEXT_I;
+    if (!finite(X1[j]) || !finite(Y1[j])) NEXT_J;
+
+    // look for pairs that are within the maximum separation
+    dX = X1[i] - X1[j];
+    if (XVERB) {
+      fprintf (stderr, "%d %d : %f\n", (int) ni, (int) nj, dX);
+    }
+
+    if (dX <= -1.02*RADIUS) NEXT_I; // negative dX: i is too small
+    if (dX >= +1.02*RADIUS) NEXT_J; // positive dX: j is too small
+
+    // within match range; look for valid matches
+    for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {     
+      if (J == i) continue;  // avoid auto-matches
+      nj = N1[J];
+
+      dX = X1[i] - X1[J];
+
+      if (slowMoving[nj]) continue;
+
+      XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
+      XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST);
+      if (XVERB) {
+	fprintf (stderr, "test object pt2 %d or %d : %f\n", (int) ni, (int) nj, dX);
+	fprintf (stderr, ".\n");
+      }
+
+      dY = Y1[i] - Y1[J];
+      dR = dX*dX + dY*dY;
+      if (dR > RADIUS2) continue;
+
+      /*** a match is found ***/
+      group[Ngroup] = nj;
+      Ngroup ++;
+      CHECK_REALLOCATE (group, int, NGROUP, Ngroup, 100);
+    }
+
+    if (Ngroup < 2) NEXT_I;
+
+    // we now have spatially associated group of objects.  now we need to see if the set of
+    // measurements can be fitted reasonably with a proper motion (& parallax?)
+
+    // the mean object will start with info from the primary object
+    // remember: testcat.Naverage = 1 -- does not change
+    testcat.average[0] = catalog[0].average[group[0]];
+    testcat.average[0].measureOffset = 0;
+    Nmatchmeas = 0;
+    I = group[0];
+    for (J = 0; J < Ngroup; J++) {
+      J = group[J];
+      m = catalog[0].average[Nj].measureOffset;
+      for (k = 0; k < catalog[0].average[Nj].Nmeasure; k++) {
+	testcat.measure[Nmatchmeas] = catalog[0].measure[m+k];
+	/* Set offset RA and Dec wrt correct average value*/
+	testcat.measure[Nmatchmeas].dR = catalog[0].measure[m+k].dR + 3600.0*(catalog[0].average[I].R - catalog[0].average[J].R);
+	testcat.measure[Nmatchmeas].dD = catalog[0].measure[m+k].dD + 3600.0*(catalog[0].average[I].D - catalog[0].average[J].D);
+	testcat.measure[Nmatchmeas].averef = 0;
+	Nmatchmeas++;
+	CHECK_REALLOCATE (testcat.measure, Measure, NMEASURE_TEST, Nmatchmeas, 1000);
+      }
+    }
+    testcat.average[0].Nmeasure = Nmatchmeas;
+
+    // we have now accumulated the measurements for this group, let's try a fit
+    // this needs to be a (fairly robust) clipped fit or we will have a hard time
+    // distinguishing a bad fit from a fit with 1 or 2 bad points
+    FIT_MODE = FIT_PM_ONLY;
+    UpdateObjects (&testcat, 1);
+
+    // logic for keeping the fit:
+    if (testcat.average[0].ChiSqPM < XXX) good = TRUE;
+    if (testcat.average[0].Npos > 0.X * testcat.average[0].Nmeasure) good = TRUE;
+    // what else?
+
+    if (good) {
+      // save the new object on catalogOut
+      {
+	catalogOut.average[Nout] = testcat.average[0];
+	catalogOut.average[Nout].measureOffset = Nmatchmeas;
+	m = testcat.average[0].measureOffset;
+	for (k = 0; k < testcat.average[0].Nmeasure; k++) {
+	  catalogOut.measure[Noutmeas] = testcat.measure[m+k];
+	  catalogOut.measure[Noutmeas].averef = Nout;
+	}
+	Nout ++;
+	CHECK_REALLOCATE (catalogOut.average, Average, NAVERAGE, Nout, 100);
+      }
+    }
+    NEXT_I;
+  }
+
+  catalogOut.Naverage=Nmatch;
+  catalogOut.Nmeasure=Nmatchmeas;
+  catalogOut.Nsecfilt=Nsecfilt;
+  catalogOut.Nsecf_mem=Nmatch*Nsecfilt;
+
+  populate_tiny_values (&catalogOut, DVO_TV_MEASURE);
+  UpdateObjects (&catalogOut, 1);
+
+  fprintf (stderr, "found %d matches\n", Nmatch);
+  dvo_catalog_save (&catalogOut, VERBOSE2);
+  dvo_catalog_unlock (&catalogOut);
+  dvo_catalog_free (&catalogOut);
+  free (slowMoving);
+  free (groupA);
+  free (groupB);
+  free (X1);
+  free (Y1);
+  free (N1);
+
+  return TRUE;
+}
+
+// return TRUE if measure->photcode match any in the photcodeSet
+// (but, return FALSE if the measurement is bad)
+int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset) {
+
+  int found, k;
+
+  if (!Nset) return TRUE;
+  // XXX require a set or not?  assert (Nset > 0);
+
+  if (!finite(measure[0].dR)) return FALSE;
+  if (!finite(measure[0].dD)) return FALSE;
+  if (!finite(measure[0].M))  return FALSE;
+  
+  float dX = GetAstromError (measure, ERROR_MODE_RA);
+  if (isnan(dX)) return FALSE;
+
+  float dY = GetAstromError (measure, ERROR_MODE_DEC);
+  if (isnan(dY)) return FALSE;
+
+  /* select measurements by photcode, or equiv photcode, if specified */
+  found = FALSE;
+  for (k = 0; (k < Nset) && !found; k++) {
+    if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE;
+    if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+  }
+
+  return found;
+}