From d499e139f08e508c232a0a88e9bccf1d9a029b0c Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Wed, 3 Jul 2024 05:48:19 -0700 Subject: [PATCH 01/11] new covariance implementation vectorized --- py/picca/pk1d/postproc_pk1d.py | 176 +++++++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index 94750b683..5bf38463f 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -862,6 +862,7 @@ def compute_average_pk_redshift( ) + def compute_cov( p1d_table, mean_p1d_table, @@ -924,6 +925,180 @@ def compute_cov( n_array (array-like): Array of the number of power spectra used to compute each covariance coefficient. + covariance_array (array-like): + Array of covariance coefficients. + """ + zbin_array = np.full(nbins_k * nbins_k, zbin_centers[izbin]) + index_zbin_array = np.full(nbins_k * nbins_k, izbin, dtype=int) + kbin_centers = (kbin_edges[1:] + kbin_edges[:-1]) / 2 + k1_array, k2_array = np.meshgrid(kbin_centers, kbin_centers, indexing="ij") + k1_array = np.ravel(k1_array) + k2_array = np.ravel(k2_array) + n_array = np.zeros(nbins_k * nbins_k) + covariance_array = np.zeros(nbins_k * nbins_k) + if n_chunks[izbin] == 0: # Fill rows with NaNs + covariance_array[:] = np.nan + return ( + zbin_array, + index_zbin_array, + n_array, + covariance_array, + k1_array, + k2_array, + ) + n_array = np.ravel( + np.outer( + mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], + mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], + ) + ) + + # Computing p1d groups, weights used in the covariance and p1d weights + p1d_groups = np.zeros((len(sub_forest_ids), nbins_k)) + covariance_weights = np.zeros((len(sub_forest_ids), nbins_k)) + p1d_weights = np.zeros((len(sub_forest_ids), nbins_k)) + + for i, sub_forest_id in enumerate(sub_forest_ids): + select_id = select_z & (p1d_table["sub_forest_id"] == sub_forest_id) + selected_pk = p1d_table["Pk"][select_id] + selected_ikbin = k_index[select_id] + if weight_method == "fit_snr": + selected_snr = p1d_table["forest_snr"][select_id] + snrfit_z = snrfit_table[izbin * nbins_k : (izbin + 1) * nbins_k, :] + selected_variance_estimated = fitfunc_variance_pk1d( + selected_snr, snrfit_z[selected_ikbin, 2], snrfit_z[selected_ikbin, 3] + ) + + for ikbin, _ in enumerate(kbin_centers): + mask_ikbin = selected_ikbin == ikbin + number_in_bins = len(mask_ikbin[mask_ikbin]) + if number_in_bins != 0: + if weight_method == "fit_snr": + weight = 1 / selected_variance_estimated[mask_ikbin][0] + p1d_weights[i][ikbin] = weight + # covariance_weights[i][ikbin] = np.sqrt(weight / number_in_bins) + covariance_weights[i][ikbin] = np.sqrt(weight) + else: + p1d_weights[i][ikbin] = 1 + covariance_weights[i][ikbin] = np.sqrt(1 / number_in_bins) + + p1d_groups[i][ikbin] = np.nansum( + selected_pk[mask_ikbin] * covariance_weights[i][ikbin] + ) + + # Calculation of covariance terms + mean_pk = mean_p1d_table["meanPk"][(nbins_k * izbin) : (nbins_k * (izbin + 1))] + mean_pk_product = np.outer(mean_pk, mean_pk) + + sum_p1d_weights = np.nansum(p1d_weights, axis=0) + weights_sum_product = np.outer(sum_p1d_weights, sum_p1d_weights) + + p1d_groups_product_sum = np.zeros((nbins_k, nbins_k)) + covariance_weights_product_sum = np.zeros((nbins_k, nbins_k)) + weights_product_sum = np.zeros((nbins_k, nbins_k)) + + for i, _ in enumerate(sub_forest_ids): + p1d_groups_product_sum = np.nansum( + [p1d_groups_product_sum, np.outer(p1d_groups[i], p1d_groups[i])], axis=0 + ) + covariance_weights_product_sum = np.nansum( + [ + covariance_weights_product_sum, + np.outer(covariance_weights[i], covariance_weights[i]), + ], + axis=0, + ) + weights_product_sum = np.nansum( + [weights_product_sum, np.outer(p1d_weights[i], p1d_weights[i])], axis=0 + ) + + del p1d_groups, covariance_weights + + covariance_matrix = (weights_product_sum / weights_sum_product) * ( + (p1d_groups_product_sum / covariance_weights_product_sum) - mean_pk_product + ) + + # For fit_snr method, due to the SNR fitting scheme used for weighting, + # the diagonal of the weigthed sample covariance matrix is not equal + # to the error in mean P1D. This is tested on Ohio mocks and data. + # We choose to renormalize the whole covariance matrix. + if weight_method == "fit_snr": + std_pk = mean_p1d_table["errorPk"][(nbins_k * izbin) : (nbins_k * (izbin + 1))] + covariance_diag = np.diag(covariance_matrix) + covariance_matrix = ( + covariance_matrix + * np.outer(std_pk, std_pk) + / np.sqrt(np.outer(covariance_diag, covariance_diag)) + ) + + covariance_array = np.ravel(covariance_matrix) + + return zbin_array, index_zbin_array, n_array, covariance_array, k1_array, k2_array + + + +def compute_cov_not_vectorized( + p1d_table, + mean_p1d_table, + zbin_centers, + n_chunks, + kbin_edges, + k_index, + nbins_k, + weight_method, + snrfit_table, + izbin, + select_z, + sub_forest_ids, +): + """Compute the covariance of a set of 1D power spectra. + + Arguments + --------- + p1d_table (array-like): + Table of 1D power spectra, with columns 'Pk' and 'sub_forest_id'. + + mean_p1d_table (array-like): + Table of mean 1D power spectra, with column 'meanPk'. + + zbin_centers (array-like): + Array of bin centers for redshift. + + n_chunks (array-like): + Array of the number of chunks in each redshift bin. + + k_index (array-like): + Array of indices for k-values, with -1 indicating values outside of the k bins. + + nbins_k (int): + Number of k bins. + + weight_method: str, + Method to weight the data. + + snrfit_table: numpy ndarray, + Table containing SNR fit infos + + izbin (int): + Current redshift bin being considered. + + select_z (array-like): + Boolean array for selecting data points based on redshift. + + sub_forest_ids (array-like): + Array of chunk ids. + + Return + ------ + zbin_array (array-like): + Array of redshift bin centers for each covariance coefficient. + + index_zbin_array (array-like): + Array of redshift bin indices for each covariance coefficient. + + n_array (array-like): + Array of the number of power spectra used to compute each covariance coefficient. + covariance_array (array-like): Array of covariance coefficients. """ @@ -1079,6 +1254,7 @@ def compute_cov( return zbin_array, index_zbin_array, n_array, covariance_array, k1_array, k2_array + def run_postproc_pk1d( data_dir, output_file, From b364ef380fc729455b4c017c2a209d9fe6266377 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Thu, 4 Jul 2024 06:28:34 -0700 Subject: [PATCH 02/11] new version of covariance computation, optimized for fast and low memory calculation --- py/picca/pk1d/postproc_pk1d.py | 399 ++++++++++++++++++--------------- 1 file changed, 216 insertions(+), 183 deletions(-) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index 5bf38463f..d66668baa 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -12,22 +12,20 @@ See the respective documentation for details """ -import os import glob +import os +from functools import partial from multiprocessing import Pool -from functools import partial -import numpy as np -from scipy.stats import binned_statistic -from scipy.optimize import curve_fit import fitsio -from astropy.table import Table, vstack +import numpy as np from astropy.stats import bootstrap - -from picca.constants import SPEED_LIGHT -from picca.constants import ABSORBER_IGM -from picca.utils import userprint +from astropy.table import Table, vstack +from picca.constants import ABSORBER_IGM, SPEED_LIGHT from picca.pk1d.utils import MEANPK_FITRANGE_SNR, fitfunc_variance_pk1d +from picca.utils import userprint +from scipy.optimize import curve_fit +from scipy.stats import binned_statistic def read_pk1d(filename, kbin_edges, snrcut=None, zbins_snrcut=None): @@ -91,9 +89,9 @@ def read_pk1d(filename, kbin_edges, snrcut=None, zbins_snrcut=None): chunk_table["forest_snr"] = float(chunk_header["MEANSNR"]) chunk_table["forest_id"] = int(chunk_header["LOS_ID"]) if "CHUNK_ID" in chunk_header: - chunk_table[ - "sub_forest_id" - ] = f"{chunk_header['LOS_ID']}_{chunk_header['CHUNK_ID']}" + chunk_table["sub_forest_id"] = ( + f"{chunk_header['LOS_ID']}_{chunk_header['CHUNK_ID']}" + ) if snrcut is not None: if len(snrcut) > 1: @@ -111,7 +109,8 @@ def read_pk1d(filename, kbin_edges, snrcut=None, zbins_snrcut=None): # Empirically remove very noisy chunks (wk,) = np.where(chunk_table["k"] < kbin_edges[-1]) if ( - np.abs(chunk_table["Pk_noise"][wk]) > 1000000 * np.abs(chunk_table["Pk_raw"][wk]) + np.abs(chunk_table["Pk_noise"][wk]) + > 1000000 * np.abs(chunk_table["Pk_raw"][wk]) ).any(): userprint( f"file {filename} hdu {i+1} has very high noise power: discarded" @@ -224,7 +223,7 @@ def compute_mean_pk1d( if "sub_forest_id" in p1d_table_cols: p1d_table_cols.remove("sub_forest_id") - p1d_table_cols = [ col for col in p1d_table_cols if "Delta_" not in col ] + p1d_table_cols = [col for col in p1d_table_cols if "Delta_" not in col] # Convert data into velocity units if velunits: @@ -345,105 +344,119 @@ def compute_mean_pk1d( nomedians, ) - if compute_covariance: - userprint("Computing covariance matrix") - params_pool = [] - # Main loop 1) z bins + if compute_covariance or compute_bootstrap: + userprint("Computation of p1d groups for covariance matrix calculation") + + p1d_groups = [] for izbin in range(nbins_z): - select_z = (p1d_table["forest_z"] < zbin_edges[izbin + 1]) & ( - p1d_table["forest_z"] > zbin_edges[izbin] + + zbin_array = np.full(nbins_k * nbins_k, zbin_centers[izbin]) + index_zbin_array = np.full(nbins_k * nbins_k, izbin, dtype=int) + kbin_centers = (kbin_edges[1:] + kbin_edges[:-1]) / 2 + k1_array, k2_array = np.meshgrid(kbin_centers, kbin_centers, indexing="ij") + k1_array = np.ravel(k1_array) + k2_array = np.ravel(k2_array) + n_array = np.ravel( + np.outer( + mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], + mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], + ) + ) + i_min = izbin * nbins_k * nbins_k + i_max = (izbin + 1) * nbins_k * nbins_k + cov_table["zbin"][i_min:i_max] = zbin_array + cov_table["index_zbin"][i_min:i_max] = index_zbin_array + cov_table["k1"][i_min:i_max] = k1_array + cov_table["k2"][i_min:i_max] = k2_array + cov_table["N"][i_min:i_max] = n_array + + if n_chunks[izbin] == 0: + p1d_weights_z, covariance_weights_z, p1d_groups_z = [], [], [] + else: + p1d_weights_z, covariance_weights_z, p1d_groups_z = compute_p1d_groups( + weight_method, + nbins_k, + zbin_edges, + izbin, + p1d_table, + k_index, + snrfit_table, + number_worker, + ) + p1d_groups.append( + [ + izbin, + p1d_weights_z, + covariance_weights_z, + p1d_groups_z, + ] ) - sub_forest_ids = np.unique(p1d_table["sub_forest_id"][select_z]) - params_pool.append([izbin, select_z, sub_forest_ids]) + + if compute_covariance: + userprint("Computation of covariance matrix") func = partial( compute_cov, - p1d_table, + weight_method, mean_p1d_table, - zbin_centers, - n_chunks, - kbin_edges, - k_index, nbins_k, - weight_method, - snrfit_table, ) if number_worker == 1: - output_cov = [func(*p) for p in params_pool] + output_cov = [func(*p1d_group) for p1d_group in p1d_groups] else: with Pool(number_worker) as pool: - output_cov = pool.starmap(func, params_pool) - # Main loop 1) z bins + output_cov = pool.starmap(func, p1d_groups) + for izbin in range(nbins_z): - ( - zbin_array, - index_zbin_array, - n_array, - covariance_array, - k1_array, - k2_array, - ) = (*output_cov[izbin],) + covariance_array = output_cov[izbin] i_min = izbin * nbins_k * nbins_k i_max = (izbin + 1) * nbins_k * nbins_k - cov_table["zbin"][i_min:i_max] = zbin_array - cov_table["index_zbin"][i_min:i_max] = index_zbin_array - cov_table["N"][i_min:i_max] = n_array cov_table["covariance"][i_min:i_max] = covariance_array - cov_table["k1"][i_min:i_max] = k1_array - cov_table["k2"][i_min:i_max] = k2_array if compute_bootstrap: userprint("Computing covariance matrix with bootstrap method") - - params_pool = [] - # Main loop 1) z bins + p1d_groups_bootstrap = [] for izbin in range(nbins_z): - select_z = (p1d_table["forest_z"] < zbin_edges[izbin + 1]) & ( - p1d_table["forest_z"] > zbin_edges[izbin] - ) - - sub_forest_ids = np.unique(p1d_table["sub_forest_id"][select_z]) - if sub_forest_ids.size > 0: + number_sub_forests = len(p1d_groups[izbin][1]) + if number_sub_forests > 0: bootid = np.array( - bootstrap(np.arange(sub_forest_ids.size), number_bootstrap) + bootstrap(np.arange(number_sub_forests), number_bootstrap) ).astype(int) else: bootid = np.full(number_bootstrap, None) - # Main loop 2) number of bootstrap samples for iboot in range(number_bootstrap): - params_pool.append([izbin, select_z, sub_forest_ids[bootid[iboot]]]) + if bootid[iboot] is None: + p1d_weights_z, covariance_weights_z, p1d_groups_z = [], [], [] + else: + p1d_weights_z = p1d_groups[izbin][1][bootid[iboot]] + covariance_weights_z = p1d_groups[izbin][2][bootid[iboot]] + p1d_groups_z = p1d_groups[izbin][3][bootid[iboot]] + p1d_groups_bootstrap.append( + [ + izbin, + p1d_weights_z, + covariance_weights_z, + p1d_groups_z, + ] + ) func = partial( compute_cov, - p1d_table, + weight_method, mean_p1d_table, - zbin_centers, - n_chunks, - kbin_edges, - k_index, nbins_k, - weight_method, - snrfit_table, ) if number_worker == 1: - output_cov = [func(*p) for p in params_pool] + output_cov = [func(*p) for p in p1d_groups_bootstrap] else: with Pool(number_worker) as pool: - output_cov = pool.starmap(func, params_pool) - # Main loop 1) z bins + output_cov = pool.starmap(func, p1d_groups_bootstrap) + for izbin in range(nbins_z): boot_cov = [] - # Main loop 2) number of bootstrap samples for iboot in range(number_bootstrap): - ( - zbin_array, - index_zbin_array, - n_array, - covariance_array, - k1_array, - k2_array, - ) = (*output_cov[izbin * number_bootstrap + iboot],) + covariance_array = (output_cov[izbin * number_bootstrap + iboot],) boot_cov.append(covariance_array) i_min = izbin * nbins_k * nbins_k @@ -829,8 +842,7 @@ def compute_average_pk_redshift( mean = np.average(p1d_table[col][select], weights=redshift_weights) # simple analytic expression: error = np.std(p1d_table[col][select]) * ( - np.sqrt(np.sum(redshift_weights**2)) - / np.sum(redshift_weights) + np.sqrt(np.sum(redshift_weights**2)) / np.sum(redshift_weights) ) else: mean = np.mean(p1d_table[col][select]) @@ -862,20 +874,127 @@ def compute_average_pk_redshift( ) +def compute_p1d_groups( + weight_method, + nbins_k, + zbin_edges, + izbin, + p1d_table, + k_index, + snrfit_table, + number_worker, +): + """Compute the P1D groups before covariance matrix calculation. + Put all the P1D in the same k grid. + + Arguments + --------- + weight_method: str, + Method to weight the data. + + nbins_k (int): + Number of k bins. + + zbin_edges (array-like): + All redshift bins. + + izbin (int): + Current redshift bin being considered. + + p1d_table (array-like): + Table of 1D power spectra, with columns 'Pk' and 'sub_forest_id'. + + k_index (array-like): + Array of indices for k-values, with -1 indicating values outside of the k bins. + + snrfit_table: numpy ndarray, + Table containing SNR fit infos + + number_worker: int + Number of workers for the parallelization + + Return + ------ + p1d_weights_z (array-like): + + covariance_weights_z (array-like): + + p1d_groups_z (array-like): + + """ + + select_z = (p1d_table["forest_z"] < zbin_edges[izbin + 1]) & ( + p1d_table["forest_z"] > zbin_edges[izbin] + ) + p1d_sub_table = Table( + [ + p1d_table["sub_forest_id"][select_z], + k_index[select_z], + p1d_table["Pk"][select_z], + ], + names=("sub_forest_id", "k_index", "pk"), + ) + if weight_method == "fit_snr": + snrfit_z = snrfit_table[izbin * nbins_k : (izbin + 1) * nbins_k, :] + p1d_sub_table["weight"] = 1 / fitfunc_variance_pk1d( + p1d_table["forest_snr"][select_z], + snrfit_z[k_index[select_z], 2], + snrfit_z[k_index[select_z], 3], + ) + else: + p1d_sub_table["weight"] = 1 + + p1d_sub_table = p1d_sub_table[p1d_sub_table["k_index"] >= 0] + + grouped_table = p1d_sub_table.group_by("sub_forest_id") + p1d_los_table = [ + group[["k_index", "pk", "weight"]] for group in grouped_table.groups + ] + + del grouped_table + + func = partial( + compute_groups, + nbins_k, + ) + if number_worker == 1: + output_cov = [func(*p1d_los) for p1d_los in p1d_los_table] + else: + with Pool(number_worker) as pool: + output_cov = np.array(pool.map(func, p1d_los_table)) + + del p1d_los_table + + return output_cov[:, 0, :], output_cov[:, 1, :], output_cov[:, 2, :] + + +def compute_groups(nbins_k, p1d_los): + + p1d_weights = np.zeros(nbins_k) + covariance_weights = np.zeros(nbins_k) + p1d_groups = np.zeros(nbins_k) + + for ikbin in range(nbins_k): + mask_ikbin = p1d_los["k_index"] == ikbin + number_in_bins = len(mask_ikbin[mask_ikbin]) + if number_in_bins != 0: + weight = p1d_los["weight"][mask_ikbin][0] + p1d_weights[ikbin] = weight + covariance_weights[ikbin] = np.sqrt(weight / number_in_bins) + p1d_groups[ikbin] = np.nansum( + p1d_los["pk"][mask_ikbin] * covariance_weights[ikbin] + ) + return p1d_weights, covariance_weights, p1d_groups + def compute_cov( - p1d_table, + weight_method, mean_p1d_table, - zbin_centers, - n_chunks, - kbin_edges, - k_index, nbins_k, - weight_method, - snrfit_table, izbin, - select_z, - sub_forest_ids, + p1d_weights, + covariance_weights, + p1d_groups, ): """Compute the covariance of a set of 1D power spectra. @@ -928,65 +1047,10 @@ def compute_cov( covariance_array (array-like): Array of covariance coefficients. """ - zbin_array = np.full(nbins_k * nbins_k, zbin_centers[izbin]) - index_zbin_array = np.full(nbins_k * nbins_k, izbin, dtype=int) - kbin_centers = (kbin_edges[1:] + kbin_edges[:-1]) / 2 - k1_array, k2_array = np.meshgrid(kbin_centers, kbin_centers, indexing="ij") - k1_array = np.ravel(k1_array) - k2_array = np.ravel(k2_array) - n_array = np.zeros(nbins_k * nbins_k) - covariance_array = np.zeros(nbins_k * nbins_k) - if n_chunks[izbin] == 0: # Fill rows with NaNs - covariance_array[:] = np.nan - return ( - zbin_array, - index_zbin_array, - n_array, - covariance_array, - k1_array, - k2_array, - ) - n_array = np.ravel( - np.outer( - mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], - mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], - ) - ) - - # Computing p1d groups, weights used in the covariance and p1d weights - p1d_groups = np.zeros((len(sub_forest_ids), nbins_k)) - covariance_weights = np.zeros((len(sub_forest_ids), nbins_k)) - p1d_weights = np.zeros((len(sub_forest_ids), nbins_k)) - - for i, sub_forest_id in enumerate(sub_forest_ids): - select_id = select_z & (p1d_table["sub_forest_id"] == sub_forest_id) - selected_pk = p1d_table["Pk"][select_id] - selected_ikbin = k_index[select_id] - if weight_method == "fit_snr": - selected_snr = p1d_table["forest_snr"][select_id] - snrfit_z = snrfit_table[izbin * nbins_k : (izbin + 1) * nbins_k, :] - selected_variance_estimated = fitfunc_variance_pk1d( - selected_snr, snrfit_z[selected_ikbin, 2], snrfit_z[selected_ikbin, 3] - ) - - for ikbin, _ in enumerate(kbin_centers): - mask_ikbin = selected_ikbin == ikbin - number_in_bins = len(mask_ikbin[mask_ikbin]) - if number_in_bins != 0: - if weight_method == "fit_snr": - weight = 1 / selected_variance_estimated[mask_ikbin][0] - p1d_weights[i][ikbin] = weight - # covariance_weights[i][ikbin] = np.sqrt(weight / number_in_bins) - covariance_weights[i][ikbin] = np.sqrt(weight) - else: - p1d_weights[i][ikbin] = 1 - covariance_weights[i][ikbin] = np.sqrt(1 / number_in_bins) - p1d_groups[i][ikbin] = np.nansum( - selected_pk[mask_ikbin] * covariance_weights[i][ikbin] - ) + if len(p1d_groups) == 0: + return np.full(nbins_k * nbins_k, np.nan) - # Calculation of covariance terms mean_pk = mean_p1d_table["meanPk"][(nbins_k * izbin) : (nbins_k * (izbin + 1))] mean_pk_product = np.outer(mean_pk, mean_pk) @@ -997,7 +1061,7 @@ def compute_cov( covariance_weights_product_sum = np.zeros((nbins_k, nbins_k)) weights_product_sum = np.zeros((nbins_k, nbins_k)) - for i, _ in enumerate(sub_forest_ids): + for i in range(len(p1d_groups)): p1d_groups_product_sum = np.nansum( [p1d_groups_product_sum, np.outer(p1d_groups[i], p1d_groups[i])], axis=0 ) @@ -1012,7 +1076,7 @@ def compute_cov( [weights_product_sum, np.outer(p1d_weights[i], p1d_weights[i])], axis=0 ) - del p1d_groups, covariance_weights + del p1d_groups, covariance_weights, p1d_weights covariance_matrix = (weights_product_sum / weights_sum_product) * ( (p1d_groups_product_sum / covariance_weights_product_sum) - mean_pk_product @@ -1033,16 +1097,13 @@ def compute_cov( covariance_array = np.ravel(covariance_matrix) - return zbin_array, index_zbin_array, n_array, covariance_array, k1_array, k2_array - + return covariance_array def compute_cov_not_vectorized( p1d_table, mean_p1d_table, - zbin_centers, n_chunks, - kbin_edges, k_index, nbins_k, weight_method, @@ -1061,9 +1122,6 @@ def compute_cov_not_vectorized( mean_p1d_table (array-like): Table of mean 1D power spectra, with column 'meanPk'. - zbin_centers (array-like): - Array of bin centers for redshift. - n_chunks (array-like): Array of the number of chunks in each redshift bin. @@ -1102,31 +1160,15 @@ def compute_cov_not_vectorized( covariance_array (array-like): Array of covariance coefficients. """ - zbin_array = np.zeros(nbins_k * nbins_k) - index_zbin_array = np.zeros(nbins_k * nbins_k, dtype=int) n_array = np.zeros(nbins_k * nbins_k) covariance_array = np.zeros(nbins_k * nbins_k) - k1_array = np.zeros(nbins_k * nbins_k) - k2_array = np.zeros(nbins_k * nbins_k) if weight_method == "fit_snr": weight_array = np.zeros(nbins_k * nbins_k) - kbin_centers = (kbin_edges[1:] + kbin_edges[:-1]) / 2 if n_chunks[izbin] == 0: # Fill rows with NaNs - zbin_array[:] = zbin_centers[izbin] - index_zbin_array[:] = izbin - n_array[:] = 0 covariance_array[:] = np.nan - k1_array[:] = np.nan - k2_array[:] = np.nan - return ( - zbin_array, - index_zbin_array, - n_array, - covariance_array, - k1_array, - k2_array, - ) + return covariance_array + # First loop 1) id sub-forest bins for sub_forest_id in sub_forest_ids: select_id = select_z & (p1d_table["sub_forest_id"] == sub_forest_id) @@ -1200,10 +1242,6 @@ def compute_cov_not_vectorized( # weight_array[index] if weights are normalized to give # sum(weight_array[index]) = n_array[index] covariance_array[index] = covariance_array[index] / n_array[index] - zbin_array[index] = zbin_centers[izbin] - index_zbin_array[index] = izbin - k1_array[index] = kbin_centers[ikbin] - k2_array[index] = kbin_centers[ikbin2] if ikbin2 != ikbin: # index of the (ikbin,ikbin2) coefficient on the bottom of the matrix @@ -1211,11 +1249,6 @@ def compute_cov_not_vectorized( covariance_array[index_2] = covariance_array[index] n_array[index_2] = n_array[index] - zbin_array[index_2] = zbin_centers[izbin] - index_zbin_array[index_2] = izbin - k1_array[index_2] = kbin_centers[ikbin2] - k2_array[index_2] = kbin_centers[ikbin] - # For fit_snr method, due to the SNR fitting scheme used for weighting, # the diagonal of the weigthed sample covariance matrix is not equal # to the error in mean P1D. This is tested on Ohio mocks. @@ -1251,8 +1284,7 @@ def compute_cov_not_vectorized( ) ) - return zbin_array, index_zbin_array, n_array, covariance_array, k1_array, k2_array - + return covariance_array def run_postproc_pk1d( @@ -1325,6 +1357,7 @@ def run_postproc_pk1d( compute_covariance=compute_covariance, compute_bootstrap=compute_bootstrap, number_bootstrap=number_bootstrap, + number_worker=ncpu, ) result = fitsio.FITS(output_file, "rw", clobber=True) From db1ce68c8d4ddc2c97d13fde0eb07c8257b59c7f Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Thu, 4 Jul 2024 08:01:24 -0700 Subject: [PATCH 03/11] improvement doc --- py/picca/pk1d/postproc_pk1d.py | 102 ++++++++++++++++----------------- 1 file changed, 48 insertions(+), 54 deletions(-) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index d66668baa..a2f3244b8 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -915,12 +915,14 @@ def compute_p1d_groups( Return ------ - p1d_weights_z (array-like): + p1d_weights (array-like): + Weights associated with p1d pixels for all subforest, used in the calculation of covariance. - covariance_weights_z (array-like): - - p1d_groups_z (array-like): + covariance_weights (array-like): + Weights for all subforest used inside the main covariance sum. + p1d_groups (array-like): + Individual p1d pixels grouped in the same wavenumber binning for all subforest """ select_z = (p1d_table["forest_z"] < zbin_edges[izbin + 1]) & ( @@ -965,26 +967,51 @@ def compute_p1d_groups( del p1d_los_table - return output_cov[:, 0, :], output_cov[:, 1, :], output_cov[:, 2, :] + p1d_weights, covariance_weights, p1d_groups = ( + output_cov[:, 0, :], + output_cov[:, 1, :], + output_cov[:, 2, :], + ) + return p1d_weights, covariance_weights, p1d_groups def compute_groups(nbins_k, p1d_los): + """Compute the P1D groups for one subforest. + + Arguments + --------- + nbins_k (int): + Number of k bins. + + p1d_los (array-like): + Table containing all p1d pixels unordered for one subforest + + Return + ------ + p1d_weights_id (array-like): + Weights associated with p1d pixels for one subforest, used in the calculation of covariance. + + covariance_weights_id (array-like): + Weights for one subforest used inside the main covariance sum. - p1d_weights = np.zeros(nbins_k) - covariance_weights = np.zeros(nbins_k) - p1d_groups = np.zeros(nbins_k) + p1d_groups_id (array-like): + Individual p1d pixels grouped in the same wavenumber binning for one subforest + """ + p1d_weights_id = np.zeros(nbins_k) + covariance_weights_id = np.zeros(nbins_k) + p1d_groups_id = np.zeros(nbins_k) for ikbin in range(nbins_k): mask_ikbin = p1d_los["k_index"] == ikbin number_in_bins = len(mask_ikbin[mask_ikbin]) if number_in_bins != 0: weight = p1d_los["weight"][mask_ikbin][0] - p1d_weights[ikbin] = weight - covariance_weights[ikbin] = np.sqrt(weight / number_in_bins) - p1d_groups[ikbin] = np.nansum( - p1d_los["pk"][mask_ikbin] * covariance_weights[ikbin] + p1d_weights_id[ikbin] = weight + covariance_weights_id[ikbin] = np.sqrt(weight / number_in_bins) + p1d_groups_id[ikbin] = np.nansum( + p1d_los["pk"][mask_ikbin] * covariance_weights_id[ikbin] ) - return p1d_weights, covariance_weights, p1d_groups + return p1d_weights_id, covariance_weights_id, p1d_groups_id def compute_cov( @@ -1000,50 +1027,26 @@ def compute_cov( Arguments --------- - p1d_table (array-like): - Table of 1D power spectra, with columns 'Pk' and 'sub_forest_id'. + weight_method: str, + Method to weight the data. mean_p1d_table (array-like): Table of mean 1D power spectra, with column 'meanPk'. - zbin_centers (array-like): - Array of bin centers for redshift. - - n_chunks (array-like): - Array of the number of chunks in each redshift bin. - - k_index (array-like): - Array of indices for k-values, with -1 indicating values outside of the k bins. - nbins_k (int): Number of k bins. - weight_method: str, - Method to weight the data. - - snrfit_table: numpy ndarray, - Table containing SNR fit infos - izbin (int): Current redshift bin being considered. - select_z (array-like): - Boolean array for selecting data points based on redshift. + p1d_weights (array-like): - sub_forest_ids (array-like): - Array of chunk ids. + covariance_weights (array-like): + + p1d_groups (array-like): Return ------ - zbin_array (array-like): - Array of redshift bin centers for each covariance coefficient. - - index_zbin_array (array-like): - Array of redshift bin indices for each covariance coefficient. - - n_array (array-like): - Array of the number of power spectra used to compute each covariance coefficient. - covariance_array (array-like): Array of covariance coefficients. """ @@ -1061,9 +1064,9 @@ def compute_cov( covariance_weights_product_sum = np.zeros((nbins_k, nbins_k)) weights_product_sum = np.zeros((nbins_k, nbins_k)) - for i in range(len(p1d_groups)): + for i, p1d_group in enumerate(p1d_groups): p1d_groups_product_sum = np.nansum( - [p1d_groups_product_sum, np.outer(p1d_groups[i], p1d_groups[i])], axis=0 + [p1d_groups_product_sum, np.outer(p1d_group, p1d_group)], axis=0 ) covariance_weights_product_sum = np.nansum( [ @@ -1148,15 +1151,6 @@ def compute_cov_not_vectorized( Return ------ - zbin_array (array-like): - Array of redshift bin centers for each covariance coefficient. - - index_zbin_array (array-like): - Array of redshift bin indices for each covariance coefficient. - - n_array (array-like): - Array of the number of power spectra used to compute each covariance coefficient. - covariance_array (array-like): Array of covariance coefficients. """ From d71221ed0c11ad46220208d7d2af0a2ed184dadc Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Thu, 4 Jul 2024 08:19:01 -0700 Subject: [PATCH 04/11] pylint fix --- py/picca/pk1d/postproc_pk1d.py | 223 ++++++++++++++++++++++----------- 1 file changed, 149 insertions(+), 74 deletions(-) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index a2f3244b8..e08a792a3 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -21,11 +21,12 @@ import numpy as np from astropy.stats import bootstrap from astropy.table import Table, vstack +from scipy.optimize import curve_fit +from scipy.stats import binned_statistic from picca.constants import ABSORBER_IGM, SPEED_LIGHT from picca.pk1d.utils import MEANPK_FITRANGE_SNR, fitfunc_variance_pk1d from picca.utils import userprint -from scipy.optimize import curve_fit -from scipy.stats import binned_statistic + def read_pk1d(filename, kbin_edges, snrcut=None, zbins_snrcut=None): @@ -392,77 +393,18 @@ def compute_mean_pk1d( ] ) - if compute_covariance: - userprint("Computation of covariance matrix") - - func = partial( - compute_cov, - weight_method, - mean_p1d_table, - nbins_k, - ) - if number_worker == 1: - output_cov = [func(*p1d_group) for p1d_group in p1d_groups] - else: - with Pool(number_worker) as pool: - output_cov = pool.starmap(func, p1d_groups) - - for izbin in range(nbins_z): - covariance_array = output_cov[izbin] - i_min = izbin * nbins_k * nbins_k - i_max = (izbin + 1) * nbins_k * nbins_k - cov_table["covariance"][i_min:i_max] = covariance_array - - if compute_bootstrap: - userprint("Computing covariance matrix with bootstrap method") - p1d_groups_bootstrap = [] - for izbin in range(nbins_z): - number_sub_forests = len(p1d_groups[izbin][1]) - if number_sub_forests > 0: - bootid = np.array( - bootstrap(np.arange(number_sub_forests), number_bootstrap) - ).astype(int) - else: - bootid = np.full(number_bootstrap, None) - - for iboot in range(number_bootstrap): - if bootid[iboot] is None: - p1d_weights_z, covariance_weights_z, p1d_groups_z = [], [], [] - else: - p1d_weights_z = p1d_groups[izbin][1][bootid[iboot]] - covariance_weights_z = p1d_groups[izbin][2][bootid[iboot]] - p1d_groups_z = p1d_groups[izbin][3][bootid[iboot]] - p1d_groups_bootstrap.append( - [ - izbin, - p1d_weights_z, - covariance_weights_z, - p1d_groups_z, - ] - ) - - func = partial( - compute_cov, - weight_method, - mean_p1d_table, - nbins_k, - ) - if number_worker == 1: - output_cov = [func(*p) for p in p1d_groups_bootstrap] - else: - with Pool(number_worker) as pool: - output_cov = pool.starmap(func, p1d_groups_bootstrap) - - for izbin in range(nbins_z): - boot_cov = [] - for iboot in range(number_bootstrap): - covariance_array = (output_cov[izbin * number_bootstrap + iboot],) - boot_cov.append(covariance_array) - - i_min = izbin * nbins_k * nbins_k - i_max = (izbin + 1) * nbins_k * nbins_k - cov_table["boot_covariance"][i_min:i_max] = np.mean(boot_cov, axis=0) - cov_table["error_boot_covariance"][i_min:i_max] = np.std(boot_cov, axis=0) + compute_and_fill_covariance( + compute_covariance, + compute_bootstrap, + weight_method, + mean_p1d_table, + nbins_k, + nbins_z, + p1d_groups, + number_worker, + number_bootstrap, + cov_table, + ) if output_snrfit is not None: np.savetxt( @@ -477,6 +419,8 @@ def compute_mean_pk1d( return mean_p1d_table, metadata_table, cov_table + + def fill_average_pk( nbins_z, nbins_k, @@ -874,6 +818,134 @@ def compute_average_pk_redshift( ) +def compute_and_fill_covariance( + compute_covariance, + compute_bootstrap, + weight_method, + mean_p1d_table, + nbins_k, + nbins_z, + p1d_groups, + number_worker, + number_bootstrap, + cov_table, +): + """Compute the covariance and bootstrap covariance and fill the corresponding + cov_table variable + + compute_covariance: Bool + If True, compute statistical covariance matrix of the mean P1D. + Requires p1d_table to contain 'sub_forest_id', since correlations are computed + within individual forests. + + compute_bootstrap: Bool + If True, compute statistical covariance using a simple bootstrap method. + + weight_method: str, + Method to weight the data. + + mean_p1d_table (array-like): + Table of mean 1D power spectra, with column 'meanPk'. + + nbins_k (int): + Number of k bins. + + nbins_z (int): + Number of z bins. + + p1d_groups (array-like): + Individual p1d pixels grouped in the same wavenumber binning for all subforest + + number_worker: int + Calculations of mean P1Ds and covariances are run parallel over redshift bins. + + number_bootstrap: int + Number of bootstrap samples used if compute_bootstrap is True. + + cov_table (array-like): + Covariance table to fill. + + Return + ------ + None + """ + + + if compute_covariance: + userprint("Computation of covariance matrix") + + func = partial( + compute_cov, + weight_method, + mean_p1d_table, + nbins_k, + ) + if number_worker == 1: + output_cov = [func(*p1d_group) for p1d_group in p1d_groups] + else: + with Pool(number_worker) as pool: + output_cov = pool.starmap(func, p1d_groups) + + for izbin in range(nbins_z): + covariance_array = output_cov[izbin] + i_min = izbin * nbins_k * nbins_k + i_max = (izbin + 1) * nbins_k * nbins_k + cov_table["covariance"][i_min:i_max] = covariance_array + + if compute_bootstrap: + userprint("Computing covariance matrix with bootstrap method") + p1d_groups_bootstrap = [] + for izbin in range(nbins_z): + number_sub_forests = len(p1d_groups[izbin][1]) + if number_sub_forests > 0: + bootid = np.array( + bootstrap(np.arange(number_sub_forests), number_bootstrap) + ).astype(int) + else: + bootid = np.full(number_bootstrap, None) + + for iboot in range(number_bootstrap): + if bootid[iboot] is None: + p1d_weights_z, covariance_weights_z, p1d_groups_z = [], [], [] + else: + p1d_weights_z = p1d_groups[izbin][1][bootid[iboot]] + covariance_weights_z = p1d_groups[izbin][2][bootid[iboot]] + p1d_groups_z = p1d_groups[izbin][3][bootid[iboot]] + p1d_groups_bootstrap.append( + [ + izbin, + p1d_weights_z, + covariance_weights_z, + p1d_groups_z, + ] + ) + + func = partial( + compute_cov, + weight_method, + mean_p1d_table, + nbins_k, + ) + if number_worker == 1: + output_cov = [func(*p) for p in p1d_groups_bootstrap] + else: + with Pool(number_worker) as pool: + output_cov = pool.starmap(func, p1d_groups_bootstrap) + + for izbin in range(nbins_z): + boot_cov = [] + for iboot in range(number_bootstrap): + covariance_array = (output_cov[izbin * number_bootstrap + iboot],) + boot_cov.append(covariance_array) + + i_min = izbin * nbins_k * nbins_k + i_max = (izbin + 1) * nbins_k * nbins_k + cov_table["boot_covariance"][i_min:i_max] = np.mean(boot_cov, axis=0) + cov_table["error_boot_covariance"][i_min:i_max] = np.std(boot_cov, axis=0) + + + + def compute_p1d_groups( weight_method, nbins_k, @@ -1039,11 +1111,14 @@ def compute_cov( izbin (int): Current redshift bin being considered. - p1d_weights (array-like): + p1d_weights (array-like): + Weights associated with p1d pixels for all subforest, used in the calculation of covariance. covariance_weights (array-like): + Weights for all subforest used inside the main covariance sum. p1d_groups (array-like): + Individual p1d pixels grouped in the same wavenumber binning for all subforest Return ------ From f2120b8e0fa941324a7386635d97f8d77cd70e72 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Thu, 4 Jul 2024 08:29:25 -0700 Subject: [PATCH 05/11] test fix --- bin/picca_Pk1D.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bin/picca_Pk1D.py b/bin/picca_Pk1D.py index 44dc9c9a3..4a22d16ac 100755 --- a/bin/picca_Pk1D.py +++ b/bin/picca_Pk1D.py @@ -192,8 +192,8 @@ def process_all_files(index_file_args): ) else: pk_noise=pk_diff=np.zeros(pk_raw.shape) - fft_delta_noise = np.zeros(pk_raw.shape,dtype=np.complex_) - fft_delta_diff = np.zeros(pk_raw.shape,dtype=np.complex_) + fft_delta_noise = np.zeros(pk_raw.shape,dtype=np.complex128) + fft_delta_diff = np.zeros(pk_raw.shape,dtype=np.complex128) # Compute resolution correction, needs uniform binning if not running_on_raw_transmission: From 47248bb634d10a93b9c084234149549a22e5ff05 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Fri, 5 Jul 2024 02:06:02 -0700 Subject: [PATCH 06/11] fixing test --- py/picca/pk1d/postproc_pk1d.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index e08a792a3..4dbe6f46c 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -393,18 +393,18 @@ def compute_mean_pk1d( ] ) - compute_and_fill_covariance( - compute_covariance, - compute_bootstrap, - weight_method, - mean_p1d_table, - nbins_k, - nbins_z, - p1d_groups, - number_worker, - number_bootstrap, - cov_table, - ) + compute_and_fill_covariance( + compute_covariance, + compute_bootstrap, + weight_method, + mean_p1d_table, + nbins_k, + nbins_z, + p1d_groups, + number_worker, + number_bootstrap, + cov_table, + ) if output_snrfit is not None: np.savetxt( From 5a69c6cf98481edcb02c1adf72497f13a7a30b2f Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Fri, 5 Jul 2024 07:20:46 -0700 Subject: [PATCH 07/11] implementing Michael comments --- py/picca/pk1d/postproc_pk1d.py | 150 +++++++++++++++++---------------- 1 file changed, 79 insertions(+), 71 deletions(-) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index 4dbe6f46c..f0fb1c24c 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -28,7 +28,6 @@ from picca.utils import userprint - def read_pk1d(filename, kbin_edges, snrcut=None, zbins_snrcut=None): """Read Pk1D data from a single file. @@ -139,6 +138,13 @@ def read_pk1d(filename, kbin_edges, snrcut=None, zbins_snrcut=None): return p1d_table, z_array +def mean_p1d_table_regular_slice(izbin, nbins_k): + return izbin * nbins_k, (izbin + 1) * nbins_k + +def cov_table_regular_slice(izbin, nbins_k): + return izbin * nbins_k * nbins_k, (izbin + 1) * nbins_k * nbins_k + + def compute_mean_pk1d( p1d_table, z_array, @@ -357,19 +363,19 @@ def compute_mean_pk1d( k1_array, k2_array = np.meshgrid(kbin_centers, kbin_centers, indexing="ij") k1_array = np.ravel(k1_array) k2_array = np.ravel(k2_array) + select = mean_p1d_table["index_zbin"] == izbin n_array = np.ravel( - np.outer( - mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], - mean_p1d_table["N"][(nbins_k * izbin) : (nbins_k * (izbin + 1))], - ) + np.outer(mean_p1d_table["N"][select], mean_p1d_table["N"][select]) ) - i_min = izbin * nbins_k * nbins_k - i_max = (izbin + 1) * nbins_k * nbins_k - cov_table["zbin"][i_min:i_max] = zbin_array - cov_table["index_zbin"][i_min:i_max] = index_zbin_array - cov_table["k1"][i_min:i_max] = k1_array - cov_table["k2"][i_min:i_max] = k2_array - cov_table["N"][i_min:i_max] = n_array + index_cov = cov_table_regular_slice(izbin, nbins_k) + cov_table["zbin"][index_cov[0]:index_cov[1]] = zbin_array + cov_table["index_zbin"][index_cov[0]:index_cov[1]] = index_zbin_array + cov_table["k1"][index_cov[0]:index_cov[1]] = k1_array + cov_table["k2"][index_cov[0]:index_cov[1]] = k2_array + cov_table["N"][index_cov[0]:index_cov[1]] = n_array + index_mean = mean_p1d_table_regular_slice(izbin, nbins_k) + mean_pk = mean_p1d_table["meanPk"][index_mean[0] : index_mean[1]] + error_pk = mean_p1d_table["errorPk"][index_mean[0] : index_mean[1]] if n_chunks[izbin] == 0: p1d_weights_z, covariance_weights_z, p1d_groups_z = [], [], [] @@ -384,9 +390,11 @@ def compute_mean_pk1d( snrfit_table, number_worker, ) + p1d_groups.append( [ - izbin, + mean_pk, + error_pk, p1d_weights_z, covariance_weights_z, p1d_groups_z, @@ -397,7 +405,6 @@ def compute_mean_pk1d( compute_covariance, compute_bootstrap, weight_method, - mean_p1d_table, nbins_k, nbins_z, p1d_groups, @@ -419,8 +426,6 @@ def compute_mean_pk1d( return mean_p1d_table, metadata_table, cov_table - - def fill_average_pk( nbins_z, nbins_k, @@ -478,21 +483,20 @@ def fill_average_pk( median_array, snrfit_array, ) = (*output_mean[izbin],) - i_min = izbin * nbins_k - i_max = (izbin + 1) * nbins_k + index_mean = mean_p1d_table_regular_slice(izbin,nbins_k) - mean_p1d_table["zbin"][i_min:i_max] = zbin_array - mean_p1d_table["index_zbin"][i_min:i_max] = index_zbin_array - mean_p1d_table["N"][i_min:i_max] = n_array + mean_p1d_table["zbin"][index_mean[0]:index_mean[1]] = zbin_array + mean_p1d_table["index_zbin"][index_mean[0]:index_mean[1]] = index_zbin_array + mean_p1d_table["N"][index_mean[0]:index_mean[1]] = n_array for icol, col in enumerate(p1d_table_cols): - mean_p1d_table["mean" + col][i_min:i_max] = mean_array[icol] - mean_p1d_table["error" + col][i_min:i_max] = error_array[icol] - mean_p1d_table["min" + col][i_min:i_max] = min_array[icol] - mean_p1d_table["max" + col][i_min:i_max] = max_array[icol] + mean_p1d_table["mean" + col][index_mean[0]:index_mean[1]] = mean_array[icol] + mean_p1d_table["error" + col][index_mean[0]:index_mean[1]] = error_array[icol] + mean_p1d_table["min" + col][index_mean[0]:index_mean[1]] = min_array[icol] + mean_p1d_table["max" + col][index_mean[0]:index_mean[1]] = max_array[icol] if not nomedians: - mean_p1d_table["median" + col][i_min:i_max] = median_array[icol] + mean_p1d_table["median" + col][index_mean[0]:index_mean[1]] = median_array[icol] if weight_method == "fit_snr": - snrfit_table[i_min:i_max, :] = snrfit_array + snrfit_table[index_mean[0]:index_mean[1], :] = snrfit_array def compute_average_pk_redshift( @@ -822,7 +826,6 @@ def compute_and_fill_covariance( compute_covariance, compute_bootstrap, weight_method, - mean_p1d_table, nbins_k, nbins_z, p1d_groups, @@ -844,9 +847,6 @@ def compute_and_fill_covariance( weight_method: str, Method to weight the data. - mean_p1d_table (array-like): - Table of mean 1D power spectra, with column 'meanPk'. - nbins_k (int): Number of k bins. @@ -870,14 +870,12 @@ def compute_and_fill_covariance( None """ - if compute_covariance: userprint("Computation of covariance matrix") func = partial( compute_cov, weight_method, - mean_p1d_table, nbins_k, ) if number_worker == 1: @@ -888,9 +886,8 @@ def compute_and_fill_covariance( for izbin in range(nbins_z): covariance_array = output_cov[izbin] - i_min = izbin * nbins_k * nbins_k - i_max = (izbin + 1) * nbins_k * nbins_k - cov_table["covariance"][i_min:i_max] = covariance_array + index_cov = cov_table_regular_slice(izbin, nbins_k) + cov_table["covariance"][index_cov[0]:index_cov[1]] = covariance_array if compute_bootstrap: userprint("Computing covariance matrix with bootstrap method") @@ -906,14 +903,23 @@ def compute_and_fill_covariance( for iboot in range(number_bootstrap): if bootid[iboot] is None: - p1d_weights_z, covariance_weights_z, p1d_groups_z = [], [], [] + ( + mean_pk, + error_pk, + p1d_weights_z, + covariance_weights_z, + p1d_groups_z, + ) = ([], [], [], [], []) else: - p1d_weights_z = p1d_groups[izbin][1][bootid[iboot]] - covariance_weights_z = p1d_groups[izbin][2][bootid[iboot]] - p1d_groups_z = p1d_groups[izbin][3][bootid[iboot]] + mean_pk = p1d_groups[izbin][0][bootid[iboot]] + error_pk = p1d_groups[izbin][1][bootid[iboot]] + p1d_weights_z = p1d_groups[izbin][2][bootid[iboot]] + covariance_weights_z = p1d_groups[izbin][3][bootid[iboot]] + p1d_groups_z = p1d_groups[izbin][4][bootid[iboot]] p1d_groups_bootstrap.append( [ - izbin, + mean_pk, + error_pk, p1d_weights_z, covariance_weights_z, p1d_groups_z, @@ -923,7 +929,6 @@ def compute_and_fill_covariance( func = partial( compute_cov, weight_method, - mean_p1d_table, nbins_k, ) if number_worker == 1: @@ -938,12 +943,9 @@ def compute_and_fill_covariance( covariance_array = (output_cov[izbin * number_bootstrap + iboot],) boot_cov.append(covariance_array) - i_min = izbin * nbins_k * nbins_k - i_max = (izbin + 1) * nbins_k * nbins_k - cov_table["boot_covariance"][i_min:i_max] = np.mean(boot_cov, axis=0) - cov_table["error_boot_covariance"][i_min:i_max] = np.std(boot_cov, axis=0) - - + index_cov = cov_table_regular_slice(izbin, nbins_k) + cov_table["boot_covariance"][index_cov[0]:index_cov[1]] = np.mean(boot_cov, axis=0) + cov_table["error_boot_covariance"][index_cov[0]:index_cov[1]] = np.std(boot_cov, axis=0) def compute_p1d_groups( @@ -1009,7 +1011,8 @@ def compute_p1d_groups( names=("sub_forest_id", "k_index", "pk"), ) if weight_method == "fit_snr": - snrfit_z = snrfit_table[izbin * nbins_k : (izbin + 1) * nbins_k, :] + index_slice = mean_p1d_table_regular_slice(izbin, nbins_k) + snrfit_z = snrfit_table[index_slice[0] : index_slice[1], :] p1d_sub_table["weight"] = 1 / fitfunc_variance_pk1d( p1d_table["forest_snr"][select_z], snrfit_z[k_index[select_z], 2], @@ -1028,7 +1031,7 @@ def compute_p1d_groups( del grouped_table func = partial( - compute_groups, + compute_groups_for_one_forest, nbins_k, ) if number_worker == 1: @@ -1047,7 +1050,7 @@ def compute_p1d_groups( return p1d_weights, covariance_weights, p1d_groups -def compute_groups(nbins_k, p1d_los): +def compute_groups_for_one_forest(nbins_k, p1d_los): """Compute the P1D groups for one subforest. Arguments @@ -1088,28 +1091,32 @@ def compute_groups(nbins_k, p1d_los): def compute_cov( weight_method, - mean_p1d_table, nbins_k, - izbin, + mean_pk, + error_pk, p1d_weights, covariance_weights, p1d_groups, ): """Compute the covariance of a set of 1D power spectra. + This is a new version of the covariance calculation. + It needs that the input data are expressed on the same + wavenumber grid. The calculation is then performed in + an entire vectorized way. Arguments --------- weight_method: str, Method to weight the data. - mean_p1d_table (array-like): - Table of mean 1D power spectra, with column 'meanPk'. - nbins_k (int): Number of k bins. - izbin (int): - Current redshift bin being considered. + mean_pk (array-like): + Mean 1D power spectra, for the considered redshift bin. + + error_pk (array-like): + Standard deviation of the 1D power spectra, for the considered redshift bin. p1d_weights (array-like): Weights associated with p1d pixels for all subforest, used in the calculation of covariance. @@ -1129,35 +1136,34 @@ def compute_cov( if len(p1d_groups) == 0: return np.full(nbins_k * nbins_k, np.nan) - mean_pk = mean_p1d_table["meanPk"][(nbins_k * izbin) : (nbins_k * (izbin + 1))] mean_pk_product = np.outer(mean_pk, mean_pk) sum_p1d_weights = np.nansum(p1d_weights, axis=0) - weights_sum_product = np.outer(sum_p1d_weights, sum_p1d_weights) + weights_sum_of_product = np.outer(sum_p1d_weights, sum_p1d_weights) p1d_groups_product_sum = np.zeros((nbins_k, nbins_k)) - covariance_weights_product_sum = np.zeros((nbins_k, nbins_k)) - weights_product_sum = np.zeros((nbins_k, nbins_k)) + covariance_weights_product_of_sum = np.zeros((nbins_k, nbins_k)) + weights_product_of_sum = np.zeros((nbins_k, nbins_k)) for i, p1d_group in enumerate(p1d_groups): p1d_groups_product_sum = np.nansum( [p1d_groups_product_sum, np.outer(p1d_group, p1d_group)], axis=0 ) - covariance_weights_product_sum = np.nansum( + covariance_weights_product_of_sum = np.nansum( [ - covariance_weights_product_sum, + covariance_weights_product_of_sum, np.outer(covariance_weights[i], covariance_weights[i]), ], axis=0, ) - weights_product_sum = np.nansum( - [weights_product_sum, np.outer(p1d_weights[i], p1d_weights[i])], axis=0 + weights_product_of_sum = np.nansum( + [weights_product_of_sum, np.outer(p1d_weights[i], p1d_weights[i])], axis=0 ) del p1d_groups, covariance_weights, p1d_weights - covariance_matrix = (weights_product_sum / weights_sum_product) * ( - (p1d_groups_product_sum / covariance_weights_product_sum) - mean_pk_product + covariance_matrix = (weights_product_of_sum / weights_sum_of_product) * ( + (p1d_groups_product_sum / covariance_weights_product_of_sum) - mean_pk_product ) # For fit_snr method, due to the SNR fitting scheme used for weighting, @@ -1165,11 +1171,10 @@ def compute_cov( # to the error in mean P1D. This is tested on Ohio mocks and data. # We choose to renormalize the whole covariance matrix. if weight_method == "fit_snr": - std_pk = mean_p1d_table["errorPk"][(nbins_k * izbin) : (nbins_k * (izbin + 1))] covariance_diag = np.diag(covariance_matrix) covariance_matrix = ( covariance_matrix - * np.outer(std_pk, std_pk) + * np.outer(error_pk, error_pk) / np.sqrt(np.outer(covariance_diag, covariance_diag)) ) @@ -1191,6 +1196,8 @@ def compute_cov_not_vectorized( sub_forest_ids, ): """Compute the covariance of a set of 1D power spectra. + This is an old implementation which is summing up each individual mode every time. + This is very slow for large data samples. Arguments --------- @@ -1251,7 +1258,8 @@ def compute_cov_not_vectorized( # adapted with weights w_i = np.sqrt(w_ipk * w_ipk2) for covariance # to obtain the right definition of variance on the diagonal selected_snr = p1d_table["forest_snr"][select_id] - snrfit_z = snrfit_table[izbin * nbins_k : (izbin + 1) * nbins_k, :] + index_slice = mean_p1d_table_regular_slice(izbin, nbins_k) + snrfit_z = snrfit_table[index_slice[0] : index_slice[1], :] selected_variance_estimated = fitfunc_variance_pk1d( selected_snr, snrfit_z[selected_ikbin, 2], snrfit_z[selected_ikbin, 3] ) From 23cd76319d5c2a4fad6c2f1904f023f6737381fb Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Fri, 5 Jul 2024 07:35:02 -0700 Subject: [PATCH 08/11] pylint fix --- py/picca/pk1d/postproc_pk1d.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index f0fb1c24c..f355bc9d1 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -139,9 +139,39 @@ def read_pk1d(filename, kbin_edges, snrcut=None, zbins_snrcut=None): def mean_p1d_table_regular_slice(izbin, nbins_k): + """Return the arguments of a slice of mean P1D table for a given redshift. + + Arguments + --------- + izbin (int): + Current redshift bin being considered. + + nbins_k (int): + Number of k bins. + + Return + ------ + Arguments of a slice of mean P1D table. + """ return izbin * nbins_k, (izbin + 1) * nbins_k + + def cov_table_regular_slice(izbin, nbins_k): + """Return the arguments of a slice of covariance table for a given redshift. + + Arguments + --------- + izbin (int): + Current redshift bin being considered. + + nbins_k (int): + Number of k bins. + + Return + ------ + Arguments of a slice of covariance table. + """ return izbin * nbins_k * nbins_k, (izbin + 1) * nbins_k * nbins_k From 060f9df2870a084165c3922fd2249f6af5006db1 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 15 Jul 2024 06:01:59 -0700 Subject: [PATCH 09/11] adding comments --- py/picca/pk1d/postproc_pk1d.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/py/picca/pk1d/postproc_pk1d.py b/py/picca/pk1d/postproc_pk1d.py index f355bc9d1..03b501658 100644 --- a/py/picca/pk1d/postproc_pk1d.py +++ b/py/picca/pk1d/postproc_pk1d.py @@ -1051,6 +1051,7 @@ def compute_p1d_groups( else: p1d_sub_table["weight"] = 1 + # Remove bins that are not associated with any wavenumber bin. p1d_sub_table = p1d_sub_table[p1d_sub_table["k_index"] >= 0] grouped_table = p1d_sub_table.group_by("sub_forest_id") @@ -1176,6 +1177,9 @@ def compute_cov( weights_product_of_sum = np.zeros((nbins_k, nbins_k)) for i, p1d_group in enumerate(p1d_groups): + # The summation is done with np.nansum instead of simple addition to not + # include the NaN that are present in the individual p1d. + # The summation is not done at the end, to prevent memory overhead. p1d_groups_product_sum = np.nansum( [p1d_groups_product_sum, np.outer(p1d_group, p1d_group)], axis=0 ) From eb28700fb0e8f09a6935e1ade0decbda4f37aee3 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 15 Jul 2024 06:37:48 -0700 Subject: [PATCH 10/11] adding covariance tests --- .../test_Pk1D/meanPk1D_covariance.fits.gz | Bin 0 -> 9172 bytes py/picca/tests/test_2_pk1d.py | 31 ++++++++++++++++++ 2 files changed, 31 insertions(+) create mode 100644 py/picca/tests/data/test_Pk1D/meanPk1D_covariance.fits.gz diff --git a/py/picca/tests/data/test_Pk1D/meanPk1D_covariance.fits.gz b/py/picca/tests/data/test_Pk1D/meanPk1D_covariance.fits.gz new file mode 100644 index 0000000000000000000000000000000000000000..745f1d3579050e07a79cf088444691d6afc6a53d GIT binary patch literal 9172 zcmY+Kby!qw*Y*YJlx_wPk&;e{0YL$g5Ts)ON$KtxX_Q7KrA0ubk?tBgq$G##8XBfP zyzl3IpZA-;_HoU=&b8KWoqNr3>}voq1o!@aaPA)l`0|4{3DB!6 zM34-m!E@;>AbZ=tZM1PL@OV*LI>ThaD1FWCg%n;dVNv`UF$vdNPo?qT;_GM*CfnoA(7J7ydVih2Ql7q0} zq{Glkh=tl8ig*0(C>m_luN6j?K%x}=JUw3Lm08v+=T~Vfq!+1+zu1cjtT+Y2A4pCe zsy?&ej;Ci*BZ(s*Auj9&SvxvAzgrKlQzd&NLCn>M7f_I#9`ykELFzFjyBLRTEcFM^ z*t?t}hx-g(pOFMnptL^u(>?t*!NASXdv>QS4CVeG82t zBaSCxU#q>ic;mLJbTZ#`wpH*n01y zRfsTdA2q$M#s6Ha0Cg>PJbq6s&h!x!)s8Rr5w?#HzWIQs$kO_h_)+4TTT*ylOxLQT zY6kB)K&2Vn$B;W)g!f~|TXa-Lz~g-1l*D^#-|h*umin}$P`HPsI>C3l7|U^&0ly5A ziX3loFYT5z1z7+4yrICC2@^33uN#x$we)MN=mrp5M6WW58Ncp(2epZyIuq-CSug3) zb1*xHo49+lxw6#6;|GVro8=+Ooa67ugG^-^v$^JBmD=Z*_B(ieIn=Z?@w!|T2G4Pm z$A@m&tW00g5BN%u(%()D3@=U%_9U?b<(#QqH_-ipiS!5RN5V-fRPU9VwaA|@dwLk# z+Pxf&`u)zUm=py-T@!Qp66ucz{5OvzIFG5V^^@f^H#)gLMT)k{L~uEuKdjY!MsQ?1hW-2V{6`z z0DkcIZ@75*`6-CC;90sdTP01A~cx&knQ$ zj86Z1muPN05Y~}L(a+od-#}-@NP;{c#ca-J!vsn%w(0(>pb?jMM(htpDkfq9YfriV zNlX!(p@e*fjY%e#%RzgCC#<`H{eH|FMbEcXqESk%d1ojX(hw2F#gGt&tQ|l`E1iX9 zb36q9L@9lWg(`xdQiH91@+?Crl$|1l)}Zr1v4Ow4WqQan0kt>`@VAG63~O1k)5(YX zD%pVA1NsEu3_8U8Pd;js4Hzr}JVjcJB>aNdC&Q)s$KDU$?GmGA*i8M;iYsdbz#n{# zcY9+w8Mb;8@$WKJD`pgEPbe}@UMBmJ9f6N~mzq5?EQO`}Ca-rMIL^zQ1uzyvcV++_ zF}W-5+Q0xzZeBlk`(4g4-RL~=1Ax3Grxv9!2T3Md>~+;8IT0)V1rt00NFegJf?BBq zfN>fG*%=Qouj)&a%%wl$v=%Y(YRsA{4Xjz)L`aI@37aVc?sRA-@yyfXD!se zcK-vsTd;`vA9#$I+@gQsMcsuKJpgd^tf+?4Zn4Im4b20&{g|g*dw^)Re$~6i=5$H=Y8z@@(RUHYS1Tbq@z+oQ4#SBO zUxH5mwjB}c_W5M*Urmy2d}II5&sL##L3)v4Pwq?@77A5F-Sz%IHJgzkqK+ydiF0vi zB*UFK@-|ZMOp1DOw>vcU{4S%WJNRb9K(qza4_x^#qwar|lv3}iQ(p(1|7*X$eSrbc z!l$H9|H1zE+p8i736*unKdu`C{gYS0Rf6UJmsi0RIqjWO{(ABjwF#hq{^G{w5U>;D z^)&$e^b(Q&8gKtk9_EZl;piZ*IyI)vFf_J;uMXCBa*sPGXQKgP&aGOSRvN>}89LnH z;Oi;2pzy-CZj@)}P4RqDUDbgJ<#p?r#RFeMCT4F+nwN2n#4Lvgm+P~@5Yg7P+OOh= zyQ}HPih4`D+)2T-r}6E{PwoYEbeYs7cdLR{t7WQLEFFT~0#f$vy4s#>F?2(D;23T+HjyOz{skCGfL3l6Uz=EF#-Hlu84!^?>Af`XlZ)gBk;@}XFrY( zV2yn>r2(;!hsgvUv^Id9Lw$H@<5j*{r}L>p!qt{XI?HpyvWbI;Lb^KXJo-5V5*zhz z!7VYw{XQopGW~M62ie*WEL1uaX<8g>r_P;D0uhV(*~IV__{V%YI>&gWu_hQUl{*f1 zj^(>rHw6Bz{vGizEMbq9_OIEmhIGSa0qWY3-xdzpoMm!+nE67La- z*x+o@(B1bfS$NzSzldWps4gL^zk_Xax%m5WwND7LVeQOc!?HiWOUQKIN8Mgntu3gI zLs>-sn&1t-N>}jNhL%%X4NFMYe0`s5Gg;%p{^gHfpnL8}=k~RhFX`K6Z^Yb;CY0Ue z6l)JG6(82L7|b-L$iHqfB-{f)CHRA;qWk9#^mmAK!fOEB?fT8`HK*BV9*15sVxNmJ zc#0<&b8SuL%qZ8k8^EMU6Be4Gn8e)BzWjsW${hPe)(r8m8P~9>(e1-mr!J%(g`@8r z%7(B$_RY?;>o4I}dv5rL4&=*57ZMbYrd&OEeYF$|)tT^aS)n1WrH?(k)t-HLRmfssj zS;Tmm8EyI3Piaw$PU6hUg`LNYo`?*|bQ!#@b03-Vb|W%S)KJxzSqn%Vr%!Rae|s(e zasOr}I@)A0tLs;!GioJi4*)+}7wAyd+D~h9F3rpTY(JiAM_Wnwfg1j7?7opF&IRnq zfABvL35m#pCQ%-pG~CFPr##*DbQ`*3XAvwRWv*hc6cszMNfBS++TirPHX359j&03P z)tjZY3hwjJm%(vY<2u3RRh6$uKT8RNMx37s1(yzAtt^_r6{97{ua^!B$!yBs_c%Z9 zyV88vS*QbrjU^cp+6+lFIkbPm+G6-XD6O^M@fj}7gpu|j9d6oPUHn8thP^VZA9XwRr~*$_LLACjYSK9nYeQrL zY{UEI%$40F=<@!D^Hq-Q95pf2hB(^6H6~0X-qv<$CL!_5j6~RRID706*&}(`WYQ*m zn)&{Z!67a2MPD~u{X(J{2z?Cjg4wlx!2M;G=@PHRSF>Q0lp7cH1dAImvR>M$i!Lms zlk3>VOD&LP|C4$+2+haf7C#%CxEJi1*Fy`FY|P(Qd@&Tl7fy<`GCoqGdpe}U8(7p9 z*jDT-x{*o$$WvMsKm55pZXF?ug#plB)Vb<&Mg!mH{mNkEO(&h3S#W{cTp3CT?&cAp zqsudOHs`K>8;9p=t~&`W_mrS}T>h2VqXl^3I{=K37@ARKkZ4xy{4f*b(xiA^CISn~ z4`kh2o@~Z8Uoewk?okONU3;;37QVSkS)xi+E)>k2+c@B*86nR-sA|@6QAyMQkWKk& zG>M%ouH&0BxD$3-K~A-2n~jHPPWC_a+G0C(_O;?;$?m?nO04L}^}=J+oc0)AXo{J>8ToGlu*kF!a+8mw`J z$nj-<#dX=1GmyeHO-f&Pk*j9#BmWXqhTG+crx-$kR;cWHs19Ab^8H2YEsOsC1qnUd)^ zzv~Bs$X{~h2nE(|vHG>Qb%>R?0{&Ebbd*@FGuopr+rms|$W8+JsxZ~vj}j|qTOcJ< zcJ~a_j!Wkpi57iahnNXwu%tn{rgnz4VcxnwP&h>UN#lq{WxxVRsM5Lv8hnMjA#zl| zv*4O)acFXk(iKZ#ixG2n`c9vwVE*{igBS*N#c*)#_bmtaEq|?fojSif20tXC#^v;r zm94)@!?G=P3_}indwI*$VlLZm=1%17OX zMMSMFm_9D+EaPVMDC=Zr-j=P9Ph|g8L^_vkIhT{SkkD%I{3o(?89OlRKy8l{uOE)z z9&S*^-ercu92xsG0&ci&<&|AVUSB^tvt1)9F;SH*--xWqVdMxe9@SYX)$%Pc(R z&zd4CI&R|((A>w{nWAU~YR)qs$nWxkFi4?VhdQw&VX6o!K5M@gr}KY)(VuUB337W$ zc)#MS$!FMi@z7uwpDSsc1Pq2RP!k*l3tYCvRNJlVj46|c?$ifkh&)HxDgo|?`tlD13 zSNN!$QpLleA3sB-YuLt!b&?k%=Pe<0+j>W7Us_n0e6W)ZqV3)MR4PKM^$pDD+HU%mN6Ac6 z%y}}k$sEcMd?DT*>T#wm87@_v8FHqGa@xa3g`2i}9gK%Ysy&H|!y;SJEX%PuWpMD* zHo~g%Op2&W4foKu_7gEkK@)sxDSB_ZKGas|I(#qQPNPlt6IQ=G>yBgDis;&Lq`VcM zcxsn}gNbI&WCPJWD<2-=0Am=-Occ?}-cz==*52LS6rujqLcV%_%Fpns3W<0K2Ufha6zi!z|b4( z`Zpaap@>8GRUtqdV+`!v7&LUA}@$dB1dX8%DJ20E2eoyr-A+fjy?`G;OhyiB*oAiWhNU&qqQ(4yX?4c1YK1-dgNUE#%u zcVCUJsA5Z+k7_|8q1&lEOZ!=nRCTa=;Cg3;phcwT&{~DZ57pWFEH&B}$jVtC^=(rv z*Pc&T9}AAB-FN)u-*tCjgA{Jj7Xu1e!e=dyCAPl%282ZI&5*9T0Mj;EUf;T;S>A|H z@J2+Yf8`cdrZm>VjJLo@>tB4Uc$JZAKDXdVLrY=YQ4X#Yv1`Gv(xg8!2XEBHC5nIb z(+0l_jO9YY6+Ls_6k#pO!2CTv(IKuhu)!P8MsNhZGx0^5GA0GgK;+f}Wg_x*hZL%g z>eLPxRqS+hj;xZ6u>#&SY?HPb9_PAnZ(e3r!g}RdCw%RGX**<@P{ZlqRP`z-Z)0zb z0bhXCGTY(yq4t(@_T0WVOc8qR-A77i)0yIHh6!3~dyrwra69aSz$=m$Y3up~B6WFn zY*+;rAydxLH)AWxCydvi@;k4+VDn2x*?(|!!!H)rn1}~4M7#>{Ug>V$+SXPcyD73) zBBn%2%eL+Gg|s*6Zy(hA3f03Es@0a}?fOOp?!!(j)@wUC-EDpD-6C z`=c&y!Q2+YpY~r@a((@FJguK*NZ=y<*^Ei-55`~_w;@4KGg_Iq#sJ*hGGYIy-`T~W zxml~YTF|f|Z`^`x$E-@~ENE|TP;QoL9cNThLE<7zUn4+p@*||nK{`2yytN#goE)R> z_p^@ig`DGST|d2hmd*|jN6RA3HMRm1-IuRIm%1XZelhvpOXprJUBt3XY&5EDFofb; z7ac`Z)eF2|s=-}@rCoHS-!CWnZbY%N1T)37%ItW|6R;WyL>qD!)2(hv|jSV&OJ#Rey|7E>sI11H}6iN#1-%D2?T5bD%Y^Z=3SoppCvt1 zo&Pf9DLuQm!Y=Orye+r2-wRW$CHe3&*Pm{=NrWF?d3q}blLqUSJwDTN?uKuLnIzx~ z+yracr1-k&MObhc3X&da$}~C*B|VCc#SpjHcix-QW0G`KGjP~^4K1?)#qr<0^oShu ztG*siCTeh7LLcG7G8J=tV5u6Jv>k6<@C5<|)T4%8eV%sIN8X?qTN}gzU!7(Xdx%qU zfLGeE(ou(ns6(2EX!CcNGE^0TpFI!)o+f=tJ z&@HV`rtoLsJ^Jx+BWWrigG&WPU=zHKpt>juAp@^!T^&7s->(9K-Hu7)JB#6MRD5?= zO0tInX;Rhu@-R_#t{O1q3tps!PR+Lh>&k#}KNI_y zdJp7?Q*#=86l79}dA7azrE}8=_~N`Xz)a!~>o}V6Gjtd4p5E-w6GnLXK21hXu4KAo zk5iq-&Y+7p^+O+nw@w^|hXOZ?ywNuDXx^3X54E#Tf`3I%Z`3SVX%LN@LtP4XLXN)( z3gVsF?c7wA0u~)AXJNbM66X%a;JyV3+`$cjY=6moDsv$H#X|)TanqU(m!fZNbH!`F znXQwrmUhl=8Z6?9Q3@}n(bb>VUL|VDjgcb)U$-Q44e&cPuUe&1Ikz!!Ck~YG* zvz#elIZaJcWR&8~hH1g?Y{@Zw_%jDP?hjoCaWZ~qY*>n8@{nC9_1SYCXVJG<@N?F% z7B!P3r@~M0tfpJHu~6a^vI^a#u>j)2iYJ)x?d(SP(X|f8ey&DDz4~aK0u}<~x>*yW zaE=3Q(3RaAK06AM+t-UAYH*hGUL83C4WGu_bAGBjm|P$hcexQ)eiN)v=jZJ?{Ynbg z7Xo2uLguq^pm?P8C2e_{Z#zp{z{kAp2ayr@2Ldu|U)+q(*`2XusU|Rsi@jV+(#0Zrf3~7~_xHsW(7b(=$}EIn zq_4jskiHFg(|&jRyRiq)_?eJ~?Z6K+D^rLMV*9hWRV9G^tH%?J%$reduM^htkbizB zTp;e=RgShoKHjaI*1*pRqlrjAT8N|n_+hd=ZyVE$r5AEb4@FJKemiX}_8vWE3RGFh zYl=M?w^2f`i(G;aApCQN9ckEaiIAz>_ZqojInP)JJ@pgyT`&S_-wKC_s_)M z2U6eUg%PE!7jgUtQX54r1}FnLabDI7DP|Ra;C@9w26K}leHg=k58nRw1nzwpQLG$< zwUOfY`Uy;neH{e~K+B$AxxctY0Dury1cIj+jA?}cdmr;>C_FRRxneL!Wy7zCF(zX< zVa3EsNZYXjscU0El>o^`xl_Q z)1~#O@}1DrhnV2|W9^0_`TaqJ^ryFM)Uxp>Oo3qFTKn6aKI&_X23BAD5ba@9S){X# zjDATd$)KmblgVZxr1jEm2*1U5sJy)4v(l@zU<*g{)Ne3LZ+@}^R#$3pE;xt=?OcIr zy;}3T8YuG+WX=(Zr&EWnlO?C$IP0AhWBEy19`oO_nQKDgMe?gGmZFz`0v1iDu8v>FNu zl8_f^pR*dcl>Z6)Lz_cl;3|AO7nU$z^UKO=oZ>q&=1Z5roX^khwec%@o||@F)qqq# zAppc3cIaDb^_(6$84GyFJOlB8j#Z)!W$!L7px1BVhUmFX7`qM@HVD@O)mKN{*sQtt z$4MsY*r0ufV}rZ<+BY+iSfAiTH7EEf#lB=+Sq4hvc*vpii(l&rXe}ftBX(U;w#EC& zieCHN6MizLsQ%LVaFOXC(DjAkfFK%G!vJN1b|;KPC3Cea&JyNkLyA!EbHXM_P*CWwDCD8J?+;0j&hRKTBmoAv1q?)DNgjM? zY|`9Z=D$Nw_ik{NP%Y%S-O1j%ZKJt|kKJOMuwb@b7Z@;E=2$7^#B883Cy=vLk{e2Q z^p!_1_>}h3wX&?a^4(RxcbL+a>1l7=x-V;~?dOae?PkF~!DT5c_D_|t8eS!=54ux~a%oM_$_CTb2_D-qYl(s@Z&VI}*=6vzviU&Tgo~UfD2F@pe~@-D?ZSXItI_ht z7?5k%or~aO8A$oggs|(P5M02U#xivt79)6N+OGSlAQh)&acCIG%W&K=s(Y^p!Dj8H4 z0TQ>I9x_Hl>^gILToRV&Ik@1v*Z#Y^xk?U1xSFUdj+4lBb734Rc9exI+yUWI* z{~&u$$5PN@$SOZ4WKA!ui_Vq(h_;O&xB?JUyNdTzIaOrM7y`(6ApX{)Jwhe+_=q}hkWyFIwLpL5eoW!;Z{)do_F1E zY`xOkV`Eks5z&39iy%i-Kf)dwUS7Xxtz&h%X-IN!2_!R5=B(x|-8po-GhE9EzLEQ} zuDKuLRHw572oR1hTGiT&Tdjn`S_x(}5vab?v{(8M1FH3b(SO z{Ji(-DgBx6cDKp%;$NwT%M9_*({|e?_idnx)+$C0gqBkF5dCzW1P>blI-zUC!S?6S zwLA5~zv`^}KoexDM(1%3qr%rVdAiK2KeRBvogSLPBhu@Dq~}k&sM_65N*jcW|IhzL S{^+kj?wKkQVic0^-TQw|G@if! literal 0 HcmV?d00001 diff --git a/py/picca/tests/test_2_pk1d.py b/py/picca/tests/test_2_pk1d.py index fd55ea201..60f3744d0 100644 --- a/py/picca/tests/test_2_pk1d.py +++ b/py/picca/tests/test_2_pk1d.py @@ -103,5 +103,36 @@ def test_meanPk1D(self): self.compare_fits(path1, path2, "picca_Pk1D_postprocess.py") + def test_meanPk1D_covariance(self): + """ + Runs a covariance test of Pk1d postprocessing + """ + import picca.bin.picca_Pk1D_postprocess as picca_Pk1D_postprocess + + self._test = True + userprint("\n") + ### Send + #- picca_Pk1D_postprocess.py takes all Pk1D*.fits.gz files in in-dir + # => copy a single file + shutil.copy(self._masterFiles + "/test_Pk1D/Pk1D.fits.gz", + self._branchFiles + "/Products/meanPk1D") + print(os.listdir(self._branchFiles + "/Products/meanPk1D")) + cmd = "picca_Pk1D_postprocess.py " + cmd += " --in-dir " + self._branchFiles + "/Products/meanPk1D" + cmd += " --output-file " + self._branchFiles + "/Products/meanPk1D/meanPk1D_covariance.fits.gz" + #- small sample => k,z-bins changed wrt default ones + cmd += " --zedge-min 2.1 --zedge-max 3.1 --zedge-bin 0.2" + cmd += " --kedge-min 0.015 --kedge-max 0.035 --kedge-bin 0.005" + cmd += " --covariance --bootstrap --nbootstrap 20" + picca_Pk1D_postprocess.main(cmd.split()[1:]) + + ### Test + if self._test: + path1 = self._masterFiles + "/test_Pk1D/meanPk1D_covariance.fits.gz" + path2 = self._branchFiles + "/Products/meanPk1D/meanPk1D_covariance.fits.gz" + self.compare_fits(path1, path2, "picca_Pk1D_postprocess.py") + + + if __name__ == '__main__': unittest.main() From a6b56580bf4465a09a9fd35da0e5dd92ebca5a64 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 15 Jul 2024 07:59:04 -0700 Subject: [PATCH 11/11] removing bootstrap from test, as it uses random --- .../test_Pk1D/meanPk1D_covariance.fits.gz | Bin 9172 -> 8395 bytes py/picca/tests/test_2_pk1d.py | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/py/picca/tests/data/test_Pk1D/meanPk1D_covariance.fits.gz b/py/picca/tests/data/test_Pk1D/meanPk1D_covariance.fits.gz index 745f1d3579050e07a79cf088444691d6afc6a53d..eb0595e03835e44554686546935105b0b514b518 100644 GIT binary patch literal 8395 zcmZ8`cRbsD*tJdVqE$0iYp-fiBQdIq+OxH5wrZ5xdn+ZjQUo1}qH52WjUBVL*n7kd zk>sWKeLtV)eez$%b*^*Hb$-b=zY7pYhIjoSBHYdKHlR11_`PZ2PspH1rmmkf|9CN* zPsPs=|5gkwEq_}?u>~=MiY_g28g(ocg~093Wc*wHs{CQ(1u+E^BxJW+IITI_)Y{Dj zothDTh-}bDCi}{DDLBv}H?zmqt>TwQ)k9v|XJpa?X?bc$3NkMG0DfEMwrIK~)qE2A zc0M8sL1RV&;SUT;8U+`>FpFk_#;Sr+mjTlxjxCAE1v^`T&?Qxq+q6%(BMOdMzKGFr zKH!LD?M#oYu+n~4>8`l$c3w4kC_a){pra=;Fgq|lJO4qsB^P72^FmvI`Oyd3XR&)z zyJC6@WTC9Lo>1o#>nKeg!aaUgO4NF*)OMqkAI5RG8?@dxxL;zG49Z+kzZ*HE#6Kn$ zrtv`d1FKy85H~0P2U-PgR5FG@sk?>fL9Ia2?I&L&bizKPkMfu3?m(%{*)=JU`9;O_ zTzra)ZMI=$`JpnLRlLn**dJ5h6Z!)9T5N5?64KKpbtAv$=Z8MY106dI@QQh*$FYYQ z5)c!SlQGcpap)-Ig(%>M65zB#b4S$}1xfws!?w%>{f(E*Og#+KDYG`Fo|C<2)+FOT z$!0jsQT`S26Yee>eD6;Fdf^j^;UoL@UmIe3>-8o^W$SPmM_^873w~Qay{R6ptZ;g` zXK#05yg&9w`Pb;RI7%dy>AZ&@eAec8C0$UE8$#^@Y!BHc7Py=I$KzN3vwrff`}@*x zYeo2`G;EvrMXP{24ZCpa4t>a;_)&4F`n?Js>&)tvbZwS7OLa0o-!TvTW$mvK(R5?q zSJD>o4UVgTiL7jD`4?I?YEgp?FHACLvg#KOnKlZ}FUrA}1l@`)ALQ?)1fAFTgguvuFMK-j}YKRzSB@wm3X1sZ&q zKb*}H4$cGzpNu;olZp(;ZxOF_ML+Ryp#<6ku?I|IClprQuoRm3wb?rX{Q;ND61Im) z1%ADK|NEt_-QX5{dHF=+z$+owWMH6{=gk3MXkiMX0{p)qwEL*A|ABnvbi!(RXW0gi zgBrlx|5Y6I@f_7%IM>0*ii{4}I5-9`$DMLqh9`l+0Wy?8y)nEKzl5{k!-G>A=O;km z;pTl*xBPNA+f6XIHrP?(@E}PW+yl#_$R`@=o+N8qUfu;ve+URz8Qov(_)RLv!47c6 z4cV_Ob3^czLBCN}Upw}Ln=pj(6)I;L;D9Jf;BVBz0BiyG?h=R)^&=(yo z$qfLU+l!yyS?-W(h10OG*TER)vH1htFc6SU5@$u+IZQ*fHgX5s)(=w!-o{-)vjJqT zu&EZeQA5@TH1h#XnD9(6_@V;+A0{;fqg)6F1@IY&3S>TMBO$w{XSKHcIS7Ll$D_U0 zw1J!gy%WRTQHTZsebA&50KiGdrH#%NmU1*4o8Pen*opBuS=|Nn2Lz0PyFfk`&`sni z&4Y+?7|7fRx8RA1BJp@dd5sDQYQSJaU6^HSj&&Bs7_Z2F3Z7-qOL<1k8cxXpuDl;RX^$`Y-Q~ z!IZJrtPihPKfk=1OL@z`s!z^UMcITSdRaMruE}Ryk)PNH5aPW>4OxlAzYVazCg^un zpLgY(qgcl8mF3Fd4Qs48A8`J!rx=s}@RY&V#p?k3l{~Rz>NR%MfBOHSws}SW7_|fV z2W;X$Vf`SIYhi!i5yG}5#28!t3w&`&*o98vByl9G%)9;G&B2|aNu61Cng);f zAX@6Ju#Pf4ZhG#01x8g^bEl{F0>i#z4)L_s&gLW0w#*AVx0e@8d@?)k4Ma$M4_7VC zJRPXw*$nwW>KAHh{()Tz{#v@B+Me=;AZs}?ryVW1&%-wo#!JQOE4L);(Av{#bAlSJ zUzMl*581(C?s-U~9$Mq?08w-uw-+DXy?X`b2d$L%I7OqXdbO1eEG)w29TF55qdZrG zQsfj3%!=PWJLTisk+yG(@5SR=nvyRE7Th4mab&yP@g%2yuD_SD#(bkpK6AQ>LOe=I zxjE_yz%?o*JxD*ml9op0b?+GFnPoVkL?GQa0Xe+xNFwcyQEyt#%$u;>z(-WP=?Gty zprD+ffOE0An}}l45&%#0Q&x+w>;&;m6F-u&a&P3;HxxC=sE5jeWPB+w^JHw}rRLN< z0&Yg>fmG2ArlIE0AfJJQ9@1X4B_G#Ogu)hcy}vjGW7sV{f4AIb+vB4_aexun#kyJ0w$BNe6)rp~{1&!ZOC(PB8Lp%VxL7yPMS(m}1> z(5nGc@jXdq^VVT(>nex(`acaYiqfYrF#R` zGi8d4-+C$Iza4SBXgcsdzW;F0Gb5SYP4o_aS(o}PP=Cpyfw$k+nIhg|#*MtKw+1J( zst;odFfOa<3b&+pTc(}jx1~0`skFO`9+=V?f$rEyG^5uY1RX7xr^=xHM8zzthL?9B zH_nWvZqovvMJ)1C5&W7}_4voenNXoA2^kqDiO|^1anEvxJ%Y_lb*_k3qePE;#_PMK zFb*fP#)69+dYzM{i?RgRa9-dzGC6}v%gQ@nW8JZke@Nhbh96vEn%{mo9rL-jYsYl8 ziRgw+?;mz=?hx%F`rBjk-1QipmViJdu8^Peb>RcXY2EP2;Vx#wrwEspl7Z4hm@Kjl8Gh)DhjPQ_hMNt zBkKz5le%iRA37E)ZZGJUyrK=1-#1eS1^`ro-;^ufGeIy3R{0HGv?;HGY|<$o(JVLZ zaLLw`h6()A3$=ti?|8xUrb)_=_30FmuZXn>rjo;LLv&f13600pLI5A*UsmWf{FM0{pon7 zG_{4-u$8F3KNqaqdr9AEbp`GbY6b^~<*kX@HX>0}cTEMYFP#eq8Mz~98C`LEi|S0_ z<21l?$?z)eKeulG&K7{W!zulo1ctqBW$08A%D621Qo>eahehLDi~a_YVcf229zlMl zWuwqL!QwGJ7KE=<=%yAnteAMb@*4dv&IW_UYLyXlqLm&s!zqlkR1};hbHjy0g{@R) z^|sP_Sn^K>QiTtIkvavQ(VxK`_@|}r7Gd|jAaR2rIif6H@VjircX&IGc5vy3A`D92 zuPQcdW%e$xY;0;rNaBwCq@R?vsGEO*gTJY*dlAq%MVn9uaGu z<+7UP$`WeL?x62!Y_UEqLU?CZ1wE2>vzeXDUJ^d8dAkpRJIxSvJ&EYEBp)MW$kI|+ z?hd@Nd|Km0@9Sr7;jshMVCXE4yzG^Wm|#mi<8QWSe1@g+#8d%Zv#I`{yPiKkeaO1i zLcNnd@u{@xE+s8+&wzPm`fMD~JEf>LQ0#V~Q*-T&x!@`#{ht&YdxUlrf zQ^qGhKH}fc2;6(4>qovjXWYj;;oF+ez0)aoOgMENUG-hx@4*ub_kDJKdQha|1V12F zTyjoTblIXT_V*r(N6<~P3>0QChlQnfK2oCTYfZ+e(&-6qDTWy>4eB38<*m$gO!o5g zw#}y#>8dtw#^^2Y@--ZJNTWsE?G2oGk`M2-9Q5nieGD>rlbJNlCLIKQ-=)FD7Aq61 z{p`eq!m$_KncC&xQPw!CJTRDMxNz{bSU82HKzckUXl-uaq045ZwRC5_mp6pvzFqf# z>Zv&(uedztv11oaVr`-U<#9Yo34k*Fjc5miDsX$3l zJQ4DoJ;&pU_{v{9PSRr-*F-xxsd~OeISzyVi4k-e!N@jiAqZDdjzCo3ikx9P%x6$6 zpD_(&^2r&xXI!<(ES*aGGTYOJ=s}~w?bZ7lcVhKGCgWMeKV8#~=u*|B?GGGJ%#qgS zxxN9aA|}v+zUp@!uN}|!aw}g>Qx7$#&HDD7bHxn&+%n@q6St+rFEjoA6xSeGr3kiegxzJP~d;Q~$biaw9n?Rr1H~S`~=iPG7 z{^*J)_BvzZLX&=nTnZ8&hW-E_g?h{_OEH=bR^x#ip&Pt%{3{0Dk2lQK^Q$%=A<>5# zyC=Q(V@7W_h^pNqyk){k&*P?Nt#tTwN7gIj&!usg+;+UQp0O z61{BgIg?Aie9K+zC8qj~Z6Ehw@2~PjowI6={e_-@`abgzRK!xU*RVr zkZrty^E!#feA(~RF}JhbYL@(_f4)QRIhIvUTfX8Ds@7XK@6m7s1v?)qMSBfl0%U#H zogeb6x?Tq9_>P@A&+=Ie;~f_l4d7kP?vh<$y7>E@P}+ZR8Nw8eG~q|eAA#PV`W z!FM@DczI1$Inn{DE&u{V>mAk-Q1DQ5wS0Tm*B^I7tm@mKws;GERljq${6!a0Y|iv+ z+=Q*Tg8Y;@BJG4cmMSwBYB5Y6Ti&1R<~rcR=M_-r*3iV8+KUQ;9%(%CLH-gFcD*p* zI}BsqufQropxU&CPPth>*;v0PKG!M4F|RP^jLtfwY2%{VHPUr1iDc$_G0on@1jR1> zF}3e&uf71P3YrHp7#PFpJq2bmo)(4)4x*{e^BdcrMpnO!0^}+vmm}maD0Df0F5@{I8pc!KWi+{Oh z&rCLAg<^`w1ozeaM-@p;UDn`YgF9W>r*XH=~U_$~` zfsbKqUZjUwY7WXT^_#ax2N?JaW|EssHvcBI?X9~qO}x+x!mJn@m#Td|;-06o?>|Ew zl%q!?#u>HDqLs?No^uIAw~iX&>J@eJSp=a(scXJr>@dKwbvgHV$=Bc!GjtD%vHJ^Iq5^Q5jWqfZ% zG_+dd-vH9<;>v#dAnsq?^pqUFXEYs}J--|PJVH?)#(EgViWW!eiE3-+eiNr}x+yFW zAe&8OqyeWbHT_BZ!Y@J6SX>}PDaV>7Tc-1$AMoc)(eux82vw0AJP;xv*N(}E#%ItxrMAG`&Q4X zo{u@}N=dt)KD|B-_Q%)1J!TrDpZ&h7YX3v6Q9J4Df|vJ_4>eF!h+WT>=J`Tg`wM-OEB z?8Tzdi^Evd^Pb}xNmFuR9@n22ihgTi+U~24b^ZtWrl;Jy3o#h=pjOXF+c-M=^#6?(0v(n18 zqIx-0@P%Itm6I7w)_hr!7H34r=SZ|sAnq4^)|9Nh{{DqhreS2tnYnal{m70}-eu#P zxCtj{_?oqoW{D_Y$iTiS_M#?Gq)J`lG;?j38YX|gStwbbx=j1EMS}|^{^RWIN2B1p|}f?Ls2DWz{uD?F4G zIL|g(N0-EUQTFV}M=I3~n`N_(V~WTCVN|^p`5E9_lgG-(_^buFe!C|PT71^>DjT|q zlB{+ca~}a2G%hlXmKf`*GjM|j86cmAMtNgXa9Y$^>^c6gt*i%oD+yn`iQ>x!rLw4} z(_ehc)sRz0cnFKUSTcLa`Xy+n=Ffti-l8a9CjZfD>-K#laslxL##~~Bs`R5;HTIXnU%c7UvBeJVmLDYW(ae0H${nls z(7TRror|V(jmQfU-jk4OC4~>ZK4v&4bCq_%B&dK~a`0gqmC5$qGX#eyHN%Gzz%VKp zugE;g>0w>LUGMeXA#%5sD5vdOlN;Yb``R>_}p&X;m~R-;TlX)qN56EpT?GO|9a0d4j#b z3wsP|9Nx(#DE-bj?Lik5(t~z3;!S;qZ~9}+`BgR@aX0VT9Xs%J?3WUDuy2Vg*;)T1 z*X2u@o98dZjQ6&d;ExyiS-zQ)n5jq#%r_do`^K&=jR`7DwZP+wYJ3c#9bPdzBc2$s zNGlFE^;%w=cf<@7?Zp@oczTp6l8iNM8fyxoNvC$^1PsJ4iDoevyY3|l<(u_zMcs>} zxT}*-6W7*RF4j}3AF=M!!<4VwQ(y(p&vNdxPHEBjq&u|jr+#$c&_!P#D`>*66QWH{ zHg7qBBcr*f=4okxCJY?`o3Ll~my#a&=LEjCzUhF%{BHtK~NA~{DL6h;<0=d%_>W3P!j~|u(pgiD%_wy;e#}1 zzz>m$3e4qLshIWs^FoA@#|-i5#?0wkgu|=h;2q|j?BWo;n~;Bhr4uHi;57&DQ4 z1ADzB3&qh9=g}S=>GQQW!O4@=v*%%b86lUs5UtY0RykeRpN@H1`{2$KU{Yjh^iuzN zF)n6wV_>3@o!S}*IZ@4)gow}*z4||Ex9dF|QWw7svyP?UJgycCXvV!kLcMUA7@P)n z8x0S`LBdXwsHdreTGUMA}EhxDM^=gjz z(WszY3+kPHa+3Equv=C$?pw>sNL4jI7#3qC6;2aZ=$(^ zpg3T&Z~ajv((1tM1R@rHCtcq`KG>i1cf97o8RIIy!B?|)juvrU`EX+QQ7F`+`K=rZ zebTm6e|TLod$P|X2^J;BFHgVjg!7NnAxl?}>CW+$Az1$X4RydlGT+p2X;|2$DMVc& z5!CdG^3uSYBdDL{*dN)vOtB2Zw8+tNz4m^9Kpk{JMcGFIpcaUO6%Gbwi~?|%|1gBg z-ns%}MX9{!>Nvaj=xQF{WsM0_4q#3CE}^m(S!iY{XB>1p@F#Y*gg*KbBWEJ_C)n+5 z!j^A%_v2WhYzE`YYX_h`PnHBaChaP7KiRyUMbl_B|HSs>KrI16xw63^n%nSd+{f6P z8io5=*5`iQ*V>^I%R_kDWBZvt*qpbP!5yf#(EV4z#OW~_YfY+%!q5|k4*;ooEoLLw zxEw;{utt=983smZr!tA(Po#Ug6$iVZagn2dw_?F8uHJWQgv#Y*-3weSaVx#9f=UeO z+X1-Ej-nQyvdfyD&gVTnHi63R=xb9q!cQE^F8+?=yfpCrA9^{Z^?bOd&4Jl@DJcDX zPQYx%ccMh}igxh-)sTj`m*vwdvZ6H}QDFjweM7;U{Y~2fUvUElT%7l|UGISQ&R)tD zJxyOK>HT7fS(ms2a%q~-ofNybW12lW+j4pJ{4)7j0j1|{6p7sUg+(W>69Mnzp~3u)b$XWJ96hKO#myWypdku6#{H zI8hGgoiCD4{aI+i_NN8EzcZZvZP>upqwzX?f9Ft*`?(R`<+19pcpa-o@ZCMO+HAR3 l=C}XsSMUQ|yD0Q-{O{%XKil$!qE;OK7}u`{3dwkQ{}1)T7YYCX literal 9172 zcmY+Kby!qw*Y*YJlx_wPk&;e{0YL$g5Ts)ON$KtxX_Q7KrA0ubk?tBgq$G##8XBfP zyzl3IpZA-;_HoU=&b8KWoqNr3>}voq1o!@aaPA)l`0|4{3DB!6 zM34-m!E@;>AbZ=tZM1PL@OV*LI>ThaD1FWCg%n;dVNv`UF$vdNPo?qT;_GM*CfnoA(7J7ydVih2Ql7q0} zq{Glkh=tl8ig*0(C>m_luN6j?K%x}=JUw3Lm08v+=T~Vfq!+1+zu1cjtT+Y2A4pCe zsy?&ej;Ci*BZ(s*Auj9&SvxvAzgrKlQzd&NLCn>M7f_I#9`ykELFzFjyBLRTEcFM^ z*t?t}hx-g(pOFMnptL^u(>?t*!NASXdv>QS4CVeG82t zBaSCxU#q>ic;mLJbTZ#`wpH*n01y zRfsTdA2q$M#s6Ha0Cg>PJbq6s&h!x!)s8Rr5w?#HzWIQs$kO_h_)+4TTT*ylOxLQT zY6kB)K&2Vn$B;W)g!f~|TXa-Lz~g-1l*D^#-|h*umin}$P`HPsI>C3l7|U^&0ly5A ziX3loFYT5z1z7+4yrICC2@^33uN#x$we)MN=mrp5M6WW58Ncp(2epZyIuq-CSug3) zb1*xHo49+lxw6#6;|GVro8=+Ooa67ugG^-^v$^JBmD=Z*_B(ieIn=Z?@w!|T2G4Pm z$A@m&tW00g5BN%u(%()D3@=U%_9U?b<(#QqH_-ipiS!5RN5V-fRPU9VwaA|@dwLk# z+Pxf&`u)zUm=py-T@!Qp66ucz{5OvzIFG5V^^@f^H#)gLMT)k{L~uEuKdjY!MsQ?1hW-2V{6`z z0DkcIZ@75*`6-CC;90sdTP01A~cx&knQ$ zj86Z1muPN05Y~}L(a+od-#}-@NP;{c#ca-J!vsn%w(0(>pb?jMM(htpDkfq9YfriV zNlX!(p@e*fjY%e#%RzgCC#<`H{eH|FMbEcXqESk%d1ojX(hw2F#gGt&tQ|l`E1iX9 zb36q9L@9lWg(`xdQiH91@+?Crl$|1l)}Zr1v4Ow4WqQan0kt>`@VAG63~O1k)5(YX zD%pVA1NsEu3_8U8Pd;js4Hzr}JVjcJB>aNdC&Q)s$KDU$?GmGA*i8M;iYsdbz#n{# zcY9+w8Mb;8@$WKJD`pgEPbe}@UMBmJ9f6N~mzq5?EQO`}Ca-rMIL^zQ1uzyvcV++_ zF}W-5+Q0xzZeBlk`(4g4-RL~=1Ax3Grxv9!2T3Md>~+;8IT0)V1rt00NFegJf?BBq zfN>fG*%=Qouj)&a%%wl$v=%Y(YRsA{4Xjz)L`aI@37aVc?sRA-@yyfXD!se zcK-vsTd;`vA9#$I+@gQsMcsuKJpgd^tf+?4Zn4Im4b20&{g|g*dw^)Re$~6i=5$H=Y8z@@(RUHYS1Tbq@z+oQ4#SBO zUxH5mwjB}c_W5M*Urmy2d}II5&sL##L3)v4Pwq?@77A5F-Sz%IHJgzkqK+ydiF0vi zB*UFK@-|ZMOp1DOw>vcU{4S%WJNRb9K(qza4_x^#qwar|lv3}iQ(p(1|7*X$eSrbc z!l$H9|H1zE+p8i736*unKdu`C{gYS0Rf6UJmsi0RIqjWO{(ABjwF#hq{^G{w5U>;D z^)&$e^b(Q&8gKtk9_EZl;piZ*IyI)vFf_J;uMXCBa*sPGXQKgP&aGOSRvN>}89LnH z;Oi;2pzy-CZj@)}P4RqDUDbgJ<#p?r#RFeMCT4F+nwN2n#4Lvgm+P~@5Yg7P+OOh= zyQ}HPih4`D+)2T-r}6E{PwoYEbeYs7cdLR{t7WQLEFFT~0#f$vy4s#>F?2(D;23T+HjyOz{skCGfL3l6Uz=EF#-Hlu84!^?>Af`XlZ)gBk;@}XFrY( zV2yn>r2(;!hsgvUv^Id9Lw$H@<5j*{r}L>p!qt{XI?HpyvWbI;Lb^KXJo-5V5*zhz z!7VYw{XQopGW~M62ie*WEL1uaX<8g>r_P;D0uhV(*~IV__{V%YI>&gWu_hQUl{*f1 zj^(>rHw6Bz{vGizEMbq9_OIEmhIGSa0qWY3-xdzpoMm!+nE67La- z*x+o@(B1bfS$NzSzldWps4gL^zk_Xax%m5WwND7LVeQOc!?HiWOUQKIN8Mgntu3gI zLs>-sn&1t-N>}jNhL%%X4NFMYe0`s5Gg;%p{^gHfpnL8}=k~RhFX`K6Z^Yb;CY0Ue z6l)JG6(82L7|b-L$iHqfB-{f)CHRA;qWk9#^mmAK!fOEB?fT8`HK*BV9*15sVxNmJ zc#0<&b8SuL%qZ8k8^EMU6Be4Gn8e)BzWjsW${hPe)(r8m8P~9>(e1-mr!J%(g`@8r z%7(B$_RY?;>o4I}dv5rL4&=*57ZMbYrd&OEeYF$|)tT^aS)n1WrH?(k)t-HLRmfssj zS;Tmm8EyI3Piaw$PU6hUg`LNYo`?*|bQ!#@b03-Vb|W%S)KJxzSqn%Vr%!Rae|s(e zasOr}I@)A0tLs;!GioJi4*)+}7wAyd+D~h9F3rpTY(JiAM_Wnwfg1j7?7opF&IRnq zfABvL35m#pCQ%-pG~CFPr##*DbQ`*3XAvwRWv*hc6cszMNfBS++TirPHX359j&03P z)tjZY3hwjJm%(vY<2u3RRh6$uKT8RNMx37s1(yzAtt^_r6{97{ua^!B$!yBs_c%Z9 zyV88vS*QbrjU^cp+6+lFIkbPm+G6-XD6O^M@fj}7gpu|j9d6oPUHn8thP^VZA9XwRr~*$_LLACjYSK9nYeQrL zY{UEI%$40F=<@!D^Hq-Q95pf2hB(^6H6~0X-qv<$CL!_5j6~RRID706*&}(`WYQ*m zn)&{Z!67a2MPD~u{X(J{2z?Cjg4wlx!2M;G=@PHRSF>Q0lp7cH1dAImvR>M$i!Lms zlk3>VOD&LP|C4$+2+haf7C#%CxEJi1*Fy`FY|P(Qd@&Tl7fy<`GCoqGdpe}U8(7p9 z*jDT-x{*o$$WvMsKm55pZXF?ug#plB)Vb<&Mg!mH{mNkEO(&h3S#W{cTp3CT?&cAp zqsudOHs`K>8;9p=t~&`W_mrS}T>h2VqXl^3I{=K37@ARKkZ4xy{4f*b(xiA^CISn~ z4`kh2o@~Z8Uoewk?okONU3;;37QVSkS)xi+E)>k2+c@B*86nR-sA|@6QAyMQkWKk& zG>M%ouH&0BxD$3-K~A-2n~jHPPWC_a+G0C(_O;?;$?m?nO04L}^}=J+oc0)AXo{J>8ToGlu*kF!a+8mw`J z$nj-<#dX=1GmyeHO-f&Pk*j9#BmWXqhTG+crx-$kR;cWHs19Ab^8H2YEsOsC1qnUd)^ zzv~Bs$X{~h2nE(|vHG>Qb%>R?0{&Ebbd*@FGuopr+rms|$W8+JsxZ~vj}j|qTOcJ< zcJ~a_j!Wkpi57iahnNXwu%tn{rgnz4VcxnwP&h>UN#lq{WxxVRsM5Lv8hnMjA#zl| zv*4O)acFXk(iKZ#ixG2n`c9vwVE*{igBS*N#c*)#_bmtaEq|?fojSif20tXC#^v;r zm94)@!?G=P3_}indwI*$VlLZm=1%17OX zMMSMFm_9D+EaPVMDC=Zr-j=P9Ph|g8L^_vkIhT{SkkD%I{3o(?89OlRKy8l{uOE)z z9&S*^-ercu92xsG0&ci&<&|AVUSB^tvt1)9F;SH*--xWqVdMxe9@SYX)$%Pc(R z&zd4CI&R|((A>w{nWAU~YR)qs$nWxkFi4?VhdQw&VX6o!K5M@gr}KY)(VuUB337W$ zc)#MS$!FMi@z7uwpDSsc1Pq2RP!k*l3tYCvRNJlVj46|c?$ifkh&)HxDgo|?`tlD13 zSNN!$QpLleA3sB-YuLt!b&?k%=Pe<0+j>W7Us_n0e6W)ZqV3)MR4PKM^$pDD+HU%mN6Ac6 z%y}}k$sEcMd?DT*>T#wm87@_v8FHqGa@xa3g`2i}9gK%Ysy&H|!y;SJEX%PuWpMD* zHo~g%Op2&W4foKu_7gEkK@)sxDSB_ZKGas|I(#qQPNPlt6IQ=G>yBgDis;&Lq`VcM zcxsn}gNbI&WCPJWD<2-=0Am=-Occ?}-cz==*52LS6rujqLcV%_%Fpns3W<0K2Ufha6zi!z|b4( z`Zpaap@>8GRUtqdV+`!v7&LUA}@$dB1dX8%DJ20E2eoyr-A+fjy?`G;OhyiB*oAiWhNU&qqQ(4yX?4c1YK1-dgNUE#%u zcVCUJsA5Z+k7_|8q1&lEOZ!=nRCTa=;Cg3;phcwT&{~DZ57pWFEH&B}$jVtC^=(rv z*Pc&T9}AAB-FN)u-*tCjgA{Jj7Xu1e!e=dyCAPl%282ZI&5*9T0Mj;EUf;T;S>A|H z@J2+Yf8`cdrZm>VjJLo@>tB4Uc$JZAKDXdVLrY=YQ4X#Yv1`Gv(xg8!2XEBHC5nIb z(+0l_jO9YY6+Ls_6k#pO!2CTv(IKuhu)!P8MsNhZGx0^5GA0GgK;+f}Wg_x*hZL%g z>eLPxRqS+hj;xZ6u>#&SY?HPb9_PAnZ(e3r!g}RdCw%RGX**<@P{ZlqRP`z-Z)0zb z0bhXCGTY(yq4t(@_T0WVOc8qR-A77i)0yIHh6!3~dyrwra69aSz$=m$Y3up~B6WFn zY*+;rAydxLH)AWxCydvi@;k4+VDn2x*?(|!!!H)rn1}~4M7#>{Ug>V$+SXPcyD73) zBBn%2%eL+Gg|s*6Zy(hA3f03Es@0a}?fOOp?!!(j)@wUC-EDpD-6C z`=c&y!Q2+YpY~r@a((@FJguK*NZ=y<*^Ei-55`~_w;@4KGg_Iq#sJ*hGGYIy-`T~W zxml~YTF|f|Z`^`x$E-@~ENE|TP;QoL9cNThLE<7zUn4+p@*||nK{`2yytN#goE)R> z_p^@ig`DGST|d2hmd*|jN6RA3HMRm1-IuRIm%1XZelhvpOXprJUBt3XY&5EDFofb; z7ac`Z)eF2|s=-}@rCoHS-!CWnZbY%N1T)37%ItW|6R;WyL>qD!)2(hv|jSV&OJ#Rey|7E>sI11H}6iN#1-%D2?T5bD%Y^Z=3SoppCvt1 zo&Pf9DLuQm!Y=Orye+r2-wRW$CHe3&*Pm{=NrWF?d3q}blLqUSJwDTN?uKuLnIzx~ z+yracr1-k&MObhc3X&da$}~C*B|VCc#SpjHcix-QW0G`KGjP~^4K1?)#qr<0^oShu ztG*siCTeh7LLcG7G8J=tV5u6Jv>k6<@C5<|)T4%8eV%sIN8X?qTN}gzU!7(Xdx%qU zfLGeE(ou(ns6(2EX!CcNGE^0TpFI!)o+f=tJ z&@HV`rtoLsJ^Jx+BWWrigG&WPU=zHKpt>juAp@^!T^&7s->(9K-Hu7)JB#6MRD5?= zO0tInX;Rhu@-R_#t{O1q3tps!PR+Lh>&k#}KNI_y zdJp7?Q*#=86l79}dA7azrE}8=_~N`Xz)a!~>o}V6Gjtd4p5E-w6GnLXK21hXu4KAo zk5iq-&Y+7p^+O+nw@w^|hXOZ?ywNuDXx^3X54E#Tf`3I%Z`3SVX%LN@LtP4XLXN)( z3gVsF?c7wA0u~)AXJNbM66X%a;JyV3+`$cjY=6moDsv$H#X|)TanqU(m!fZNbH!`F znXQwrmUhl=8Z6?9Q3@}n(bb>VUL|VDjgcb)U$-Q44e&cPuUe&1Ikz!!Ck~YG* zvz#elIZaJcWR&8~hH1g?Y{@Zw_%jDP?hjoCaWZ~qY*>n8@{nC9_1SYCXVJG<@N?F% z7B!P3r@~M0tfpJHu~6a^vI^a#u>j)2iYJ)x?d(SP(X|f8ey&DDz4~aK0u}<~x>*yW zaE=3Q(3RaAK06AM+t-UAYH*hGUL83C4WGu_bAGBjm|P$hcexQ)eiN)v=jZJ?{Ynbg z7Xo2uLguq^pm?P8C2e_{Z#zp{z{kAp2ayr@2Ldu|U)+q(*`2XusU|Rsi@jV+(#0Zrf3~7~_xHsW(7b(=$}EIn zq_4jskiHFg(|&jRyRiq)_?eJ~?Z6K+D^rLMV*9hWRV9G^tH%?J%$reduM^htkbizB zTp;e=RgShoKHjaI*1*pRqlrjAT8N|n_+hd=ZyVE$r5AEb4@FJKemiX}_8vWE3RGFh zYl=M?w^2f`i(G;aApCQN9ckEaiIAz>_ZqojInP)JJ@pgyT`&S_-wKC_s_)M z2U6eUg%PE!7jgUtQX54r1}FnLabDI7DP|Ra;C@9w26K}leHg=k58nRw1nzwpQLG$< zwUOfY`Uy;neH{e~K+B$AxxctY0Dury1cIj+jA?}cdmr;>C_FRRxneL!Wy7zCF(zX< zVa3EsNZYXjscU0El>o^`xl_Q z)1~#O@}1DrhnV2|W9^0_`TaqJ^ryFM)Uxp>Oo3qFTKn6aKI&_X23BAD5ba@9S){X# zjDATd$)KmblgVZxr1jEm2*1U5sJy)4v(l@zU<*g{)Ne3LZ+@}^R#$3pE;xt=?OcIr zy;}3T8YuG+WX=(Zr&EWnlO?C$IP0AhWBEy19`oO_nQKDgMe?gGmZFz`0v1iDu8v>FNu zl8_f^pR*dcl>Z6)Lz_cl;3|AO7nU$z^UKO=oZ>q&=1Z5roX^khwec%@o||@F)qqq# zAppc3cIaDb^_(6$84GyFJOlB8j#Z)!W$!L7px1BVhUmFX7`qM@HVD@O)mKN{*sQtt z$4MsY*r0ufV}rZ<+BY+iSfAiTH7EEf#lB=+Sq4hvc*vpii(l&rXe}ftBX(U;w#EC& zieCHN6MizLsQ%LVaFOXC(DjAkfFK%G!vJN1b|;KPC3Cea&JyNkLyA!EbHXM_P*CWwDCD8J?+;0j&hRKTBmoAv1q?)DNgjM? zY|`9Z=D$Nw_ik{NP%Y%S-O1j%ZKJt|kKJOMuwb@b7Z@;E=2$7^#B883Cy=vLk{e2Q z^p!_1_>}h3wX&?a^4(RxcbL+a>1l7=x-V;~?dOae?PkF~!DT5c_D_|t8eS!=54ux~a%oM_$_CTb2_D-qYl(s@Z&VI}*=6vzviU&Tgo~UfD2F@pe~@-D?ZSXItI_ht z7?5k%or~aO8A$oggs|(P5M02U#xivt79)6N+OGSlAQh)&acCIG%W&K=s(Y^p!Dj8H4 z0TQ>I9x_Hl>^gILToRV&Ik@1v*Z#Y^xk?U1xSFUdj+4lBb734Rc9exI+yUWI* z{~&u$$5PN@$SOZ4WKA!ui_Vq(h_;O&xB?JUyNdTzIaOrM7y`(6ApX{)Jwhe+_=q}hkWyFIwLpL5eoW!;Z{)do_F1E zY`xOkV`Eks5z&39iy%i-Kf)dwUS7Xxtz&h%X-IN!2_!R5=B(x|-8po-GhE9EzLEQ} zuDKuLRHw572oR1hTGiT&Tdjn`S_x(}5vab?v{(8M1FH3b(SO z{Ji(-DgBx6cDKp%;$NwT%M9_*({|e?_idnx)+$C0gqBkF5dCzW1P>blI-zUC!S?6S zwLA5~zv`^}KoexDM(1%3qr%rVdAiK2KeRBvogSLPBhu@Dq~}k&sM_65N*jcW|IhzL S{^+kj?wKkQVic0^-TQw|G@if! diff --git a/py/picca/tests/test_2_pk1d.py b/py/picca/tests/test_2_pk1d.py index 60f3744d0..f4f509311 100644 --- a/py/picca/tests/test_2_pk1d.py +++ b/py/picca/tests/test_2_pk1d.py @@ -123,7 +123,7 @@ def test_meanPk1D_covariance(self): #- small sample => k,z-bins changed wrt default ones cmd += " --zedge-min 2.1 --zedge-max 3.1 --zedge-bin 0.2" cmd += " --kedge-min 0.015 --kedge-max 0.035 --kedge-bin 0.005" - cmd += " --covariance --bootstrap --nbootstrap 20" + cmd += " --covariance" picca_Pk1D_postprocess.main(cmd.split()[1:]) ### Test