diff --git a/README.md b/README.md index 8513389..3cabf9d 100644 --- a/README.md +++ b/README.md @@ -20,19 +20,20 @@ git clone git@github.com:sweverett/Balrog-GalSim.git ## Running Balrog -Balrog is run with a python call to the script `balrog_injection.py`, with a few extra (required and optional) input arguments: +(Updated as of 4/5/18) Balrog is run with a python call to the script `balrog_injection.py`, with a few extra (required and optional) input arguments that can either be passed by command line: ``` -python balrog_injection.py [config_file] [tile_list] [geom_file] [-t tile_dir] [-c config_dir] \ +python balrog_injection.py [config_file] [-l tile_list] [-g geom_file] [-t tile_dir] [-c config_dir] \ [-p psf_dir] [-o output_dir] [-v --verbose] ``` +or explicitly in the `config_file` described below: -* `config_file`: A GalSim config file (for now .yaml, but more types in future) that defines universal parameters for all injections. Chip-specific config parameters are appended to this file during processing. Each DES tile will produce a separate yaml config that handles all injections for all chips in all bands that overlap with the tile ([see GalSim's Demo 6](https://github.com/GalSim-developers/GalSim/blob/master/examples/demo6.yaml) to see how appended configs work / look). An example balrog config file is given in `configs/bal_config.yaml`. +* `config_file`: A GalSim config file (for now .yaml, but more types in future) that defines universal parameters for all injections. Chip-specific config parameters are appended to this file during processing. Each DES tile will produce a separate yaml config that handles all injections for all chips in all bands that overlap with the tile ([see GalSim's Demo 6](https://github.com/GalSim-developers/GalSim/blob/master/examples/demo6.yaml) to see how appended configs work / look). An example balrog config file is given in `configs/bal_config.yaml`. It also includes a description of how to pass these command-line input arguments in the config itself. * `tile_list`: A txt or csv file that contains a list of all tiles that are to be processed. Delimiter can be commas or newlines. * `geom_file`: The fits file containing tile geometry (e.g. `Y3A2_COADDTILE_GEOM.fits`). * `tile_dir`: Directory location that contains all desired DES tile files/folders (i.e. for N tiles, the location would contain at least the N directories `/DES0001-0001`, `/DES0001-0002`, etc.). Set to `.` by default. -* `config_dir`: Directory location of the global GalSim config file. tile list file, and geometry file if not given in inputted filenames (the files only must be in the same directory if this is passed). Set to `.` by default. +* `config_dir`: Directory location of the global GalSim config file. tile list file, and geometry file if not given in inputted filenames (the files only must be in the same directory if this is passed). Set to `.` by default. Cannot be set in `config_file` (for obvious reasons). * `psf_dir`: Relative directory path of a tile's psf *from* a given tile location (e.g. `{tile_dir}/{TILE}/{psf_dir}`). Set to `psfs` by default. * `output_dir`: Location of parent output directory for Balrog images and config files. Images are saved as `{output_dir}/balrog_images/{realization}/{tilename}/{band}/{chipname}_balrog_inj.fits`, and tile configs are saved to `{output_dir}/configs/bal_config_{realization}_{tilename}.yaml`. * `verbose`: Use -v for more verbose messages. @@ -106,6 +107,12 @@ python balrog/balrog_injection.py config/bal_config.yaml inputs/tilelist.csv inp -t inputs/tiles -o outputs/ ``` +Alternatively, these values could have been set in the `image` field of the `config_file` (except for `output_dir`; see example file) in which case you would simply type + +``` +python balrog/balrog_injection.py config/bal_config.yaml +``` + ## Balrog Config The global balrog config file is very similar to yaml config files used for GalSim. Check out the [GalSim demo page](https://github.com/GalSim-developers/GalSim/wiki/Tutorials) for lots of examples of how to use the yaml files and how they are translated to python scripts. The important parts are summarized in the comments of the example config [bal_config.yaml](https://github.com/sweverett/Balrog-GalSim/blob/master/config/bal_config.yaml). (**Note**: You will need to change a few path names in the config file to work for your local machine!) @@ -113,11 +120,11 @@ The global balrog config file is very similar to yaml config files used for GalS However, note that a Balrog config will **not** run successfully if called by the `galsim` executable; it only houses the global simulation variables while the individual chip injections parameters are set during the `balrog_injection.py` script. Each tile produces its own complete multi-output yaml config file that is sent to `galsim` at the end of processing. There are a few config inputs specific to the `image` field of Balrog configs that are worth highlighting here: -* `n_galaxies`: The total number of galaxies to be injected per DES tile. -* `gal_density`: The final injected Balrog galaxy density in the tile field. +* `n_galaxies`: The total number of galaxies to be injected per DES tile **per realization**. +* `gal_density`: The injected Balrog galaxy density in the tile field **per realization**. * `n_realizations`: The number of injection realizations used to reach the desired galaxy count or density. -Two things to note: (1) **Only one** of `n_galaxies` or `gal_density` is allowed as an input; not both! (2) Either input should give the desired **final** count or density. The script will calculate the injection number per realization given the input and the total number of realizations to be computed. +Two things to note: (1) **Only one** of `n_galaxies` or `gal_density` is allowed as an input; not both! (2) Either input should give the desired count or density **per realization**!. An older version of the code had this set to the desired final count/desnity, but this was counter-intuitive for users. ## Input Catalogs diff --git a/balrog/balrog_injection.py b/balrog/balrog_injection.py index d91f891..d9e89f1 100644 --- a/balrog/balrog_injection.py +++ b/balrog/balrog_injection.py @@ -26,11 +26,7 @@ import argparse import datetime import itertools -# TODO: Have automatic check for astropy vs. fitsio! -# try: import fitsio import fitsio -# except: import astropy -# from fitsio import FITS, FITSHDR, etc. from astropy.io import fits from astropy import wcs from astropy.table import Table @@ -43,9 +39,9 @@ #------------------------------------------------------------------------------- # Urgent todo's: -# TODO: Implement error handling for galaxy injections / gsparams! (I think fixed now) +# TODO: Implement error handling for galaxy injections / gsparams! (Working solution, but need +# to look into more detail per Erin) # TODO: Clean up evals in add_gs_injection()! -# TODO: Change n_galaxies to be divided up per realization! # TODO: Figure out injector.py parameter parsing issue! # TODO: Check pixel origin for transformations! # TODO: Remove all stars not contained in unique region! @@ -53,18 +49,12 @@ # Some extra todo's: # TODO: Allow grid parameters to be passed! # TODO: Have grid correctly handle chips w/ no injections on them -# TODO: Use fitsio when available! -# TODO: Get rid of extra '/' on config file parsing! +# TODO: Have code automatically check fitsio vs astropy # TODO: Implement some print statements as warnings -# TODO: Redistribute config.galaxies_remainder among first m delta_mag_bds[i]) & (np.abs(delta_mag) <= delta_mag_bds[i+1]) + ax1.plot(meas_matched['cm_T'][these],delta_mag[these],'.',alpha=0.5,color='k') + ax1.set_xscale('log') + ax1.set_xlabel('cm_T (r)') + ax1.set_ylabel('cm_mag Measured - True (r)') + + outfile = os.path.join(outdir, 't_vs_mag_dif'+bands[bandind]) + fig.savefig(outfile) + plt.close(fig) # Fun color-coded magnitude plot. fig,ax1 = plt.subplots(nrows=1,ncols=1,figsize=(7,7)) @@ -251,7 +279,7 @@ def make_plots(truth = None, meas_matched = None, truth_matched = None, sep = No outfile = os.path.join(outdir, 't_vs_mag_color-coded-by-delta-mag'+bands[bandind]) fig.savefig(outfile) - fig.close() + plt.close(fig) # Do matching errors matter? fig,ax1 = plt.subplots(figsize=(7,7)) @@ -271,7 +299,7 @@ def make_plots(truth = None, meas_matched = None, truth_matched = None, sep = No outfile = os.path.join(outdir, 'sep_vs_t') fig.savefig(outfile) - fig.close() + plt.close(fig) fig,ax1 = plt.subplots(figsize=(7,7)) #small = (truth_matched['cm_T'] / truth_matched['cm_T_err'] ) > 10. @@ -290,7 +318,7 @@ def make_plots(truth = None, meas_matched = None, truth_matched = None, sep = No outfile = os.path.join(outdir, 'Terr') fig.savefig(outfile) - fig.close() + plt.close(fig) fig,ax1 = plt.subplots(figsize=(7,7)) for i in xrange(n_delta_bins): @@ -305,7 +333,136 @@ def make_plots(truth = None, meas_matched = None, truth_matched = None, sep = No outfile = os.path.join(outdir, 'logsb_vs_cm_T') fig.savefig(outfile) - fig.close() + plt.close(fig) + + # Shape and Correlation plots + # fig, ax = ... + # pudb.set_trace() + + if bandind == 1: + # only do once! + + # min/max separation and bin size + mins, maxs = 0.1, 10 + bs = 0.075 + + # First do the truth catalog + g1, g2 = truth_matched['cm_g'][:,0], truth_matched['cm_g'][:,1] + e1, e2 = ngmix.shape.g1g2_to_e1e2(g1, g2) + + # Easiest to use treecorr by making a new temporary catalog + truth_outfile = os.path.join(outdir,'treecorr_temp_file_truth.fits') + delete_file(truth_outfile) + fits = fitsio.FITS(truth_outfile,'rw') + data = [truth_matched['ra'], truth_matched['dec'], g1, g2, e1, e2] + names = ['ra', 'dec', 'g1', 'g2', 'e1', 'e2'] + fits.write(data, names=names) + fits.close() + + cat = treecorr.Catalog(truth_outfile, ra_col='ra', dec_col='dec', + ra_units='degrees', dec_units='degrees', + g1_col='e1', g2_col='e2') + gg = treecorr.GGCorrelation(min_sep=mins, max_sep=maxs, bin_size=bs, + sep_units='arcmin') + gg.process(cat) + fig = plot_gg_corr(gg, plt_type='Truth') + outfile = os.path.join(outdir, 'gg_corr_truth.png') + fig.savefig(outfile) + plt.close(fig) + + # Now for measured catalog + g1, g2 = meas_matched['cm_g'][:,0], meas_matched['cm_g'][:,1] + e1, e2 = ngmix.shape.g1g2_to_e1e2(g1, g2) + + # Easiest to use treecorr by making a new temporary catalog + meas_outfile = os.path.join(outdir,'treecorr_temp_file_meas.fits') + delete_file(meas_outfile) + fits = fitsio.FITS(meas_outfile,'rw', clobber=True) + data = [meas_matched['ra'], meas_matched['dec'], g1, g2, e1, e2] + names = ['ra', 'dec', 'g1', 'g2', 'e1', 'e2'] + fits.write(data, names=names) + fits.close() + + cat = treecorr.Catalog(meas_outfile, ra_col='ra', dec_col='dec', + ra_units='degrees', dec_units='degrees', + g1_col='e1', g2_col='e2') + gg = treecorr.GGCorrelation(min_sep=mins, max_sep=maxs, bin_size=bs, + sep_units='arcmin') + gg.process(cat) + fig = plot_gg_corr(gg, plt_type='Measured') + outfile = os.path.join(outdir, 'gg_corr_meas.png') + fig.savefig(outfile) + plt.close(fig) + + # Now for differences + g1t, g2t = truth_matched['cm_g'][:,0], truth_matched['cm_g'][:,1] + e1t, e2t = ngmix.shape.g1g2_to_e1e2(g1t, g2t) + g1m, g2m = meas_matched['cm_g'][:,0], meas_matched['cm_g'][:,1] + e1m, e2m = ngmix.shape.g1g2_to_e1e2(g1m, g2m) + g1d, g2d = g1m-g1t, g2m-g2t + e1d, e2d = e1m-e1t, e2m-e2t + + # Easiest to use treecorr by making a new temporary catalog + diff_outfile = os.path.join(outdir,'treecorr_temp_file_diff.fits') + delete_file(diff_outfile) + fits = fitsio.FITS(diff_outfile,'rw', clobber=True) + data = [truth_matched['ra'], truth_matched['dec'], g1d, g2d, e1d, e2d] + names = ['ra', 'dec', 'g1', 'g2', 'e1', 'e2'] + fits.write(data, names=names) + fits.close() + + cat = treecorr.Catalog(diff_outfile, ra_col='ra', dec_col='dec', + ra_units='degrees', dec_units='degrees', + g1_col='e1', g2_col='e2') + gg = treecorr.GGCorrelation(min_sep=mins, max_sep=maxs, bin_size=bs, + sep_units='arcmin') + gg.process(cat) + fig = plot_gg_corr(gg, plt_type='Measured-True') + outfile = os.path.join(outdir, 'gg_corr_diff.png') + fig.savefig(outfile) + plt.close(fig) + +def delete_file(filename): + try: + os.remove(filename) + except OSError as e: # this would be "except OSError, e:" before Python 2.6 + if e.errno != errno.ENOENT: # errno.ENOENT = no such file or directory + raise +def plot_gg_corr(gg, plt_type=None): + r = np.exp(gg.meanlogr) + xip = gg.xip + xim = gg.xim + sig = np.sqrt(gg.varxi) + + plt.plot(r, xip, color='blue') + # plt.plot(r, -xip, color='blue', ls=':') + plt.errorbar(r[xip>0], xip[xip>0], yerr=sig[xip>0], color='blue', lw=1, ls='') + # plt.errorbar(r[xip<0], -xip[xip<0], yerr=sig[xip<0], color='blue', lw=0.1, ls='') + lp = plt.errorbar(-r, xip, yerr=sig, color='blue') + + plt.plot(r, xim, color='green') + # plt.plot(r, -xim, color='green', ls=':') + plt.errorbar(r[xim>0], xim[xim>0], yerr=sig[xim>0], color='green', lw=1, ls='') + # plt.errorbar(r[xim<0], -xim[xim<0], yerr=sig[xim<0], color='green', lw=0.1, ls='') + lm = plt.errorbar(-r, xim, yerr=sig, color='green') + + # Reference line + plt.axhline(0, linestyle='--', c='k') + + plt.xscale('log') + # if plt_type == 'Measured' or plt_type=='Truth': + # plt.yscale('log', nonposy='clip') + # plt.yscale('log', nonposy='clip') + plt.xlabel(r'$\theta$ (arcmin)') + + plt.legend([lp, lm], [r'$\xi_+(\theta)$', r'$\xi_-(\theta)$']) + # plt.xlim( [1,100] ) + plt.ylabel(r'$\xi_{+,-}$') + plt.title(plt_type, fontsize=16) + + plt.gcf().set_size_inches(7,7) + + return plt.gcf() def make_star_plots(truth = None, meas_matched = None, truth_matched = None, filetag = ''): tag = 'cm_flux' @@ -316,8 +473,11 @@ def make_star_plots(truth = None, meas_matched = None, truth_matched = None, fil pass def main(argv): - # path = '../Data/sof_stars/' - run_name = 'grid_with_noise' + run_name = 'shear_test_ps' + # run_name = 'sof_stars' + # run_name = 'grid_with_noise' + # run_name = 'grid_test_ngmix' + # run_name = 'grid_test_shear_sof' path = os.path.join('/home/spencer/research/balrog/outputs/' + run_name) tilename = 'DES0347-5540' re_id = '0' diff --git a/plots/plot_YZ.py b/plots/plot_YZ.py index 28869ba..b917aac 100644 --- a/plots/plot_YZ.py +++ b/plots/plot_YZ.py @@ -5,6 +5,7 @@ ##this file also contains the followign functions readfiles(), match_func(), running_medians() ## dm_m_plot(), dm_T_plot(), dm_dT_plot() +from sys import argv import numpy as np import astropy.io.fits as pyfits import matplotlib.pyplot as plt @@ -587,13 +588,19 @@ def dm_dT_plot(t_gm, o_gm, t_sm, o_sm, up_perc=1, lo_perc=99, figname=None): plt.savefig(figname) return figname -def make_all(basepath=None, tile_list=None, realizations=None): +def make_all(basepath=None, tile_list=None, realizations=None, outdir=None): if basepath is None: basepath='/data/des71.a/data/kuropat/blank_test/y3v02/balrog_images/' if realizations is None: realizations=os.listdir(basepath) if tile_list is None: tile_list=os.listdir( os.path.join(basepath, realizations[0]) ) + if outdir is None: + outdir = os.getcwd() + + # Make output directory if it does not yet exist + if not os.path.exists(outdir): + os.makedirs(outdir) ##read in files tg, ts, oo=read_files(basepath, tile_list, realizations) @@ -604,13 +611,33 @@ def make_all(basepath=None, tile_list=None, realizations=None): ### make plots names=[] ##Diff_m vs True_m plots - names=np.append(names, dm_m_plot(truth_gm, obs_gm, truth_sm, obs_sm, figname='dm_m_YZ.png')) + fn1 = os.path.join(outdir, 'dm_m_YZ.png') + names=np.append(names, dm_m_plot(truth_gm, obs_gm, truth_sm, obs_sm, figname=fn1)) ##Diff_m vs Obs_T plots - names=np.append(names, dm_T_plot(truth_gm, obs_gm, truth_sm, obs_sm, figname='dm_T_YZ.png')) + fn2 = os.path.join(outdir, 'dm_T_YZ.png') + names=np.append(names, dm_T_plot(truth_gm, obs_gm, truth_sm, obs_sm, figname=fn2)) ##Diff_m vs diff T plots - names=np.append(names, dm_dT_plot(truth_gm, obs_gm, truth_sm, obs_sm, figname='dm_dT_gals_YZ.png')) + fn3 = os.path.join(outdir, 'dm_dT_gals_YZ.png') + names=np.append(names, dm_dT_plot(truth_gm, obs_gm, truth_sm, obs_sm, figname=fn3)) print 'genearted plots: ', names return names if __name__ == "__main__": - make_all() + + ## We can make this fancier, but for now this is simple enough. Could use argparse instead. + + # First argument is basepath + try: + print(argv[1]) + basepath = argv[1] + except IndexError: + basepath = None + + # Second argument is output directory + try: + print(argv[2]) + outdir = argv[2] + except IndexError: + outdir = None + + make_all(basepath=basepath, outdir=outdir)