Changeset 39488 for trunk

Timestamp:

Mar 28, 2016, 10:07:15 AM (10 years ago)

Author:

eugene

Message:

running slow, trying to speed things up

Location:

trunk/Ohana/src/dvolens/src

Files:

• WarpImageMaps.c (modified) (4 diffs)
• myIndex.c (modified) (1 diff)
• select_images.c (modified) (1 diff)
• update_objects_catalog.c (modified) (7 diffs)

trunk/Ohana/src/dvolens/src/WarpImageMaps.c

@@ -11 +11 @@
   off_t *warps;
   unsigned short photcode;
+  double *Rmin;
+  double *Rmax;
+  double *Dmin;
+  double *Dmax;
+  int    *onBoundary;
 } WarpGroup;
 
@@ -112 +117 @@
   }
 
+  // we now have the warps assigned to groups.  now we need to generate the grid data to find the warp assignments
+  for (int i = 0; i < Nwarpgroup; i++) {
+
+    myAssert (warpgroup[i].Nwarps == uniqcount[i], "failure");
+    ALLOCATE (warpgroup[i].Rmin, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].Rmax, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].onBoundary, int, warpgroup[i].Nwarps);
+
+    for (int j = 0; j < warpgroup[i].Nwarps; j++) {
+
+      off_t N = warpgroup[i].warps[j];
+      
+      int Nx = image[N].NX;
+      int Ny = image[N].NY;
+
+      warpgroup[i].onBoundary[j] = FALSE;
+
+      double R, D;
+      XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = R;
+      warpgroup[i].Rmax[j] = R;
+      warpgroup[i].Dmin[j] = D;
+      warpgroup[i].Dmax[j] = D;
+
+      // XXX need to worry about 0,360 boundary
+      XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
+      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
+      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
+      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
+
+      // XXX need to worry about 0,360 boundary
+      XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
+      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
+      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
+      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
+
+      // XXX need to worry about 0,360 boundary
+      XY_to_RD (&R, &D, Nx, Ny, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
+      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
+      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
+      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
+
+      // bump Rmin,Rmax,Dmin,Dmax 10 arcsec worth of padding
+      double dR = 10.0/3600.0 / cos (RAD_DEG*warpgroup[i].Dmin[j]);
+      double dD = 10.0/3600.0;
+
+      // at north pole, force test of all nearby skycells
+      if (warpgroup[i].Dmax[j] > 89.8) {
+        warpgroup[i].Rmin[j] =   0.0;
+        warpgroup[i].Rmax[j] = 360.0;
+        warpgroup[i].Dmax[j] =  90.0;
+        warpgroup[i].Dmin[j] -=  dD;
+        continue;
+      } 
+
+      if (warpgroup[i].Rmax[j] - warpgroup[i].Rmin[j] > 270.0) {
+        // Rmin and Rmax are in the range 0 - 360.  For images at the 0,360 boundary,
+        // "Rmax" is the lower edge, and "Rmin" is the upper edge.  we need to flip them 
+        // and then break the 0-360 range:
+        double tmp = warpgroup[i].Rmin[j];
+        warpgroup[i].Rmin[j] = warpgroup[i].Rmax[j] - 360.0;
+        warpgroup[i].Rmax[j] = tmp;
+        warpgroup[i].onBoundary[j] = TRUE;
+      }
+
+      warpgroup[i].Rmin[j] -=  dR;
+      warpgroup[i].Rmax[j] +=  dR;
+      warpgroup[i].Dmin[j] -=  dD;
+      warpgroup[i].Dmax[j] +=  dD;
+    }
+  }
+
   free (uniqtimes); 
   free (uniqcount); 
@@ -131 +217 @@
 
   // can we find an image in this set which matches one of our measures?
+  double Rave = average->R;
+  double Dave = average->D;
 
   // we now have the warp group, but which is the correct warp?
@@ -138 +226 @@
     // average->R,D could pin-point more quickly
 
-    // this is a possible image: check the coordinates:
+    if (warpgroup[seq].onBoundary[i]) {
+      int inRange1 = (Rave >= warpgroup[seq].Rmin[i]) && (Rave <= warpgroup[seq].Rmax[i]);
+      int inRange2 = (Rave >= warpgroup[seq].Rmin[i] + 360.0) && (Rave <= warpgroup[seq].Rmax[i] + 360.0);
+      if (!inRange1 && !inRange2) continue;
+    } else {
+      if (Rave < warpgroup[seq].Rmin[i]) continue;
+      if (Rave > warpgroup[seq].Rmax[i]) continue;
+    }
+    if (Dave < warpgroup[seq].Dmin[i]) continue;
+    if (Dave > warpgroup[seq].Dmax[i]) continue;
+
+    // this is a possible image based on coordinates: does it match a warp?
     int N = warpgroup[seq].warps[i];
 

trunk/Ohana/src/dvolens/src/myIndex.c

@@ -47 +47 @@
 
   myIndex->Nindex = myIndex->maxID - myIndex->minID + 1;
-  myIndex->NINDEX = myIndex->Nindex;
-  
-  REALLOCATE (myIndex->index, int, myIndex->NINDEX);
+  if (myIndex->Nindex > myIndex->NINDEX) {
+    myIndex->NINDEX = myIndex->Nindex;
+    REALLOCATE (myIndex->index, int, myIndex->NINDEX);
+  }
 
   int i;

trunk/Ohana/src/dvolens/src/select_images.c

@@ -43 +43 @@
   ALLOCATE (image, Image, NIMAGE);
   
+  // I need to ensure we get the images from any skycells in a warp
+  double dD = 2.0;
+  double dR = dD / cos(RAD_DEG*DminSkyRegion);
+  RminSkyRegion -= dR;
+  RmaxSkyRegion += dR;
+  DminSkyRegion -= dD;
+  DmaxSkyRegion += dD;
+
   // go through the complete list of images, selecting ones which overlap any region
   for (i = 0; i < Ntimage; i++) {

trunk/Ohana/src/dvolens/src/update_objects_catalog.c

@@ -2 +2 @@
 # define SCALE 0.001
 float MagToFlux (float Mag); // in libdvo, but not exposed?
+
+# define TIMESTAMP(TIME){                               \
+  gettimeofday (&stopTimer, (void *) NULL);             \
+  double dtime = DTIME (stopTimer, startTimer);         \
+  TIME += dtime;                                        \
+  gettimeofday (&startTimer, (void *) NULL); }
 
 int update_objects_catalog (Catalog *catalog) {
@@ -31 +37 @@
   myIndexType *measureIndex = myIndexAlloc();
 
+  int Nfixed = 0;
+  int Nfailed = 0;
+
+  INITTIME;
+
+  double time1 = 0.0;
+  double time2 = 0.0;
+  double time3 = 0.0;
+  double time4 = 0.0;
+
   for (i = 0; i < catalog->Naverage; i++) {
+    
+    if (i % 1000 == 0) {
+      MARKTIME ("done with %d aves: %f sec (%f %f %f %f)\n", (int) i, dtime, time1, time2, time3, time4);
+      time1 = 0.0;
+      time2 = 0.0;
+      time3 = 0.0;
+      time4 = 0.0;
+      gettimeofday (&startTimer, NULL);
+    }
+
     Average *average = &catalog->average[i];
     if (average->Nlensing == 0) continue;
@@ -62 +88 @@
       myIndexUpdateLimits (measureIndex, catalog->measure[Moff + Mj].imageID);
     }
+
+    TIMESTAMP (time1);
     myIndexSetRange (measureIndex);
+    TIMESTAMP (time2);
+
     for (Mj = 0; Mj < average->Nmeasure; Mj++) {
       // if there are duplicates (multiple measures from the same image), 
@@ -89 +119 @@
       // 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");
+      if (Mj < 0) {
+        if (REPAIR_LENSING_IDS) {
+          Mj = RecoverLensingIndex (average, measureIndex, lensing);
+          if (Mj == -1) {
+            Nfailed ++;
+            continue;
+          }
+          Nfixed ++;
+        } else {
+          Nfailed ++;
+          continue;
+        }
+      }
 
       Measure *measure = &catalog->measure[Moff + Mj];
@@ -196 +235 @@
     }
     
+    TIMESTAMP (time3);
+
     for (j = 0; j < Nsecfilt; j++) {
       if (!lensobj[j].Nmeas) {
@@ -273 +314 @@
     }
 
+    TIMESTAMP (time4);
+
     average->Nlensobj = Nsecfilt;
     average->lensobjOffset = Nlensobj;
@@ -290 +333 @@
   myIndexFree (measureIndex);
 
+  fprintf (stderr, "corrected %d lensing IDs, failed on %d lensing IDs\n", Nfixed, Nfailed);
+
   catalog->Nlensobj = Nlensobj;
   catalog->Nlensobj_disk = 0;