Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/include/relphot.h (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/include/relphot.h (revision 31258) @@ -90,4 +90,7 @@ PhotCode *refPhotcode; + +int MaxDensityUse; +double MaxDensityValue; int AreaSelect; @@ -223,2 +226,5 @@ void relphot_usage (void); void relphot_help (int argc, char **argv); + +int rationalize_mosaics (); +int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog); Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/MosaicOps.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/MosaicOps.c (revision 31258) @@ -302,5 +302,6 @@ return; } - fprintf (stderr, "missed measurement\n"); + // if we miss a measurement, it might be a REF type of value + // fprintf (stderr, "missed measurement\n"); return; } @@ -403,5 +404,5 @@ mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]); if (mark) { - if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]); } + if (1 || VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]); } mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW; Nfew ++; @@ -433,4 +434,144 @@ } +// When we rationalize the images/mosaics, we are driving the negative cloud images back +// to 0.0. At the same time, we make a guess to the effective impact on all other images, +// driven by the coupling of common stars. +int rationalize_mosaics (Catalog *catalog, int Ncatalog) { + + double *imageOffset, *starOffset; + int *starNcount; + int **Slist, *NSlist, *NSLIST; + int i, j, k, m, nMos; + + off_t Nimage; + Image *image; + + // set a test value for now + float CLOUD_TOLERANCE = 0.025; + + if (!MOSAIC_ZEROPT) return (FALSE); + if (FREEZE_MOSAICS) return (FALSE); + + image = getimages (&Nimage); + + 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 + ALLOCATE (NSLIST, int, Nmosaic); // number of Slist entries allocated + + // find the images / mosaics with negative clouds and save their offset + for (i = 0; i < Nmosaic; i++) { + + NSlist[i] = 0; + NSLIST[i] = 100; + ALLOCATE (Slist[i], int, NSLIST[i]); + + imageOffset[i] = 0.0; + + if (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE) { + fprintf (stderr, "cloud-free: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + } + if (mosaic[i].Mcal < -CLOUD_TOLERANCE) { + imageOffset[i] = -mosaic[i].Mcal; + // NOTE the negative sign: down below, we are going to add in the negative of Mcal + // to this image, and the propagated mean values for other images + fprintf (stderr, "anti-clouds: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + } + if (mosaic[i].Mcal > CLOUD_TOLERANCE) { + fprintf (stderr, "cloudy : %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + } + } + + // allocate an array for star offset + int Nstars = 0; + int NSTARS = 1000; + ALLOCATE (starOffset, double, NSTARS); + ALLOCATE (starNcount, int, NSTARS); + + // find the mean offset for each star + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + starOffset[Nstars] = 0.0; + starNcount[Nstars] = 0; + + // skip bad stars + if (catalog[i].average[j].flags & STAR_BAD) continue; + m = catalog[i].average[j].measureOffset; + + // determine the mosaic for each measurement + for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + + // skip unused measurements + if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; + + // skip REF measurements (not tied to an image) + if (getImageEntry (m, i) < 0) continue; + + // find the source of this measurement (skip unassigned measurements) + nMos = bin[i][m]; + if (nMos == -1) continue; + + // accumulate the offsets from the negative cloud images (others have 0.0 value) + starOffset[Nstars] += imageOffset[nMos]; + starNcount[Nstars] ++; + + // record the mosaic->star reference + Slist[nMos][NSlist[nMos]] = Nstars; + NSlist[nMos] ++; + if (NSlist[nMos] == NSLIST[nMos]) { + NSLIST[nMos] += 100; + REALLOCATE (Slist[nMos], int, NSLIST[nMos]); + } + } + Nstars ++; + if (Nstars == NSTARS) { + NSTARS += 1000; + REALLOCATE (starOffset, double, NSTARS); + REALLOCATE (starNcount, int, NSTARS); + } + } + } + + // find the mean offset of the images without negative clouds + for (i = 0; i < Nmosaic; i++) { + + // a negative cloud image (cloud: Mcal > 0; anti-clouds: Mcal < 0; imageOffset = -Mcal) + if (imageOffset[i] > 0.0) continue; + + // we need to actually have cross-references to count + if (NSlist[i] < 2) continue; + + float dM = 0.0; + for (j = 0; j < NSlist[i]; j++) { + Nstars = Slist[i][j]; + if (starNcount[Nstars] > 1) { + dM += (starOffset[Nstars] / starNcount[Nstars]); + } + } + imageOffset[i] = dM / NSlist[i]; + fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[imlist[i][0]].name, dM, NSlist[i], imageOffset[i]); + } + + for (i = 0; i < Nmosaic; i++) { + fprintf (stderr, "correction: %s : %f\n", image[imlist[i][0]].name, imageOffset[i]); + } + + // apply all offset values to the mosaics + // find the images / mosaics with negative clouds and save their offset + for (i = 0; i < Nmosaic; i++) { + mosaic[i].Mcal += imageOffset[i]; + } + + for (i = 0; i < Nmosaic; i++){ + free (Slist[i]); + } + free (NSlist); + free (NSLIST); + free (imageOffset); + + return (TRUE); +} + StatType statsMosaicM (Catalog *catalog) { @@ -563,5 +704,5 @@ void clean_mosaics () { - off_t i, N, mark, Nmark; + off_t i, N, mark, Nmark, Nscatter, Noffset; double *mlist, *slist, *dlist; double MaxOffset, MedOffset, MaxScatter; @@ -592,8 +733,15 @@ fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); - Nmark = 0; + Nmark = Nscatter = Noffset = 0; for (i = 0; i < Nmosaic; i++) { mark = FALSE; - mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset); + if (mosaic[i].dMcal > MaxScatter) { + mark = TRUE; + Nscatter ++; + } + if (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset) { + mark = TRUE; + Noffset ++; + } if (mark) { Nmark ++; @@ -604,5 +752,5 @@ } - fprintf (stderr, OFF_T_FMT" mosaics marked poor\n", Nmark); + fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Nscatter, Noffset); initstats (STATMODE); free (mlist); Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c (revision 31258) @@ -5,4 +5,10 @@ static double *dlist; +// When we rationalize the images/mosaics, we are driving the negative cloud images back +// to 0.0. At the same time, we make a guess to the effective impact on all other images, +// driven by the coupling of common stars. This array carries the impact of those offsets +// on each star +static double *Moffset; + void initMrel (Catalog *catalog, int Ncatalog) { @@ -16,6 +22,7 @@ } - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); + ALLOCATE (list, double, MAX (1, Nmax)); + ALLOCATE (dlist, double, MAX (1, Nmax)); + ALLOCATE (Moffset, double, MAX (1, Nmax)); } Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/args.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/args.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/args.c (revision 31258) @@ -128,4 +128,12 @@ } + MaxDensityUse = FALSE; + if ((N = get_argument (argc, argv, "-max-density"))) { + remove_argument (N, &argc, argv); + MaxDensityValue = atof(argv[N]); + remove_argument (N, &argc, argv); + MaxDensityUse = TRUE; + } + SHOW_PARAMS = FALSE; if ((N = get_argument (argc, argv, "-params"))) { Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/bcatalog.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/bcatalog.c (revision 31258) @@ -1,3 +1,5 @@ # include "relphot.h" + +extern double drand48(); int bcatalog (Catalog *subcatalog, Catalog *catalog) { @@ -70,9 +72,9 @@ // check for galaxies if (!isnan(catalog[0].measure[offset].Map)) { - if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) { - nEXT ++; - } else { - nPSF ++; - } + if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) { + nEXT ++; + } else { + nPSF ++; + } } @@ -139,4 +141,9 @@ assert (PhotNsec == catalog[0].Nsecfilt); + // limit the total number of stars in the catalog + if (MaxDensityUse) { + LimitDensityCatalog (subcatalog, catalog); + } + if (VERBOSE) { fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", @@ -147,2 +154,90 @@ return (TRUE); } + +int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) { + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + gfits_scan (&catalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&catalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax); + + if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); + + float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin)); + assert (AREA > 0); + + off_t Nmax = MaxDensityValue * AREA; + if (subcatalog[0].Naverage <= Nmax) { + if (VERBOSE) { + fprintf (stderr, "subcatalog has less than the max density\n"); + } + return (TRUE); + } + + off_t Naverage = subcatalog[0].Naverage; + + // select a random subset of Nmax stars from subcatalog using Fisher-Yates + + // we are going to select Nmax entries by generating a random-sorted index list + off_t *index, tmp, i, j, ave; + ALLOCATE (index, off_t, Naverage); + for (i = 0; i < Naverage; i++) { + index[i] = i; + } + for (i = 0; i < Naverage; i++) { + j = (Naverage - i) * drand48() + i; // a number between i and Naverage + tmp = index[j]; + index[j] = index[i]; + index[i] = tmp; + } + + // count the number of measurements this selection will yield + off_t NMEASURE = 0; + for (i = 0; i < Nmax; i++) { + ave = index[i]; + NMEASURE += subcatalog[0].average[ave].Nmeasure; + } + + // allocate the output data + ALLOCATE (tmpcatalog.average, Average, Nmax); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * PhotNsec); + ALLOCATE (tmpcatalog.measure, Measure, NMEASURE); + + off_t Nmeasure = 0; + + // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links) + for (i = 0; i < Nmax; i++) { + ave = index[i]; + tmpcatalog.average[i] = subcatalog[0].average[ave]; + tmpcatalog.average[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { + off_t offset = subcatalog[0].average[ave].measureOffset + j; + tmpcatalog.measure[Nmeasure] = subcatalog[0].measure[offset]; + tmpcatalog.measure[Nmeasure].averef = i; + Nmeasure ++; + } + } + + if (VERBOSE) { + fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", + Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, catalog[0].filename); + } + + free (subcatalog[0].average); + free (subcatalog[0].measure); + free (subcatalog[0].secfilt); + + subcatalog[0].average = tmpcatalog.average; + subcatalog[0].measure = tmpcatalog.measure; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nmax; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt; + + return (TRUE); +} Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/initialize.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/initialize.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/initialize.c (revision 31258) @@ -55,4 +55,9 @@ exit (0); } + + // init the random seed + long A, B; + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); } - Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/load_images.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/load_images.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/load_images.c (revision 31258) @@ -44,5 +44,5 @@ // select the images which overlap the selected sky regions subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); - MARKTIME("selected images: %f sec\n", dtime); + MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); // generate db->vtable from db->ftable based on the selection Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/relphot.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/relphot.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/relphot.c (revision 31258) @@ -140,19 +140,27 @@ plot_chisq (catalog, Ncatalog); } - if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); - if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); - if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); - if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); + if (i % 6 == 1) rationalize_mosaics (catalog, Ncatalog); + // if (i % 6 == 1) rationalize_images (); + if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE); + if (i % 6 == 3) clean_stars (catalog, Ncatalog); + if (i % 6 == 5) clean_mosaics (); + if (i % 6 == 5) clean_images (); + + // if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); + // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); global_stats (catalog, Ncatalog); } - SAVEPLOT = TRUE; - plot_scatter (catalog, Ncatalog); - plot_grid (catalog); - plot_mosaics (); - plot_images (); - plot_stars (catalog, Ncatalog); - plot_chisq (catalog, Ncatalog); - + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog); + plot_grid (catalog); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (catalog, Ncatalog); + } + if (USE_GRID) dump_grid (); if (!UPDATE) exit (0); Index: /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/select_images.c =================================================================== --- /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/select_images.c (revision 31257) +++ /branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/select_images.c (revision 31258) @@ -64,4 +64,6 @@ DmaxSkyRegion = -90.0; + // FILE *ftest = fopen ("relphot.dump.dat", "w"); + /* compare with each region file */ for (i = 0; i < skylist[0].Nregions; i++) { @@ -131,4 +133,7 @@ } + // XXX temporary test : record center coords for each accepted and rejected chip/mosaic + // double RAo, DECo; + /* define image corners - note the DIS images (mosaic phu) are special */ if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) { @@ -138,4 +143,5 @@ Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY; Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY; + // XY_to_RD(&RAo, &DECo, 0.0, 0.0, &timage[i].coords); } else { Xi[0] = 0; Yi[0] = 0; @@ -144,10 +150,11 @@ Xi[3] = 0; Yi[3] = timage[i].NY; Xi[4] = 0; Yi[4] = 0; + // XY_to_RD(&RAo, &DECo, 0.5*timage[i].NX, 0.5*timage[i].NY, &timage[i].coords); } found = FALSE; /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */ - double RminImage = 360.0; - double RmaxImage = 0.0; + double RminImage = RmidSkyRegion + 180.0; + double RmaxImage = RmidSkyRegion - 180.0; double DminImage = +90.0; double DmaxImage = -90.0; @@ -218,4 +225,7 @@ found_it: + // XXX We claim this is a good image: write to a test file + // fprintf (ftest, "%s : %lf %lf : %f %f %d %x\n", timage[i].name, RAo, DECo, timage[i].Mcal, timage[i].dMcal, timage[i].nstar, timage[i].flags); + image[nimage] = timage[i]; /* always allow 'few' images to succeed, if possible */ @@ -239,4 +249,6 @@ } MARKTIME("finish image selection: %f sec\n", dtime); + + // fclose (ftest); if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);