From c751e8fa7fe5584d4c6b22936a7621c6329dc110 Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 1 Apr 2024 23:00:57 -0700 Subject: [PATCH 1/5] ana_cov_comp branch and GP diagonal smoother --- pspipe_utils/covariance.py | 31 +++++++++++++++++++++++++++++++ setup.py | 1 + 2 files changed, 32 insertions(+) diff --git a/pspipe_utils/covariance.py b/pspipe_utils/covariance.py index 7a1e19f..3c918ea 100644 --- a/pspipe_utils/covariance.py +++ b/pspipe_utils/covariance.py @@ -7,6 +7,9 @@ from scipy.optimize import curve_fit import pylab as plt +from sklearn.gaussian_process import GaussianProcessRegressor +from sklearn.gaussian_process.kernels import RBF, WhiteKernel + from pixell import utils from itertools import combinations_with_replacement as cwr @@ -339,7 +342,35 @@ def skew(cov, dir=1): return S, corrected_cov else: return corrected_cov + + +def smooth_gp_diag(lb, res, ell_cut, ell_cut_transition=50, length_scale=100.0, + length_scale_bounds=(10, 1e4), noise_level=0.1, noise_level_bounds=(1e-6, 1e1)): + kernel = 1.0 * RBF(length_scale=length_scale, + length_scale_bounds=length_scale_bounds) + WhiteKernel( + noise_level=noise_level, noise_level_bounds=noise_level_bounds + ) + # fit the first GP on the bins above the ell_cut + gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0, normalize_y=True) # sample mean + i_cut = np.argmax(lb > ell_cut) + X_train = lb[i_cut:,np.newaxis] + y_train = np.diag(res)[i_cut:] + gpr.fit(X_train, y_train) + y_mean_high = gpr.predict(lb[:,np.newaxis], return_std=False) + + # fit the second GP on the residual, for the bins below the ell_cut + gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0, normalize_y=False) # mean 0 + i_cut = np.argmax(lb > ell_cut) + X_train = np.array([lb[:i_cut]]).T + y_train = (np.diag(res) - y_mean_high)[:i_cut] + gpr.fit(X_train, y_train) + y_mean_low = gpr.predict(lb[:,np.newaxis], return_std=False) + + # add the second GP but smoothly transition it to zero at the ell cut + step = 1 - 1 / (1 + np.exp(-2 * (lb - ell_cut) / (ell_cut_transition/2))) + return y_mean_high + step * y_mean_low + def correct_analytical_cov_keep_res_diag(an_full_cov, mc_full_cov, return_diag=False): sqrt_an_full_cov = utils.eigpow(an_full_cov, 0.5) diff --git a/setup.py b/setup.py index 660415f..9b16b7a 100644 --- a/setup.py +++ b/setup.py @@ -17,6 +17,7 @@ install_requires=[ "pspy>=1.5.3", "mflike>=0.8.2", + "scikit-learn>=1" ], package_data={"": ["data/**"]}, ) From 967798c6b1ca02e3db9b93224013a82b22eddc46 Mon Sep 17 00:00:00 2001 From: xzackli Date: Mon, 1 Apr 2024 23:04:57 -0700 Subject: [PATCH 2/5] remove smooth step --- pspipe_utils/covariance.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pspipe_utils/covariance.py b/pspipe_utils/covariance.py index 3c918ea..a2aab73 100644 --- a/pspipe_utils/covariance.py +++ b/pspipe_utils/covariance.py @@ -344,7 +344,7 @@ def skew(cov, dir=1): return corrected_cov -def smooth_gp_diag(lb, res, ell_cut, ell_cut_transition=50, length_scale=100.0, +def smooth_gp_diag(lb, res, ell_cut, length_scale=100.0, length_scale_bounds=(10, 1e4), noise_level=0.1, noise_level_bounds=(1e-6, 1e1)): kernel = 1.0 * RBF(length_scale=length_scale, @@ -367,9 +367,9 @@ def smooth_gp_diag(lb, res, ell_cut, ell_cut_transition=50, length_scale=100.0, gpr.fit(X_train, y_train) y_mean_low = gpr.predict(lb[:,np.newaxis], return_std=False) - # add the second GP but smoothly transition it to zero at the ell cut - step = 1 - 1 / (1 + np.exp(-2 * (lb - ell_cut) / (ell_cut_transition/2))) - return y_mean_high + step * y_mean_low + # add together the two GPs but only include the second GP below the ell_cut + y_mean_high[:i_cut] += y_mean_low[:i_cut] + return y_mean_high def correct_analytical_cov_keep_res_diag(an_full_cov, mc_full_cov, return_diag=False): From 3614422594354a39a0d9a0dfe902b40084708b2d Mon Sep 17 00:00:00 2001 From: Xavier Garrido Date: Tue, 2 Apr 2024 08:38:27 +0200 Subject: [PATCH 3/5] Fix missing comma --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 41e0e2c..9698c7d 100644 --- a/setup.py +++ b/setup.py @@ -16,7 +16,7 @@ python_requires=">=3.9", install_requires=[ "pspy>=1.5.3", - "scikit-learn>=1" + "scikit-learn>=1", "mflike>=0.9.5", ], package_data={"": ["data/**"]}, From b3fa51f6d4dad3abc85ff92a5f491812466fb394 Mon Sep 17 00:00:00 2001 From: xzackli Date: Tue, 2 Apr 2024 08:20:45 -0700 Subject: [PATCH 4/5] fit exponential at low ell instead --- pspipe_utils/covariance.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/pspipe_utils/covariance.py b/pspipe_utils/covariance.py index 3f605ca..f298d94 100644 --- a/pspipe_utils/covariance.py +++ b/pspipe_utils/covariance.py @@ -347,31 +347,32 @@ def skew(cov, dir=1): return corrected_cov -def smooth_gp_diag(lb, res, ell_cut, length_scale=100.0, - length_scale_bounds=(10, 1e4), noise_level=0.1, noise_level_bounds=(1e-6, 1e1)): +def smooth_gp_diag(lb, arr_diag, ell_cut, length_scale=500.0, + length_scale_bounds=(100, 1e4), noise_level=0.01, + noise_level_bounds=(1e-6, 1e1), low_ell_scale=100, n_restarts_optimizer=20): kernel = 1.0 * RBF(length_scale=length_scale, length_scale_bounds=length_scale_bounds) + WhiteKernel( noise_level=noise_level, noise_level_bounds=noise_level_bounds ) # fit the first GP on the bins above the ell_cut - gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0, normalize_y=True) # sample mean + gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0, normalize_y=True, + n_restarts_optimizer=n_restarts_optimizer) i_cut = np.argmax(lb > ell_cut) X_train = lb[i_cut:,np.newaxis] - y_train = np.diag(res)[i_cut:] + y_train = arr_diag[i_cut:] gpr.fit(X_train, y_train) y_mean_high = gpr.predict(lb[:,np.newaxis], return_std=False) - # fit the second GP on the residual, for the bins below the ell_cut - gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0, normalize_y=False) # mean 0 + # fit an exponential at the low end i_cut = np.argmax(lb > ell_cut) - X_train = np.array([lb[:i_cut]]).T - y_train = (np.diag(res) - y_mean_high)[:i_cut] - gpr.fit(X_train, y_train) - y_mean_low = gpr.predict(lb[:,np.newaxis], return_std=False) - - # add together the two GPs but only include the second GP below the ell_cut - y_mean_high[:i_cut] += y_mean_low[:i_cut] + X_train = lb[:i_cut] + y_train = (arr_diag - y_mean_high)[:i_cut] + pos_el = y_train > 0 + X_train, y_train = X_train[pos_el], y_train[pos_el] + z = np.polyfit(X_train, np.log(y_train), 1) + f = np.poly1d(z) + y_mean_high[:i_cut] += np.exp(f(lb[:i_cut])) return y_mean_high From 347032982cb08cdf65e82e26f76f76302adf48b7 Mon Sep 17 00:00:00 2001 From: xzackli Date: Tue, 2 Apr 2024 08:42:52 -0700 Subject: [PATCH 5/5] remove extraneous parameter and smoooth diag --- pspipe_utils/covariance.py | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/pspipe_utils/covariance.py b/pspipe_utils/covariance.py index f298d94..cc42207 100644 --- a/pspipe_utils/covariance.py +++ b/pspipe_utils/covariance.py @@ -349,7 +349,7 @@ def skew(cov, dir=1): def smooth_gp_diag(lb, arr_diag, ell_cut, length_scale=500.0, length_scale_bounds=(100, 1e4), noise_level=0.01, - noise_level_bounds=(1e-6, 1e1), low_ell_scale=100, n_restarts_optimizer=20): + noise_level_bounds=(1e-6, 1e1), n_restarts_optimizer=20): kernel = 1.0 * RBF(length_scale=length_scale, length_scale_bounds=length_scale_bounds) + WhiteKernel( @@ -367,9 +367,7 @@ def smooth_gp_diag(lb, arr_diag, ell_cut, length_scale=500.0, # fit an exponential at the low end i_cut = np.argmax(lb > ell_cut) X_train = lb[:i_cut] - y_train = (arr_diag - y_mean_high)[:i_cut] - pos_el = y_train > 0 - X_train, y_train = X_train[pos_el], y_train[pos_el] + y_train = np.abs(arr_diag - y_mean_high)[:i_cut] z = np.polyfit(X_train, np.log(y_train), 1) f = np.poly1d(z) y_mean_high[:i_cut] += np.exp(f(lb[:i_cut])) @@ -401,6 +399,25 @@ def correct_analytical_cov_keep_res_diag(an_full_cov, mc_full_cov, return_diag=F else: return corrected_cov +def correct_analytical_cov_keep_res_diag_gp(an_full_cov, mc_full_cov, lb, ell_cut, return_diag=False): + d_an, O_an = np.linalg.eigh(an_full_cov) + sqrt_an_full_cov = O_an @ np.diag(d_an**.5) + inv_sqrt_an_full_cov = np.diag(d_an**-.5) @ O_an.T + res = inv_sqrt_an_full_cov @ mc_full_cov @ inv_sqrt_an_full_cov.T # res should be close to the identity if an_full_cov is good + res_diag = np.diag(res) + + n_spec = len(res_diag) // len(lb) + diags = np.array_split(res_diag, n_spec) + smoothed_diags = [smooth_gp_diag(lb, r, ell_cut) for r in diags] + smooth_res = np.hstack(smoothed_diags) + + corrected_cov = sqrt_an_full_cov @ np.diag(smooth_res) @ sqrt_an_full_cov.T + + if return_diag: + return corrected_cov, res_diag + else: + return corrected_cov + def canonize_connected_2pt(leg1, leg2, all_legs): """A connected 2-point term has two legs but is invariant to their