IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 30, 2012, 1:46:12 PM (14 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20120405

Location:
trunk
Files:
22 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/Ohana/src/relphot/Makefile

    r33651 r33963  
    4646$(SRC)/relphot_objects.$(ARCH).o         \
    4747$(SRC)/select_images.$(ARCH).o   \
     48$(SRC)/client_logger.$(ARCH).o   \
    4849$(SRC)/setExclusions.$(ARCH).o   \
    4950$(SRC)/setMrelFinal.$(ARCH).o    \
    … …  
    7374$(SRC)/relphot_objects.$(ARCH).o         \
    7475$(SRC)/relphot_client.$(ARCH).o  \
     76$(SRC)/client_logger.$(ARCH).o   \
    7577$(SRC)/setExclusions.$(ARCH).o   \
    7678$(SRC)/setMrelFinal.$(ARCH).o    \

  • trunk/Ohana/src/relphot/include/relphot.h

    r33651 r33963  
    33# include <kapa.h>
    44# include <signal.h>
     5# include <pthread.h>
    56
    67/* # define GRID_V1 */
    … …  
    2627  float dMcal;
    2728  float dMsys;
    28   short nFitPhotom;
     29  unsigned short nFitPhotom;
    2930  short Xm;
    3031  float secz;
    3132  float ubercalDist;
    3233  unsigned int flags;
     34  char skipCal;               // if TRUE, this mosaic is incomplete and should not be calibrated
    3335  Coords coords;
    3436} Mosaic;
     37
     38typedef enum {
     39  STATS_NONE,
     40  STATS_MEAN,
     41  STATS_MEDIAN,
     42  STATS_WT_MEAN,
     43  STATS_INNER_MEAN,
     44  STATS_INNER_WTMEAN,
     45  STATS_CHI_INNER_MEAN,
     46  STATS_CHI_INNER_WTMEAN
     47} ListStatsMode;
    3548
    3649typedef struct {
    … …  
    4659  double total;
    4760  int    Nmeas;
     61  ListStatsMode statmode;
    4862} StatType;
    4963
    … …  
    7993
    8094/* global variables set in parameter file */
    81 # define MAX_PATH_LENGTH 1024
    82 char   ImageCat[MAX_PATH_LENGTH];
    83 char   ImageTemplate[MAX_PATH_LENGTH];
    84 char   CatTemplate[MAX_PATH_LENGTH];
    85 char   GSCFILE[MAX_PATH_LENGTH];
     95# define DVO_MAX_PATH 1024
     96char   ImageCat[DVO_MAX_PATH];
     97char   ImageTemplate[DVO_MAX_PATH];
     98char   CatTemplate[DVO_MAX_PATH];
     99char   GSCFILE[DVO_MAX_PATH];
    86100char  *CATDIR;
    87101char   CATMODE[16];    /* raw, mef, split, mysql */
    88102char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
    89 char   CameraConfig[MAX_PATH_LENGTH];
     103char   CameraConfig[DVO_MAX_PATH];
    90104char   CAMERA[64];    /* eg, gpc1 */
    91 char   SKY_TABLE[MAX_PATH_LENGTH];
     105char   SKY_TABLE[DVO_MAX_PATH];
    92106int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
    93107
    … …  
    111125int    PARALLEL_MANUAL;
    112126int    PARALLEL_SERIAL;
     127
     128int    NTHREADS;
    113129
    114130int    VERBOSE;
    … …  
    141157int    PLOTDELAY;
    142158int    UpdateAverages;
     159int    ApplyOffsets;
    143160
    144161char  *PhotcodeList;
    … …  
    194211int           args                PROTO((int argc, char **argv));
    195212int           args_client         PROTO((int argc, char **argv));
    196 int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
     213int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog, int Ncat));
    197214void          clean_images        PROTO((void));
    198215void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr));
    … …  
    236253void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
    237254void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
    238 void          initMosaics         PROTO((Image *image, off_t Nimage));
     255void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
    239256void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
    240257void          initialize          PROTO((int argc, char **argv));
    241258void          initialize_client   PROTO((int argc, char **argv));
    242 void          initstats           PROTO((char *mode));
    243 int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
     259void          liststats_setmode   PROTO((StatType *stats, char *strmode));
     260int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
    244261Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath));
    245262Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog));
    246263
    247264SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region));
    248 Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
     265Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage, SkyRegion *region));
    249266
    250267int           main                PROTO((int argc, char **argv));
    … …  
    273290int           setMmos             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
    274291int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    275 void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
     292void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage));
    276293int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr));
    277294int           setMave             PROTO((Catalog *catalog, int Ncatalog));
    … …  
    319336int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage);
    320337ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage);
     338
     339int client_logger_init ();
     340int client_logger_message (char *format,...);
     341

  • trunk/Ohana/src/relphot/src/ConfigInit.c

    r33651 r33963  
    3535  // force CATDIR to be absolute (so parallel mode will work)
    3636  char *tmpcatdir = NULL;
    37   ALLOCATE (tmpcatdir, char, MAX_PATH_LENGTH);
     37  ALLOCATE (tmpcatdir, char, DVO_MAX_PATH);
    3838  GetConfig (config, "CATDIR",                 "%s",  0, tmpcatdir);
    39   CATDIR = abspath (tmpcatdir, MAX_PATH_LENGTH);
     39  CATDIR = abspath (tmpcatdir, DVO_MAX_PATH);
    4040  free (tmpcatdir);
    4141

  • trunk/Ohana/src/relphot/src/GridOps.c

    r33651 r33963  
    527527  double *list, *dlist;
    528528  float Msys, Mrel, Mcal, Mmos;
     529
    529530  StatType stats;
     531  liststats_setmode (&stats, STATMODE);
    530532 
    531533  if (!USE_GRID) return;
    … …  
    595597    }
    596598
    597     liststats (list, dlist, N, &stats);
     599    liststats (list, dlist, NULL, N, &stats);
    598600    gridM[i] = stats.mean;
    599601    gridS[i] = stats.sigma;

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r33651 r33963  
    4141
    4242// relationships between the measure,catalog set and the images:
    43 static IDX_T       **MeasureToImage; // image index from measure,catalog   : MeasureToImage[cat][meas] = ImageIndex
    44 static IDX_T       **ImageToCatalog; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
    45 static IDX_T       **ImageToMeasure; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
     43static IDX_T       **MeasureToImage = NULL; // image index from measure,catalog   : MeasureToImage[cat][meas] = ImageIndex
     44static IDX_T       **ImageToCatalog = NULL; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
     45static IDX_T       **ImageToMeasure = NULL; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
    4646
    4747// MeasureToImage was 'bin'
    … …  
    322322  off_t i;
    323323
     324  if (!MeasureToImage) return -1;
     325
    324326  i = MeasureToImage[cat][meas];
    325327  return (i);
    … …  
    355357  short distance;
    356358
     359  if (!MeasureToImage) return -1;
     360
    357361  i = MeasureToImage[cat][meas];
    358362  if (i == -1) return (1000);
    … …  
    395399  float Msys, Mrel, Mmos, Mgrid, Mflat;
    396400  double *list, *dlist, *Mlist, *dMlist;
     401
    397402  StatType stats;
     403  liststats_setmode (&stats, STATMODE);
    398404
    399405  if (FREEZE_IMAGES) return;
    … …  
    514520    }
    515521
    516     liststats (list, dlist, N, &stats);
     522    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
     523    liststats (list, dlist, NULL, N, &stats);
    517524    image[i].Mcal  = stats.mean;
    518525    image[i].dMcal = stats.error;
    … …  
    523530
    524531    // bright end scatter
    525     liststats (Mlist, dMlist, Nbright, &stats);
     532    liststats (Mlist, dMlist, NULL, Nbright, &stats);
    526533    image[i].dMagSys = stats.sigma;
    527534
    … …  
    556563  double *mlist, *slist, *dlist;
    557564  double MaxOffset, MaxScatter, MedOffset;
    558   StatType stats;
    559565
    560566  if (FREEZE_IMAGES) return;
    … …  
    574580    N++;
    575581  }
    576   initstats ("MEAN");
    577   liststats (mlist, dlist, N, &stats);
     582
     583  // use a straight mean to find the global image statistics (no weighting)
     584  StatType stats;
     585  liststats_setmode (&stats, "MEAN");
     586
     587  liststats (mlist, dlist, NULL, N, &stats);
    578588  MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
    579589  MedOffset = stats.median;
    580   liststats (slist, dlist, N, &stats);
     590
     591  liststats (slist, dlist, NULL, N, &stats);
    581592  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
    582593  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
    … …  
    599610
    600611  fprintf (stderr, "%d images marked poor\n", Nmark);
    601   initstats (STATMODE);
    602612  free (mlist);
    603613  free (slist);
    … …  
    720730  double *list, *dlist;
    721731  float Mcal, Mmos, Mgrid;
     732
    722733  StatType stats;
    723 
    724734  bzero (&stats, sizeof (StatType));
    725735  if (FREEZE_IMAGES) return (stats);
    … …  
    751761  }
    752762
    753   liststats (list, dlist, n, &stats);
     763  liststats_setmode (&stats, "MEAN");
     764
     765  liststats (list, dlist, NULL, n, &stats);
    754766  free (list);
    755767  free (dlist);
    … …  
    779791  }
    780792
    781   liststats (list, dlist, n, &stats);
     793  liststats_setmode (&stats, "MEAN");
     794
     795  liststats (list, dlist, NULL, n, &stats);
    782796  free (list);
    783797  free (dlist);
    … …  
    807821  }
    808822
    809   liststats (list, dlist, n, &stats);
     823  liststats_setmode (&stats, "MEAN");
     824
     825  liststats (list, dlist, NULL, n, &stats);
    810826  free (list);
    811827  free (dlist);
    … …  
    835851  }
    836852
    837   liststats (list, dlist, n, &stats);
     853  liststats_setmode (&stats, "MEAN");
     854
     855  liststats (list, dlist, NULL, n, &stats);
    838856  free (list);
    839857  free (dlist);

  • trunk/Ohana/src/relphot/src/ImageSubset.c

    r33651 r33963  
    9292}
    9393
     94// STATUS is value expected for success
     95# define CHECK_STATUS(STATUS,MSG,...)                                   \
     96  if (!(STATUS)) {                                                      \
     97    fprintf (stderr, MSG, __VA_ARGS__);                                 \
     98    return FALSE;                                                       \
     99  }
     100
    94101int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage) {
    95102
    … …  
    167174  }
    168175
    169   gfits_fwrite_header  (f, &header);
    170   gfits_fwrite_matrix  (f, &matrix);
    171   gfits_fwrite_Theader (f, &theader);
    172   gfits_fwrite_table  (f, &ftable);
     176  int status;
     177  status = gfits_fwrite_header  (f, &header);
     178  CHECK_STATUS (status, "ERROR: cannot write header for image subset %s\n", filename);
     179
     180  status = gfits_fwrite_matrix  (f, &matrix);
     181  CHECK_STATUS (status, "ERROR: cannot write matrix for image subset %s\n", filename);
     182
     183  status = gfits_fwrite_Theader (f, &theader);
     184  CHECK_STATUS (status, "ERROR: cannot write table header for image subset %s\n", filename);
     185
     186  status = gfits_fwrite_table  (f, &ftable);
     187  CHECK_STATUS (status, "ERROR: cannot write table data for image subset %s\n", filename);
    173188
    174189  gfits_free_header (&header);
    … …  
    177192  gfits_free_table (&ftable);
    178193
    179   fclose (f);
    180   fflush (f);
     194  int fd = fileno (f);
     195
     196  status = fflush (f);
     197  CHECK_STATUS (!status, "ERROR: cannot flush file image subset %s\n", filename);
     198
     199  status = fsync (fd);
     200  CHECK_STATUS (!status, "ERROR: cannot flush file image subset %s\n", filename);
     201
     202  status = fclose (f);
     203  CHECK_STATUS (!status, "ERROR: problem closing image subset file %s\n", filename);
    181204
    182205  return TRUE;

  • trunk/Ohana/src/relphot/src/MosaicOps.c

    r33651 r33963  
    5050  ALLOCATE (mosaic, Mosaic, NMOSAIC);
    5151
    52   ALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
     52  ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
    5353  ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
    5454  ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
    … …  
    9595    if (found) continue;
    9696   
    97     /* a new mosaic, define ranges */
     97    /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
    9898    mosaic[Nmosaic].start = start;
    9999    mosaic[Nmosaic].stop  = stop;
    100     mosaic[Nmosaic].Mcal  = 0.0;
    101     mosaic[Nmosaic].dMcal = 0.0;
     100    mosaic[Nmosaic].Mcal  = 0.0; // note : mosaic stores only offsets relative to the original image values
     101    mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
    102102    mosaic[Nmosaic].dMsys = 0.0;
    103103    mosaic[Nmosaic].Xm    = 0.0;
    … …  
    152152}
    153153
    154 /* find mosaic frames (unique time periods & photcode name matches mosaic) */
    155 void initMosaics (Image *image, off_t Nimage) {
     154# define MARKTIME(MSG,...) { \
     155  float dtime; \
     156  gettimeofday (&stopTimer, (void *) NULL); \
     157  dtime = DTIME (stopTimer, startTimer); \
     158  fprintf (stderr, MSG, __VA_ARGS__); }
     159
     160/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...)
     161   this function will also identify the images NOT in the subset which belong to a selected mosaic
     162 */
     163void initMosaics (Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage) {
    156164
    157165  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
    … …  
    159167  char *pname;
    160168
     169  struct timeval startTimer, stopTimer;
     170
    161171  if (!MOSAIC_ZEROPT) return;
     172
     173  gettimeofday (&startTimer, (void *) NULL);
    162174
    163175  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no
    … …  
    166178   */
    167179
    168   // generate a list of unique start times (these define the mosaics)
    169   ALLOCATE (startTimes, unsigned int, Nimage);
    170   for (i = 0; i < Nimage; i++) {
    171     startTimes[i] = image[i].tzero;
    172   }
    173   sort_times (startTimes, Nimage);
     180  // generate a list of all subset image start times
     181  ALLOCATE (startTimes, unsigned int, Nsubset);
     182  for (i = 0; i < Nsubset; i++) {
     183    startTimes[i] = subset[i].tzero;
     184  }
     185  sort_times (startTimes, Nsubset);
     186  MARKTIME("create array of all image obstimes: %f sec\n", dtime);
    174187 
    175188  Nmosaic = 0;
    … …  
    178191  startTimesMosaic[0] = startTimes[0];
    179192
    180   for (i = 0; i < Nimage; i++) {
     193  // generate a list of the unique start times (these define the mosaics)
     194  for (i = 0; i < Nsubset; i++) {
    181195    if (startTimes[i] < startTimesMosaic[Nmosaic]) {
    182196      fprintf (stderr, "error?\n");
    … …  
    192206  }
    193207  Nmosaic ++;
     208  MARKTIME("create subset array of mosaic obstimes: %f sec\n", dtime);
    194209
    195210  // now I have a list of uniq start times, and they are in order
    … …  
    212227    mosaic[i].secz  = NAN;
    213228    mosaic[i].photcode = 0;
     229    mosaic[i].skipCal = FALSE;
     230   
     231    memset (&mosaic[i].coords, 0, sizeof(Coords));
    214232
    215233    MosaicN_IMAGE[i] = 10;
    … …  
    219237  }
    220238
    221   ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
    222 
    223   // assign each image to a mosaic
     239  int Nskip, Nmark;
     240  Nskip = Nmark = 0;
     241  // find any mosaics (startTimesMosaic) which match unselected images
    224242  for (i = 0; i < Nimage; i++) {
    225     ImageToMosaic[i] = -1;
    226 
     243    if (inSubset[i]) {
     244      Nskip ++;
     245      continue;
     246    }
     247   
    227248    /* select valid mosaic images by photcode */
    228249    pname = GetPhotcodeNamebyCode (image[i].photcode);
     250    if (!pname) continue;
    229251    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
    230252    if (status) continue;
    231253
    232254    /* set image time range */
    233     // start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
    234     // stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    235 
    236255    start = image[i].tzero;
    237256    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    238257
     258    /* find a matching mosaic */
     259    j = findMosaic(startTimesMosaic, Nmosaic, start);
     260    if (j != -1) {
     261      // mark this mosaic as bad
     262      mosaic[j].skipCal = TRUE;
     263      Nmark ++;
     264    }
     265  }
     266  fprintf (stderr, "%d total images, %d overlap skyregion, %d do not overlap, but match overlapping mosaics\n", (int) Nimage, (int) Nskip, (int) Nmark);
     267  MARKTIME("find unselected images matching selected mosaics: %f sec\n", dtime);
     268
     269  ALLOCATE (ImageToMosaic, off_t, Nsubset); // mosaic to which image belongs
     270
     271  // assign each image to a mosaic
     272  for (i = 0; i < Nsubset; i++) {
     273    ImageToMosaic[i] = -1;
     274
     275    /* select valid mosaic images by photcode */
     276    pname = GetPhotcodeNamebyCode (subset[i].photcode);
     277    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
     278    if (status) continue;
     279
     280    start = subset[i].tzero;
     281    stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);
     282
    239283    j = findMosaic(startTimesMosaic, Nmosaic, start);
    240284    if (j == -1) {
    241       fprintf (stderr, "error?\n");
     285      fprintf (stderr, "programming error? all subset images should belong to a mosaic\n");
    242286      abort();
    243287    }
    … …  
    268312    mosaic[j].dMsys = 0.0;
    269313    mosaic[j].Xm    = 0.0;
    270     mosaic[j].flags  = image[i].flags;
    271     mosaic[j].secz  = image[i].secz;
    272     mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
     314    mosaic[j].flags  = subset[i].flags;
     315    mosaic[j].secz  = subset[i].secz;
     316    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode);
    273317  }
    274318
    … …  
    278322  free (startTimesMosaic);
    279323
    280   initMosaicGrid (image, Nimage);
    281 
    282   fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
     324  initMosaicGrid (subset, Nsubset);
     325
     326  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic);
    283327  return;
    284328}
    285329
    286 // use bisection to find the overlapping mosaic
     330// use bisection to find the overlapping mosaic (returns exact match)
     331// startTimes is a sorted, unique list of times
     332// start might not be in the list
    287333off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start) {
    288334
    … …  
    290336
    291337  // find the last mosaic before start
    292   Nlo = 0; // startTimes[Nlo] guaranteed to be <= start
    293   Nhi = Nmosaic - 1; // startTimes[Nhi] guaranteed to be >= start
    294   while (Nhi - Nlo > 10) {
     338  Nlo = 0; // first valid startTimes value
     339  Nhi = Nmosaic - 1; // last valid startTimes value
     340
     341  // if start is not in this range, return -1
     342  if (start < startTimes[Nlo]) return (-1);
     343  if (start > startTimes[Nhi]) return (-1);
     344
     345  while (Nhi - Nlo > 4) {
    295346    N = 0.5*(Nlo + Nhi);
    296347    if (startTimes[N] < start) {
    … …  
    300351    }
    301352  }
    302 
    303   // check for the matched mosaic starting from Nlo
    304   for (N = Nlo; N < Nmosaic; N++) {
    305     if (start > startTimes[N])  continue;
    306     return (N);
     353  // we now have : startTimes[Nlo] < start <= starTimes[Nhi]
     354
     355  // find a matched mosaic starting from Nlo, or return -1
     356  for (N = Nlo; N <= Nhi; N++) {
     357    if (startTimes[N] == start) return N;
    307358  }
    308359  return (-1);
    … …  
    394445}
    395446
     447// XXX : what about mosaics with skipCal == TRUE?
    396448void setMcalFinal () {
    397449
    … …  
    569621}
    570622
     623typedef struct {
     624  int Nfew;
     625  int Nbad;
     626  int Ncal;
     627  int Ngrid;
     628  int Nrel;
     629  int Nsys;
     630  int Nskip;
     631  off_t Nmax;
     632  int PoorImages;
     633  double *list;
     634  double *dlist;
     635  double *Mlist;
     636  double *dMlist;
     637} SetMmosInfo;
     638
     639enum {THREAD_RUN, THREAD_DONE};
     640
     641typedef struct {
     642  int entry;
     643  int state;
     644  Catalog *catalog;
     645  Image *image;
     646  FlatCorrectionTable *flatcorr;
     647  SetMmosInfo info;
     648} ThreadInfo;
     649
     650int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr);
     651void *setMmos_worker (void *data);
     652int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);
     653
     654void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) {
     655  info->Nfew = 0;
     656  info->Nbad = 0;
     657  info->Ncal = 0;
     658  info->Nrel = 0;
     659  info->Ngrid = 0;
     660  info->Nskip = 0;
     661  info->Nsys = 0;
     662
     663  info->Nmax = Nmax;
     664  info->PoorImages = PoorImages;
     665
     666  if (allocLists) {
     667    ALLOCATE (info->list, double, Nmax);
     668    ALLOCATE (info->dlist, double, Nmax);
     669    ALLOCATE (info->Mlist, double, Nmax);
     670    ALLOCATE (info->dMlist, double, Nmax);
     671  }
     672}
     673
     674void SetMmosInfoFree (SetMmosInfo *info) {
     675  free (info->list);
     676  free (info->dlist);
     677  free (info->Mlist);
     678  free (info->dMlist);
     679}
     680
     681void SetMmosInfoAccum (SetMmosInfo *summary, SetMmosInfo *results) {
     682  summary->Nfew  += results->Nfew ;
     683  summary->Nsys  += results->Nsys ;
     684  summary->Nbad  += results->Nbad ;
     685  summary->Ncal  += results->Ncal ;
     686  summary->Nrel  += results->Nrel ;
     687  summary->Ngrid += results->Ngrid;
     688  summary->Nskip += results->Nskip;
     689}
     690
     691static int npass_output = 0;
     692
     693// mutex to lock setMmos_worker operations
     694static pthread_mutex_t setMmos_mutex = PTHREAD_MUTEX_INITIALIZER;
     695static int nextMosaic = 0;
     696
     697// we have an array of mosaics (mosaic, Nmosaic).  we need to hand out mosaics one at a time to
     698// the worker threads as they need
     699off_t getNextMosaicForThread () {
     700
     701  pthread_mutex_lock (&setMmos_mutex);
     702  if (nextMosaic >= Nmosaic) {
     703    pthread_mutex_unlock (&setMmos_mutex);
     704    return (-1);
     705  }
     706  int thisMosaic = nextMosaic;
     707  nextMosaic ++;
     708
     709  pthread_mutex_unlock (&setMmos_mutex);
     710  return (thisMosaic);
     711}
     712
    571713int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    572714
    573   off_t i, j, m, c, n, N, Nmax;
    574   int mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nbright;
    575   float Msys, Mrel, Mcal, Mgrid, Mflat;
    576   double *list, *dlist, *Mlist, *dMlist;
    577   StatType stats;
     715  off_t i, N, Nmax;
    578716  Image *image;
    579717
    … …  
    581719  if (FREEZE_MOSAICS) return (FALSE);
    582720
     721  if (NTHREADS) {
     722    int status = setMmos_threaded (catalog, PoorImages, flatcorr);
     723    return status;
     724  }
     725
    583726  image = getimages (&N, NULL);
    584727
    585728  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
    586 
    587   int Nsecfilt = GetPhotcodeNsecfilt ();
    588729
    589730  if (PoorImages) {
    … …  
    597738    Nmax = MAX (Nmax, N_onMosaic[i]);
    598739  }
    599   ALLOCATE (list, double, Nmax);
    600   ALLOCATE (dlist, double, Nmax);
    601   ALLOCATE (Mlist, double, Nmax);
    602   ALLOCATE (dMlist, double, Nmax);
    603 
    604   Nfew = Nbad = Ncal = Nrel = Ngrid = Nsys = 0;
    605 
    606   for (i = 0; i < Nmosaic; i++) {
    607    
    608     /* on PoorImages run, skip good images */
    609     if (PoorImages) {
    610       bad = mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
    611       if (!bad) continue;
    612     }     
    613 
    614     // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
    615     // we optionally do not recalibrate images with UBERCAL zero points
    616     if (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL) {
    617       mosaic[i].ubercalDist = 0;
    618       // propagate ubercalDist to the images
    619       for (j = 0; j < MosaicN_Image[i]; j++) {
    620         off_t im = MosaicToImage[i][j];
    621         image[im].ubercalDist = mosaic[i].ubercalDist;
    622         // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
    623       }
    624       if (KEEP_UBERCAL) continue;
    625     }
    626 
    627     int minUbercalDist = 1000;
    628 
    629     // number of stars to measure the bright-end scatter
    630     Nbright = 0;
    631 
    632     N = 0;
    633     for (j = 0; j < N_onMosaic[i]; j++) {
    634      
    635       m = MosaicToMeasure[i][j];
    636       c = MosaicToCatalog[i][j];
    637      
    638       if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    639           Nbad ++;
    640           continue;
    641       }
    642       Mcal  = getMcal  (m, c, flatcorr, catalog);
    643       if (isnan(Mcal)) {
    644           Ncal++;
    645           continue;
    646       }
    647       Mgrid = getMgrid (m, c);
    648       if (isnan(Mgrid)) {
    649           Ngrid ++;
    650           continue;
    651       }
    652       Mrel  = getMrel  (catalog, m, c);
    653       if (isnan(Mrel)) {
    654           Nrel ++;
    655           continue;
    656       }
    657      
    658       // image.Mcal is not supposed to include the flat-field correction, so we need to
    659       // apply that offset as well here for this image (in other words, each detection is
    660       // being compared to the model, excluding the zero point, Mcal.  The model includes
    661       // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
    662 
    663       Mflat = getMflat (m, c, flatcorr, catalog);
    664 
    665       n = catalog[c].measureT[m].averef;
    666       Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
    667       if (isnan(Msys)) {
    668         Nsys++;
    669         continue;
    670       }
    671 
    672       PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
    673       if (!code) goto skip;
    674       if (code->equiv < 1) goto skip;
    675       int Nsec = GetPhotcodeNsec (code->equiv);
    676       if (Nsec == -1) goto skip;
    677       minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);
    678 
    679     skip:
    680       list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
    681       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    682       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
    683         Mlist[Nbright] = list[N];
    684         dMlist[Nbright] = dlist[N];
    685         Nbright ++;
    686       }
    687       N++;
    688     }
    689     /* N_onMosaic[i] is all measurements, N is good measurements */
    690 
    691     /* too few good measurements or too many bad measurements (skip in PoorImages run) */
    692     if (!PoorImages) {
    693       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[i]);
    694       if (mark) {
    695         if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[MosaicToImage[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  N_onMosaic[i]); }
    696         mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW;
    697         Nfew ++;
    698       } else {
    699         mosaic[i].flags &= ~ID_IMAGE_PHOTOM_FEW;
    700       }
    701     }
    702     liststats (list, dlist, N, &stats);
    703     if (VERBOSE2 && PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
    704 
    705     mosaic[i].Mcal  = stats.mean;
    706     mosaic[i].dMcal = stats.error;
    707     mosaic[i].nFitPhotom = N;
    708     mosaic[i].Xm    = 100.0*log10(stats.chisq);
    709 
    710     plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
    711 
    712     // bright end scatter
    713     liststats (Mlist, dMlist, Nbright, &stats);
    714     mosaic[i].dMsys = stats.sigma;
    715 
    716     if (mosaic[i].Mcal < -CLOUD_TOLERANCE) {
    717       mosaic[i].Mcal = 0.0;
    718     }
    719 
    720     // minUbercalDist calculated here is the min value for any star owned by this image
    721     // since this particular image is tied to that star, bump its distance by 1
    722     mosaic[i].ubercalDist = minUbercalDist + 1;
    723 
    724     // propagate ubercalDist to the images
    725     for (j = 0; j < MosaicN_Image[i]; j++) {
    726       off_t im = MosaicToImage[i][j];
    727       image[im].ubercalDist = mosaic[i].ubercalDist;
    728       // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
    729     }
    730   }
    731   free (list);
    732   free (dlist);
    733   free (Mlist);
    734   free (dMlist);
    735 
    736   fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Ngrid, Nrel, Nsys);
     740
     741  SetMmosInfo info;
     742  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages);
     743
     744  for (i = 0; i < Nmosaic; i++) {
     745    setMmos_mosaic (&mosaic[i], i, image, catalog, &info, flatcorr);
     746  }
     747  SetMmosInfoFree (&info);
     748
     749  npass_output ++;
     750
     751  fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys);
    737752
    738753  if (PoorImages) {
    … …  
    744759}
    745760 
     761// 'mosaic' is a pointer to the current mosaic of interest (Nmos)
     762int setMmos_mosaic (Mosaic *myMosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) {
     763
     764  Image *imageReal;
     765  off_t j, NimageReal;
     766
     767  StatType stats;
     768  liststats_setmode (&stats, STATMODE);
     769
     770  double *list   = info->list;
     771  double *dlist  = info->dlist;
     772  double *Mlist  = info->Mlist;
     773  double *dMlist = info->dMlist;
     774
     775  assert (Nmos >= 0);
     776  assert (Nmos < Nmosaic);
     777
     778  imageReal = getimages (&NimageReal, NULL);
     779
     780  /* on PoorImages run, skip good images */
     781  if (info->PoorImages) {
     782    int bad = myMosaic[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
     783    if (!bad) return TRUE;
     784  }     
     785
     786  // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
     787  // we optionally do not recalibrate images with UBERCAL zero points
     788  if (myMosaic[0].flags & ID_IMAGE_PHOTOM_UBERCAL) {
     789    myMosaic[0].ubercalDist = 0;
     790    // propagate ubercalDist to the images
     791    for (j = 0; j < MosaicN_Image[Nmos]; j++) {
     792      off_t im = MosaicToImage[Nmos][j];
     793      assert (im < NimageReal);
     794      assert (im >= 0);
     795      image[im].ubercalDist = myMosaic[0].ubercalDist;
     796      // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
     797    }
     798    if (KEEP_UBERCAL) return TRUE;
     799  }
     800
     801  int Nsecfilt = GetPhotcodeNsecfilt ();
     802
     803  int minUbercalDist = 1000;
     804
     805  // do not modify the calibration for mosaics with partial images loaded (skipCal TRUE)
     806  if (myMosaic[0].skipCal) {
     807    info->Nskip ++;
     808    return TRUE;
     809  }
     810
     811  int testImage = FALSE;
     812  // testImage |= (abs(myMosaic[0].start - 1324104046) < 10);
     813  // testImage |= (abs(myMosaic[0].start - 1324103823) < 10);
     814  // testImage |= (abs(myMosaic[0].start - 1323003245) < 10);
     815  // testImage |= (abs(myMosaic[0].start - 1323003069) < 10);
     816  // testImage |= (abs(myMosaic[0].start - 1323003125) < 10);
     817  // testImage |= (abs(myMosaic[0].start - 1323003300) < 10);
     818  // testImage |= (abs(myMosaic[0].start - 1323003365) < 10);
     819  // testImage |= (abs(myMosaic[0].start - 1323003191) < 10);
     820  // testImage |= (abs(myMosaic[0].start - 1323003014) < 10);
     821  // testImage |= (abs(myMosaic[0].start - 1323003484) < 10);
     822  // testImage |= (abs(myMosaic[0].start - 1323003419) < 10);
     823  // testImage |= (abs(myMosaic[0].start - 1323002949) < 10);
     824
     825  FILE *fout = NULL;
     826  if (testImage) {
     827    char filename[64];
     828    snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output);
     829    fout = fopen (filename, "w");
     830  }
     831
     832  // number of stars to measure the bright-end scatter
     833  int Nbright = 0;
     834
     835  int N = 0;
     836  for (j = 0; j < N_onMosaic[Nmos]; j++) {
     837    float Msys, Mrel, Mcal, Mgrid, Mflat;
     838     
     839    off_t m = MosaicToMeasure[Nmos][j];
     840    off_t c = MosaicToCatalog[Nmos][j];
     841     
     842    if (fout) {
     843      Mcal  = getMcal  (m, c, flatcorr, catalog);
     844      Mgrid = getMgrid (m, c);
     845      Mrel  = getMrel  (catalog, m, c);
     846      Mflat = getMflat (m, c, flatcorr, catalog);
     847
     848      off_t n = catalog[c].measureT[m].averef;
     849      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
     850
     851      float delta = Msys - Mrel - Mcal - Mgrid + Mflat;
     852
     853      int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
     854
     855      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mcal, Mgrid, Mflat, delta, isBad);
     856    }
     857
     858    if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
     859      info->Nbad ++;
     860      continue;
     861    }
     862    Mcal  = getMcal  (m, c, flatcorr, catalog);
     863    if (isnan(Mcal)) {
     864      info->Ncal++;
     865      continue;
     866    }
     867    Mgrid = getMgrid (m, c);
     868    if (isnan(Mgrid)) {
     869      info->Ngrid ++;
     870      continue;
     871    }
     872    Mrel  = getMrel  (catalog, m, c);
     873    if (isnan(Mrel)) {
     874      info->Nrel ++;
     875      continue;
     876    }
     877     
     878    // image.Mcal is not supposed to include the flat-field correction, so we need to
     879    // apply that offset as well here for this image (in other words, each detection is
     880    // being compared to the model, excluding the zero point, Mcal.  The model includes
     881    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
     882
     883    Mflat = getMflat (m, c, flatcorr, catalog);
     884
     885    off_t n = catalog[c].measureT[m].averef;
     886    Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
     887    if (isnan(Msys)) {
     888      info->Nsys++;
     889      continue;
     890    }
     891
     892    PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
     893    if (!code) goto skip;
     894    if (code->equiv < 1) goto skip;
     895    int Nsec = GetPhotcodeNsec (code->equiv);
     896    if (Nsec == -1) goto skip;
     897    minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);
     898
     899  skip:
     900    assert (N < info->Nmax);
     901    assert (N >= 0);
     902    assert (Nbright < info->Nmax);
     903    assert (Nbright >= 0);
     904
     905    list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
     906    dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
     907    if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
     908      Mlist[Nbright] = list[N];
     909      dMlist[Nbright] = dlist[N];
     910      Nbright ++;
     911    }
     912    N++;
     913  }
     914  /* N_onMosaic[Nmos] is all measurements, N is good measurements */
     915
     916  if (fout) {
     917    fclose (fout);
     918  }
     919
     920  /* too few good measurements or too many bad measurements (skip in PoorImages run) */
     921  if (!info->PoorImages) {
     922    int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]);
     923    if (mark) {
     924      if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[MosaicToImage[Nmos][0]].name,  Nmos,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  N_onMosaic[Nmos]); }
     925      myMosaic[0].flags |= ID_IMAGE_PHOTOM_FEW;
     926      info->Nfew ++;
     927      if (testImage) {
     928        fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]));
     929      }
     930    } else {
     931      myMosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW;
     932    }
     933  }
     934
     935  liststats (list, dlist, NULL, N, &stats);
     936  if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
     937
     938  myMosaic[0].Mcal  = stats.mean;
     939  myMosaic[0].dMcal = stats.error;
     940  myMosaic[0].nFitPhotom = N;
     941  myMosaic[0].Xm    = 100.0*log10(stats.chisq);
     942
     943  if (testImage) {
     944    fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic[0].start, stats.mean, stats.error, myMosaic[0].nFitPhotom);
     945  }
     946
     947  plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
     948
     949  // bright end scatter
     950  liststats (Mlist, dMlist, NULL, Nbright, &stats);
     951  myMosaic[0].dMsys = stats.sigma;
     952
     953  if (myMosaic[0].Mcal < -CLOUD_TOLERANCE) {
     954    myMosaic[0].Mcal = 0.0;
     955  }
     956
     957  if (testImage) {
     958    fprintf (stderr, "%f %f  :  %d %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Xm, pow(10.0, 0.01*myMosaic[0].Xm));
     959  }
     960
     961  // minUbercalDist calculated here is the min value for any star owned by this image
     962  // since this particular image is tied to that star, bump its distance by 1
     963  myMosaic[0].ubercalDist = minUbercalDist + 1;
     964
     965  // propagate ubercalDist to the images
     966  for (j = 0; j < MosaicN_Image[Nmos]; j++) {
     967    off_t im = MosaicToImage[Nmos][j];
     968    image[im].ubercalDist = myMosaic[0].ubercalDist;
     969    // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
     970  }
     971  return TRUE;
     972}
     973 
     974int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
     975
     976  int i;
     977  off_t N;
     978
     979  Image *image = getimages (&N, NULL);
     980
     981  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
     982
     983  if (PoorImages) {
     984    // XXX use bad stars and measurements for PoorImages? or not?
     985    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
     986    IMAGE_BAD = 0;
     987  }
     988
     989  off_t Nmax = 0;
     990  for (i = 0; i < Nmosaic; i++) {
     991    Nmax = MAX (Nmax, N_onMosaic[i]);
     992  }
     993
     994  SetMmosInfo summary;
     995  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages);
     996
     997  pthread_attr_t attr;
     998  pthread_attr_init (&attr);
     999  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
     1000 
     1001  pthread_t *threads;
     1002  ALLOCATE (threads, pthread_t, NTHREADS);
     1003
     1004  ThreadInfo *threadinfo;
     1005  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
     1006
     1007  // each time this function is called, we cycle through the available mosaics
     1008  // make sure we start at 0
     1009  nextMosaic = 0;;
     1010
     1011  // launch N worker threads
     1012  for (i = 0; i < NTHREADS; i++) {
     1013    threadinfo[i].entry = i;
     1014    threadinfo[i].state = THREAD_RUN;
     1015    threadinfo[i].catalog  =  catalog;
     1016    threadinfo[i].image    =    image;
     1017    threadinfo[i].flatcorr = flatcorr;
     1018    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
     1019    pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]);
     1020  }
     1021  pthread_attr_destroy (&attr);
     1022
     1023  // wait until all threads have finished
     1024  while (1) {
     1025    int allDone = TRUE;
     1026    for (i = 0; i < NTHREADS; i++) {
     1027      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
     1028    }
     1029    if (allDone) {
     1030      break;
     1031    }
     1032    usleep (500000);
     1033  }
     1034
     1035  // all threads are done, free the threads array and grab the info
     1036  free (threads);
     1037 
     1038  // report stats & summary from the threads
     1039  for (i = 0; i < NTHREADS; i++) {
     1040    fprintf (stderr, "setMmos thread %d : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n",
     1041             i,
     1042             threadinfo[i].info.Nfew,
     1043             threadinfo[i].info.Nbad,
     1044             threadinfo[i].info.Ncal,
     1045             threadinfo[i].info.Ngrid,
     1046             threadinfo[i].info.Nrel,
     1047             threadinfo[i].info.Nsys,
     1048             threadinfo[i].info.Nskip);
     1049    SetMmosInfoAccum (&summary, &threadinfo[i].info);
     1050  }
     1051  fprintf (stderr, "total : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n",
     1052           summary.Nfew,
     1053           summary.Nbad,
     1054           summary.Ncal,
     1055           summary.Ngrid,
     1056           summary.Nrel,
     1057           summary.Nsys,
     1058           summary.Nskip);
     1059  free (threadinfo);
     1060
     1061  npass_output ++;
     1062
     1063  if (PoorImages) {
     1064    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
     1065    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
     1066    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
     1067  }
     1068  return TRUE;
     1069}
     1070
     1071void *setMmos_worker (void *data) {
     1072
     1073  ThreadInfo *threadinfo = data;
     1074
     1075  SetMmosInfo results;
     1076  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here
     1077
     1078  while (1) {
     1079
     1080    off_t i = getNextMosaicForThread();
     1081    if (i == -1) {
     1082      threadinfo->state = THREAD_DONE;
     1083      return NULL;
     1084    }
     1085
     1086    Catalog *catalog = threadinfo->catalog;
     1087    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
     1088    Image *image = threadinfo->image;
     1089
     1090    setMmos_mosaic (&mosaic[i], i, image, catalog, &results, flatcorr);
     1091    SetMmosInfoAccum (&threadinfo->info, &results);
     1092  }
     1093
     1094  SetMmosInfoFree (&results);
     1095  return NULL;
     1096}
     1097
    7461098// When we rationalize the images/mosaics, we are driving the negative cloud images back
    7471099// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
    … …  
    9551307  for (i = 0; i < Nmosaic; i++) {
    9561308    if (mosaic[i].flags & IMAGE_BAD) continue;
     1309    if (mosaic[i].skipCal) continue;
    9571310    list[n] = mosaic[i].Mcal;
    9581311    dlist[n] = 1;
    … …  
    9601313  }
    9611314
    962   liststats (list, dlist, n, &stats);
     1315  liststats_setmode (&stats, "MEAN");
     1316
     1317  liststats (list, dlist, NULL, n, &stats);
    9631318  free (list);
    9641319  free (dlist);
    … …  
    9811336  n = 0;
    9821337  for (i = 0; i < Nmosaic; i++) {
    983 
    9841338    if (mosaic[i].flags & IMAGE_BAD) continue;
     1339    if (mosaic[i].skipCal) continue;
    9851340    list[n] = mosaic[i].dMcal;
    9861341    dlist[n] = 1;
    … …  
    9881343  }
    9891344
    990   liststats (list, dlist, n, &stats);
     1345  liststats_setmode (&stats, "MEAN");
     1346
     1347  liststats (list, dlist, NULL, n, &stats);
    9911348  free (list);
    9921349  free (dlist);
    … …  
    10111368  for (i = 0; i < Nmosaic; i++) {
    10121369    if (mosaic[i].flags & IMAGE_BAD)  continue;
     1370    if (mosaic[i].skipCal) continue;
    10131371
    10141372    N = 0;
    … …  
    10321390  // fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n",  Nmosaic,  n);
    10331391
    1034   liststats (list, dlist, n, &stats);
     1392  liststats_setmode (&stats, "MEAN");
     1393
     1394  liststats (list, dlist, NULL, n, &stats);
    10351395  free (list);
    10361396  free (dlist);
    … …  
    10531413  n = 0;
    10541414  for (i = 0; i < Nmosaic; i++) {
    1055 
    10561415    if (mosaic[i].flags & IMAGE_BAD) continue;
     1416    if (mosaic[i].skipCal) continue;
    10571417    list[n] = pow(10.0, 0.01*mosaic[i].Xm);
    10581418    dlist[n] = 1;
    … …  
    10601420  }
    10611421
    1062   liststats (list, dlist, n, &stats);
     1422  liststats_setmode (&stats, "MEAN");
     1423
     1424  liststats (list, dlist, NULL, n, &stats);
    10631425  free (list);
    10641426  free (dlist);
    … …  
    10721434  double *mlist, *slist, *dlist;
    10731435  double MaxOffset, MedOffset, MaxScatter;
    1074   StatType stats;
    10751436
    10761437  if (!MOSAIC_ZEROPT) return;
    … …  
    10851446  for (i = N = 0; i < Nmosaic; i++) {
    10861447    if (mosaic[i].flags & IMAGE_BAD) continue;
     1448    if (mosaic[i].skipCal) continue;
    10871449    mlist[N] = mosaic[i].Mcal;
    10881450    slist[N] = mosaic[i].dMcal;
    … …  
    10901452    N++;
    10911453  }
    1092   initstats ("MEAN");
    1093   liststats (mlist, dlist, N, &stats);
     1454
     1455  StatType stats;
     1456  liststats_setmode (&stats, "MEAN");
     1457
     1458  liststats (mlist, dlist, NULL, N, &stats);
    10941459  MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
    10951460  MedOffset = stats.median;
    1096   liststats (slist, dlist, N, &stats);
     1461
     1462  liststats (slist, dlist, NULL, N, &stats);
    10971463  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
    10981464  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
    … …  
    11021468    // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
    11031469    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
     1470
     1471    if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue;
     1472    if (mosaic[i].skipCal) continue;
    11041473
    11051474    mark = FALSE;
    … …  
    11211490
    11221491  fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);
    1123   initstats (STATMODE);
     1492
    11241493  free (mlist);
    11251494  free (slist);

  • trunk/Ohana/src/relphot/src/StarOps.c

    r33651 r33963  
    22
    33static int Nmax;
    4 static double *list;
    5 static double *dlist;
    6 
    7 // When we rationalize the images/mosaics, we are driving the negative cloud images back
    8 // to 0.0.  At the same time, we make a guess to the effective impact on all other images,
    9 // driven by the coupling of common stars.  This array carries the impact of those offsets
    10 // on each star
    11 static double *Moffset;
    12 
     4
     5typedef struct {
     6  int Nfew;
     7  int Ncode;
     8  int Nsys;
     9  int Nbad;
     10  int Ncal;
     11  int Nmos;
     12  int Ngrid;
     13  double *list;
     14  double *dlist;
     15  double *wlist;
     16  double *aplist;
     17  double *daplist;
     18} SetMrelInfo;
     19
     20enum {THREAD_RUN, THREAD_DONE};
     21
     22typedef struct {
     23  int entry;
     24  int state;
     25  Catalog *catalog;
     26  int Ncatalog;
     27  FlatCorrectionTable *flatcorr;
     28  SetMrelInfo summary;
     29} ThreadInfo;
     30
     31// ** internal functions:
     32void *setMrel_worker (void *data);
     33int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr);
     34int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
     35int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
     36
     37// we want to allocate the stats list,dlist arrays only once (or once per thread).
     38// this function finds the largest array so we can allocate that max size when needed
    1339void initMrel (Catalog *catalog, int Ncatalog) {
    1440
    … …  
    2147    }
    2248  }
    23 
    24   ALLOCATE (list,    double, MAX (1, Nmax));
    25   ALLOCATE (dlist,   double, MAX (1, Nmax));
    26   ALLOCATE (Moffset, double, MAX (1, Nmax));
    2749} 
    2850
    … …  
    4769}
    4870
    49 int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    50 
    51   off_t j, k, m;
    52   int i, N, Nfew, Nsys, Nbad, Ncal, Nmos, Ngrid;
    53   float Msys, Mcal, Mmos, Mgrid;
    54   StatType stats;
    55 
    56   int Nsecfilt = GetPhotcodeNsecfilt ();
    57   Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
    58 
    59   for (i = 0; i < Ncatalog; i++) {
    60     for (j = 0; j < catalog[i].Naverage; j++) {
    61       // XXX accumulate all secfilt values in a single pass?
    62 
    63       int minUbercalDist = 1000;
    64 
    65       int Ns;
    66       for (Ns = 0; Ns < Nphotcodes; Ns++) {
    67 
    68         int thisCode = photcodes[Ns][0].code;
    69         int Nsec = GetPhotcodeNsec(thisCode);
    70 
    71         /* calculate the average mag in this SEC photcode for a single star */
    72 
    73         // skip bad stars
    74         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
    75 
    76         N = 0;
    77         m = catalog[i].averageT[j].measureOffset;
    78         for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    79 
    80           // skip measurements that do not match the current photcode
    81           PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode);
    82           if (!code) continue;
    83           if (code->equiv != thisCode) { continue; }
    84 
    85           if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; }
    86 
    87           if (getImageEntry (m, i) < 0) {
    88             Mcal = Mmos = Mgrid = 0;
    89           } else {
    90             Mcal  = getMcal  (m, i, flatcorr, catalog);
    91             if (isnan(Mcal)) { Ncal ++; continue; }
    92             Mmos  = getMmos  (m, i);
    93             if (isnan(Mmos)) { Nmos ++; continue; }
    94             Mgrid = getMgrid (m, i);
    95             if (isnan(Mgrid)) { Ngrid++; continue; }
    96           }
    97 
    98           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
    99           if (isnan(Msys)) { Nsys++; continue; }
    100           list[N] = Msys - Mcal - Mmos - Mgrid;
    101 
    102           int myUbercalDist = getUbercalDist(m,i);
    103           minUbercalDist = MIN(minUbercalDist, myUbercalDist);
    104 
    105           // dlist gives the error, which is used as the weight in WT_MEAN.
    106           // we can modify the resulting weight in a few ways:
    107           // 1) MIN_ERROR guarantees a floor
    108           // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    109           // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    110           // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    111           dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR);
    112 
    113           // up-weight the ubercal values (or convergence can take a long time...)
    114           if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    115             dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR);
    116           }
    117 
    118           // tie down reference photometry if the -refcode (code) option is selected
    119           // eg, -refcode g_SDSS
    120           // this probably makes no sense in the context of multifilter analysis
    121           if (refPhotcode) {
    122             if (code->code == refPhotcode->code) {
    123               // tiny error -> large weight
    124               // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR);
    125               dlist[N] = 0.0001;
    126             }
    127           }
    128           N++;
    129         }
    130 
    131         // when performing the grid analysis, STAR_TOOFEW will be set to 1;
    132         if (N <= STAR_TOOFEW) { /* too few measurements */
    133           // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].measureT[catalog[i].averageT[j].measureOffset].imageID, N, STAR_TOOFEW);
    134           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
    135           Nfew ++;
    136         } else {
    137           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
    138         }       
    139 
    140         liststats (list, dlist, N, &stats);
    141        
    142         catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
    143         catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
    144         catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
    145        
    146         catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
    147       }
    148     }
    149   }
    150   fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys);
    151 
    152   return (TRUE);
    153 }
    154 
    155 int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
    156 
    157   off_t k;
    158 
    159   int Nsecfilt = GetPhotcodeNsecfilt ();
    160 
    161   off_t m = average[0].measureOffset;
    162 
    163   for (k = 0; k < average[0].Nmeasure; k++, m++) {
    164     fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
    165   }
    166 
    167   int Ns;
    168   for (Ns = 0; Ns < Nsecfilt; Ns++) {
    169     fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
    170   }
    171   return 1;
    172 }
    173 
    174 # define MARK_SKIP_MEAS \
    175   catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \
    176   catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
     71void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
     72  results->Nfew  = 0;
     73  results->Ncode  = 0;
     74  results->Nsys  = 0;
     75  results->Nbad  = 0;
     76  results->Ncal  = 0;
     77  results->Nmos  = 0;
     78  results->Ngrid = 0;
     79  if (allocLists) {
     80    ALLOCATE (results->list,  double, Nmax);
     81    ALLOCATE (results->dlist, double, Nmax);
     82    ALLOCATE (results->wlist, double, Nmax);
     83  }
     84}
     85
     86void SetMrelInfoFree (SetMrelInfo *results) {
     87  free (results->list);
     88  free (results->dlist);
     89  free (results->wlist);
     90}
     91
     92void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
     93  summary->Nfew  += results->Nfew ;
     94  summary->Ncode  += results->Ncode ;
     95  summary->Nsys  += results->Nsys ;
     96  summary->Nbad  += results->Nbad ;
     97  summary->Ncal  += results->Ncal ;
     98  summary->Nmos  += results->Nmos ;
     99  summary->Ngrid += results->Ngrid;
     100}
     101
     102// mutex to lock setMrel_worker operations
     103static pthread_mutex_t setMrel_mutex = PTHREAD_MUTEX_INITIALIZER;
     104static int nextCatalog = 0;
     105
     106// we have an array of catalogs (catalog, Ncatalog).  we need to hand out catalogs one at a time to
     107// the worker threads as they need
     108off_t getNextCatalogForThread (int Ncatalog) {
     109
     110  pthread_mutex_lock (&setMrel_mutex);
     111  if (nextCatalog >= Ncatalog) {
     112    pthread_mutex_unlock (&setMrel_mutex);
     113    return (-1);
     114  }
     115  int thisCatalog = nextCatalog;
     116  nextCatalog ++;
     117
     118  pthread_mutex_unlock (&setMrel_mutex);
     119  return (thisCatalog);
     120}
    177121
    178122// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
    … …  
    180124// * setMrelOutput skips stars for which there are too few good measurements
    181125// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass
    182 
    183 // setMave is also similar to the above.  but:
    184 // * setMave is called by relphot -update-objects
    185 // * setMave excludes all detections with (PSF_QF < 0.85), setMrelOutput allows these for PASS > 2
    186 // * setMave updates 2MASS average flags
    187 // * setMave updates average EXT flags (PS1 and 2MASS)
     126// * setMrelOutput updates 2MASS average flags
     127// * setMrelOutput updates average EXT flags (PS1 and 2MASS)
     128
     129int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
     130
     131  int i;
     132
     133  if (NTHREADS) {
     134    int status = setMrel_threaded (catalog, Ncatalog, flatcorr);
     135    return status;
     136  }
     137
     138  int Nsecfilt = GetPhotcodeNsecfilt ();
     139
     140  SetMrelInfo summary, results;
     141  SetMrelInfoInit (&summary, FALSE);
     142  SetMrelInfoInit (&results, TRUE);
     143
     144  for (i = 0; i < Ncatalog; i++) {
     145    // pass == -1 for anything other than the final pass
     146    setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
     147    SetMrelInfoAccum (&summary, &results);
     148  }
     149  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     150
     151  SetMrelInfoFree (&results);
     152
     153  return (TRUE);
     154}
     155
    188156int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {
    189157
    190   int i, N;
    191   off_t j, k, m, Nmax;
    192   float Msys, Mcal, Mmos, Mgrid, Map;
    193   double *list, *dlist, *aplist, *daplist;
     158  int i;
     159
     160  int Nsecfilt = GetPhotcodeNsecfilt ();
     161
     162  SetMrelInfo summary, results;
     163  SetMrelInfoInit (&summary, FALSE);
     164  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
     165  ALLOCATE (results.aplist, double, Nmax);
     166  ALLOCATE (results.daplist, double, Nmax);
     167
     168  for (i = 0; i < Ncatalog; i++) {
     169    setMrel_catalog  (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options
     170    SetMrelInfoAccum (&summary, &results);
     171  }
     172  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     173
     174  SetMrelInfoFree (&results);
     175  free (results.aplist);
     176  free (results.daplist);
     177  return (TRUE);
     178}
     179
     180int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
     181
     182  int i;
     183
     184  SetMrelInfo summary;
     185  SetMrelInfoInit (&summary, FALSE);
     186
     187  pthread_attr_t attr;
     188  pthread_attr_init (&attr);
     189  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
     190 
     191  pthread_t *threads;
     192  ALLOCATE (threads, pthread_t, NTHREADS);
     193
     194  ThreadInfo *threadinfo;
     195  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
     196
     197  // each time this function is called, we cycle through the available catalogs.
     198  // make sure we start at 0
     199  nextCatalog = 0;;
     200
     201  // launch N worker threads
     202  for (i = 0; i < NTHREADS; i++) {
     203    threadinfo[i].entry = i;
     204    threadinfo[i].state = THREAD_RUN;
     205    threadinfo[i].catalog  =  catalog;
     206    threadinfo[i].Ncatalog = Ncatalog;
     207    threadinfo[i].flatcorr = flatcorr;
     208    SetMrelInfoInit (&threadinfo[i].summary, FALSE);
     209    pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]);
     210  }
     211  pthread_attr_destroy (&attr);
     212
     213  // wait until all threads have finished
     214  while (1) {
     215    int allDone = TRUE;
     216    for (i = 0; i < NTHREADS; i++) {
     217      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
     218    }
     219    if (allDone) {
     220      break;
     221    }
     222    usleep (500000);
     223  }
     224
     225  // all threads are done, free the threads array and grab the info
     226  free (threads);
     227 
     228  // report stats & summary from the threads
     229  for (i = 0; i < NTHREADS; i++) {
     230    if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n",
     231             i,
     232             threadinfo[i].summary.Ncode,
     233             threadinfo[i].summary.Nfew,
     234             threadinfo[i].summary.Nbad,
     235             threadinfo[i].summary.Ncal,
     236             threadinfo[i].summary.Nmos,
     237             threadinfo[i].summary.Ngrid,
     238             threadinfo[i].summary.Nsys);
     239    SetMrelInfoAccum (&summary, &threadinfo[i].summary);
     240  }
     241  if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     242  free (threadinfo);
     243
     244  return TRUE;
     245}
     246
     247void *setMrel_worker (void *data) {
     248
     249  ThreadInfo *threadinfo = data;
     250
     251  int Nsecfilt = GetPhotcodeNsecfilt ();
     252
     253  SetMrelInfo results;
     254  SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here
     255
     256  while (1) {
     257
     258    off_t i = getNextCatalogForThread(threadinfo->Ncatalog);
     259    if (i == -1) {
     260      threadinfo->state = THREAD_DONE;
     261      return NULL;
     262    }
     263
     264    Catalog *catalog = threadinfo->catalog;
     265    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
     266
     267    // pass == -1 for anything other than the final pass
     268    setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
     269    SetMrelInfoAccum (&threadinfo->summary, &results);
     270  }
     271
     272  SetMrelInfoFree (&results);
     273  return NULL;
     274}
     275
     276# define SKIP_THIS_MEAS(REASON) {                               \
     277    catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;     \
     278    if (catalog[Nc].measure) {                                  \
     279      catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;   \
     280    }                                                           \
     281    results->REASON ++;                                         \
     282    continue; }
     283
     284
     285// setMrel_catalog is used in 3 different contexts:
     286
     287// * during the main relphot iterations:
     288// ** operations only apply to measureT / averageT
     289// ** we are applying the analysis to the bright subset catalog
     290// ** we skip any stars found to be bad (STAR_BAD)
     291
     292// * during the final pass
     293// ** operations are applied to all objects
     294// ** we have special tests for PS1 extended data, 2MASS extended & good data, & synthetic photometry
     295
     296// * during the stand-alone '-averages' mode (relphot_objects)
     297// ** no image data is loaded
     298// ** we only use the provided calibration (measure.Mcal)
     299// ** we skip outlier rejection of measurements
     300
     301// set the Mrel values for the specified catalog
     302// NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array
     303int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
     304
     305  off_t j, k, m;
     306  int N;
     307  float Msys, Mcal, Mmos, Mgrid;
     308
    194309  StatType stats, apstats;
    195   int Nsec;
    196 
    197   int Nsecfilt = GetPhotcodeNsecfilt ();
    198 
    199   /* Nmeasure is now different, need to reallocate */
    200   Nmax = 0;
    201   for (i = 0; i < Ncatalog; i++) {
    202     for (j = 0; j < catalog[i].Naverage; j++) {
    203       Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
    204     }
    205   }
    206   ALLOCATE (list, double, MAX (1, Nmax));
    207   ALLOCATE (dlist, double, MAX (1, Nmax));
    208 
    209   ALLOCATE (aplist, double, MAX (1, Nmax));
    210   ALLOCATE (daplist, double, MAX (1, Nmax));
    211 
    212   for (i = 0; i < Ncatalog; i++) {
    213     for (j = 0; j < catalog[i].Naverage; j++) {
    214 
    215       if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) {
    216         fprintf (stderr, "test obj\n");
    217         print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure);
    218       }
     310  liststats_setmode (&stats, STATMODE);
     311  liststats_setmode (&apstats, STATMODE);
     312
     313  double *list    = results->list;
     314  double *dlist   = results->dlist;
     315  double *wlist   = results->wlist;
     316  double *aplist  = results->aplist;
     317  double *daplist = results->daplist;
     318
     319  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
     320
     321  int isSetMrelFinal = (pass >= 0);
     322
     323  for (j = 0; j < catalog[Nc].Naverage; j++) {
     324    // XXX accumulate all secfilt values in a single pass?
     325
     326    // option for a test print
     327    if (FALSE && (catalog[Nc].average[j].objID == 0x46a4) && (catalog[Nc].average[j].catID == 0xf40e)) {
     328      fprintf (stderr, "test obj\n");
     329      print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure);
     330    }
     331
     332    int GoodPS1 = FALSE;
     333    int Good2MASS = FALSE;
     334    int Galaxy2MASS = FALSE;
     335
     336    int NextPS1 = 0;
     337    int NpsfPS1 = 0;
     338
     339    int Ns;
     340    for (Ns = 0; Ns < Nphotcodes; Ns++) {
     341
     342      int thisCode = photcodes[Ns][0].code;
     343      int Nsec = GetPhotcodeNsec(thisCode);
     344
     345      /* calculate the average mag in this SEC photcode for a single star */
     346
     347      /* star/photcodes already calibrated */
     348      if ( isSetMrelFinal && catalog[Nc].found[Nsecfilt*j+Nsec]) continue; 
     349     
     350      // skip bad stars
     351      if (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue;
     352
     353      int Ncode = 0;
     354      int Next = 0;
     355      int haveSynth = FALSE;
     356
     357      int forceSynth = FALSE;
     358      int forceSynthEntry = -1;
     359
     360      int haveUbercal = FALSE;
    219361
    220362      int minUbercalDist = 1000;
    221 
    222       int Ns;
    223       for (Ns = 0; Ns < Nphotcodes; Ns++) {
    224 
    225         int thisCode = photcodes[Ns][0].code;
    226         Nsec = GetPhotcodeNsec(thisCode);
    227 
    228         /* star/photcodes already calibrated */
    229         if (catalog[i].found[Nsecfilt*j+Nsec]) continue; 
    230 
    231         int Ncode = 0;
    232         int Next = 0;
    233         int haveSynth = FALSE;
    234         int haveUbercal = FALSE;
    235         int isRefPhot = FALSE;
    236 
    237         N = 0;
    238         m = catalog[i].averageT[j].measureOffset;
    239         for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    240 
    241           // skip measurements that do not match the current photcode
    242           PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode);
    243           if (!code) continue;
    244           if (code->equiv != thisCode) { continue; }
    245           Ncode ++;
    246 
    247           if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
    248 
    249           if (getImageEntry (m, i) < 0) {
    250             // these detetions have no image (eg, ref values such as 2MASS)
    251             Mcal = Mmos = Mgrid = 0;
    252             isRefPhot = TRUE;
    253           } else {
    254             Mcal  = getMcal  (m, i, flatcorr, catalog);
    255             if (isnan(Mcal)) {
    256               MARK_SKIP_MEAS;
    257               continue;
    258             }
    259             Mmos  = getMmos  (m, i);
    260             if (isnan(Mmos)) {
    261               MARK_SKIP_MEAS;
    262               continue;
    263             }
    264             Mgrid = getMgrid (m, i);
    265             if (isnan(Mgrid)) {
    266               MARK_SKIP_MEAS;
    267               continue;
    268             }
    269           }
    270 
    271           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
    272           if (isnan(Msys)) {
    273             if (!isRefPhot) {
    274               MARK_SKIP_MEAS;
    275             }
    276             continue;
    277           }
    278           list[N] = Msys - Mcal - Mmos - Mgrid;
    279 
    280           Map = PhotAper (&catalog[i].measure[m]);
     363   
     364      N = 0;
     365      m = catalog[Nc].averageT[j].measureOffset;
     366      for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) {
     367
     368        // skip measurements that do not match the current photcode
     369        PhotCode *code = GetPhotcodebyCode (catalog[Nc].measureT[m].photcode);
     370        if (!code) continue;
     371        if (code->equiv != thisCode) { continue; }
     372        Ncode ++;
     373
     374        if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
     375
     376        if (getImageEntry (m, Nc) < 0) {
     377          // measurements without an image are either external reference photometry or
     378          // data for which the associated image has not been loaded (probably because of
     379          // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
     380          Mmos = Mgrid = 0;
     381          Mcal = catalog[Nc].measureT[m].Mcal; // check that this is zero for loaded REF value
     382        } else {
     383          Mcal  = getMcal  (m, Nc, flatcorr, catalog);
     384          if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
     385          Mmos  = getMmos  (m, Nc);
     386          if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
     387          Mgrid = getMgrid (m, Nc);
     388          if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
     389        }
     390
     391        Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]);
     392        if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
     393
     394        list[N] = Msys - Mcal - Mmos - Mgrid;
     395
     396        int myUbercalDist = getUbercalDist(m,Nc);
     397        minUbercalDist = MIN(minUbercalDist, myUbercalDist);
     398
     399        if (isSetMrelFinal) {
     400          float Map = PhotAper (&catalog[Nc].measure[m]);
    281401          aplist[N] = Map - Mcal - Mmos - Mgrid;
    282402
    283           int myUbercalDist = getUbercalDist(m,i);
    284           minUbercalDist = MIN(minUbercalDist, myUbercalDist);
    285 
    286           // count the extended detections
    287           if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
    288             if (!isnan(catalog[i].measure[m].Map)) {
    289               if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
     403          // special options for PS1 data
     404          if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
     405            // count the extended detections
     406            if (!isnan(catalog[Nc].measure[m].Map)) {
     407              if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) {
    290408                Next ++;
     409                NextPS1 ++;
     410              } else {
     411                NpsfPS1 ++;
    291412              }
    292413            }
    293414          }
    294415          // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
    295           if ((catalog[i].measure[m].photcode >= 2011) && (catalog[i].measure[m].photcode <= 2013)) {
    296             if (catalog[i].measure[m].photFlags & 0x00c00000) {
     416          if ((catalog[Nc].measure[m].photcode >= 2011) && (catalog[Nc].measure[m].photcode <= 2013)) {
     417            if (catalog[Nc].measure[m].photFlags & 0x00c00000) {
    297418              Next ++;
     419              Galaxy2MASS = TRUE;
    298420            }
    299             if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) {
    300               // detections without one of these bits should only be used in PASS_1
    301               MARK_SKIP_MEAS;
    302               continue;
    303             }
    304           }
    305           // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars)
    306           if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) {
    307             if (pass < 4) {
    308               MARK_SKIP_MEAS;
    309               continue;
    310             }
    311             haveSynth = TRUE;
    312           }
    313 
    314           // dlist gives the error, which is used as the weight in WT_MEAN.
    315           // we can modify the resulting weight in a few ways:
    316           // 1) MIN_ERROR guarantees a floor
    317           // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    318           // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    319           // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    320           dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR);
    321 
    322           // up-weight the ubercal values (or convergence can take a long time...) (XXX make this optional?)
    323           if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    324             dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR);
    325             haveUbercal = TRUE;
    326           }
    327 
    328           // tie down reference photometry if the -refcode (code) option is selected
    329           // eg, -refcode g_SDSS
    330           // this probably makes no sense in the context of multifilter analysis
    331           if (refPhotcode) {
    332             if (code->code == refPhotcode->code) {
    333               // tiny error -> large weight
    334               // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR);
    335               dlist[N] = 0.0001;
    336             }
    337           }
    338           N++;
    339         }
    340         if (N < 1) continue;
    341 
    342         for (m = 0; m < N; m++) {
    343           daplist[m] = dlist[m];
    344         }
    345 
    346         // XXX force WT_MEAN or MEAN here?
    347         liststats (list, dlist, N, &stats);
    348         liststats (aplist, daplist, N, &apstats);
    349         catalog[i].found[Nsecfilt*j+Nsec] = TRUE;
    350         switch (pass) {
    351           case 0:
    352             catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
    353             break;
    354           case 1:
    355             catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
    356             break;
    357           case 2:
    358             catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
    359             break;
    360           case 3:
    361             catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
    362             break;
    363           case 4:
    364             catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
    365             break;
    366         }
    367 
    368         if (haveSynth) {
    369           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
    370         }       
    371         if (haveUbercal) {
    372           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL;
    373         }       
    374 
    375         /* use sigma or error in dM for output? */
    376         catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
    377         catalog[i].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean;
    378         catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma);
    379         catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
    380         catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
    381         catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
    382 
    383         catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
    384         catalog[i].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
    385         catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
    386 
    387         if ((Next > 0) && (Next > 0.5*N)) {
    388           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT;
    389         }
    390       }
    391     }
    392   }
    393 
    394   free (list);
    395   free (dlist);
    396 
    397   free (aplist);
    398   free (daplist);
    399   return (TRUE);
    400 }
    401 
    402 /* grab Nsec for named photcode */
    403 # define NAMED_PHOTCODE_NSEC(MY_NSEC,NAME)      \
    404   short MY_NSEC = -1;                           \
    405   {                                             \
    406   PhotCode *code = GetPhotcodebyName (NAME);    \
    407   if (code) {                                   \
    408     MY_NSEC = GetPhotcodeNsec (code->equiv);    \
    409   } }
    410 
    411 // For each average object, set the average mags based on existing equiv photometry.
    412 // NOTE: this function operates on the real Measure & Average structures, not the
    413 // MeasureTiny & AverageTiny structures
    414 // NOTE: this function is called on the remote machine -- make sure any (global) options
    415 // are passed to the relphot_client program
    416 int setMave (Catalog *catalog, int Ncatalog) {
    417 
    418   off_t j, k, m, Nmax;
    419   int i, Ns, Nsecfilt, N, Nc;
    420   float Msys;
    421   double *list, *dlist;
    422   StatType stats;
    423   DVOAverageFlags flagBits;
    424 
    425   flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
    426 
    427   // pre-allocate a list for stats purposes
    428   Nmax = 0;
    429   for (i = 0; i < Ncatalog; i++) {
    430     for (j = 0; j < catalog[i].Naverage; j++) {
    431       Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
    432     }
    433   }
    434   ALLOCATE (list, double, MAX (1, Nmax));
    435   ALLOCATE (dlist, double, MAX (1, Nmax));
    436 
    437   Nsecfilt = GetPhotcodeNsecfilt ();
    438 
    439   // we want to raise some bits on the 2MASS (JHK) secfilt flags, if we have 2MASS data
    440   NAMED_PHOTCODE_NSEC (Nsec_J, "2MASS_J");
    441   NAMED_PHOTCODE_NSEC (Nsec_H, "2MASS_H");
    442   NAMED_PHOTCODE_NSEC (Nsec_K, "2MASS_K");
    443 
    444   for (i = 0; i < Ncatalog; i++) {
    445     for (j = 0; j < catalog[i].Naverage; j++) {
    446 
    447       // update average photometry for each of the average filters
    448 
    449       // XXX Note that this would be faster if we had an array of results and accumulated
    450       // them in a single pass
    451 
    452       int Next = 0;
    453 
    454       for (Ns = 0; Ns < Nsecfilt; Ns++) {
    455 
    456         PhotCode *code = GetPhotcodebyNsec (Ns);
    457         Nc = code[0].code;
    458        
    459         N = 0;
    460         m = catalog[i].average[j].measureOffset;
    461         for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    462           if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
    463           if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
    464 
    465           Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    466           if (isnan(Msys)) continue;
    467 
    468           // reject POOR detections (PSF_QF < 0.85) and count extended detections
    469           // XXX only apply this filter for psphot data from GPC1 for now...
    470           if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) {
    471             if (catalog[i].measure[m].psfQual < 0.85) continue;
    472             if (!isnan(catalog[i].measure[m].Map)) {
    473               if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
    474                 Next ++;
     421            if (pass == 0) {
     422              if (catalog[Nc].measure[m].photFlags & 0x00000007) {
     423                Good2MASS = TRUE;
     424              } else {
     425                // detections without one of these bits should only be used in PASS_1
     426                SKIP_THIS_MEAS(Nbad);
     427                continue;
    475428              }
    476429            }
    477430          }
    478431
    479           // XXX make it optional to apply Mcal (do I need an extra field for advisory Mcal?)
    480           list[N] = Msys - catalog[i].measure[m].Mcal;
    481           dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    482           N++;
    483         }
    484         if (N < 1) continue;
    485 
    486         liststats (list, dlist, N, &stats);
    487 
    488         /* use sigma or error in dM for output? */
    489         catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
    490         catalog[i].secfilt[Nsecfilt*j+Ns].dM = MAX (stats.error, stats.sigma);
    491         catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
    492         catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
    493         catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
    494 
    495         catalog[i].secfilt[Nsecfilt*j+Ns].M_80 = 1000 * stats.Upper80;
    496         catalog[i].secfilt[Nsecfilt*j+Ns].M_20 = 1000 * stats.Lower20;
    497 
    498         if ((Next > 0) && (Next > 0.5*N)) {
    499           catalog[i].secfilt[Nsecfilt*j+Ns].flags |= ID_SECF_OBJ_EXT;
    500         }
    501       }
    502 
    503       // update average flags based on the detection stats. 
    504 
    505       int Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000])
    506       int goodPS1 = FALSE;     // true if any PS1 measurements have psfQual > 0.85
    507       int good2MASS = FALSE;   // true if 2MASS J measurements have significant detections
    508       int nEXT = 0;
    509       int nPSF = 0;     // number of PS1 PSF vs EXT measurements
    510       int have2MASS = FALSE;
    511      
    512       // count, flag good and extended detections
    513       m = catalog[i].average[j].measureOffset;
    514       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    515 
    516         // PS1 data :
    517         if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
    518           if (catalog[i].measure[m].psfQual > 0.85) {
    519             goodPS1 = TRUE;
    520             if (!isnan(catalog[i].measure[m].Map)) {
    521               if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
    522                 nEXT ++;
    523               } else {
    524                 nPSF ++;
     432          // Blindly accepth the SYNTH mags if we are above saturation, otherwise,
     433          // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
     434          if ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) {
     435            // something of a hack: force object to use synth values if synth mags >>
     436            // saturation (3pi instrumental mags < -15)
     437            float MaxMagForceSynth = NAN;
     438            switch (catalog[Nc].measure[m].photcode) {
     439              case 3001:
     440                MaxMagForceSynth = 13.64;
     441                break;
     442              case 3002:
     443                MaxMagForceSynth = 13.76;
     444                break;
     445              case 3003:
     446                MaxMagForceSynth = 13.74;
     447                break;
     448              case 3004:
     449                MaxMagForceSynth = 12.94;
     450                break;
     451              case 3005:
     452                MaxMagForceSynth = 12.01;
     453                break;
     454            }
     455            if (catalog[Nc].measureT[m].M < MaxMagForceSynth) {
     456              forceSynth = TRUE;
     457              forceSynthEntry = N;
     458            } else {
     459              if (pass < 4) {
     460                SKIP_THIS_MEAS(Nbad);
     461                continue;
    525462              }
     463              haveSynth = TRUE;
    526464            }
    527465          }
    528466        }
    529          
    530         // 2MASS data J-band flags
    531         if (catalog[i].measure[m].photcode == 2011) {
    532           // only need to do this once (always have JHK triplet; galaxy flag is same for all 3)
    533           have2MASS = TRUE;
    534           if (catalog[i].measure[m].photFlags & 0x00c00000) {
    535             Galaxy2MASS = TRUE;
     467
     468        // dlist gives the error per measurement, wlist gives the weight
     469        // we can modify the error and weight in a few ways:
     470        // 1) MIN_ERROR guarantees a floor
     471        // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
     472        // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
     473        // 4) some reference photcode of some kind can be specified as fixed and have a high weight
     474        dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
     475        wlist[N] = 1.0;
     476
     477        // up-weight the ubercal values (or convergence can take a long time...)
     478        if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     479          wlist[N] = 10.0;
     480        }
     481
     482        // tie down reference photometry if the -refcode (code) option is selected
     483        // eg, -refcode g_SDSS
     484        // this probably makes no sense in the context of multifilter analysis
     485        if (refPhotcode) {
     486          if (code->code == refPhotcode->code) {
     487            wlist[N] = 100.0;
    536488          }
    537           if (catalog[i].measure[m].photFlags & 0x00000007) {
    538             good2MASS = TRUE;
    539             if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_0;
    540           } else {
    541             if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_1;
    542           }
    543         } 
    544         // 2MASS data H-band flags
    545         if (catalog[i].measure[m].photcode == 2012) {
    546           if (catalog[i].measure[m].photFlags & 0x00000007) {
    547             good2MASS = TRUE;
    548             if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_0;
    549           } else {
    550             if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_1;
    551           }
    552         } 
    553         // 2MASS data K-band flags
    554         if (catalog[i].measure[m].photcode == 2013) {
    555           if (catalog[i].measure[m].photFlags & 0x00000007) {
    556             good2MASS = TRUE;
    557             if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_0;
    558           } else {
    559             if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_1;
    560           }
    561         } 
    562       }
     489        }
     490        N++;
     491      }
     492
     493      int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
     494
     495      // when performing the grid analysis, STAR_TOOFEW should be set to 1;
     496      if (N < Nminmeas) { /* too few measurements */
     497        // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D, catalog[Nc].measureT[catalog[Nc].averageT[j].measureOffset].imageID, N, STAR_TOOFEW);
     498        catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
     499        if (Ncode == 0) {
     500          results->Ncode ++;
     501        } else {
     502          results->Nfew ++;
     503        }
     504        continue;
     505      } else {
     506        catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
     507      }
     508
     509      if (forceSynth) {
     510        // use the single SYNTH value instead of the other mags here
     511        myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
     512        list[0]  = list[forceSynthEntry];
     513        dlist[0] = dlist[forceSynthEntry];
     514        wlist[0] = wlist[forceSynthEntry];
     515        N = 1;
     516      }
     517      liststats (list, dlist, wlist, N, &stats);
     518
     519      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
     520      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error;
     521      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
     522
     523      // when running -averages, we have no information about the images, so we cannot set this
     524      if (minUbercalDist > -1) {
     525        catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
     526      }
     527
     528      if (isSetMrelFinal) {
     529        catalog[Nc].found[Nsecfilt*j+Nsec] = TRUE;
     530
     531        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
     532        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
     533        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
     534
     535        catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
     536        catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
     537
     538        // NOTE : use the modified weight for apmags as well as psf mags
     539        liststats (aplist, daplist, wlist, N, &apstats);
     540
     541        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean;
     542
     543        // NOTE: for 2MASS measurements, Next should be 1, as should N
     544        if ((Next > 0) && (Next > 0.5*N)) {
     545          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT;
     546        }
     547
     548        switch (pass) {
     549          case 0:
     550            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
     551            GoodPS1 = TRUE;
     552            break;
     553          case 1:
     554            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
     555            GoodPS1 = TRUE;
     556            break;
     557          case 2:
     558            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
     559            GoodPS1 = TRUE;
     560            break;
     561          case 3:
     562            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
     563            break;
     564          case 4:
     565            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
     566            break;
     567        }
     568        if (haveSynth) {
     569          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
     570        }       
     571        if (haveUbercal) {
     572          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL;
     573        }       
     574      }
     575    }
     576
     577    if (isSetMrelFinal) {
     578      DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
    563579
    564580      // we attempt to set a few flags here; reset those bits before trying:
    565       catalog[i].average[j].flags &= ~flagBits;
    566 
    567       // XXX set the secfilt bits?
    568       if (nEXT && (nEXT > nPSF)) {
    569         catalog[i].average[j].flags |= ID_OBJ_EXT;
    570       }
    571       if (goodPS1) {
    572         catalog[i].average[j].flags |= ID_OBJ_GOOD;
     581      catalog[Nc].average[j].flags &= ~flagBits;
     582
     583      if (NextPS1 && (NextPS1 > NpsfPS1)) {
     584        catalog[Nc].average[j].flags |= ID_OBJ_EXT;
     585      }
     586      if (GoodPS1) {
     587        catalog[Nc].average[j].flags |= ID_OBJ_GOOD;
    573588      }
    574589      if (Galaxy2MASS) {
    575         catalog[i].average[j].flags |= ID_OBJ_EXT_ALT;
    576         if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_SECF_OBJ_EXT;
    577         if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_SECF_OBJ_EXT;
    578         if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_SECF_OBJ_EXT;
    579       }
    580       if (good2MASS) {
    581         catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT;
    582       }
    583     }
    584   }
    585 
    586   free (list);
    587   free (dlist);
     590        catalog[Nc].average[j].flags |= ID_OBJ_EXT_ALT;
     591      }
     592      if (Good2MASS) {
     593        catalog[Nc].average[j].flags |= ID_OBJ_GOOD_ALT;
     594      }
     595    }
     596  }
    588597  return (TRUE);
     598}
     599
     600int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
     601
     602  off_t k;
     603
     604  int Nsecfilt = GetPhotcodeNsecfilt ();
     605
     606  off_t m = average[0].measureOffset;
     607
     608  for (k = 0; k < average[0].Nmeasure; k++, m++) {
     609    fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
     610  }
     611
     612  int Ns;
     613  for (Ns = 0; Ns < Nsecfilt; Ns++) {
     614    fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
     615  }
     616  return 1;
    589617}
    590618
    … …  
    612640        if (isnan(Mgrid)) continue;
    613641
     642        // note that measurements for which the image is not selected will not be modified
     643        // (that is a good thing! -- we keep a prior calibration)
     644
    614645        // set the output calibration
    615646        catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
    … …  
    630661  double Chisq, MaxScatter, MaxChisq;
    631662  double *xlist, *slist, *dlist;
     663
    632664  StatType stats;
     665  liststats_setmode (&stats, "MEAN");
    633666
    634667  if (VERBOSE) fprintf (stderr, "marking poor stars\n");
    … …  
    644677  int Nsecfilt = GetPhotcodeNsecfilt ();
    645678
    646   // XX int oldPLOTSTUFF = PLOTSTUFF;
    647   // XX PLOTSTUFF = TRUE;
    648   // XX plot_chisq (catalog, Ncatalog);
    649   // XX PLOTSTUFF = oldPLOTSTUFF;
    650 
    651679  // eliminate bad stars using the stats for a single secfilt at a time
    652   // XXX DEP replace average.flags with secfilt flags
    653680  for (Ns = 0; Ns < Nphotcodes; Ns ++) {
    654681   
    … …  
    669696    }
    670697 
    671     initstats ("MEAN");
    672     liststats (xlist, dlist, Ntot, &stats);
     698    liststats (xlist, dlist, NULL, Ntot, &stats);
    673699    MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
    674     liststats (slist, dlist, Ntot, &stats);
     700
     701    liststats (slist, dlist, NULL, Ntot, &stats);
    675702    MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
    676703    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
    … …  
    696723    }
    697724    fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
    698     initstats (STATMODE);
    699725  }
    700726  free (xlist);
    … …  
    702728  free (dlist);
    703729}
     730
     731// clean_measures examines the stats for a single star.  It measures the INNER 50% mean
     732// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
     733// sigma of the INNER 50% mean.  it then flags any measurements which are more than
     734// NSIGMA_REJECT (5) sigma of the mean
    704735
    705736# define NSIGMA_CLIP 3.0
    … …  
    712743  double *tlist, *list, *dlist;
    713744  float Msys, Mcal, Mmos, Mgrid;
    714   StatType stats;
    715745  int Ncal, Nmos, Ngrid, Nfew;
    716746
    … …  
    732762  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
    733763  TOOFEW = MAX (5, STAR_TOOFEW);
     764
     765  // stats structures for inner and full stats
     766  StatType instats, stats;
     767  liststats_setmode (&instats, "INNER_MEAN");
     768  liststats_setmode (&stats, "MEAN");
    734769
    735770  Ndel = Nave = 0;
    … …  
    779814
    780815        // calculated mean of inner 50%
    781         initstats ("INNER_MEAN");
    782         liststats (list, dlist, N, &stats);
    783         stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
     816        liststats (list, dlist, NULL, N, &instats);
     817        instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */
    784818
    785819        // ignore entries > 3sigma from inner mean
    786820        for (k = m = 0; k < N; k++) {
    787           if (fabs (list[k] - stats.median) < NSIGMA_CLIP*stats.sigma) {
     821          if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
    788822            list[m] = list[k];
    789823            m++;
    … …  
    791825        }
    792826        // recalculate the mean & sigma of the accepted measurements
    793         initstats ("MEAN");
    794         liststats (list, dlist, m, &stats);
     827        liststats (list, dlist, NULL, m, &stats);
    795828        stats.sigma = MAX (MIN_ERROR, stats.sigma);
    796829
    … …  
    841874    }
    842875  }
    843   initstats (STATMODE);
    844876  if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
     877
    845878  free (list);
    846879  free (dlist);
    847880  free (ilist);
    848881  free (tlist);
    849 
    850 
    851882}
    852883
    … …  
    898929  }
    899930
    900   // fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0);
    901   liststats (list, dlist, n, &stats);
     931  liststats_setmode (&stats, "MEAN");
     932
     933  liststats (list, dlist, NULL, n, &stats);
    902934  free (list);
    903935  free (dlist);
    … …  
    938970  }
    939971
    940   liststats (list, dlist, n, &stats);
     972  liststats_setmode (&stats, "MEAN");
     973
     974  liststats (list, dlist, NULL, n, &stats);
    941975  free (list);
    942976  free (dlist);
    … …  
    9771011  }
    9781012
    979   liststats (list, dlist, n, &stats);
     1013  liststats_setmode (&stats, "MEAN");
     1014
     1015  liststats (list, dlist, NULL, n, &stats);
    9801016  free (list);
    9811017  free (dlist);

  • trunk/Ohana/src/relphot/src/args.c

    r33651 r33963  
    7474  }
    7575
     76  NTHREADS = 0;
     77  if ((N = get_argument (argc, argv, "-threads"))) {
     78    remove_argument (N, &argc, argv);
     79    NTHREADS = atof (argv[N]);
     80    remove_argument (N, &argc, argv);
     81  }
     82
    7683  // XXX for the moment, make this selection manual.  it needs to be automatic
    7784  // based on the state of the SkyTable
    … …  
    300307  }
    301308
     309  ApplyOffsets = FALSE;
     310  if ((N = get_argument (argc, argv, "-apply-offsets"))) {
     311    remove_argument (N, &argc, argv);
     312    ApplyOffsets = TRUE;
     313  }
     314
    302315  if (UpdateAverages && (argc == 1)) return TRUE;
    303316  if (argc != 2) relphot_usage ();
    … …  
    331344  }
    332345  if (!HOSTDIR) relphot_client_usage();
    333 
    334   // CATDIR[0] = 0;
    335   // if ((N = get_argument (argc, argv, "-catdir"))) {
    336   //   remove_argument (N, &argc, argv);
    337   //   strcpy (CATDIR, argv[N]);
    338   //   remove_argument (N, &argc, argv);
    339   // }
    340   // if (!CATDIR[0]) relphot_client_usage();
    341346
    342347  IMAGES = NULL; // used in -update mode

  • trunk/Ohana/src/relphot/src/bcatalog.c

    r33651 r33963  
    33int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
    44
    5 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
     5int bcatalog (Catalog *subcatalog, Catalog *catalog, int Ncat) {
    66 
    77  off_t i, j, offset;
    … …  
    7474        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map    = NAN;
    7575        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
     76        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT;
    7677        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm     = NAN;
    7778        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
    … …  
    158159      subcatalog[0].measureT[Nmeasure].averef = Naverage;
    159160      if (RESET) {
    160         if (!KEEP_UBERCAL) {
    161           subcatalog[0].measureT[Nmeasure].Mcal = 0;
    162           subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     161        // only reset Mcal for measures with a matching image
     162        // do not reset Mcal for ubercal images unless explicitly requested
     163        if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     164          if (!KEEP_UBERCAL) {
     165            subcatalog[0].measureT[Nmeasure].Mcal = 0.0;
     166            subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     167          }
     168        } else {
     169          if (getImageEntry (Nmeasure, Ncat) >= 0) {
     170            subcatalog[0].measureT[Nmeasure].Mcal = 0.0;
     171          }
    163172        }
    164173        subcatalog[0].measureT[Nmeasure].dbFlags &= ~PHOTOM_FLAGS;

  • trunk/Ohana/src/relphot/src/global_stats.c

    r33651 r33963  
    11# include "relphot.h"
     2
     3# define MARKTIME(MSG,...) { \
     4  float dtime; \
     5  gettimeofday (&stopTimer, (void *) NULL); \
     6  dtime = DTIME (stopTimer, startTimer); \
     7  fprintf (stderr, MSG, __VA_ARGS__); }
    28
    39void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    … …  
    511  StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD;
    612
    7   initstats ("MEAN");
     13  // struct timeval startTimer, stopTimer;
    814
    915  fprintf (stderr, "\n");
    1016  fprintf (stderr, "STATS            median     mean    sigma      min      max   Nmeas\n");
     17
     18  // gettimeofday (&startTimer, (void *) NULL);
    1119
    1220  int Ns;
    … …  
    1725    int seccode = photcodes[Ns][0].code;
    1826
     27    // MARKTIME("start star stats: %f sec\n", dtime);
    1928    stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr);
     29    // MARKTIME("done star N: %f sec\n", dtime);
    2030    stX = statsStarX (catalog, Ncatalog, Nsec);
     31    // MARKTIME("done star X: %f sec\n", dtime);
    2132    stS = statsStarS (catalog, Ncatalog, Nsec);
     33    // MARKTIME("done star S: %f sec\n", dtime);
    2234 
    2335    fprintf (stderr, "   --- stats for %s ---\n", photcodes[Ns][0].name);
    … …  
    4658  fprintf (stderr, "dMcal mosaic:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas);
    4759  fprintf (stderr, "chisq mosaic:   %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas);
    48  
    49 
    50   initstats (STATMODE);
    51 
    5260}
    5361

  • trunk/Ohana/src/relphot/src/initialize.c

    r33651 r33963  
    1212    PhotcodeList = strcreate (argv[1]);
    1313    photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes
    14   }
     14  } else {
     15    char tmpline1[256];
     16    int Ns;
     17    Nphotcodes = GetPhotcodeNsecfilt ();
     18    ALLOCATE (photcodes, PhotCode *, Nphotcodes);
     19    ALLOCATE (PhotcodeList, char, 256);
     20    for (Ns = 0; Ns < Nphotcodes; Ns++) {
     21      photcodes[Ns] = GetPhotcodebyNsec (Ns);
     22      if (Ns > 0) {
     23        snprintf (tmpline1, 256, "%s,%s", PhotcodeList, photcodes[Ns][0].name);
     24      } else {
     25        snprintf (tmpline1, 256, "%s", photcodes[Ns][0].name);
     26      }
     27      strcpy (PhotcodeList, tmpline1);
     28    }
     29  }   
     30   
    1531  if (USE_GRID && (Nphotcodes > 1)) {
    1632    fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n");
    1733    exit (1);
    1834  }
    19 
    20   initstats (STATMODE);
    2135
    2236  IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;

  • trunk/Ohana/src/relphot/src/liststats.c

    r33651 r33963  
    11# include "relphot.h"
    22
    3 enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
    4       M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
     3void liststats_setmode (StatType *stats, char *strmode) {
    54
    6 static int statmode;
     5  stats->statmode = -1;
     6  if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
     7  if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
     8  if (!strcmp (strmode, "WT_MEAN"))          { stats->statmode = STATS_WT_MEAN; return; }
     9  if (!strcmp (strmode, "INNER_MEAN"))       { stats->statmode = STATS_INNER_MEAN; return; }
     10  if (!strcmp (strmode, "INNER_WTMEAN"))     { stats->statmode = STATS_INNER_WTMEAN; return; }
     11  if (!strcmp (strmode, "CHI_INNER_MEAN"))   { stats->statmode = STATS_CHI_INNER_MEAN; return; }
     12  if (!strcmp (strmode, "CHI_INNER_WTMEAN")) { stats->statmode = STATS_CHI_INNER_WTMEAN; return; }
    713
    8 void initstats (char *mode) {
    9 
    10   statmode = -1;
    11   if (!strcmp (mode, "MEAN")) statmode = M_MEAN;
    12   if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN;
    13   if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN;
    14   if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN;
    15   if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN;
    16   if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN;
    17   if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN;
    18 
    19   if (statmode == -1) {
    20     fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode);
    21     exit (1);
    22   }
     14  fprintf (stderr, "ERROR: invalid stats mode: %s\n", strmode);
     15  exit (1);
    2316}
    2417
    25 int liststats (double *value, double *dvalue, int N, StatType *stats) {
     18int liststats (double *value, double *dvalue, double *weight, int N, StatType *stats) {
    2619 
    27   int i, ks, ke, Nm;
    28   double Mo, dMo, M, dM, X2, dS, *chi;
     20  int i, ks, ke;
     21  double Mo, dMo, M, dM, Nm, X2, dS, R, W, *chi;
     22
     23  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
    2924
    3025  ke = ks = dMo = 0;
    … …  
    3732  if (N < 1) return (FALSE);
    3833
    39   dsortpair (value, dvalue, N);
    40   stats[0].median = value[(int)(0.5*N)];
    41   stats[0].min    = value[0];
    42   stats[0].max    = value[N-1];
    4334  int N80 = MIN (N-1, 0.8*N);
    4435  int N20 = MAX (0, 0.2*N);
     36
     37  if (weight) {
     38    dsortthree (value, dvalue, weight, N);
     39  } else {
     40    dsortpair (value, dvalue, N);
     41  }
     42
     43  // these values do not depend on the errors or weighting scheme
     44  stats[0].median  = value[(int)(0.5*N)];
     45  stats[0].min     = value[0];
     46  stats[0].max     = value[N-1];
    4547  stats[0].Upper80 = value[N80];
    4648  stats[0].Lower20 = value[N20];
    4749
    48   switch (statmode) {
    49   case M_MEDIAN:
    50     ks = 0;
    51     ke = N;
    52     Mo = stats[0].median;
    53     Nm = N;
    54     goto chisq;
    55     break;
    56   case M_MEAN:
    57   case M_WT_MEAN:
    58     ks = 0;
    59     ke = N;
    60     break;
    61   case M_INNER_MEAN:
    62   case M_INNER_WTMEAN:
    63   case M_CHI_INNER_MEAN:
    64   case M_CHI_INNER_WTMEAN:
    65     ks = 0.25*N + 0.50;
    66     ke = 0.75*N + 0.25;
    67     if (N <= 3) {
     50  switch (stats->statmode) {
     51    case STATS_MEDIAN:
    6852      ks = 0;
    6953      ke = N;
    70     }
    71     break;
     54      Mo = stats[0].median;
     55      Nm = N;
     56      goto chisq;
     57      break;
     58    case STATS_MEAN:
     59    case STATS_WT_MEAN:
     60      ks = 0;
     61      ke = N;
     62      break;
     63    case STATS_INNER_MEAN:
     64    case STATS_INNER_WTMEAN:
     65    case STATS_CHI_INNER_MEAN:
     66    case STATS_CHI_INNER_WTMEAN:
     67      ks = 0.25*N + 0.50;
     68      ke = 0.75*N + 0.25;
     69      if (N <= 3) {
     70        ks = 0;
     71        ke = N;
     72      }
     73      break;
     74    case STATS_NONE:
     75      myAbort ("undefined stats");
    7276  }   
    7377
    74   if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
     78  // for these two modes, I need a vector of the chi-square contribution
     79  // I'm actually just using chisq to get the correct sorting order
     80  if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) {
    7581    ALLOCATE (chi, double, N);
    7682    for (i = 0; i < N; i++) {
    7783      chi[i] = (value[i] - stats[0].median) / dvalue[i];
    7884    }
    79     dsortthree (chi, value, dvalue, N);
     85    if (weight) {
     86      dsortfour (chi, value, dvalue, weight, N);
     87    } else {
     88      dsortthree (chi, value, dvalue, N);
     89    }
    8090    free (chi);
    8191  }
    8292
    83   /* calculating the per-star offset based on the weighted average */
    84   M = dM = Nm = 0;
    85   if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
    86     for (i = ks; i < ke; i++) {
    87       M   += value[i] / SQ (dvalue[i]);
    88       dM  += 1.0 / SQ (dvalue[i]);
    89       Nm  ++; 
    90     }   
    91     Mo = M / dM;
    92     dMo = sqrt (1.0 / dM);
     93  int WeightedMean = FALSE;
     94  WeightedMean |= (stats->statmode == STATS_WT_MEAN);
     95  WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
     96  WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
     97
     98  /* calculating the per-star offset based on the desired weighting scheme */
     99  M = dM = Nm = W = R = 0;
     100  if (weight) {
     101    // the weight value is multiplied by whichever nominal weighting scheme is provided
     102    // thus the user should set weight to 1.0 for nominal weight, and 100 for heavy weight (or so)
     103    // and 0.01 for under-weight
     104    if (WeightedMean) {
     105      for (i = ks; i < ke; i++) {
     106        M  += value[i] * weight[i] / SQ(dvalue[i]);
     107        W  +=            weight[i] / SQ(dvalue[i]);
     108        dM += SQ (weight[i] / dvalue[i]);
     109        R  += weight[i] / SQ(dvalue[i]);
     110        Nm += 1.0; 
     111      }
     112      Mo  = M / W;
     113      dMo = sqrt (dM) / R;
     114    } else {
     115      for (i = ks; i < ke; i++) {
     116        M  += value[i] * weight[i];
     117        W  +=            weight[i];
     118        dM += SQ (weight[i] * dvalue[i]);
     119        R  += weight[i];
     120        Nm += 1.0; 
     121      }
     122      Mo  = M / W;
     123      dMo = sqrt (dM) / R;
     124    }
    93125  } else {
    94     for (i = ks; i < ke; i++) {
    95       M   += value[i];
    96       dM  += SQ (dvalue[i]);
    97       Nm  ++; 
    98     }   
    99     Mo = M / (double) Nm;
    100     dMo = sqrt (dM / (double) Nm);
     126    // NULL weight vector is supplied, revert to standard form (above reverts if weight[i] == 1)
     127    if (WeightedMean) {
     128      // weighted by inverse-variance
     129      for (i = ks; i < ke; i++) {
     130        M   += value[i] / SQ (dvalue[i]);
     131        dM  += 1.0 / SQ (dvalue[i]);
     132        Nm  += 1.0; 
     133      }
     134      Mo = M / dM;
     135      dMo = sqrt (1.0 / dM);
     136    } else {
     137      // pure un-weighted
     138      for (i = ks; i < ke; i++) {
     139        M   += value[i];
     140        dM  += SQ (dvalue[i]);
     141        Nm  += 1.0; 
     142      }
     143      Mo = M / Nm;
     144      dMo = sqrt (dM) / Nm;
     145    }
    101146  }
    102147
    … …  
    111156  }
    112157  X2 = X2 / (Nm - 1);
    113   dS = sqrt (dS / Nm);
     158  dS = sqrt (dS / (Nm - 1));
    114159
    115160  stats[0].mean  = Mo;
    … …  
    122167}
    123168
     169// we could define the weight to be the only scale factor:
     170// \mu      = \sum (value_i * weight_i) / \sum (weight_i)
     171// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2)
     172// R = \sum (weight_i^2)
     173
     174// or, we could define the weight to be a scale factor times the inverse error:
     175// \mu      = \sum (value_i * weight_i / sigma_i) / \sum (weight_i)
     176// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2)
     177// R = \sum (weight_i^2)

  • trunk/Ohana/src/relphot/src/load_catalogs.c

    r33651 r33963  
    6868
    6969    // select only the brighter stars
    70     bcatalog (&catalog[i], &tcatalog);
     70    bcatalog (&catalog[i], &tcatalog, i);
    7171    dvo_catalog_unlock (&tcatalog);
    7272    dvo_catalog_free (&tcatalog);
    … …  
    123123
    124124    // ensure that the paths are absolute path names
    125     char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
     125    char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH);
    126126    free (table->hosts[i].pathname);
    127127    table->hosts[i].pathname = tmppath;
    … …  
    190190  }
    191191  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
    192     HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     192    // if any one of the remote jobs fails, we should fail and exit
     193    int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     194    if (!status) {
     195      fprintf (stderr, "at least one remote client job failed to load data, exiting\n");
     196      exit (3);
     197    }
    193198  }
    194199

  • trunk/Ohana/src/relphot/src/load_images.c

    r33651 r33963  
    1515  off_t     *LineNumber;
    1616  struct timeval start, stop;
     17  char *inSubset;
    1718
    1819  SkyTable *sky = NULL;
    … …  
    3435  MARKTIME("read image table: %f sec\n", dtime);
    3536
     37  // allocate and init an array to identify the images included in the subset
     38  ALLOCATE (inSubset, char, Nimage);
     39  memset (inSubset, 0, Nimage);
     40
    3641  // determine the populated SkyRegions overlapping the requested area
    3742
    … …  
    4550    // select the images which overlap the selected sky regions
    4651    // 'subset' points to a new copy of the data (different from 'image')
    47     subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     52    subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset, NULL);
    4853    MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    4954  } else {
    … …  
    5459      skylist = SkyListByPatch (sky, -1, &subregion);
    5560
    56       subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0;
    57       SkyList *extraList = SkyListByPatch (sky, -1, &subregion);
    58 
    5961      // select the images which overlap the selected sky regions
    6062      // 'subset' points to a new copy of the data (different from 'image')
    61       subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     63      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset, &subregion);
    6264      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    6365
    … …  
    6668      off_t *LineNumberExtra, i;
    6769
    68       subsetExtra = select_images (extraList, image, Nimage, &LineNumberExtra, &NsubsetExtra);
     70      subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0;
     71      SkyList *extraList = SkyListByPatch (sky, -1, &subregion);
     72
     73      subsetExtra = select_images (extraList, image, Nimage, inSubset, &LineNumberExtra, &NsubsetExtra, &subregion);
    6974      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    7075
    … …  
    8388      // select the images which overlap the selected sky regions
    8489      // 'subset' points to a new copy of the data (different from 'image')
    85       subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     90      subset = select_images (skylist, image, Nimage, inSubset, &LineNumber, &Nsubset, region);
    8691      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    8792    }
    8893  }
    8994
     95# if (1)
     96  FILE *ftest = fopen ("skydump.dat", "w");
     97  int i;
     98  for (i = 0; i < skylist[0].Nregions; i++) {
     99    fprintf (ftest, "%s : %f %f : %f %f\n", skylist[0].regions[i][0].name, skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, skylist[0].regions[i][0].Dmax);
     100  }
     101  fclose (ftest);
     102# endif
     103
     104  // match chips to mosaics (if applicable)
     105  initMosaics (subset, Nsubset, image, inSubset, Nimage);
     106  MARKTIME("init mosaics: %f sec\n", dtime);
     107 
    90108  // save the subset of images in the static reference in ImageOps, set up indexes
    91109  initImages (subset, LineNumber, Nsubset);
    92110  MARKTIME("init images: %f sec\n", dtime);
    93111
    94   // match chips to mosaics (if applicable)
    95   initMosaics (subset, Nsubset);
    96   MARKTIME("init mosaics: %f sec\n", dtime);
    97  
    98112  return (skylist);
    99113}

  • trunk/Ohana/src/relphot/src/reload_catalogs.c

    r33651 r33963  
    8686    TIMESTAMP(time4);
    8787
    88     setMrelFinal (&catalog, flatcorr);
     88    initMrel (&catalog, 1);
     89    setMrelFinal (&catalog, flatcorr, FALSE);
    8990    TIMESTAMP(time5);
    9091
    … …  
    104105    freeGridBins (1);
    105106    TIMESTAMP(time7);
     107
     108    if (hostID) {
     109      client_logger_message ("updated catalog file %s\n", catalog.filename);
     110    }
    106111  }
    107112
    … …  
    119124int reload_catalogs_parallel (SkyList *sky) {
    120125
    121   off_t Nimage;
    122   ImageSubset *image = getimages_subset (&Nimage);
    123 
    124   // write out the subset table of image information
     126  // name of image subset file:
    125127  char imageFile[512];
    126128  snprintf (imageFile, 512, "%s/Images.subset.dat", CATDIR);
    127   if (!ImageSubsetSave (imageFile, image, Nimage)) {
    128     fprintf (stderr, "failed to write image subset\n");
    129     exit (1);
    130   }
    131   free (image);
     129
     130  // the ApplyOffsets option re-uses an existing Images.subset.dat file
     131  if (!ApplyOffsets) {
     132    off_t Nimage;
     133    ImageSubset *image = getimages_subset (&Nimage);
     134    if (!ImageSubsetSave (imageFile, image, Nimage)) {
     135      fprintf (stderr, "failed to write image subset\n");
     136      exit (1);
     137    }
     138    free (image);
     139  } else {
     140    // load the current sky table (layout of all SkyRegions)
     141    SkyTable *skytable = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
     142    // XXX kind of hackish...
     143    ALLOCATE (sky, SkyList, 1);
     144    strcpy (sky->hosts, skytable->hosts);
     145  }
    132146
    133147  // now launch the relphot_client jobs to the parallel hosts
    … …  
    140154
    141155    // ensure that the paths are absolute path names
    142     char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
     156    char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH);
    143157    free (table->hosts[i].pathname);
    144158    table->hosts[i].pathname = tmppath;
    … …  
    201215  }
    202216  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
    203     HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     217    int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     218    if (!status) {
     219      fprintf (stderr, "at least one remote client job failed to load data, exiting\n");
     220      exit (3);
     221    }
    204222  }
    205223  return (TRUE);

  • trunk/Ohana/src/relphot/src/relphot.c

    r33651 r33963  
    2323  if (UpdateAverages) {
    2424    relphot_objects (0, NULL);
     25    exit (0);
     26  }
     27
     28  if (ApplyOffsets) {
     29    // re-run this step from a previous attempt (assumes an existing Images.subset.dat file)
     30    if (!PARALLEL) {
     31      fprintf (stderr, "-apply-offsets only makes sense in an parallel context\n");
     32      exit (2);
     33    }
     34    reload_catalogs (NULL, NULL, 0, NULL);
    2535    exit (0);
    2636  }
    … …  
    121131        plot_scatter (catalog, Ncatalog, flatcorr);
    122132      }
    123       setMrel  (catalog, Ncatalog, flatcorr);
     133      setMrel  (catalog, Ncatalog, flatcorr); // threaded
    124134      if (PLOTSTUFF) {
    125135        plot_scatter (catalog, Ncatalog, flatcorr);
    … …  
    128138      setMmos  (catalog, FALSE, flatcorr);
    129139      setMgrid (catalog, flatcorr);
     140      MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    130141   
    131142      if (PLOTSTUFF) {
    … …  
    139150      // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
    140151      // if (i % 6 == 1) rationalize_images ();
    141       if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
    142       if (i % 6 == 3) clean_stars (catalog, Ncatalog);
    143       if (i % 6 == 5) clean_mosaics ();
    144       if (i % 6 == 5) clean_images ();
     152
     153      // NOTE : in the past, I was not iterating enough before cleaning.  make sure we do
     154      // at least 8 loops first -- that should get the systematic errors down to the ~1%
     155      // level, even in cases where we have an even split between photometric data and
     156      // data with 1 mag of extinction.
     157
     158      if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
     159      if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
     160      if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
     161      if ((i > 8) && (i % 8 == 5)) clean_images ();
    145162
    146163      // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
    … …  
    148165      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
    149166      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
    150       global_stats (catalog, Ncatalog, flatcorr);
     167      if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr);
    151168      MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    152169    }

  • trunk/Ohana/src/relphot/src/relphot_client.c

    r33651 r33963  
    1919  // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc)
    2020  initialize_client (argc, argv);
     21  client_logger_init ();
    2122
    2223  // load the current sky table (layout of all SkyRegions)
    … …  
    3334      exit (2);
    3435  }
     36  client_logger_message ("loaded sky table and defined patch\n");
    3537 
    3638  switch (MODE) {
    … …  
    4244          exit (2);
    4345      }
     46      client_logger_message ("loaded catalogs\n");
     47
    4448      BrightCatalog *bcatalog = BrightCatalogMerge (catalog, Ncatalog);
     49      client_logger_message ("generated subset data\n");
     50
    4551      if (!BrightCatalogSave (BCATALOG, bcatalog)) {
    4652          fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR);
    4753          exit (2);
    4854      }
     55      client_logger_message ("generated subset table\n");
    4956      break;
    5057    }
    … …  
    5865          exit (2);
    5966      }
     67      client_logger_message ("loaded image subset data\n");
    6068     
    6169      // save the available image information in the static array in ImageOps.c
    … …  
    6674      snprintf (flatcorrFile, 1024, "%s/flatcorr.fits", CATDIR);
    6775      FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrFile, VERBOSE);
     76      client_logger_message ("loaded flat-field correction data\n");
    6877
    6978      reload_catalogs (skylist, flatcorr, HOST_ID, HOSTDIR);
     79      client_logger_message ("updated catalogs\n");
    7080      break;
    7181    }
    … …  
    7383    case MODE_UPDATE_OBJECTS: {
    7484      relphot_objects (HOST_ID, HOSTDIR);
     85      client_logger_message ("updated objects\n");
    7586      break;
    7687    }
    … …  
    8091      abort();
    8192  }
     93  client_logger_message ("done with relphot_client\n");
    8294
    8395  exit (0);
    8496}
     97

  • trunk/Ohana/src/relphot/src/relphot_objects.c

    r33651 r33963  
    7878          catalog.secfilt[j*Nsecfilt + k].Map    = NAN;
    7979          catalog.secfilt[j*Nsecfilt + k].dM     = NAN;
     80          catalog.secfilt[j*Nsecfilt + k].Mstdev = NAN_S_SHORT;
    8081          catalog.secfilt[j*Nsecfilt + k].Xm     = NAN;
    8182          catalog.secfilt[j*Nsecfilt + k].M_20   = NAN_S_SHORT;
    … …  
    8990    }
    9091
    91     setMave (&catalog, 1);
     92    populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE);
     93
     94    initMrel (&catalog, 1);
     95    setMrelFinal (&catalog, NULL, TRUE);
    9296
    9397    if (!UPDATE) {
    9498      dvo_catalog_unlock (&catalog);
     99      free_tiny_values(&catalog);
    95100      dvo_catalog_free (&catalog);
    96101      continue;
    … …  
    104109    dvo_catalog_save (&catalog, VERBOSE);
    105110    dvo_catalog_unlock (&catalog);
     111    free_tiny_values(&catalog);
    106112    dvo_catalog_free (&catalog);
    107113  }
    … …  
    127133
    128134    // ensure that the paths are absolute path names
    129     char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
     135    char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH);
    130136    free (table->hosts[i].pathname);
    131137    table->hosts[i].pathname = tmppath;

  • trunk/Ohana/src/relphot/src/select_images.c

    r33651 r33963  
    1313} SkyRegionCoords;
    1414
     15void MakeSkyCoordIndex (SkyRegionCoords **skycoords_out, double **RmaxSky_out, off_t **index_out, SkyList *skylist);
    1516void dsortindex (double *X, off_t *Y, int N);
    1617off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
    … …  
    2223  fprintf (stderr, MSG, __VA_ARGS__); }
    2324
    24 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
     25Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage, SkyRegion *region) {
    2526 
    2627  Image *image;
    … …  
    2829  off_t *line_number;
    2930  int InRange, ecode, found;
    30   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
     31  double Ri[5], Di[5], Xi[5], Yi[5];
    3132  Coords tcoords;
    32   SkyRegionCoords *skycoords;
    3333  struct timeval start, stop;
    3434 
    35   double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion;
    36 
    37   double *RmaxSky;
    38   off_t *index;
    39 
    4035  if (skylist[0].Nregions < 1) {
    4136    *Nimage = 0;
    … …  
    5449  strcpy (tcoords.ctype, "RA---TAN");
    5550
    56   ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
    57 
    58   ALLOCATE (RmaxSky, double, skylist[0].Nregions);
    59   ALLOCATE (index, off_t, skylist[0].Nregions);
    60 
    61   RminSkyRegion = +360.0;
    62   RmaxSkyRegion = -360.0;
    63   DminSkyRegion = +90.0;
    64   DmaxSkyRegion = -90.0;
    65 
    66   D_NIMAGE = 6000;
    67 
    68   // FILE *ftest = fopen ("relphot.dump.dat", "w");
    69 
    70   /* compare with each region file */
     51  double RminSkyRegion = region[0].Rmin;
     52  double RmaxSkyRegion = region[0].Rmax;
     53  double DminSkyRegion = region[0].Dmin;
     54  double DmaxSkyRegion = region[0].Dmax;
     55
     56  double RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
     57
     58  // generate a set of Rmin/Rmax values using this RmidSkyRegion for DEC bands of width
     59  // dDecBand (0.1 deg).  DecIndex = (Dec + 90) / dDecBand
     60  double *RminBand, *RmaxBand;
     61  float dDecBand = 0.1;
     62  int NDecBands = 180 / dDecBand + 1;
     63  ALLOCATE (RminBand, double, NDecBands);
     64  ALLOCATE (RmaxBand, double, NDecBands);
     65  for (i = 0; i < NDecBands; i++) {
     66    RminBand[i] = +360.0;
     67    RmaxBand[i] = -360.0;
     68  }
     69 
     70   /* compare with each region file */
    7171  for (i = 0; i < skylist[0].Nregions; i++) {
    72 
    73     /* we make positional comparisons in the projection of catalog */
    74     skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin);
    75     skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin);
    76     tcoords.crval1 = skycoords[i].Rc;
    77     tcoords.crval2 = skycoords[i].Dc;
    78 
    79     /* define catalog corners */
    80     RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords);
    81     RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords);
    82     RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords);
    83     RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords);
    84     skycoords[i].Xc[4] = skycoords[i].Xc[0];   
    85     skycoords[i].Yc[4] = skycoords[i].Yc[0];   
    86 
    87     RmaxSky[i] = skylist[0].regions[i][0].Rmax;
    88     index[i] = i;
    89 
    90     dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
    91     dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
    92     skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy;
    93     skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy;
    94     skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy;
    95     skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
    96     skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
    97 
    98     RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
    99     RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
    100     DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
    101     DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
    102   }
    103   RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
    104   MARKTIME("create sky region coords: %f sec\n", dtime);
    105 
    106   dsortindex (RmaxSky, index, skylist[0].Nregions);
    107   MARKTIME("sort sky coords: %f sec\n", dtime);
     72    int iDecBandMin = (skylist[0].regions[i][0].Dmin + 90.0) / dDecBand;
     73    int iDecBandMax = (skylist[0].regions[i][0].Dmax + 90.0) / dDecBand;
     74
     75    double RminAlt = ohana_normalize_angle_to_midpoint (skylist[0].regions[i][0].Rmin, RmidSkyRegion);
     76    double RmaxAlt = ohana_normalize_angle_to_midpoint (skylist[0].regions[i][0].Rmax, RmidSkyRegion);
     77
     78    for (j = iDecBandMin; j <= iDecBandMax; j++) {
     79      RminBand[j] = MIN(RminBand[j], RminAlt);
     80      RmaxBand[j] = MAX(RmaxBand[j], RmaxAlt);
     81    }
     82  }
     83
     84# if (1)
     85  // XXX quick test of the dec bands:
     86  FILE *fout = fopen ("dec.bands.dat", "w");
     87  for (i = 0; i < NDecBands; i++) {
     88    fprintf (fout, "%d %f : %f - %f\n", (int) i, dDecBand * i - 90.0, RminBand[i], RmaxBand[i]);
     89  }
     90  fclose (fout);
     91# endif
     92
     93  SkyRegionCoords *skycoords = NULL;
     94  double *RmaxSky = NULL;
     95  off_t *index = NULL;
     96
     97  if (!USE_BASIC_CHECK) {
     98    MakeSkyCoordIndex (&skycoords, &RmaxSky, &index, skylist);
     99  }
    108100
    109101  if (VERBOSE) fprintf (stderr, "finding images\n");
    110102  BuildChipMatch (timage, Ntimage);
    111   MARKTIME("build chip match: %f sec\n", dtime);
     103  MARKTIME("build chip match for %d images: %f sec\n", (int) Ntimage, dtime);
     104
     105  D_NIMAGE = 6000;
    112106
    113107  nimage = 0;
    … …  
    185179    if (RminImage > RmaxSkyRegion) continue;
    186180    if (RmaxImage < RminSkyRegion) continue;
     181
     182    // the above checks are only valid for the outermost region.  however, at a given
     183    // declination, the range of RA is limited by the actual catalog boundaries
     184    // check if this image is in range for the Dmax and Dmin locations
     185
     186    int iDecBandMin = (DminImage + 90.0) / dDecBand;
     187    int iDecBandMax = (DmaxImage + 90.0) / dDecBand;
     188
     189    // the sky region RA is defined to be 0 - 360.0
     190    if (RminImage < RminBand[iDecBandMin]) {
     191      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     192      continue;
     193    }
     194    if (RminImage < RminBand[iDecBandMax]) {
     195      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     196      continue;
     197    }
     198    if (RmaxImage > RmaxBand[iDecBandMin]) {
     199      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     200      continue;
     201    }
     202    if (RmaxImage > RmaxBand[iDecBandMax]) {
     203      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     204      continue;
     205    }
    187206
    188207    // image overlaps region, keep it
    … …  
    230249  found_it:
    231250    image[nimage] = timage[i];
     251    inSubset[i] = TRUE;
    232252    /* always allow 'few' images to succeed, if possible (new images / detections may have
    233253     * been added) */
    … …  
    255275  MARKTIME("finish image selection: %f sec\n", dtime);
    256276
    257   // fclose (ftest);
    258 
    259277  if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);
    260278
    261279  REALLOCATE (image, Image, MAX (nimage, 1));
    262280  REALLOCATE (line_number, off_t, MAX (nimage, 1));
    263   free (skycoords);
    264   free (RmaxSky);
    265   free (index);
     281
     282  if (!USE_BASIC_CHECK) {
     283    free (skycoords);
     284    free (RmaxSky);
     285    free (index);
     286  }
    266287
    267288  *Nimage  = nimage;
    … …  
    379400  return (Nlo);
    380401}
     402
     403void MakeSkyCoordIndex (SkyRegionCoords **skycoords_out, double **RmaxSky_out, off_t **index_out, SkyList *skylist) {
     404
     405  off_t i;
     406  SkyRegionCoords *skycoords;
     407  Coords tcoords;
     408  double *RmaxSky;
     409  off_t *index;
     410  double dx, dy;
     411
     412  // the comparison is made in the catalog local projection. below we set crval1,2
     413  tcoords.crpix1 = tcoords.crpix2 = 0.0;
     414  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
     415  tcoords.pc1_1  = tcoords.pc2_2 = 1.0;
     416  tcoords.pc1_2  = tcoords.pc2_1 = 0.0;
     417  strcpy (tcoords.ctype, "RA---TAN");
     418
     419  /* compare with each region file */
     420  ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
     421
     422  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
     423  ALLOCATE (index, off_t, skylist[0].Nregions);
     424
     425  for (i = 0; i < skylist[0].Nregions; i++) {
     426
     427    /* we make positional comparisons in the projection of catalog */
     428    skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin);
     429    skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin);
     430    tcoords.crval1 = skycoords[i].Rc;
     431    tcoords.crval2 = skycoords[i].Dc;
     432
     433    /* define catalog corners */
     434    RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords);
     435    RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords);
     436    RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords);
     437    RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords);
     438    skycoords[i].Xc[4] = skycoords[i].Xc[0];   
     439    skycoords[i].Yc[4] = skycoords[i].Yc[0];   
     440
     441    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
     442    index[i] = i;
     443
     444    dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
     445    dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
     446    skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy;
     447    skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy;
     448    skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy;
     449    skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
     450    skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
     451  }
     452
     453  dsortindex (RmaxSky, index, skylist[0].Nregions);
     454
     455  *RmaxSky_out = RmaxSky;
     456  *index_out = index;
     457  *skycoords_out = skycoords;
     458
     459  return;
     460}

  • trunk/Ohana/src/relphot/src/setMrelFinal.c

    r33651 r33963  
    55// output dbFlags values
    66
    7 void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) {
     7void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) {
    88
    99  off_t i;
    … …  
    4343        catalog[0].secfilt[N].Map    = NAN;
    4444        catalog[0].secfilt[N].dM     = NAN;
     45        catalog[0].secfilt[N].Mstdev = NAN_S_SHORT;
    4546        catalog[0].secfilt[N].Xm     = NAN;
    4647        catalog[0].secfilt[N].M_20   = NAN_S_SHORT;
    … …  
    6566          }
    6667       
    67           if (!KEEP_UBERCAL) {
    68             catalog[0].measure[m].Mcal = 0;
    69             catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     68          // only reset Mcal for measures with a matching image
     69          // do not reset Mcal for ubercal images unless explicitly requested
     70
     71          if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     72            if (!KEEP_UBERCAL) {
     73              catalog[0].measure[m].Mcal = 0.0;
     74              catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     75            }
     76          } else {
     77            if (getImageEntry (m, 0) >= 0) {
     78              catalog[0].measure[m].Mcal = 0.0;
     79            }
    7080          }
    7181          catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS;
    … …  
    8494  }
    8595
    86   clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
     96  // XXX make this optional? (do not clean for -averages?)
     97  if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
    8798  for (i = 0; i < 5; i++) {
    8899    skip_measurements (catalog, i, flatcorr);       /* set ID_MEAS_SKIP for measures to be skipped */
    89100    setMrelOutput  (catalog, 1, i, flatcorr);    /* set Mrel using remaining measures */
    90101  }
    91   setMcalOutput (catalog, 1, flatcorr);
     102  if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr);
    92103
    93104  /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */
    … …  
    248259    }
    249260  }
    250   if (VERBOSE2) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
     261  if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
    251262                        pass, Ntot, Ntry, Nskip, Nkeep);
    252263}
Note: See TracChangeset for help on using the changeset viewer.