Index: trunk/Ohana/src/relphot/doc/allsky.txt =================================================================== --- trunk/Ohana/src/relphot/doc/allsky.txt (revision 31667) +++ trunk/Ohana/src/relphot/doc/allsky.txt (revision 31668) @@ -1,2 +1,34 @@ + +2011.05.10 + + Looking into them memory usage of relphot: + AverageTiny = 32 bytes (including 8 byte padding) + MeasureTiny = 56 bytes (including 8 byte padding) + Image Indexes: + bin : Nmeasure * off_t + clist : Nmeasure * off_t + mlist : Nmeasure * off_t + Mosaic Indexes: + bin : Nmeasure * off_t + clist : Nmeasure * off_t + mlist : Nmeasure * off_t + Total: + 32 bytes * Naverage (values) + 56 bytes * Nmeasure (values) + 64 bytes * Nmeasure (indexes) + + In my test run, I have 3.3M average + 37.6M measure: + 2.1GB values + 2.2GB indexes + + This matches the memory usage graph. But, after the indexes have + been loaded, the memory footprint grows (slowly) by another 1GB as + the processing runs along. + + +using 23153493 of 640663697 stars (192585953 of 289672671 measurements) +using 3306896 of 48399614 stars ( 37582679 of 193863612 measurements) + + 2011.04.12 Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 31667) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 31668) @@ -63,4 +63,5 @@ int RESET; int UPDATE; +int SAVE_IMAGE_UPDATES; int PLOTSTUFF; int SAVEPLOT; @@ -139,6 +140,6 @@ void dump_grid PROTO((void)); int edge_check PROTO((double *x1, double *y1, double *x2, double *y2)); -void findImages PROTO((Catalog *catalog, int Ncatalog)); -int findMosaics PROTO((Catalog *catalog, int Ncatalog)); +void findImages PROTO((Catalog *catalog, int Ncatalog, int doImageList)); +int findMosaics PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); void set_db (FITS_DB *in); @@ -149,6 +150,6 @@ void freeGridBins PROTO((int Ncatalog)); -void freeImageBins PROTO((int Ncatalog)); -void freeMosaicBins PROTO((int Ncatalog)); +void freeImageBins PROTO((int Ncatalog, int doImageList)); +void freeMosaicBins PROTO((int Ncatalog, int doMosaicList)); void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); int gcatalog PROTO((Catalog *catalog, int FINAL)); @@ -160,11 +161,11 @@ float getMrel PROTO((Catalog *catalog, off_t meas, int cat)); Image *getimage PROTO((off_t N)); -Image *getimages PROTO((off_t *N)); +Image *getimages PROTO((off_t *N, off_t **LineNumber)); void global_stats PROTO((Catalog *catalog, int Ncatalog)); void initGrid PROTO((int dX, int dY)); void initGridBins PROTO((Catalog *catalog, int Ncatalog)); -void initImageBins PROTO((Catalog *catalog, int Ncatalog)); -void initImages PROTO((Image *input, off_t N)); -void initMosaicBins PROTO((Catalog *catalog, int Ncatalog)); +void initImageBins PROTO((Catalog *catalog, int Ncatalog, int doImageList)); +void initImages PROTO((Image *input, off_t *LineNumber, off_t N)); +void initMosaicBins PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); void initMosaicGrid PROTO((Image *image, off_t Nimage)); void initMosaics PROTO((Image *image, off_t Nimage)); @@ -179,6 +180,6 @@ int main PROTO((int argc, char **argv)); void mark_images PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage)); -void matchImage PROTO((Catalog *catalog, off_t meas, int cat)); -void matchMosaics PROTO((Catalog *catalog, off_t meas, int cat)); +void matchImage PROTO((Catalog *catalog, off_t meas, int cat, int doImageList)); +void matchMosaics PROTO((Catalog *catalog, off_t meas, int cat, int doMosaicList)); double opening_angle PROTO((double x1, double y1, double x2, double y2, double x3, double y3)); void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 31667) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 31668) @@ -9,4 +9,5 @@ static Image *image; // list of available images static off_t Nimage; // number of available images +static off_t *LineNumber; // match of subset to full image table // to search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted index @@ -14,7 +15,8 @@ static off_t *imageIdx; // list of index for image IDs -Image *getimages (off_t *N) { +Image *getimages (off_t *N, off_t **line_number) { *N = Nimage; + if (line_number) *line_number = LineNumber; return (image); } @@ -24,9 +26,10 @@ } -void initImages (Image *input, off_t N) { +void initImages (Image *input, off_t *line_number, off_t N) { off_t i; image = input; + LineNumber = line_number; Nimage = N; @@ -66,5 +69,5 @@ } -void initImageBins (Catalog *catalog, int Ncatalog) { +void initImageBins (Catalog *catalog, int Ncatalog, int doImageList) { off_t i, j; @@ -76,18 +79,20 @@ } - ALLOCATE (Nlist, off_t, Nimage); - ALLOCATE (NLIST, off_t, Nimage); - ALLOCATE (clist, off_t *, Nimage); - ALLOCATE (mlist, off_t *, Nimage); - - for (i = 0; i < Nimage; i++) { - Nlist[i] = 0; - NLIST[i] = 100; - ALLOCATE (clist[i], off_t, NLIST[i]); - ALLOCATE (mlist[i], off_t, NLIST[i]); - } -} - -void freeImageBins (int Ncatalog) { + if (doImageList) { + ALLOCATE (Nlist, off_t, Nimage); + ALLOCATE (NLIST, off_t, Nimage); + ALLOCATE (clist, off_t *, Nimage); + ALLOCATE (mlist, off_t *, Nimage); + + for (i = 0; i < Nimage; i++) { + Nlist[i] = 0; + NLIST[i] = 100; + ALLOCATE (clist[i], off_t, NLIST[i]); + ALLOCATE (mlist[i], off_t, NLIST[i]); + } + } +} + +void freeImageBins (int Ncatalog, int doImageList) { off_t i; @@ -97,16 +102,19 @@ } free (bin); - for (i = 0; i < Nimage; i++) { - free (clist[i]); - free (mlist[i]); - } - free (clist); - free (mlist); - free (Nlist); - free (NLIST); + + if (doImageList) { + for (i = 0; i < Nimage; i++) { + free (clist[i]); + free (mlist[i]); + } + free (clist); + free (mlist); + free (Nlist); + free (NLIST); + } } /* select all image equivalent to the active photcode set */ -void findImages (Catalog *catalog, int Ncatalog) { +void findImages (Catalog *catalog, int Ncatalog, int doImageList) { off_t j; @@ -122,5 +130,5 @@ } if (!found) continue; - matchImage (catalog, j, i); + matchImage (catalog, j, i, doImageList); Nmatch ++; } @@ -178,5 +186,5 @@ } -void matchImage (Catalog *catalog, off_t meas, int cat) { +void matchImage (Catalog *catalog, off_t meas, int cat, int doImageList) { off_t idx, ID; @@ -204,16 +212,18 @@ bin[cat][meas] = idx; - // index for image, Nentry -> catalog - clist[idx][Nlist[idx]] = cat; - - // index for image, Nentry -> measure - mlist[idx][Nlist[idx]] = meas; - Nlist[idx] ++; - - if (Nlist[idx] == NLIST[idx]) { - NLIST[idx] += 100; - REALLOCATE (clist[idx], off_t, NLIST[idx]); - REALLOCATE (mlist[idx], off_t, NLIST[idx]); - } + if (doImageList) { + // index for image, Nentry -> catalog + clist[idx][Nlist[idx]] = cat; + + // index for image, Nentry -> measure + mlist[idx][Nlist[idx]] = meas; + Nlist[idx] ++; + + if (Nlist[idx] == NLIST[idx]) { + NLIST[idx] += 100; + REALLOCATE (clist[idx], off_t, NLIST[idx]); + REALLOCATE (mlist[idx], off_t, NLIST[idx]); + } + } return; @@ -343,4 +353,6 @@ liststats (list, dlist, N, &stats); + + float CLOUD_TOLERANCE = 0.01; image[i].Mcal = stats.mean; image[i].dMcal = stats.error; @@ -348,4 +360,8 @@ image[i].nFitPhotom = N; image[i].Xm = 100.0*log10(stats.chisq); + + if (image[i].Mcal < -CLOUD_TOLERANCE) { + image[i].Mcal = 0.0; + } } free (list); Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 31667) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 31668) @@ -195,10 +195,12 @@ if (!MOSAIC_ZEROPT) return; - image = getimages (&Nimage); - + image = getimages (&Nimage, NULL); + + // copy the mosaic results to the images. set the mosaic Mcal to 0.0 since we have moved its + // impact to the images for (i = 0; i < Nmosaic; i++) { for (j = 0; j < Nimlist[i]; j++) { im = imlist[i][j]; - image[im].Mcal = mosaic[i].Mcal; + image[im].Mcal += mosaic[i].Mcal; image[im].dMcal = mosaic[i].dMcal; image[im].Xm = mosaic[i].Xm; @@ -208,8 +210,9 @@ image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR); } + mosaic[i].Mcal = 0.0; } } -void initMosaicBins (Catalog *catalog, int Ncatalog) { +void initMosaicBins (Catalog *catalog, int Ncatalog, int doMosaicList) { off_t i, j; @@ -224,19 +227,21 @@ } - /* mosaic -> measure */ - ALLOCATE (Nlist, off_t, Nmosaic); - ALLOCATE (NLIST, off_t, Nmosaic); - ALLOCATE (clist, int *, Nmosaic); - ALLOCATE (mlist, off_t *, Nmosaic); - - for (i = 0; i < Nmosaic; i++) { - Nlist[i] = 0; - NLIST[i] = 100; - ALLOCATE (clist[i], int, NLIST[i]); - ALLOCATE (mlist[i], off_t, NLIST[i]); - } -} - -void freeMosaicBins (int Ncatalog) { + if (doMosaicList) { + /* mosaic -> measure */ + ALLOCATE (Nlist, off_t, Nmosaic); + ALLOCATE (NLIST, off_t, Nmosaic); + ALLOCATE (clist, int *, Nmosaic); + ALLOCATE (mlist, off_t *, Nmosaic); + + for (i = 0; i < Nmosaic; i++) { + Nlist[i] = 0; + NLIST[i] = 100; + ALLOCATE (clist[i], int, NLIST[i]); + ALLOCATE (mlist[i], off_t, NLIST[i]); + } + } +} + +void freeMosaicBins (int Ncatalog, int doMosaicList) { off_t i; @@ -250,16 +255,18 @@ free (bin); - /* mosaic -> measure */ - for (i = 0; i < Nmosaic; i++) { - free (clist[i]); - free (mlist[i]); - } - free (Nlist); - free (NLIST); - free (clist); - free (mlist); -} - -int findMosaics (Catalog *catalog, int Ncatalog) { + if (doMosaicList) { + /* mosaic -> measure */ + for (i = 0; i < Nmosaic; i++) { + free (clist[i]); + free (mlist[i]); + } + free (Nlist); + free (NLIST); + free (clist); + free (mlist); + } +} + +int findMosaics (Catalog *catalog, int Ncatalog, int doMosaicList) { int i, ecode, found, Ns; @@ -281,5 +288,5 @@ } if (!found) continue; - matchMosaics (catalog, j, i); + matchMosaics (catalog, j, i, doMosaicList); Nmatch ++; } @@ -289,5 +296,5 @@ } -void matchMosaics (Catalog *catalog, off_t meas, int cat) { +void matchMosaics (Catalog *catalog, off_t meas, int cat, int doMosaicList) { off_t idx, ID, mosID; @@ -322,13 +329,15 @@ bin[cat][meas] = mosID; - clist[mosID][Nlist[mosID]] = cat; - mlist[mosID][Nlist[mosID]] = meas; - Nlist[mosID] ++; + if (doMosaicList) { + clist[mosID][Nlist[mosID]] = cat; + mlist[mosID][Nlist[mosID]] = meas; + Nlist[mosID] ++; - if (Nlist[mosID] == NLIST[mosID]) { - NLIST[mosID] += 100; - REALLOCATE (clist[mosID], int, NLIST[mosID]); - REALLOCATE (mlist[mosID], off_t, NLIST[mosID]); - } + if (Nlist[mosID] == NLIST[mosID]) { + NLIST[mosID] += 100; + REALLOCATE (clist[mosID], int, NLIST[mosID]); + REALLOCATE (mlist[mosID], off_t, NLIST[mosID]); + } + } return; } @@ -361,5 +370,5 @@ if (FREEZE_MOSAICS) return (FALSE); - image = getimages (&N); + image = getimages (&N, NULL); Nsecfilt = GetPhotcodeNsecfilt (); @@ -442,4 +451,6 @@ liststats (list, dlist, N, &stats); if (PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas, N); + + float CLOUD_TOLERANCE = 0.01; mosaic[i].Mcal = stats.mean; mosaic[i].dMcal = stats.error; @@ -447,4 +458,8 @@ mosaic[i].nFitPhotom = N; mosaic[i].Xm = 100.0*log10(stats.chisq); + + if (mosaic[i].Mcal < -CLOUD_TOLERANCE) { + mosaic[i].Mcal = 0.0; + } } free (list); @@ -482,8 +497,7 @@ if (FREEZE_MOSAICS) return (FALSE); - image = getimages (&Nimage); + image = getimages (&Nimage, NULL); ALLOCATE (imageOffset, double, Nmosaic); - ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic ALLOCATE (NSlist, int, Nmosaic); // number of stars on mosaic @@ -621,7 +635,7 @@ float dM = 0.0; for (j = 0; j < NSlist[i]; j++) { - Nstars = Slist[i][j]; - if (starNcount[Nstars][Ns] > 1) { - dM += (starOffset[Nstars][Ns] / starNcount[Nstars][Ns]); + k = Slist[i][j]; + if (starNcount[k][Ns] > 1) { + dM += (starOffset[k][Ns] / starNcount[k][Ns]); } } @@ -647,11 +661,12 @@ free (starNcount); + free (seclist); + free (NSlist); + free (NSLIST); for (i = 0; i < Nmosaic; i++){ free (Slist[i]); } - free (NSlist); - free (NSLIST); + free (Slist); free (imageOffset); - free (seclist); return (TRUE); } @@ -853,5 +868,5 @@ if (!MOSAIC_ZEROPT) return; - image = getimages (&Nimage); + image = getimages (&Nimage, NULL); N = 0; Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 31667) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 31668) @@ -587,4 +587,6 @@ initstats (STATMODE); if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave); + free (list); + free (dlist); free (ilist); free (tlist); Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 31667) +++ trunk/Ohana/src/relphot/src/args.c (revision 31668) @@ -128,4 +128,10 @@ } + SAVE_IMAGE_UPDATES = TRUE; + if ((N = get_argument (argc, argv, "-skip-image-updates"))) { + remove_argument (N, &argc, argv); + SAVE_IMAGE_UPDATES = FALSE; + } + MaxDensityUse = FALSE; if ((N = get_argument (argc, argv, "-max-density"))) { Index: trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/bcatalog.c (revision 31667) +++ trunk/Ohana/src/relphot/src/bcatalog.c (revision 31668) @@ -5,7 +5,7 @@ int CopyAverageTiny (AverageTiny *averageT, Average *average) { - averageT[0].R = average[0].R; - averageT[0].D = average[0].D; - averageT[0].flags = average[0].flags; + averageT[0].R = average[0].R; + averageT[0].D = average[0].D; + averageT[0].flags = average[0].flags; averageT[0].Nmeasure = average[0].Nmeasure; averageT[0].measureOffset = average[0].measureOffset; @@ -293,13 +293,44 @@ off_t i; + AverageTiny *averageT; + Average *average; + + MeasureTiny *measureT; + Measure *measure; + ALLOCATE (catalog[0].measureT, MeasureTiny, catalog[0].Nmeasure); ALLOCATE (catalog[0].averageT, AverageTiny, catalog[0].Naverage); + average = catalog[0].average; + averageT = catalog[0].averageT; + + measure = catalog[0].measure; + measureT = catalog[0].measureT; + for (i = 0; i < catalog[0].Naverage; i++) { - CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]); + // CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]); + averageT[i].R = average[i].R; + averageT[i].D = average[i].D; + averageT[i].flags = average[i].flags; + averageT[i].Nmeasure = average[i].Nmeasure; + averageT[i].measureOffset = average[i].measureOffset; } for (i = 0; i < catalog[0].Nmeasure; i++) { - CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]); + // CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]); + measureT[i].dR = measure[i].dR; + measureT[i].dD = measure[i].dD; + measureT[i].M = measure[i].M; + measureT[i].Mcal = measure[i].Mcal; + measureT[i].dM = measure[i].dM; + measureT[i].airmass = measure[i].airmass; + measureT[i].Xccd = measure[i].Xccd; + measureT[i].Yccd = measure[i].Yccd; + measureT[i].t = measure[i].t; + measureT[i].dt = measure[i].dt; + measureT[i].averef = measure[i].averef; + measureT[i].imageID = measure[i].imageID; + measureT[i].dbFlags = measure[i].dbFlags; + measureT[i].photcode = measure[i].photcode; } @@ -309,7 +340,7 @@ int free_tiny_values (Catalog *catalog) { - free (catalog[0].averageT); - free (catalog[0].measureT); - return (TRUE); -} - + if (catalog[0].averageT) free (catalog[0].averageT); + if (catalog[0].measureT) free (catalog[0].measureT); + return (TRUE); +} + Index: trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/load_catalogs.c (revision 31667) +++ trunk/Ohana/src/relphot/src/load_catalogs.c (revision 31668) @@ -3,5 +3,5 @@ Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) { - int i, Nmeas, Nstar, Nmeas_total, Nstar_total; + off_t i, Nmeas, Nstar, Nmeas_total, Nstar_total; Catalog *catalog, tcatalog; @@ -14,4 +14,5 @@ // load data from each region file, only use bright stars for (i = 0; i < skylist[0].Nregions; i++) { + dvo_catalog_init (&catalog[i], TRUE); // set up the basic catalog info @@ -49,8 +50,8 @@ } if (Nstar < 2) { - fprintf (stderr, "warning: insufficient stars %d\n", Nstar); + fprintf (stderr, "warning: insufficient stars "OFF_T_FMT"\n", Nstar); } - fprintf (stderr, "using %d of %d stars (%d of %d measurements)\n", Nstar, Nstar_total, Nmeas, Nmeas_total); + 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); if (Nstar < 1) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n"); Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 31667) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 31668) @@ -25,14 +25,6 @@ SkyTableSetFilenames (sky, CATDIR, "cpt"); - // determine the populated SkyRegions overlapping the requested area - if (RegionSelect) { - skylist = SkyListByPatch (sky, -1, region); - } else { - Nchar = strlen(regionName); - if (!strcmp (®ionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0; - skylist = SkyListByName (sky, regionName); - } - // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); if (!image) { @@ -42,16 +34,67 @@ MARKTIME("read image table: %f sec\n", dtime); - // select the images which overlap the selected sky regions - subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); - MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); + // determine the populated SkyRegions overlapping the requested area + if (RegionSelect) { + if (region[0].Rmin > region[0].Rmax) { + SkyRegion subregion; + subregion = *region; + subregion.Rmin = 0.0; subregion.Rmax = region[0].Rmax; + skylist = SkyListByPatch (sky, -1, &subregion); + + subregion.Rmin = region[0].Rmin; subregion.Rmax = 360.0; + SkyList *extraList = SkyListByPatch (sky, -1, &subregion); + + // select the images which overlap the selected sky regions + // 'subset' points to a new copy of the data (different from 'image') + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); + + Image *subsetExtra; + off_t NsubsetExtra; + off_t *LineNumberExtra, i; + + subsetExtra = select_images (extraList, image, Nimage, &LineNumberExtra, &NsubsetExtra); + MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); + + REALLOCATE (subset, Image, MAX (Nsubset + NsubsetExtra, 1)); + REALLOCATE (LineNumber, off_t, MAX (Nsubset + NsubsetExtra, 1)); + + for (i = 0; i < NsubsetExtra; i++) { + subset[i+Nsubset] = subsetExtra[i]; + LineNumber[i+Nsubset] = LineNumberExtra[i]; + } + Nsubset += NsubsetExtra; + MARKTIME("selected %d overlapping images: %f sec\n", (int) NsubsetExtra, dtime); + } else { + skylist = SkyListByPatch (sky, -1, region); + + // select the images which overlap the selected sky regions + // 'subset' points to a new copy of the data (different from 'image') + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); + } + } else { + Nchar = strlen(regionName); + if (!strcmp (®ionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0; + skylist = SkyListByName (sky, regionName); + + // select the images which overlap the selected sky regions + // 'subset' points to a new copy of the data (different from 'image') + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); + } + // generate db->vtable from db->ftable based on the selection // XXX does this simply duplicate the memory needlessly? we recreate these lines // in reload_images. If we had saved the line numbers, we could avoid this - gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset); - MARKTIME("converted ftable to vtable: %f sec\n", dtime); + // vtable points *another* copy of the subset rows + // (the later call to 'reload_images' copies the subset elements back on top of + // the rows of the vtable) + // gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset); + // MARKTIME("converted ftable to vtable: %f sec\n", dtime); // save the subset of images in the static reference in ImageOps, set up indexes - initImages (subset, Nsubset); + initImages (subset, LineNumber, Nsubset); MARKTIME("init images: %f sec\n", dtime); @@ -66,9 +109,10 @@ Image *image; - off_t Nimage, Nx, i; + off_t Nimage, Nx, i, *LineNumber; VTable *vtable; - image = getimages (&Nimage); + image = getimages (&Nimage, &LineNumber); + gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nimage); vtable = &db[0].vtable; Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 31667) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 31668) @@ -18,5 +18,6 @@ double time1 = 0.0; double time2 = 0.0; - double time3 = 0.0; + double time3a = 0.0; + double time3b = 0.0; double time4 = 0.0; double time5 = 0.0; @@ -57,12 +58,13 @@ populate_tiny_values(&catalog); + TIMESTAMP(time3a); - initImageBins (&catalog, 1); - initMosaicBins (&catalog, 1); + initImageBins (&catalog, 1, FALSE); + initMosaicBins (&catalog, 1, FALSE); initGridBins (&catalog, 1); - TIMESTAMP(time3); + TIMESTAMP(time3b); - findImages (&catalog, 1); // FX - findMosaics (&catalog, 1); // + findImages (&catalog, 1, FALSE); + findMosaics (&catalog, 1, FALSE); TIMESTAMP(time4); @@ -77,16 +79,17 @@ TIMESTAMP(time6); - freeImageBins (1); - freeMosaicBins (1); + freeImageBins (1, FALSE); + freeMosaicBins (1, FALSE); freeGridBins (1); TIMESTAMP(time7); } - fprintf (stderr, "time1 %f : find catalog\n", time1); - fprintf (stderr, "time2 %f : load catalog\n", time2); - fprintf (stderr, "time3 %f : init imbins\n", time3); - fprintf (stderr, "time4 %f : find images\n", time4); - fprintf (stderr, "time5 %f : set Mrel\n", time5); - fprintf (stderr, "time6 %f : save catalog\n", time6); - fprintf (stderr, "time7 %f : free catalog\n", time7); + fprintf (stderr, "time1 %f : find catalog\n", time1); + fprintf (stderr, "time2 %f : load catalog\n", time2); + fprintf (stderr, "time3a %f : init imbins\n", time3a); + fprintf (stderr, "time3b %f : init imbins\n", time3b); + fprintf (stderr, "time4 %f : find images\n", time4); + fprintf (stderr, "time5 %f : set Mrel\n", time5); + fprintf (stderr, "time6 %f : save catalog\n", time6); + fprintf (stderr, "time7 %f : free catalog\n", time7); } Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 31667) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 31668) @@ -1,4 +1,10 @@ # include "relphot.h" # define USE_DIRECT 0 + +# define TIMESTAMP(TIME) \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + TIME += dtime; \ + gettimeofday (&start, (void *) NULL); # define MARKTIME(MSG,...) { \ @@ -17,4 +23,82 @@ gettimeofday (&start, (void *) NULL); + + // XXX quick and stupid test: + if (0) { + int i, j; + off_t *Nlist, *NLIST, **mlist; + off_t *NlistI, *NLISTI, **mlistI; + int **clist; + int **clistI; + float dtime; + float time1 = 0.0; + float time2 = 0.0; + + int Nmosaic = 10000; + int Nimage = 600000; + + for (j = 0; j < 19000; j++) { + + // mosaic indexes + ALLOCATE (Nlist, off_t, Nmosaic); + ALLOCATE (NLIST, off_t, Nmosaic); + ALLOCATE (clist, int *, Nmosaic); + ALLOCATE (mlist, off_t *, Nmosaic); + + for (i = 0; i < Nmosaic; i++) { + Nlist[i] = 0; + NLIST[i] = 100; + ALLOCATE (clist[i], int, NLIST[i]); + ALLOCATE (mlist[i], off_t, NLIST[i]); + } + + // image indexes + ALLOCATE (NlistI, off_t, Nimage); + ALLOCATE (NLISTI, off_t, Nimage); + ALLOCATE (clistI, int *, Nimage); + ALLOCATE (mlistI, off_t *, Nimage); + + for (i = 0; i < Nimage; i++) { + NlistI[i] = 0; + NLISTI[i] = 100; + ALLOCATE (clistI[i], int, NLISTI[i]); + ALLOCATE (mlistI[i], off_t, NLISTI[i]); + } + + TIMESTAMP(time1); + + // free mosaic indexes + for (i = 0; i < Nmosaic; i++) { + free (clist[i]); + free (mlist[i]); + } + free (Nlist); + free (NLIST); + free (clist); + free (mlist); + + // free image indexes + for (i = 0; i < Nimage; i++) { + free (clistI[i]); + free (mlistI[i]); + } + free (NlistI); + free (NLISTI); + free (clistI); + free (mlistI); + + TIMESTAMP(time2); + + if (j % 100 == 0) { + fprintf (stderr, "done with %d\n", j); + fprintf (stderr, "time1 %f : init mosaic\n", time1); + fprintf (stderr, "time2 %f : free mosaic\n", time2); + } + } + + fprintf (stderr, "time1 %f : init mosaic\n", time1); + fprintf (stderr, "time2 %f : free mosaic\n", time2); + exit (1); + } /* get configuration info, args */ @@ -66,15 +150,15 @@ /* match measurements with images, mosaics */ - initImageBins (catalog, Ncatalog); + initImageBins (catalog, Ncatalog, TRUE); MARKTIME("-- make image bins: %f sec\n", dtime); - initMosaicBins (catalog, Ncatalog); + initMosaicBins (catalog, Ncatalog, TRUE); initGridBins (catalog, Ncatalog); initMrel (catalog, Ncatalog); - findImages (catalog, Ncatalog); + findImages (catalog, Ncatalog, TRUE); MARKTIME("-- set up image indexes: %f sec\n", dtime); - findMosaics (catalog, Ncatalog); /* also sets Grid values */ + findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); @@ -123,5 +207,5 @@ plot_chisq (catalog, Ncatalog); } - if (i % 6 == 1) rationalize_mosaics (catalog, Ncatalog); + // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); // if (i % 6 == 1) rationalize_images (); if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE); @@ -148,5 +232,4 @@ if (USE_GRID) dump_grid (); - if (!UPDATE) exit (0); /* set Mcal & Mmos for bad images */ @@ -155,4 +238,65 @@ MARKTIME("-- finalize Mcal values: %f sec\n", dtime); + setMcalFinal (); // copy per-mosaic calibrations to the images + + if (SAVE_IMAGE_UPDATES) { + + FITS_DB dbX; + off_t *LineNumber; + Image *image, *subset; + off_t Nimage, Nsubset, i, Nx; + char filename[1024]; + + gfits_db_init (&dbX); + dbX.lockstate = LCK_XCLD; + dbX.timeout = 60.0; + dbX.mode = db.mode; + dbX.format = db.format; + + snprintf (filename, 1024, "%s.bck", ImageCat); + if (!gfits_db_lock (&dbX, filename)) { + fprintf (stderr, "can't lock backup image image catalog\n"); + return (FALSE); + } + + // apply the changes from the image subset to the full table: + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1, &Nx); + + subset = getimages (&Nsubset, &LineNumber); + for (i = 0; i < Nsubset; i++) { + if (LineNumber[i] == -1) continue; + memcpy (&image[LineNumber[i]], &subset[i], Nx); + } + + // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump + // the entire Image table, but with the updates in place I think this says: re-work the + // ftable/vtable usage in this program to be more sensible... + gfits_copy_header (&db.header, &dbX.header); + gfits_copy_matrix (&db.matrix, &dbX.matrix); + gfits_copy_header (&db.theader, &dbX.theader); + gfits_copy_ftable (&db.ftable, &dbX.ftable); + + dbX.ftable.header = &dbX.theader; + dbX.virtual = FALSE; + + // copy Images.dat to another structure + dvo_image_save (&dbX, VERBOSE); + dvo_image_unlock (&dbX); + MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); + } + + // only change the real database if -update is requested + if (!UPDATE) exit (0); + + reload_images (&db); + /* at this point, we have correct cal coeffs in the image/mosaic structures */ for (i = 0; i < Ncatalog; i++) { @@ -160,6 +304,6 @@ dvo_catalog_free (&catalog[i]); } - freeImageBins (Ncatalog); - freeMosaicBins (Ncatalog); + freeImageBins (Ncatalog, TRUE); + freeMosaicBins (Ncatalog, TRUE); freeGridBins (Ncatalog); @@ -168,7 +312,4 @@ MARKTIME("-- updated all catalogs: %f sec\n", dtime); - setMcalFinal (); - reload_images (&db); - dvo_image_update (&db, VERBOSE); dvo_image_unlock (&db);