IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31668


Ignore:
Timestamp:
Jun 22, 2011, 12:43:18 AM (15 years ago)
Author:
eugene
Message:

merged from eam/ipp20110505: various memory and speed improvements; fixed 0 cloud issue

Location:
trunk/Ohana/src/relphot
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/doc/allsky.txt

    r31450 r31668  
     1
     22011.05.10
     3
     4  Looking into them memory usage of relphot:
     5  AverageTiny = 32 bytes (including 8 byte padding)
     6  MeasureTiny = 56 bytes (including 8 byte padding)
     7  Image Indexes:
     8    bin   : Nmeasure * off_t
     9    clist : Nmeasure * off_t
     10    mlist : Nmeasure * off_t
     11  Mosaic Indexes:
     12    bin   : Nmeasure * off_t
     13    clist : Nmeasure * off_t
     14    mlist : Nmeasure * off_t
     15  Total:
     16    32 bytes * Naverage (values)
     17    56 bytes * Nmeasure (values)
     18    64 bytes * Nmeasure (indexes)
     19
     20  In my test run, I have 3.3M average + 37.6M measure:
     21    2.1GB values
     22    2.2GB indexes
     23
     24  This matches the memory usage graph.  But, after the indexes have
     25  been loaded, the memory footprint grows (slowly) by another 1GB as
     26  the processing runs along.
     27
     28
     29using 23153493 of 640663697 stars (192585953 of 289672671 measurements)
     30using  3306896 of  48399614 stars ( 37582679 of 193863612 measurements)
     31
     32
    133
    2342011.04.12

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

    r31639 r31668  
    6363int    RESET;
    6464int    UPDATE;
     65int    SAVE_IMAGE_UPDATES;
    6566int    PLOTSTUFF;
    6667int    SAVEPLOT;
    … …  
    139140void          dump_grid           PROTO((void));
    140141int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
    141 void          findImages          PROTO((Catalog *catalog, int Ncatalog));
    142 int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
     142void          findImages          PROTO((Catalog *catalog, int Ncatalog, int doImageList));
     143int           findMosaics         PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
    143144
    144145void set_db (FITS_DB *in);
    … …  
    149150
    150151void          freeGridBins        PROTO((int Ncatalog));
    151 void          freeImageBins       PROTO((int Ncatalog));
    152 void          freeMosaicBins      PROTO((int Ncatalog));
     152void          freeImageBins       PROTO((int Ncatalog, int doImageList));
     153void          freeMosaicBins      PROTO((int Ncatalog, int doMosaicList));
    153154void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
    154155int           gcatalog            PROTO((Catalog *catalog, int FINAL));
    … …  
    160161float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
    161162Image        *getimage            PROTO((off_t N));
    162 Image        *getimages           PROTO((off_t *N));
     163Image        *getimages           PROTO((off_t *N, off_t **LineNumber));
    163164void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
    164165void          initGrid            PROTO((int dX, int dY));
    165166void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
    166 void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
    167 void          initImages          PROTO((Image *input, off_t N));
    168 void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
     167void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int doImageList));
     168void          initImages          PROTO((Image *input, off_t *LineNumber, off_t N));
     169void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
    169170void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
    170171void          initMosaics         PROTO((Image *image, off_t Nimage));
    … …  
    179180int           main                PROTO((int argc, char **argv));
    180181void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
    181 void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat));
    182 void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat));
     182void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat, int doImageList));
     183void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat, int doMosaicList));
    183184double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
    184185void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));

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

    r31450 r31668  
    99static Image        *image;   // list of available images 
    1010static off_t        Nimage;   // number of available images
     11static off_t       *LineNumber; // match of subset to full image table
    1112
    1213// to search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted index
    … …  
    1415static off_t        *imageIdx; // list of index for image IDs
    1516
    16 Image *getimages (off_t *N) {
     17Image *getimages (off_t *N, off_t **line_number) {
    1718
    1819  *N = Nimage;
     20  if (line_number) *line_number = LineNumber;
    1921  return (image);
    2022}
    … …  
    2426}
    2527
    26 void initImages (Image *input, off_t N) {
     28void initImages (Image *input, off_t *line_number, off_t N) {
    2729
    2830  off_t i;
    2931
    3032  image = input;
     33  LineNumber = line_number;
    3134  Nimage = N;
    3235
    … …  
    6669}
    6770
    68 void initImageBins (Catalog *catalog, int Ncatalog) {
     71void initImageBins (Catalog *catalog, int Ncatalog, int doImageList) {
    6972
    7073  off_t i, j;
    … …  
    7679  }
    7780
    78   ALLOCATE (Nlist, off_t, Nimage);
    79   ALLOCATE (NLIST, off_t, Nimage);
    80   ALLOCATE (clist, off_t *, Nimage);
    81   ALLOCATE (mlist, off_t *, Nimage);
    82 
    83   for (i = 0; i < Nimage; i++) {
    84     Nlist[i] = 0;
    85     NLIST[i] = 100;
    86     ALLOCATE (clist[i], off_t, NLIST[i]);
    87     ALLOCATE (mlist[i], off_t, NLIST[i]);
    88   }
    89 }
    90 
    91 void freeImageBins (int Ncatalog) {
     81  if (doImageList) {
     82    ALLOCATE (Nlist, off_t, Nimage);
     83    ALLOCATE (NLIST, off_t, Nimage);
     84    ALLOCATE (clist, off_t *, Nimage);
     85    ALLOCATE (mlist, off_t *, Nimage);
     86
     87    for (i = 0; i < Nimage; i++) {
     88      Nlist[i] = 0;
     89      NLIST[i] = 100;
     90      ALLOCATE (clist[i], off_t, NLIST[i]);
     91      ALLOCATE (mlist[i], off_t, NLIST[i]);
     92    }
     93  }
     94}
     95
     96void freeImageBins (int Ncatalog, int doImageList) {
    9297
    9398  off_t i;
    … …  
    97102  }
    98103  free (bin);
    99   for (i = 0; i < Nimage; i++) {
    100     free (clist[i]);
    101     free (mlist[i]);
    102   }
    103   free (clist);
    104   free (mlist);
    105   free (Nlist);
    106   free (NLIST);
     104
     105  if (doImageList) {
     106    for (i = 0; i < Nimage; i++) {
     107      free (clist[i]);
     108      free (mlist[i]);
     109    }
     110    free (clist);
     111    free (mlist);
     112    free (Nlist);
     113    free (NLIST);
     114  }
    107115}
    108116
    109117/* select all image equivalent to the active photcode set */
    110 void findImages (Catalog *catalog, int Ncatalog) {
     118void findImages (Catalog *catalog, int Ncatalog, int doImageList) {
    111119
    112120  off_t j;
    … …  
    122130      }
    123131      if (!found) continue;
    124       matchImage (catalog, j, i);
     132      matchImage (catalog, j, i, doImageList);
    125133      Nmatch ++;
    126134    }
    … …  
    178186}
    179187
    180 void matchImage (Catalog *catalog, off_t meas, int cat) {
     188void matchImage (Catalog *catalog, off_t meas, int cat, int doImageList) {
    181189
    182190  off_t idx, ID;
    … …  
    204212  bin[cat][meas] = idx;
    205213
    206   // index for image, Nentry -> catalog
    207   clist[idx][Nlist[idx]] = cat;
    208 
    209   // index for image, Nentry -> measure
    210   mlist[idx][Nlist[idx]] = meas;
    211   Nlist[idx] ++;
    212 
    213   if (Nlist[idx] == NLIST[idx]) {
    214     NLIST[idx] += 100;
    215     REALLOCATE (clist[idx], off_t, NLIST[idx]);
    216     REALLOCATE (mlist[idx], off_t, NLIST[idx]);
    217   }     
     214  if (doImageList) {
     215    // index for image, Nentry -> catalog
     216    clist[idx][Nlist[idx]] = cat;
     217
     218    // index for image, Nentry -> measure
     219    mlist[idx][Nlist[idx]] = meas;
     220    Nlist[idx] ++;
     221
     222    if (Nlist[idx] == NLIST[idx]) {
     223      NLIST[idx] += 100;
     224      REALLOCATE (clist[idx], off_t, NLIST[idx]);
     225      REALLOCATE (mlist[idx], off_t, NLIST[idx]);
     226    }   
     227  }
    218228
    219229  return;
    … …  
    343353   
    344354    liststats (list, dlist, N, &stats);
     355
     356    float CLOUD_TOLERANCE = 0.01;
    345357    image[i].Mcal  = stats.mean;
    346358    image[i].dMcal = stats.error;
    … …  
    348360    image[i].nFitPhotom = N;
    349361    image[i].Xm    = 100.0*log10(stats.chisq);
     362
     363    if (image[i].Mcal < -CLOUD_TOLERANCE) {
     364      image[i].Mcal = 0.0;
     365    }
    350366  }
    351367  free (list);

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

    r31450 r31668  
    195195  if (!MOSAIC_ZEROPT) return;
    196196
    197   image = getimages (&Nimage);
    198 
     197  image = getimages (&Nimage, NULL);
     198
     199  // copy the mosaic results to the images.  set the mosaic Mcal to 0.0 since we have moved its
     200  // impact to the images
    199201  for (i = 0; i < Nmosaic; i++) {
    200202    for (j = 0; j < Nimlist[i]; j++) {
    201203      im = imlist[i][j];
    202       image[im].Mcal = mosaic[i].Mcal;
     204      image[im].Mcal += mosaic[i].Mcal;
    203205      image[im].dMcal = mosaic[i].dMcal;
    204206      image[im].Xm = mosaic[i].Xm;
    … …  
    208210      image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR);
    209211    }
     212    mosaic[i].Mcal = 0.0;
    210213  }     
    211214}
    212215
    213 void initMosaicBins (Catalog *catalog, int Ncatalog) {
     216void initMosaicBins (Catalog *catalog, int Ncatalog, int doMosaicList) {
    214217
    215218  off_t i, j;
    … …  
    224227  }
    225228
    226   /* mosaic -> measure */
    227   ALLOCATE (Nlist, off_t,   Nmosaic);
    228   ALLOCATE (NLIST, off_t,   Nmosaic);
    229   ALLOCATE (clist, int *,   Nmosaic);
    230   ALLOCATE (mlist, off_t *, Nmosaic);
    231 
    232   for (i = 0; i < Nmosaic; i++) {
    233     Nlist[i] = 0;
    234     NLIST[i] = 100;
    235     ALLOCATE (clist[i], int,   NLIST[i]);
    236     ALLOCATE (mlist[i], off_t, NLIST[i]);
    237   }
    238 }
    239 
    240 void freeMosaicBins (int Ncatalog) {
     229  if (doMosaicList) {
     230    /* mosaic -> measure */
     231    ALLOCATE (Nlist, off_t,   Nmosaic);
     232    ALLOCATE (NLIST, off_t,   Nmosaic);
     233    ALLOCATE (clist, int *,   Nmosaic);
     234    ALLOCATE (mlist, off_t *, Nmosaic);
     235
     236    for (i = 0; i < Nmosaic; i++) {
     237      Nlist[i] = 0;
     238      NLIST[i] = 100;
     239      ALLOCATE (clist[i], int,   NLIST[i]);
     240      ALLOCATE (mlist[i], off_t, NLIST[i]);
     241    }
     242  }
     243}
     244
     245void freeMosaicBins (int Ncatalog, int doMosaicList) {
    241246
    242247  off_t i;
    … …  
    250255  free (bin);
    251256
    252   /* mosaic -> measure */
    253   for (i = 0; i < Nmosaic; i++) {
    254     free (clist[i]);
    255     free (mlist[i]);
    256   }
    257   free (Nlist);
    258   free (NLIST);
    259   free (clist);
    260   free (mlist);
    261 }
    262 
    263 int findMosaics (Catalog *catalog, int Ncatalog) {
     257  if (doMosaicList) {
     258    /* mosaic -> measure */
     259    for (i = 0; i < Nmosaic; i++) {
     260      free (clist[i]);
     261      free (mlist[i]);
     262    }
     263    free (Nlist);
     264    free (NLIST);
     265    free (clist);
     266    free (mlist);
     267  }
     268}
     269
     270int findMosaics (Catalog *catalog, int Ncatalog, int doMosaicList) {
    264271 
    265272  int i, ecode, found, Ns;
    … …  
    281288      }
    282289      if (!found) continue;
    283       matchMosaics (catalog, j, i);
     290      matchMosaics (catalog, j, i, doMosaicList);
    284291      Nmatch ++;
    285292    }
    … …  
    289296}
    290297
    291 void matchMosaics (Catalog *catalog, off_t meas, int cat) {
     298void matchMosaics (Catalog *catalog, off_t meas, int cat, int doMosaicList) {
    292299
    293300  off_t idx, ID, mosID;
    … …  
    322329  bin[cat][meas] = mosID;
    323330
    324   clist[mosID][Nlist[mosID]] = cat;
    325   mlist[mosID][Nlist[mosID]] = meas;
    326   Nlist[mosID] ++;
     331  if (doMosaicList) {
     332    clist[mosID][Nlist[mosID]] = cat;
     333    mlist[mosID][Nlist[mosID]] = meas;
     334    Nlist[mosID] ++;
    327335   
    328   if (Nlist[mosID] == NLIST[mosID]) {
    329     NLIST[mosID] += 100;
    330     REALLOCATE (clist[mosID], int,   NLIST[mosID]);
    331     REALLOCATE (mlist[mosID], off_t, NLIST[mosID]);
    332   }     
     336    if (Nlist[mosID] == NLIST[mosID]) {
     337      NLIST[mosID] += 100;
     338      REALLOCATE (clist[mosID], int,   NLIST[mosID]);
     339      REALLOCATE (mlist[mosID], off_t, NLIST[mosID]);
     340    }   
     341  }
    333342  return;
    334343}
    … …  
    361370  if (FREEZE_MOSAICS) return (FALSE);
    362371
    363   image = getimages (&N);
     372  image = getimages (&N, NULL);
    364373
    365374  Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    442451    liststats (list, dlist, N, &stats);
    443452    if (PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
     453
     454    float CLOUD_TOLERANCE = 0.01;
    444455    mosaic[i].Mcal  = stats.mean;
    445456    mosaic[i].dMcal = stats.error;
    … …  
    447458    mosaic[i].nFitPhotom = N;
    448459    mosaic[i].Xm    = 100.0*log10(stats.chisq);
     460
     461    if (mosaic[i].Mcal < -CLOUD_TOLERANCE) {
     462      mosaic[i].Mcal = 0.0;
     463    }
    449464  }
    450465  free (list);
    … …  
    482497  if (FREEZE_MOSAICS) return (FALSE);
    483498
    484   image = getimages (&Nimage);
     499  image = getimages (&Nimage, NULL);
    485500
    486501  ALLOCATE (imageOffset, double, Nmosaic);
    487 
    488502  ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic
    489503  ALLOCATE (NSlist, int,   Nmosaic); // number of stars on mosaic
    … …  
    621635    float dM = 0.0;
    622636    for (j = 0; j < NSlist[i]; j++) {
    623       Nstars = Slist[i][j];
    624       if (starNcount[Nstars][Ns] > 1) {
    625         dM += (starOffset[Nstars][Ns] / starNcount[Nstars][Ns]);
     637      k = Slist[i][j];
     638      if (starNcount[k][Ns] > 1) {
     639        dM += (starOffset[k][Ns] / starNcount[k][Ns]);
    626640      }
    627641    }
    … …  
    647661  free (starNcount);
    648662
     663  free (seclist);
     664  free (NSlist);
     665  free (NSLIST);
    649666  for (i = 0; i < Nmosaic; i++){
    650667    free (Slist[i]);
    651668  }
    652   free (NSlist);
    653   free (NSLIST);
     669  free (Slist);
    654670  free (imageOffset);
    655   free (seclist);
    656671  return (TRUE);
    657672}
    … …  
    853868  if (!MOSAIC_ZEROPT) return;
    854869
    855   image = getimages (&Nimage);
     870  image = getimages (&Nimage, NULL);
    856871
    857872  N = 0;

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

    r31450 r31668  
    587587  initstats (STATMODE);
    588588  if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
     589  free (list);
     590  free (dlist);
    589591  free (ilist);
    590592  free (tlist);

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

    r31450 r31668  
    128128  }
    129129
     130  SAVE_IMAGE_UPDATES = TRUE;
     131  if ((N = get_argument (argc, argv, "-skip-image-updates"))) {
     132    remove_argument (N, &argc, argv);
     133    SAVE_IMAGE_UPDATES = FALSE;
     134  }
     135
    130136  MaxDensityUse = FALSE;
    131137  if ((N = get_argument (argc, argv, "-max-density"))) {

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

    r31450 r31668  
    55int CopyAverageTiny (AverageTiny *averageT, Average *average) {
    66
    7   averageT[0].R     = average[0].R;
    8   averageT[0].D     = average[0].D;
    9   averageT[0].flags = average[0].flags;
     7  averageT[0].R             = average[0].R;
     8  averageT[0].D             = average[0].D;
     9  averageT[0].flags         = average[0].flags;
    1010  averageT[0].Nmeasure      = average[0].Nmeasure;
    1111  averageT[0].measureOffset = average[0].measureOffset;
    … …  
    293293  off_t i;
    294294
     295  AverageTiny *averageT;
     296  Average *average;
     297
     298  MeasureTiny *measureT;
     299  Measure *measure;
     300
    295301  ALLOCATE (catalog[0].measureT, MeasureTiny, catalog[0].Nmeasure);
    296302  ALLOCATE (catalog[0].averageT, AverageTiny, catalog[0].Naverage);
    297303
     304  average  = catalog[0].average;
     305  averageT = catalog[0].averageT;
     306
     307  measure  = catalog[0].measure;
     308  measureT = catalog[0].measureT;
     309
    298310  for (i = 0; i < catalog[0].Naverage; i++) {
    299     CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]);
     311    // CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]);
     312    averageT[i].R             = average[i].R;
     313    averageT[i].D             = average[i].D;
     314    averageT[i].flags         = average[i].flags;
     315    averageT[i].Nmeasure      = average[i].Nmeasure;
     316    averageT[i].measureOffset = average[i].measureOffset;
    300317  }
    301318
    302319  for (i = 0; i < catalog[0].Nmeasure; i++) {
    303     CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]);
     320    // CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]);
     321    measureT[i].dR       = measure[i].dR;
     322    measureT[i].dD       = measure[i].dD;
     323    measureT[i].M        = measure[i].M;
     324    measureT[i].Mcal     = measure[i].Mcal;
     325    measureT[i].dM       = measure[i].dM;
     326    measureT[i].airmass  = measure[i].airmass;
     327    measureT[i].Xccd     = measure[i].Xccd;
     328    measureT[i].Yccd     = measure[i].Yccd;
     329    measureT[i].t        = measure[i].t;
     330    measureT[i].dt       = measure[i].dt;
     331    measureT[i].averef   = measure[i].averef;
     332    measureT[i].imageID  = measure[i].imageID;
     333    measureT[i].dbFlags  = measure[i].dbFlags;
     334    measureT[i].photcode = measure[i].photcode;
    304335  }
    305336
    … …  
    309340int free_tiny_values (Catalog *catalog) {
    310341
    311   free (catalog[0].averageT);
    312   free (catalog[0].measureT);
    313   return (TRUE);
    314 }
    315 
     342  if (catalog[0].averageT) free (catalog[0].averageT);
     343  if (catalog[0].measureT) free (catalog[0].measureT);
     344  return (TRUE);
     345}
     346

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

    r30616 r31668  
    33Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) {
    44
    5   int i, Nmeas, Nstar, Nmeas_total, Nstar_total;
     5  off_t i, Nmeas, Nstar, Nmeas_total, Nstar_total;
    66  Catalog *catalog, tcatalog;
    77
    … …  
    1414  // load data from each region file, only use bright stars
    1515  for (i = 0; i < skylist[0].Nregions; i++) {
     16    dvo_catalog_init (&catalog[i], TRUE);
    1617
    1718    // set up the basic catalog info
    … …  
    4950  }
    5051  if (Nstar < 2) {
    51     fprintf (stderr, "warning: insufficient stars %d\n", Nstar);
     52    fprintf (stderr, "warning: insufficient stars "OFF_T_FMT"\n", Nstar);
    5253  }
    5354
    54   fprintf (stderr, "using %d of %d stars (%d of %d measurements)\n", Nstar, Nstar_total, Nmeas, Nmeas_total);
     55  fprintf (stderr, "using "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measurements)\n", Nstar, Nstar_total, Nmeas, Nmeas_total);
    5556  if (Nstar < 1) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n");
    5657   

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

    r31450 r31668  
    2525  SkyTableSetFilenames (sky, CATDIR, "cpt");
    2626 
    27   // determine the populated SkyRegions overlapping the requested area
    28   if (RegionSelect) {
    29     skylist = SkyListByPatch (sky, -1, region);
    30   } else {
    31     Nchar = strlen(regionName);
    32     if (!strcmp (&regionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0;
    33     skylist = SkyListByName (sky, regionName);
    34   }
    35 
    3627  // convert database table to internal structure (binary to Image)
     28  // 'image' points to the same memory as db->ftable->buffer
    3729  image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
    3830  if (!image) {
    … …  
    4234  MARKTIME("read image table: %f sec\n", dtime);
    4335
    44   // select the images which overlap the selected sky regions
    45   subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
    46   MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
     36  // determine the populated SkyRegions overlapping the requested area
     37  if (RegionSelect) {
     38    if (region[0].Rmin > region[0].Rmax) {
     39      SkyRegion subregion;
     40      subregion = *region;
     41      subregion.Rmin = 0.0; subregion.Rmax = region[0].Rmax;
     42      skylist = SkyListByPatch (sky, -1, &subregion);
     43
     44      subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0;
     45      SkyList *extraList = SkyListByPatch (sky, -1, &subregion);
     46
     47      // select the images which overlap the selected sky regions
     48      // 'subset' points to a new copy of the data (different from 'image')
     49      subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     50      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
     51
     52      Image *subsetExtra;
     53      off_t NsubsetExtra;
     54      off_t *LineNumberExtra, i;
     55
     56      subsetExtra = select_images (extraList, image, Nimage, &LineNumberExtra, &NsubsetExtra);
     57      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
     58
     59      REALLOCATE (subset, Image, MAX (Nsubset + NsubsetExtra, 1));
     60      REALLOCATE (LineNumber, off_t, MAX (Nsubset + NsubsetExtra, 1));
     61     
     62      for (i = 0; i < NsubsetExtra; i++) {
     63        subset[i+Nsubset] = subsetExtra[i];
     64        LineNumber[i+Nsubset] = LineNumberExtra[i];
     65      }
     66      Nsubset += NsubsetExtra;
     67      MARKTIME("selected %d overlapping images: %f sec\n", (int) NsubsetExtra, dtime);
     68    } else {
     69      skylist = SkyListByPatch (sky, -1, region);
     70
     71      // select the images which overlap the selected sky regions
     72      // 'subset' points to a new copy of the data (different from 'image')
     73      subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     74      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
     75    }
     76  } else {
     77    Nchar = strlen(regionName);
     78    if (!strcmp (&regionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0;
     79    skylist = SkyListByName (sky, regionName);
     80
     81    // select the images which overlap the selected sky regions
     82    // 'subset' points to a new copy of the data (different from 'image')
     83    subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     84    MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
     85  }
     86
    4787
    4888  // generate db->vtable from db->ftable based on the selection
    4989  // XXX does this simply duplicate the memory needlessly?  we recreate these lines
    5090  // in reload_images.  If we had saved the line numbers, we could avoid this
    51   gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
    52   MARKTIME("converted ftable to vtable: %f sec\n", dtime);
     91  // vtable points *another* copy of the subset rows
     92  // (the later call to 'reload_images' copies the subset elements back on top of
     93  // the rows of the vtable)
     94  // gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
     95  // MARKTIME("converted ftable to vtable: %f sec\n", dtime);
    5396
    5497  // save the subset of images in the static reference in ImageOps, set up indexes
    55   initImages (subset, Nsubset);
     98  initImages (subset, LineNumber, Nsubset);
    5699  MARKTIME("init images: %f sec\n", dtime);
    57100
    … …  
    66109
    67110  Image     *image;
    68   off_t     Nimage, Nx, i;
     111  off_t     Nimage, Nx, i, *LineNumber;
    69112  VTable    *vtable;
    70113
    71   image = getimages (&Nimage);
     114  image = getimages (&Nimage, &LineNumber);
    72115
     116  gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nimage);
    73117  vtable = &db[0].vtable;
    74118

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

    r31450 r31668  
    1818  double time1 = 0.0;
    1919  double time2 = 0.0;
    20   double time3 = 0.0;
     20  double time3a = 0.0;
     21  double time3b = 0.0;
    2122  double time4 = 0.0;
    2223  double time5 = 0.0;
    … …  
    5758
    5859    populate_tiny_values(&catalog);
     60    TIMESTAMP(time3a);
    5961
    60     initImageBins  (&catalog, 1);
    61     initMosaicBins (&catalog, 1);
     62    initImageBins  (&catalog, 1, FALSE);
     63    initMosaicBins (&catalog, 1, FALSE);
    6264    initGridBins   (&catalog, 1);
    63     TIMESTAMP(time3);
     65    TIMESTAMP(time3b);
    6466
    65     findImages (&catalog, 1);  // FX
    66     findMosaics (&catalog, 1); //
     67    findImages (&catalog, 1, FALSE);
     68    findMosaics (&catalog, 1, FALSE);
    6769    TIMESTAMP(time4);
    6870
    … …  
    7779    TIMESTAMP(time6);
    7880
    79     freeImageBins (1);
    80     freeMosaicBins (1);
     81    freeImageBins (1, FALSE);
     82    freeMosaicBins (1, FALSE);
    8183    freeGridBins (1);
    8284    TIMESTAMP(time7);
    8385  }
    8486
    85   fprintf (stderr, "time1 %f : find catalog\n", time1);
    86   fprintf (stderr, "time2 %f : load catalog\n", time2);
    87   fprintf (stderr, "time3 %f : init imbins\n",  time3);
    88   fprintf (stderr, "time4 %f : find images\n",  time4);
    89   fprintf (stderr, "time5 %f : set Mrel\n",     time5);
    90   fprintf (stderr, "time6 %f : save catalog\n", time6);
    91   fprintf (stderr, "time7 %f : free catalog\n", time7);
     87  fprintf (stderr, "time1  %f : find catalog\n", time1);
     88  fprintf (stderr, "time2  %f : load catalog\n", time2);
     89  fprintf (stderr, "time3a %f : init imbins\n",  time3a);
     90  fprintf (stderr, "time3b %f : init imbins\n",  time3b);
     91  fprintf (stderr, "time4  %f : find images\n",  time4);
     92  fprintf (stderr, "time5  %f : set Mrel\n",     time5);
     93  fprintf (stderr, "time6  %f : save catalog\n", time6);
     94  fprintf (stderr, "time7  %f : free catalog\n", time7);
    9295}

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

    r31450 r31668  
    11# include "relphot.h"
    22# define USE_DIRECT 0
     3
     4# define TIMESTAMP(TIME) \
     5    gettimeofday (&stop, (void *) NULL);        \
     6    dtime = DTIME (stop, start);                \
     7    TIME += dtime;                              \
     8    gettimeofday (&start, (void *) NULL);
    39
    410# define MARKTIME(MSG,...) { \
    … …  
    1723
    1824  gettimeofday (&start, (void *) NULL);
     25
     26  // XXX quick and stupid test:
     27  if (0) {
     28    int i, j;
     29    off_t *Nlist, *NLIST, **mlist;
     30    off_t *NlistI, *NLISTI, **mlistI;
     31    int **clist;
     32    int **clistI;
     33    float dtime;
     34    float time1 = 0.0;
     35    float time2 = 0.0;
     36
     37    int Nmosaic = 10000;
     38    int Nimage = 600000;
     39
     40    for (j = 0; j < 19000; j++) {
     41
     42      // mosaic indexes
     43      ALLOCATE (Nlist, off_t,   Nmosaic);
     44      ALLOCATE (NLIST, off_t,   Nmosaic);
     45      ALLOCATE (clist, int *,   Nmosaic);
     46      ALLOCATE (mlist, off_t *, Nmosaic);
     47   
     48      for (i = 0; i < Nmosaic; i++) {
     49        Nlist[i] = 0;
     50        NLIST[i] = 100;
     51        ALLOCATE (clist[i], int,   NLIST[i]);
     52        ALLOCATE (mlist[i], off_t, NLIST[i]);
     53      }
     54
     55      // image indexes
     56      ALLOCATE (NlistI, off_t,   Nimage);
     57      ALLOCATE (NLISTI, off_t,   Nimage);
     58      ALLOCATE (clistI, int *,   Nimage);
     59      ALLOCATE (mlistI, off_t *, Nimage);
     60   
     61      for (i = 0; i < Nimage; i++) {
     62        NlistI[i] = 0;
     63        NLISTI[i] = 100;
     64        ALLOCATE (clistI[i], int,   NLISTI[i]);
     65        ALLOCATE (mlistI[i], off_t, NLISTI[i]);
     66      }
     67
     68      TIMESTAMP(time1);
     69
     70      // free mosaic indexes
     71      for (i = 0; i < Nmosaic; i++) {
     72        free (clist[i]);
     73        free (mlist[i]);
     74      }
     75      free (Nlist);
     76      free (NLIST);
     77      free (clist);
     78      free (mlist);
     79
     80      // free image indexes
     81      for (i = 0; i < Nimage; i++) {
     82        free (clistI[i]);
     83        free (mlistI[i]);
     84      }
     85      free (NlistI);
     86      free (NLISTI);
     87      free (clistI);
     88      free (mlistI);
     89
     90      TIMESTAMP(time2);
     91
     92      if (j % 100 == 0) {
     93        fprintf (stderr, "done with %d\n", j);
     94        fprintf (stderr, "time1  %f : init mosaic\n", time1);
     95        fprintf (stderr, "time2  %f : free mosaic\n", time2);
     96      }
     97    }
     98   
     99    fprintf (stderr, "time1  %f : init mosaic\n", time1);
     100    fprintf (stderr, "time2  %f : free mosaic\n", time2);
     101    exit (1);
     102  }
    19103
    20104  /* get configuration info, args */
    … …  
    66150
    67151  /* match measurements with images, mosaics */
    68   initImageBins  (catalog, Ncatalog);
     152  initImageBins  (catalog, Ncatalog, TRUE);
    69153  MARKTIME("-- make image bins: %f sec\n", dtime);
    70154
    71   initMosaicBins (catalog, Ncatalog);
     155  initMosaicBins (catalog, Ncatalog, TRUE);
    72156  initGridBins   (catalog, Ncatalog);
    73157  initMrel (catalog, Ncatalog);
    74158
    75   findImages (catalog, Ncatalog);
     159  findImages (catalog, Ncatalog, TRUE);
    76160  MARKTIME("-- set up image indexes: %f sec\n", dtime);
    77161
    78   findMosaics (catalog, Ncatalog);  /* also sets Grid values */
     162  findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
    79163  MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
    80164
    … …  
    123207      plot_chisq (catalog, Ncatalog);
    124208    }
    125     if (i % 6 == 1) rationalize_mosaics (catalog, Ncatalog);
     209    // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
    126210    // if (i % 6 == 1) rationalize_images ();
    127211    if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE);
    … …  
    148232 
    149233  if (USE_GRID) dump_grid ();
    150   if (!UPDATE) exit (0);
    151234
    152235  /* set Mcal & Mmos for bad images */
    … …  
    155238  MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    156239
     240  setMcalFinal (); // copy per-mosaic calibrations to the images
     241
     242  if (SAVE_IMAGE_UPDATES) {
     243
     244    FITS_DB dbX;
     245    off_t *LineNumber;
     246    Image *image, *subset;
     247    off_t Nimage, Nsubset, i, Nx;
     248    char filename[1024];
     249
     250    gfits_db_init (&dbX);
     251    dbX.lockstate = LCK_XCLD;
     252    dbX.timeout   = 60.0;
     253    dbX.mode      = db.mode;
     254    dbX.format    = db.format;
     255
     256    snprintf (filename, 1024, "%s.bck", ImageCat);
     257    if (!gfits_db_lock (&dbX, filename)) {
     258      fprintf (stderr, "can't lock backup image image catalog\n");
     259      return (FALSE);
     260    }
     261   
     262    // apply the changes from the image subset to the full table:
     263
     264    // convert database table to internal structure (binary to Image)
     265    // 'image' points to the same memory as db->ftable->buffer
     266    image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
     267    if (!image) {
     268        fprintf (stderr, "ERROR: failed to read images\n");
     269        exit (2);
     270    }
     271    gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
     272
     273    subset = getimages (&Nsubset, &LineNumber);
     274    for (i = 0; i < Nsubset; i++) {
     275      if (LineNumber[i] == -1) continue;
     276      memcpy (&image[LineNumber[i]], &subset[i], Nx);
     277    }
     278
     279    // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
     280    // the entire Image table, but with the updates in place I think this says: re-work the
     281    // ftable/vtable usage in this program to be more sensible...
     282    gfits_copy_header (&db.header,  &dbX.header);
     283    gfits_copy_matrix (&db.matrix,  &dbX.matrix);
     284    gfits_copy_header (&db.theader, &dbX.theader);
     285    gfits_copy_ftable (&db.ftable,  &dbX.ftable);
     286
     287    dbX.ftable.header = &dbX.theader;
     288    dbX.virtual = FALSE;
     289
     290    // copy Images.dat to another structure
     291    dvo_image_save (&dbX, VERBOSE);
     292    dvo_image_unlock (&dbX);
     293    MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
     294  }
     295
     296  // only change the real database if -update is requested
     297  if (!UPDATE) exit (0);
     298 
     299  reload_images (&db);
     300
    157301  /* at this point, we have correct cal coeffs in the image/mosaic structures */
    158302  for (i = 0; i < Ncatalog; i++) {
    … …  
    160304    dvo_catalog_free (&catalog[i]);
    161305  }
    162   freeImageBins (Ncatalog);
    163   freeMosaicBins (Ncatalog);
     306  freeImageBins (Ncatalog, TRUE);
     307  freeMosaicBins (Ncatalog, TRUE);
    164308  freeGridBins (Ncatalog);
    165309
    … …  
    168312  MARKTIME("-- updated all catalogs: %f sec\n", dtime);
    169313
    170   setMcalFinal ();
    171   reload_images (&db);
    172  
    173314  dvo_image_update (&db, VERBOSE);
    174315  dvo_image_unlock (&db);
Note: See TracChangeset for help on using the changeset viewer.