Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 33650) +++ trunk/Ohana/src/relphot/Makefile (revision 33651) @@ -1,3 +1,3 @@ -default: relphot +default: relphot relphot_client help: @echo "make options: relphot (default)" @@ -18,5 +18,7 @@ relphot: $(BIN)/relphot.$(ARCH) -install: $(DESTBIN)/relphot +relphot_client: $(BIN)/relphot_client.$(ARCH) + +install: $(DESTBIN)/relphot $(DESTBIN)/relphot_client RELPHOT = \ @@ -24,4 +26,6 @@ $(SRC)/GridOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/ImageSubset.$(ARCH).o \ +$(SRC)/BrightCatalog.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ @@ -48,2 +52,29 @@ $(RELPHOT): $(INC)/relphot.h $(BIN)/relphot.$(ARCH): $(RELPHOT) + +RELPHOT_CLIENT = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/GridOps.$(ARCH).o \ +$(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/ImageSubset.$(ARCH).o \ +$(SRC)/BrightCatalog.$(ARCH).o \ +$(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/StarOps.$(ARCH).o \ +$(SRC)/bcatalog.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/help.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/liststats.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/reload_catalogs.$(ARCH).o \ +$(SRC)/relphot_objects.$(ARCH).o \ +$(SRC)/relphot_client.$(ARCH).o \ +$(SRC)/setExclusions.$(ARCH).o \ +$(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/write_coords.$(ARCH).o + +$(RELPHOT_CLIENT): $(INC)/relphot.h +$(BIN)/relphot_client.$(ARCH): $(RELPHOT_CLIENT) Index: trunk/Ohana/src/relphot/doc/flatcorr.txt =================================================================== --- trunk/Ohana/src/relphot/doc/flatcorr.txt (revision 33651) +++ trunk/Ohana/src/relphot/doc/flatcorr.txt (revision 33651) @@ -0,0 +1,90 @@ + +2012.01.18 + +I have updated libdvo to define a flatcorr table, and have the APIs to +read and write the table. I have also updated setphot to load the +ubercal measurements and to define the flatcorr table and links (from +image -> flatcorr). Now I need to understand where to include this in +relphot. + +One option would be to include it as part of the function 'getMcal'. +This function goes from (cat,meas) -> (image) to find the image zero +point (Mcal). It could also check for a photom_map_id value for that +image, and include that in the calculation (if it exists). + + +Another option would be to add a new function (getMmap or something) +and to include this essentially everywhere getMcal is called. + +2012.01.13 + +Thinking it over, I've decided that the easiest way to move forward on +the implementation is to define a second grid-like correction to +represent the flat-field, and defer unifying the two (conceptually +similar) things to the future. + +So, we start with infrastructure: + +* APIs to read and write the flat-field corrections + + * how to represent the flat-field correction? + * each pixel as a row in a table? + * easy to do the file I/O (one table for all corrections) + * need to unwrap the values into correction[photcode][x][y] + * + * each photcode as a FITS image extention? + * + * each camera as a full image (like Grid)? + + +2012.01.12 + +The Ubercal analysis gives me zero points along with flat-field +corrections defined for certain time ranges as a function of +fractional cell location. I need to include equivalent information +into the DVO tables, both to apply the ubercal values and to do +equivalent fits. This needs to be done in both relphot and setphot. +(Other programs refer to the value of Mcal for each detection; only +relphot and setphot modify or apply the image Mcal to the detection +value). + +In setphot, the needed modification is in update_catalog_setphot, +where the value of Mcal is set from the image. This function needs to +also look up the corresponding flat-correction and include it in that +calculation. It is also necessary to load the externally-defined +description of the flat-correction and to determine which images go +which which corrections. + +In relphot, we would like to be able to fit for the value of interest, +as well as to apply the values loaded from the database. + +In relphot: + +The survey-wide flat-field correction information is similar to the +Grid that is currently used to calculate the flat-field correction +terms, though there are some subtleties. + +Currently, the grid analysis is a bit complicated by the way we relate +the chip coordinates to the camera coordinates. The main problem is +that we do not do a good job of tracking which chip a detection comes +from. + +GridOps loads a configuration table which describes the layout of the +camera. This means the relative orientations of each of the chips. +Each chip is define by its ccdnum value and by its ccdname. + +The code assigns the Grid bin by using the camera layout information +to relate the chip X,Y position to the camera X,Y coordinate. The +grid cells are then defined as a function of camera coordinate. + +This is not a terrible solution, but it conflates two things: the chip +and position in the chip with the layout of the focal plane. To solve +or apply a flat-field term, we don't actually need to know anything +about the location of the chip in the camera. + +In other words, Grid[Xbin][Ybin] is currently a single array assuming +a monotonic relationship between Xbin and Xchip, etc. We could have +Grid[chip][Xbin][Ybin], where Xbin and Ybin come directl from +Xchip,Ychip and have a different function which converts the resulting +image to a camera-coordinate system image. + Index: trunk/Ohana/src/relphot/doc/parallel.txt =================================================================== --- trunk/Ohana/src/relphot/doc/parallel.txt (revision 33651) +++ trunk/Ohana/src/relphot/doc/parallel.txt (revision 33651) @@ -0,0 +1,35 @@ + +2012.02.13 + +Image fields needed by relphot_client (saved by ImageSubsetSave): + +imageID +photom_map_id +flags +Mcal +dMcal +tzero +trate + +(if we have USE_GRID == T): +photcode + + +2012.02.10 + +I have put together the start of parallel dvo examples. I have +'dvodist' to distribute the data (and pull it back), and parallel +modifications to 'setphot'. Here is what I need to do to relphot to +make it work in parallel: + +* load_catalogs.c : a parallel version of this would spawn the remote + jobs in -load mode. these should read from their owned tables, + collect the subset of detections (subset naturally generated by + bcatalog), and write the results to local temp fits files. the main + program would wait until the sub programs have been finished and + then load in all of the subset tables. + +* reload_catalogs : a parallel version needs to write the subset image + table and then run the remote jobs. these would read the image + data, and load and update their local database tables. This is + nearly identical to the code in setphot / setphot_client. Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 33650) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 33651) @@ -7,4 +7,15 @@ # define GRID_V2 # define NO_IMAGE -100 + +// choose off_t or int depending on full-scale relphot analysis resources +// # define IDX_T off_t +# define IDX_T int + +typedef enum { + MODE_NONE = 0, + MODE_LOAD = 1, + MODE_UPDATE = 2, + MODE_UPDATE_OBJECTS = 3, +} ModeType; typedef struct { @@ -18,5 +29,6 @@ short Xm; float secz; - char flags; + float ubercalDist; + unsigned int flags; Coords coords; } Mosaic; @@ -30,19 +42,59 @@ double min; double max; + double Upper80; + double Lower20; double total; int Nmeas; } StatType; +typedef struct { + AverageTiny *average; // array of (minimal) average data + MeasureTiny *measure; // array of (minimal) measure data + SecFilt *secfilt; // array of secfilt data (matched to average by Nsecfilt) + off_t Naverage; + off_t Nmeasure; +} BrightCatalog; + +typedef struct { + float Mcal; + float dMcal; + unsigned int imageID; + unsigned int photom_map_id; + unsigned int flags; + unsigned int tzero; + unsigned char trate; +} ImageSubset; + +typedef struct { + Catalog *catalog; // array of catalogs generated + int NCATALOG; // number of catalogs allocated + int Ncatalog; // number of catalogs generated + int Nsecfilt; // number of catalogs generated + off_t *NAVERAGE; // allocated Averages per catalog + off_t *NMEASURE; // allocated Measures per catalog + int *index; // lookup table catID -> catalog[i] + int *catIDs; // lookup table catID <- catalog[i] + int maxID; // max catID value to date +} CatalogSplitter; + /* global variables set in parameter file */ -char ImageCat[256]; -char ImageTemplate[256]; -char CatTemplate[256]; -char GSCFILE[256]; -char CATDIR[256]; +# define MAX_PATH_LENGTH 1024 +char ImageCat[MAX_PATH_LENGTH]; +char ImageTemplate[MAX_PATH_LENGTH]; +char CatTemplate[MAX_PATH_LENGTH]; +char GSCFILE[MAX_PATH_LENGTH]; +char *CATDIR; char CATMODE[16]; /* raw, mef, split, mysql */ char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ -char CameraConfig[256]; -char SKY_TABLE[256]; +char CameraConfig[MAX_PATH_LENGTH]; +char CAMERA[64]; /* eg, gpc1 */ +char SKY_TABLE[MAX_PATH_LENGTH]; int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ + +int HOST_ID; +char *HOSTDIR; +char *IMAGES; +char *BCATALOG; +ModeType MODE; double MAG_LIM; @@ -54,4 +106,9 @@ double MIN_ERROR; double IMFIT_SYS_SIGMA_LIM; +double CLOUD_TOLERANCE; + +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; int VERBOSE; @@ -79,7 +136,11 @@ int FREEZE_MOSAICS; int USE_GRID; +int KEEP_UBERCAL; char *OUTROOT; +char *UPDATE_CATFORMAT; int PLOTDELAY; int UpdateAverages; + +char *PhotcodeList; int RELPHOT_GRID_X; @@ -101,5 +162,6 @@ double AreaXmin, AreaXmax, AreaYmin, AreaYmax; -int ImagSelect, ImagMin, ImagMax; +int ImagSelect; +double ImagMin, ImagMax; int DophotSelect, DophotValue; @@ -112,5 +174,5 @@ SkyRegion UserPatch; -int UserPatchSelect; +char *UserCatalog; int USE_BASIC_CHECK; @@ -131,7 +193,8 @@ off_t *SelectRefMosaic PROTO((Mosaic **refmosaic, off_t *Nimage)); int args PROTO((int argc, char **argv)); +int args_client PROTO((int argc, char **argv)); int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); void clean_images PROTO((void)); -void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final)); +void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr)); void clean_mosaics PROTO((void)); void clean_stars PROTO((Catalog *catalog, int Ncatalog)); @@ -156,14 +219,18 @@ Coords *getCoords PROTO((off_t meas, int cat)); off_t getImageEntry PROTO((off_t meas, int cat)); -float getMcal PROTO((off_t meas, int cat)); +float getMcal PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); +float getMflat PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); float getMgrid PROTO((off_t meas, int cat)); float getMmos PROTO((off_t meas, int cat)); float getMrel PROTO((Catalog *catalog, off_t meas, int cat)); +short getUbercalDist PROTO((off_t meas, int cat)); Image *getimage PROTO((off_t N)); Image *getimages PROTO((off_t *N, off_t **LineNumber)); -void global_stats PROTO((Catalog *catalog, int Ncatalog)); +ImageSubset *getimages_subset PROTO((off_t *N)); +void global_stats PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void initGrid PROTO((int dX, int dY)); void initGridBins PROTO((Catalog *catalog, int Ncatalog)); void initImageBins PROTO((Catalog *catalog, int Ncatalog, int doImageList)); +void initImagesSubset PROTO((ImageSubset *input, off_t *line_number, off_t N)); void initImages PROTO((Image *input, off_t *LineNumber, off_t N)); void initMosaicBins PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); @@ -172,9 +239,12 @@ void initMrel PROTO((Catalog *catalog, int Ncatalog)); void initialize PROTO((int argc, char **argv)); +void initialize_client PROTO((int argc, char **argv)); void initstats PROTO((char *mode)); int liststats PROTO((double *value, double *dvalue, int N, StatType *stats)); -Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog)); -SkyList *load_images PROTO((FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect)); -Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage)); +Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath)); +Catalog *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog)); + +SkyList *load_images PROTO((FITS_DB *db, char *regionName, SkyRegion *region)); +Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage)); int main PROTO((int argc, char **argv)); @@ -185,27 +255,28 @@ void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); void plot_defaults PROTO((Graphdata *graphdata)); -void plot_grid PROTO((Catalog *catalog)); +void plot_grid PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); void plot_images PROTO((void)); void plot_list PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) OHANA_FORMAT(printf, 6, 7) ); void plot_mosaic_fields PROTO((Catalog *catalog)); void plot_mosaics PROTO((void)); -void plot_scatter PROTO((Catalog *catalog, int Ncatalog)); +void plot_scatter PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void plot_star_coords PROTO((Catalog *catalog, int Ncatalog)); void plot_stars PROTO((Catalog *catalog, int Ncatalog)); -void reload_catalogs PROTO((SkyList *skylist)); +void reload_catalogs PROTO((SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath)); +int reload_catalogs_parallel PROTO((SkyList *sky)); int reload_images PROTO((FITS_DB *db)); -int setExclusions PROTO((Catalog *catalog, int Ncatalog)); -void setMcal PROTO((Catalog *catalog, int Poor)); +int setExclusions PROTO((Catalog *catalog, int Ncatalog, int verbose)); +void setMcal PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); void setMcalFinal PROTO((void)); -int setMcalOutput PROTO((Catalog *catalog, int Ncatalog)); -void setMgrid PROTO((Catalog *catalog)); -int setMmos PROTO((Catalog *catalog, int Poor)); -int setMrel PROTO((Catalog *catalog, int Ncatalog)); -void setMrelFinal PROTO((Catalog *catalog)); -int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int mark)); +int setMcalOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); +void setMgrid PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +int setMmos PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); +int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); +void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr)); int setMave PROTO((Catalog *catalog, int Ncatalog)); void set_ZP PROTO((double ZERO)); int setrefcode PROTO((Image *image, off_t Nimage)); -void skip_measurements PROTO((Catalog *catalog, int pass)); +void skip_measurements PROTO((Catalog *catalog, int pass, FlatCorrectionTable *flatcorr)); void sortA PROTO((double *X, int N)); void sortB PROTO((double *X, double *Y, int N)); @@ -213,12 +284,12 @@ void sortD PROTO((double *X, double *Y, double *Z, int N)); StatType statsImageM PROTO((Catalog *catalog)); -StatType statsImageN PROTO((Catalog *catalog)); +StatType statsImageN PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); StatType statsImageX PROTO((Catalog *catalog)); StatType statsImagedM PROTO((Catalog *catalog)); StatType statsMosaicM PROTO((Catalog *catalog)); -StatType statsMosaicN PROTO((Catalog *catalog)); +StatType statsMosaicN PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); StatType statsMosaicX PROTO((Catalog *catalog)); StatType statsMosaicdM PROTO((Catalog *catalog)); -StatType statsStarN PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode)); +StatType statsStarN PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr)); StatType statsStarS PROTO((Catalog *catalog, int Ncatalog, int Nsec)); StatType statsStarX PROTO((Catalog *catalog, int Ncatalog, int Nsec)); @@ -226,9 +297,12 @@ void wimages PROTO((void)); void write_coords PROTO((Header *header, Coords *coords)); -int relphot_objects (void); +int relphot_objects (int hostID, char *hostpath); void relphot_usage (void); void relphot_help (int argc, char **argv); +void relphot_client_usage (void); +void relphot_client_help (int argc, char **argv); + off_t getImageByID (off_t ID); @@ -236,4 +310,11 @@ int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog); -int populate_tiny_values (Catalog *catalog); -int free_tiny_values (Catalog *catalog); +BrightCatalog *BrightCatalogLoad(char *filename); +int BrightCatalogSave(char *filename, BrightCatalog *catalog); +BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog); +CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt); +int BrightCatalogSplitFree (CatalogSplitter *catalogs); +int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog); + +int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage); +ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage); Index: trunk/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 33651) +++ trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 33651) @@ -0,0 +1,713 @@ +# include "relphot.h" + +// relphot_client reads from the local catalogs and generates a set of bright, subset catalogs +// we are going to save these as a single big FITS file, with only 3 extensions: + +// measure, average, secfilt. these only need to contain the fields relevant to the +// MeasureTiny, AverageTiny, and Secfilt tables, but with cat_id appended to each row so we can +// reassign correctly on read + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// XXX double check the header extname values for each table +// XXX make sure we free things as we can +// XXX make sure we close files as we can +// XXX handle free and close on error return as well +BrightCatalog *BrightCatalogLoad(char *filename) { + + int i, Ncol; + off_t Nrow; + char type[16]; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + header.buffer = NULL; + matrix.buffer = NULL; + ftable.buffer = NULL; + theader.buffer = NULL; + BrightCatalog *catalog = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + goto escape; + } + + ALLOCATE (catalog, BrightCatalog, 1); + catalog->Naverage = 0; + catalog->Nmeasure = 0; + + ftable.header = &theader; + + // *** Measure Data *** + { + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN(dR, "RA_OFF", float); + GET_COLUMN(dD, "DEC_OFF", float); + GET_COLUMN(M, "MAG_SYS", float); + GET_COLUMN(Mcal, "MAG_CAL", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(airmass, "AIRMASS", float); + GET_COLUMN(Xccd, "X_CCD", float); + GET_COLUMN(Yccd, "Y_CCD", float); + GET_COLUMN(dt, "EXPTIME", float); + GET_COLUMN(t, "TIME", int); + GET_COLUMN(averef, "AVE_REF", int); // XXX signed vs unsigned? + GET_COLUMN(imageID, "IMAGE_ID", int); + GET_COLUMN(dbFlags, "DB_FLAGS", int); + GET_COLUMN(photFlags, "PHOT_FLAGS", int); + GET_COLUMN(catID, "CAT_ID", int); + GET_COLUMN(photcode, "PHOTCODE", short); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + MeasureTiny *measure = NULL; + ALLOCATE (measure, MeasureTiny, Nrow); + for (i = 0; i < Nrow; i++) { + measure[i].dR = dR[i]; + measure[i].dD = dD[i]; + measure[i].M = M[i]; + measure[i].Mcal = Mcal[i]; + measure[i].dM = dM[i]; + measure[i].airmass = airmass[i]; + measure[i].Xccd = Xccd[i]; + measure[i].Yccd = Yccd[i]; + measure[i].dt = dt[i]; + measure[i].t = t[i]; + measure[i].averef = averef[i]; + measure[i].imageID = imageID[i]; + measure[i].dbFlags = dbFlags[i]; + measure[i].photFlags = photFlags[i]; + measure[i].catID = catID[i]; + measure[i].photcode = photcode[i]; + } + fprintf (stderr, "loaded data for %lld measure\n", (long long) Nrow); + + free (dR ); + free (dD ); + free (M ); + free (Mcal ); + free (dM ); + free (airmass ); + free (Xccd ); + free (Yccd ); + free (dt ); + free (t ); + free (averef ); + free (imageID ); + free (dbFlags ); + free (photFlags ); + free (catID ); + free (photcode); + + catalog->measure = measure; + catalog->Nmeasure = Nrow; + } + + /*** load Average values ***/ + { + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN(R, "RA", double); + GET_COLUMN(D, "DEC", double); + GET_COLUMN(Nmeasure, "NMEAS", int); + GET_COLUMN(measureOffset, "MEASURE_OFF", int); + GET_COLUMN(flags, "FLAGS", int); + GET_COLUMN(catID, "CAT_ID", int); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + AverageTiny *average = NULL; + ALLOCATE (average, AverageTiny, Nrow); + for (i = 0; i < Nrow; i++) { + average[i].R = R[i]; + average[i].D = D[i]; + average[i].Nmeasure = Nmeasure[i]; + average[i].measureOffset = measureOffset[i]; + average[i].flags = flags[i]; + average[i].catID = catID[i]; + } + fprintf (stderr, "loaded data for %lld average\n", (long long) Nrow); + + free (R ); + free (D ); + free (Nmeasure ); + free (measureOffset ); + free (flags ); + free (catID ); + + catalog->average = average; + catalog->Naverage = Nrow; + } + + /*** load Secfilt values ***/ + { + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN(M, "MAG", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(Xm, "MAG_CHI", float); + GET_COLUMN(flags, "FLAGS", int); + GET_COLUMN(Ncode, "NCODE", short); + GET_COLUMN(Nused, "NUSED", short); + GET_COLUMN(M_20, "MAG_20", short); + GET_COLUMN(M_80, "MAG_80", short); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + SecFilt *secfilt = NULL; + ALLOCATE (secfilt, SecFilt, Nrow); + for (i = 0; i < Nrow; i++) { + secfilt[i].M = M[i]; + secfilt[i].dM = dM[i]; + secfilt[i].Xm = Xm[i]; + secfilt[i].flags = flags[i]; + secfilt[i].Ncode = Ncode[i]; + secfilt[i].Nused = Nused[i]; + secfilt[i].M_20 = M_20[i]; + secfilt[i].M_80 = M_80[i]; + } + fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow); + + free (M ); + free (dM ); + free (Xm ); + free (flags); + free (Ncode); + free (Nused); + free (M_20 ); + free (M_80 ); + catalog->secfilt = secfilt; + // assert Nsecfilt * Naverage = Nrow? + } + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + return catalog; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + if (catalog) free (catalog); + + fclose (f); + return NULL; +} + +// we are passed a BrightCatalog structure, write it to a FITS table (3 ext) +int BrightCatalogSave(char *filename, BrightCatalog *catalog) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename); + return FALSE; + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + /*** MeasureTiny ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY"); + + gfits_define_bintable_column (&theader, "E", "RA_OFF", "ra offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DEC_OFF", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "magnitude (err)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "AIRMASS", "airmass", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "X_CCD", "ccd x coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Y_CCD", "ccd y coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "EXPTIME", "-2.5 * log (exposure time)", "sec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "TIME", "time of exp", "sec", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "AVE_REF", "pointer to average table", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode", NULL, 1.0, 1.0*0x80); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + float *dR ; ALLOCATE (dR , float, catalog->Nmeasure); + float *dD ; ALLOCATE (dD , float, catalog->Nmeasure); + float *M ; ALLOCATE (M , float, catalog->Nmeasure); + float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); + float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); + float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); + float *Xccd ; ALLOCATE (Xccd , float, catalog->Nmeasure); + float *Yccd ; ALLOCATE (Yccd , float, catalog->Nmeasure); + float *dt ; ALLOCATE (dt , float, catalog->Nmeasure); + int *t ; ALLOCATE (t , int , catalog->Nmeasure); + int *averef ; ALLOCATE (averef , int , catalog->Nmeasure); + int *imageID ; ALLOCATE (imageID , int , catalog->Nmeasure); + int *dbFlags ; ALLOCATE (dbFlags , int , catalog->Nmeasure); + int *photFlags ; ALLOCATE (photFlags, int , catalog->Nmeasure); + int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); + short *photcode ; ALLOCATE (photcode , short, catalog->Nmeasure); + + // assign the storage arrays + MeasureTiny *measure = catalog->measure; + for (i = 0; i < catalog->Nmeasure; i++) { + dR[i] = measure[i].dR ; + dD[i] = measure[i].dD ; + M[i] = measure[i].M ; + Mcal[i] = measure[i].Mcal ; + dM[i] = measure[i].dM ; + airmass[i] = measure[i].airmass ; + Xccd[i] = measure[i].Xccd ; + Yccd[i] = measure[i].Yccd ; + dt[i] = measure[i].dt ; + t[i] = measure[i].t ; + averef[i] = measure[i].averef ; + catID[i] = measure[i].catID ; + imageID[i] = measure[i].imageID ; + dbFlags[i] = measure[i].dbFlags ; + photFlags[i]= measure[i].photFlags; + photcode[i] = measure[i].photcode ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA_OFF", dR, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DEC_OFF", dD, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "X_CCD", Xccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "Y_CCD", Yccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "EXPTIME", dt, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "TIME", t, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AVE_REF", averef, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS", dbFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOT_FLAGS", photFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, catalog->Nmeasure); + + free (dR ); + free (dD ); + free (M ); + free (Mcal ); + free (dM ); + free (airmass ); + free (Xccd ); + free (Yccd ); + free (dt ); + free (t ); + free (averef ); + free (imageID ); + free (dbFlags ); + free (photFlags ); + free (catID ); + free (photcode); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** AverageTiny ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "AVERAGE_TINY"); + + gfits_define_bintable_column (&theader, "D", "RA", "ra (J2000)", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "dec (J2000)", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NMEAS", "number of measures", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "MEASURE_OFF", "index to measurements", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ref", NULL, 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + double *R ; ALLOCATE (R, double, catalog->Naverage); + double *D ; ALLOCATE (D, double, catalog->Naverage); + int *Nmeasure ; ALLOCATE (Nmeasure, int, catalog->Naverage); + int *measureOffset ; ALLOCATE (measureOffset, int, catalog->Naverage); + int *flags ; ALLOCATE (flags, int, catalog->Naverage); + int *catID ; ALLOCATE (catID, int, catalog->Naverage); + + // assign the storage arrays + AverageTiny *average = catalog->average; + for (i = 0; i < catalog->Naverage; i++) { + R[i] = average[i].R ; + D[i] = average[i].D ; + Nmeasure[i] = average[i].Nmeasure; + measureOffset[i] = average[i].measureOffset; + flags[i] = average[i].flags; + catID[i] = average[i].catID; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "NMEAS", Nmeasure, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "MEASURE_OFF", measureOffset, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Naverage); + + free (R ); + free (D ); + free (Nmeasure ); + free (measureOffset ); + free (flags ); + free (catID ); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** SecFilt ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "SECFILT"); + + gfits_define_bintable_column (&theader, "E", "MAG", "ra offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CHI", "magnitude (sys)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NCODE", "magnitude (err)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NUSED", "airmass", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "MAG_20", "ccd x coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "MAG_80", "ccd y coord", "pix", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // Nsecfilt is number of average filters, Nsec is total number of values + int Nsecfilt = GetPhotcodeNsecfilt(); + int Nsec = Nsecfilt * catalog->Naverage; + + // create intermediate storage arrays + float *M ; ALLOCATE (M , float, Nsec); + float *dM ; ALLOCATE (dM , float, Nsec); + float *Xm ; ALLOCATE (Xm , float, Nsec); + int *flags ; ALLOCATE (flags , int, Nsec); + short *Ncode ; ALLOCATE (Ncode , short, Nsec); + short *Nused ; ALLOCATE (Nused , short, Nsec); + short *M_20 ; ALLOCATE (M_20 , short, Nsec); + short *M_80 ; ALLOCATE (M_80 , short, Nsec); + + // assign the storage arrays + SecFilt *secfilt = catalog->secfilt; + for (i = 0; i < Nsec; i++) { + M [i] = secfilt[i]. M ; + dM [i] = secfilt[i]. dM ; + Xm [i] = secfilt[i]. Xm ; + flags [i] = secfilt[i]. flags ; + Ncode [i] = secfilt[i]. Ncode ; + Nused [i] = secfilt[i]. Nused ; + M_20 [i] = secfilt[i]. M_20 ; + M_80 [i] = secfilt[i]. M_80 ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Xm , Nsec); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nsec); + gfits_set_bintable_column (&theader, &ftable, "NCODE", Ncode, Nsec); + gfits_set_bintable_column (&theader, &ftable, "NUSED", Nused, Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_20", M_20 , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_80", M_80 , Nsec); + + free (M ); + free (dM ); + free (Xm ); + free (flags ); + free (Ncode ); + free (Nused ); + free (M_20 ); + free (M_80 ); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + return TRUE; +} + +// merge a list of catalogs into a single BrightCatalog array set +BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog) { + + off_t i, j, k; + + BrightCatalog *bcatalog = NULL; + ALLOCATE (bcatalog, BrightCatalog, 1); + + int Nmeas = 0; + int Naves = 0; + for (i = 0; i < Ncatalog; i++) { + Nmeas += catalog[i].Nmeasure; + Naves += catalog[i].Naverage; + } + + // XXX this prevents different catalogs from having different Nsecfilt values + int Nsecfilt = GetPhotcodeNsecfilt(); + + ALLOCATE (bcatalog[0].measure, MeasureTiny, Nmeas); + ALLOCATE (bcatalog[0].average, AverageTiny, Naves); + ALLOCATE (bcatalog[0].secfilt, SecFilt, Naves*Nsecfilt); + + int Nm = 0; + int Na = 0; + for (i = 0; i < Ncatalog; i++) { + if (!catalog[i].Naverage) continue; + for (j = 0; j < catalog[i].Naverage; j++) { + bcatalog[0].average[Na] = catalog[i].averageT[j]; + for (k = 0; k < Nsecfilt; k++) { + bcatalog[0].secfilt[Nsecfilt*Na + k] = catalog[i].secfilt[Nsecfilt*j + k]; + } + Na++; + assert (Na <= Naves); + } + for (j = 0; j < catalog[i].Nmeasure; j++) { + bcatalog[0].measure[Nm] = catalog[i].measureT[j]; + Nm++; + assert (Nm <= Nmeas); + } + } + bcatalog->Naverage = Na; + bcatalog->Nmeasure = Nm; + return bcatalog; +} + +// distribute a bright catalog across separate catalogs +CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt) { + + int i; + + CatalogSplitter *catalogs = NULL; + + ALLOCATE (catalogs, CatalogSplitter, 1); + + // as we see new BrightCatalogs, we will update maxID and extend this array + catalogs->maxID = 0; + ALLOCATE (catalogs->index, int, catalogs->maxID + 1); + for (i = 0; i <= catalogs->maxID; i++) catalogs->index[i] = -1; + + catalogs->Nsecfilt = Nsecfilt; + + catalogs->Ncatalog = 0; + catalogs->NCATALOG = 16; + catalogs->catalog = NULL; + ALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG); + + ALLOCATE (catalogs->catIDs, int, catalogs->NCATALOG); + ALLOCATE (catalogs->NAVERAGE, off_t, catalogs->NCATALOG); + ALLOCATE (catalogs->NMEASURE, off_t, catalogs->NCATALOG); + + for (i = 0; i < catalogs->NCATALOG; i++) { + dvo_catalog_init (&catalogs->catalog[i], TRUE); + catalogs->catIDs[i] = 0; + catalogs->NAVERAGE[i] = 256; + catalogs->NMEASURE[i] = 256; + catalogs->catalog[i].Naverage = 0; + catalogs->catalog[i].Nmeasure = 0; + ALLOCATE (catalogs->catalog[i].averageT, AverageTiny, catalogs->NAVERAGE[i]); + ALLOCATE (catalogs->catalog[i].measureT, MeasureTiny, catalogs->NMEASURE[i]); + ALLOCATE (catalogs->catalog[i].secfilt, SecFilt, catalogs->NAVERAGE[i]*Nsecfilt); + } + return catalogs; +} + +// distribute a bright catalog across separate catalogs +int BrightCatalogSplitFree (CatalogSplitter *catalogs) { + + // XXX don't free the catalogs : + // for (i = 0; i < NCATALOG; i++) { + // free (catalogs->catalog[i].averageT); + // free (catalogs->catalog[i].measureT); + // free (catalogs->catalog[i].secfilt); + // } + // free (catalogs->catalog); + + free (catalogs->catIDs); + free (catalogs->NAVERAGE); + free (catalogs->NMEASURE); + free (catalogs->index); + free (catalogs); + return TRUE; +} + +// distribute a bright catalog across separate catalogs +int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog) { + + int i; + + int D_NCATALOG = 16; + + int Nsecfilt = catalogs->Nsecfilt; + + // find the max value of catID in this BrightCatalog + int catIDmax = 0; + for (i = 0; i < bcatalog->Naverage; i++) { + catIDmax = MAX(catIDmax, bcatalog->average[i].catID); + } + // XXX validate the measure ID range here + + int maxIDold = catalogs->maxID; + catalogs->maxID = MAX (maxIDold, catIDmax); + + // extend the index array and init + REALLOCATE (catalogs->index, int, catalogs->maxID + 1); + for (i = maxIDold + 1; i <= catalogs->maxID; i++) catalogs->index[i] = -1; + + // identify the new catID values + for (i = 0; i < bcatalog->Naverage; i++) { + int catID = bcatalog->average[i].catID; + assert (catID > 0); + assert (catID < catalogs->maxID + 1); + int idx = catalogs->index[catID]; + if (idx != -1) continue; // already have seen this one + + // a new catID value + int Ncat = catalogs->Ncatalog; // the next available slot + catalogs->catIDs[Ncat] = catID; + assert (Ncat >= 0); + assert (Ncat < catalogs->NCATALOG); + + catalogs->catalog[Ncat].Nsecfilt = Nsecfilt; + catalogs->catalog[Ncat].catID = catID; + + catalogs->index[catID] = Ncat; + assert (catID > 0); + assert (catID < catalogs->maxID + 1); + + catalogs->Ncatalog ++; + + if (catalogs->Ncatalog >= catalogs->NCATALOG) { + catalogs->NCATALOG += D_NCATALOG; + + // fprintf (stderr, "realloc catalogs->catalog: old: %llx ", (long long) catalogs->catalog); + REALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->catalog, (long long) (catalogs->catalog + sizeof(Catalog)*catalogs->NCATALOG)); + + // fprintf (stderr, "realloc catalogs->NAVERAGE: old: %llx ", (long long) catalogs->NAVERAGE); + REALLOCATE (catalogs->NAVERAGE, off_t, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->NAVERAGE, (long long) (catalogs->NAVERAGE + sizeof(off_t)*catalogs->NCATALOG)); + + // fprintf (stderr, "realloc catalogs->NMEASURE: old: %llx ", (long long) catalogs->NMEASURE); + REALLOCATE (catalogs->NMEASURE, off_t, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->NMEASURE, (long long) (catalogs->NMEASURE + sizeof(off_t)*catalogs->NCATALOG)); + + REALLOCATE (catalogs->catIDs, int, catalogs->NCATALOG); + + int j; + for (j = catalogs->NCATALOG - D_NCATALOG; j < catalogs->NCATALOG; j++) { + dvo_catalog_init (&catalogs->catalog[j], TRUE); + catalogs->catIDs[j] = 0; + catalogs->NAVERAGE[j] = 256; + catalogs->NMEASURE[j] = 256; + catalogs->catalog[j].Naverage = 0; + catalogs->catalog[j].Nmeasure = 0; + ALLOCATE (catalogs->catalog[j].averageT, AverageTiny, catalogs->NAVERAGE[j]); + ALLOCATE (catalogs->catalog[j].measureT, MeasureTiny, catalogs->NMEASURE[j]); + ALLOCATE (catalogs->catalog[j].secfilt, SecFilt, catalogs->NAVERAGE[j]*Nsecfilt); + } + D_NCATALOG = MAX(2000, 2*D_NCATALOG); + } + } + + // each bcatalog (from each host) has a different set of catID ranges + // they also (probably) have contiguous ranges. But, it is a bit tricky to be sure I + // can use that info, so perhaps ignore it + + // assign the averages to the corresponding catalog + for (i = 0; i < bcatalog->Naverage; i++) { + int ID = bcatalog->average[i].catID; + int Nc = catalogs->index[ID]; + assert (Nc > -1); + assert (Nc < catalogs->NCATALOG); + + int Na = catalogs->catalog[Nc].Naverage; + catalogs->catalog[Nc].averageT[Na] = bcatalog->average[i]; + + // secfilt entries are grouped in blocks for each average + int k; + for (k = 0; k < Nsecfilt; k++) { + catalogs->catalog[Nc].secfilt[Na*Nsecfilt + k] = bcatalog->secfilt[i*Nsecfilt + k]; + } + + catalogs->catalog[Nc].Naverage ++; + if (catalogs->catalog[Nc].Naverage >= catalogs->NAVERAGE[Nc]) { + catalogs->NAVERAGE[Nc] += MAX(catalogs->NAVERAGE[Nc], 1024); + REALLOCATE (catalogs->catalog[Nc].averageT, AverageTiny, catalogs->NAVERAGE[Nc]); + REALLOCATE (catalogs->catalog[Nc].secfilt, SecFilt, catalogs->NAVERAGE[Nc]*Nsecfilt); + } + } + + // assign the measures to the corresponding catalog + // XXX what about averef and related links? Do I need them? + for (i = 0; i < bcatalog->Nmeasure; i++) { + int ID = bcatalog->measure[i].catID; + int Nc = catalogs->index[ID]; + assert (Nc > -1); + assert (Nc < catalogs->NCATALOG); + int Na = catalogs->catalog[Nc].Nmeasure; + catalogs->catalog[Nc].measureT[Na] = bcatalog->measure[i]; + catalogs->catalog[Nc].Nmeasure ++; + if (catalogs->catalog[Nc].Nmeasure >= catalogs->NMEASURE[Nc]) { + catalogs->NMEASURE[Nc] += MAX(catalogs->NMEASURE[Nc], 4096); + REALLOCATE (catalogs->catalog[Nc].measureT, MeasureTiny, catalogs->NMEASURE[Nc]); + } + } + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relphot/src/ConfigInit.c (revision 33650) +++ trunk/Ohana/src/relphot/src/ConfigInit.c (revision 33651) @@ -33,6 +33,14 @@ GetConfig (config, "IMAGE_GOOD_FRACTION", "%lf", 0, &IMAGE_GOOD_FRACTION); - GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); - GetConfig (config, "CATDIR", "%s", 0, CATDIR); + // force CATDIR to be absolute (so parallel mode will work) + char *tmpcatdir = NULL; + ALLOCATE (tmpcatdir, char, MAX_PATH_LENGTH); + GetConfig (config, "CATDIR", "%s", 0, tmpcatdir); + CATDIR = abspath (tmpcatdir, MAX_PATH_LENGTH); + free (tmpcatdir); + + GetConfig (config, "CAMERA", "%s", 0, CAMERA); + + GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); ScanConfig (config, "CATMODE", "%s", 0, CATMODE); ScanConfig (config, "CATFORMAT", "%s", 0, CATFORMAT); Index: trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- trunk/Ohana/src/relphot/src/GridOps.c (revision 33650) +++ trunk/Ohana/src/relphot/src/GridOps.c (revision 33651) @@ -284,5 +284,5 @@ /* direct (non-iterative) solution for Mgrid values for all grid bins */ -void setMgridDirect (Catalog *catalog, int Ncatalog) { +void setMgridDirect (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { int **gotstar, **gridmeas; @@ -399,5 +399,5 @@ // skip images marked as BAD - Mcal = getMcal (m, c); + Mcal = getMcal (m, c, flatcorr, catalog); if (isnan(Mcal)) { Ncal ++; @@ -522,5 +522,5 @@ /* determine Mgrid values for all grid bins */ -void setMgrid (Catalog *catalog) { +void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) { int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit; @@ -554,5 +554,5 @@ continue; } - Mcal = getMcal (m, c); + Mcal = getMcal (m, c, flatcorr, catalog); if (isnan(Mcal)) { Ncal ++; @@ -609,5 +609,5 @@ } -void plot_grid (Catalog *catalog) { +void plot_grid (Catalog *catalog, FlatCorrectionTable *flatcorr) { int i, j, m, c, n, N, Narea; @@ -641,5 +641,5 @@ continue; } - Mcal = getMcal (m, c); + Mcal = getMcal (m, c, flatcorr, catalog); if (isnan(Mcal)) continue; Mmos = getMmos (m, c); Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 33650) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 33651) @@ -1,12 +1,32 @@ # include "relphot.h" - -static off_t **bin; // link from catalog,measure to image -static off_t **clist; // catalog which supplied measurement on image -static off_t **mlist; // measure reference for measurement on image -static off_t *Nlist; // number of measurements on image -static off_t *NLIST; // allocated number of measurements on image - -static Image *image; // list of available images +void plot_setMcal (double *list, int Npts, StatType *stats, float clouds); + +// the MeasureToImage, ImageToCatalog, and ImageToMeasure arrays are a substantial part of the memory footprint. +// total data volume is currently: +// a) Naverage*sizeof(AverageTiny) [averages] : 32 +// b) Naverage*sizeof(Secfilt)*Nsecfilt [secfilt values] : 8 * 32 +// b) Nmeasure*sizeof(Measure) [measurements] : 72 +// c) Nmeasure*sizeof(IDX_T)*3 [image idx] : 3 * 16 +// d) Nmeasure*sizeof(IDX_T)*3 [mosaic idx] : 3 * 16 +// e) Nimage*sizeof(Image) [image data] : 360 + +// for 3pi analysis in 2012, we have Nmeasure ~ 20 x Naverage +// so, each measurement costs ~14.5B (Ave + Sec) + 72B (Meas) + 96B (idx)! +// we are using off_t (64bit) to avoid the 32bit limit of an int, but +// if we really had 2^31 measurements in a single analysis, we would be using > 350GB of ram... +// until we reach that point, it is sort of silly to use IDX_T = off_t here. + +// in fact, we are safer than this, because the number of detections per table is rarely so large. + +// with IDX_T = int, each measurement costs ~14.5B (Ave + Sec) + 72B (Meas) + 48B (idx)! + +// elsewhere (not in relphot_images), we need to use off_t because a single catalog + +// we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage) static off_t Nimage; // number of available images +static Image *image; // array of available images + +// if we read only a subset of the rows from the Image FITS, LineNumber tells us to which row +// each image belongs static off_t *LineNumber; // match of subset to full image table @@ -15,4 +35,21 @@ static off_t *imageIdx; // list of index for image IDs +// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas]) +// each image has N_onImage[ImageIndex] measurements +static off_t *N_onImage; // actual number of measurements on image +static off_t *N_ONIMAGE; // allocated number of measurements on image + +// relationships between the measure,catalog set and the images: +static IDX_T **MeasureToImage; // image index from measure,catalog : MeasureToImage[cat][meas] = ImageIndex +static IDX_T **ImageToCatalog; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) +static IDX_T **ImageToMeasure; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex]) + +// MeasureToImage was 'bin' +// ImageToCatalog was 'clist' +// ImageToMeasure was 'mlist' + +// N_onImage was 'Nlist' +// N_ONIMAGE was 'NLIST' + Image *getimages (off_t *N, off_t **line_number) { @@ -42,4 +79,55 @@ } llsortpair (imageIDs, imageIdx, Nimage); +} + +void initImagesSubset (ImageSubset *input, off_t *line_number, off_t N) { + + off_t i; + + // we have been given an ImageSubset array, containing a reduced set of image fields + // create full a Image array and save the needed values + ALLOCATE (image, Image, N); + for (i = 0; i < N; i++) { + image[i].imageID = input[i].imageID ; + image[i].photom_map_id = input[i].photom_map_id; + image[i].flags = input[i].flags ; + image[i].Mcal = input[i].Mcal ; + image[i].dMcal = input[i].dMcal ; + image[i].tzero = input[i].tzero ; + image[i].trate = input[i].trate ; + } + LineNumber = line_number; + Nimage = N; + + ALLOCATE (imageIDs, off_t, Nimage); + ALLOCATE (imageIdx, off_t, Nimage); + + for (i = 0; i < Nimage; i++) { + imageIdx[i] = i; + imageIDs[i] = image[i].imageID; + } + llsortpair (imageIDs, imageIdx, Nimage); +} + +ImageSubset *getimages_subset (off_t *N) { + + *N = Nimage; + + off_t i; + + // we have been given an ImageSubset array, containing a reduced set of image fields + // create full a Image array and save the needed values + ImageSubset *subset = NULL; + ALLOCATE (subset, ImageSubset, Nimage); + for (i = 0; i < Nimage; i++) { + subset[i].imageID = image[i].imageID ; + subset[i].photom_map_id = image[i].photom_map_id; + subset[i].flags = image[i].flags ; + subset[i].Mcal = image[i].Mcal ; + subset[i].dMcal = image[i].dMcal ; + subset[i].tzero = image[i].tzero ; + subset[i].trate = image[i].trate ; + } + return subset; } @@ -71,23 +159,23 @@ void initImageBins (Catalog *catalog, int Ncatalog, int doImageList) { - off_t i, j; - - ALLOCATE (bin, off_t *, Ncatalog); + IDX_T i, j; + + ALLOCATE (MeasureToImage, IDX_T *, Ncatalog); for (i = 0; i < Ncatalog; i++) { - ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1)); - for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1; + ALLOCATE (MeasureToImage[i], IDX_T, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToImage[i][j] = -1; } if (doImageList) { - ALLOCATE (Nlist, off_t, Nimage); - ALLOCATE (NLIST, off_t, Nimage); - ALLOCATE (clist, off_t *, Nimage); - ALLOCATE (mlist, off_t *, Nimage); + ALLOCATE (N_onImage, off_t, Nimage); + ALLOCATE (N_ONIMAGE, off_t, Nimage); + ALLOCATE (ImageToCatalog, IDX_T *, Nimage); + ALLOCATE (ImageToMeasure, IDX_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]); + N_onImage[i] = 0; + N_ONIMAGE[i] = 30; + ALLOCATE (ImageToCatalog[i], IDX_T, N_ONIMAGE[i]); + ALLOCATE (ImageToMeasure[i], IDX_T, N_ONIMAGE[i]); } } @@ -99,21 +187,21 @@ for (i = 0; i < Ncatalog; i++) { - free (bin[i]); - } - free (bin); + free (MeasureToImage[i]); + } + free (MeasureToImage); 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 */ + free (ImageToCatalog[i]); + free (ImageToMeasure[i]); + } + free (ImageToCatalog); + free (ImageToMeasure); + free (N_onImage); + free (N_ONIMAGE); + } +} + +/* select all images equivalent to the active photcode set */ void findImages (Catalog *catalog, int Ncatalog, int doImageList) { @@ -210,18 +298,18 @@ // index for (catalog, measure) -> image - bin[cat][meas] = idx; + MeasureToImage[cat][meas] = idx; if (doImageList) { // index for image, Nentry -> catalog - clist[idx][Nlist[idx]] = cat; + ImageToCatalog[idx][N_onImage[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]); + ImageToMeasure[idx][N_onImage[idx]] = meas; + N_onImage[idx] ++; + + if (N_onImage[idx] == N_ONIMAGE[idx]) { + N_ONIMAGE[idx] += 30; + REALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]); + REALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]); } } @@ -234,26 +322,66 @@ off_t i; - i = bin[cat][meas]; + i = MeasureToImage[cat][meas]; return (i); } -float getMcal (off_t meas, int cat) { +// returns image.Mcal - ff(x,y) +float getMcal (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) { off_t i; - float value; - - i = bin[cat][meas]; + float value, offset; + + i = MeasureToImage[cat][meas]; if (i == -1) return (NAN); if (image[i].flags & IMAGE_BAD) return (NAN); value = image[i].Mcal; + offset = 0.0; + + // to do this, I need to pass in the catalog and flatcorr pointers + int flat_id = image[i].photom_map_id; + if (flat_id) { + offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd); + } + value -= offset; + return (value); } +// returns image.Mcal - ff(x,y) +short getUbercalDist (off_t meas, int cat) { + + off_t i; + short distance; + + i = MeasureToImage[cat][meas]; + if (i == -1) return (1000); + + if (image[i].flags & IMAGE_BAD) return (1000); + distance = image[i].ubercalDist; // was dummy3 in structure + return (distance); +} + +float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) { + + off_t i = MeasureToImage[cat][meas]; + if (i == -1) return (NAN); + + float offset = 0.0; + + // to do this, I need to pass in the catalog and flatcorr pointers + int flat_id = image[i].photom_map_id; + if (flat_id) { + offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd); + } + + return (offset); +} + Coords *getCoords (off_t meas, int cat) { off_t i; - i = bin[cat][meas]; + i = MeasureToImage[cat][meas]; if (i == -1) return (NULL); return (&image[i].coords); @@ -261,32 +389,34 @@ /* determine Mcal values for all images */ -void setMcal (Catalog *catalog, int PoorImages) { - - off_t i, j, m, c, n, N, Nmax, mark, bad; - int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys, Nbright; - float Msys, Mrel, Mmos, Mgrid, McalBright, McalBright2; - double *list, *dlist; +void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + off_t i, j, m, c, n, N, Nmax; + int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys, Nbright; + float Msys, Mrel, Mmos, Mgrid, Mflat; + double *list, *dlist, *Mlist, *dMlist; StatType stats; if (FREEZE_IMAGES) return; + fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); + int Nsecfilt = GetPhotcodeNsecfilt (); if (PoorImages) { - IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + // XXX use bad stars and measurements for PoorImages? or not? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + IMAGE_BAD = 0; } Nmax = 0; for (i = 0; i < Nimage; i++) { - Nmax = MAX (Nmax, Nlist[i]); + Nmax = MAX (Nmax, N_onImage[i]); } ALLOCATE (list, double, Nmax); ALLOCATE (dlist, double, Nmax); + ALLOCATE (Mlist, double, Nmax); + ALLOCATE (dMlist, double, Nmax); Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0; - - // counters to measure the bright-end scatter - McalBright = McalBright2 = 0.0; - Nbright = 0; for (i = 0; i < Nimage; i++) { @@ -298,9 +428,21 @@ } + // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0: + // we optionally do not recalibrate images with UBERCAL zero points + if (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL) { + image[i].ubercalDist = 0; // was dummy3 + if (KEEP_UBERCAL) continue; + } + + int minUbercalDist = 1000; + + // number of stars to measure the bright-end scatter + Nbright = 0; + N = 0; - for (j = 0; j < Nlist[i]; j++) { + for (j = 0; j < N_onImage[i]; j++) { - m = mlist[i][j]; - c = clist[i][j]; + m = ImageToMeasure[i][j]; + c = ImageToCatalog[i][j]; if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { @@ -324,4 +466,11 @@ } + // image.Mcal is not supposed to include the flat-field correction, so we need to + // apply that offset as well here for this image (in other words, each detection is + // being compared to the model, excluding the zero point, Mcal. The model includes + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) + + Mflat = getMflat (m, c, flatcorr, catalog); + n = catalog[c].measureT[m].averef; Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); @@ -330,18 +479,31 @@ continue; } - list[N] = Msys - Mrel - Mmos - Mgrid; + + PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); + if (!code) goto skip; + if (code->equiv < 1) goto skip; + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec == -1) goto skip; + minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist); + + // fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f %6.3f %6.3f %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat); + + skip: + list[N] = Msys - Mrel - Mmos - Mgrid + Mflat; dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { - McalBright += list[N]; - McalBright2 += SQ(list[N]); + Mlist[Nbright] = list[N]; + dMlist[Nbright] = dlist[N]; Nbright ++; } N++; } - /* Nlist[i] is all measurements, N is good measurements */ + /* N_onImage[i] is all measurements, N is good measurements */ + + // if (VERBOSE2) fprintf (stderr, "meas skipped: (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nbad, Nmos, Ngrid, Nrel, Nsys); /* too few good measurements or too many bad measurements */ if (!PoorImages) { - mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]); + mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onImage[i]); if (mark) { image[i].flags |= ID_IMAGE_PHOTOM_FEW; @@ -351,20 +513,29 @@ } } - + liststats (list, dlist, N, &stats); - - float CLOUD_TOLERANCE = 0.01; image[i].Mcal = stats.mean; image[i].dMcal = stats.error; - image[i].dMagSys = stats.sigma; image[i].nFitPhotom = N; image[i].Xm = 100.0*log10(stats.chisq); + plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); + + // bright end scatter + liststats (Mlist, dMlist, Nbright, &stats); + image[i].dMagSys = stats.sigma; + if (image[i].Mcal < -CLOUD_TOLERANCE) { image[i].Mcal = 0.0; } + + // minUbercalDist calculated here is the min value for any star owned by this image + // since this particular image is tied to that star, bump its distance by 1 + image[i].ubercalDist = minUbercalDist + 1; } free (list); free (dlist); + free (Mlist); + free (dMlist); fprintf (stderr, "%d images marked having too few measurements (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys); @@ -395,7 +566,8 @@ ALLOCATE (dlist, double, Nimage); + // measure stats for Mcal and dMcal for (i = N = 0; i < Nimage; i++) { if (image[i].flags & IMAGE_BAD) continue; - mlist[N] = fabs (image[i].Mcal); + mlist[N] = image[i].Mcal; slist[N] = image[i].dMcal; dlist[N] = 1; @@ -412,4 +584,7 @@ Nmark = 0; for (i = 0; i < Nimage; i++) { + // if we are keeping ubercal sacrosanct, then we should not be allowed to break them... + if (KEEP_UBERCAL && (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + mark = FALSE; image[i].flags &= ~ID_IMAGE_PHOTOM_POOR; @@ -430,4 +605,36 @@ } +void plot_setMcal (double *list, int Npts, StatType *stats, float clouds) { + + off_t i; + double *xlist; + Graphdata graphdata; + + if (!PLOTSTUFF) return; + + ALLOCATE (xlist, double, Npts); + + /**** Mcal vs seq ****/ + float minMcal = +100.0; + float maxMcal = -100.0; + for (i = 0; i < Npts; i++) { + xlist[i] = i; + minMcal = MIN (list[i], minMcal); + maxMcal = MAX (list[i], maxMcal); + } + + float McalRange = MAX(1.2*(maxMcal - minMcal), 0.21); + float McalCenter = 0.5*(maxMcal + minMcal); + + plot_defaults (&graphdata); + graphdata.xmin = -1; + graphdata.xmax = Npts + 1; + graphdata.ymin = McalCenter - 0.5*McalRange; + graphdata.ymax = McalCenter + 0.5*McalRange; + plot_list (&graphdata, xlist, list, Npts, "sequence vs Mcal", "%s.seq.png", OUTROOT); + + free (xlist); +} + void plot_images () { @@ -443,15 +650,43 @@ /**** dMcal vs airmass ****/ + float minAirmass = 1000.0; + float maxAirmass = 0.0; + float minMcal = +100.0; + float maxMcal = -100.0; + float mindMcal = +100.0; + float maxdMcal = -100.0; for (i = 0; i < Nimage; i++) { Mlist[i] = image[i].Mcal; dlist[i] = image[i].dMcal; xlist[i] = image[i].secz; - } + minAirmass = MIN (image[i].secz, minAirmass); + maxAirmass = MAX (image[i].secz, maxAirmass); + minMcal = MIN (image[i].Mcal, minMcal); + maxMcal = MAX (image[i].Mcal, maxMcal); + mindMcal = MIN (image[i].dMcal, mindMcal); + maxdMcal = MAX (image[i].dMcal, maxdMcal); + } + + float AirmassRange = MAX(1.2*(maxAirmass - minAirmass), 0.25); + float AirmassCenter = 0.5*(maxAirmass + minAirmass); + + float McalRange = MAX(1.2*(maxMcal - minMcal), 0.21); + float McalCenter = 0.5*(maxMcal + minMcal); + + float dMcalRange = MAX(1.2*(maxdMcal - mindMcal), 0.21); + float dMcalCenter = 0.5*(maxdMcal + mindMcal); plot_defaults (&graphdata); - graphdata.ymin = PlotdMmin; - graphdata.ymax = PlotdMmax; + graphdata.xmin = AirmassCenter - 0.5*AirmassRange; + graphdata.xmax = AirmassCenter + 0.5*AirmassRange; + graphdata.ymin = McalCenter - 0.5*McalRange; + graphdata.ymax = McalCenter + 0.5*McalRange; plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "%s.airmass.png", OUTROOT); + plot_defaults (&graphdata); + graphdata.xmin = McalCenter - 0.5*McalRange; + graphdata.xmax = McalCenter + 0.5*McalRange; + graphdata.ymin = dMcalCenter - 0.5*dMcalRange; + graphdata.ymax = dMcalCenter + 0.5*dMcalRange; plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", "%s.Mcal.dMcal.png", OUTROOT); @@ -471,4 +706,6 @@ plot_defaults (&graphdata); graphdata.style = 1; + graphdata.xmin = dMcalCenter - 0.5*dMcalRange; + graphdata.xmax = dMcalCenter + 0.5*dMcalRange; plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT); @@ -478,5 +715,5 @@ } -StatType statsImageN (Catalog *catalog) { +StatType statsImageN (Catalog *catalog, FlatCorrectionTable *flatcorr) { off_t i, j, m, c, n, N; @@ -496,10 +733,10 @@ N = 0; - for (j = 0; j < Nlist[i]; j++) { - - m = mlist[i][j]; - c = clist[i][j]; - - Mcal = getMcal (m, c); + for (j = 0; j < N_onImage[i]; j++) { + + m = ImageToMeasure[i][j]; + c = ImageToCatalog[i][j]; + + Mcal = getMcal (m, c, flatcorr, catalog); if (isnan(Mcal)) continue; Mmos = getMmos (m, c); Index: trunk/Ohana/src/relphot/src/ImageSubset.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageSubset.c (revision 33651) +++ trunk/Ohana/src/relphot/src/ImageSubset.c (revision 33651) @@ -0,0 +1,183 @@ +# include "relphot.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nimage = 0; + ImageSubset *image = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // a bit annoying : we read the entire block of data, then extract the columns, then set the image structure values. + // this means I need 3 copies in memory at some point. ugh. + + char type[16]; + + GET_COLUMN (Mcal, "MCAL", float); + GET_COLUMN (dMcal, "MCAL_ERR", float); + GET_COLUMN (imageID, "IMAGE_ID", int); + GET_COLUMN (map, "PHOTOM_MAP", int); + GET_COLUMN (flags, "FLAGS", int); + GET_COLUMN (tzero, "TZERO", int); + GET_COLUMN (trate, "TRATE", short); + + // XXX free the fits table data here + + ALLOCATE (image, ImageSubset, Nrow); + for (i = 0; i < Nrow; i++) { + image[i].Mcal = Mcal[i]; + image[i].dMcal = dMcal[i]; + image[i].imageID = imageID[i]; + image[i].photom_map_id = map[i]; + image[i].flags = flags[i]; + image[i].tzero = tzero[i]; + image[i].trate = trate[i]; + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (Mcal); + free (dMcal); + free (imageID); + free (map); + free (flags); + free (tzero); + free (trate); + + *nimage = Nrow; + return image; +} + +int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "IMAGE_SUBSET"); + + gfits_define_bintable_column (&theader, "E", "MCAL", "zero point offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "zero point error", "magnitudes", 1.0, 0.0); + + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "PHOTOM_MAP", "map", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 1.0*0x8000); + + gfits_define_bintable_column (&theader, "J", "TZERO", "exposure start", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "I", "TRATE", "tti rate", NULL, 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + float *Mcal, *dMcal; + unsigned int *imageID, *map, *flags, *tzero; + unsigned short *trate; + + // create intermediate storage arrays + ALLOCATE (Mcal, float, Nimage); + ALLOCATE (dMcal, float, Nimage); + ALLOCATE (imageID, unsigned int, Nimage); + ALLOCATE (map, unsigned int, Nimage); + ALLOCATE (flags, unsigned int, Nimage); + ALLOCATE (tzero, unsigned int, Nimage); + ALLOCATE (trate, unsigned short, Nimage); + + // assign the storage arrays + for (i = 0; i < Nimage; i++) { + Mcal[i] = image[i].Mcal; + dMcal[i] = image[i].dMcal; + imageID[i] = image[i].imageID; + map[i] = image[i].photom_map_id; + flags[i] = image[i].flags; + tzero[i] = image[i].tzero; + trate[i] = image[i].trate; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MCAL", Mcal, Nimage); + gfits_set_bintable_column (&theader, &ftable, "MCAL_ERR", dMcal, Nimage); + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, Nimage); + gfits_set_bintable_column (&theader, &ftable, "PHOTOM_MAP", map, Nimage); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nimage); + gfits_set_bintable_column (&theader, &ftable, "TZERO", tzero, Nimage); + gfits_set_bintable_column (&theader, &ftable, "TRATE", trate, Nimage); + + free (Mcal); + free (dMcal); + free (imageID); + free (map); + free (flags); + free (tzero); + free (trate); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename); + return FALSE; + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + fflush (f); + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 33650) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 33651) @@ -1,3 +1,6 @@ # include "relphot.h" +void plot_setMcal (double *list, int Npts, StatType *stats, float clouds); + +// see discussion in ImagesOps.c re: IDX_T // array of mosaic definition structures @@ -5,22 +8,37 @@ static Mosaic *mosaic; -// list of all images associated with a mosaic -static off_t *Nimlist; -static off_t **imlist; /* mosaic -> image[] */ -static off_t **bin; /* catalog, measure -> mosaic */ - -// list of mosaics associated with an image -static off_t *mosimage; - -// list of mosaic associated with each image -static int **clist; /* mosaic -> catalog[] */ -static off_t **mlist; /* mosiac -> measure[] */ -static off_t *Nlist; -static off_t *NLIST; +// relationships between the mosaics and their associated images +static off_t *MosaicN_Image; // number of images associated with the given mosaic +static off_t **MosaicToImage; // list of imagesa associated with the given mosaic + +// mosaic index for given image : ImageToMosaic[ImageIndex] = MosaicIndex (ImageIndex : 0 < Nimage) +static off_t *ImageToMosaic; + +// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas]) +// each mosaic has N_onMosaic[MosaicIndex] measurements +static off_t *N_onMosaic; // actual number of measurements on mosaic +static off_t *N_ONMOSAIC; // allocated number of measurements on mosaic + +// relationships between the measure,catalog set and the mosaics: +static off_t **MeasureToMosaic; // Mosaic index from measure,catalog : MeasureToMosaic[cat][meas] = MosaicIndex +static off_t **MosaicToCatalog; // catalog for given measure on mosaic : MosaicCatalog[MosaicIndex][i] = cat (i : 0 < NonMosaic[MosaicIndex]) +static off_t **MosaicToMeasure; // measure for given measure on mosaic : MosaicMeasure[MosaicIndex][i] = cat (i : 0 < NonMosaic[MosaicIndex]) + +// MeasureToMosaic was 'bin' +// MosaicToCatalog was 'clist' +// MosaicToMeasure was 'mlist' + +// N_onMosaic was 'Nlist' +// N_ONMOSAIC was 'NLIST' + +// ImageToMosaic was 'mosimage' + +// MosaicN_image was 'Nimlist' +// MosaicToImage was 'imlist' /* find mosaic frames (unique time periods & photcode name matches mosaic) */ -void initMosaics (Image *image, off_t Nimage) { - - off_t i, j, status, found, NMOSAIC, *NIMLIST; +void initMosaics_old (Image *image, off_t Nimage) { + + off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; unsigned int start, stop; char *pname; @@ -32,9 +50,9 @@ ALLOCATE (mosaic, Mosaic, NMOSAIC); - ALLOCATE (imlist, off_t *, NMOSAIC); - ALLOCATE (Nimlist, off_t, NMOSAIC); - ALLOCATE (NIMLIST, off_t, NMOSAIC); - - ALLOCATE (mosimage, off_t, Nimage); // mosaic to which image belongs + ALLOCATE (MosaicToImage, off_t *, NMOSAIC); + ALLOCATE (MosaicN_Image, off_t, NMOSAIC); + ALLOCATE (MosaicN_IMAGE, off_t, NMOSAIC); + + ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no @@ -45,5 +63,5 @@ /* generate list of unique mosaics */ for (i = 0; i < Nimage; i++) { - mosimage[i] = -1; + ImageToMosaic[i] = -1; /* select valid mosaic images by photcode */ @@ -64,12 +82,12 @@ // add reference from image to mosaic - mosimage[i] = j; + ImageToMosaic[i] = j; /* add image to mosaic image list */ - imlist[j][Nimlist[j]] = i; - Nimlist[j] ++; - if (Nimlist[j] == NIMLIST[j]) { - NIMLIST[j] += 10; - REALLOCATE (imlist[j], off_t, NIMLIST[j]); + MosaicToImage[j][MosaicN_Image[j]] = i; + MosaicN_Image[j] ++; + if (MosaicN_Image[j] == MosaicN_IMAGE[j]) { + MosaicN_IMAGE[j] += 10; + REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]); } @@ -88,12 +106,14 @@ mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + // XXX do we need to do something about flag consistency across a mosaic? + /* add image to mosaic image list */ - NIMLIST[Nmosaic] = 10; - Nimlist[Nmosaic] = 1; - ALLOCATE (imlist[Nmosaic], off_t, NIMLIST[Nmosaic]); - imlist[Nmosaic][0] = i; + MosaicN_IMAGE[Nmosaic] = 10; + MosaicN_Image[Nmosaic] = 1; + ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]); + MosaicToImage[Nmosaic][0] = i; // add reference from image to mosaic - mosimage[i] = Nmosaic; + ImageToMosaic[i] = Nmosaic; Nmosaic ++; @@ -101,12 +121,190 @@ NMOSAIC += 10; REALLOCATE (mosaic, Mosaic, NMOSAIC); - REALLOCATE (imlist, off_t *, NMOSAIC); - REALLOCATE (Nimlist, off_t, NMOSAIC); - REALLOCATE (NIMLIST, off_t, NMOSAIC); - } - } + REALLOCATE (MosaicToImage, off_t *, NMOSAIC); + REALLOCATE (MosaicN_Image, off_t, NMOSAIC); + REALLOCATE (MosaicN_IMAGE, off_t, NMOSAIC); + } + } + + // free this or not? + free (MosaicN_IMAGE); initMosaicGrid (image, Nimage); + + fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic); return; +} + +off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start); + +void sort_times (unsigned int *T, int N) { + +# define SWAPFUNC(A,B){ unsigned int tmp; \ + tmp = T[A]; T[A] = T[B]; T[B] = tmp; \ +} +# define COMPARE(A,B)(T[A] < T[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +/* find mosaic frames (unique time periods & photcode name matches mosaic) */ +void initMosaics (Image *image, off_t Nimage) { + + off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; + unsigned int start, stop, *startTimes, *startTimesMosaic; + char *pname; + + if (!MOSAIC_ZEROPT) return; + + /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no + * entry in the image table that represents this mosaic. Instead, it is an + * internal construct that defines a group of related images + */ + + // generate a list of unique start times (these define the mosaics) + ALLOCATE (startTimes, unsigned int, Nimage); + for (i = 0; i < Nimage; i++) { + startTimes[i] = image[i].tzero; + } + sort_times (startTimes, Nimage); + + Nmosaic = 0; + NMOSAIC = 1000; + ALLOCATE (startTimesMosaic, unsigned int, NMOSAIC); + startTimesMosaic[0] = startTimes[0]; + + for (i = 0; i < Nimage; i++) { + if (startTimes[i] < startTimesMosaic[Nmosaic]) { + fprintf (stderr, "error?\n"); + abort(); + } + if (startTimes[i] == startTimesMosaic[Nmosaic]) continue; + Nmosaic ++; + if (Nmosaic >= NMOSAIC) { + NMOSAIC += 1000; + REALLOCATE (startTimesMosaic, unsigned int, NMOSAIC); + } + startTimesMosaic[Nmosaic] = startTimes[i]; + } + Nmosaic ++; + + // now I have a list of uniq start times, and they are in order + // create the mosaic arrays for these times + ALLOCATE (mosaic, Mosaic, Nmosaic); + + ALLOCATE (MosaicToImage, off_t *, Nmosaic); + ALLOCATE (MosaicN_Image, off_t, Nmosaic); + ALLOCATE (MosaicN_IMAGE, off_t, Nmosaic); + + for (i = 0; i < Nmosaic; i++) { + /* a new mosaic, define ranges */ + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].Mcal = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].dMsys = 0.0; + mosaic[i].Xm = 0.0; + mosaic[i].flags = 0; + mosaic[i].secz = NAN; + mosaic[i].photcode = 0; + + MosaicN_IMAGE[i] = 10; + MosaicN_Image[i] = 0; + ALLOCATE (MosaicToImage[i], off_t, MosaicN_IMAGE[i]); + MosaicToImage[i][0] = -1; + } + + ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs + + // assign each image to a mosaic + for (i = 0; i < Nimage; i++) { + ImageToMosaic[i] = -1; + + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (image[i].photcode); + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + + /* set image time range */ + // start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); + // stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + + start = image[i].tzero; + stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + + j = findMosaic(startTimesMosaic, Nmosaic, start); + if (j == -1) { + fprintf (stderr, "error?\n"); + abort(); + } + + // add reference from image to mosaic + ImageToMosaic[i] = j; + + // have we already found this mosaic? + found = (MosaicN_Image[j] > 0); + + /* add image to mosaic image list */ + MosaicToImage[j][MosaicN_Image[j]] = i; + MosaicN_Image[j] ++; + if (MosaicN_Image[j] == MosaicN_IMAGE[j]) { + MosaicN_IMAGE[j] += 10; + REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]); + } + if (found) continue; + + /* a new mosaic, define ranges */ + if (mosaic[j].start != start) { + fprintf (stderr, "error?\n"); + abort(); + } + mosaic[j].stop = stop; + mosaic[j].Mcal = 0.0; + mosaic[j].dMcal = 0.0; + mosaic[j].dMsys = 0.0; + mosaic[j].Xm = 0.0; + mosaic[j].flags = image[i].flags; + mosaic[j].secz = image[i].secz; + mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + } + + // free this or not? + free (MosaicN_IMAGE); + free (startTimes); + free (startTimesMosaic); + + initMosaicGrid (image, Nimage); + + fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic); + return; +} + +// use bisection to find the overlapping mosaic +off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start) { + + off_t Nlo, Nhi, N; + + // find the last mosaic before start + Nlo = 0; // startTimes[Nlo] guaranteed to be <= start + Nhi = Nmosaic - 1; // startTimes[Nhi] guaranteed to be >= start + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (startTimes[N] < start) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N, Nmosaic - 1); + } + } + + // check for the matched mosaic starting from Nlo + for (N = Nlo; N < Nmosaic; N++) { + if (start > startTimes[N]) continue; + return (N); + } + return (-1); } @@ -129,6 +327,6 @@ dS = 0.0; Mcal = dMcal = Xm = 0; - for (j = 0; j < Nimlist[i]; j++) { - m = imlist[i][j]; + for (j = 0; j < MosaicN_Image[i]; j++) { + m = MosaicToImage[i][j]; NX = image[m].NX; NY = image[m].NY; @@ -157,4 +355,12 @@ dMcal += image[m].dMcal; Xm += image[m].Xm; + + // XXX: how does this work with UBERCAL? We want to keep the Mcal values supplied by ubercal, but + // solve for a single offset for each exposure (Mosaic.Mcal). + + // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not) + + // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value + /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ image[m].Mcal = 0.0; @@ -162,5 +368,5 @@ image[m].Xm = NAN_S_SHORT; } - dS /= Nimlist[i]; + dS /= MosaicN_Image[i]; strcpy (mosaic[i].coords.ctype, "RA---TAN"); mosaic[i].coords.crval1 = Rmin; @@ -176,7 +382,7 @@ RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords); - mosaic[i].Mcal = Mcal / Nimlist[i]; - mosaic[i].dMcal = dMcal / Nimlist[i]; - mosaic[i].Xm = Xm / Nimlist[i]; + mosaic[i].Mcal = Mcal / MosaicN_Image[i]; + mosaic[i].dMcal = dMcal / MosaicN_Image[i]; + mosaic[i].Xm = Xm / MosaicN_Image[i]; } if (!USE_GRID) return; @@ -197,9 +403,12 @@ image = getimages (&Nimage, NULL); + // XXX I think this is OK in the ubercal context, but probably need to skip UBERCAL + // images? (no need to update them) + // 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]; + for (j = 0; j < MosaicN_Image[i]; j++) { + im = MosaicToImage[i][j]; image[im].Mcal += mosaic[i].Mcal; image[im].dMcal = mosaic[i].dMcal; @@ -221,22 +430,22 @@ if (!MOSAIC_ZEROPT) return; - ALLOCATE (bin, off_t *, Ncatalog); + ALLOCATE (MeasureToMosaic, off_t *, Ncatalog); for (i = 0; i < Ncatalog; i++) { - ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1)); - for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1; + ALLOCATE (MeasureToMosaic[i], off_t, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToMosaic[i][j] = -1; } if (doMosaicList) { /* mosaic -> measure */ - ALLOCATE (Nlist, off_t, Nmosaic); - ALLOCATE (NLIST, off_t, Nmosaic); - ALLOCATE (clist, int *, Nmosaic); - ALLOCATE (mlist, off_t *, Nmosaic); + ALLOCATE (N_onMosaic, off_t, Nmosaic); + ALLOCATE (N_ONMOSAIC, off_t, Nmosaic); + ALLOCATE (MosaicToCatalog, off_t *, Nmosaic); + ALLOCATE (MosaicToMeasure, 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]); + N_onMosaic[i] = 0; + N_ONMOSAIC[i] = 100; + ALLOCATE (MosaicToCatalog[i], off_t, N_ONMOSAIC[i]); + ALLOCATE (MosaicToMeasure[i], off_t, N_ONMOSAIC[i]); } } @@ -251,18 +460,18 @@ for (i = 0; i < Ncatalog; i++) { - free (bin[i]); - } - free (bin); + free (MeasureToMosaic[i]); + } + free (MeasureToMosaic); if (doMosaicList) { /* mosaic -> measure */ for (i = 0; i < Nmosaic; i++) { - free (clist[i]); - free (mlist[i]); - } - free (Nlist); - free (NLIST); - free (clist); - free (mlist); + free (MosaicToCatalog[i]); + free (MosaicToMeasure[i]); + } + free (N_onMosaic); + free (N_ONMOSAIC); + free (MosaicToCatalog); + free (MosaicToMeasure); } } @@ -310,5 +519,5 @@ } - mosID = mosimage[idx]; + mosID = ImageToMosaic[idx]; if (mosID < 0) { Image *image = getimage(idx); @@ -320,5 +529,7 @@ { Image *image = getimage(idx); - unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1); + // XXX we are now matching with just tzero. be careful for cameras with drift + // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1); + unsigned int imageStart = image[0].tzero; if (imageStart != mosaic[mosID].start) { fprintf (stderr, "error in image to mosaic match\n"); @@ -327,15 +538,15 @@ } - bin[cat][meas] = mosID; + MeasureToMosaic[cat][meas] = mosID; if (doMosaicList) { - clist[mosID][Nlist[mosID]] = cat; - mlist[mosID][Nlist[mosID]] = meas; - Nlist[mosID] ++; + MosaicToCatalog[mosID][N_onMosaic[mosID]] = cat; + MosaicToMeasure[mosID][N_onMosaic[mosID]] = meas; + N_onMosaic[mosID] ++; - if (Nlist[mosID] == NLIST[mosID]) { - NLIST[mosID] += 100; - REALLOCATE (clist[mosID], int, NLIST[mosID]); - REALLOCATE (mlist[mosID], off_t, NLIST[mosID]); + if (N_onMosaic[mosID] == N_ONMOSAIC[mosID]) { + N_ONMOSAIC[mosID] += 100; + REALLOCATE (MosaicToCatalog[mosID], off_t, N_ONMOSAIC[mosID]); + REALLOCATE (MosaicToMeasure[mosID], off_t, N_ONMOSAIC[mosID]); } } @@ -350,5 +561,5 @@ if (!MOSAIC_ZEROPT) return (0); - i = bin[cat][meas]; + i = MeasureToMosaic[cat][meas]; if (i == -1) return (NAN); @@ -358,9 +569,9 @@ } -int setMmos (Catalog *catalog, int PoorImages) { - - off_t i, j, m, c, N, Nmax; - int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nsecfilt; - float Msys, Mrel, Mcal, Mgrid; +int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + off_t i, j, m, c, n, N, Nmax; + int mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nbright; + float Msys, Mrel, Mcal, Mgrid, Mflat; double *list, *dlist, *Mlist, *dMlist; StatType stats; @@ -372,7 +583,11 @@ image = getimages (&N, NULL); - Nsecfilt = GetPhotcodeNsecfilt (); + fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); + + int Nsecfilt = GetPhotcodeNsecfilt (); if (PoorImages) { + // XXX use bad stars and measurements for PoorImages? or not? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; IMAGE_BAD = 0; } @@ -380,5 +595,5 @@ Nmax = 0; for (i = 0; i < Nmosaic; i++) { - Nmax = MAX (Nmax, Nlist[i]); + Nmax = MAX (Nmax, N_onMosaic[i]); } ALLOCATE (list, double, Nmax); @@ -397,9 +612,27 @@ } + // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0: + // we optionally do not recalibrate images with UBERCAL zero points + if (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL) { + mosaic[i].ubercalDist = 0; + // propagate ubercalDist to the images + for (j = 0; j < MosaicN_Image[i]; j++) { + off_t im = MosaicToImage[i][j]; + image[im].ubercalDist = mosaic[i].ubercalDist; + // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist); + } + if (KEEP_UBERCAL) continue; + } + + int minUbercalDist = 1000; + + // number of stars to measure the bright-end scatter + Nbright = 0; + N = 0; - for (j = 0; j < Nlist[i]; j++) { + for (j = 0; j < N_onMosaic[i]; j++) { - m = mlist[i][j]; - c = clist[i][j]; + m = MosaicToMeasure[i][j]; + c = MosaicToCatalog[i][j]; if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { @@ -407,5 +640,5 @@ continue; } - Mcal = getMcal (m, c); + Mcal = getMcal (m, c, flatcorr, catalog); if (isnan(Mcal)) { Ncal++; @@ -423,4 +656,11 @@ } + // image.Mcal is not supposed to include the flat-field correction, so we need to + // apply that offset as well here for this image (in other words, each detection is + // being compared to the model, excluding the zero point, Mcal. The model includes + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) + + Mflat = getMflat (m, c, flatcorr, catalog); + n = catalog[c].measureT[m].averef; Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); @@ -429,18 +669,29 @@ continue; } - list[N] = Msys - Mrel - Mcal - Mgrid; + + PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); + if (!code) goto skip; + if (code->equiv < 1) goto skip; + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec == -1) goto skip; + minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist); + + skip: + list[N] = Msys - Mrel - Mcal - Mgrid + Mflat; dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); - Mlist[N] = Msys; - dMlist[N] = list[N]; + if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { + Mlist[Nbright] = list[N]; + dMlist[Nbright] = dlist[N]; + Nbright ++; + } N++; } - /* Nlist[i] is all measurements, N is good measurements */ + /* N_onMosaic[i] is all measurements, N is good measurements */ /* too few good measurements or too many bad measurements (skip in PoorImages run) */ - if (!PoorImages) { - mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]); + mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[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 (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]); } mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW; Nfew ++; @@ -450,15 +701,30 @@ } 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; + if (VERBOSE2 && PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas, N); + mosaic[i].Mcal = stats.mean; mosaic[i].dMcal = stats.error; - mosaic[i].dMsys = stats.sigma; mosaic[i].nFitPhotom = N; mosaic[i].Xm = 100.0*log10(stats.chisq); + plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); + + // bright end scatter + liststats (Mlist, dMlist, Nbright, &stats); + mosaic[i].dMsys = stats.sigma; + if (mosaic[i].Mcal < -CLOUD_TOLERANCE) { mosaic[i].Mcal = 0.0; + } + + // minUbercalDist calculated here is the min value for any star owned by this image + // since this particular image is tied to that star, bump its distance by 1 + mosaic[i].ubercalDist = minUbercalDist + 1; + + // propagate ubercalDist to the images + for (j = 0; j < MosaicN_Image[i]; j++) { + off_t im = MosaicToImage[i][j]; + image[im].ubercalDist = mosaic[i].ubercalDist; + // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist); } } @@ -481,4 +747,5 @@ // 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. +// XXX this function is no longer used because we force significantly negative clouds to 0.0 int rationalize_mosaics (Catalog *catalog, int Ncatalog) { @@ -515,5 +782,5 @@ if (VERBOSE2 && (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE)) { - fprintf (stderr, "cloud-free: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal); } if (VERBOSE2 && (mosaic[i].Mcal < -CLOUD_TOLERANCE)) { @@ -521,8 +788,8 @@ // 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); + fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal); } if (VERBOSE2 && (mosaic[i].Mcal > CLOUD_TOLERANCE)) { - fprintf (stderr, "cloudy : %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + fprintf (stderr, "cloudy : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal); } } @@ -579,5 +846,5 @@ // find the source of this measurement (skip unassigned measurements) - nMos = bin[i][m]; + nMos = MeasureToMosaic[i][m]; if (nMos == -1) continue; @@ -641,9 +908,9 @@ } imageOffset[i] = dM / NSlist[i]; - // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[imlist[i][0]].name, dM, NSlist[i], imageOffset[i]); + // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[MosaicToImage[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]); + // fprintf (stderr, "correction: %s : %f\n", image[MosaicToImage[i][0]].name, imageOffset[i]); // } @@ -727,5 +994,5 @@ } -StatType statsMosaicN (Catalog *catalog) { +StatType statsMosaicN (Catalog *catalog, FlatCorrectionTable *flatcorr) { off_t i, j, m, c, n, N; @@ -746,10 +1013,10 @@ N = 0; - for (j = 0; j < Nlist[i]; j++) { - - m = mlist[i][j]; - c = clist[i][j]; - - Mcal = getMcal (m, c); + for (j = 0; j < N_onMosaic[i]; j++) { + + m = MosaicToMeasure[i][j]; + c = MosaicToCatalog[i][j]; + + Mcal = getMcal (m, c, flatcorr, catalog); if (isnan(Mcal)) continue; Mgrid = getMgrid (m, c); @@ -833,4 +1100,7 @@ Nmark = Nscatter = Noffset = 0; for (i = 0; i < Nmosaic; i++) { + // if we are keeping ubercal sacrosanct, then we should not be allowed to break them... + if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + mark = FALSE; if (mosaic[i].dMcal > MaxScatter) { @@ -872,5 +1142,5 @@ N = 0; for (i = 0; i < Nmosaic; i++) - N = MAX (N, Nlist[i]); + N = MAX (N, N_onMosaic[i]); ALLOCATE (xlist, double, N); @@ -881,8 +1151,8 @@ Xmin = Ymin = +360.0; Xmax = Ymax = -360.0; - for (j = 0; j < Nlist[i]; j++) { + for (j = 0; j < N_onMosaic[i]; j++) { - m = mlist[i][j]; - c = clist[i][j]; + m = MosaicToMeasure[i][j]; + c = MosaicToCatalog[i][j]; if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; @@ -959,9 +1229,9 @@ Imax = 0; - Nmax = Nimlist[0]; - for (i = 0; i < Nmosaic; i++) { - if (Nimlist[i] > Nmax) { + Nmax = MosaicN_Image[0]; + for (i = 0; i < Nmosaic; i++) { + if (MosaicN_Image[i] > Nmax) { Imax = i; - Nmax = Nimlist[i]; + Nmax = MosaicN_Image[i]; } } @@ -969,4 +1239,4 @@ *refmosaic = &mosaic[Imax]; *Nimage = Nmax; - return (imlist[Imax]); -} + return (MosaicToImage[Imax]); +} Index: trunk/Ohana/src/relphot/src/SetSignals.c =================================================================== --- trunk/Ohana/src/relphot/src/SetSignals.c (revision 33650) +++ trunk/Ohana/src/relphot/src/SetSignals.c (revision 33651) @@ -30,7 +30,8 @@ for (i = 0; i < 36; i++) { switch (i) { - /* can't redirect this signals */ + /* can't redirect these signals */ case SIGKILL: /* kill -9: cannot be caught or ignored */ case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */ + /* ignore these signals */ case SIGCHLD: /* child halted: ignore */ Index: trunk/Ohana/src/relphot/src/Shutdown.c =================================================================== --- trunk/Ohana/src/relphot/src/Shutdown.c (revision 33650) +++ trunk/Ohana/src/relphot/src/Shutdown.c (revision 33651) @@ -1,5 +1,5 @@ # include "relphot.h" -static FITS_DB *db; +static FITS_DB *db = NULL; void set_db (FITS_DB *in) { @@ -22,5 +22,5 @@ SetProtect (TRUE); - gfits_db_close (db); + if (db) gfits_db_close (db); fprintf (stderr, "ERROR: addstar halted\n"); exit (1); Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 33650) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 33651) @@ -47,5 +47,5 @@ } -int setMrel (Catalog *catalog, int Ncatalog) { +int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { off_t j, k, m; @@ -59,4 +59,7 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { + // XXX accumulate all secfilt values in a single pass? + + int minUbercalDist = 1000; int Ns; @@ -66,56 +69,59 @@ int Nsec = GetPhotcodeNsec(thisCode); - /* calculate the average value for a single star */ + /* calculate the average mag in this SEC photcode for a single star */ // skip bad stars if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + + N = 0; m = catalog[i].averageT[j].measureOffset; - - N = 0; for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode - int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); - if (ecode != thisCode) { continue; } - - if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { - Nbad ++; - continue; - } - // XXX allow REF stars (no Image Entry) to be included in the calculation this - // should be optionally set, and should allow for REF stars to be downweighted by - // more than their reported errors. how such info is carried is unclear... + PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode); + if (!code) continue; + if (code->equiv != thisCode) { continue; } + + if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; } + if (getImageEntry (m, i) < 0) { Mcal = Mmos = Mgrid = 0; } else { - Mcal = getMcal (m, i); - if (isnan(Mcal)) { - Ncal ++; - continue; - } + Mcal = getMcal (m, i, flatcorr, catalog); + if (isnan(Mcal)) { Ncal ++; continue; } Mmos = getMmos (m, i); - if (isnan(Mmos)) { - Nmos ++; - continue; - } + if (isnan(Mmos)) { Nmos ++; continue; } Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) { - Ngrid++; - continue; - } + if (isnan(Mgrid)) { Ngrid++; continue; } } Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); - if (isnan(Msys)) { - Nsys++; - continue; - } + if (isnan(Msys)) { Nsys++; continue; } list[N] = Msys - Mcal - Mmos - Mgrid; - dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); + + int myUbercalDist = getUbercalDist(m,i); + minUbercalDist = MIN(minUbercalDist, myUbercalDist); + + // dlist gives the error, which is used as the weight in WT_MEAN. + // we can modify the resulting weight in a few ways: + // 1) MIN_ERROR guarantees a floor + // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode + // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages + // 4) some reference photcode of some kind can be specified as fixed and have a high weight + dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR); + + // up-weight the ubercal values (or convergence can take a long time...) + if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR); + } // tie down reference photometry if the -refcode (code) option is selected + // eg, -refcode g_SDSS + // this probably makes no sense in the context of multifilter analysis if (refPhotcode) { - if (GetPhotcodeEquivCodebyCode(catalog[i].measureT[m].photcode) == refPhotcode[0].equiv) { - // increase the weight by a factor of 100: - dlist[N] = 0.01*catalog[i].measureT[m].dM; + if (code->code == refPhotcode->code) { + // tiny error -> large weight + // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR); + dlist[N] = 0.0001; } } @@ -125,4 +131,5 @@ // when performing the grid analysis, STAR_TOOFEW will be set to 1; if (N <= STAR_TOOFEW) { /* too few measurements */ + // 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); catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW; Nfew ++; @@ -132,8 +139,10 @@ liststats (list, dlist, N, &stats); - + catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma; - catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; + catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; + + catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; } } @@ -144,14 +153,47 @@ } -int setMrelOutput (Catalog *catalog, int Ncatalog, int mark) { +int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) { + + off_t k; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + off_t m = average[0].measureOffset; + + for (k = 0; k < average[0].Nmeasure; k++, m++) { + fprintf (stderr, "meas: %08x\n", measure[m].dbFlags); + } + + int Ns; + for (Ns = 0; Ns < Nsecfilt; Ns++) { + fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags); + } + return 1; +} + +# define MARK_SKIP_MEAS \ + catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; + +// setMrel and setMrelOutput are extremely similar, but have slightly different implications: +// * setMrel uses the internal Tiny structures only +// * setMrelOutput skips stars for which there are too few good measurements +// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass + +// setMave is also similar to the above. but: +// * setMave is called by relphot -update-objects +// * setMave excludes all detections with (PSF_QF < 0.85), setMrelOutput allows these for PASS > 2 +// * setMave updates 2MASS average flags +// * setMave updates average EXT flags (PS1 and 2MASS) +int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) { int i, N; off_t j, k, m, Nmax; - float Msys, Mcal, Mmos, Mgrid; - double *list, *dlist; - StatType stats; - int Nsec, Nsecfilt, ecode; - - Nsecfilt = GetPhotcodeNsecfilt (); + float Msys, Mcal, Mmos, Mgrid, Map; + double *list, *dlist, *aplist, *daplist; + StatType stats, apstats; + int Nsec; + + int Nsecfilt = GetPhotcodeNsecfilt (); /* Nmeasure is now different, need to reallocate */ @@ -165,51 +207,185 @@ ALLOCATE (dlist, double, MAX (1, Nmax)); - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - /* skip stars already calibrated */ - if (catalog[i].found[j]) continue; + ALLOCATE (aplist, double, MAX (1, Nmax)); + ALLOCATE (daplist, double, MAX (1, Nmax)); + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) { + fprintf (stderr, "test obj\n"); + print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure); + } + + int minUbercalDist = 1000; int Ns; for (Ns = 0; Ns < Nphotcodes; Ns++) { + int thisCode = photcodes[Ns][0].code; Nsec = GetPhotcodeNsec(thisCode); + + /* star/photcodes already calibrated */ + if (catalog[i].found[Nsecfilt*j+Nsec]) continue; + + int Ncode = 0; + int Next = 0; + int haveSynth = FALSE; + int haveUbercal = FALSE; + int isRefPhot = FALSE; N = 0; m = catalog[i].averageT[j].measureOffset; for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode - ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); - if (ecode != thisCode) { continue; } + PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode); + if (!code) continue; + if (code->equiv != thisCode) { continue; } + Ncode ++; if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - // XXX allow REF stars (no Image Entry) to be included in the calculation this - // should be optionally set, and should allow for REF stars to be downweighted by - // more than their reported errors. how such info is carried is unclear... if (getImageEntry (m, i) < 0) { + // these detetions have no image (eg, ref values such as 2MASS) Mcal = Mmos = Mgrid = 0; + isRefPhot = TRUE; } else { - Mcal = getMcal (m, i); - if (isnan(Mcal)) continue; + Mcal = getMcal (m, i, flatcorr, catalog); + if (isnan(Mcal)) { + MARK_SKIP_MEAS; + continue; + } Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; + if (isnan(Mmos)) { + MARK_SKIP_MEAS; + continue; + } Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; + if (isnan(Mgrid)) { + MARK_SKIP_MEAS; + continue; + } } Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); + if (isnan(Msys)) { + if (!isRefPhot) { + MARK_SKIP_MEAS; + } + continue; + } list[N] = Msys - Mcal - Mmos - Mgrid; - dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); + + Map = PhotAper (&catalog[i].measure[m]); + aplist[N] = Map - Mcal - Mmos - Mgrid; + + int myUbercalDist = getUbercalDist(m,i); + minUbercalDist = MIN(minUbercalDist, myUbercalDist); + + // count the extended detections + if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) { + if (!isnan(catalog[i].measure[m].Map)) { + if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { + Next ++; + } + } + } + // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013) + if ((catalog[i].measure[m].photcode >= 2011) && (catalog[i].measure[m].photcode <= 2013)) { + if (catalog[i].measure[m].photFlags & 0x00c00000) { + Next ++; + } + if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) { + // detections without one of these bits should only be used in PASS_1 + MARK_SKIP_MEAS; + continue; + } + } + // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars) + if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) { + if (pass < 4) { + MARK_SKIP_MEAS; + continue; + } + haveSynth = TRUE; + } + + // dlist gives the error, which is used as the weight in WT_MEAN. + // we can modify the resulting weight in a few ways: + // 1) MIN_ERROR guarantees a floor + // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode + // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages + // 4) some reference photcode of some kind can be specified as fixed and have a high weight + dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR); + + // up-weight the ubercal values (or convergence can take a long time...) (XXX make this optional?) + if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR); + haveUbercal = TRUE; + } + + // tie down reference photometry if the -refcode (code) option is selected + // eg, -refcode g_SDSS + // this probably makes no sense in the context of multifilter analysis + if (refPhotcode) { + if (code->code == refPhotcode->code) { + // tiny error -> large weight + // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR); + dlist[N] = 0.0001; + } + } N++; } if (N < 1) continue; + for (m = 0; m < N; m++) { + daplist[m] = dlist[m]; + } + + // XXX force WT_MEAN or MEAN here? liststats (list, dlist, N, &stats); - if (mark) catalog[i].found[j] = TRUE; + liststats (aplist, daplist, N, &apstats); + catalog[i].found[Nsecfilt*j+Nsec] = TRUE; + switch (pass) { + case 0: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0; + break; + case 1: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1; + break; + case 2: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2; + break; + case 3: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3; + break; + case 4: + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4; + break; + } + + if (haveSynth) { + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH; + } + if (haveUbercal) { + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL; + } /* use sigma or error in dM for output? */ catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; + catalog[i].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean; catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma); catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; + catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode; + catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas; + + catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80; + catalog[i].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20; + catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist; + + if ((Next > 0) && (Next > 0.5*N)) { + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT; + } } } @@ -218,10 +394,24 @@ free (list); free (dlist); + + free (aplist); + free (daplist); return (TRUE); } + +/* grab Nsec for named photcode */ +# define NAMED_PHOTCODE_NSEC(MY_NSEC,NAME) \ + short MY_NSEC = -1; \ + { \ + PhotCode *code = GetPhotcodebyName (NAME); \ + if (code) { \ + MY_NSEC = GetPhotcodeNsec (code->equiv); \ + } } // For each average object, set the average mags based on existing equiv photometry. // NOTE: this function operates on the real Measure & Average structures, not the // MeasureTiny & AverageTiny structures +// NOTE: this function is called on the remote machine -- make sure any (global) options +// are passed to the relphot_client program int setMave (Catalog *catalog, int Ncatalog) { @@ -231,5 +421,4 @@ double *list, *dlist; StatType stats; - PhotCode *code; DVOAverageFlags flagBits; @@ -248,12 +437,22 @@ Nsecfilt = GetPhotcodeNsecfilt (); -# define PSFQUALSTATS 1 + // we want to raise some bits on the 2MASS (JHK) secfilt flags, if we have 2MASS data + NAMED_PHOTCODE_NSEC (Nsec_J, "2MASS_J"); + NAMED_PHOTCODE_NSEC (Nsec_H, "2MASS_H"); + NAMED_PHOTCODE_NSEC (Nsec_K, "2MASS_K"); + for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { // update average photometry for each of the average filters + + // XXX Note that this would be faster if we had an array of results and accumulated + // them in a single pass + + int Next = 0; + for (Ns = 0; Ns < Nsecfilt; Ns++) { - code = GetPhotcodebyNsec (Ns); + PhotCode *code = GetPhotcodebyNsec (Ns); Nc = code[0].code; @@ -267,7 +466,13 @@ if (isnan(Msys)) continue; + // reject POOR detections (PSF_QF < 0.85) and count extended detections // XXX only apply this filter for psphot data from GPC1 for now... - if (PSFQUALSTATS && (catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) { - if (catalog[i].measure[m].psfQual < 0.85) continue; + if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) { + if (catalog[i].measure[m].psfQual < 0.85) continue; + if (!isnan(catalog[i].measure[m].Map)) { + if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { + Next ++; + } + } } @@ -287,60 +492,92 @@ catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N; catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas; + + catalog[i].secfilt[Nsecfilt*j+Ns].M_80 = 1000 * stats.Upper80; + catalog[i].secfilt[Nsecfilt*j+Ns].M_20 = 1000 * stats.Lower20; + + if ((Next > 0) && (Next > 0.5*N)) { + catalog[i].secfilt[Nsecfilt*j+Ns].flags |= ID_SECF_OBJ_EXT; + } } // update average flags based on the detection stats. - // XXX we need to clean this up -- this is a serious set of hacks... - if (PSFQUALSTATS) { - int Galaxy2MASS, GalaxySDSS, goodPS1, nEXT, nPSF, good2MASS; - - Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000]) - GalaxySDSS = FALSE; // best guess for galaxy based on SDSS measurements (XXX need to fix SDSS flags) - goodPS1 = FALSE; // true if any PS1 measurements have psfQual > 0.85 - good2MASS = FALSE; // true if 2MASS J measurements have significant detections - nEXT = nPSF = 0; // number of PS1 PSF vs EXT measurements - - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - - // PS1 data : - if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) { - if (catalog[i].measure[m].psfQual > 0.85) { - goodPS1 = TRUE; - if (!isnan(catalog[i].measure[m].Map)) { - if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { - nEXT ++; - } else { - nPSF ++; - } + + int Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000]) + int goodPS1 = FALSE; // true if any PS1 measurements have psfQual > 0.85 + int good2MASS = FALSE; // true if 2MASS J measurements have significant detections + int nEXT = 0; + int nPSF = 0; // number of PS1 PSF vs EXT measurements + int have2MASS = FALSE; + + // count, flag good and extended detections + m = catalog[i].average[j].measureOffset; + for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + + // PS1 data : + if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) { + if (catalog[i].measure[m].psfQual > 0.85) { + goodPS1 = TRUE; + if (!isnan(catalog[i].measure[m].Map)) { + if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) { + nEXT ++; + } else { + nPSF ++; } } } + } - // 2MASS data: - if (catalog[i].measure[m].photcode == 2011) { - if (catalog[i].measure[m].photFlags & 0x00c00000) { - Galaxy2MASS = TRUE; - } - if (catalog[i].measure[m].photFlags & 0x00000007) { - good2MASS = TRUE; - } - } - } - - // we attempt to set a few flags here; reset those bits before trying: - catalog[i].average[j].flags &= ~flagBits; - - if (nEXT && (nEXT > nPSF)) { - catalog[i].average[j].flags |= ID_OBJ_EXT; - } - if (goodPS1) { - catalog[i].average[j].flags |= ID_OBJ_GOOD; - } - if (Galaxy2MASS) { - catalog[i].average[j].flags |= ID_OBJ_EXT_ALT; - } - if (good2MASS) { - catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT; - } + // 2MASS data J-band flags + if (catalog[i].measure[m].photcode == 2011) { + // only need to do this once (always have JHK triplet; galaxy flag is same for all 3) + have2MASS = TRUE; + if (catalog[i].measure[m].photFlags & 0x00c00000) { + Galaxy2MASS = TRUE; + } + if (catalog[i].measure[m].photFlags & 0x00000007) { + good2MASS = TRUE; + if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_0; + } else { + if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_1; + } + } + // 2MASS data H-band flags + if (catalog[i].measure[m].photcode == 2012) { + if (catalog[i].measure[m].photFlags & 0x00000007) { + good2MASS = TRUE; + if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_0; + } else { + if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_1; + } + } + // 2MASS data K-band flags + if (catalog[i].measure[m].photcode == 2013) { + if (catalog[i].measure[m].photFlags & 0x00000007) { + good2MASS = TRUE; + if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_0; + } else { + if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_1; + } + } + } + + // we attempt to set a few flags here; reset those bits before trying: + catalog[i].average[j].flags &= ~flagBits; + + // XXX set the secfilt bits? + if (nEXT && (nEXT > nPSF)) { + catalog[i].average[j].flags |= ID_OBJ_EXT; + } + if (goodPS1) { + catalog[i].average[j].flags |= ID_OBJ_GOOD; + } + if (Galaxy2MASS) { + catalog[i].average[j].flags |= ID_OBJ_EXT_ALT; + if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_SECF_OBJ_EXT; + if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_SECF_OBJ_EXT; + if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_SECF_OBJ_EXT; + } + if (good2MASS) { + catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT; } } @@ -353,5 +590,5 @@ /* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */ -int setMcalOutput (Catalog *catalog, int Ncatalog) { +int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { int i; @@ -368,5 +605,5 @@ for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - Mcal = getMcal (m, i); + Mcal = getMcal (m, i, flatcorr, catalog); if (isnan(Mcal)) continue; Mmos = getMmos (m, i); @@ -377,4 +614,8 @@ // set the output calibration catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid; + + if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + myAssert (isfinite(catalog[i].measure[m].Mcal), "oops, broke an ubercal mag"); + } } } @@ -385,5 +626,5 @@ void clean_stars (Catalog *catalog, int Ncatalog) { - int i, j, Ndel, Nave, Ntot, mark, Ns; + int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan; float dM, Xm; double Chisq, MaxScatter, MaxChisq; @@ -402,4 +643,9 @@ int Nsecfilt = GetPhotcodeNsecfilt (); + + // XX int oldPLOTSTUFF = PLOTSTUFF; + // XX PLOTSTUFF = TRUE; + // XX plot_chisq (catalog, Ncatalog); + // XX PLOTSTUFF = oldPLOTSTUFF; // eliminate bad stars using the stats for a single secfilt at a time @@ -430,5 +676,5 @@ fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq); - Ndel = Nave = 0; + Ndel = Nave = Nscat = Nnan = Nchi = 0; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { @@ -440,4 +686,7 @@ catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR; Ndel ++; + if (dM > MaxScatter) { Nscat ++; } + if (Xm == NAN_S_SHORT) { Nnan ++; } + if (Chisq > MaxChisq) { Nchi ++; } } else { catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR; @@ -446,5 +695,5 @@ } } - fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave); + fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave); initstats (STATMODE); } @@ -456,5 +705,5 @@ # define NSIGMA_CLIP 3.0 # define NSIGMA_REJECT 5.0 -void clean_measures (Catalog *catalog, int Ncatalog, int final) { +void clean_measures (Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr) { off_t j, k, m, Nmax, Ndel, Nave; @@ -492,10 +741,10 @@ for (Ns = 0; Ns < Nphotcodes; Ns++) { - /* on final processing, skip stars already measured */ - if (final && catalog[i].found[j]) continue; - int thisCode = photcodes[Ns][0].code; int Nsec = GetPhotcodeNsec(thisCode); + /* on final processing, skip stars already measured */ + if (final && catalog[i].found[Nsecfilt*j + Nsec]) continue; + /* skip bad stars to prevent them from becoming good (on inner sample) */ if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; @@ -509,6 +758,6 @@ if (ecode != thisCode) { continue; } - /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */ - Mcal = getMcal (m, i); + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers + Mcal = getMcal (m, i, flatcorr, catalog); if (isnan(Mcal)) { Ncal ++; continue; } Mmos = getMmos (m, i); @@ -522,5 +771,8 @@ N++; } - if (N <= TOOFEW) { Nfew ++; continue; } + if (N <= TOOFEW) { + Nfew ++; + continue; + } /* 3-sigma clip based on stats of inner 50% */ @@ -553,6 +805,6 @@ if (ecode != thisCode) { continue; } - /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */ - Mcal = getMcal (m, i); + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers + Mcal = getMcal (m, i, flatcorr, catalog); if (isnan(Mcal)) continue; Mmos = getMmos (m, i); @@ -572,5 +824,9 @@ /* mark bad measures (> 3 sigma deviant) */ for (k = 0; k < N; k++) { - if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) { + // treat the scatter of the star as a systematic term. this is a bit of an + // over-estimage (a perfect Gauss distribution with perfect errors would have + // mySigma = sqrt(2) too large) + float mySigma = hypot (stats.sigma, dlist[k]); + if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) { catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM; if (final) { @@ -591,7 +847,9 @@ free (ilist); free (tlist); -} - -StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) { + + +} + +StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr) { off_t j, k, m, Ntot; @@ -625,5 +883,5 @@ int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); if (ecode != seccode) { continue;} - Mcal = getMcal (m, i); + Mcal = getMcal (m, i, flatcorr, catalog); if (isnan(Mcal)) { continue;} Mmos = getMmos (m, i); @@ -712,4 +970,5 @@ dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + if (isnan(dM)) continue; list[n] = dM; dlist[n] = 1; Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 33650) +++ trunk/Ohana/src/relphot/src/args.c (revision 33651) @@ -38,5 +38,4 @@ UserPatch.Dmin = -90; UserPatch.Dmax = +90; - UserPatchSelect = FALSE; if ((N = get_argument (argc, argv, "-region"))) { remove_argument (N, &argc, argv); @@ -49,5 +48,12 @@ UserPatch.Dmax = atof (argv[N]); remove_argument (N, &argc, argv); - UserPatchSelect = TRUE; + } + + /* specify region file by name (eg n0000/0000.00) */ + UserCatalog = NULL; + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserCatalog = strcreate (argv[N]); + remove_argument (N, &argc, argv); } @@ -68,4 +74,32 @@ } + // XXX for the moment, make this selection manual. it needs to be automatic + // based on the state of the SkyTable + PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + PARALLEL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the remote jobs and waiting for completion, + // relphot will simply list the remote command and wait for the user to signal completion + PARALLEL_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-manual"))) { + PARALLEL = TRUE; // -parallel-manual implies -parallel + PARALLEL_MANUAL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the relphot_client jobs remotely, they are + // run in serial via 'system' + PARALLEL_SERIAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-serial"))) { + if (PARALLEL_MANUAL) { + fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); + exit (1); + } + PARALLEL = TRUE; // -parallel-serial implies -parallel + PARALLEL_SERIAL = TRUE; + remove_argument (N, &argc, argv); + } + PLOTSTUFF = FALSE; if ((N = get_argument (argc, argv, "-plot"))) { @@ -90,5 +124,5 @@ } - strcpy (STATMODE, "CHI_INNER_WTMEAN"); + strcpy (STATMODE, "WT_MEAN"); if ((N = get_argument (argc, argv, "-statmode"))) { remove_argument (N, &argc, argv); @@ -128,4 +162,11 @@ } + UPDATE_CATFORMAT = NULL; + if ((N = get_argument (argc, argv, "-update-catformat"))) { + remove_argument (N, &argc, argv); + UPDATE_CATFORMAT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + SAVE_IMAGE_UPDATES = TRUE; if ((N = get_argument (argc, argv, "-skip-image-updates"))) { @@ -140,4 +181,11 @@ remove_argument (N, &argc, argv); MaxDensityUse = TRUE; + } + + CLOUD_TOLERANCE = 0.02; + if ((N = get_argument (argc, argv, "-cloud-limit"))) { + remove_argument (N, &argc, argv); + CLOUD_TOLERANCE = atof(argv[N]); + remove_argument (N, &argc, argv); } @@ -189,4 +237,10 @@ } + KEEP_UBERCAL = TRUE; + if ((N = get_argument (argc, argv, "-reset-ubercal"))) { + remove_argument (N, &argc, argv); + KEEP_UBERCAL = FALSE; + } + MIN_ERROR = 0.001; if ((N = get_argument (argc, argv, "-minerror"))) { @@ -247,7 +301,211 @@ if (UpdateAverages && (argc == 1)) return TRUE; - if (UserPatchSelect && (argc == 2)) return TRUE; - if (argc != 3) relphot_usage (); + if (argc != 2) relphot_usage (); return TRUE; } + +int args_client (int argc, char **argv) { + + int N; + double trange; + + // by definition, the client is not parallel + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + HOST_ID = 0; + if ((N = get_argument (argc, argv, "-hostID"))) { + remove_argument (N, &argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOST_ID) relphot_client_usage(); + + HOSTDIR = NULL; + if ((N = get_argument (argc, argv, "-hostdir"))) { + remove_argument (N, &argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOSTDIR) relphot_client_usage(); + + // CATDIR[0] = 0; + // if ((N = get_argument (argc, argv, "-catdir"))) { + // remove_argument (N, &argc, argv); + // strcpy (CATDIR, argv[N]); + // remove_argument (N, &argc, argv); + // } + // if (!CATDIR[0]) relphot_client_usage(); + + IMAGES = NULL; // used in -update mode + BCATALOG = NULL; // used in -load mode + MODE = MODE_NONE; + if ((N = get_argument (argc, argv, "-load"))) { + MODE = MODE_LOAD; + remove_argument (N, &argc, argv); + BCATALOG = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-update-catalogs"))) { + if (MODE) { + fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); + relphot_client_usage(); + } + MODE = MODE_UPDATE; + remove_argument (N, &argc, argv); + IMAGES = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-update-objects"))) { + if (MODE) { + fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); + relphot_client_usage(); + } + MODE = MODE_UPDATE_OBJECTS; + remove_argument (N, &argc, argv); + } + if (!MODE) relphot_client_usage(); + + strcpy (STATMODE, "WT_MEAN"); + if ((N = get_argument (argc, argv, "-statmode"))) { + remove_argument (N, &argc, argv); + strcpy (STATMODE, argv[N]); + remove_argument (N, &argc, argv); + } + + /* specify portion of the sky */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE2 = VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + RESET = FALSE; + if ((N = get_argument (argc, argv, "-reset"))) { + remove_argument (N, &argc, argv); + RESET = TRUE; + } + + KEEP_UBERCAL = TRUE; + if ((N = get_argument (argc, argv, "-reset-ubercal"))) { + remove_argument (N, &argc, argv); + KEEP_UBERCAL = FALSE; + } + + /* define time */ + TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &TSTART)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_dtime (argv[N], &trange)) { + if (!ohana_str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + } + + + MIN_ERROR = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = atof (argv[N]); + remove_argument (N, &argc, argv); + /* require MIN_ERROR > 0 */ + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + UPDATE_CATFORMAT = NULL; + if ((N = get_argument (argc, argv, "-update-catformat"))) { + remove_argument (N, &argc, argv); + UPDATE_CATFORMAT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + AreaSelect = FALSE; + if ((N = get_argument (argc, argv, "-area"))) { + remove_argument (N, &argc, argv); + AreaXmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaXmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaSelect = TRUE; + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; + } + + 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; + } + + if ((MODE == MODE_UPDATE_OBJECTS) && (argc == 1)) return TRUE; + if (argc != 2) relphot_client_usage (); + + return TRUE; +} + Index: trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/bcatalog.c (revision 33650) +++ trunk/Ohana/src/relphot/src/bcatalog.c (revision 33651) @@ -1,46 +1,13 @@ # include "relphot.h" -extern double drand48(); - -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].Nmeasure = average[0].Nmeasure; - averageT[0].measureOffset = average[0].measureOffset; - - // make Nmeasure & measureOffset optional? - - return (TRUE); -} - -int CopyMeasureTiny (MeasureTiny *measureT, Measure *measure) { - - measureT[0].dR = measure[0].dR; - measureT[0].dD = measure[0].dD; - measureT[0].M = measure[0].M; - measureT[0].Mcal = measure[0].Mcal; - measureT[0].dM = measure[0].dM; - measureT[0].airmass = measure[0].airmass; - measureT[0].Xccd = measure[0].Xccd; - measureT[0].Yccd = measure[0].Yccd; - measureT[0].t = measure[0].t; - measureT[0].dt = measure[0].dt; - measureT[0].averef = measure[0].averef; - measureT[0].imageID = measure[0].imageID; - measureT[0].dbFlags = measure[0].dbFlags; - measureT[0].photcode = measure[0].photcode; - - return (TRUE); -} +int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog); int bcatalog (Catalog *subcatalog, Catalog *catalog) { off_t i, j, offset; - int ecode, found, Ns; + int found, Ns, *Nvalid; off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; float mag; - int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf; + int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor; int Nsecfilt = GetPhotcodeNsecfilt (); @@ -53,8 +20,22 @@ ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); ALLOCATE (subcatalog[0].averageT, AverageTiny, NAVERAGE); - ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); Nmeasure = Naverage = 0; - Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = 0; + Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0; + + // copy the following fields to the subcatalog: + subcatalog[0].catID = catalog[0].catID; + subcatalog[0].filename = catalog[0].filename; + + // array to count measures with photcodes->equiv to each of the secfilts + ALLOCATE (Nvalid, int, Nsecfilt); + + // flags used by the photometry analysis (excluding UBERCAL) + unsigned int PHOTOM_FLAGS = + ID_MEAS_NOCAL | // detection ignored for this analysis (photcode, time range) + ID_MEAS_POOR_PHOTOM | // detection is photometry outlier + ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement + ID_MEAS_AREA; // detetion was outside acceptable area of device /* exclude stars not in range or with too few measurements */ @@ -63,13 +44,26 @@ /* start with all stars good */ - CopyAverageTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]); + CopyAverageToTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]); subcatalog[0].averageT[Naverage].measureOffset = Nmeasure; for (j = 0; j < Nsecfilt; j++) { subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j]; - } - + Nvalid[j] = 0; // reset the Nvalid array for this star + } + + // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images) if (RESET) { int Ns; + DVOSecfiltFlags secfiltBits = + ID_SECF_STAR_FEW | + ID_SECF_STAR_POOR | + ID_PHOTOM_PASS_0 | + ID_PHOTOM_PASS_1 | + ID_PHOTOM_PASS_2 | + ID_PHOTOM_PASS_3 | + ID_PHOTOM_PASS_4 | + ID_SECF_USE_SYNTH | + ID_SECF_USE_UBERCAL | + ID_SECF_OBJ_EXT; for (Ns = 0; Ns < Nphotcodes; Ns++) { @@ -77,8 +71,14 @@ int Nsec = GetPhotcodeNsec(thisCode); - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M = NAN; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_FEW; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_POOR; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M = NAN; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map = NAN; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm = NAN; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20 = NAN_S_SHORT; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80 = NAN_S_SHORT; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode = 0; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused = 0; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000; } } @@ -92,8 +92,12 @@ /* select measurements by photcode */ - ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode); + PhotCode *code = GetPhotcodebyCode (catalog[0].measure[offset].photcode); + if (!code) continue; found = FALSE; + int Nsec = -1; // Nsec equivalent for the measurement for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { - if (ecode == photcodes[Ns][0].code) found = TRUE; + if (code->equiv != photcodes[Ns][0].code) continue; + found = TRUE; + Nsec = GetPhotcodeNsec(code->equiv); } if (!found) { @@ -113,7 +117,14 @@ // skip garbage measurements if (catalog[0].measure[offset].psfQual < 0.85) { Npsfqf ++; continue; } - - // check for galaxies - if (!isnan(catalog[0].measure[offset].Map)) { + if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; } + + if (catalog[0].measure[offset].photFlags & code->photomBadMask) { + Nbad++; + continue; + } + if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;} + + // check for galaxies (XXX skip for now) + if (FALSE && !isnan(catalog[0].measure[offset].Map)) { if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) { nEXT ++; @@ -137,13 +148,19 @@ } - CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); + // count this measurement as valid for this secfilt entry + if (Nsec > -1) { + assert (Nsec < Nsecfilt); + Nvalid[Nsec] ++; + } + + CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; subcatalog[0].measureT[Nmeasure].averef = Naverage; if (RESET) { - subcatalog[0].measureT[Nmeasure].Mcal = 0; - subcatalog[0].measureT[Nmeasure].dbFlags &= 0xff00; - subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_PHOTOM; - subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA; - subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; + if (!KEEP_UBERCAL) { + subcatalog[0].measureT[Nmeasure].Mcal = 0; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + } + subcatalog[0].measureT[Nmeasure].dbFlags &= ~PHOTOM_FLAGS; } Nmeasure ++; @@ -156,5 +173,5 @@ // skip object if it is likely to be a galaxy - if (nEXT >= nPSF) { + if (FALSE && (nEXT >= nPSF)) { Nmeasure -= Nm; Ngalaxy ++; @@ -163,9 +180,20 @@ // XXXX test : what checks do I need to make elsewhere to avoid problems here? - if (Nm <= STAR_TOOFEW) { /* enough measurements in band? */ + if (Nm <= STAR_TOOFEW) { /* enough measurements total? */ Nmeasure -= Nm; Nfew ++; continue; } + int anySecfiltGood = FALSE; + for (Ns = 0; Ns < Nsecfilt; Ns++) { /* enough measurements in at least one band? */ + if (Nvalid[Ns] <= STAR_TOOFEW) continue; + anySecfiltGood = TRUE; + } + if (!anySecfiltGood) { + Nmeasure -= Nm; + Nfew ++; + continue; + } + subcatalog[0].averageT[Naverage].Nmeasure = Nm; Naverage ++; @@ -186,5 +214,5 @@ // limit the total number of stars in the catalog if (MaxDensityUse) { - LimitDensityCatalog (subcatalog, catalog); + LimitDensityCatalog_ByNmeasure (subcatalog, catalog); } @@ -192,11 +220,27 @@ fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", subcatalog[0].Naverage, subcatalog[0].Nmeasure, i, catalog[0].filename); - fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d galaxies\n", - Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Ngalaxy); + fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", + Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor); } return (TRUE); } -int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) { +// sort by decreasing Nmeasure (X) +void sort_by_Nmeasure (int *X, off_t *Y, off_t N) { + +# define SWAPFUNC(A,B){ int tmpI; off_t tmpT; \ + tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \ + tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \ +} +# define COMPARE(A,B)(X[A] > X[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) { Catalog tmpcatalog; @@ -206,8 +250,8 @@ int Nsecfilt = GetPhotcodeNsecfilt (); - 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); + gfits_scan (&oldcatalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&oldcatalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&oldcatalog[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); @@ -229,15 +273,13 @@ // we are going to select Nmax entries by generating a random-sorted index list - off_t *index, tmp, i, j, ave; + int *value; + off_t *index, i, j, ave; ALLOCATE (index, off_t, Naverage); + ALLOCATE (value, int, 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; - } + value[i] = subcatalog[0].averageT[i].Nmeasure; + } + sort_by_Nmeasure (value, index, Naverage); // count the number of measurements this selection will yield @@ -266,4 +308,104 @@ Nmeasure ++; } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } + } + + if (VERBOSE) { + char *basename = filebasename (oldcatalog[0].filename); + fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", + Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, oldcatalog[0].Nmeasure, basename); + free (basename); + } + + free (index); + free (value); + free (subcatalog[0].averageT); + free (subcatalog[0].measureT); + free (subcatalog[0].secfilt); + + subcatalog[0].averageT = tmpcatalog.averageT; + subcatalog[0].measureT = tmpcatalog.measureT; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nmax; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt; + + return (TRUE); +} + +int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) { + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + 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].averageT[ave].Nmeasure; + } + + // allocate the output data + ALLOCATE (tmpcatalog.averageT, AverageTiny, Nmax); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); + + 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.averageT[i] = subcatalog[0].averageT[ave]; + tmpcatalog.averageT[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.averageT[i].Nmeasure; j++) { + off_t offset = subcatalog[0].averageT[ave].measureOffset + j; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; + Nmeasure ++; + } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } } @@ -287,60 +429,2 @@ return (TRUE); } - -// for the cases where we are not using a subset of the data, we still need to have a copy of these fields -int populate_tiny_values (Catalog *catalog) { - - 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]); - 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]); - 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; - } - - return (TRUE); -} - -int free_tiny_values (Catalog *catalog) { - - if (catalog[0].averageT) free (catalog[0].averageT); - if (catalog[0].measureT) free (catalog[0].measureT); - return (TRUE); -} - Index: trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- trunk/Ohana/src/relphot/src/global_stats.c (revision 33650) +++ trunk/Ohana/src/relphot/src/global_stats.c (revision 33651) @@ -1,5 +1,5 @@ # include "relphot.h" -void global_stats (Catalog *catalog, int Ncatalog) { +void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD; @@ -17,5 +17,5 @@ int seccode = photcodes[Ns][0].code; - stN = statsStarN (catalog, Ncatalog, Nsec, seccode); + stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr); stX = statsStarX (catalog, Ncatalog, Nsec); stS = statsStarS (catalog, Ncatalog, Nsec); @@ -27,10 +27,10 @@ } - imN = statsImageN (catalog); + imN = statsImageN (catalog, flatcorr); imX = statsImageX (catalog); imM = statsImageM (catalog); imD = statsImagedM (catalog); - msN = statsMosaicN (catalog); + msN = statsMosaicN (catalog, flatcorr); msM = statsMosaicM (catalog); msD = statsMosaicdM (catalog); Index: trunk/Ohana/src/relphot/src/help.c =================================================================== --- trunk/Ohana/src/relphot/src/help.c (revision 33650) +++ trunk/Ohana/src/relphot/src/help.c (revision 33651) @@ -2,6 +2,6 @@ void relphot_usage (void) { - fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n"); - fprintf (stderr, " or: relphot (photcode) -region RA RA DEC DEC\n"); + fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -region RA RA DEC DEC\n"); + fprintf (stderr, " or: relphot (photcodes) -catalog (name)\n"); fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); fprintf (stderr, " use -h for more usage information\n"); @@ -18,6 +18,6 @@ show_help: - fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n"); - fprintf (stderr, " or: relphot (photcode) -region RA RA DEC DEC\n"); + fprintf (stderr, "ERROR: USAGE: relphot (photcode) -region RA RA DEC DEC\n"); + fprintf (stderr, " or: relphot (photcode) -catalog (name)\n"); fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); fprintf (stderr, " options: \n"); @@ -37,4 +37,5 @@ fprintf (stderr, " -imfreeze\n"); fprintf (stderr, " -grid\n"); + fprintf (stderr, " -reset-ubercal : also reset ubercal-ed zero points (otherwise they are sacrosanct)\n"); fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); fprintf (stderr, " -instmag min max\n"); @@ -43,2 +44,33 @@ } +void relphot_client_usage (void) { + fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -load (filename) -hostID (hostID) -hostdir (hostdir) [options]\n"); + fprintf (stderr, " or: relphot -update (filename) -hostID (hostID) -hostdir (hostdir) [options]\n"); + fprintf (stderr, " or: relphot -update-objects -hostID (hostID) -hostdir (hostdir) [options]\n"); + fprintf (stderr, " use -h for more usage information\n"); + exit (2); +} + +void relphot_client_help (int argc, char **argv) { + + /* check for help request */ + if (get_argument (argc, argv, "-help")) goto show_help; + if (get_argument (argc, argv, "-h")) goto show_help; + if (argc == 1) relphot_client_usage(); + return; + +show_help: + fprintf (stderr, "USAGE: relphot_client [-load / -update] (db info)\n\n"); + fprintf (stderr, " relphot_client -load (bcatalog) : extract the bright catalog subset from client's tables\n"); + fprintf (stderr, " (bcatalog) : location where the bright subset is saved\n"); + fprintf (stderr, " relphot_client -update (images) : apply calculated zero points to the client's tables\n\n"); + fprintf (stderr, " (images) : location of the table with the image zero points\n"); + fprintf (stderr, " relphot_client -update-objects : determine average magnitudes for objects\n\n"); + fprintf (stderr, " db info : -hostID (hostID) -hostdir (hostdir) -catdir (catdir)\n"); + fprintf (stderr, "other options:\n"); + fprintf (stderr, " -v : verbose output\n"); + fprintf (stderr, " -vv : extra verbose output\n"); + fprintf (stderr, " \n"); + exit (2); +} + Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 33650) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 33651) @@ -2,6 +2,4 @@ void initialize (int argc, char **argv) { - - int N; relphot_help (argc, argv); @@ -10,48 +8,9 @@ if (!UpdateAverages) { - - N = UserPatchSelect ? 1 : 2; - -# if (0) - // XXX DEP - if ((photcode = GetPhotcodebyName (argv[N])) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); - exit (1); - } - if (photcode[0].type != PHOT_SEC) { - fprintf (stderr, "photcode %s is not a primary or secondary filter\n", argv[N]); - exit (1); - } - // PhotSec is used to select the single average photcode being processed - PhotSec = GetPhotcodeNsec (photcode[0].code); -# endif - - Nphotcodes = 0; - photcodes = NULL; - int NPHOTCODES = 10; - ALLOCATE (photcodes, PhotCode *, NPHOTCODES); - - /* parse the comma-separated list of photcodesKeep */ - char *myList = strcreate(argv[N]); - char *list = myList; - char *codename = NULL; - char *ptr = NULL; - while ((codename = strtok_r (list, ",", &ptr)) != NULL) { - list = NULL; // pass NULL on successive strtok_r calls - fprintf (stderr, "PHOTCODE LIST: %s\n", myList); - fprintf (stderr, "codename: %s\n", codename); - if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); - exit (1); - } - if (photcodes[Nphotcodes][0].type != PHOT_SEC) { - fprintf (stderr, "photcode %s is not an filter type (SEC)\n", codename); - exit (1); - } - Nphotcodes ++; - CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10); - } + // load the list of photcodes into the globals (photcodes, Nphotcodes) + // only a single remaining argument in this mode (the list of photcodes, eg g,r,i) + PhotcodeList = strcreate (argv[1]); + photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes } - // XXX DEP PhotNsec = GetPhotcodeNsecfilt (); if (USE_GRID && (Nphotcodes > 1)) { fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n"); @@ -95,2 +54,45 @@ srand48(B); } + +void initialize_client (int argc, char **argv) { + + // XXX need to determine which globals can affect relphot_client in either mode and pass appropriately + + relphot_client_help (argc, argv); + ConfigInit (&argc, argv); + args_client (argc, argv); + + if (MODE == MODE_UPDATE_OBJECTS) return; + + // load the list of photcodes into the globals (photcodes, Nphotcodes) + PhotcodeList = strcreate (argv[1]); + photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes +} + +void ParsePhotcodeList_old (char *word) { + + Nphotcodes = 0; + photcodes = NULL; + int NPHOTCODES = 10; + ALLOCATE (photcodes, PhotCode *, NPHOTCODES); + + /* parse the comma-separated list of photcodesKeep */ + char *myList = strcreate(word); + char *list = myList; + char *codename = NULL; + char *ptr = NULL; + while ((codename = strtok_r (list, ",", &ptr)) != NULL) { + list = NULL; // pass NULL on successive strtok_r calls + if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); + exit (1); + } + if (photcodes[Nphotcodes][0].type != PHOT_SEC) { + fprintf (stderr, "photcode %s is not an filter type (SEC)\n", codename); + exit (1); + } + Nphotcodes ++; + CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10); + } +} + Index: trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- trunk/Ohana/src/relphot/src/liststats.c (revision 33650) +++ trunk/Ohana/src/relphot/src/liststats.c (revision 33651) @@ -41,4 +41,8 @@ stats[0].min = value[0]; stats[0].max = value[N-1]; + int N80 = MIN (N-1, 0.8*N); + int N20 = MAX (0, 0.2*N); + stats[0].Upper80 = value[N80]; + stats[0].Lower20 = value[N20]; switch (statmode) { @@ -106,5 +110,5 @@ dS += M; } - X2 = X2 / Nm; + X2 = X2 / (Nm - 1); dS = sqrt (dS / Nm); Index: trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/load_catalogs.c (revision 33650) +++ trunk/Ohana/src/relphot/src/load_catalogs.c (revision 33651) @@ -1,11 +1,29 @@ # include "relphot.h" - -Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) { +# include +# include + +// this function loops over the catalogs, loads the data, and extracts a bright subset +// the returned array (catalog, Ncatalog) has the same layout as the full database, but +// only a subset of the detetions & objects + +// if this function is called in parallel mode, it in turn calls load_catalogs_parallel, +// which distributes the work to the remote hosts and loads their results + +// if this function is called with a specified hostID, then only the fraction of the +// database hosted by that hostID is loaded +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int hostID, char *hostpath) { off_t i, Nmeas, Nstar, Nmeas_total, Nstar_total; Catalog *catalog, tcatalog; + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL && !hostID) { + catalog = load_catalogs_parallel (skylist, Ncatalog); + return catalog; + } + if (VERBOSE2) fprintf (stderr, "loading catalog data\n"); + // a bit of an over-alloc, since we don't load all catalogs for a region ALLOCATE (catalog, Catalog, skylist[0].Nregions); @@ -14,8 +32,15 @@ // load data from each region file, only use bright stars for (i = 0; i < skylist[0].Nregions; i++) { + // XXX keep in mind that not all catalogs are loaded dvo_catalog_init (&catalog[i], TRUE); + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + // set up the basic catalog info - tcatalog.filename = skylist[0].filename[i]; + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + tcatalog.filename = hostID ? hostfile : skylist[0].filename[i]; + tcatalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data tcatalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data @@ -34,4 +59,9 @@ } + if (!tcatalog.sorted) { + fprintf (stderr, "this database is not sorted. please sort using addstar -resort\n"); + exit (1); + } + Nstar_total += tcatalog.Naverage; Nmeas_total += tcatalog.Nmeasure; @@ -48,4 +78,9 @@ Nstar += catalog[i].Naverage; Nmeas += catalog[i].Nmeasure; + if ((sizeof(IDX_T) == 8) && (catalog[i].Nmeasure > 0xffffffff)) { + fprintf (stderr, "ERROR: using small-sized IDX_T on data with more than 2G detections per table will cause errors\n"); + fprintf (stderr, " If you need to do this, and you can afford the RAM, rebuild with IDX_T set to off_t (see relphot.h, ImagesOps.c)\n"); + exit (3); + } } if (Nstar < 2) { @@ -54,5 +89,6 @@ 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"); + if (!hostID && (Nstar < 1)) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n"); + // in regular relphot, we shutdown here; in relphot_client, we generate and return an empty table (for consistency) @@ -67,2 +103,136 @@ need to use an XCLD lock here. */ + +// CATDIR is supplied globally +# define DEBUG 1 +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { + + int Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog + + // launch the setphot_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.subset.dat", table->hosts[i].pathname); + + // options / arguments that can affect relphot_client -load: + // VERBOSE, VERBOSE2 + // KEEP_UBERCAL + // RESET (-reset) + // TimeSelect -time + // DophotSelect + // (note that psfQual is applied rigidly at 0.85, as is the galaxy test) + // MAG_LIM + // SIGMA_LIM + // ImagSelect, ImagMin, ImagMax + // MaxDensityUse, MaxDensityValue + + char command[1024]; + snprintf (command, 1024, "relphot_client %s -load %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -D CAMERA %s -D MAG_LIM %f -D SIGMA_LIM %f", + PhotcodeList, table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, CAMERA, MAG_LIM, SIGMA_LIM); + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal", command); strcpy (command, tmpline); } + if (DophotSelect) { snprintf (tmpline, 1024, "%s -dophot %d", command, DophotValue); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); + free (tstart); + free (tstop); + strcpy (command, tmpline); + } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running relphot_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + continue; + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the relphot_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + // each host generates a BrightCatalog structure, with the measure, average, etc value + // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt). I need to + // split out the per-catalog measurements into separate catalog entries. + + // set up an initial array of catalogs + CatalogSplitter *catalogs = BrightCatalogSplitInit (Nsecfilt); + + for (i = 0; i < table->Nhosts; i++) { + + BrightCatalog *bcatalog = NULL; + while ((bcatalog = BrightCatalogLoad (table->hosts[i].results)) == NULL) { + // failed to get the data from this host. This can happen for various reasons. Give the user a chance to try again... + fprintf (stderr, "failed to read data from %s, stopping operations until this can be fixed\n", table->hosts[i].hostname); + fprintf (stderr, "you may run the command manually and send this process the CONT signal\n"); + int pid = getpid(); + kill (pid, SIGSTOP); + fprintf (stderr, "retrying %s\n", table->hosts[i].results); + } + free (table->hosts[i].results); + table->hosts[i].results = NULL; + + BrightCatalogSplit (catalogs, bcatalog); + + free (bcatalog->average); + free (bcatalog->measure); + free (bcatalog->secfilt); + free (bcatalog); + } + + Catalog *catalog = catalogs->catalog; + *Ncatalog = catalogs->Ncatalog; + BrightCatalogSplitFree (catalogs); + + int Nmeasure = 0; + int Naverage = 0; + for (i = 0; i < catalogs->Ncatalog; i++) { + Nmeasure += catalogs->catalog[i].Nmeasure; + Naverage += catalogs->catalog[i].Naverage; + } + + fprintf (stderr, "loaded %d catalogs, using a total of %d stars (%d measures)\n", catalogs->Ncatalog, Naverage, Nmeasure); + + return (catalog); +} Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 33650) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 33651) @@ -9,5 +9,5 @@ // This function generates a subset of the images based on selections. Input db has already // been loaded with the raw fits table data -SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) { +SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region) { Image *image, *subset; @@ -35,5 +35,17 @@ // determine the populated SkyRegions overlapping the requested area - if (RegionSelect) { + + // if the user selects a specific catalog, do that one; + // otherwise limit to the selection region + if (regionName) { + 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); + } else { if (region[0].Rmin > region[0].Rmax) { SkyRegion subregion; @@ -74,14 +86,5 @@ 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); - } + } // save the subset of images in the static reference in ImageOps, set up indexes Index: trunk/Ohana/src/relphot/src/plot_scatter.c =================================================================== --- trunk/Ohana/src/relphot/src/plot_scatter.c (revision 33650) +++ trunk/Ohana/src/relphot/src/plot_scatter.c (revision 33651) @@ -1,5 +1,5 @@ # include "relphot.h" -void plot_scatter (Catalog *catalog, int Ncatalog) { +void plot_scatter (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { off_t i, j, k, m, N, Ntot; @@ -11,5 +11,5 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - Ntot += catalog[i].average[j].Nmeasure; + Ntot += catalog[i].averageT[j].Nmeasure; } } @@ -32,7 +32,7 @@ /* calculate the average value for a single star */ if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - m = catalog[i].average[j].measureOffset; + m = catalog[i].averageT[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { // skip measurements that do not match the current photcode int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); @@ -40,5 +40,5 @@ if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - Mcal = getMcal (m, i); + Mcal = getMcal (m, i, flatcorr, catalog); if (isnan(Mcal)) continue; Mmos = getMmos (m, i); Index: trunk/Ohana/src/relphot/src/plotstuff.c =================================================================== --- trunk/Ohana/src/relphot/src/plotstuff.c (revision 33650) +++ trunk/Ohana/src/relphot/src/plotstuff.c (revision 33651) @@ -165,5 +165,5 @@ graphdata[0].color = black; graphdata[0].lweight = 0; - graphdata[0].size = 0.5; + graphdata[0].size = 1.0; graphdata[0].xmin = dUNDEF; Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 33650) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 33651) @@ -1,3 +1,5 @@ # include "relphot.h" + +# define DEBUG 1 # define TIMESTAMP(TIME) \ @@ -7,5 +9,5 @@ gettimeofday (&start, (void *) NULL); -void reload_catalogs (SkyList *skylist) { +void reload_catalogs (SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath) { int i; @@ -25,10 +27,24 @@ double time7 = 0.0; + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL && !hostID) { + reload_catalogs_parallel (skylist); + return; + } + if (VERBOSE) fprintf (stderr, "re-loading catalog data\n"); /* load data from each region file */ for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + gettimeofday (&start, (void *) NULL); - catalog.filename = skylist[0].filename[i]; + + // set up the basic catalog info + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + catalog.filename = hostID ? hostfile : skylist[0].filename[i]; // only update existing db tables @@ -57,7 +73,8 @@ TIMESTAMP(time2); - populate_tiny_values(&catalog); + populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE); TIMESTAMP(time3a); + // XXX need to worry about the image subset data initImageBins (&catalog, 1, FALSE); initMosaicBins (&catalog, 1, FALSE); @@ -69,7 +86,11 @@ TIMESTAMP(time4); - setMrelFinal (&catalog); + setMrelFinal (&catalog, flatcorr); TIMESTAMP(time5); + // modify the output format as desired (ignore current format on disk) + if (UPDATE_CATFORMAT) { + catalog.catformat = dvo_catalog_catformat (UPDATE_CATFORMAT); + } dvo_catalog_save (&catalog, VERBOSE); dvo_catalog_unlock (&catalog); @@ -85,11 +106,101 @@ } - 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); + fprintf (stderr, "time step 1 %10.3f sec : find catalog\n", time1); + fprintf (stderr, "time step 2 %10.3f sec : load catalog\n", time2); + fprintf (stderr, "time step 3 %10.3f sec : make tiny values\n", time3a); + fprintf (stderr, "time step 4 %10.3f sec : init imbins\n", time3b); + fprintf (stderr, "time step 5 %10.3f sec : find images\n", time4); + fprintf (stderr, "time step 6 %10.3f sec : set Mrel\n", time5); + fprintf (stderr, "time step 7 %10.3f sec : save catalog\n", time6); + fprintf (stderr, "time step 8 %10.3f sec : free catalog\n", time7); } + +// XXX Image to Image Subset +int reload_catalogs_parallel (SkyList *sky) { + + off_t Nimage; + ImageSubset *image = getimages_subset (&Nimage); + + // write out the subset table of image information + char imageFile[512]; + snprintf (imageFile, 512, "%s/Images.subset.dat", CATDIR); + if (!ImageSubsetSave (imageFile, image, Nimage)) { + fprintf (stderr, "failed to write image subset\n"); + exit (1); + } + free (image); + + // now launch the relphot_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + char command[1024]; + snprintf (command, 1024, "relphot_client %s -update-catalogs %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f", + PhotcodeList, imageFile, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR); + + // options & configs which affect relphot_client -update-catalogs + // VERBOSE, VERBOSE2 + // RESET + // TimeSelect + // AreaSelect + // STATMODE + // STAR_TOOFEW + // MIN_ERROR + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (UPDATE) { snprintf (tmpline, 1024, "%s -update", command); strcpy (command, tmpline); } + if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal", command); strcpy (command, tmpline); } + if (UPDATE_CATFORMAT) { snprintf (tmpline, 1024, "%s -update-catformat %s", command, UPDATE_CATFORMAT); strcpy (command, tmpline); } + if (AreaSelect) { snprintf (tmpline, 1024, "%s -area %f %f %f %f", command, AreaXmin, AreaXmax, AreaYmin, AreaYmax); strcpy (command, tmpline); } + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); + free (tstart); + free (tstop); + strcpy (command, tmpline); + } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running relphot_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the relphot_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + return (TRUE); +} Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 33650) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 33651) @@ -1,21 +1,14 @@ # include "relphot.h" -# define USE_DIRECT 0 - -# define TIMESTAMP(TIME) \ + +# define MARKTIME(MSG,...) { \ + float dtime; \ gettimeofday (&stop, (void *) NULL); \ dtime = DTIME (stop, start); \ - TIME += dtime; \ - gettimeofday (&start, (void *) NULL); - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } + fprintf (stderr, MSG, __VA_ARGS__); } int main (int argc, char **argv) { int i, status, Ncatalog; - Catalog *catalog; + Catalog *catalog = NULL; FITS_DB db; struct timeval start, stop; @@ -24,82 +17,4 @@ 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 */ initialize (argc, argv); @@ -107,5 +22,5 @@ /* the object analysis is a separate process iterating over catalogs */ if (UpdateAverages) { - relphot_objects (); + relphot_objects (0, NULL); exit (0); } @@ -125,7 +40,13 @@ // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { - // XXX get ZERO_POINT from config - dvo_image_create (&db, 25.0); - // Shutdown ("ERROR: No images in catalog %s (1)", db.filename); + + // 25.0 is the nominal zero point for measurements in the database + // for all measurements except PHOT_REF: + // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0 + // for measurements using PHOT_REF: + // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known) + dvo_image_create (&db, 25.0); + + // XXX why are we running relphot on an empty database? } else { if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); @@ -133,9 +54,6 @@ MARKTIME("-- load image data: %f sec\n", dtime); - /* load regions and images based on specified sky patch */ - // XXX need to mimic old-style load by passing patch name - // XXX need to reduce number of global variables in use. - // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true - skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect); + /* load regions and images based on specified sky patch and/or catalog */ + skylist = load_images (&db, UserCatalog, &UserPatch); MARKTIME("-- load images: %f sec\n", dtime); @@ -143,37 +61,43 @@ if (!UPDATE) dvo_image_unlock (&db); - /* load catalog data from region files */ - catalog = load_catalogs (skylist, &Ncatalog); - MARKTIME("-- load catalog data: %f sec\n", dtime); + // load the flat correction table (if defined) + char flatcorrfile[256]; + sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); + + if (NLOOP > 0) { + /* load catalog data from region files (hostID is 0 since we are not a client */ + catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); + MARKTIME("-- load catalog data: %f sec\n", dtime); - /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ - - /* match measurements with images, mosaics */ - initImageBins (catalog, Ncatalog, TRUE); - MARKTIME("-- make image bins: %f sec\n", dtime); - - initMosaicBins (catalog, Ncatalog, TRUE); - initGridBins (catalog, Ncatalog); - initMrel (catalog, Ncatalog); - - findImages (catalog, Ncatalog, TRUE); - MARKTIME("-- set up image indexes: %f sec\n", dtime); - - findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ - MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); - - SAVEPLOT = FALSE; - - setExclusions (catalog, Ncatalog); - - global_stats (catalog, Ncatalog); - - if (PLOTSTUFF) { - plot_star_coords (catalog, Ncatalog); - // plot_mosaic_fields (catalog); - } - - // if we are measuring the flat-field correction grid, we need to perform a number of iterations first: - if (USE_GRID) { + /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("-- make image bins: %f sec\n", dtime); + + initMosaicBins (catalog, Ncatalog, TRUE); + initGridBins (catalog, Ncatalog); + initMrel (catalog, Ncatalog); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("-- set up image indexes: %f sec\n", dtime); + + findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); + + SAVEPLOT = FALSE; + + setExclusions (catalog, Ncatalog, TRUE); + + global_stats (catalog, Ncatalog, flatcorr); + + if (PLOTSTUFF) { + plot_star_coords (catalog, Ncatalog); + // plot_mosaic_fields (catalog); + } + + // if we are measuring the flat-field correction grid, we need to perform a number of iterations first: + if (USE_GRID) { int star_toofew; @@ -183,23 +107,52 @@ STAR_TOOFEW = 0; for (i = 0; i < NGRID; i++) { - STAR_BAD = ID_STAR_POOR; - setMrel (catalog, Ncatalog); - STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; - setMgrid (catalog); + STAR_BAD = ID_STAR_POOR; + setMrel (catalog, Ncatalog, flatcorr); + STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + setMgrid (catalog, flatcorr); } STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; STAR_TOOFEW = star_toofew; - } - - /* determine fit values */ - for (i = 0; i < NLOOP; i++) { - setMrel (catalog, Ncatalog); - setMcal (catalog, FALSE); - setMmos (catalog, FALSE); - setMgrid (catalog); + } + + /* determine fit values */ + for (i = 0; i < NLOOP; i++) { + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } + setMrel (catalog, Ncatalog, flatcorr); + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } + setMcal (catalog, FALSE, flatcorr); + setMmos (catalog, FALSE, flatcorr); + setMgrid (catalog, flatcorr); + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + plot_grid (catalog, flatcorr); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (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, flatcorr); + 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, flatcorr); + MARKTIME("-- finished loop %d: %f sec\n", i, dtime); + } + if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog); - plot_grid (catalog); + plot_scatter (catalog, Ncatalog, flatcorr); + plot_grid (catalog, flatcorr); plot_mosaics (); plot_images (); @@ -207,109 +160,99 @@ plot_chisq (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); - 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); - MARKTIME("-- finished loop %d: %f sec\n", i, dtime); - } - - 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 (); - - /* set Mcal & Mmos for bad images */ - setMcal (catalog, TRUE); - setMmos (catalog, TRUE); - 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); - } + if (USE_GRID) dump_grid (); + + /* set Mcal & Mmos for bad images */ + setMcal (catalog, TRUE, flatcorr); + setMmos (catalog, TRUE, flatcorr); + 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) { + // 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 + } + 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; + + // save Images.dat using the copied structure + if (UPDATE_CATFORMAT) { + // ensure the db format is updated + dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + char photcodeFile[1024]; + sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR); + SavePhotcodesFITS (photcodeFile); + } + dvo_image_save (&dbX, VERBOSE); + dvo_image_unlock (&dbX); + MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); + } + + /* at this point, we have correct cal coeffs in the image/mosaic structures */ + for (i = 0; i < Ncatalog; i++) { + free_tiny_values (&catalog[i]); + dvo_catalog_free (&catalog[i]); + } + freeImageBins (Ncatalog, TRUE); + freeMosaicBins (Ncatalog, TRUE); + freeGridBins (Ncatalog); + + } // (NLOOP > 0) : this loop determines the offsets per chip + + reload_images (&db); + + // only change the real database files 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++) { - free_tiny_values (&catalog[i]); - dvo_catalog_free (&catalog[i]); - } - freeImageBins (Ncatalog, TRUE); - freeMosaicBins (Ncatalog, TRUE); - freeGridBins (Ncatalog); - /* load catalog data from region files, update Mrel include all data */ - reload_catalogs (skylist); + reload_catalogs (skylist, flatcorr, 0, NULL); MARKTIME("-- updated all catalogs: %f sec\n", dtime); + if (UPDATE_CATFORMAT) { + // ensure the db format is updated + db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + } dvo_image_update (&db, VERBOSE); dvo_image_unlock (&db); Index: trunk/Ohana/src/relphot/src/relphot_client.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_client.c (revision 33651) +++ trunk/Ohana/src/relphot/src/relphot_client.c (revision 33651) @@ -0,0 +1,84 @@ +# include "relphot.h" + +// relphot_client is run on a remote host and is responsible for updating the catalogs +// owned by that host. + +// there are two modes: + +// relphot_client -load : extract the bright catalog subset from the client's tables and +// save to a local FITS table for relphot to load and used + +// (relphot loads the bright catalog subsets from all clients, then uses this to determine +// the per-image zero points (and potentially the flat-field corrections) + +// relphot_client -update : load image table containing the update zero point information +// and apply this to the client's tables + +int main (int argc, char **argv) { + + // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc) + initialize_client (argc, argv); + + // load the current sky table (layout of all SkyRegions) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); + if (!sky) { + fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); + exit (2); + } + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + if (!skylist) { + fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); + exit (2); + } + + switch (MODE) { + case MODE_LOAD: { + int Ncatalog; + Catalog *catalog = load_catalogs (skylist, &Ncatalog, HOST_ID, HOSTDIR); + if (!catalog) { + fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); + exit (2); + } + BrightCatalog *bcatalog = BrightCatalogMerge (catalog, Ncatalog); + if (!BrightCatalogSave (BCATALOG, bcatalog)) { + fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); + exit (2); + } + break; + } + + case MODE_UPDATE: { + // load the image subset table from the specified location + off_t Nimage; + ImageSubset *image = ImageSubsetLoad (IMAGES, &Nimage); + if (!image) { + fprintf (stderr, "ERROR loading image subset %s\n", CATDIR); + exit (2); + } + + // save the available image information in the static array in ImageOps.c + initImagesSubset (image, NULL, Nimage); + + // load the flat-field correction table from CATDIR + char flatcorrFile[1024]; + snprintf (flatcorrFile, 1024, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrFile, VERBOSE); + + reload_catalogs (skylist, flatcorr, HOST_ID, HOSTDIR); + break; + } + + case MODE_UPDATE_OBJECTS: { + relphot_objects (HOST_ID, HOSTDIR); + break; + } + + default: + fprintf (stderr, "impossible!"); + abort(); + } + + exit (0); +} Index: trunk/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_objects.c (revision 33650) +++ trunk/Ohana/src/relphot/src/relphot_objects.c (revision 33651) @@ -1,5 +1,12 @@ # include "relphot.h" -int relphot_objects () { +// apply the current calibration values to the measurements to determine the average +// magnitudes, applying the desired clipping analysis. Note that this does not depend on +// the form of the image zero point measurement. If there is a grid or map term, and this +// has been correctly propagated to the measurement, the average will respect that value. + +int relphot_objects_parallel (SkyList *sky); + +int relphot_objects (int hostID, char *hostpath) { off_t i, j, k; @@ -17,9 +24,21 @@ skylist = SkyListByPatch (sky, -1, &UserPatch); + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL && !hostID) { + relphot_objects_parallel (skylist); + return TRUE; + } + // load data from each region file, only use bright stars for (i = 0; i < skylist[0].Nregions; i++) { + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + // set up the basic catalog info - catalog.filename = skylist[0].filename[i]; + // set up the basic catalog info + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + catalog.filename = hostID ? hostfile : skylist[0].filename[i]; catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data @@ -38,18 +57,32 @@ } - // reset + // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images) if (RESET) { Nsecfilt = catalog.Nsecfilt; + DVOSecfiltFlags secfiltBits = + ID_SECF_STAR_FEW | + ID_SECF_STAR_POOR | + ID_PHOTOM_PASS_0 | + ID_PHOTOM_PASS_1 | + ID_PHOTOM_PASS_2 | + ID_PHOTOM_PASS_3 | + ID_PHOTOM_PASS_4 | + ID_SECF_USE_SYNTH | + ID_SECF_USE_UBERCAL | + ID_SECF_OBJ_EXT; + for (j = 0; j < catalog.Naverage; j++) { catalog.average[j].flags = 0; for (k = 0; k < Nsecfilt; k++) { - catalog.secfilt[j*Nsecfilt + k].M = NAN; - catalog.secfilt[j*Nsecfilt + k].M_20 = NAN_S_SHORT; - catalog.secfilt[j*Nsecfilt + k].M_80 = NAN_S_SHORT; - catalog.secfilt[j*Nsecfilt + k].dM = NAN; - catalog.secfilt[j*Nsecfilt + k].Xm = NAN_S_SHORT; - catalog.secfilt[j*Nsecfilt + k].Ncode = 0; - catalog.secfilt[j*Nsecfilt + k].Nused = 0; - // XXX reset the photometry flags for secfilt entries? + catalog.secfilt[j*Nsecfilt + k].M = NAN; + catalog.secfilt[j*Nsecfilt + k].Map = NAN; + catalog.secfilt[j*Nsecfilt + k].dM = NAN; + catalog.secfilt[j*Nsecfilt + k].Xm = NAN; + catalog.secfilt[j*Nsecfilt + k].M_20 = NAN_S_SHORT; + catalog.secfilt[j*Nsecfilt + k].M_80 = NAN_S_SHORT; + catalog.secfilt[j*Nsecfilt + k].Ncode = 0; + catalog.secfilt[j*Nsecfilt + k].Nused = 0; + catalog.secfilt[j*Nsecfilt + k].flags &= ~secfiltBits; + catalog.secfilt[j*Nsecfilt + k].ubercalDist = 1000; } } @@ -76,2 +109,79 @@ return (TRUE); } + +// CATDIR is supplied globally +# define DEBUG 1 +int relphot_objects_parallel (SkyList *sky) { + + // launch the setphot_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + // options / arguments that can affect relphot_client -update-objects: + // VERBOSE, VERBOSE2 + // KEEP_UBERCAL + // RESET (-reset) + // TimeSelect -time + // DophotSelect + // (note that psfQual is applied rigidly at 0.85, as is the galaxy test) + // MAG_LIM + // SIGMA_LIM + // ImagSelect, ImagMin, ImagMax + // MaxDensityUse, MaxDensityValue + + char command[1024]; + snprintf (command, 1024, "relphot_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR); + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (UPDATE) { snprintf (tmpline, 1024, "%s -update", command); strcpy (command, tmpline); } + if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal", command); strcpy (command, tmpline); } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running relphot_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the relphot_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- trunk/Ohana/src/relphot/src/select_images.c (revision 33650) +++ trunk/Ohana/src/relphot/src/select_images.c (revision 33651) @@ -25,5 +25,5 @@ Image *image; - off_t i, j, k, m, nStart, iSky, nimage, NIMAGE; + off_t i, j, k, m, nStart, iSky, nimage, NIMAGE, D_NIMAGE; off_t *line_number; int InRange, ecode, found; @@ -63,4 +63,6 @@ DminSkyRegion = +90.0; DmaxSkyRegion = -90.0; + + D_NIMAGE = 6000; // FILE *ftest = fopen ("relphot.dump.dat", "w"); @@ -110,5 +112,5 @@ nimage = 0; - NIMAGE = 100; + NIMAGE = D_NIMAGE; ALLOCATE (image, Image, NIMAGE); ALLOCATE (line_number, off_t, NIMAGE); @@ -117,4 +119,6 @@ for (i = 0; i < Ntimage; i++) { + if (!(i % 300000)) fprintf (stderr, "."); + /* exclude images by photcode */ ecode = GetPhotcodeEquivCodebyCode (timage[i].photcode); @@ -137,7 +141,4 @@ continue; } - - // 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 */ @@ -148,5 +149,4 @@ 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; @@ -155,5 +155,4 @@ 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; @@ -230,23 +229,24 @@ 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 */ + /* always allow 'few' images to succeed, if possible (new images / detections may have + * been added) */ if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR); } if (RESET) { - // XXX assignMcal (&image[nimage], (double *) NULL, -1); - // XXX this needs to be thought through a bit more - image[nimage].Mcal = 0.0; - image[nimage].dMcal = NAN; + if (!KEEP_UBERCAL) { + image[nimage].Mcal = 0.0; + image[nimage].dMcal = NAN; + image[nimage].flags &= ~ID_IMAGE_PHOTOM_UBERCAL; + } image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR; + image[nimage].ubercalDist = 1000; } line_number[nimage] = i; nimage ++; if (nimage == NIMAGE) { - NIMAGE += 100; + NIMAGE += D_NIMAGE; + D_NIMAGE = MAX (100000, D_NIMAGE * 1.5); REALLOCATE (image, Image, NIMAGE); REALLOCATE (line_number, off_t, NIMAGE); @@ -262,4 +262,6 @@ REALLOCATE (line_number, off_t, MAX (nimage, 1)); free (skycoords); + free (RmaxSky); + free (index); *Nimage = nimage; Index: trunk/Ohana/src/relphot/src/setExclusions.c =================================================================== --- trunk/Ohana/src/relphot/src/setExclusions.c (revision 33650) +++ trunk/Ohana/src/relphot/src/setExclusions.c (revision 33651) @@ -6,5 +6,5 @@ // and uses the bits read from disk as the test -int setExclusions (Catalog *catalog, int Ncatalog) { +int setExclusions (Catalog *catalog, int Ncatalog, int verbose) { off_t i, j, k, m, Narea, Nnocal, Ngood; @@ -59,7 +59,7 @@ } } - if (VERBOSE) fprintf (stderr, OFF_T_FMT" measurements marked by area\n", Narea); - if (VERBOSE) fprintf (stderr, OFF_T_FMT" measurements marked nocal\n", Nnocal); - if (VERBOSE) fprintf (stderr, OFF_T_FMT" measurements kept for analysis\n", Ngood); + if (verbose) fprintf (stderr, OFF_T_FMT" measurements marked by area\n", Narea); + if (verbose) fprintf (stderr, OFF_T_FMT" measurements marked nocal\n", Nnocal); + if (verbose) fprintf (stderr, OFF_T_FMT" measurements kept for analysis\n", Ngood); return (TRUE); } Index: trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 33650) +++ trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 33651) @@ -5,14 +5,31 @@ // output dbFlags values -void setMrelFinal (Catalog *catalog) { - - off_t i, j, m; +void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) { + + off_t i; int ecode; + + int Nsecfilt = GetPhotcodeNsecfilt (); /* if we reset the catalog, reset all the current measurements */ if (RESET) { - - int Nsecfilt = GetPhotcodeNsecfilt (); - + // flags used by the photometry analysis (excluding UBERCAL) + unsigned int PHOTOM_FLAGS = + ID_MEAS_NOCAL | // detection ignored for this analysis (photcode, time range) + ID_MEAS_POOR_PHOTOM | // detection is photometry outlier + ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement + ID_MEAS_AREA; // detetion was outside acceptable area of device + + // flags used by the photometry analysis (excluding UBERCAL) + // unsigned int secfiltFlags = + // ID_PHOTOM_PASS_0 | // average measured at pass 0 + // ID_PHOTOM_PASS_1 | // average measured at pass 1 + // ID_PHOTOM_PASS_2 | // average measured at pass 2 + // ID_PHOTOM_PASS_3 | // average measured at pass 3 + // ID_PHOTOM_PASS_4 | // average measured at pass 3 + // ID_SECF_USE_SYNTH | // average measured at pass 3 + // ID_SECF_USE_UBERCAL | // average measured at pass 3 + // ID_SECF_OBJ_EXT; // average measured at pass 3 + int Ns; for (Ns = 0; Ns < Nphotcodes; Ns++) { @@ -22,9 +39,18 @@ for (i = 0; i < catalog[0].Naverage; i++) { - catalog[0].secfilt[Nsecfilt*i+Nsec].M = NAN; - catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN; - catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT; - - m = catalog[0].average[i].measureOffset; + off_t N = Nsecfilt*i+Nsec; + catalog[0].secfilt[N].M = NAN; + catalog[0].secfilt[N].Map = NAN; + catalog[0].secfilt[N].dM = NAN; + catalog[0].secfilt[N].Xm = NAN; + catalog[0].secfilt[N].M_20 = NAN_S_SHORT; + catalog[0].secfilt[N].M_80 = NAN_S_SHORT; + catalog[0].secfilt[N].Ncode = 0; + catalog[0].secfilt[N].Nused = 0; + catalog[0].secfilt[N].flags = 0; // XXX there are no *astrometry* bits in secfilt.flags + catalog[0].secfilt[N].ubercalDist = 1000; + + off_t m = catalog[0].average[i].measureOffset; + off_t j; for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { @@ -39,10 +65,9 @@ } - catalog[0].measure[m].Mcal = 0; - catalog[0].measure[m].dbFlags &= 0xff00; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_PHOTOM; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_AREA; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL; + if (!KEEP_UBERCAL) { + catalog[0].measure[m].Mcal = 0; + catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + } + catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS; } } @@ -51,32 +76,31 @@ // this sets flags in the measureT element, not the measure element - setExclusions (catalog, 1); /* mark by area */ + setExclusions (catalog, 1, VERBOSE); /* mark by area */ /* set catalog[0].found[i] = FALSE */ - ALLOCATE (catalog[0].found, off_t, MAX (1, catalog[0].Naverage)); - for (i = 0; i < catalog[0].Naverage; i++) { + ALLOCATE (catalog[0].found, off_t, MAX (1, Nsecfilt*catalog[0].Naverage)); + for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) { catalog[0].found[i] = FALSE; } + clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ for (i = 0; i < 5; i++) { - skip_measurements (catalog, i); /* set ID_MEAS_SKIP for measures to be skipped */ - setMrelOutput (catalog, 1, FALSE); /* set Mrel using all measures not skipped */ - clean_measures (catalog, 1, TRUE); /* mark outliers ID_MEAS_POOR_PHOTOM */ - setMrelOutput (catalog, 1, TRUE); /* set Mrel using remaining measures */ - } - setMcalOutput (catalog, 1); - - int Nsecfilt = GetPhotcodeNsecfilt (); + skip_measurements (catalog, i, flatcorr); /* set ID_MEAS_SKIP for measures to be skipped */ + setMrelOutput (catalog, 1, i, flatcorr); /* set Mrel using remaining measures */ + } + setMcalOutput (catalog, 1, flatcorr); /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */ /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */ + // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found) for (i = 0; i < catalog[0].Naverage; i++) { catalog[0].average[i].flags &= ~ID_STAR_FEW; catalog[0].average[i].flags &= ~ID_STAR_POOR; + off_t j; for (j = 0; j < Nsecfilt; j++) { catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_FEW; catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_POOR; } - m = catalog[0].average[i].measureOffset; + off_t m = catalog[0].average[i].measureOffset; for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL; @@ -86,83 +110,143 @@ /* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */ -void skip_measurements (Catalog *catalog, int pass) { +// For each star & photcode, we are selecting the measurements to use in the average +// photometry value. This function is called several times, lowering the bar to +// acceptance on each pass. The rejections at the different passes are: + +// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD) +// 1 : accept measurements thought to be POOR based on photflags +// 2 : accept the measurements marked as outliers; accept images marked as outliers +// 3 : accept measurements thought to be BAD based on photflags +// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT) + +// XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all +// non-outliers with Imag range + +// int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) { +// +// off_t k; +// +// int Nsecfilt = GetPhotcodeNsecfilt (); +// +// off_t m = average[0].measureOffset; +// +// for (k = 0; k < average[0].Nmeasure; k++, m++) { +// fprintf (stderr, "meas: %08x\n", measure[m].dbFlags); +// } +// +// int Ns; +// for (Ns = 0; Ns < Nsecfilt; Ns++) { +// fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags); +// } +// return 1; +// } + +int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure); + +void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) { off_t i, k, m; - int ecode, found, Ns; off_t Ntot, Ntry, Nkeep, Nskip; float mag; + int Nsecfilt = GetPhotcodeNsecfilt (); + Ntot = Ntry = Nskip = Nkeep = 0; /* allow measures from images marked POOR and FEW */ - if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; + if (pass >= 2) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; /* allow measures marked as outliers (POOR) and off image region (AREA) */ - if (pass >= 3) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM; - - /* mark measures which should be ignored on second pass */ + if (pass >= 2) { + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM; + } else { + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; + } + + /* mark measures which should be ignored on this pass */ for (i = 0; i < catalog[0].Naverage; i++) { Ntot += catalog[0].average[i].Nmeasure; - if (catalog[0].found[i]) continue; - - m = catalog[0].average[i].measureOffset; - for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) { - Ntry++; - - /* clear SKIP for all measures at first */ - catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - - /** never use these measurements (wrong photcode, bad time range) */ - /* skipped via NOCAL, don't mark as skipped */ - ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode); - found = FALSE; - for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { - if (ecode == photcodes[Ns][0].code) found = TRUE; + + if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) { + fprintf (stderr, "test obj\n"); + print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure); + } + + // mark measurements for each secfilt separately + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + /* star/photcodes already calibrated */ + if (catalog[0].found[Nsecfilt*i+Nsec]) continue; + + m = catalog[0].average[i].measureOffset; + for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) { + Ntry++; + + // skip measurements not related to this photcode + PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode); + if (!code) continue; + if (code->equiv != thisCode) continue; + + // skip measurements by time range (mark as skipped for this secfilt) + // XXX note that this is a bit dangerous : some objects may never get calibrated + if (TimeSelect) { + if (catalog[0].measure[m].t < TSTART) goto skip; + if (catalog[0].measure[m].t > TSTOP) goto skip; + } + + // we now have a measurement relavant to this photcode (and time range) + + // clear SKIP for all measures at first + catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + + // skip measurements from BAD images and mosaics (not REF mags) + // do NOT skip measurements without a matching image (REF mags) + off_t Nim = getImageEntry (m, 0); + if (Nim > -1) { + if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip; + if (isnan(getMmos (m, 0))) goto skip; + } + + // PASS 4 : skip measurements by inst mag limit (not REF mags) + if ((pass < 4) && ImagSelect) { + if (Nim > -1) { + mag = PhotInst (&catalog[0].measure[m]); + if (mag < ImagMin) goto skip; + if (mag > ImagMax) goto skip; + } + } + + // PASS 3 : accept bad measurements (eg, SAT, CR) + if (pass < 3) { + if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip; + if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) { + if (catalog[0].measure[m].psfQual < 0.85) goto skip; + } + } + + // PASS 2 : internal outliers accepted + + // PASS 1 : accept poor measurements as well (eg, POOR FIT, etc) + if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip; + + // PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc) + + Nkeep ++; + + continue; + + skip: + catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; + catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; + Nskip ++; } - if (!found) continue; - - /* skip measurements by time range */ - if (TimeSelect) { - if (catalog[0].measure[m].t < TSTART) goto skip; - if (catalog[0].measure[m].t > TSTOP) goto skip; - } - - /* skip measurements with sat. dophot values */ - // XXX we need to fix this: use the photFlags to ID saturation... - // if ((pass < 4) && (catalog[0].measure[m].dophot == 10)) goto skip; - - /* skip measurements from BAD images and mosaics */ - /* do NOT skip measurements without a matching image */ - if (isnan(getMcal (m, 0))) goto skip; - if (isnan(getMmos (m, 0))) goto skip; - - /* skip measurements by inst mag limit */ - if ((pass < 4) && ImagSelect) { - mag = PhotInst (&catalog[0].measure[m]); - if (mag < ImagMin) goto skip; - if (mag > ImagMax) goto skip; - } - - // d1 = (catalog[0].measure[m].dophot == 1); - // d2 = (catalog[0].measure[m].dophot == 2); - // fprintf (stderr, "fix the dophot / type selections!"); - // abort (); - - /* skip measurements by measurement error */ - // if ((pass < 2) && !(d1 || d2)) goto skip; - - /* if ((pass < 1) && (catalog[0].measure[m].dM > SIGMA_LIM)) goto skip; */ - // if ((pass < 1) && !d1) goto skip; - Nkeep ++; - - continue; - - skip: - catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; - catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; - Nskip ++; - } - } - if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n", + } + } + if (VERBOSE2) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n", pass, Ntot, Ntry, Nskip, Nkeep); }