Index: /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/doc/notes.txt
===================================================================
--- /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/doc/notes.txt	(revision 27424)
+++ /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/doc/notes.txt	(revision 27425)
@@ -1,2 +1,44 @@
+
+2010.03.23
+
+  I've been working on optimization.  I've updated ImageOps to make
+  the measure->image relationship generation use the image IDs along
+  with a bracket search -- this goes much faster as a result.  More
+  work is needed:
+
+  * select_images needs to apply astrometry to all image corners: 
+    -> cache the image center and radius and use this to narrow down
+       the searches.
+
+  * MosaicOps uses the time to match mosaics.  I've added sorting to
+    speed this up, but this should be done based on the image and
+    parent index.
+
+    -> make a tool to create the parent ID for existing DBs.
+    -> make sure addstar is populating the parent IDs
+    -> use the parent IDs to make the link.
+
+load image data: 0.006347 sec
+  setup sky: 0.091475 sec
+  convert image table: 0.091498 sec
+  select images: 2.925978 sec
+  init images: 2.926957 sec
+  init mosaics: 2.927319 sec
+load images: 2.933762 sec
+load catalog data: 9.772259 sec
+make image bins: 9.780618 sec
+set up image indexes: 10.196519 sec
+
+load image data: 0.006746 sec
+  setup sky: 0.091440 sec
+  convert image table: 0.091467 sec
+  select images: 0.374193 sec
+  init images: 0.375355 sec
+  init mosaics: 0.375753 sec
+load images: 0.382606 sec
+load catalog data: 7.260094 sec
+make image bins: 7.268564 sec
+set up image indexes: 7.687485 sec
+
 
 2008.03.01
Index: /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/ImageOps.c	(revision 27424)
+++ /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/ImageOps.c	(revision 27425)
@@ -1,7 +1,5 @@
 # include "relastro.h"
-#define TESTING
-
-static unsigned int *start;
-static unsigned int *stop;
+# define TESTING
+# define USE_IMAGE_ID 1
 
 static off_t       **bin;     // link from catalog,measure to image
@@ -14,4 +12,13 @@
 static off_t        Nimage;   // number of available images
 
+# if USE_IMAGE_ID
+// we sort (imageIDs, imageIdx) by imageIDs to get a sorted index
+static off_t        *imageIDs; // list of all image IDs
+static off_t        *imageIdx; // list of index for image IDs 
+# else
+static unsigned int *start;
+static unsigned int *stop;
+# endif
+
 Image *getimages (off_t *N) {
   *N = Nimage;
@@ -30,6 +37,16 @@
   Nimage = N;
 
-  ALLOCATE (start,   unsigned, Nimage);
-  ALLOCATE (stop,    unsigned, Nimage);
+# if USE_IMAGE_ID
+  ALLOCATE (imageIDs, off_t, Nimage);
+  ALLOCATE (imageIdx, off_t, Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+    imageIdx[i] = i;
+    imageIDs[i] = image[i].imageID;
+  }
+  llsortpair (imageIDs, imageIdx, Nimage);
+# else
+  ALLOCATE (start, unsigned, Nimage);
+  ALLOCATE (stop, unsigned, Nimage);
 
   for (i = 0; i < Nimage; i++) {
@@ -37,4 +54,33 @@
     stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
   }
+# endif
+
+}
+
+off_t getImageByID (off_t ID) {
+
+  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
+  // use bisection to find the specified image ID
+
+# if USE_IMAGE_ID
+  off_t Nlo, Nhi, N;
+
+  Nlo = 0; Nhi = Nimage;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (imageIDs[N] < ID) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nimage);
+    }
+  }
+
+  for (N = Nlo; N < Nhi; N++) {
+    if (imageIDs[N] == ID)
+      return (imageIdx[N]);
+  }
+# endif
+
+  return (-1);
 }
 
@@ -97,7 +143,42 @@
 }
 
+# if USE_IMAGE_ID
 /* modify this function to use the measure->imageID field */
 void matchImage (Catalog *catalog, off_t meas, int cat) {
 
+  off_t idx, ID;
+  Measure *measure;
+
+  measure = &catalog[cat].measure[meas];
+
+  ID = measure[0].imageID;
+  idx = getImageByID (ID);
+  if (idx == -1) {
+    fprintf (stderr, "can't match detection to image?\n");
+    abort();
+  }
+
+  // index for (catalog, measure) -> image
+  bin[cat][meas] = idx;
+
+  // index for image, Nentry -> catalog
+  clist[idx][Nlist[idx]] = cat;
+
+  // index for image, Nentry -> measure
+  mlist[idx][Nlist[idx]] = meas;
+  Nlist[idx] ++;
+
+  if (Nlist[idx] == NLIST[idx]) {
+    NLIST[idx] += 100;
+    REALLOCATE (clist[idx], int,   NLIST[idx]);
+    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
+  }
+  return;
+}
+
+# else
+/* modify this function to use the measure->imageID field */
+void matchImage (Catalog *catalog, off_t meas, int cat) {
+
   off_t i;
   Measure *measure;
@@ -107,4 +188,6 @@
   /* find the image that supplied this measurement */
   for (i = 0; i < Nimage; i++) {
+    // let's try the very slow method first before adding a bisection search
+    // if (image[i].imageID != measure[0].imageID) continue;
     if (image[0].photcode == -1) continue;
     if (measure[0].photcode != image[i].photcode) continue;
@@ -129,5 +212,8 @@
     return;
   }
-}
+  fprintf (stderr, "can't match detection to image?\n");
+  abort();
+}
+# endif
 
 Coords *getCoords (off_t meas, int cat) {
@@ -199,5 +285,5 @@
   mosaic = getMosaicForImage (im);
   if (mosaic != NULL) {
-      moscoords = &mosaic[0].coords;
+    moscoords = &mosaic[0].coords;
   }
   imcoords = &image[im].coords;
@@ -270,10 +356,10 @@
   moscoords = NULL;
   if (mode == MODE_MOSAIC) {
-      mosaic = getMosaicForImage (im);
-      if (mosaic == NULL) {
-        fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
-        exit (1);
-      }
-      moscoords = &mosaic[0].coords;
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) {
+      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
+      exit (1);
+    }
+    moscoords = &mosaic[0].coords;
   }
 
@@ -319,7 +405,7 @@
         LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
         break;
-    default:
-      fprintf (stderr, "error: invalid mode in getImageRaw");
-      abort ();
+      default:
+	fprintf (stderr, "error: invalid mode in getImageRaw");
+	abort ();
     }
   }
@@ -371,14 +457,14 @@
     switch (mode) {
       case MODE_SIMPLE:
-      RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
-      ref[i].L = ref[i].P;
-      ref[i].M = ref[i].Q;
-      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
-      break;
+	RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
+	ref[i].L = ref[i].P;
+	ref[i].M = ref[i].Q;
+	LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
+	break;
       case MODE_MOSAIC:
-      RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
-      LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
-      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
-      break;
+	RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
+	LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
+	LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
+	break;
       default:
         fprintf (stderr, "invalid case");
Index: /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/MosaicOps.c
===================================================================
--- /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/MosaicOps.c	(revision 27424)
+++ /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/MosaicOps.c	(revision 27425)
@@ -19,9 +19,61 @@
 }
 
+off_t getMosaicByTimes (unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nmosaic;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (startMos[N] < start) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nmosaic);
+    }
+  }
+
+  // check for the matched mosaic starting from Nlo 
+  // we may have to go much beyond Nlo since stop is not sorted
+  // can we use a sorted version of stop to check when we are beyond the valid range??
+  for (N = Nlo; N < Nmosaic; N++) { 
+    if (stop  < stopMos[N]) continue;
+    if (start > startMos[N])  continue;
+    if (stop  < startMos[N]) return (-1);
+    return (indexMos[N]);
+  }
+
+  return (-1);
+}
+
+// sort two times vectors and an index by first time vector
+void sort_mosaic_times (unsigned int *S, unsigned int *E, off_t *I, off_t N) {
+
+# define SWAPFUNC(A,B){ unsigned int tmp_t; off_t tmp_i; \
+  tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \
+  tmp_t = E[A]; E[A] = E[B]; E[B] = tmp_t; \
+  tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \
+}
+# define COMPARE(A,B)(S[A] < S[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+// first, let's continue to use the time to make the match, but use bracketing to make it faster:
+
 // find mosaic frames (unique time periods & photcode name matches mosaic) 
 void initMosaics (Image *image, off_t Nimage) {
 
-  off_t i, j, found, NMOSAIC;
+  off_t i, Nmos, NMOSAIC;
   unsigned int start, stop;
+
+  unsigned int *startMos, *stopMos;
+  off_t *indexMos;
 
   Nmosaic = 0;
@@ -32,4 +84,7 @@
   ALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
   ALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+  ALLOCATE (startMos, unsigned int, NMOSAIC);
+  ALLOCATE (stopMos, unsigned int, NMOSAIC);
+  ALLOCATE (indexMos, off_t, NMOSAIC);
 
   /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */
@@ -56,4 +111,8 @@
     ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]);
 
+    startMos[Nmosaic] = start;
+    stopMos[Nmosaic] = stop;
+    indexMos[Nmosaic] = Nmosaic;
+
     Nmosaic ++;
     if (Nmosaic == NMOSAIC) {
@@ -63,7 +122,13 @@
       REALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
       REALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
-    }
-  }
-
+      REALLOCATE (startMos, unsigned int, NMOSAIC);
+      REALLOCATE (stopMos, unsigned int, NMOSAIC);
+      REALLOCATE (indexMos, off_t, NMOSAIC);
+    }
+  }
+
+  // sort the index, start, and stop by the start times:
+  sort_mosaic_times (startMos, stopMos, indexMos, Nmosaic);
+  
   // array to store image->mosaic index
   Nmosaic_for_images = Nimage;
@@ -80,26 +145,18 @@
     stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
 
-    /* find existing mosaic with this time range */
-    found = FALSE;
-    for (j = 0; !found && (j < Nmosaic); j++) { 
-      if (stop  < mosaic[j].start) continue;
-      if (start > mosaic[j].stop)  continue;
-      found = TRUE;
-      break;
-    }
-    /* if no matching mosaic exists, skip this image */
-    if (!found) continue;
+    Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos);
+    if (Nmos == -1) continue;
 
     // mosaic corresponding to this image
-    mosaic_for_images[i] = j;
+    mosaic_for_images[i] = Nmos;
 
     // add image to mosaic_own_image list 
-    mosaic_own_images[j][Nmosaic_own_images[j]] = i;
-    Nmosaic_own_images[j] ++;
-    if (Nmosaic_own_images[j] == Amosaic_own_images[j]) {
-      Amosaic_own_images[j] += 10;
-      REALLOCATE (mosaic_own_images[j], off_t, Amosaic_own_images[j]);
-    }
-    assert (Nmosaic_own_images[j] < Amosaic_own_images[j]);
+    mosaic_own_images[Nmos][Nmosaic_own_images[Nmos]] = i;
+    Nmosaic_own_images[Nmos] ++;
+    if (Nmosaic_own_images[Nmos] == Amosaic_own_images[Nmos]) {
+      Amosaic_own_images[Nmos] += 10;
+      REALLOCATE (mosaic_own_images[Nmos], off_t, Amosaic_own_images[Nmos]);
+    }
+    assert (Nmosaic_own_images[Nmos] < Amosaic_own_images[Nmos]);
   }
 
Index: /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/load_images.c
===================================================================
--- /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/load_images.c	(revision 27424)
+++ /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/load_images.c	(revision 27425)
@@ -1,3 +1,9 @@
 # include "relastro.h"
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 SkyList *load_images (FITS_DB *db, SkyRegion *region) {
@@ -6,7 +12,10 @@
   off_t      Nimage, Nsubset;
   off_t     *LineNumber;
+  struct timeval start, stop;
 
   SkyTable *sky = NULL;
   SkyList *skylist = NULL;
+
+  gettimeofday (&start, (void *) NULL);
 
   // load the current sky table (layout of all SkyRegions) 
@@ -16,15 +25,21 @@
   // determine the populated SkyRegions overlapping the requested area
   skylist = SkyListByPatch (sky, -1, region);
+  MARKTIME("  setup sky: %f sec\n", dtime);
 
   // convert database table to internal structure
   image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
+  MARKTIME("  convert image table: %f sec\n", dtime);
 
   // select the images which overlap the selected sky regions
   subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
+  MARKTIME("  select images: %f sec\n", dtime);
 
   gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
 
   initImages (subset, Nsubset);
+  MARKTIME("  init images: %f sec\n", dtime);
+
   initMosaics (subset, Nsubset);
+  MARKTIME("  init mosaics: %f sec\n", dtime);
   
   /* unlock, if we can (else, unlocked below) */
Index: /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/relastro.c
===================================================================
--- /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/relastro.c	(revision 27424)
+++ /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/relastro.c	(revision 27425)
@@ -1,3 +1,9 @@
 # include "relastro.h"
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 int main (int argc, char **argv) {
@@ -6,6 +12,9 @@
   Catalog *catalog;
   FITS_DB db;
+  struct timeval start, stop;
 
   SkyList *skylist = NULL;
+
+  gettimeofday (&start, (void *) NULL);
 
   /* get configuration info, args */
@@ -26,14 +35,20 @@
   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);
+  MARKTIME("load image data: %f sec\n", dtime);
 
   /* load regions and images based on specified sky patch (default depth) */
   skylist = load_images (&db, &UserPatch);
+  MARKTIME("load images: %f sec\n", dtime);
 
   /* load catalog data from region files : subselect high-quality measurements */
   catalog = load_catalogs (skylist, &Ncatalog, TRUE);
+  MARKTIME("load catalog data: %f sec\n", dtime);
 
   /* match measurements with images */
   initImageBins (catalog, Ncatalog);
+  MARKTIME("make image bins: %f sec\n", dtime);
+
   findImages (catalog, Ncatalog);
+  MARKTIME("set up image indexes: %f sec\n", dtime);
 
   if (PLOTSTUFF) {
Index: /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/select_images.c
===================================================================
--- /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/select_images.c	(revision 27424)
+++ /branches/eam_branches/largefiles.20100314/Ohana/src/relastro/src/select_images.c	(revision 27425)
@@ -13,8 +13,11 @@
 } SkyRegionCoords;
 
+void dsortindex (double *X, off_t *Y, int N);
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
   
   Image *image;
-  off_t i, j, k, m, nimage, NIMAGE;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
   int InRange, found;
   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
@@ -23,4 +26,7 @@
   SkyRegionCoords *skycoords;
   
+  double *RmaxSky;
+  off_t *index;
+
   if (skylist[0].Nregions < 1) {
     *Nimage = 0;
@@ -38,4 +44,7 @@
 
   ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
+
+  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
+  ALLOCATE (index, off_t, skylist[0].Nregions);
 
   /* compare with each region file */
@@ -56,4 +65,7 @@
     skycoords[i].Yc[4] = skycoords[i].Yc[0];    
 
+    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
+    index[i] = i;
+
     dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
     dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
@@ -64,4 +76,6 @@
     skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
   }
+
+  dsortindex (RmaxSky, index, skylist[0].Nregions);
 
   if (VERBOSE) fprintf (stderr, "finding images\n");
@@ -100,5 +114,8 @@
     }
     
-    if (!FindMosaicForImage (timage, Ntimage, i)) continue;
+    if (!FindMosaicForImage (timage, Ntimage, i)) {
+      fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i);
+      continue;
+    }
 
     /* define image corners */
@@ -111,10 +128,22 @@
 
     /* transform to ra,dec */
+    double RminImage = 360.0;
+    // double RmaxImage =   0.0;
     for (j = 0; j < 5; j++) {
       XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
-    }
+      RminImage = MIN(RminImage, Ri[j]);
+      // RmaxImage = MIN(RmaxImage, Ri[j]);
+    }
+
+    // RA(nStart) is guaranteed to be < RminImage:
+    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
 
     /* compare with each region file */
-    for (m = 0; (m < skylist[0].Nregions) && !found; m++) { 
+    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) { 
+
+      m = index[iSky];
+      // if (RmaxImage < skylist[0].regions[m][0].Rmin) {
+      // break;
+      // }
 
       /* we make positional comparisons in the projection of catalog */
@@ -122,5 +151,5 @@
       tcoords.crval2 = skycoords[m].Dc;
 
-      /* transform to ra,dec */
+      /* transform corner coords to X,Y in this catalog system */
       InRange = TRUE;
       for (j = 0; (j < 5) && InRange; j++) {
@@ -253,2 +282,35 @@
 }
 
+void dsortindex (double *X, off_t *Y, int N) {
+
+# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \
+  tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \
+  tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \
+}
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nregions;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (Rref[N] < R) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N, Nregions);
+    }
+  }
+  return (Nlo);
+}