Changeset 33652 for trunk/Ohana/src/relastro/src/UpdateChips.c

Timestamp:

Apr 1, 2012, 3:00:19 PM (14 years ago)

Author:

eugene

Message:

updates to relastro to support distributed dvo; some optmimization (eg, image matching to area); cleanup names of static helper arrays in ImageOps, MosaicOps; cleanup rules for -reset options (relastro does not reset photom bits); skip objects with NAN error bars (can be used to exclude some photcodes from astrometry); enum for mask bits in FitChip and related; choose max order using same rules as psastro; plug leaks; reduce verbosity; big re-org of main (for parallel code); only use accepted measurements for Tmin,Tmax,Trange,Tmean; require a min parallax factor to fit parallaxes; use MeasureTiny where possible; use new ParsePhotcodeList function; testparallax function; threaded version of UpdateChips; bcatalog limit density selecting objects by number of measurements; high-speed now copies support tables to output database; high-speed may write to a parallel database; test for unsorted database (& fail if true)

File:

• trunk/Ohana/src/relastro/src/UpdateChips.c (modified) (7 diffs)

trunk/Ohana/src/relastro/src/UpdateChips.c

@@ -1 +1 @@
 # include "relastro.h"
+
 int plotChipFits (double *Ro, double *Do, char *mode, int Nimage);
 int saveCenter (Image *image, double *Ro, double *Do, int im);
-
+off_t getNextImageForThread ();
+int UpdateChips_threaded (Catalog *catalog, int Ncatalog);
+void *UpdateChips_worker (void *data);
+
+enum {THREAD_RUN, THREAD_DONE};
+
+typedef struct {
+  int entry;
+  int state;
+  double *Ro;
+  double *Do;
+  char *mode;
+  off_t Nskip;
+  off_t Nmosaic;
+  off_t NnewFit;
+  off_t NoldFit;
+  Catalog *catalog;
+  int Ncatalog;
+} ThreadInfo;
+
+static Image *image = NULL;
+static off_t Nimage = 0;
+static off_t nextImage = 0;
+
+// update astrometry of all chips relative to the average positions
+// if NTHREADS is non-zero, call the threaded version of this function
 int UpdateChips (Catalog *catalog, int Ncatalog) {
 
-  int Nskip, Nmosaic, NnewFit, NoldFit;
+  off_t Nskip, Nmosaic, NnewFit, NoldFit;
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  off_t i, Nimage, Nraw, Nref, nFitAstr;
-  Image *image;
+  off_t i, Nraw, Nref, nFitAstr;
   StarData *raw, *ref;
   Coords *oldCoords;
@@ -16 +41 @@
   char *mode;
 
+  if (NTHREADS) {
+    UpdateChips_threaded (catalog, Ncatalog);
+    return TRUE;
+  }
+
   Nskip = Nmosaic = NnewFit = NoldFit = 0;
 
@@ -25 +55 @@
   ALLOCATE (mode, char, Nimage);
 
+  // XXX for faster processing in the future, this can be easily run in parallel
+  // each chip is fitted independently, so we could do N at once in parallel
   for (i = 0; i < Nimage; i++) {
+
+    // XXX looks like everything below is thread safe : we can unroll this into a set of
+    // helper functions that grab the next available chip....
 
     /* skip all except WRP images */
@@ -47 +82 @@
       Nskip ++;
       mode[i] = 0;
+      free (raw);
       continue;
     }
@@ -62 +98 @@
     // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
     if (!FitChip (raw, ref, Nraw, &image[i])) {
-      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
 
       // restore status quo ante
@@ -105 +141 @@
 
   plotChipFits (Ro, Do, mode, Nimage);
-
-  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic);
+  free (Ro);
+  free (Do);
+  free (mode);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+           (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic);
   return (TRUE);
 }
 
+int UpdateChips_threaded (Catalog *catalog, int Ncatalog) {
+
+  int i;
+  off_t Nskip, Nmosaic, NnewFit, NoldFit;
+  double *Ro, *Do;
+  char *mode;
+
+  image = getimages (&Nimage, NULL);
+  nextImage = 0;
+
+  // save fit results for summary plot
+  ALLOCATE (Ro, double, Nimage);
+  ALLOCATE (Do, double, Nimage);
+  ALLOCATE (mode, char, Nimage);
+
+  pthread_attr_t attr;
+  pthread_attr_init (&attr);
+  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
+  
+  pthread_t *threads;
+  ALLOCATE (threads, pthread_t, NTHREADS);
+
+  ThreadInfo *threadinfo;
+  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
+
+  // launch N worker threads
+  for (i = 0; i < NTHREADS; i++) {
+    threadinfo[i].entry = i;
+    threadinfo[i].state = THREAD_RUN;
+    threadinfo[i].Ro = Ro;
+    threadinfo[i].Do = Do;
+    threadinfo[i].mode = mode;
+    threadinfo[i].Nskip = 0;
+    threadinfo[i].Nmosaic = 0;
+    threadinfo[i].NnewFit = 0;
+    threadinfo[i].NoldFit = 0;
+    threadinfo[i].catalog  =  catalog;
+    threadinfo[i].Ncatalog = Ncatalog;
+    pthread_create (&threads[i], NULL, UpdateChips_worker, &threadinfo[i]);
+  }
+  pthread_attr_destroy (&attr);
+
+  // wait until all threads have finished
+  while (1) {
+    int allDone = TRUE;
+    for (i = 0; i < NTHREADS; i++) {
+      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
+    }
+    if (allDone) {
+      break;
+    }
+    usleep (500000);
+  }
+
+  // all threads are done, free the threads array and grab the info
+  free (threads);
+  
+  plotChipFits (Ro, Do, mode, Nimage);
+  free (Ro);
+  free (Do);
+  free (mode);
+
+  // report stats & summary from the threads
+  Nskip = Nmosaic = NnewFit = NoldFit = 0;
+  for (i = 0; i < NTHREADS; i++) {
+    fprintf (stderr, "UpdateChips thread %d : %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+             i, (int) threadinfo[i].NnewFit, (int) threadinfo[i].NoldFit, (int) threadinfo[i].Nskip, (int) threadinfo[i].Nmosaic);
+
+    NnewFit += threadinfo[i].NnewFit;
+    NoldFit += threadinfo[i].NoldFit;
+    Nskip   += threadinfo[i].Nskip;
+    Nmosaic += threadinfo[i].Nmosaic;
+  }
+  free (threadinfo);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+           (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic);
+  return (TRUE);
+}
+
+void *UpdateChips_worker (void *data) {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+  off_t Nraw, Nref, nFitAstr;
+  StarData *raw, *ref;
+  Coords *oldCoords;
+  float dXpixSys, dYpixSys;
+
+  ThreadInfo *threadinfo = data;
+
+  while (1) {
+
+    off_t i = getNextImageForThread();
+    if (i == -1) {
+      threadinfo->state = THREAD_DONE;
+      return NULL;
+    }
+
+    /* skip all except WRP images */
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
+      threadinfo->Nmosaic ++;
+      threadinfo->mode[i] = 0;
+      continue;
+    }
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (threadinfo->catalog, threadinfo->Ncatalog, i, &Nraw, MODE_MOSAIC);
+    if (!raw) {
+      threadinfo->Nskip ++;
+      threadinfo->mode[i] = 0;
+      continue;
+    }
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (threadinfo->catalog, threadinfo->Ncatalog, i, &Nref, MODE_MOSAIC);
+    if (!ref) {
+      threadinfo->Nskip ++;
+      threadinfo->mode[i] = 0;
+      free (raw);
+      continue;
+    }
+
+    // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
+    assert (Nraw == Nref);
+
+    // save these in case of failure
+    saveCoords (&image[i].coords, i);
+    dXpixSys = image[i].dXpixSys;
+    dYpixSys = image[i].dYpixSys;
+    nFitAstr = image[i].nFitAstrom;
+
+    // FitChip does iterative, clipped fitting
+    // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
+    if (!FitChip (raw, ref, Nraw, &image[i])) {
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
+
+      // restore status quo ante
+      oldCoords = getCoords (i);
+      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+      threadinfo->mode[i] = 1;
+      threadinfo->NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    }
+
+    if (!checkStarMap (i)) {
+      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      // restore status quo ante
+      oldCoords = getCoords (i);
+      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+      threadinfo->mode[i] = 2;
+      image[i].flags |= ID_IMAGE_ASTROM_POOR;
+      threadinfo->NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    } 
+
+    saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+    threadinfo->mode[i] = 3;
+    threadinfo->NnewFit ++;
+    free (raw);
+    free (ref);
+  }
+  
+  // we should never reach here...
+  return NULL;
+}
+
+// mutex to lock UpdateChips_worker operations 
+static pthread_mutex_t UpdateChips_mutex = PTHREAD_MUTEX_INITIALIZER;
+
+// we have an array of chips (image, Nimage).  we need to hand out images one at a time to
+// the worker threads as they need
+off_t getNextImageForThread () {
+
+  pthread_mutex_lock (&UpdateChips_mutex);
+  if (nextImage >= Nimage) {
+    pthread_mutex_unlock (&UpdateChips_mutex);
+    return (-1);
+  }
+  off_t thisImage = nextImage;
+  nextImage ++;
+
+  pthread_mutex_unlock (&UpdateChips_mutex);
+  return (thisImage);
+}
+
+// This function uses the mosaic array in MosaicOps.c to associate PHU + CHIP this is
+// thread safe : the lookup (getMosaicForImage) is selecting an element of a previously
+// allocated and assigned static array.  This function does NOT use the thread-unsafe (and
+// somewhat slower) functions in libdvo/src/mosaic_astrom.c to associate CHIP and PHU
 int saveCenter (Image *image, double *Ro, double *Do, int im) {
 
@@ -136 +379 @@
     // note that for a Simple image, L,M = P,Q
     XY_to_LM (&L, &M, X, Y, imcoords);
-    LM_to_RD (&R, &D, L, M, imcoords);
+    LM_to_RD (&R, &D, L, M, imcoords); // because of the block above, ctype is not -WRP 
   } else {
     XY_to_LM (&L, &M, X, Y, imcoords);
     XY_to_LM (&P, &Q, L, M, moscoords);
-    LM_to_RD (&R, &D, P, Q, moscoords);
-  }
-
-  double Rmid;
-  if (UserCatalog) {
-      Rmid = UserCatalogRA;
-  } else {
-      Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax);
-  }
+    LM_to_RD (&R, &D, P, Q, moscoords); // moscoords.ctype is -DIS
+  }
+
+  double Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax);
 
   R = ohana_normalize_angle_to_midpoint (R, Rmid);