Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Paper figs #223

Merged
merged 6 commits into from
Jun 17, 2021
Merged

Paper figs #223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5f7ef99
Paper figure updates
evevkovacs Apr 27, 2021
9110882
Updates to tests and config files required for validation paper
evevkovacs Jun 8, 2021
51e6985
Update descqa/ColorRedshiftTest.py
yymao Jun 9, 2021
c6f4720
Update descqa/ColorRedshiftTest.py
yymao Jun 9, 2021
a615102
Update descqa/BiasVersusRedshift.py
yymao Jun 9, 2021
ffbba36
Update descqa/BiasVersusRedshift.py
yymao Jun 9, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 84 additions & 11 deletions descqa/BiasVersusRedshift.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,40 @@

__all__ = ['BiasValidation']



def neglnlike(b, x, y, yerr):
return 0.5*np.sum((b**2*x-y)**2/yerr**2) # We ignore the covariance

def wtheta(x, b):
return b**2*x

class BiasValidation(CorrelationsAngularTwoPoint):
"""
Validation test of 2pt correlation function
"""

possible_observations = {'SRD':{
'filename_template': 'galaxy_bias/bias_SRD.txt',
'label': 'SRD ($i<25.3$)',
'colnames': ('z', 'bias'),
'skip':2,
},
'nicola_27':{
'filename_template': 'galaxy_bias/bias_nicola_mlim27.txt',
'label': 'Nicola et al. \n($i<27$)',
'colnames': ('z', 'bias'),
'skip':2,
},
'nicola_25.3':{
'filename_template': 'galaxy_bias/bias_nicola_mlim25.3.txt',
'label': 'Nicola et al. \n($i<25.3$)',
'colnames': ('z', 'bias'),
'skip':2,
},
}


def __init__(self, **kwargs): #pylint: disable=W0231
super().__init__(**kwargs)
self.data_label = kwargs['data_label']
Expand All @@ -39,19 +65,45 @@ def __init__(self, **kwargs): #pylint: disable=W0231
self.validation_data = np.loadtxt(validation_filepath, skiprows=2)
self.truncate_cat_name = kwargs.get('truncate_cat_name', False)
self.title_in_legend = kwargs.get('title_in_legend', True)
self.observations = kwargs.get('observations', [])

self.validation_data = self.get_validation_data(self.observations)

def get_validation_data(self, observations):

validation_data = {}
if observations:
for obs in observations:
print(obs)
data_args = self.possible_observations[obs]
fn = os.path.join(self.data_dir, data_args['filename_template'])
validation_data[obs] = dict(zip(data_args['colnames'], np.loadtxt(fn, skiprows=data_args['skip'], unpack=True)))
validation_data[obs]['label'] = data_args['label']
validation_data[obs]['colnames'] = data_args['colnames']

print(validation_data)
return validation_data


def plot_bias_results(self, corr_data, corr_theory, bias, z, catalog_name, output_dir):
fig, ax = plt.subplots(1, 2, gridspec_kw={'width_ratios': [5, 2]})
def plot_bias_results(self, corr_data, corr_theory, bias, z, catalog_name, output_dir, err=None):
fig, ax = plt.subplots(1, 2, gridspec_kw={'width_ratios': [5, 3]})
colors = plt.cm.plasma_r(np.linspace(0.1, 1, len(self.test_samples))) # pylint: disable=no-member

for sample_name, color in zip(self.test_samples, colors):
sample_corr = corr_data[sample_name]
sample_label = self.test_sample_labels.get(sample_name)
sample_th = corr_theory[sample_name]
ax[0].loglog(sample_corr[0], sample_th, c=color)
ax[0].errorbar(sample_corr[0], sample_corr[1], sample_corr[2], marker='o', ls='', c=color,
label=sample_label)

_, caps, bars = ax[0].errorbar(sample_corr[0], sample_corr[1], sample_corr[2], marker='o', ls='', c=color,
label=sample_label)
# add transparency for error bars
[bar.set_alpha(0.2) for bar in bars]
[cap.set_alpha(0.2) for cap in caps]
#add shaded band for fit range
ax[0].fill_between([self.fit_range[sample_name]['min_theta'], self.fit_range[sample_name]['max_theta']],
[self.fig_ylim[0], self.fig_ylim[0]], [self.fig_ylim[1], self.fig_ylim[1]],
alpha=0.07, color='grey')

if self.title_in_legend:
lgnd_title = catalog_name
title = self.data_label
Expand All @@ -63,10 +115,21 @@ def plot_bias_results(self, corr_data, corr_theory, bias, z, catalog_name, outpu
ax[0].set_ylim(*self.fig_ylim)
ax[0].set_ylabel(self.fig_ylabel, size=self.font_size)
ax[0].set_title(title, fontsize='medium')
ax[1].plot(z,bias)

ax[1].errorbar(z, bias, err, marker='o', ls='', label=catalog_name)
if not self.observations:
ax[1].plot(z, bias)
#add validation data
for v in self.validation_data.values():
colz = v['colnames'][0]
colb = v['colnames'][1]
zmask = (v[colz] < np.max(z)*1.25)
ax[1].plot(v[colz][zmask], v[colb][zmask], label=v['label'])

ax[1].set_title('Bias vs redshift', fontsize='medium')
ax[1].set_xlabel('$z$', size=self.font_size)
ax[1].set_ylabel('$b(z)$', size=self.font_size)
ax[1].legend(loc='upper right', framealpha=0.5, frameon=True, fontsize=self.legend_size)
plt.subplots_adjust(wspace=.05)
fig.tight_layout()
fig.savefig(os.path.join(output_dir, '{:s}.png'.format(self.test_name)), bbox_inches='tight')
Expand Down Expand Up @@ -102,6 +165,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
correlation_theory = dict()
best_fit_bias = []
z_mean = []
best_fit_err = []
for sample_name, sample_conditions in self.test_samples.items():
tmp_catalog_data = self.create_test_sample(
catalog_data, sample_conditions)
Expand Down Expand Up @@ -138,16 +202,25 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
args=(w_th[angles], xi[angles],
xi_sig[angles]), bounds=[(0.1, 10)])
best_bias = result['x']
#extract covariance matrix
#use curve_fit to get error on fit which has documented normalization for covariance matrix
cfit = op.curve_fit(wtheta, w_th[angles], xi[angles], p0=1.0, sigma=xi_sig[angles], bounds=(0.1, 10))
#best_bias_obj = result.hess_inv*np.identity(1)[0] #unknown relative normalization
best_bias_err = np.sqrt(cfit[1][0][0])
correlation_theory[sample_name] = best_bias**2*w_th
best_fit_bias.append(best_bias)
best_fit_err.append(best_bias_err)

z_mean = np.array(z_mean)
best_fit_bias = np.array(best_fit_bias)
best_fit_err = np.array(best_fit_err)
self.plot_bias_results(corr_data=correlation_data,
catalog_name=catalog_name,
corr_theory=correlation_theory,
bias=best_fit_bias,
z=z_mean,
output_dir=output_dir)
catalog_name=catalog_name,
corr_theory=correlation_theory,
bias=best_fit_bias,
z=z_mean,
output_dir=output_dir,
err=best_fit_err)

passed = np.all((best_fit_bias[1:]-best_fit_bias[:-1]) > 0)
score = np.count_nonzero((best_fit_bias[:-1]-best_fit_bias[1:])>0)*1.0/(len(best_fit_bias)-1.0)
Expand Down
2 changes: 1 addition & 1 deletion descqa/CheckColors.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -364,7 +364,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
print(xbins,ybins)
cntr1 = ax.contour(counts.transpose(), extent=[xmin,xmax,ymin,ymax],
colors='black',linestyles='solid',levels=self.levels)

ax.clabel(cntr1, inline=True, fmt='%1.1f', fontsize=10)
h1,_ = cntr1.legend_elements()

### CompareDensity block (Wasserstein metric)
Expand Down
94 changes: 88 additions & 6 deletions descqa/ColorRedshiftTest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
import os
import re
import numpy as np
from scipy.stats import binned_statistic as bs
import matplotlib.colors as clr
from .base import BaseValidationTest, TestResult
from .plotting import plt
Expand All @@ -19,15 +20,41 @@ class ColorRedshiftTest(BaseValidationTest):
"""
This test plots various color-redshfit diagnostics
"""

possible_observations = {
'des_fit': {'filename_template':'red_sequence/des/rykoff_et_al_1026',
'keys': (0),
'coefficients': (1, 2, 3, 4),
'skip': 7,
'label': 'DES fit',
'zmin': 0.2,
'zmax': 0.9,
'format': 'fit',
},
'des_y1':{'filename_template':'red_sequence/des/des_y1_redshift_ri_color.txt',
'skip': 1,
'usecols': (0,1),
'colnames': ('z', 'r-i'),
'label':'DES Y1',
'format': 'data',
'bins': [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.9],
},
}


def __init__(self, **kwargs):
super(ColorRedshiftTest, self).__init__()
# load test config options
self.kwargs = kwargs
self.truncate_cat_name = kwargs.get('truncate_cat_name', False)
self.title_in_legend = kwargs.get('title_in_legend', True)
self.font_size = kwargs.get('font_size', 16)
with open(os.path.join(self.data_dir, 'README.md')) as f:
self.validation_data = f.readline().strip()
self.text_size = kwargs.get('text_size', 12)
self.legend_size = kwargs.get('legend_size', 13)
self.observations = kwargs.get('observations', None)

#with open(os.path.join(self.data_dir, 'README.md')) as f:
# self.validation_data = f.readline().strip()
self.plot_list = kwargs.get("plot_list", [])
for plot_param in self.plot_list:
color = plot_param['color']
Expand All @@ -52,6 +79,34 @@ def __init__(self, **kwargs):
plot_param["redshift_block_limit"] = plot_param.get("redshift_block_limit", 1)
assert plot_param['redshift_block_limit'] in [1, 2, 3], "redshift_block_limit must be set to 1,2 or 3. It is set to: {}".format(plot_param['redshift_block_limit'])

#read in validation data
self.validation_data = self.get_validation_data(self.observations)

def get_validation_data(self, observations):
validation_data = {}
if observations:
for obs in observations:
data_args = self.possible_observations[obs]
fn = os.path.join(self.data_dir, data_args['filename_template'])
if 'keys' in data_args.keys():
keys = np.genfromtxt(fn, skip_header=data_args['skip'], usecols=data_args['keys'], dtype=str)
coefficients = np.genfromtxt(fn, skip_header=data_args['skip'], usecols=data_args['coefficients'])
validation_data[obs] = dict(zip(keys, coefficients))
else:
validation_data[obs] = dict(zip(data_args['colnames'],
np.loadtxt(fn, skiprows=data_args['skip'],
unpack=True, usecols=data_args['usecols'])))
validation_data[obs]['label'] = data_args['label']
validation_data[obs]['format'] = data_args['format']
if 'zmin' in data_args.keys():
validation_data[obs]['zmin'] = data_args['zmin']
if 'zmax' in data_args.keys():
validation_data[obs]['zmax'] = data_args['zmax']
if 'bins' in data_args.keys():
validation_data[obs]['bins'] = data_args['bins']

return validation_data

def post_process_plot(self, ax):
pass
# ax.text(0.05, 0.95, self.validation_data)
Expand Down Expand Up @@ -88,7 +143,7 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
slct, title = self._get_selection_and_title(catalog_instance, title, plot_param,
redshift_limit=plot_param['redshift_limit'],
redshift_block_limit=plot_param['redshift_block_limit'])
print(title)

fig, ax = plt.subplots()
# for ax_this in (ax, self.summary_ax):
if plot_param['redshift_limit'] is not None:
Expand All @@ -108,15 +163,41 @@ def run_on_single_catalog(self, catalog_instance, catalog_name, output_dir):
fig.colorbar(pc, ax=ax).set_label("Population Density")
mag1 = re.split('_', mag1_str)[1] #get filter
mag2 = re.split('_', mag2_str)[1] #get filter
# plot observations
for v in self.validation_data.values():
color=mag1 + '-' + mag2
if v['format'] == 'fit':
coeffs = v[color]
zmask = (redshift_bins >= v['zmin']) & (redshift_bins <= v['zmax'])
obs = np.zeros(len(redshift_bins[zmask]))
for n, coeff in enumerate(coeffs):
obs += coeff*redshift_bins[zmask]**(len(coeffs)-1-n)

ax.plot(redshift_bins[zmask], obs, color='r', label=v['label'])
elif v['format'] == 'data':
if color in v.keys():
zbins = np.asarray(v['bins'])
mean, _, num = bs(v['z'], v[color], bins=zbins)
std, _, num = bs(v['z'], v[color], bins=zbins, statistic='std')
fmask = np.isfinite(mean)
z_cen = 0.5*(zbins[1:]+zbins[:-1])
ax.errorbar(z_cen[fmask], mean[fmask], ls='', marker='o',
yerr=np.sqrt(std[fmask]), c='orange', label=v['label'])
counts = [np.sum(num==i+1) for i in range(len(zbins))]
print(z_cen[fmask], mean[fmask], std[fmask], counts)

legend = ax.legend(loc='lower right', fontsize=self.legend_size)
plt.setp(legend.get_texts(), color='w')

ax.set_ylabel('{} - {}'.format(mag1, mag2), size=self.font_size)
ax.set_xlabel('z', size=self.font_size)
ax.set_xlabel('Redshift $z$', size=self.font_size)
if self.title_in_legend:
title = '{}\n{}'.format(catalog_name, title)
else:
ax.set_title(catalog_name)
ax.text(0.05, 0.95, title, transform=ax.transAxes,
verticalalignment='top', color='white',
fontsize='small')
fontsize=self.text_size)
fig.savefig(os.path.join(output_dir, 'plot_{}.png'.format(plot_num)))
plt.close(fig)

Expand Down Expand Up @@ -221,7 +302,8 @@ def _get_selection_and_title(self, catalog_instance, title, plot_param,
redshift_limit=redshift_limit,
redshift_block_limit=redshift_block_limit)
slct = slct & (rs == plot_param["red_sequence_cut"])
title += "red sequence = {}, ".format(plot_param["red_sequence_cut"])
if plot_param["red_sequence_cut"]:
title += "red sequence galaxies, "
title_elem += 1
if title_elem % title_elem_per_line == 0:
title += "\n"
Expand Down
25 changes: 20 additions & 5 deletions descqa/CorrelationsTwoPoint.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,12 +66,17 @@ def __init__(self, **kwargs):
self.fig_subplot_groups = kwargs.get('fig_subplot_groups', [None])
self.fig_xlim = kwargs.get('fig_xlim', None)
self.tick_size = kwargs.get('tick_size', 12)

self.mask_large_errors = kwargs.get('mask_large_errors', False)

self.treecorr_config = {
'min_sep': kwargs['min_sep'],
'max_sep': kwargs['max_sep'],
'bin_size': kwargs['bin_size'],
}
if kwargs.get('var_method', None):
self.treecorr_config['var_method'] = kwargs['var_method']
self.npatch = kwargs.get('npatch', 1)

self.random_nside = kwargs.get('random_nside', 1024)
self.random_mult = kwargs.get('random_mult', 3)

Expand Down Expand Up @@ -224,7 +229,8 @@ def plot_data_comparison(self, corr_data, catalog_name, output_dir):
fig_xsize = 5 if self.fig_subplots_ncols==1 else 7 #widen figure for subplots
fig_ysize = 5 if self.fig_subplots_ncols==1 else 4 #narrow y-axis for subplots
fig, ax_all = plt.subplots(self.fig_subplots_nrows, self.fig_subplots_ncols, squeeze=False,
figsize=(min(2, self.fig_subplots_ncols)*fig_xsize, min(2, self.fig_subplots_nrows)*fig_ysize))
figsize=(min(2, self.fig_subplots_ncols)*fig_xsize,
min(2, self.fig_subplots_nrows)*fig_ysize))

for nx, (ax, this_group) in enumerate(zip(ax_all.flat, self.fig_subplot_groups)):
if this_group is None:
Expand Down Expand Up @@ -253,11 +259,15 @@ def plot_data_comparison(self, corr_data, catalog_name, output_dir):
self.validation_data[:, sample_data['data_err_col']])
y2 = (self.validation_data[:, sample_data['data_col']] -
self.validation_data[:, sample_data['data_err_col']])
if self.fig_xlim is not None:
if self.fig_ylim is not None:
y2[y2 <= 0] = self.fig_ylim[0]*0.9
ax.fill_between(self.validation_data[:, 0], y1, y2, lw=0, color=color, alpha=0.25)
if cat_data:
ax.errorbar(sample_corr[0], sample_corr[1], sample_corr[2],
if self.mask_large_errors and self.fig_ylim is not None:
mask = (sample_corr[1] - sample_corr[2]) > min(self.fig_ylim)
else:
mask = np.ones(len(sample_corr[1]), dtype=bool)
ax.errorbar(sample_corr[0][mask], sample_corr[1][mask], sample_corr[2][mask],
label=' '.join([catalog_name, sample_label]),
marker='o', ls='', c=color)

Expand Down Expand Up @@ -463,7 +473,9 @@ def generate_processed_randoms(self, catalog_data):
get_healpixel_footprint(catalog_data['ra'], catalog_data['dec'], self.random_nside),
self.random_nside,
)
rand_cat = treecorr.Catalog(ra=rand_ra, dec=rand_dec, ra_units='deg', dec_units='deg')
rand_cat = treecorr.Catalog(ra=rand_ra, dec=rand_dec, ra_units='deg', dec_units='deg',
npatch= self.npatch,
)
rr = treecorr.NNCorrelation(**self.treecorr_config)
rr.process(rand_cat)

Expand Down Expand Up @@ -497,6 +509,7 @@ def run_treecorr(self, catalog_data, treecorr_rand_cat, rr, output_file_name):
dec=catalog_data['dec'],
ra_units='deg',
dec_units='deg',
npatch= self.npatch,
)

dd = treecorr.NNCorrelation(**self.treecorr_config)
Expand Down Expand Up @@ -672,6 +685,7 @@ def run_treecorr_projected(self, catalog_data, rand_ra, rand_dec,
dec=catalog_data['dec'],
ra_units='deg',
dec_units='deg',
npatch=self.npatch,
r=redshift2dist(catalog_data['z'], cosmology),
)

Expand All @@ -683,6 +697,7 @@ def run_treecorr_projected(self, catalog_data, rand_ra, rand_dec,
dec=rand_dec,
ra_units='deg',
dec_units='deg',
npatch=self.npatch,
r=generate_uniform_random_dist(
rand_ra.size, *redshift2dist(np.array([z_min, z_max]), cosmology)),
)
Expand Down
Loading