Index: trunk/Ohana/src/relastro/doc/notes.txt
===================================================================
--- trunk/Ohana/src/relastro/doc/notes.txt	(revision 25757)
+++ trunk/Ohana/src/relastro/doc/notes.txt	(revision 27435)
@@ -1,2 +1,61 @@
+
+2010.03.24
+
+  I have finished only part of the optimizations listed below: I have
+  made the detection->image and image->mosaic lookups use bracket
+  searches.  I have NOT converted the image->mosaic lookups to use the
+  parentID, so these matches are still based on time ranges.  I have
+  also added a bit of positional constraint to the image match
+  function in 'select_images' (comparing Rmax of the sky regions to
+  Rmin of the images).  These have had a huge improvement.  Here is
+  the execution time for the SAS data (~450 images) with the old vs
+  the new code on pikake:
+
+  old : 39min
+  new :  7min
+
+  The bulk of the time for the new code is now disk I/O.
+
+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: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 25757)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 27435)
@@ -140,5 +140,5 @@
 char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
 void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
-int          *SelectRefMosaic     PROTO((Mosaic **refmosaic, int *Nimage));
+off_t          *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
 int           args                PROTO((int argc, char **argv));
 int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
@@ -153,12 +153,12 @@
 void          findImages          PROTO((Catalog *catalog, int Ncatalog));
 int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
-Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, int Nregion, int *Nimage, int **LineNum));
+Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum));
 void set_db (FITS_DB *in);
-int Shutdown (char *format, ...);
+int Shutdown (char *format, ...) OHANA_FORMAT(printf, 1, 2);
 void TrapSignal (int sig);
 void SetProtect (int mode);
 int SetSignals (void);
 
-GSCRegion    *find_regions        PROTO((Image *image, int Nimage, int *Nregions, GSCRegion *fullregion));
+GSCRegion    *find_regions        PROTO((Image *image, off_t Nimage, int *Nregions, GSCRegion *fullregion));
 void          freeGridBins        PROTO((int Ncatalog));
 void          freeImageBins       PROTO((int Ncatalog));
@@ -166,21 +166,21 @@
 void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
 int           gcatalog            PROTO((Catalog *catalog, int FINAL));
-Coords       *getCoords           PROTO((int meas, int cat));
-float         getMcal             PROTO((int meas, int cat));
-float         getMgrid            PROTO((int meas, int cat));
-float         getMmos             PROTO((int meas, int cat));
-float         getMrel             PROTO((Catalog *catalog, int meas, int cat));
-GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, int *Nregions));
-void          getfullregion       PROTO((Image *image, int Nimage, GSCRegion *fullregion));
-Image        *getimage            PROTO((int N));
-Image        *getimages           PROTO((int *N));
+Coords       *getCoords           PROTO((off_t meas, int cat));
+float         getMcal             PROTO((off_t meas, int cat));
+float         getMgrid            PROTO((off_t meas, int cat));
+float         getMmos             PROTO((off_t meas, int cat));
+float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
+GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, off_t *Nregions));
+void          getfullregion       PROTO((Image *image, off_t Nimage, GSCRegion *fullregion));
+Image        *getimage            PROTO((off_t N));
+Image        *getimages           PROTO((off_t *N));
 void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
 void          initGrid            PROTO((int dX, int dY));
 void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
 void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
-void          initImages          PROTO((Image *input, int N));
+void          initImages          PROTO((Image *input, off_t N));
 void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
-void          initMosaicGrid      PROTO((Image *image, int Nimage));
-void          initMosaics         PROTO((Image *image, int Nimage));
+void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
+void          initMosaics         PROTO((Image *image, off_t Nimage));
 void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
 void          initialize          PROTO((int argc, char **argv));
@@ -189,5 +189,5 @@
 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
 SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
-Image        *select_images       PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage));
+Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
 
 void check_permissions (char *basefile);
@@ -198,8 +198,8 @@
 
 int           main                PROTO((int argc, char **argv));
-void          mark_images         PROTO((Image *image, int Nimage, Image *timage, int Ntimage));
-void          matchImage          PROTO((Catalog *catalog, int meas, int cat));
-void          matchMosaics        PROTO((Catalog *catalog, int meas, int cat));
-GSCRegion    *name_region         PROTO((char *name, int *Nregions));
+void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
+void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat));
+void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat));
+GSCRegion    *name_region         PROTO((char *name, off_t *Nregions));
 double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
 void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
@@ -224,5 +224,5 @@
 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
 void          set_ZP              PROTO((double ZERO));
-int           setrefcode          PROTO((Image *image, int Nimage));
+int           setrefcode          PROTO((Image *image, off_t Nimage));
 void          skip_measurements   PROTO((Catalog *catalog, int pass));
 void          sortA               PROTO((double *X, int N));
@@ -268,5 +268,5 @@
 int UpdateMosaic (Catalog *catalog, int Ncatalog);
 int UpdateMeasures (Catalog *catalog, int Ncatalog);
-void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
+void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im);
 void FlagOutliers(Catalog *catalog);
 int MeasFilterTest(Measure *measure);
@@ -278,14 +278,14 @@
 int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);
 
-Mosaic *getMosaicForImage (int N);
-
-StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode);
-StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode);
-
-Mosaic *getmosaics (int *N);
-void initMosaics (Image *image, int Nimage);
-StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars);
-StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars);
-Mosaic *getMosaicForImage (int im);
+Mosaic *getMosaicForImage (off_t N);
+
+StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
+
+Mosaic *getmosaics (off_t *N);
+void initMosaics (Image *image, off_t Nimage);
+StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
+StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
+Mosaic *getMosaicForImage (off_t im);
 
 double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
@@ -297,2 +297,8 @@
 int relastro_objects (void);
 int UpdateObjectOffsets (SkyList *skylist);
+
+int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
+int relastroVisualPlotScatter(double values[], double thresh, int npts);
+int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, 
+			       StatType statsR, StatType statsD, double thresh);
+
Index: trunk/Ohana/src/relastro/src/FitSimple.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitSimple.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/FitSimple.c	(revision 27435)
@@ -17,5 +17,5 @@
 
     // XXX the value needs to be set in a more intelligent way
-    if (dR > 0.15) continue;
+    if (dR > 0.50) continue;
     
     fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0);
Index: trunk/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/ImageOps.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/ImageOps.c	(revision 27435)
@@ -1,35 +1,52 @@
 # include "relastro.h"
-#define TESTING
-
+# define TESTING
+# define USE_IMAGE_ID 1
+
+static off_t       **bin;     // link from catalog,measure to image
+static int         **clist;   // catalog which supplied measurement on image
+static off_t       **mlist;   // measure reference for measurement on image
+static off_t        *Nlist;   // number of measurements on image
+static off_t        *NLIST;   // allocated number of measurements on image
+
+static Image        *image;   // list of available images
+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;
-static int         **bin;
-
-static int         **clist;
-static int         **mlist;
-static int          *Nlist;
-static int          *NLIST;
-
-static Image        *image;
-static int          Nimage;
-
-Image *getimages (int *N) {
+# endif
+
+Image *getimages (off_t *N) {
   *N = Nimage;
   return (image);
 }
 
-Image *getimage (int N) {
+Image *getimage (off_t N) {
   return (&image[N]);
 }
 
-void initImages (Image *input, int N) {
-
-  int i;
+void initImages (Image *input, off_t N) {
+
+  off_t i;
 
   image = input;
   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,20 +54,49 @@
     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);
 }
 
 void initImageBins (Catalog *catalog, int Ncatalog) {
 
-  int i, j;
-
-  ALLOCATE (bin, int *, Ncatalog);
+  off_t i, j;
+
+  ALLOCATE (bin, off_t *, Ncatalog);
   for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
+    ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
     for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
   }
 
-  ALLOCATE (Nlist, int, Nimage);
-  ALLOCATE (NLIST, int, Nimage);
-  ALLOCATE (clist, int *, Nimage);
-  ALLOCATE (mlist, int *, Nimage);
+  ALLOCATE (Nlist, off_t,   Nimage);
+  ALLOCATE (NLIST, off_t,   Nimage);
+  ALLOCATE (clist, int *,   Nimage);
+  ALLOCATE (mlist, off_t *, Nimage);
 
   for (i = 0; i < Nimage; i++) {
@@ -58,5 +104,5 @@
     NLIST[i] = 100;
     ALLOCATE (clist[i], int, NLIST[i]);
-    ALLOCATE (mlist[i], int, NLIST[i]);
+    ALLOCATE (mlist[i], off_t, NLIST[i]);
   }
 }
@@ -64,5 +110,5 @@
 void freeImageBins (int Ncatalog) {
 
-  int i;
+  off_t i;
 
   for (i = 0; i < Ncatalog; i++) {
@@ -81,5 +127,5 @@
 void findImages (Catalog *catalog, int Ncatalog) {
 
-  int i, j;
+  off_t i, j;
   char *name;
 
@@ -92,13 +138,48 @@
   for (i = 0; VERBOSE && (i < Nimage); i++) {
     name = GetPhotcodeNamebyCode (image[i].photcode);
-    fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i],
+    fprintf (stderr, "image %lld has %lld measures (%s, %s)\n", (long long) i, (long long) Nlist[i],
              ohana_sec_to_date(image[i].tzero), name);
   }
 }
 
+# if USE_IMAGE_ID
 /* modify this function to use the measure->imageID field */
-void matchImage (Catalog *catalog, int meas, int cat) {
-
-  int i;
+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;
@@ -124,14 +207,17 @@
     if (Nlist[i] == NLIST[i]) {
       NLIST[i] += 100;
-      REALLOCATE (clist[i], int, NLIST[i]);
-      REALLOCATE (mlist[i], int, NLIST[i]);
+      REALLOCATE (clist[i], int,   NLIST[i]);
+      REALLOCATE (mlist[i], off_t, NLIST[i]);
     }
     return;
   }
-}
-
-Coords *getCoords (int meas, int cat) {
-
-  int i;
+  fprintf (stderr, "can't match detection to image?\n");
+  abort();
+}
+# endif
+
+Coords *getCoords (off_t meas, int cat) {
+
+  off_t i;
 
   i = bin[cat][meas];
@@ -142,5 +228,5 @@
 void plot_images () {
 
-  int i, bin;
+  off_t i, bin;
   double *xlist, *Mlist, *dlist;
   Graphdata graphdata;
@@ -188,7 +274,7 @@
 // return StarData values for detections in the specified image, converting coordinates from the
 // chip positions: X,Y -> L,M -> P,Q -> R,D
-void fixImageRaw (Catalog *catalog, int Ncatalog, int im) {
-
-  int i, m, c, n;
+void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
+
+  off_t i, m, c, n;
   double X, Y, L, M, P, Q, R, D, dR, dD;
 
@@ -199,5 +285,5 @@
   mosaic = getMosaicForImage (im);
   if (mosaic != NULL) {
-      moscoords = &mosaic[0].coords;
+    moscoords = &mosaic[0].coords;
   }
   imcoords = &image[im].coords;
@@ -257,7 +343,7 @@
 // chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
 // which the detections have already been filtered when they were loaded (bcatalog)
-StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
-
-  int i, m, c, n;
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
+
+  off_t i, m, c, n;
 
   Mosaic *mosaic;
@@ -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 ();
     }
   }
@@ -332,7 +418,7 @@
 // the sky positions: R,D -> P,Q -> L,M -> X,Y
 
-StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
-
-  int i, m, c, n;
+StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
+
+  off_t i, m, c, n;
 
   Mosaic *mosaic;
@@ -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");
@@ -399,5 +485,6 @@
   return;
 
-  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
+  int Ndel, Nave;
+  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
   double Ns, theta, x, y;
   double *R, *D, *dR, *dD;
@@ -490,5 +577,6 @@
 void FlagOutliers2D (Catalog *catalog) {
 
-  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
+  int Ndel, Nave;
+  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
   double *index;
   double Ns, theta, x, y;
@@ -569,5 +657,5 @@
     // recalculate image center, sigma based on closest 50% of points
     for(k = 0;  k < N; k++) {
-      int ind = (int) index[k];
+      off_t ind = (off_t) index[k];
       R[k] = catalog[0].measure[ind].dR;
       D[k] = catalog[0].measure[ind].dD;
Index: trunk/Ohana/src/relastro/src/MosaicOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 27435)
@@ -2,26 +2,78 @@
 
 // array of mosaic definition structures
-static int    Nmosaic;
+static off_t    Nmosaic;
 static Mosaic *mosaic;
 
 // list of all images associated with a mosaic
-static int   *Nmosaic_own_images; // number of images for this mosaic
-static int   *Amosaic_own_images; // size of allocated array
-static int   **mosaic_own_images; // array of arrays: mosaic -> images
+static off_t   *Nmosaic_own_images; // number of images for this mosaic
+static off_t   *Amosaic_own_images; // size of allocated array
+static off_t   **mosaic_own_images; // array of arrays: mosaic -> images
 
 // list of mosaic associated with each image  
-static int    Nmosaic_for_images; // number of images (for internal checks)
-static int    *mosaic_for_images; // array of: image -> mosaic
-
-Mosaic *getmosaics (int *N) {
+static off_t    Nmosaic_for_images; // number of images (for off_ternal checks)
+static off_t    *mosaic_for_images; // array of: image -> mosaic
+
+Mosaic *getmosaics (off_t *N) {
   *N = Nmosaic;
   return (mosaic);
 }
 
+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, int Nimage) {
-
-  int i, j, found, NMOSAIC;
+void initMosaics (Image *image, off_t Nimage) {
+
+  off_t i, Nmos, NMOSAIC;
   unsigned int start, stop;
+
+  unsigned int *startMos, *stopMos;
+  off_t *indexMos;
 
   Nmosaic = 0;
@@ -29,7 +81,10 @@
   ALLOCATE (mosaic, Mosaic, NMOSAIC);
 
-  ALLOCATE (Nmosaic_own_images, int, NMOSAIC);
-  ALLOCATE (Amosaic_own_images, int, NMOSAIC);
-  ALLOCATE (mosaic_own_images, int *, NMOSAIC);
+  ALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
+  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 */
@@ -54,5 +109,9 @@
     Nmosaic_own_images[Nmosaic] = 0;
     Amosaic_own_images[Nmosaic] = 10;
-    ALLOCATE (mosaic_own_images[Nmosaic], int, Amosaic_own_images[Nmosaic]);
+    ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]);
+
+    startMos[Nmosaic] = start;
+    stopMos[Nmosaic] = stop;
+    indexMos[Nmosaic] = Nmosaic;
 
     Nmosaic ++;
@@ -60,13 +119,19 @@
       NMOSAIC += 10;
       REALLOCATE (mosaic, Mosaic, NMOSAIC);
-      REALLOCATE (mosaic_own_images, int *, NMOSAIC);
-      REALLOCATE (Nmosaic_own_images, int, NMOSAIC);
-      REALLOCATE (Amosaic_own_images, int, NMOSAIC);
-    }
-  }
-
+      REALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+      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;
-  ALLOCATE (mosaic_for_images, int, Nmosaic_for_images);
+  ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images);
 
   /* now assign the WRP images to these mosaics */
@@ -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], int, 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]);
   }
 
@@ -109,7 +166,7 @@
 // return StarData values for detections in the specified image, converting coordinates from the
 // chip positions: X,Y -> L,M -> P,Q -> R,D
-StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
-
-  int i, j, im, Nraw, Nnew;
+StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nraw, Nnew;
   StarData *raw, *new;
 
@@ -142,7 +199,7 @@
 // return StarData values for averages positions in the specified image, converting coordinates from
 // the sky positions: R,D -> P,Q -> L,M -> X,Y
-StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
-
-  int i, j, im, Nref, Nnew;
+StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nref, Nnew;
   StarData *ref, *new;
   
@@ -172,7 +229,7 @@
 }
 
-Mosaic *getMosaicForImage (int im) {
-
-  int mos;
+Mosaic *getMosaicForImage (off_t im) {
+
+  off_t mos;
 
   if (im < 0) abort();
Index: trunk/Ohana/src/relastro/src/UpdateChips.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 27435)
@@ -4,5 +4,5 @@
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  int i, Nimage, Nraw, Nref;
+  off_t i, Nimage, Nraw, Nref;
   Image *image;
   StarData *raw, *ref;
@@ -25,5 +25,5 @@
 
     // FitChip does iterative, clipped fitting
-    fprintf (stderr, "image %d : Nstars: %d\n", i, Nraw);
+    fprintf (stderr, "image %lld : Nstars: %lld\n", (long long) i, (long long) Nraw);
     FitChip (raw, ref, Nraw, &image[i].coords);
 
Index: trunk/Ohana/src/relastro/src/UpdateMeasures.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateMeasures.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/UpdateMeasures.c	(revision 27435)
@@ -3,5 +3,5 @@
 int UpdateMeasures (Catalog *catalog, int Ncatalog) {
 
-  int i, Nimage;
+  off_t i, Nimage;
   Image *image;
 
Index: trunk/Ohana/src/relastro/src/UpdateMosaic.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateMosaic.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/UpdateMosaic.c	(revision 27435)
@@ -4,5 +4,5 @@
 
   /* we can measure new image parameters for each mosaic independently */
-  int i, Nmosaic, Nstars;
+  off_t i, Nmosaic, Nstars;
   Mosaic *mosaic;
   StarData *raw, *ref;
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 27435)
@@ -1,5 +1,5 @@
 # include "relastro.h"
 
-static int Nmax;
+static off_t   Nmax;
 static double *X, *dX;
 static double *Y, *dY;
@@ -13,5 +13,5 @@
 void initObjectData (Catalog *catalog, int Ncatalog) {
 
-  int i, j;
+  off_t i, j;
   
   Nmax = 0;
@@ -40,5 +40,5 @@
 int UpdateObjects (Catalog *catalog, int Ncatalog) {
 
-  int i, j, k, m, N, Nsecfilt, found, kp;
+  int i, j, k, m, N, Nsecfilt;
   StatType statsR, statsD;
   Coords coords;
@@ -47,7 +47,4 @@
   int mode, Nave, Npm, Npar, Nskip;
   double Tmin, Tmax;
-  float mag;
-  int mask;
-  PhotCode *code;
 
   initObjectData (catalog, Ncatalog);
@@ -74,5 +71,5 @@
   for (i = 0; i < Ncatalog; i++) {
 
-    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure);
+    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %lld ave, %lld meas\n", i, (long long) catalog[i].Naverage, (long long) catalog[i].Nmeasure);
 
     Nskip = 0;
Index: trunk/Ohana/src/relastro/src/UpdateSimple.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateSimple.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/UpdateSimple.c	(revision 27435)
@@ -4,5 +4,5 @@
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  int i, Nimage, Nstars;
+  off_t i, Nimage, Nstars;
   Image *image;
   StarData *raw, *ref;
Index: trunk/Ohana/src/relastro/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relastro/src/bcatalog.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/bcatalog.c	(revision 27435)
@@ -3,9 +3,7 @@
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
 
-  int i, j, k, offset, found;
-  int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
-  float mag;
-  int mask;
-  PhotCode *code;
+  off_t i, j, offset;
+  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
+  int Nsecfilt;
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -106,6 +104,5 @@
 
   if (VERBOSE) {
-    fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,
-             subcatalog[0].Naverage, subcatalog[0].Nmeasure);
+    fprintf (stderr, "%lld: using %lld stars (%lld measures) for catalog\n", (long long) i, (long long) subcatalog[0].Naverage, (long long) subcatalog[0].Nmeasure);
    }
   return (TRUE);
Index: trunk/Ohana/src/relastro/src/load_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/load_images.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/load_images.c	(revision 27435)
@@ -1,12 +1,21 @@
 # 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) {
 
   Image     *image, *subset;
-  int        Nimage, Nsubset;
-  int       *LineNumber;
+  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: trunk/Ohana/src/relastro/src/relastro.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/relastro.c	(revision 27435)
@@ -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: trunk/Ohana/src/relastro/src/select_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/select_images.c	(revision 25757)
+++ trunk/Ohana/src/relastro/src/select_images.c	(revision 27435)
@@ -13,14 +13,20 @@
 } SkyRegionCoords;
 
-Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) {
+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;
-  int i, j, k, m, found, nimage, NIMAGE;
-  int InRange;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
+  int InRange, found;
   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
-  int *line_number;
+  off_t *line_number;
   Coords tcoords;
   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]);
@@ -65,4 +77,6 @@
   }
 
+  dsortindex (RmaxSky, index, skylist[0].Nregions);
+
   if (VERBOSE) fprintf (stderr, "finding images\n");
   BuildChipMatch (timage, Ntimage);
@@ -71,5 +85,5 @@
   NIMAGE = 100;
   ALLOCATE (image, Image, NIMAGE);
-  ALLOCATE (line_number, int, NIMAGE);
+  ALLOCATE (line_number, off_t, NIMAGE);
   
   // go through the complete list of images, selecting ones which overlap any region
@@ -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++) {
@@ -161,13 +190,13 @@
 	NIMAGE += 100;
 	REALLOCATE (image, Image, NIMAGE);
-	REALLOCATE (line_number, int, NIMAGE);
+	REALLOCATE (line_number, off_t, NIMAGE);
       }
     }
   }
       
-  if (VERBOSE) fprintf (stderr, "found %d images\n", nimage);
+  if (VERBOSE) fprintf (stderr, "found %lld images\n", (long long) nimage);
 
   REALLOCATE (image, Image, MAX (nimage, 1));
-  REALLOCATE (line_number, int, MAX (nimage, 1));
+  REALLOCATE (line_number, off_t, MAX (nimage, 1));
   free (skycoords);
 
@@ -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);
+}