Index: /trunk/Ohana/src/dvolens/Makefile
===================================================================
--- /trunk/Ohana/src/dvolens/Makefile	(revision 39486)
+++ /trunk/Ohana/src/dvolens/Makefile	(revision 39487)
@@ -23,4 +23,7 @@
 
 DVOLENS = \
+$(SRC)/load_images.$(ARCH).o 	 \
+$(SRC)/select_images.$(ARCH).o 	 \
+$(SRC)/WarpImageMaps.$(ARCH).o 	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
 $(SRC)/SetSignals.$(ARCH).o 	 \
@@ -39,4 +42,7 @@
 
 DVOLENS_CLIENT = \
+$(SRC)/load_images.$(ARCH).o 	 \
+$(SRC)/select_images.$(ARCH).o 	 \
+$(SRC)/WarpImageMaps.$(ARCH).o 	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
 $(SRC)/SetSignals.$(ARCH).o 	 \
Index: /trunk/Ohana/src/dvolens/include/dvolens.h
===================================================================
--- /trunk/Ohana/src/dvolens/include/dvolens.h	(revision 39486)
+++ /trunk/Ohana/src/dvolens/include/dvolens.h	(revision 39487)
@@ -27,4 +27,5 @@
 char   CATMODE[16];    /* raw, mef, split, mysql */
 char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
+char   ImageCat[DVO_MAX_PATH];
 
 int    HOST_ID;
@@ -39,4 +40,6 @@
 int    UPDATE;
 int    NTHREADS;
+
+int    REPAIR_LENSING_IDS;
 
 DvoLensMode MODE;
@@ -82,2 +85,16 @@
 int myIndexSetEntry (myIndexType *myIndex, int value, int entry);
 int myIndexGetEntry (myIndexType *myIndex, int value);
+
+int FindWarpGroups (void);
+int RecoverLensingIndex (Average *average, myIndexType *measureIndex, Lensing *lensing);
+void FreeWarpGroups (void);
+
+int isGPC1warp (int photcode);
+
+Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t *Nimage);
+
+int load_images (SkyList *skylist);
+void free_images (void);
+Image *getimages (off_t *N);
+off_t getImageByID (off_t ID);
+
Index: /trunk/Ohana/src/dvolens/src/ConfigInit.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/ConfigInit.c	(revision 39486)
+++ /trunk/Ohana/src/dvolens/src/ConfigInit.c	(revision 39487)
@@ -21,4 +21,6 @@
   CATDIR = abspath (tmpcatdir, DVO_MAX_PATH);
   free (tmpcatdir);
+
+  snprintf (ImageCat, DVO_MAX_PATH, "%s/Images.dat", CATDIR);
 
   ScanConfig (config, "CATMODE",                "%s",  0, CATMODE);
Index: /trunk/Ohana/src/dvolens/src/WarpImageMaps.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/WarpImageMaps.c	(revision 39487)
+++ /trunk/Ohana/src/dvolens/src/WarpImageMaps.c	(revision 39487)
@@ -0,0 +1,170 @@
+# include "dvolens.h"
+
+// we have a list of warp images, each with obstime & photcode.  I need to make a list of the obstime groups
+
+# define D_NOBSTIMES 10000
+
+typedef struct {
+  int obstime;
+  int Nwarps;
+  int NWARPS;
+  off_t *warps;
+  unsigned short photcode;
+} WarpGroup;
+
+WarpGroup *warpgroup = NULL;
+int Nwarpgroup = 0;
+int NWARPGROUP = 0;
+
+myIndexType *warpObstimeIndex = NULL;
+
+static off_t Nimage = 0;
+static Image *image = NULL;
+
+int FindWarpGroups (void) {
+
+  // keep a local static copy of the image pointer
+  image = getimages (&Nimage);
+
+  int *obstimes = NULL;
+  int Nobstimes = 0;
+  int NOBSTIMES = D_NOBSTIMES;
+  ALLOCATE (obstimes, int, NOBSTIMES);
+
+  // first generate a list of obstime values
+  for (off_t i = 0; i < Nimage; i++) {
+    if (!isGPC1warp(image[i].photcode)) continue;
+
+    obstimes[Nobstimes] = image[i].tzero;
+    Nobstimes ++;
+
+    CHECK_REALLOCATE (obstimes, int, NOBSTIMES, Nobstimes, D_NOBSTIMES);
+  }
+
+  // sort the list:
+  isort (obstimes, Nobstimes);
+
+  // find the unique obstimes
+  int *uniqtimes = NULL;
+  int *uniqcount = NULL;
+  int Nuniqtimes = 0;
+  ALLOCATE (uniqtimes, int, Nobstimes);
+  ALLOCATE (uniqcount, int, Nobstimes);
+
+  // generate a uniq set of obstimes
+  for (int i = 0; i < Nobstimes; Nuniqtimes ++) {
+    uniqtimes[Nuniqtimes] = obstimes[i];
+    int Ndup = 0;
+    int lastValue = uniqtimes[Nuniqtimes];
+    while ((i < Nobstimes) && (obstimes[i] == lastValue)) {
+      i++;
+      Ndup ++;
+      uniqcount[Nuniqtimes] = Ndup;
+    }
+  }
+  REALLOCATE (uniqtimes, int, Nuniqtimes);
+  REALLOCATE (uniqcount, int, Nuniqtimes);
+  
+  // now I need to assign all warp images to a warpgroup
+  // I have 2 options to go from warp obstime to warpgroup:
+  // a) use a bisection lookup [may be slow]
+  // b) use an index on obstime [require ~5*3.14*1e7*4 bytes ~ 600MB for index]
+
+  // generate an index on obstime
+  warpObstimeIndex = myIndexAlloc ();
+  myIndexInit (warpObstimeIndex);
+  warpObstimeIndex->minID = obstimes[0];
+  warpObstimeIndex->maxID = obstimes[Nobstimes-1];
+
+  myIndexSetRange (warpObstimeIndex);
+
+  // generate the warp groups (create the index as we go)
+  Nwarpgroup = Nuniqtimes;
+  ALLOCATE (warpgroup, WarpGroup, Nwarpgroup);
+  
+  for (int i = 0; i < Nuniqtimes; i++) {
+    warpgroup[i].obstime = uniqtimes[i];
+    warpgroup[i].photcode = 0; // not yet set
+    warpgroup[i].Nwarps = 0;
+    warpgroup[i].NWARPS = 100;
+    ALLOCATE (warpgroup[i].warps, off_t, warpgroup[i].NWARPS);
+    myIndexSetEntry (warpObstimeIndex, uniqtimes[i], i);
+  }
+
+  // assign all warps to one of the warpgroups
+  for (off_t i = 0; i < Nimage; i++) {
+    if (!isGPC1warp(image[i].photcode)) continue;
+
+    int seq = myIndexGetEntry (warpObstimeIndex, image[i].tzero);
+    myAssert (warpgroup[seq].obstime == image[i].tzero, "oops");
+
+    if (!warpgroup[seq].photcode) {
+      warpgroup[seq].photcode = image[i].photcode;
+    } else {
+      myAssert (warpgroup[seq].photcode == image[i].photcode, "oops");
+    }
+
+    int N = warpgroup[seq].Nwarps;
+    warpgroup[seq].warps[N] = i;
+    warpgroup[seq].Nwarps ++;
+    
+    CHECK_REALLOCATE (warpgroup[seq].warps, off_t, warpgroup[seq].NWARPS, warpgroup[seq].Nwarps, 100);
+  }
+
+  free (uniqtimes); 
+  free (uniqcount); 
+  free (obstimes);
+  return TRUE;
+}
+
+int RecoverLensingIndex (Average *average, myIndexType *measureIndex, Lensing *lensing) {
+
+  // Although the imageID is wrong (no matching measures), the image should be one of the
+  // correct set of warps with the same obstime. Find that set so we can find a matching image
+  int im = getImageByID (lensing->imageID);
+  if (im < 0) {
+    return -1;
+  }
+
+  int seq = myIndexGetEntry (warpObstimeIndex, image[im].tzero);
+  myAssert (seq > -1, "oops");
+
+  // can we find an image in this set which matches one of our measures?
+
+  // we now have the warp group, but which is the correct warp?
+  for (int i = 0; i < warpgroup[seq].Nwarps; i++) {
+
+    // XXX check on Rmin,Rmax,Dmin,Dmax
+    // average->R,D could pin-point more quickly
+
+    // this is a possible image: check the coordinates:
+    int N = warpgroup[seq].warps[i];
+
+    int Mj = myIndexGetEntry (measureIndex, image[N].imageID);
+    if (Mj < 0) continue;
+
+    // We have a match between a warp measure and an image in this obstime group.  this
+    // may not be the only possible solution, but any other solution will have to be one
+    // of the matching overlapping images.  I'm not sure I can distinguish these.
+    lensing->oldImID = lensing->imageID;
+    lensing->imageID = image[N].imageID;
+    return Mj;
+  }
+  
+  return -1;
+}
+
+void FreeWarpGroups (void) {
+
+  for (int i = 0; i < Nwarpgroup; i++) {
+    free (warpgroup[i].warps);
+    // free (warpgroup[i].Rmin);
+    // free (warpgroup[i].Rmax);
+    // free (warpgroup[i].Dmin);
+    // free (warpgroup[i].Dmax);
+    // free (warpgroup[i].onBoundary);
+  }    
+
+  free (warpgroup);
+  myIndexFree (warpObstimeIndex);
+}
Index: /trunk/Ohana/src/dvolens/src/args.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/args.c	(revision 39486)
+++ /trunk/Ohana/src/dvolens/src/args.c	(revision 39487)
@@ -37,4 +37,10 @@
   if ((N = get_argument (argc, argv, "-vv"))) {
     VERBOSE2 = VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  REPAIR_LENSING_IDS = FALSE;
+  if ((N = get_argument (argc, argv, "-repair-lensing-ids"))) {
+    REPAIR_LENSING_IDS = TRUE;
     remove_argument (N, &argc, argv);
   }
@@ -162,4 +168,10 @@
   }
 
+  REPAIR_LENSING_IDS = FALSE;
+  if ((N = get_argument (argc, argv, "-repair-lensing-ids"))) {
+    REPAIR_LENSING_IDS = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
   UPDATE = FALSE;
   if ((N = get_argument (argc, argv, "-update"))) {
Index: /trunk/Ohana/src/dvolens/src/dvolens.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/dvolens.c	(revision 39486)
+++ /trunk/Ohana/src/dvolens/src/dvolens.c	(revision 39487)
@@ -13,4 +13,6 @@
       FREE (UserCatalog);
       FREE (CATDIR);
+      free_images();
+      FreeWarpGroups();
       ohana_memcheck (VERBOSE);
       ohana_memdump (VERBOSE);
Index: /trunk/Ohana/src/dvolens/src/dvolens_client.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/dvolens_client.c	(revision 39486)
+++ /trunk/Ohana/src/dvolens/src/dvolens_client.c	(revision 39487)
@@ -19,4 +19,6 @@
       FREE (HOSTDIR);
       FREE (CATDIR);
+      free_images();
+      FreeWarpGroups();
       ohana_memcheck (VERBOSE);
       ohana_memdump (VERBOSE);
Index: /trunk/Ohana/src/dvolens/src/load_images.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/load_images.c	(revision 39487)
+++ /trunk/Ohana/src/dvolens/src/load_images.c	(revision 39487)
@@ -0,0 +1,127 @@
+# include "dvolens.h"
+
+static off_t Nsubset = 0;
+static Image *subset = NULL;
+
+// MAX_ID requires 512M to store the image index
+# define MAX_ID      0x8000000
+static off_t         minImageID = MAX_ID;
+static off_t         maxImageID = 0;
+
+// as an alternative, we generate imageSeq, which directly maps imageID -> seq
+static off_t *imageIDs = NULL; // list of all image IDs
+static off_t *imageIdx = NULL; // list of index for image IDs 
+static off_t *imageSeq = NULL; // list of index for image IDs 
+
+void initImages (Image *image, off_t Nimage);
+
+int load_images (SkyList *skylist) {
+
+  FITS_DB db;
+  
+  gfits_db_init (&db);
+  
+  /* lock and load the image db table */
+  int status = dvo_image_lock (&db, ImageCat, 60.0, LCK_SOFT);
+  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
+  if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
+  if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
+  // the raw FITS data is freed by dvo_image_load, leaving only the Image structure data
+  
+  INITTIME;
+
+  // convert database table to internal structure
+  off_t Nimage;
+  Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.scaledValue, &db.nativeOrder);
+  if (!image) {
+      fprintf (stderr, "ERROR: failed to read images\n");
+      exit (2);
+  }
+  MARKTIME("  convert image table: %f sec\n", dtime);
+
+  subset = select_images (skylist, image, Nimage, &Nsubset);
+  MARKTIME("  select images: %f sec\n", dtime);
+
+  // if we have generated an image subset and we are running UPDATE_OFFSETS, the we can free images here
+  dvo_image_unlock (&db); 
+  gfits_db_free (&db);
+
+  initImages (subset, Nsubset);
+
+  return TRUE;
+}
+
+void initImages (Image *image, off_t Nimage) {
+
+  ALLOCATE (imageIDs, off_t, Nimage);
+  ALLOCATE (imageIdx, off_t, Nimage);
+
+  off_t i;
+  for (i = 0; i < Nimage; i++) {
+    imageIdx[i] = i;
+    imageIDs[i] = image[i].imageID;
+    minImageID = MIN(minImageID, image[i].imageID);
+    maxImageID = MAX(maxImageID, image[i].imageID);
+    myAssert (image[i].imageID < MAX_ID, "image IDs too large for index memory");
+  }
+  llsortpair (imageIDs, imageIdx, Nimage);
+
+  ALLOCATE (imageSeq, off_t, maxImageID + 1);
+  for (i = 0; i < maxImageID + 1; i++)  {
+    imageSeq[i] = -1; // not yet assigned
+  }
+
+  for (i = 0; i < Nimage; i++) {
+    off_t N = image[i].imageID;
+    myAssert (imageSeq[N] == -1, "previously assigned");
+    imageSeq[N] = i;
+  }
+}
+
+Image *getimages (off_t *N) {
+
+  *N = Nsubset;
+  return subset;
+}
+
+void free_images (void) {
+
+  if (!subset) return;
+
+  FREE (subset);
+  FREE (imageIDs);
+  FREE (imageIdx);
+  FREE (imageSeq);
+}
+
+off_t getImageByID (off_t ID) {
+  if (imageSeq) {
+    if (ID < minImageID) return (-1);
+    if (ID > maxImageID) return (-1);
+    off_t N = imageSeq[ID];
+    return N;
+  }
+
+  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
+  // use bisection to find the specified image ID
+
+  off_t Nlo, Nhi, N;
+
+  Nlo = 0; Nhi = Nsubset;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (imageIDs[N] < ID) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nsubset);
+    }
+  }
+
+  for (N = Nlo; N < Nhi; N++) {
+    if (imageIDs[N] == ID)
+      return (imageIdx[N]);
+  }
+
+  return (-1);
+}
+
Index: /trunk/Ohana/src/dvolens/src/select_images.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/select_images.c	(revision 39487)
+++ /trunk/Ohana/src/dvolens/src/select_images.c	(revision 39487)
@@ -0,0 +1,110 @@
+# include "dvolens.h"
+
+// for the purposes of dvolens / RepairLensing, we are only keeping warp images
+int isGPC1warp (int photcode);
+
+Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t *Nimage) {
+  
+  Image *image;
+  off_t i, j;
+  double Ri[5], Di[5], Xi[5], Yi[5];
+  
+  double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion;
+
+  if (skylist[0].Nregions < 1) {
+    *Nimage = 0;
+    if (VERBOSE) fprintf (stderr, "no matching sky regions\n");
+    return NULL;
+  }
+
+  INITTIME;
+
+  RminSkyRegion = +360.0;
+  RmaxSkyRegion = -360.0;
+  DminSkyRegion = +90.0;
+  DmaxSkyRegion = -90.0;
+
+  /* compare with each region file */
+  for (i = 0; i < skylist[0].Nregions; i++) { 
+    RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
+    RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
+    DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
+    DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
+  }
+  RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
+  MARKTIME("create sky region coords: %f sec\n", dtime);
+
+  if (VERBOSE) fprintf (stderr, "finding images\n");
+
+  off_t D_NIMAGE = 6000;
+
+  off_t nimage = 0;
+  off_t NIMAGE = D_NIMAGE;
+  ALLOCATE (image, Image, NIMAGE);
+  
+  // go through the complete list of images, selecting ones which overlap any region
+  for (i = 0; i < Ntimage; i++) {
+      
+    if (!isGPC1warp(timage[i].photcode)) continue;
+   
+    /* define image corners - only Warps (TAN) /*/
+    Xi[0] = 0;            Yi[0] = 0;
+    Xi[1] = timage[i].NX; Yi[1] = 0;
+    Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
+    Xi[3] = 0;            Yi[3] = timage[i].NY;
+    Xi[4] = 0;            Yi[4] = 0;
+
+    /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */
+    double RminImage = RmidSkyRegion + 180.0;
+    double RmaxImage = RmidSkyRegion - 180.0;
+    double DminImage = +90.0;
+    double DmaxImage = -90.0;
+    for (j = 0; j < 5; j++) {
+      XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
+      Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion);
+
+      RminImage = MIN(RminImage, Ri[j]);
+      RmaxImage = MAX(RmaxImage, Ri[j]);
+      DminImage = MIN(DminImage, Di[j]);
+      DmaxImage = MAX(DmaxImage, Di[j]);
+    }
+
+    // check that this image is even in range of the searched region
+    if (DminImage > DmaxSkyRegion) continue;
+    if (DmaxImage < DminSkyRegion) continue;
+    
+    // the sky region RA is defined to be 0 - 360.0
+    if (RminImage > RmaxSkyRegion) continue;
+    if (RmaxImage < RminSkyRegion) continue;
+
+    image[nimage] = timage[i]; 
+
+    nimage ++;
+    if (nimage == NIMAGE) {
+      NIMAGE += D_NIMAGE;
+      D_NIMAGE = MAX (100000, D_NIMAGE * 1.5);
+      REALLOCATE (image, Image, NIMAGE);
+    }
+  }
+  MARKTIME("finish image selection: %f sec\n", dtime);
+
+  if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);
+
+  REALLOCATE (image, Image, MAX (nimage, 1));
+
+  *Nimage  = nimage;
+  return (image);
+}
+
+// for now (20140710) I need to identify gpc1 stacks explicitly.  generalize in the future
+int isGPC1warp (int photcode) {
+
+  if (photcode == 12000) return TRUE; // g-band
+  if (photcode == 12100) return TRUE; // r-band
+  if (photcode == 12200) return TRUE; // i-band
+  if (photcode == 12300) return TRUE; // z-band
+  if (photcode == 12400) return TRUE; // y-band
+  if (photcode == 12500) return TRUE; // w-band
+
+  return FALSE;
+}
Index: /trunk/Ohana/src/dvolens/src/update_objects.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/update_objects.c	(revision 39486)
+++ /trunk/Ohana/src/dvolens/src/update_objects.c	(revision 39487)
@@ -37,4 +37,11 @@
 
   if (VERBOSE) fprintf (stderr, "re-loading catalog data\n");
+
+  // XXX I need to load enough of the images to load the warp groups
+
+  if (REPAIR_LENSING_IDS) {
+    load_images (skylist);
+    FindWarpGroups();
+  }
 
   /* load data from each region file */
@@ -150,20 +157,23 @@
     table->hosts[i].pathname = tmppath;
 
-    char command[1024];
-    snprintf (command, 1024, "dvolens_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", 
-	      table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, 
-	      UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
-
-    // options & configs which affect dvolens_client -update-catalogs
+    char *command = NULL;
+    strextend (&command, "dvolens_client -update-objects");
+    strextend (&command, "-hostID %d", table->hosts[i].hostID);
+    strextend (&command, "-D CATDIR %s", CATDIR);
+    strextend (&command, "-hostdir %s", table->hosts[i].pathname);
+    strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+
+    // options & configs which affect dvolens_client -update-catalogs:
     // VERBOSE, VERBOSE2, UPDATE
-    
-    char tmpline[1024];
-    if (VERBOSE)          { snprintf (tmpline, 1024, "%s -v",           	command);                    			  strcpy (command, tmpline); }
-    if (VERBOSE2)         { snprintf (tmpline, 1024, "%s -vv",          	command); 		     			  strcpy (command, tmpline); }
-    if (UPDATE)           { snprintf (tmpline, 1024, "%s -update",       	command); 		     			  strcpy (command, tmpline); }
+
+    if (VERBOSE)            strextend (&command, "-v");
+    if (VERBOSE2)           strextend (&command, "-vv");
+    if (UPDATE)             strextend (&command, "-update");
+    if (REPAIR_LENSING_IDS) strextend (&command, "-repair-lensing-ids");
 
     fprintf (stderr, "command: %s\n", command);
 
     if (PARALLEL_MANUAL) { 
+      free (command);
       continue;
     }
@@ -185,4 +195,5 @@
       table->hosts[i].pid = pid; // save for future reference
     }
+    free (command);
   }
 
@@ -198,4 +209,6 @@
     }
   }
+
+  FreeHostTable (table);
   return (TRUE);
 }      
Index: /trunk/Ohana/src/dvolens/src/update_objects_catalog.c
===================================================================
--- /trunk/Ohana/src/dvolens/src/update_objects_catalog.c	(revision 39486)
+++ /trunk/Ohana/src/dvolens/src/update_objects_catalog.c	(revision 39487)
@@ -83,4 +83,13 @@
 
       Mj = myIndexGetEntry (measureIndex, lensing->imageID);
+
+      // if we are unable to find a match in the measure table, we may have the wrong 
+      // imageID.  there was a bug in which we only matched to the correct warp obstime,
+      // but got the wrong skycell.  if this is the case, we can try to recover by finding
+      // the set of warps which match our obstime, then choosing the one from that set
+      // which matches one of our warps.
+      if ((Mj < 0) && REPAIR_LENSING_IDS) {
+	Mj = RecoverLensingIndex (average, measureIndex, lensing);
+      }
       myAssert (Mj > -1, "missing index");