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

new option per_r_par for covariance smoothing #1030

Merged
merged 2 commits into from
Jul 28, 2023
Merged

new option per_r_par for covariance smoothing #1030

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9970550
new option per_r_par for covariance smoothing
julienguy Jul 26, 2023
15c390a
dont compute index if not needed
julienguy Jul 27, 2023
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
10 changes: 9 additions & 1 deletion bin/picca_export.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,11 @@ def main(cmdargs):
action='store_true',
default=False,
help='Do not smooth the covariance matrix')
parser.add_argument(
'--smooth-per-r-par',
action='store_true',
default=False,
help='Consider different correlation coefficients per r_par')

parser.add_argument(
'--blind-corr-type',
Expand Down Expand Up @@ -149,12 +154,15 @@ def main(cmdargs):
delta_r_trans = (r_trans_max - 0.) / num_bins_r_trans
if not args.do_not_smooth_cov:
userprint("INFO: The covariance will be smoothed")
if args.smooth_per_r_par :
userprint("INFO: with different correlation coefficients per r_par")
covariance = smooth_cov(xi,
weights,
r_par,
r_trans,
delta_r_trans=delta_r_trans,
delta_r_par=delta_r_par)
delta_r_par=delta_r_par,
per_r_par=args.smooth_per_r_par)
else:
userprint("INFO: The covariance will not be smoothed")
covariance = compute_cov(xi, weights)
Expand Down
44 changes: 32 additions & 12 deletions py/picca/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,8 @@ def smooth_cov(xi,
r_trans,
delta_r_trans=4.0,
delta_r_par=4.0,
covariance=None):
covariance=None,
per_r_par=False):
"""Smoothes the covariance matrix

Args:
Expand All @@ -83,6 +84,8 @@ def smooth_cov(xi,
covariance: array of floats or None - defautl: None
Covariance matrix. If None, it will be computed using the
subsampling technique
per_r_par: if true, average the correlation coef per
(r_par,delta_r_par,delta_r_trans)

Returns:
The smooth covariance matrix. If data is not correct, then print a
Expand All @@ -107,30 +110,48 @@ def smooth_cov(xi,
sum_correlation = {}
counts_correlation = {}
for index in range(num_bins):
if per_r_par :
index_r_par = int(r_par[index]/delta_r_par)
userprint("\rsmoothing {}".format(index), end="")
for index2 in range(index + 1, num_bins):
index_delta_r_par = round(
abs(r_par[index2] - r_par[index]) / delta_r_par)
index_delta_r_trans = round(
abs(r_trans[index] - r_trans[index2]) / delta_r_trans)
if (index_delta_r_par, index_delta_r_trans) not in sum_correlation:
sum_correlation[(index_delta_r_par, index_delta_r_trans)] = 0.
counts_correlation[(index_delta_r_par, index_delta_r_trans)] = 0

sum_correlation[(index_delta_r_par,
index_delta_r_trans)] += correlation[index, index2]
counts_correlation[(index_delta_r_par, index_delta_r_trans)] += 1
if per_r_par :
if (index_r_par, index_delta_r_par, index_delta_r_trans) not in sum_correlation:
sum_correlation[(index_r_par, index_delta_r_par, index_delta_r_trans)] = 0.
counts_correlation[(index_r_par, index_delta_r_par, index_delta_r_trans)] = 0

sum_correlation[(index_r_par, index_delta_r_par,
index_delta_r_trans)] += correlation[index, index2]
counts_correlation[(index_r_par, index_delta_r_par, index_delta_r_trans)] += 1
else :
if (index_delta_r_par, index_delta_r_trans) not in sum_correlation:
sum_correlation[(index_delta_r_par, index_delta_r_trans)] = 0.
counts_correlation[(index_delta_r_par, index_delta_r_trans)] = 0

sum_correlation[(index_delta_r_par,
index_delta_r_trans)] += correlation[index, index2]
counts_correlation[(index_delta_r_par, index_delta_r_trans)] += 1

for index in range(num_bins):
correlation_smooth[index, index] = 1.
if per_r_par :
index_r_par = int(r_par[index]/delta_r_par)
for index2 in range(index + 1, num_bins):
index_delta_r_par = round(
abs(r_par[index2] - r_par[index]) / delta_r_par)
index_delta_r_trans = round(
abs(r_trans[index] - r_trans[index2]) / delta_r_trans)
correlation_smooth[index, index2] = (
sum_correlation[(index_delta_r_par, index_delta_r_trans)] /
counts_correlation[(index_delta_r_par, index_delta_r_trans)])
if per_r_par :
correlation_smooth[index, index2] = (
sum_correlation[(index_r_par, index_delta_r_par, index_delta_r_trans)] /
counts_correlation[(index_r_par , index_delta_r_par, index_delta_r_trans)])
else :
correlation_smooth[index, index2] = (
sum_correlation[(index_delta_r_par, index_delta_r_trans)] /
counts_correlation[(index_delta_r_par, index_delta_r_trans)])
correlation_smooth[index2, index] = correlation_smooth[index,
index2]

Expand Down Expand Up @@ -470,4 +491,3 @@ def unred(wave, ebv, R_V=3.1, LMC2=False, AVGLMC=False):
corr = 1. / (10.**(0.4 * curve))

return corr