From b64cedbcddba9d0e6255365ce2505d5285d2c99c Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Fri, 1 Sep 2023 16:19:31 -0700 Subject: [PATCH 01/13] fast computation of metaldmat using stack of delta --- bin/picca_fast_metal_dmat.py | 570 +++++++++++++++++++++++++++++++++++ 1 file changed, 570 insertions(+) create mode 100755 bin/picca_fast_metal_dmat.py diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py new file mode 100755 index 000000000..1d091addc --- /dev/null +++ b/bin/picca_fast_metal_dmat.py @@ -0,0 +1,570 @@ +#!/usr/bin/env python +"""Compute the auto and cross-correlation of delta fields for a list of IGM +absorption. + +This module follow the procedure described in sections 4.3 of du Mas des +Bourboux et al. 2020 (In prep) to compute the distortion matrix +""" +import sys,os +import time +import argparse +import numpy as np +import fitsio + +from picca import constants, cf, utils, io +from picca.utils import userprint + + +def read_stack_deltas_table(filename) : + return fitsio.read(filename,"STACK_DELTAS") + +def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, + out_lambda_abs_1, out_lambda_abs_2, + stack_table_1, stack_table_2, rebin_factor=None): + """Computes the metal distortion matrix. + + Args: + in_lambda_abs_1 : str + Name of absorption in picca.constants in forest pixels from stack 1 (input, i.e. 'true' absorber) + in_lambda_abs_2 : str + Name of absorption in picca.constants in forest pixels from stack 2 (input, i.e. 'true' absorber) + out_lambda_abs_1 : str + Name of absorption in picca.constants in forest pixels from stack 1 (output, i.e. 'assumed' absorber, usually LYA) + out_lambda_abs_2 : str + Name of absorption in picca.constants in forest pixels from stack 2 (output, i.e. 'assumed' absorber, usually LYA) + stack_table_1: table + table with cols LOGLAM and WEIGHT for first series of deltas + stack_table_2: table + table with cols LOGLAM and WEIGHT for second series of deltas + Optionnal : rebin_factor + rebin loglam and weights + Returns: + The distortion matrix data + Note the global picca.cf contains the cosmology and the rp grid + """ + + loglam1=stack_table_1["LOGLAM"] + weight1=stack_table_1["WEIGHT"] + loglam2=stack_table_2["LOGLAM"] + weight2=stack_table_2["WEIGHT"] + if rebin_factor is not None : + n1=loglam1.size + loglam1 = loglam1[:(n1//rebin_factor)*rebin_factor].reshape((n1//rebin_factor),rebin_factor).mean(-1) + weight1 = weight1[:(n1//rebin_factor)*rebin_factor].reshape((n1//rebin_factor),rebin_factor).mean(-1) + n2=loglam2.size + loglam2 = loglam2[:(n2//rebin_factor)*rebin_factor].reshape((n2//rebin_factor),rebin_factor).mean(-1) + weight2 = weight2[:(n2//rebin_factor)*rebin_factor].reshape((n2//rebin_factor),rebin_factor).mean(-1) + + # input + iz1 = (10**loglam1)/constants.ABSORBER_IGM[in_lambda_abs_1] - 1. + iz2 = (10**loglam2)/constants.ABSORBER_IGM[in_lambda_abs_2] - 1. + ir1 = cf.cosmo.get_r_comov(iz1) + ir2 = cf.cosmo.get_r_comov(iz2) + # all pairs + irp = (ir1[:,None]-ir2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) + if not cf.x_correlation: + irp = np.abs(irp) + # output + oz1 = (10**loglam1)/constants.ABSORBER_IGM[out_lambda_abs_1] - 1. + oz2 = (10**loglam2)/constants.ABSORBER_IGM[out_lambda_abs_2] - 1. + or1 = cf.cosmo.get_r_comov(oz1) + or2 = cf.cosmo.get_r_comov(oz2) + # all pairs + orp = (or1[:,None]-or2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) + if not cf.x_correlation: + orp = np.abs(orp) + + # weights + sw = ((weight1*((1+iz1)**(cf.alpha-1)))[:,None]*(weight2*((1+iz2)**(cf.alpha2-1)))[None,:]).ravel() + + # distortion matrix + rpbins = cf.r_par_min + (cf.r_par_max-cf.r_par_min)/cf.num_bins_r_par*np.arange(cf.num_bins_r_par+1) + + # I checked the orientation of the matrix + dmat,_,_ = np.histogram2d(orp,irp,bins=(rpbins,rpbins),weights=sw) + + # normalize (sum of weight should be one for each input rp,rt) + sum_in_weight,_ = np.histogram(irp,bins=rpbins,weights=sw) + dmat *= ((sum_in_weight>0)/(sum_in_weight+(sum_in_weight==0)))[None,:] + + # mean outputs + sum_out_weight,_ = np.histogram(orp,bins=rpbins,weights=sw) + sum_out_weight_rp,_ = np.histogram(orp,bins=rpbins,weights=sw*(orp[None,:].ravel())) + sum_out_weight_z,_ = np.histogram(orp,bins=rpbins,weights=sw*(((oz1[:,None]+oz2[None,:])/2.).ravel())) + r_par_eff = sum_out_weight_rp/(sum_out_weight+(sum_out_weight==0)) + z_eff = sum_out_weight_z/(sum_out_weight+(sum_out_weight==0)) + + r_trans_eff = np.zeros(r_par_eff.shape) + + if True : + # now we will return the full dmat to be consistent with the other computation + # it consists in duplicating the result found to all rt, with output_rt = input_rt + ntot = cf.num_bins_r_par*cf.num_bins_r_trans + + full_dmat = np.zeros((ntot,ntot)) + full_r_par_eff = np.zeros(ntot) + full_r_trans_eff = np.zeros(ntot) + full_z_eff = np.zeros(ntot) + ii = np.arange(cf.num_bins_r_par) + r_trans = (0.5+np.arange(cf.num_bins_r_trans))*cf.r_trans_max/cf.num_bins_r_trans + for j in range(cf.num_bins_r_trans) : + indices = j + cf.num_bins_r_trans * ii + for k,i in zip(indices,ii) : + full_dmat[indices,k] = dmat[ii,i] + full_r_par_eff[indices] = r_par_eff + full_z_eff[indices] = z_eff + full_r_trans_eff[indices] = r_trans[j] + + return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff + + return dmat, r_par_eff, r_trans_eff, z_eff + +def main(cmdargs): + # pylint: disable-msg=too-many-locals,too-many-branches,too-many-statements + """Compute the auto and cross-correlation of delta fields for a list of IGM + absorption.""" + parser = argparse.ArgumentParser( + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + description=('Computes metal matrices ')) + + parser.add_argument('--out', + type=str, + default=None, + required=True, + help='Output file name') + + parser.add_argument('-i','--in-attributes', + type=str, + default=None, + required=True, + help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS containing table with at least rows "LOGLAM" and "WEIGHT"') + + parser.add_argument('--in-attributes2', + type=str, + default=None, + required=False, + help='Path to 2nd delta_attributes.fits.gz file') + + parser.add_argument('--delta-dir', + type=str, + default=None, + required=False, + help='Path to directory with delta*.gz to get the blinding info (default is trying to guess from attributes file)') + + parser.add_argument('--rp-min', + type=float, + default=0., + required=False, + help='Min r-parallel [h^-1 Mpc]') + + parser.add_argument('--rp-max', + type=float, + default=200., + required=False, + help='Max r-parallel [h^-1 Mpc]') + + parser.add_argument('--rt-max', + type=float, + default=200., + required=False, + help='Max r-transverse [h^-1 Mpc]') + + parser.add_argument('--np', + type=int, + default=50, + required=False, + help='Number of r-parallel bins') + + parser.add_argument('--nt', + type=int, + default=50, + required=False, + help='Number of r-transverse bins') + + parser.add_argument( + '--coef-binning-model', + type=int, + default=1, + required=False, + help=('Coefficient multiplying np and nt to get finner binning for the ' + 'model')) + + parser.add_argument( + '--z-cut-min', + type=float, + default=0., + required=False, + help=('Use only pairs of forest x object with the mean of the last ' + 'absorber redshift and the object redshift larger than ' + 'z-cut-min')) + + parser.add_argument( + '--z-cut-max', + type=float, + default=10., + required=False, + help=('Use only pairs of forest x object with the mean of the last ' + 'absorber redshift and the object redshift smaller than ' + 'z-cut-max')) + + parser.add_argument( + '--z-min-sources', + type=float, + default=0., + required=False, + help=('Limit the minimum redshift of the quasars ' + 'used as sources for spectra')) + + parser.add_argument( + '--z-max-sources', + type=float, + default=10., + required=False, + help=('Limit the maximum redshift of the quasars ' + 'used as sources for spectra')) + + parser.add_argument( + '--lambda-abs', + type=str, + default='LYA', + required=False, + help=('Name of the absorption in picca.constants defining the redshift ' + 'of the delta')) + + parser.add_argument( + '--lambda-abs2', + type=str, + default=None, + required=False, + help=('Name of the absorption in picca.constants defining the redshift ' + 'of the 2nd delta')) + + parser.add_argument( + '--abs-igm', + type=str, + default=[], + required=True, + nargs='*', + help=('List of names of metal absorption in picca.constants present in ' + 'forest')) + + parser.add_argument( + '--abs-igm2', + type=str, + default=[], + required=False, + nargs='*', + help=('List of names of metal absorption in picca.constants present in ' + '2nd forest')) + + parser.add_argument('--z-ref', + type=float, + default=2.25, + required=False, + help='Reference redshift') + + parser.add_argument( + '--z-evol', + type=float, + default=2.9, + required=False, + help='Exponent of the redshift evolution of the delta field') + + parser.add_argument( + '--z-evol2', + type=float, + default=2.9, + required=False, + help='Exponent of the redshift evolution of the 2nd delta field') + + parser.add_argument( + '--metal-alpha', + type=float, + default=1., + required=False, + help='Exponent of the redshift evolution of the metal delta field') + + parser.add_argument( + '--fid-Om', + type=float, + default=0.315, + required=False, + help='Omega_matter(z=0) of fiducial LambdaCDM cosmology') + + parser.add_argument( + '--fid-Or', + type=float, + default=0., + required=False, + help='Omega_radiation(z=0) of fiducial LambdaCDM cosmology') + + parser.add_argument('--fid-Ok', + type=float, + default=0., + required=False, + help='Omega_k(z=0) of fiducial LambdaCDM cosmology') + + parser.add_argument( + '--fid-wl', + type=float, + default=-1., + required=False, + help='Equation of state of dark energy of fiducial LambdaCDM cosmology') + + parser.add_argument( + '--unfold-cf', + action='store_true', + required=False, + help=('rp can be positive or negative depending on the relative ' + 'position between absorber1 and absorber2')) + + parser.add_argument('--rebin-factor', + type=int, + default=None, + required=False, + help='Rebin factor for deltas. If not None, deltas will ' + 'be rebinned by that factor') + + args = parser.parse_args(cmdargs) + + # setup variables in module cf + cf.r_par_max = args.rp_max + cf.r_trans_max = args.rt_max + cf.r_par_min = args.rp_min + cf.z_cut_max = args.z_cut_max + cf.z_cut_min = args.z_cut_min + cf.num_bins_r_par = args.np * args.coef_binning_model + cf.num_bins_r_trans = args.nt * args.coef_binning_model + cf.num_model_bins_r_par = args.np * args.coef_binning_model + cf.num_model_bins_r_trans = args.nt * args.coef_binning_model + cf.z_ref = args.z_ref + cf.alpha = args.z_evol + cf.alpha2 = args.z_evol2 + cf.lambda_abs = constants.ABSORBER_IGM[args.lambda_abs] + cf.x_correlation = False # I guess I have to specify this! + + cf.alpha_abs = {} + cf.alpha_abs[args.lambda_abs] = cf.alpha + for metal in args.abs_igm: + cf.alpha_abs[metal] = args.metal_alpha + + # read blinding keyword + if args.delta_dir is None : + args.delta_dir = os.path.dirname(args.in_attributes)+"/../Delta/" + if not os.path.isdir(args.delta_dir) : + print(f"Tried to guess the delta directory (containing the delta*.gz files) but '{args.delta_dir}' is not valid") + print("Please specify the directory with option --delta-dir") + sys.exit(1) + blinding = io.read_blinding(args.delta_dir) + + # load fiducial cosmology + cf.cosmo = constants.Cosmo(Om=args.fid_Om, + Or=args.fid_Or, + Ok=args.fid_Ok, + wl=args.fid_wl, + blinding=blinding) + + + + t0 = time.time() + + ### Read data + stack_table1 = read_stack_deltas_table(args.in_attributes) + + if args.in_attributes2 is not None : + stack_table2 = read_stack_deltas_table(args.in_attributes2) + else : + stack_table2 = stack_table1 # reference to first one + + t1 = time.time() + userprint(f'picca_fast_metal_dmat.py - Time reading data: {(t1-t0)/60:.3f} minutes') + + abs_igm = [args.lambda_abs] + args.abs_igm + userprint("abs_igm = {}".format(abs_igm)) + + if args.lambda_abs2 is None: + args.lambda_abs2 = args.lambda_abs + args.abs_igm2 = args.abs_igm + + abs_igm_2 = [args.lambda_abs2] + args.abs_igm2 + + if cf.x_correlation: + userprint("abs_igm2 = {}".format(abs_igm_2)) + + # intitialize arrays to store the results for the different metal absorption + dmat_all = [] + r_par_all = [] + r_trans_all = [] + z_all = [] + names = [] + #weights_dmat_all = [] + #num_pairs_all = [] + #num_pairs_used_all = [] + + # loop over metals + for index1, abs_igm1 in enumerate(abs_igm): + index0 = index1 + if args.lambda_abs != args.lambda_abs2: + index0 = 0 + for index2, abs_igm2 in enumerate(abs_igm_2[index0:]): + if index1 == 0 and index2 == 0: + continue + + print("Computing",abs_igm1,abs_igm2) + + # this a matrix as a function of rp only + dmat, r_par_eff, r_trans_eff, z_eff = calc_fast_metal_dmat(abs_igm1,abs_igm2, + args.lambda_abs,args.lambda_abs2, + stack_table1, stack_table2, + rebin_factor = args.rebin_factor) + + # add these results to the list ofor the different metal absorption + dmat_all.append(dmat) + r_par_all.append(r_par_eff) + r_trans_all.append(r_trans_eff) + z_all.append(z_eff) + names.append(abs_igm1 + "_" + abs_igm2) + #weights_dmat_all.append(weights_dmat) + #r_trans_all.append(r_trans) NA + #num_pairs_all.append(num_pairs) NA + #num_pairs_used_all.append(num_pairs_used) NA + + t2 = time.time() + userprint(f'picca_fast_metal_dmat.py - Time computing all metal matrices : {(t2-t1)/60:.3f} minutes') + + # save the results + results = fitsio.FITS(args.out, 'rw', clobber=True) + header = [ + { + 'name': 'RPMIN', + 'value': cf.r_par_min, + 'comment': 'Minimum r-parallel [h^-1 Mpc]' + }, + { + 'name': 'RPMAX', + 'value': cf.r_par_max, + 'comment': 'Maximum r-parallel [h^-1 Mpc]' + }, + { + 'name': 'RTMAX', + 'value': cf.r_trans_max, + 'comment': 'Maximum r-transverse [h^-1 Mpc]' + }, + { + 'name': 'NP', + 'value': cf.num_bins_r_par, + 'comment': 'Number of bins in r-parallel' + }, + { + 'name': 'NT', + 'value': cf.num_bins_r_trans, + 'comment': ' Number of bins in r-transverse' + }, + { + 'name': 'COEFMOD', + 'value': args.coef_binning_model, + 'comment': 'Coefficient for model binning' + }, + { + 'name': 'ZCUTMIN', + 'value': cf.z_cut_min, + 'comment': 'Minimum redshift of pairs' + }, + { + 'name': 'ZCUTMAX', + 'value': cf.z_cut_max, + 'comment': 'Maximum redshift of pairs' + }, + { + 'name': 'REJ', + 'value': cf.reject, + 'comment': 'Rejection factor' + }, + { + 'name': 'ALPHAMET', + 'value': args.metal_alpha, + 'comment': 'Evolution of metal bias' + }, { + 'name': 'OMEGAM', + 'value': args.fid_Om, + 'comment': 'Omega_matter(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAR', + 'value': args.fid_Or, + 'comment': 'Omega_radiation(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAK', + 'value': args.fid_Ok, + 'comment': 'Omega_k(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'WL', + 'value': args.fid_wl, + 'comment': 'Equation of state of dark energy of fiducial LambdaCDM cosmology' + }, { + 'name': "BLINDING", + 'value': blinding, + 'comment': 'String specifying the blinding strategy' + } + ] + len_names = np.array([len(name) for name in names]).max() + names = np.array(names, dtype='S' + str(len_names)) + results.write( + [ + #np.array(num_pairs_all), + #np.array(num_pairs_used_all), + np.array(names) + ], + #names=['NPALL', 'NPUSED', 'ABS_IGM'], + names=['ABS_IGM'], + header=header, + #comment=['Number of pairs', 'Number of used pairs', 'Absorption name'], + comment=['Absorption name'], + extname='ATTRI') + + dmat_name = "DM_" + if blinding != "none": + dmat_name += "BLIND_" + names = names.astype(str) + out_list = [] + out_names = [] + out_comment = [] + out_units = [] + for index, name in enumerate(names): + out_names += ['RP_' + name] + out_list += [r_par_all[index]] + out_comment += ['R-parallel'] + out_units += ['h^-1 Mpc'] + + out_names += ['RT_' + name] + out_list += [r_trans_all[index]] + out_comment += ['R-transverse'] + out_units += ['h^-1 Mpc'] + + out_names += ['Z_' + name] + out_list += [z_all[index]] + out_comment += ['Redshift'] + out_units += [''] + + out_names += [dmat_name + name] + out_list += [dmat_all[index]] + out_comment += ['Distortion matrix'] + out_units += [''] + + #out_names += ['WDM_' + name] + #out_list += [weights_dmat_all[index]] + #out_comment += ['Sum of weight'] + #out_units += [''] + + results.write(out_list, + names=out_names, + comment=out_comment, + units=out_units, + extname='MDMAT') + results.close() + + t3 = time.time() + userprint(f'picca_fast_metal_dmat.py - Time total : {(t3-t0)/60:.3f} minutes') + +if __name__ == '__main__': + cmdargs=sys.argv[1:] + main(cmdargs) From a40f18da9f004e6f7cafb981747bb7f892e5d851 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Tue, 5 Sep 2023 10:32:18 -0700 Subject: [PATCH 02/13] remove 'if True' statement used for tests --- bin/picca_fast_metal_dmat.py | 44 ++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index 1d091addc..261a2d02a 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -96,28 +96,28 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, r_trans_eff = np.zeros(r_par_eff.shape) - if True : - # now we will return the full dmat to be consistent with the other computation - # it consists in duplicating the result found to all rt, with output_rt = input_rt - ntot = cf.num_bins_r_par*cf.num_bins_r_trans - - full_dmat = np.zeros((ntot,ntot)) - full_r_par_eff = np.zeros(ntot) - full_r_trans_eff = np.zeros(ntot) - full_z_eff = np.zeros(ntot) - ii = np.arange(cf.num_bins_r_par) - r_trans = (0.5+np.arange(cf.num_bins_r_trans))*cf.r_trans_max/cf.num_bins_r_trans - for j in range(cf.num_bins_r_trans) : - indices = j + cf.num_bins_r_trans * ii - for k,i in zip(indices,ii) : - full_dmat[indices,k] = dmat[ii,i] - full_r_par_eff[indices] = r_par_eff - full_z_eff[indices] = z_eff - full_r_trans_eff[indices] = r_trans[j] - - return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff - - return dmat, r_par_eff, r_trans_eff, z_eff + # we could return the quantities computed as a function of rp only (and not rt): + # return dmat, r_par_eff, r_trans_eff, z_eff + # but for now we will return the full dmat to be consistent with the other computation + # it consists in duplicating the result found to all rt, with output_rt = input_rt + ntot = cf.num_bins_r_par*cf.num_bins_r_trans + + full_dmat = np.zeros((ntot,ntot)) + full_r_par_eff = np.zeros(ntot) + full_r_trans_eff = np.zeros(ntot) + full_z_eff = np.zeros(ntot) + ii = np.arange(cf.num_bins_r_par) + r_trans = (0.5+np.arange(cf.num_bins_r_trans))*cf.r_trans_max/cf.num_bins_r_trans + for j in range(cf.num_bins_r_trans) : + indices = j + cf.num_bins_r_trans * ii + for k,i in zip(indices,ii) : + full_dmat[indices,k] = dmat[ii,i] + full_r_par_eff[indices] = r_par_eff + full_z_eff[indices] = z_eff + full_r_trans_eff[indices] = r_trans[j] + + return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff + def main(cmdargs): # pylint: disable-msg=too-many-locals,too-many-branches,too-many-statements From a7365d87b32efc41f20e3f3c8bcdbb4fc155346f Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Tue, 5 Sep 2023 10:39:56 -0700 Subject: [PATCH 03/13] more explicit variable names --- bin/picca_fast_metal_dmat.py | 58 ++++++++++++++++++------------------ 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index 261a2d02a..566a4ae9b 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -48,49 +48,49 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, loglam2=stack_table_2["LOGLAM"] weight2=stack_table_2["WEIGHT"] if rebin_factor is not None : - n1=loglam1.size - loglam1 = loglam1[:(n1//rebin_factor)*rebin_factor].reshape((n1//rebin_factor),rebin_factor).mean(-1) - weight1 = weight1[:(n1//rebin_factor)*rebin_factor].reshape((n1//rebin_factor),rebin_factor).mean(-1) - n2=loglam2.size - loglam2 = loglam2[:(n2//rebin_factor)*rebin_factor].reshape((n2//rebin_factor),rebin_factor).mean(-1) - weight2 = weight2[:(n2//rebin_factor)*rebin_factor].reshape((n2//rebin_factor),rebin_factor).mean(-1) + size1=loglam1.size + loglam1 = loglam1[:(size1//rebin_factor)*rebin_factor].reshape((size1//rebin_factor),rebin_factor).mean(-1) + weight1 = weight1[:(size1//rebin_factor)*rebin_factor].reshape((size1//rebin_factor),rebin_factor).mean(-1) + size2=loglam2.size + loglam2 = loglam2[:(size2//rebin_factor)*rebin_factor].reshape((size2//rebin_factor),rebin_factor).mean(-1) + weight2 = weight2[:(size2//rebin_factor)*rebin_factor].reshape((size2//rebin_factor),rebin_factor).mean(-1) # input - iz1 = (10**loglam1)/constants.ABSORBER_IGM[in_lambda_abs_1] - 1. - iz2 = (10**loglam2)/constants.ABSORBER_IGM[in_lambda_abs_2] - 1. - ir1 = cf.cosmo.get_r_comov(iz1) - ir2 = cf.cosmo.get_r_comov(iz2) + input_z1 = (10**loglam1)/constants.ABSORBER_IGM[in_lambda_abs_1] - 1. + input_z2 = (10**loglam2)/constants.ABSORBER_IGM[in_lambda_abs_2] - 1. + input_r1 = cf.cosmo.get_r_comov(input_z1) + input_r2 = cf.cosmo.get_r_comov(input_z2) # all pairs - irp = (ir1[:,None]-ir2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) + input_rp = (input_r1[:,None]-input_r2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) if not cf.x_correlation: - irp = np.abs(irp) + input_rp = np.abs(input_rp) # output - oz1 = (10**loglam1)/constants.ABSORBER_IGM[out_lambda_abs_1] - 1. - oz2 = (10**loglam2)/constants.ABSORBER_IGM[out_lambda_abs_2] - 1. - or1 = cf.cosmo.get_r_comov(oz1) - or2 = cf.cosmo.get_r_comov(oz2) + output_z1 = (10**loglam1)/constants.ABSORBER_IGM[out_lambda_abs_1] - 1. + output_z2 = (10**loglam2)/constants.ABSORBER_IGM[out_lambda_abs_2] - 1. + output_r1 = cf.cosmo.get_r_comov(output_z1) + output_r2 = cf.cosmo.get_r_comov(output_z2) # all pairs - orp = (or1[:,None]-or2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) + output_rp = (output_r1[:,None]-output_r2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) if not cf.x_correlation: - orp = np.abs(orp) + output_rp = np.abs(output_rp) # weights - sw = ((weight1*((1+iz1)**(cf.alpha-1)))[:,None]*(weight2*((1+iz2)**(cf.alpha2-1)))[None,:]).ravel() + weights = ((weight1*((1+input_z1)**(cf.alpha-1)))[:,None]*(weight2*((1+input_z2)**(cf.alpha2-1)))[None,:]).ravel() # distortion matrix rpbins = cf.r_par_min + (cf.r_par_max-cf.r_par_min)/cf.num_bins_r_par*np.arange(cf.num_bins_r_par+1) # I checked the orientation of the matrix - dmat,_,_ = np.histogram2d(orp,irp,bins=(rpbins,rpbins),weights=sw) + dmat,_,_ = np.histogram2d(output_rp,input_rp,bins=(rpbins,rpbins),weights=weights) # normalize (sum of weight should be one for each input rp,rt) - sum_in_weight,_ = np.histogram(irp,bins=rpbins,weights=sw) + sum_in_weight,_ = np.histogram(input_rp,bins=rpbins,weights=weights) dmat *= ((sum_in_weight>0)/(sum_in_weight+(sum_in_weight==0)))[None,:] # mean outputs - sum_out_weight,_ = np.histogram(orp,bins=rpbins,weights=sw) - sum_out_weight_rp,_ = np.histogram(orp,bins=rpbins,weights=sw*(orp[None,:].ravel())) - sum_out_weight_z,_ = np.histogram(orp,bins=rpbins,weights=sw*(((oz1[:,None]+oz2[None,:])/2.).ravel())) + sum_out_weight,_ = np.histogram(output_rp,bins=rpbins,weights=weights) + sum_out_weight_rp,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(output_rp[None,:].ravel())) + sum_out_weight_z,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(((output_z1[:,None]+output_z2[None,:])/2.).ravel())) r_par_eff = sum_out_weight_rp/(sum_out_weight+(sum_out_weight==0)) z_eff = sum_out_weight_z/(sum_out_weight+(sum_out_weight==0)) @@ -100,12 +100,12 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, # return dmat, r_par_eff, r_trans_eff, z_eff # but for now we will return the full dmat to be consistent with the other computation # it consists in duplicating the result found to all rt, with output_rt = input_rt - ntot = cf.num_bins_r_par*cf.num_bins_r_trans + num_bins_total = cf.num_bins_r_par*cf.num_bins_r_trans - full_dmat = np.zeros((ntot,ntot)) - full_r_par_eff = np.zeros(ntot) - full_r_trans_eff = np.zeros(ntot) - full_z_eff = np.zeros(ntot) + full_dmat = np.zeros((num_bins_total,num_bins_total)) + full_r_par_eff = np.zeros(num_bins_total) + full_r_trans_eff = np.zeros(num_bins_total) + full_z_eff = np.zeros(num_bins_total) ii = np.arange(cf.num_bins_r_par) r_trans = (0.5+np.arange(cf.num_bins_r_trans))*cf.r_trans_max/cf.num_bins_r_trans for j in range(cf.num_bins_r_trans) : From 4b12b8ab68be4929d628d04fb0a1a8f3f1c930b8 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Tue, 5 Sep 2023 10:41:36 -0700 Subject: [PATCH 04/13] remove commented lines --- bin/picca_fast_metal_dmat.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index 566a4ae9b..b3e18764b 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -397,9 +397,6 @@ def main(cmdargs): r_trans_all = [] z_all = [] names = [] - #weights_dmat_all = [] - #num_pairs_all = [] - #num_pairs_used_all = [] # loop over metals for index1, abs_igm1 in enumerate(abs_igm): @@ -424,10 +421,6 @@ def main(cmdargs): r_trans_all.append(r_trans_eff) z_all.append(z_eff) names.append(abs_igm1 + "_" + abs_igm2) - #weights_dmat_all.append(weights_dmat) - #r_trans_all.append(r_trans) NA - #num_pairs_all.append(num_pairs) NA - #num_pairs_used_all.append(num_pairs_used) NA t2 = time.time() userprint(f'picca_fast_metal_dmat.py - Time computing all metal matrices : {(t2-t1)/60:.3f} minutes') @@ -510,14 +503,10 @@ def main(cmdargs): names = np.array(names, dtype='S' + str(len_names)) results.write( [ - #np.array(num_pairs_all), - #np.array(num_pairs_used_all), np.array(names) ], - #names=['NPALL', 'NPUSED', 'ABS_IGM'], names=['ABS_IGM'], header=header, - #comment=['Number of pairs', 'Number of used pairs', 'Absorption name'], comment=['Absorption name'], extname='ATTRI') From 3a8ace9d08ca6efc343bbfe016ab436a5429f2d9 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Wed, 6 Sep 2023 10:12:30 -0700 Subject: [PATCH 05/13] fast metal matrix computation for cross-correlation --- bin/picca_fast_metal_xdmat.py | 544 ++++++++++++++++++++++++++++++++++ 1 file changed, 544 insertions(+) create mode 100755 bin/picca_fast_metal_xdmat.py diff --git a/bin/picca_fast_metal_xdmat.py b/bin/picca_fast_metal_xdmat.py new file mode 100755 index 000000000..c5e03a7e1 --- /dev/null +++ b/bin/picca_fast_metal_xdmat.py @@ -0,0 +1,544 @@ +#!/usr/bin/env python +"""Compute the distortion matrix of the cross-correlation delta x object for a +list of IGM absorption. + +This module follow the procedure described in sections 4.3 of du Mas des +Bourboux et al. 2020 (In prep) to compute the distortion matrix +""" +import sys,os +import time +import argparse +import numpy as np +import fitsio + +from picca import constants, xcf, io, utils +from picca.utils import userprint + +def read_stack_deltas_table(filename) : + return fitsio.read(filename,"STACK_DELTAS") + +def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, + z_qso, weight_qso, + rebin_factor=None): + """Computes the metal distortion matrix. + + Args: + in_lambda_abs : str + Name of absorption in picca.constants in forest pixels from stack (input, i.e. 'true' absorber) + out_lambda_abs : str + Name of absorption in picca.constants in forest pixels from stack (output, i.e. 'assumed' absorber, usually LYA) + stack_table: table + table with cols LOGLAM and WEIGHT for first series of deltas + z_qso : float 1D array + QSO redshifts + weight_qso : float 1D array + QSO weights (as computed in picca.io.read_objects) + + Optionnal : rebin_factor + rebin loglam and weights + + Returns: + The distortion matrix data + Note the global picca.xcf contains the cosmology, the rp grid, and the QSO catalog + """ + + loglam=stack_table["LOGLAM"] + weight_forest=stack_table["WEIGHT"] + if rebin_factor is not None : + size=loglam.size + loglam = loglam[:(size//rebin_factor)*rebin_factor].reshape((size//rebin_factor),rebin_factor).mean(-1) + weight_forest = weight_forest[:(size//rebin_factor)*rebin_factor].reshape((size//rebin_factor),rebin_factor).mean(-1) + + # input + input_zf = (10**loglam)/constants.ABSORBER_IGM[in_lambda_abs] - 1. # redshift in the forest + input_rf = xcf.cosmo.get_r_comov(input_zf) + + rq = xcf.cosmo.get_r_comov(z_qso) # comoving distance for qso + + # all pairs + input_rp = (input_rf[:,None]-rq[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) + + # output + output_zf = (10**loglam)/constants.ABSORBER_IGM[out_lambda_abs] - 1. + output_rf = xcf.cosmo.get_r_comov(output_zf) + + # all pairs + output_rp = (output_rf[:,None]-rq[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) + + # weights + alpha = xcf.alpha_abs[out_lambda_abs] # this is how the correlation is computed + weights = ((weight_forest*((1+input_zf)**(alpha-1)))[:,None]*weight_qso[None,:]).ravel() # qso weights have already been scale with (1+z) in picca.io.read_objects + + # distortion matrix + rpbins = xcf.r_par_min + (xcf.r_par_max-xcf.r_par_min)/xcf.num_bins_r_par*np.arange(xcf.num_bins_r_par+1) + + # I checked the orientation of the matrix + dmat,_,_ = np.histogram2d(output_rp,input_rp,bins=(rpbins,rpbins),weights=weights) + + # normalize (sum of weight should be one for each input rp,rt) + sum_in_weight,_ = np.histogram(input_rp,bins=rpbins,weights=weights) + dmat *= ((sum_in_weight>0)/(sum_in_weight+(sum_in_weight==0)))[None,:] + + # mean outputs + sum_out_weight,_ = np.histogram(output_rp,bins=rpbins,weights=weights) + sum_out_weight_rp,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(output_rp[None,:].ravel())) + sum_out_weight_z,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(((output_zf[:,None]+z_qso[None,:])/2.).ravel())) + r_par_eff = sum_out_weight_rp/(sum_out_weight+(sum_out_weight==0)) + z_eff = sum_out_weight_z/(sum_out_weight+(sum_out_weight==0)) + + r_trans_eff = np.zeros(r_par_eff.shape) + + # we could return the quantities computed as a function of rp only (and not rt): + # return dmat, r_par_eff, r_trans_eff, z_eff + # but for now we will return the full dmat to be consistent with the other computation + # it consists in duplicating the result found to all rt, with output_rt = input_rt + num_bins_total = xcf.num_bins_r_par*xcf.num_bins_r_trans + + full_dmat = np.zeros((num_bins_total,num_bins_total)) + full_r_par_eff = np.zeros(num_bins_total) + full_r_trans_eff = np.zeros(num_bins_total) + full_z_eff = np.zeros(num_bins_total) + ii = np.arange(xcf.num_bins_r_par) + r_trans = (0.5+np.arange(xcf.num_bins_r_trans))*xcf.r_trans_max/xcf.num_bins_r_trans + for j in range(xcf.num_bins_r_trans) : + indices = j + xcf.num_bins_r_trans * ii + for k,i in zip(indices,ii) : + full_dmat[indices,k] = dmat[ii,i] + full_r_par_eff[indices] = r_par_eff + full_z_eff[indices] = z_eff + full_r_trans_eff[indices] = r_trans[j] + + return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff + +def main(cmdargs): + """Compute the metal matrix of the cross-correlation delta x object for + a list of IGM absorption.""" + parser = argparse.ArgumentParser( + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + description=('Computes metal matrices for the cross-correlation ' + 'delta x object for a list of IGM absorption.')) + + parser.add_argument('--out', + type=str, + default=None, + required=True, + help='Output file name') + + parser.add_argument('-i','--in-attributes', + type=str, + default=None, + required=True, + help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS containing table with at least rows "LOGLAM" and "WEIGHT"') + + parser.add_argument('--delta-dir', + type=str, + default=None, + required=False, + help='Path to directory with delta*.gz to get the blinding info (default is trying to guess from attributes file)') + + parser.add_argument('--drq', + type=str, + default=None, + required=True, + help='Catalog of objects in DRQ format') + + parser.add_argument('--mode', + type=str, + default='sdss', + choices=['sdss','desi','desi_mocks','desi_healpix'], + required=False, + help='type of catalog supplied, default sdss') + + parser.add_argument('--rp-min', + type=float, + default=-200., + required=False, + help='Min r-parallel [h^-1 Mpc]') + + parser.add_argument('--rp-max', + type=float, + default=200., + required=False, + help='Max r-parallel [h^-1 Mpc]') + + parser.add_argument('--rt-max', + type=float, + default=200., + required=False, + help='Max r-transverse [h^-1 Mpc]') + + parser.add_argument('--np', + type=int, + default=100, + required=False, + help='Number of r-parallel bins') + + parser.add_argument('--nt', + type=int, + default=50, + required=False, + help='Number of r-transverse bins') + + parser.add_argument( + '--coef-binning-model', + type=int, + default=1, + required=False, + help=('Coefficient multiplying np and nt to get finner binning for the ' + 'model')) + + parser.add_argument('--z-min-obj', + type=float, + default=0, + required=False, + help='Min redshift for object field') + + parser.add_argument('--z-max-obj', + type=float, + default=10, + required=False, + help='Max redshift for object field') + + parser.add_argument( + '--z-cut-min', + type=float, + default=0., + required=False, + help=('Use only pairs of forest x object with the mean of the last ' + 'absorber redshift and the object redshift larger than ' + 'z-cut-min')) + + parser.add_argument( + '--z-cut-max', + type=float, + default=10., + required=False, + help=('Use only pairs of forest x object with the mean of the last ' + 'absorber redshift and the object redshift smaller than ' + 'z-cut-max')) + + parser.add_argument( + '--z-min-sources', + type=float, + default=0., + required=False, + help=('Limit the minimum redshift of the quasars ' + 'used as sources for spectra')) + + parser.add_argument( + '--z-max-sources', + type=float, + default=10., + required=False, + help=('Limit the maximum redshift of the quasars ' + 'used as sources for spectra')) + + parser.add_argument( + '--lambda-abs', + type=str, + default='LYA', + required=False, + help=('Name of the absorption in picca.constants defining the redshift ' + 'of the delta')) + + parser.add_argument('--obj-name', + type=str, + default='QSO', + required=False, + help='Name of the object tracer') + + parser.add_argument( + '--abs-igm', + type=str, + default=None, + required=False, + nargs='*', + help='List of names of metal absorption in picca.constants') + + parser.add_argument('--z-ref', + type=float, + default=2.25, + required=False, + help='Reference redshift') + + parser.add_argument( + '--z-evol-del', + type=float, + default=2.9, + required=False, + help='Exponent of the redshift evolution of the delta field') + + parser.add_argument( + '--z-evol-obj', + type=float, + default=1.44, + required=False, + help='Exponent of the redshift evolution of the object field') + + parser.add_argument( + '--metal-alpha', + type=float, + default=1., + required=False, + help='Exponent of the redshift evolution of the metal delta field') + + parser.add_argument( + '--fid-Om', + type=float, + default=0.315, + required=False, + help='Omega_matter(z=0) of fiducial LambdaCDM cosmology') + + parser.add_argument( + '--fid-Or', + type=float, + default=0., + required=False, + help='Omega_radiation(z=0) of fiducial LambdaCDM cosmology') + + parser.add_argument('--fid-Ok', + type=float, + default=0., + required=False, + help='Omega_k(z=0) of fiducial LambdaCDM cosmology') + + parser.add_argument( + '--fid-wl', + type=float, + default=-1., + required=False, + help='Equation of state of dark energy of fiducial LambdaCDM cosmology') + + parser.add_argument('--rebin-factor', + type=int, + default=None, + required=False, + help='Rebin factor for deltas. If not None, deltas will ' + 'be rebinned by that factor') + + args = parser.parse_args(cmdargs) + + # setup variables in module xcf + xcf.r_par_max = args.rp_max + xcf.r_par_min = args.rp_min + xcf.r_trans_max = args.rt_max + xcf.z_cut_max = args.z_cut_max + xcf.z_cut_min = args.z_cut_min + xcf.num_bins_r_par = args.np * args.coef_binning_model + xcf.num_bins_r_trans = args.nt * args.coef_binning_model + xcf.num_model_bins_r_par = args.np * args.coef_binning_model + xcf.num_model_bins_r_trans = args.nt * args.coef_binning_model + xcf.z_ref = args.z_ref + xcf.lambda_abs = constants.ABSORBER_IGM[args.lambda_abs] + + xcf.alpha_abs = {} + xcf.alpha_abs[args.lambda_abs] = args.z_evol_del + for metal in args.abs_igm: + xcf.alpha_abs[metal] = args.metal_alpha + + # read blinding keyword + if args.delta_dir is None : + args.delta_dir = os.path.dirname(args.in_attributes)+"/../Delta/" + if not os.path.isdir(args.delta_dir) : + userprint(f"Tried to guess the delta directory (containing the delta*.gz files) but '{args.delta_dir}' is not valid") + userprint("Please specify the directory with option --delta-dir") + sys.exit(1) + blinding = io.read_blinding(args.delta_dir) + + # load fiducial cosmology + cosmo = constants.Cosmo(Om=args.fid_Om, + Or=args.fid_Or, + Ok=args.fid_Ok, + wl=args.fid_wl, + blinding=blinding) + xcf.cosmo = cosmo + + t0 = time.time() + + ### Read data + stack_table = read_stack_deltas_table(args.in_attributes) + + # read objets + if 0 : + # dummy nside, we do not use this pixelization + # but we still want to use the existing picca I/O routines + nside = 4 + objs, z_min2 = io.read_objects(args.drq, nside, args.z_min_obj, + args.z_max_obj, args.z_evol_obj, args.z_ref, + cosmo, mode=args.mode) + # collect redshifts and weights + z_qso = [] + weight_qso = [] + for qsos in objs.values() : # objs is a dictionary of healpixels + z_qso += [qso.z_qso for qso in qsos] + weight_qso += [qso.weights for qso in qsos] + # convert to numpy array + z_qso = np.hstack(z_qso) + weight_qso = np.hstack(weight_qso) + else : + catalog = io.read_drq(args.drq, z_min=args.z_min_obj, + z_max=args.z_max_obj, keep_bal=True, mode=args.mode) + z_qso = catalog['Z'] + weight_qso = ((1. + z_qso) / (1. + args.z_ref))**(args.z_evol_obj - 1.) + + + zbins=100 + userprint(f"Use histogram of QSO redshifts with {zbins} bins") + hw,zbins = np.histogram(z_qso,bins=zbins,weights=weight_qso) + hwz,_ = np.histogram(z_qso,bins=zbins,weights=weight_qso*z_qso) + ok = (hw>0) + z_qso = hwz[ok]/hw[ok] # weighted mean in bins + weight_qso = hw[ok] + + t1 = time.time() + userprint(f'picca_metal_xdmat.py - Time reading data: {(t1-t0)/60:.3f} minutes') + + # intitialize arrays to store the results for the different metal absorption + dmat_all = [] + r_par_all = [] + r_trans_all = [] + z_all = [] + names = [] + + # loop over metals + for index, abs_igm in enumerate(args.abs_igm): + + userprint("Computing",abs_igm) + + # this a matrix as a function of rp only + dmat, r_par_eff, r_trans_eff, z_eff = calc_fast_metal_dmat(abs_igm, + args.lambda_abs, + stack_table, + z_qso, + weight_qso, + rebin_factor = args.rebin_factor) + + # add these results to the list ofor the different metal absorption + dmat_all.append(dmat) + r_par_all.append(r_par_eff) + r_trans_all.append(r_trans_eff) + z_all.append(z_eff) + names.append(abs_igm) + + t2 = time.time() + userprint(f'picca_metal_xdmat.py - Time computing all metal matrix: {(t2-t1)/60:.3f} minutes') + + # save the results + results = fitsio.FITS(args.out, 'rw', clobber=True) + header = [ + { + 'name': 'RPMIN', + 'value': xcf.r_par_min, + 'comment': 'Minimum r-parallel [h^-1 Mpc]' + }, + { + 'name': 'RPMAX', + 'value': xcf.r_par_max, + 'comment': 'Maximum r-parallel [h^-1 Mpc]' + }, + { + 'name': 'RTMAX', + 'value': xcf.r_trans_max, + 'comment': 'Maximum r-transverse [h^-1 Mpc]' + }, + { + 'name': 'NP', + 'value': xcf.num_bins_r_par, + 'comment': 'Number of bins in r-parallel' + }, + { + 'name': 'NT', + 'value': xcf.num_bins_r_trans, + 'comment': 'Number of bins in r-transverse' + }, + { + 'name': 'COEFMOD', + 'value': args.coef_binning_model, + 'comment': 'Coefficient for model binning' + }, + { + 'name': 'ZCUTMIN', + 'value': xcf.z_cut_min, + 'comment': 'Minimum redshift of pairs' + }, + { + 'name': 'ZCUTMAX', + 'value': xcf.z_cut_max, + 'comment': 'Maximum redshift of pairs' + }, + { + 'name': 'OMEGAM', + 'value': args.fid_Om, + 'comment': 'Omega_matter(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAR', + 'value': args.fid_Or, + 'comment': 'Omega_radiation(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAK', + 'value': args.fid_Ok, + 'comment': 'Omega_k(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'WL', + 'value': args.fid_wl, + 'comment': 'Equation of state of dark energy of fiducial LambdaCDM cosmology' + }, { + 'name': "BLINDING", + 'value': blinding, + 'comment': 'String specifying the blinding strategy' + } + ] + len_names = np.array([len(name) for name in names]).max() + names = np.array(names, dtype='S' + str(len_names)) + results.write( + [ + np.array(names) + ], + names=['ABS_IGM'], + header=header, + comment=['Number of pairs', 'Number of used pairs', 'Absorption name'], + extname='ATTRI') + + dmat_name = "DM_" + if blinding != "none": + dmat_name += "BLIND_" + names = names.astype(str) + out_list = [] + out_names = [] + out_comment = [] + out_units = [] + for index, name in enumerate(names): + out_names += ['RP_' + args.obj_name + '_' + name] + out_list += [r_par_all[index]] + out_comment += ['R-parallel'] + out_units += ['h^-1 Mpc'] + + out_names += ['RT_' + args.obj_name + '_' + name] + out_list += [r_trans_all[index]] + out_comment += ['R-transverse'] + out_units += ['h^-1 Mpc'] + + out_names += ['Z_' + args.obj_name + '_' + name] + out_list += [z_all[index]] + out_comment += ['Redshift'] + out_units += [''] + + out_names += [dmat_name + args.obj_name + '_' + name] + out_list += [dmat_all[index]] + out_comment += ['Distortion matrix'] + out_units += [''] + + results.write(out_list, + names=out_names, + comment=out_comment, + units=out_units, + extname='MDMAT') + results.close() + + t3 = time.time() + userprint(f'picca_metal_xdmat.py - Time total: {(t3-t0)/60:.3f} minutes') + + +if __name__ == '__main__': + cmdargs=sys.argv[1:] + main(cmdargs) From 4961ef30002cc5c4184e5a485a85fa5462a73afb Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Wed, 6 Sep 2023 10:12:51 -0700 Subject: [PATCH 06/13] userprint instead of print --- bin/picca_fast_metal_dmat.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index b3e18764b..4d91b1f2f 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -352,8 +352,8 @@ def main(cmdargs): if args.delta_dir is None : args.delta_dir = os.path.dirname(args.in_attributes)+"/../Delta/" if not os.path.isdir(args.delta_dir) : - print(f"Tried to guess the delta directory (containing the delta*.gz files) but '{args.delta_dir}' is not valid") - print("Please specify the directory with option --delta-dir") + userprint(f"Tried to guess the delta directory (containing the delta*.gz files) but '{args.delta_dir}' is not valid") + userprint("Please specify the directory with option --delta-dir") sys.exit(1) blinding = io.read_blinding(args.delta_dir) @@ -407,7 +407,7 @@ def main(cmdargs): if index1 == 0 and index2 == 0: continue - print("Computing",abs_igm1,abs_igm2) + userprint("Computing",abs_igm1,abs_igm2) # this a matrix as a function of rp only dmat, r_par_eff, r_trans_eff, z_eff = calc_fast_metal_dmat(abs_igm1,abs_igm2, From 1b4a704dd10467eaacab45448a143422fc53eb6c Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Wed, 6 Sep 2023 10:45:59 -0700 Subject: [PATCH 07/13] qso-z-bins is a parameter with default = 1000 --- bin/picca_fast_metal_xdmat.py | 35 +++++++++++------------------------ 1 file changed, 11 insertions(+), 24 deletions(-) diff --git a/bin/picca_fast_metal_xdmat.py b/bin/picca_fast_metal_xdmat.py index c5e03a7e1..708368b37 100755 --- a/bin/picca_fast_metal_xdmat.py +++ b/bin/picca_fast_metal_xdmat.py @@ -315,6 +315,11 @@ def main(cmdargs): required=False, help='Rebin factor for deltas. If not None, deltas will ' 'be rebinned by that factor') + parser.add_argument('--qso-z-bins', + type=int, + default=1000, + required=False, + help='Bins for the distribution of QSO redshifts') args = parser.parse_args(cmdargs) @@ -359,30 +364,12 @@ def main(cmdargs): stack_table = read_stack_deltas_table(args.in_attributes) # read objets - if 0 : - # dummy nside, we do not use this pixelization - # but we still want to use the existing picca I/O routines - nside = 4 - objs, z_min2 = io.read_objects(args.drq, nside, args.z_min_obj, - args.z_max_obj, args.z_evol_obj, args.z_ref, - cosmo, mode=args.mode) - # collect redshifts and weights - z_qso = [] - weight_qso = [] - for qsos in objs.values() : # objs is a dictionary of healpixels - z_qso += [qso.z_qso for qso in qsos] - weight_qso += [qso.weights for qso in qsos] - # convert to numpy array - z_qso = np.hstack(z_qso) - weight_qso = np.hstack(weight_qso) - else : - catalog = io.read_drq(args.drq, z_min=args.z_min_obj, - z_max=args.z_max_obj, keep_bal=True, mode=args.mode) - z_qso = catalog['Z'] - weight_qso = ((1. + z_qso) / (1. + args.z_ref))**(args.z_evol_obj - 1.) - - - zbins=100 + catalog = io.read_drq(args.drq, z_min=args.z_min_obj, + z_max=args.z_max_obj, keep_bal=True, mode=args.mode) + z_qso = catalog['Z'] + weight_qso = ((1. + z_qso) / (1. + args.z_ref))**(args.z_evol_obj - 1.) + + zbins=args.qso_z_bins userprint(f"Use histogram of QSO redshifts with {zbins} bins") hw,zbins = np.histogram(z_qso,bins=zbins,weights=weight_qso) hwz,_ = np.histogram(z_qso,bins=zbins,weights=weight_qso*z_qso) From a849666ba2b03c375efba6b19d01263762891955 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Wed, 6 Sep 2023 15:49:43 -0700 Subject: [PATCH 08/13] change in a comment --- bin/picca_fast_metal_xdmat.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bin/picca_fast_metal_xdmat.py b/bin/picca_fast_metal_xdmat.py index 708368b37..019206a02 100755 --- a/bin/picca_fast_metal_xdmat.py +++ b/bin/picca_fast_metal_xdmat.py @@ -67,7 +67,7 @@ def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, # weights alpha = xcf.alpha_abs[out_lambda_abs] # this is how the correlation is computed - weights = ((weight_forest*((1+input_zf)**(alpha-1)))[:,None]*weight_qso[None,:]).ravel() # qso weights have already been scale with (1+z) in picca.io.read_objects + weights = ((weight_forest*((1+input_zf)**(alpha-1)))[:,None]*weight_qso[None,:]).ravel() # qso weights have already been scale with (1+z) # distortion matrix rpbins = xcf.r_par_min + (xcf.r_par_max-xcf.r_par_min)/xcf.num_bins_r_par*np.arange(xcf.num_bins_r_par+1) From bc8256f802573eaefbfafaaa62f5b65cb794dd36 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Wed, 6 Sep 2023 17:11:59 -0700 Subject: [PATCH 09/13] fix scaling of weights with (1+z). still mismatch for lya x qso with std. computation --- bin/picca_fast_metal_dmat.py | 10 +++++++++- bin/picca_fast_metal_xdmat.py | 9 +++++++-- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index 4d91b1f2f..e759d59e9 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -75,7 +75,15 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, output_rp = np.abs(output_rp) # weights - weights = ((weight1*((1+input_z1)**(cf.alpha-1)))[:,None]*(weight2*((1+input_z2)**(cf.alpha2-1)))[None,:]).ravel() + # alpha_in: in (1+z)^(alpha_in-1) is a scaling used to model how the metal contribution evolves with redshift (by default alpha=1 so that this has no effect) + # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and used for the correlation function. It also has to be applied here. + # we have them for both forests (1 and 2) + alpha_in_1 = cf.alpha_abs[in_lambda_abs_1] + alpha_out_1 = cf.alpha # = args.z_evol in main function + alpha_in_2 = cf.alpha_abs[in_lambda_abs_2] + alpha_out_2 = cf.alpha2 # = args.z_evol2 in main function + # so here we have to apply both scalings (in the original code : alpha_in is applied in cf.calc_metal_dmat and alpha_out in io.read_deltas) + weights = ((weight1*((1+input_z1)**(alpha_in_1+alpha_out_1-2)))[:,None]*(weight2*((1+input_z2)**(alpha_in_2+alpha_out_2-2)))[None,:]).ravel() # distortion matrix rpbins = cf.r_par_min + (cf.r_par_max-cf.r_par_min)/cf.num_bins_r_par*np.arange(cf.num_bins_r_par+1) diff --git a/bin/picca_fast_metal_xdmat.py b/bin/picca_fast_metal_xdmat.py index 019206a02..f77f4b67a 100755 --- a/bin/picca_fast_metal_xdmat.py +++ b/bin/picca_fast_metal_xdmat.py @@ -66,8 +66,13 @@ def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, output_rp = (output_rf[:,None]-rq[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) # weights - alpha = xcf.alpha_abs[out_lambda_abs] # this is how the correlation is computed - weights = ((weight_forest*((1+input_zf)**(alpha-1)))[:,None]*weight_qso[None,:]).ravel() # qso weights have already been scale with (1+z) + # alpha_in: in (1+z)^(alpha_in-1) is a scaling used to model how the metal contribution evolves with redshift (by default alpha=1 so that this has no effect) + alpha_in = xcf.alpha_abs[in_lambda_abs] + # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and used for the correlation function. It also has to be applied here. + alpha_out = xcf.alpha_abs[out_lambda_abs] + # so here we have to apply both scalings (in the original code : alpha_in is applied in xcf.calc_metal_dmat and alpha_out in io.read_deltas) + # qso weights have already been scaled with (1+z)^alpha_obj + weights = ((weight_forest*((1+input_zf)**(alpha_in+alpha_out-2)))[:,None]*weight_qso[None,:]).ravel() # distortion matrix rpbins = xcf.r_par_min + (xcf.r_par_max-xcf.r_par_min)/xcf.num_bins_r_par*np.arange(xcf.num_bins_r_par+1) From 98b6c27255a31c56f6fa61beb3a2dd1e71e1a428 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Tue, 12 Sep 2023 10:29:26 -0700 Subject: [PATCH 10/13] making pylint happier --- bin/picca_fast_metal_dmat.py | 70 ++++++++++++++++++++++-------------- 1 file changed, 43 insertions(+), 27 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index e759d59e9..09a6b6a85 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -5,17 +5,25 @@ This module follow the procedure described in sections 4.3 of du Mas des Bourboux et al. 2020 (In prep) to compute the distortion matrix """ -import sys,os +import sys +import os import time import argparse import numpy as np import fitsio -from picca import constants, cf, utils, io +from picca import constants, cf, io from picca.utils import userprint - def read_stack_deltas_table(filename) : + """ + Read stack. + + Args: + filename : std , path + Returns: + table as numpy.ndarray + """ return fitsio.read(filename,"STACK_DELTAS") def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, @@ -25,13 +33,17 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, Args: in_lambda_abs_1 : str - Name of absorption in picca.constants in forest pixels from stack 1 (input, i.e. 'true' absorber) + Name of absorption in picca.constants in forest pixels from stack 1 + (input, i.e. 'true' absorber) in_lambda_abs_2 : str - Name of absorption in picca.constants in forest pixels from stack 2 (input, i.e. 'true' absorber) + Name of absorption in picca.constants in forest pixels from stack 2 + (input, i.e. 'true' absorber) out_lambda_abs_1 : str - Name of absorption in picca.constants in forest pixels from stack 1 (output, i.e. 'assumed' absorber, usually LYA) + Name of absorption in picca.constants in forest pixels from stack 1 + (output, i.e. 'assumed' absorber, usually LYA) out_lambda_abs_2 : str - Name of absorption in picca.constants in forest pixels from stack 2 (output, i.e. 'assumed' absorber, usually LYA) + Name of absorption in picca.constants in forest pixels from stack 2 + (output, i.e. 'assumed' absorber, usually LYA) stack_table_1: table table with cols LOGLAM and WEIGHT for first series of deltas stack_table_2: table @@ -75,18 +87,23 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, output_rp = np.abs(output_rp) # weights - # alpha_in: in (1+z)^(alpha_in-1) is a scaling used to model how the metal contribution evolves with redshift (by default alpha=1 so that this has no effect) - # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and used for the correlation function. It also has to be applied here. + # alpha_in: in (1+z)^(alpha_in-1) is a scaling used to model how the metal contribution + # evolves with redshift (by default alpha_in=1 so that this has no effect) + # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and + # used for the correlation function. It also has to be applied here. # we have them for both forests (1 and 2) alpha_in_1 = cf.alpha_abs[in_lambda_abs_1] alpha_out_1 = cf.alpha # = args.z_evol in main function alpha_in_2 = cf.alpha_abs[in_lambda_abs_2] alpha_out_2 = cf.alpha2 # = args.z_evol2 in main function - # so here we have to apply both scalings (in the original code : alpha_in is applied in cf.calc_metal_dmat and alpha_out in io.read_deltas) - weights = ((weight1*((1+input_z1)**(alpha_in_1+alpha_out_1-2)))[:,None]*(weight2*((1+input_z2)**(alpha_in_2+alpha_out_2-2)))[None,:]).ravel() + # so here we have to apply both scalings (in the original code : + # alpha_in is applied in cf.calc_metal_dmat and alpha_out in io.read_deltas) + weights = ((weight1*((1+input_z1)**(alpha_in_1+alpha_out_1-2)))[:,None] + *(weight2*((1+input_z2)**(alpha_in_2+alpha_out_2-2)))[None,:]).ravel() # distortion matrix - rpbins = cf.r_par_min + (cf.r_par_max-cf.r_par_min)/cf.num_bins_r_par*np.arange(cf.num_bins_r_par+1) + rpbins = cf.r_par_min \ + + (cf.r_par_max-cf.r_par_min)/cf.num_bins_r_par*np.arange(cf.num_bins_r_par+1) # I checked the orientation of the matrix dmat,_,_ = np.histogram2d(output_rp,input_rp,bins=(rpbins,rpbins),weights=weights) @@ -102,8 +119,6 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, r_par_eff = sum_out_weight_rp/(sum_out_weight+(sum_out_weight==0)) z_eff = sum_out_weight_z/(sum_out_weight+(sum_out_weight==0)) - r_trans_eff = np.zeros(r_par_eff.shape) - # we could return the quantities computed as a function of rp only (and not rt): # return dmat, r_par_eff, r_trans_eff, z_eff # but for now we will return the full dmat to be consistent with the other computation @@ -114,12 +129,12 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, full_r_par_eff = np.zeros(num_bins_total) full_r_trans_eff = np.zeros(num_bins_total) full_z_eff = np.zeros(num_bins_total) - ii = np.arange(cf.num_bins_r_par) + r_par_indices = np.arange(cf.num_bins_r_par) r_trans = (0.5+np.arange(cf.num_bins_r_trans))*cf.r_trans_max/cf.num_bins_r_trans for j in range(cf.num_bins_r_trans) : - indices = j + cf.num_bins_r_trans * ii - for k,i in zip(indices,ii) : - full_dmat[indices,k] = dmat[ii,i] + indices = j + cf.num_bins_r_trans * r_par_indices + for k,i in zip(indices,r_par_indices) : + full_dmat[indices,k] = dmat[r_par_indices,i] full_r_par_eff[indices] = r_par_eff full_z_eff[indices] = z_eff full_r_trans_eff[indices] = r_trans[j] @@ -127,13 +142,13 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff -def main(cmdargs): +def main(): # pylint: disable-msg=too-many-locals,too-many-branches,too-many-statements """Compute the auto and cross-correlation of delta fields for a list of IGM absorption.""" parser = argparse.ArgumentParser( formatter_class=argparse.ArgumentDefaultsHelpFormatter, - description=('Computes metal matrices ')) + description='Computes metal matrices') parser.add_argument('--out', type=str, @@ -145,7 +160,8 @@ def main(cmdargs): type=str, default=None, required=True, - help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS containing table with at least rows "LOGLAM" and "WEIGHT"') + help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS' + ' containing table with at least rows "LOGLAM" and "WEIGHT"') parser.add_argument('--in-attributes2', type=str, @@ -157,7 +173,8 @@ def main(cmdargs): type=str, default=None, required=False, - help='Path to directory with delta*.gz to get the blinding info (default is trying to guess from attributes file)') + help='Path to directory with delta*.gz to get the blinding info' + ' (default is trying to guess from attributes file)') parser.add_argument('--rp-min', type=float, @@ -323,8 +340,8 @@ def main(cmdargs): '--unfold-cf', action='store_true', required=False, - help=('rp can be positive or negative depending on the relative ' - 'position between absorber1 and absorber2')) + help='rp can be positive or negative depending on the relative ' + 'position between absorber1 and absorber2') parser.add_argument('--rebin-factor', type=int, @@ -333,7 +350,7 @@ def main(cmdargs): help='Rebin factor for deltas. If not None, deltas will ' 'be rebinned by that factor') - args = parser.parse_args(cmdargs) + args = parser.parse_args() # setup variables in module cf cf.r_par_max = args.rp_max @@ -563,5 +580,4 @@ def main(cmdargs): userprint(f'picca_fast_metal_dmat.py - Time total : {(t3-t0)/60:.3f} minutes') if __name__ == '__main__': - cmdargs=sys.argv[1:] - main(cmdargs) + main() From 0ed478376a08e433cae3c7cf437894bbcb1b7635 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Tue, 12 Sep 2023 10:40:38 -0700 Subject: [PATCH 11/13] making pylint happier --- bin/picca_fast_metal_xdmat.py | 69 +++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 28 deletions(-) diff --git a/bin/picca_fast_metal_xdmat.py b/bin/picca_fast_metal_xdmat.py index f77f4b67a..36145a78e 100755 --- a/bin/picca_fast_metal_xdmat.py +++ b/bin/picca_fast_metal_xdmat.py @@ -5,16 +5,25 @@ This module follow the procedure described in sections 4.3 of du Mas des Bourboux et al. 2020 (In prep) to compute the distortion matrix """ -import sys,os +import sys +import os import time import argparse import numpy as np import fitsio -from picca import constants, xcf, io, utils +from picca import constants, xcf, io from picca.utils import userprint def read_stack_deltas_table(filename) : + """ + Read stack. + + Args: + filename : std , path + Returns: + table as numpy.ndarray + """ return fitsio.read(filename,"STACK_DELTAS") def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, @@ -24,9 +33,11 @@ def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, Args: in_lambda_abs : str - Name of absorption in picca.constants in forest pixels from stack (input, i.e. 'true' absorber) + Name of absorption in picca.constants in forest pixels from stack + (input, i.e. 'true' absorber) out_lambda_abs : str - Name of absorption in picca.constants in forest pixels from stack (output, i.e. 'assumed' absorber, usually LYA) + Name of absorption in picca.constants in forest pixels from stack + (output, i.e. 'assumed' absorber, usually LYA) stack_table: table table with cols LOGLAM and WEIGHT for first series of deltas z_qso : float 1D array @@ -53,24 +64,27 @@ def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, input_zf = (10**loglam)/constants.ABSORBER_IGM[in_lambda_abs] - 1. # redshift in the forest input_rf = xcf.cosmo.get_r_comov(input_zf) - rq = xcf.cosmo.get_r_comov(z_qso) # comoving distance for qso + r_qso = xcf.cosmo.get_r_comov(z_qso) # comoving distance for qso # all pairs - input_rp = (input_rf[:,None]-rq[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) + input_rp = (input_rf[:,None]-r_qso[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) # output output_zf = (10**loglam)/constants.ABSORBER_IGM[out_lambda_abs] - 1. output_rf = xcf.cosmo.get_r_comov(output_zf) # all pairs - output_rp = (output_rf[:,None]-rq[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) + output_rp = (output_rf[:,None]-r_qso[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) # weights - # alpha_in: in (1+z)^(alpha_in-1) is a scaling used to model how the metal contribution evolves with redshift (by default alpha=1 so that this has no effect) + # alpha_in: in (1+z)^(alpha_in-1) is a scaling used to model how the metal contribution + # evolves with redshift (by default alpha=1 so that this has no effect) alpha_in = xcf.alpha_abs[in_lambda_abs] - # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and used for the correlation function. It also has to be applied here. + # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and + # used for the correlation function. It also has to be applied here. alpha_out = xcf.alpha_abs[out_lambda_abs] - # so here we have to apply both scalings (in the original code : alpha_in is applied in xcf.calc_metal_dmat and alpha_out in io.read_deltas) + # so here we have to apply both scalings (in the original code : alpha_in is applied in + # xcf.calc_metal_dmat and alpha_out in io.read_deltas) # qso weights have already been scaled with (1+z)^alpha_obj weights = ((weight_forest*((1+input_zf)**(alpha_in+alpha_out-2)))[:,None]*weight_qso[None,:]).ravel() @@ -91,8 +105,6 @@ def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, r_par_eff = sum_out_weight_rp/(sum_out_weight+(sum_out_weight==0)) z_eff = sum_out_weight_z/(sum_out_weight+(sum_out_weight==0)) - r_trans_eff = np.zeros(r_par_eff.shape) - # we could return the quantities computed as a function of rp only (and not rt): # return dmat, r_par_eff, r_trans_eff, z_eff # but for now we will return the full dmat to be consistent with the other computation @@ -103,19 +115,19 @@ def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, full_r_par_eff = np.zeros(num_bins_total) full_r_trans_eff = np.zeros(num_bins_total) full_z_eff = np.zeros(num_bins_total) - ii = np.arange(xcf.num_bins_r_par) + r_par_indices = np.arange(xcf.num_bins_r_par) r_trans = (0.5+np.arange(xcf.num_bins_r_trans))*xcf.r_trans_max/xcf.num_bins_r_trans for j in range(xcf.num_bins_r_trans) : - indices = j + xcf.num_bins_r_trans * ii - for k,i in zip(indices,ii) : - full_dmat[indices,k] = dmat[ii,i] + indices = j + xcf.num_bins_r_trans * r_par_indices + for k,i in zip(indices,r_par_indices) : + full_dmat[indices,k] = dmat[r_par_indices,i] full_r_par_eff[indices] = r_par_eff full_z_eff[indices] = z_eff full_r_trans_eff[indices] = r_trans[j] return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff -def main(cmdargs): +def main(): """Compute the metal matrix of the cross-correlation delta x object for a list of IGM absorption.""" parser = argparse.ArgumentParser( @@ -133,13 +145,15 @@ def main(cmdargs): type=str, default=None, required=True, - help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS containing table with at least rows "LOGLAM" and "WEIGHT"') + help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS' + ' containing table with at least rows "LOGLAM" and "WEIGHT"') parser.add_argument('--delta-dir', type=str, default=None, required=False, - help='Path to directory with delta*.gz to get the blinding info (default is trying to guess from attributes file)') + help='Path to directory with delta*.gz to get the blinding info' + ' (default is trying to guess from attributes file)') parser.add_argument('--drq', type=str, @@ -236,7 +250,7 @@ def main(cmdargs): default=10., required=False, help=('Limit the maximum redshift of the quasars ' - 'used as sources for spectra')) + 'used as sources for spectra')) parser.add_argument( '--lambda-abs', @@ -326,7 +340,7 @@ def main(cmdargs): required=False, help='Bins for the distribution of QSO redshifts') - args = parser.parse_args(cmdargs) + args = parser.parse_args() # setup variables in module xcf xcf.r_par_max = args.rp_max @@ -376,11 +390,11 @@ def main(cmdargs): zbins=args.qso_z_bins userprint(f"Use histogram of QSO redshifts with {zbins} bins") - hw,zbins = np.histogram(z_qso,bins=zbins,weights=weight_qso) - hwz,_ = np.histogram(z_qso,bins=zbins,weights=weight_qso*z_qso) - ok = (hw>0) - z_qso = hwz[ok]/hw[ok] # weighted mean in bins - weight_qso = hw[ok] + histo_w,zbins = np.histogram(z_qso,bins=zbins,weights=weight_qso) + histo_wz,_ = np.histogram(z_qso,bins=zbins,weights=weight_qso*z_qso) + selection = histo_w>0 + z_qso = histo_wz[selection]/histo_w[selection] # weighted mean in bins + weight_qso = histo_w[selection] t1 = time.time() userprint(f'picca_metal_xdmat.py - Time reading data: {(t1-t0)/60:.3f} minutes') @@ -532,5 +546,4 @@ def main(cmdargs): if __name__ == '__main__': - cmdargs=sys.argv[1:] - main(cmdargs) + main() From 8bce2254a0dce07e3dba4f06ec52419468c9ab02 Mon Sep 17 00:00:00 2001 From: Julien Guy Date: Tue, 12 Sep 2023 10:43:41 -0700 Subject: [PATCH 12/13] changed description --- bin/picca_fast_metal_dmat.py | 6 +----- bin/picca_fast_metal_xdmat.py | 6 +----- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index 09a6b6a85..9b9d79119 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -1,9 +1,5 @@ #!/usr/bin/env python -"""Compute the auto and cross-correlation of delta fields for a list of IGM -absorption. - -This module follow the procedure described in sections 4.3 of du Mas des -Bourboux et al. 2020 (In prep) to compute the distortion matrix +"""Compute the metal matrices """ import sys import os diff --git a/bin/picca_fast_metal_xdmat.py b/bin/picca_fast_metal_xdmat.py index 36145a78e..781608e8d 100755 --- a/bin/picca_fast_metal_xdmat.py +++ b/bin/picca_fast_metal_xdmat.py @@ -1,9 +1,5 @@ #!/usr/bin/env python -"""Compute the distortion matrix of the cross-correlation delta x object for a -list of IGM absorption. - -This module follow the procedure described in sections 4.3 of du Mas des -Bourboux et al. 2020 (In prep) to compute the distortion matrix +"""Compute the metal matrices """ import sys import os From d5da456d43c6b5790af5194ab2e78254281e7e40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ignasi=20P=C3=A9rez-R=C3=A0fols?= Date: Thu, 14 Sep 2023 14:26:33 +0200 Subject: [PATCH 13/13] yapdfed and linted --- bin/picca_fast_metal_dmat.py | 396 ++++++++++++++++++---------------- bin/picca_fast_metal_xdmat.py | 346 +++++++++++++++-------------- 2 files changed, 403 insertions(+), 339 deletions(-) diff --git a/bin/picca_fast_metal_dmat.py b/bin/picca_fast_metal_dmat.py index 9b9d79119..4be8de3ab 100755 --- a/bin/picca_fast_metal_dmat.py +++ b/bin/picca_fast_metal_dmat.py @@ -11,7 +11,8 @@ from picca import constants, cf, io from picca.utils import userprint -def read_stack_deltas_table(filename) : + +def read_stack_deltas_table(filename): """ Read stack. @@ -20,11 +21,16 @@ def read_stack_deltas_table(filename) : Returns: table as numpy.ndarray """ - return fitsio.read(filename,"STACK_DELTAS") + return fitsio.read(filename, "STACK_DELTAS") + -def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, - out_lambda_abs_1, out_lambda_abs_2, - stack_table_1, stack_table_2, rebin_factor=None): +def calc_fast_metal_dmat(in_lambda_abs_1, + in_lambda_abs_2, + out_lambda_abs_1, + out_lambda_abs_2, + stack_table_1, + stack_table_2, + rebin_factor=None): """Computes the metal distortion matrix. Args: @@ -51,34 +57,42 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, Note the global picca.cf contains the cosmology and the rp grid """ - loglam1=stack_table_1["LOGLAM"] - weight1=stack_table_1["WEIGHT"] - loglam2=stack_table_2["LOGLAM"] - weight2=stack_table_2["WEIGHT"] - if rebin_factor is not None : - size1=loglam1.size - loglam1 = loglam1[:(size1//rebin_factor)*rebin_factor].reshape((size1//rebin_factor),rebin_factor).mean(-1) - weight1 = weight1[:(size1//rebin_factor)*rebin_factor].reshape((size1//rebin_factor),rebin_factor).mean(-1) - size2=loglam2.size - loglam2 = loglam2[:(size2//rebin_factor)*rebin_factor].reshape((size2//rebin_factor),rebin_factor).mean(-1) - weight2 = weight2[:(size2//rebin_factor)*rebin_factor].reshape((size2//rebin_factor),rebin_factor).mean(-1) + loglam1 = stack_table_1["LOGLAM"] + weight1 = stack_table_1["WEIGHT"] + loglam2 = stack_table_2["LOGLAM"] + weight2 = stack_table_2["WEIGHT"] + if rebin_factor is not None: + size1 = loglam1.size + loglam1 = loglam1[:(size1 // rebin_factor) * rebin_factor].reshape( + (size1 // rebin_factor), rebin_factor).mean(-1) + weight1 = weight1[:(size1 // rebin_factor) * rebin_factor].reshape( + (size1 // rebin_factor), rebin_factor).mean(-1) + size2 = loglam2.size + loglam2 = loglam2[:(size2 // rebin_factor) * rebin_factor].reshape( + (size2 // rebin_factor), rebin_factor).mean(-1) + weight2 = weight2[:(size2 // rebin_factor) * rebin_factor].reshape( + (size2 // rebin_factor), rebin_factor).mean(-1) # input - input_z1 = (10**loglam1)/constants.ABSORBER_IGM[in_lambda_abs_1] - 1. - input_z2 = (10**loglam2)/constants.ABSORBER_IGM[in_lambda_abs_2] - 1. + input_z1 = (10**loglam1) / constants.ABSORBER_IGM[in_lambda_abs_1] - 1. + input_z2 = (10**loglam2) / constants.ABSORBER_IGM[in_lambda_abs_2] - 1. input_r1 = cf.cosmo.get_r_comov(input_z1) input_r2 = cf.cosmo.get_r_comov(input_z2) # all pairs - input_rp = (input_r1[:,None]-input_r2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) + input_rp = ( + input_r1[:, None] - + input_r2[None, :]).ravel() # same sign as line 676 of cf.py (1-2) if not cf.x_correlation: input_rp = np.abs(input_rp) # output - output_z1 = (10**loglam1)/constants.ABSORBER_IGM[out_lambda_abs_1] - 1. - output_z2 = (10**loglam2)/constants.ABSORBER_IGM[out_lambda_abs_2] - 1. + output_z1 = (10**loglam1) / constants.ABSORBER_IGM[out_lambda_abs_1] - 1. + output_z2 = (10**loglam2) / constants.ABSORBER_IGM[out_lambda_abs_2] - 1. output_r1 = cf.cosmo.get_r_comov(output_z1) output_r2 = cf.cosmo.get_r_comov(output_z2) # all pairs - output_rp = (output_r1[:,None]-output_r2[None,:]).ravel() # same sign as line 676 of cf.py (1-2) + output_rp = ( + output_r1[:, None] - + output_r2[None, :]).ravel() # same sign as line 676 of cf.py (1-2) if not cf.x_correlation: output_rp = np.abs(output_rp) @@ -88,51 +102,65 @@ def calc_fast_metal_dmat(in_lambda_abs_1, in_lambda_abs_2, # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and # used for the correlation function. It also has to be applied here. # we have them for both forests (1 and 2) - alpha_in_1 = cf.alpha_abs[in_lambda_abs_1] - alpha_out_1 = cf.alpha # = args.z_evol in main function - alpha_in_2 = cf.alpha_abs[in_lambda_abs_2] - alpha_out_2 = cf.alpha2 # = args.z_evol2 in main function + alpha_in_1 = cf.alpha_abs[in_lambda_abs_1] + alpha_out_1 = cf.alpha # = args.z_evol in main function + alpha_in_2 = cf.alpha_abs[in_lambda_abs_2] + alpha_out_2 = cf.alpha2 # = args.z_evol2 in main function # so here we have to apply both scalings (in the original code : # alpha_in is applied in cf.calc_metal_dmat and alpha_out in io.read_deltas) - weights = ((weight1*((1+input_z1)**(alpha_in_1+alpha_out_1-2)))[:,None] - *(weight2*((1+input_z2)**(alpha_in_2+alpha_out_2-2)))[None,:]).ravel() + weights = ( + (weight1 * ((1 + input_z1)**(alpha_in_1 + alpha_out_1 - 2)))[:, None] * + (weight2 * + ((1 + input_z2)**(alpha_in_2 + alpha_out_2 - 2)))[None, :]).ravel() # distortion matrix rpbins = cf.r_par_min \ + (cf.r_par_max-cf.r_par_min)/cf.num_bins_r_par*np.arange(cf.num_bins_r_par+1) # I checked the orientation of the matrix - dmat,_,_ = np.histogram2d(output_rp,input_rp,bins=(rpbins,rpbins),weights=weights) + dmat, _, _ = np.histogram2d(output_rp, + input_rp, + bins=(rpbins, rpbins), + weights=weights) # normalize (sum of weight should be one for each input rp,rt) - sum_in_weight,_ = np.histogram(input_rp,bins=rpbins,weights=weights) - dmat *= ((sum_in_weight>0)/(sum_in_weight+(sum_in_weight==0)))[None,:] + sum_in_weight, _ = np.histogram(input_rp, bins=rpbins, weights=weights) + dmat *= ((sum_in_weight > 0) / (sum_in_weight + + (sum_in_weight == 0)))[None, :] # mean outputs - sum_out_weight,_ = np.histogram(output_rp,bins=rpbins,weights=weights) - sum_out_weight_rp,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(output_rp[None,:].ravel())) - sum_out_weight_z,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(((output_z1[:,None]+output_z2[None,:])/2.).ravel())) - r_par_eff = sum_out_weight_rp/(sum_out_weight+(sum_out_weight==0)) - z_eff = sum_out_weight_z/(sum_out_weight+(sum_out_weight==0)) + sum_out_weight, _ = np.histogram(output_rp, bins=rpbins, weights=weights) + sum_out_weight_rp, _ = np.histogram(output_rp, + bins=rpbins, + weights=weights * + (output_rp[None, :].ravel())) + sum_out_weight_z, _ = np.histogram( + output_rp, + bins=rpbins, + weights=weights * + (((output_z1[:, None] + output_z2[None, :]) / 2.).ravel())) + r_par_eff = sum_out_weight_rp / (sum_out_weight + (sum_out_weight == 0)) + z_eff = sum_out_weight_z / (sum_out_weight + (sum_out_weight == 0)) # we could return the quantities computed as a function of rp only (and not rt): # return dmat, r_par_eff, r_trans_eff, z_eff # but for now we will return the full dmat to be consistent with the other computation # it consists in duplicating the result found to all rt, with output_rt = input_rt - num_bins_total = cf.num_bins_r_par*cf.num_bins_r_trans + num_bins_total = cf.num_bins_r_par * cf.num_bins_r_trans - full_dmat = np.zeros((num_bins_total,num_bins_total)) + full_dmat = np.zeros((num_bins_total, num_bins_total)) full_r_par_eff = np.zeros(num_bins_total) full_r_trans_eff = np.zeros(num_bins_total) - full_z_eff = np.zeros(num_bins_total) + full_z_eff = np.zeros(num_bins_total) r_par_indices = np.arange(cf.num_bins_r_par) - r_trans = (0.5+np.arange(cf.num_bins_r_trans))*cf.r_trans_max/cf.num_bins_r_trans - for j in range(cf.num_bins_r_trans) : - indices = j + cf.num_bins_r_trans * r_par_indices - for k,i in zip(indices,r_par_indices) : - full_dmat[indices,k] = dmat[r_par_indices,i] - full_r_par_eff[indices] = r_par_eff - full_z_eff[indices] = z_eff + r_trans = (0.5 + np.arange( + cf.num_bins_r_trans)) * cf.r_trans_max / cf.num_bins_r_trans + for j in range(cf.num_bins_r_trans): + indices = j + cf.num_bins_r_trans * r_par_indices + for k, i in zip(indices, r_par_indices): + full_dmat[indices, k] = dmat[r_par_indices, i] + full_r_par_eff[indices] = r_par_eff + full_z_eff[indices] = z_eff full_r_trans_eff[indices] = r_trans[j] return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff @@ -152,12 +180,14 @@ def main(): required=True, help='Output file name') - parser.add_argument('-i','--in-attributes', - type=str, - default=None, - required=True, - help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS' - ' containing table with at least rows "LOGLAM" and "WEIGHT"') + parser.add_argument( + '-i', + '--in-attributes', + type=str, + default=None, + required=True, + help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS' + ' containing table with at least rows "LOGLAM" and "WEIGHT"') parser.add_argument('--in-attributes2', type=str, @@ -165,12 +195,13 @@ def main(): required=False, help='Path to 2nd delta_attributes.fits.gz file') - parser.add_argument('--delta-dir', - type=str, - default=None, - required=False, - help='Path to directory with delta*.gz to get the blinding info' - ' (default is trying to guess from attributes file)') + parser.add_argument( + '--delta-dir', + type=str, + default=None, + required=False, + help='Path to directory with delta*.gz to get the blinding info' + ' (default is trying to guess from attributes file)') parser.add_argument('--rp-min', type=float, @@ -228,21 +259,19 @@ def main(): 'absorber redshift and the object redshift smaller than ' 'z-cut-max')) - parser.add_argument( - '--z-min-sources', - type=float, - default=0., - required=False, - help=('Limit the minimum redshift of the quasars ' - 'used as sources for spectra')) + parser.add_argument('--z-min-sources', + type=float, + default=0., + required=False, + help=('Limit the minimum redshift of the quasars ' + 'used as sources for spectra')) - parser.add_argument( - '--z-max-sources', - type=float, - default=10., - required=False, - help=('Limit the maximum redshift of the quasars ' - 'used as sources for spectra')) + parser.add_argument('--z-max-sources', + type=float, + default=10., + required=False, + help=('Limit the maximum redshift of the quasars ' + 'used as sources for spectra')) parser.add_argument( '--lambda-abs', @@ -337,14 +366,15 @@ def main(): action='store_true', required=False, help='rp can be positive or negative depending on the relative ' - 'position between absorber1 and absorber2') + 'position between absorber1 and absorber2') - parser.add_argument('--rebin-factor', - type=int, - default=None, - required=False, - help='Rebin factor for deltas. If not None, deltas will ' - 'be rebinned by that factor') + parser.add_argument( + '--rebin-factor', + type=int, + default=None, + required=False, + help='Rebin factor for deltas. If not None, deltas will ' + 'be rebinned by that factor') args = parser.parse_args() @@ -362,7 +392,7 @@ def main(): cf.alpha = args.z_evol cf.alpha2 = args.z_evol2 cf.lambda_abs = constants.ABSORBER_IGM[args.lambda_abs] - cf.x_correlation = False # I guess I have to specify this! + cf.x_correlation = False # I guess I have to specify this! cf.alpha_abs = {} cf.alpha_abs[args.lambda_abs] = cf.alpha @@ -370,38 +400,41 @@ def main(): cf.alpha_abs[metal] = args.metal_alpha # read blinding keyword - if args.delta_dir is None : - args.delta_dir = os.path.dirname(args.in_attributes)+"/../Delta/" - if not os.path.isdir(args.delta_dir) : - userprint(f"Tried to guess the delta directory (containing the delta*.gz files) but '{args.delta_dir}' is not valid") + if args.delta_dir is None: + args.delta_dir = os.path.dirname(args.in_attributes) + "/../Delta/" + if not os.path.isdir(args.delta_dir): + userprint( + "Tried to guess the delta directory (containing the delta*.gz " + f"files) but '{args.delta_dir}' is not valid" + ) userprint("Please specify the directory with option --delta-dir") sys.exit(1) blinding = io.read_blinding(args.delta_dir) # load fiducial cosmology cf.cosmo = constants.Cosmo(Om=args.fid_Om, - Or=args.fid_Or, - Ok=args.fid_Ok, - wl=args.fid_wl, - blinding=blinding) - - + Or=args.fid_Or, + Ok=args.fid_Ok, + wl=args.fid_wl, + blinding=blinding) t0 = time.time() ### Read data stack_table1 = read_stack_deltas_table(args.in_attributes) - if args.in_attributes2 is not None : + if args.in_attributes2 is not None: stack_table2 = read_stack_deltas_table(args.in_attributes2) - else : - stack_table2 = stack_table1 # reference to first one + else: + stack_table2 = stack_table1 # reference to first one t1 = time.time() - userprint(f'picca_fast_metal_dmat.py - Time reading data: {(t1-t0)/60:.3f} minutes') + userprint( + f'picca_fast_metal_dmat.py - Time reading data: {(t1-t0)/60:.3f} minutes' + ) abs_igm = [args.lambda_abs] + args.abs_igm - userprint("abs_igm = {}".format(abs_igm)) + userprint(f"abs_igm = {abs_igm}") if args.lambda_abs2 is None: args.lambda_abs2 = args.lambda_abs @@ -410,7 +443,7 @@ def main(): abs_igm_2 = [args.lambda_abs2] + args.abs_igm2 if cf.x_correlation: - userprint("abs_igm2 = {}".format(abs_igm_2)) + userprint(f"abs_igm2 = {abs_igm_2}") # intitialize arrays to store the results for the different metal absorption dmat_all = [] @@ -428,13 +461,17 @@ def main(): if index1 == 0 and index2 == 0: continue - userprint("Computing",abs_igm1,abs_igm2) + userprint("Computing", abs_igm1, abs_igm2) # this a matrix as a function of rp only - dmat, r_par_eff, r_trans_eff, z_eff = calc_fast_metal_dmat(abs_igm1,abs_igm2, - args.lambda_abs,args.lambda_abs2, - stack_table1, stack_table2, - rebin_factor = args.rebin_factor) + dmat, r_par_eff, r_trans_eff, z_eff = calc_fast_metal_dmat( + abs_igm1, + abs_igm2, + args.lambda_abs, + args.lambda_abs2, + stack_table1, + stack_table2, + rebin_factor=args.rebin_factor) # add these results to the list ofor the different metal absorption dmat_all.append(dmat) @@ -444,92 +481,83 @@ def main(): names.append(abs_igm1 + "_" + abs_igm2) t2 = time.time() - userprint(f'picca_fast_metal_dmat.py - Time computing all metal matrices : {(t2-t1)/60:.3f} minutes') + userprint( + f'picca_fast_metal_dmat.py - Time computing all metal matrices : {(t2-t1)/60:.3f} minutes' + ) # save the results results = fitsio.FITS(args.out, 'rw', clobber=True) - header = [ - { - 'name': 'RPMIN', - 'value': cf.r_par_min, - 'comment': 'Minimum r-parallel [h^-1 Mpc]' - }, - { - 'name': 'RPMAX', - 'value': cf.r_par_max, - 'comment': 'Maximum r-parallel [h^-1 Mpc]' - }, - { - 'name': 'RTMAX', - 'value': cf.r_trans_max, - 'comment': 'Maximum r-transverse [h^-1 Mpc]' - }, - { - 'name': 'NP', - 'value': cf.num_bins_r_par, - 'comment': 'Number of bins in r-parallel' - }, - { - 'name': 'NT', - 'value': cf.num_bins_r_trans, - 'comment': ' Number of bins in r-transverse' - }, - { - 'name': 'COEFMOD', - 'value': args.coef_binning_model, - 'comment': 'Coefficient for model binning' - }, - { - 'name': 'ZCUTMIN', - 'value': cf.z_cut_min, - 'comment': 'Minimum redshift of pairs' - }, - { - 'name': 'ZCUTMAX', - 'value': cf.z_cut_max, - 'comment': 'Maximum redshift of pairs' - }, - { - 'name': 'REJ', - 'value': cf.reject, - 'comment': 'Rejection factor' - }, - { - 'name': 'ALPHAMET', - 'value': args.metal_alpha, - 'comment': 'Evolution of metal bias' - }, { - 'name': 'OMEGAM', - 'value': args.fid_Om, - 'comment': 'Omega_matter(z=0) of fiducial LambdaCDM cosmology' - }, { - 'name': 'OMEGAR', - 'value': args.fid_Or, - 'comment': 'Omega_radiation(z=0) of fiducial LambdaCDM cosmology' - }, { - 'name': 'OMEGAK', - 'value': args.fid_Ok, - 'comment': 'Omega_k(z=0) of fiducial LambdaCDM cosmology' - }, { - 'name': 'WL', - 'value': args.fid_wl, - 'comment': 'Equation of state of dark energy of fiducial LambdaCDM cosmology' - }, { - 'name': "BLINDING", - 'value': blinding, - 'comment': 'String specifying the blinding strategy' - } - ] + header = [{ + 'name': 'RPMIN', + 'value': cf.r_par_min, + 'comment': 'Minimum r-parallel [h^-1 Mpc]' + }, { + 'name': 'RPMAX', + 'value': cf.r_par_max, + 'comment': 'Maximum r-parallel [h^-1 Mpc]' + }, { + 'name': 'RTMAX', + 'value': cf.r_trans_max, + 'comment': 'Maximum r-transverse [h^-1 Mpc]' + }, { + 'name': 'NP', + 'value': cf.num_bins_r_par, + 'comment': 'Number of bins in r-parallel' + }, { + 'name': 'NT', + 'value': cf.num_bins_r_trans, + 'comment': ' Number of bins in r-transverse' + }, { + 'name': 'COEFMOD', + 'value': args.coef_binning_model, + 'comment': 'Coefficient for model binning' + }, { + 'name': 'ZCUTMIN', + 'value': cf.z_cut_min, + 'comment': 'Minimum redshift of pairs' + }, { + 'name': 'ZCUTMAX', + 'value': cf.z_cut_max, + 'comment': 'Maximum redshift of pairs' + }, { + 'name': 'REJ', + 'value': cf.reject, + 'comment': 'Rejection factor' + }, { + 'name': 'ALPHAMET', + 'value': args.metal_alpha, + 'comment': 'Evolution of metal bias' + }, { + 'name': 'OMEGAM', + 'value': args.fid_Om, + 'comment': 'Omega_matter(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAR', + 'value': args.fid_Or, + 'comment': 'Omega_radiation(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAK', + 'value': args.fid_Ok, + 'comment': 'Omega_k(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': + 'WL', + 'value': + args.fid_wl, + 'comment': + 'Equation of state of dark energy of fiducial LambdaCDM cosmology' + }, { + 'name': "BLINDING", + 'value': blinding, + 'comment': 'String specifying the blinding strategy' + }] len_names = np.array([len(name) for name in names]).max() names = np.array(names, dtype='S' + str(len_names)) - results.write( - [ - np.array(names) - ], - names=['ABS_IGM'], - header=header, - comment=['Absorption name'], - extname='ATTRI') + results.write([np.array(names)], + names=['ABS_IGM'], + header=header, + comment=['Absorption name'], + extname='ATTRI') dmat_name = "DM_" if blinding != "none": @@ -573,7 +601,9 @@ def main(): results.close() t3 = time.time() - userprint(f'picca_fast_metal_dmat.py - Time total : {(t3-t0)/60:.3f} minutes') + userprint( + f'picca_fast_metal_dmat.py - Time total : {(t3-t0)/60:.3f} minutes') + if __name__ == '__main__': main() diff --git a/bin/picca_fast_metal_xdmat.py b/bin/picca_fast_metal_xdmat.py index 781608e8d..a933270ae 100755 --- a/bin/picca_fast_metal_xdmat.py +++ b/bin/picca_fast_metal_xdmat.py @@ -11,7 +11,8 @@ from picca import constants, xcf, io from picca.utils import userprint -def read_stack_deltas_table(filename) : + +def read_stack_deltas_table(filename): """ Read stack. @@ -20,10 +21,14 @@ def read_stack_deltas_table(filename) : Returns: table as numpy.ndarray """ - return fitsio.read(filename,"STACK_DELTAS") + return fitsio.read(filename, "STACK_DELTAS") + -def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, - z_qso, weight_qso, +def calc_fast_metal_dmat(in_lambda_abs, + out_lambda_abs, + stack_table, + z_qso, + weight_qso, rebin_factor=None): """Computes the metal distortion matrix. @@ -49,80 +54,106 @@ def calc_fast_metal_dmat(in_lambda_abs, out_lambda_abs, stack_table, Note the global picca.xcf contains the cosmology, the rp grid, and the QSO catalog """ - loglam=stack_table["LOGLAM"] - weight_forest=stack_table["WEIGHT"] - if rebin_factor is not None : - size=loglam.size - loglam = loglam[:(size//rebin_factor)*rebin_factor].reshape((size//rebin_factor),rebin_factor).mean(-1) - weight_forest = weight_forest[:(size//rebin_factor)*rebin_factor].reshape((size//rebin_factor),rebin_factor).mean(-1) + loglam = stack_table["LOGLAM"] + weight_forest = stack_table["WEIGHT"] + if rebin_factor is not None: + size = loglam.size + loglam = loglam[:(size // rebin_factor) * rebin_factor].reshape( + (size // rebin_factor), rebin_factor).mean(-1) + weight_forest = weight_forest[:(size // rebin_factor) * + rebin_factor].reshape( + (size // rebin_factor), + rebin_factor).mean(-1) # input - input_zf = (10**loglam)/constants.ABSORBER_IGM[in_lambda_abs] - 1. # redshift in the forest + input_zf = (10**loglam) / constants.ABSORBER_IGM[ + in_lambda_abs] - 1. # redshift in the forest input_rf = xcf.cosmo.get_r_comov(input_zf) - r_qso = xcf.cosmo.get_r_comov(z_qso) # comoving distance for qso + r_qso = xcf.cosmo.get_r_comov(z_qso) # comoving distance for qso # all pairs - input_rp = (input_rf[:,None]-r_qso[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) + input_rp = ( + input_rf[:, None] - + r_qso[None, :]).ravel() # same sign as line 528 of xcf.py (forest-qso) # output - output_zf = (10**loglam)/constants.ABSORBER_IGM[out_lambda_abs] - 1. + output_zf = (10**loglam) / constants.ABSORBER_IGM[out_lambda_abs] - 1. output_rf = xcf.cosmo.get_r_comov(output_zf) # all pairs - output_rp = (output_rf[:,None]-r_qso[None,:]).ravel() # same sign as line 528 of xcf.py (forest-qso) + output_rp = ( + output_rf[:, None] - + r_qso[None, :]).ravel() # same sign as line 528 of xcf.py (forest-qso) # weights # alpha_in: in (1+z)^(alpha_in-1) is a scaling used to model how the metal contribution # evolves with redshift (by default alpha=1 so that this has no effect) - alpha_in = xcf.alpha_abs[in_lambda_abs] + alpha_in = xcf.alpha_abs[in_lambda_abs] # alpha_out: (1+z)^(alpha_out-1) is applied to the delta weights in io.read_deltas and # used for the correlation function. It also has to be applied here. alpha_out = xcf.alpha_abs[out_lambda_abs] # so here we have to apply both scalings (in the original code : alpha_in is applied in # xcf.calc_metal_dmat and alpha_out in io.read_deltas) # qso weights have already been scaled with (1+z)^alpha_obj - weights = ((weight_forest*((1+input_zf)**(alpha_in+alpha_out-2)))[:,None]*weight_qso[None,:]).ravel() + weights = ((weight_forest * + ((1 + input_zf)**(alpha_in + alpha_out - 2)))[:, None] * + weight_qso[None, :]).ravel() # distortion matrix - rpbins = xcf.r_par_min + (xcf.r_par_max-xcf.r_par_min)/xcf.num_bins_r_par*np.arange(xcf.num_bins_r_par+1) + rpbins = xcf.r_par_min + ( + xcf.r_par_max - + xcf.r_par_min) / xcf.num_bins_r_par * np.arange(xcf.num_bins_r_par + 1) # I checked the orientation of the matrix - dmat,_,_ = np.histogram2d(output_rp,input_rp,bins=(rpbins,rpbins),weights=weights) + dmat, _, _ = np.histogram2d(output_rp, + input_rp, + bins=(rpbins, rpbins), + weights=weights) # normalize (sum of weight should be one for each input rp,rt) - sum_in_weight,_ = np.histogram(input_rp,bins=rpbins,weights=weights) - dmat *= ((sum_in_weight>0)/(sum_in_weight+(sum_in_weight==0)))[None,:] + sum_in_weight, _ = np.histogram(input_rp, bins=rpbins, weights=weights) + dmat *= ((sum_in_weight > 0) / (sum_in_weight + + (sum_in_weight == 0)))[None, :] # mean outputs - sum_out_weight,_ = np.histogram(output_rp,bins=rpbins,weights=weights) - sum_out_weight_rp,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(output_rp[None,:].ravel())) - sum_out_weight_z,_ = np.histogram(output_rp,bins=rpbins,weights=weights*(((output_zf[:,None]+z_qso[None,:])/2.).ravel())) - r_par_eff = sum_out_weight_rp/(sum_out_weight+(sum_out_weight==0)) - z_eff = sum_out_weight_z/(sum_out_weight+(sum_out_weight==0)) + sum_out_weight, _ = np.histogram(output_rp, bins=rpbins, weights=weights) + sum_out_weight_rp, _ = np.histogram(output_rp, + bins=rpbins, + weights=weights * + (output_rp[None, :].ravel())) + sum_out_weight_z, _ = np.histogram( + output_rp, + bins=rpbins, + weights=weights * + (((output_zf[:, None] + z_qso[None, :]) / 2.).ravel())) + r_par_eff = sum_out_weight_rp / (sum_out_weight + (sum_out_weight == 0)) + z_eff = sum_out_weight_z / (sum_out_weight + (sum_out_weight == 0)) # we could return the quantities computed as a function of rp only (and not rt): # return dmat, r_par_eff, r_trans_eff, z_eff # but for now we will return the full dmat to be consistent with the other computation # it consists in duplicating the result found to all rt, with output_rt = input_rt - num_bins_total = xcf.num_bins_r_par*xcf.num_bins_r_trans + num_bins_total = xcf.num_bins_r_par * xcf.num_bins_r_trans - full_dmat = np.zeros((num_bins_total,num_bins_total)) + full_dmat = np.zeros((num_bins_total, num_bins_total)) full_r_par_eff = np.zeros(num_bins_total) full_r_trans_eff = np.zeros(num_bins_total) - full_z_eff = np.zeros(num_bins_total) + full_z_eff = np.zeros(num_bins_total) r_par_indices = np.arange(xcf.num_bins_r_par) - r_trans = (0.5+np.arange(xcf.num_bins_r_trans))*xcf.r_trans_max/xcf.num_bins_r_trans - for j in range(xcf.num_bins_r_trans) : - indices = j + xcf.num_bins_r_trans * r_par_indices - for k,i in zip(indices,r_par_indices) : - full_dmat[indices,k] = dmat[r_par_indices,i] - full_r_par_eff[indices] = r_par_eff - full_z_eff[indices] = z_eff + r_trans = (0.5 + np.arange( + xcf.num_bins_r_trans)) * xcf.r_trans_max / xcf.num_bins_r_trans + for j in range(xcf.num_bins_r_trans): + indices = j + xcf.num_bins_r_trans * r_par_indices + for k, i in zip(indices, r_par_indices): + full_dmat[indices, k] = dmat[r_par_indices, i] + full_r_par_eff[indices] = r_par_eff + full_z_eff[indices] = z_eff full_r_trans_eff[indices] = r_trans[j] return full_dmat, full_r_par_eff, full_r_trans_eff, full_z_eff + def main(): """Compute the metal matrix of the cross-correlation delta x object for a list of IGM absorption.""" @@ -137,19 +168,22 @@ def main(): required=True, help='Output file name') - parser.add_argument('-i','--in-attributes', - type=str, - default=None, - required=True, - help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS' - ' containing table with at least rows "LOGLAM" and "WEIGHT"') + parser.add_argument( + '-i', + '--in-attributes', + type=str, + default=None, + required=True, + help='Path to delta_attributes.fits.gz file with hdu STACK_DELTAS' + ' containing table with at least rows "LOGLAM" and "WEIGHT"') - parser.add_argument('--delta-dir', - type=str, - default=None, - required=False, - help='Path to directory with delta*.gz to get the blinding info' - ' (default is trying to guess from attributes file)') + parser.add_argument( + '--delta-dir', + type=str, + default=None, + required=False, + help='Path to directory with delta*.gz to get the blinding info' + ' (default is trying to guess from attributes file)') parser.add_argument('--drq', type=str, @@ -160,7 +194,7 @@ def main(): parser.add_argument('--mode', type=str, default='sdss', - choices=['sdss','desi','desi_mocks','desi_healpix'], + choices=['sdss', 'desi', 'desi_mocks', 'desi_healpix'], required=False, help='type of catalog supplied, default sdss') @@ -232,21 +266,19 @@ def main(): 'absorber redshift and the object redshift smaller than ' 'z-cut-max')) - parser.add_argument( - '--z-min-sources', - type=float, - default=0., - required=False, - help=('Limit the minimum redshift of the quasars ' - 'used as sources for spectra')) + parser.add_argument('--z-min-sources', + type=float, + default=0., + required=False, + help=('Limit the minimum redshift of the quasars ' + 'used as sources for spectra')) - parser.add_argument( - '--z-max-sources', - type=float, - default=10., - required=False, - help=('Limit the maximum redshift of the quasars ' - 'used as sources for spectra')) + parser.add_argument('--z-max-sources', + type=float, + default=10., + required=False, + help=('Limit the maximum redshift of the quasars ' + 'used as sources for spectra')) parser.add_argument( '--lambda-abs', @@ -324,12 +356,13 @@ def main(): required=False, help='Equation of state of dark energy of fiducial LambdaCDM cosmology') - parser.add_argument('--rebin-factor', - type=int, - default=None, - required=False, - help='Rebin factor for deltas. If not None, deltas will ' - 'be rebinned by that factor') + parser.add_argument( + '--rebin-factor', + type=int, + default=None, + required=False, + help='Rebin factor for deltas. If not None, deltas will ' + 'be rebinned by that factor') parser.add_argument('--qso-z-bins', type=int, default=1000, @@ -357,10 +390,13 @@ def main(): xcf.alpha_abs[metal] = args.metal_alpha # read blinding keyword - if args.delta_dir is None : - args.delta_dir = os.path.dirname(args.in_attributes)+"/../Delta/" - if not os.path.isdir(args.delta_dir) : - userprint(f"Tried to guess the delta directory (containing the delta*.gz files) but '{args.delta_dir}' is not valid") + if args.delta_dir is None: + args.delta_dir = os.path.dirname(args.in_attributes) + "/../Delta/" + if not os.path.isdir(args.delta_dir): + userprint( + "Tried to guess the delta directory (containing the delta*.gz " + f"files) but '{args.delta_dir}' is not valid" + ) userprint("Please specify the directory with option --delta-dir") sys.exit(1) blinding = io.read_blinding(args.delta_dir) @@ -379,21 +415,25 @@ def main(): stack_table = read_stack_deltas_table(args.in_attributes) # read objets - catalog = io.read_drq(args.drq, z_min=args.z_min_obj, - z_max=args.z_max_obj, keep_bal=True, mode=args.mode) - z_qso = catalog['Z'] + catalog = io.read_drq(args.drq, + z_min=args.z_min_obj, + z_max=args.z_max_obj, + keep_bal=True, + mode=args.mode) + z_qso = catalog['Z'] weight_qso = ((1. + z_qso) / (1. + args.z_ref))**(args.z_evol_obj - 1.) - zbins=args.qso_z_bins + zbins = args.qso_z_bins userprint(f"Use histogram of QSO redshifts with {zbins} bins") - histo_w,zbins = np.histogram(z_qso,bins=zbins,weights=weight_qso) - histo_wz,_ = np.histogram(z_qso,bins=zbins,weights=weight_qso*z_qso) - selection = histo_w>0 - z_qso = histo_wz[selection]/histo_w[selection] # weighted mean in bins + histo_w, zbins = np.histogram(z_qso, bins=zbins, weights=weight_qso) + histo_wz, _ = np.histogram(z_qso, bins=zbins, weights=weight_qso * z_qso) + selection = histo_w > 0 + z_qso = histo_wz[selection] / histo_w[selection] # weighted mean in bins weight_qso = histo_w[selection] t1 = time.time() - userprint(f'picca_metal_xdmat.py - Time reading data: {(t1-t0)/60:.3f} minutes') + userprint( + f'picca_metal_xdmat.py - Time reading data: {(t1-t0)/60:.3f} minutes') # intitialize arrays to store the results for the different metal absorption dmat_all = [] @@ -405,15 +445,16 @@ def main(): # loop over metals for index, abs_igm in enumerate(args.abs_igm): - userprint("Computing",abs_igm) + userprint("Computing", abs_igm) # this a matrix as a function of rp only - dmat, r_par_eff, r_trans_eff, z_eff = calc_fast_metal_dmat(abs_igm, - args.lambda_abs, - stack_table, - z_qso, - weight_qso, - rebin_factor = args.rebin_factor) + dmat, r_par_eff, r_trans_eff, z_eff = calc_fast_metal_dmat( + abs_igm, + args.lambda_abs, + stack_table, + z_qso, + weight_qso, + rebin_factor=args.rebin_factor) # add these results to the list ofor the different metal absorption dmat_all.append(dmat) @@ -423,79 +464,72 @@ def main(): names.append(abs_igm) t2 = time.time() - userprint(f'picca_metal_xdmat.py - Time computing all metal matrix: {(t2-t1)/60:.3f} minutes') + userprint( + f'picca_metal_xdmat.py - Time computing all metal matrix: {(t2-t1)/60:.3f} minutes' + ) # save the results results = fitsio.FITS(args.out, 'rw', clobber=True) - header = [ - { - 'name': 'RPMIN', - 'value': xcf.r_par_min, - 'comment': 'Minimum r-parallel [h^-1 Mpc]' - }, - { - 'name': 'RPMAX', - 'value': xcf.r_par_max, - 'comment': 'Maximum r-parallel [h^-1 Mpc]' - }, - { - 'name': 'RTMAX', - 'value': xcf.r_trans_max, - 'comment': 'Maximum r-transverse [h^-1 Mpc]' - }, - { - 'name': 'NP', - 'value': xcf.num_bins_r_par, - 'comment': 'Number of bins in r-parallel' - }, - { - 'name': 'NT', - 'value': xcf.num_bins_r_trans, - 'comment': 'Number of bins in r-transverse' - }, - { - 'name': 'COEFMOD', - 'value': args.coef_binning_model, - 'comment': 'Coefficient for model binning' - }, - { - 'name': 'ZCUTMIN', - 'value': xcf.z_cut_min, - 'comment': 'Minimum redshift of pairs' - }, - { - 'name': 'ZCUTMAX', - 'value': xcf.z_cut_max, - 'comment': 'Maximum redshift of pairs' - }, - { - 'name': 'OMEGAM', - 'value': args.fid_Om, - 'comment': 'Omega_matter(z=0) of fiducial LambdaCDM cosmology' - }, { - 'name': 'OMEGAR', - 'value': args.fid_Or, - 'comment': 'Omega_radiation(z=0) of fiducial LambdaCDM cosmology' - }, { - 'name': 'OMEGAK', - 'value': args.fid_Ok, - 'comment': 'Omega_k(z=0) of fiducial LambdaCDM cosmology' - }, { - 'name': 'WL', - 'value': args.fid_wl, - 'comment': 'Equation of state of dark energy of fiducial LambdaCDM cosmology' - }, { - 'name': "BLINDING", - 'value': blinding, - 'comment': 'String specifying the blinding strategy' - } - ] + header = [{ + 'name': 'RPMIN', + 'value': xcf.r_par_min, + 'comment': 'Minimum r-parallel [h^-1 Mpc]' + }, { + 'name': 'RPMAX', + 'value': xcf.r_par_max, + 'comment': 'Maximum r-parallel [h^-1 Mpc]' + }, { + 'name': 'RTMAX', + 'value': xcf.r_trans_max, + 'comment': 'Maximum r-transverse [h^-1 Mpc]' + }, { + 'name': 'NP', + 'value': xcf.num_bins_r_par, + 'comment': 'Number of bins in r-parallel' + }, { + 'name': 'NT', + 'value': xcf.num_bins_r_trans, + 'comment': 'Number of bins in r-transverse' + }, { + 'name': 'COEFMOD', + 'value': args.coef_binning_model, + 'comment': 'Coefficient for model binning' + }, { + 'name': 'ZCUTMIN', + 'value': xcf.z_cut_min, + 'comment': 'Minimum redshift of pairs' + }, { + 'name': 'ZCUTMAX', + 'value': xcf.z_cut_max, + 'comment': 'Maximum redshift of pairs' + }, { + 'name': 'OMEGAM', + 'value': args.fid_Om, + 'comment': 'Omega_matter(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAR', + 'value': args.fid_Or, + 'comment': 'Omega_radiation(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': 'OMEGAK', + 'value': args.fid_Ok, + 'comment': 'Omega_k(z=0) of fiducial LambdaCDM cosmology' + }, { + 'name': + 'WL', + 'value': + args.fid_wl, + 'comment': + 'Equation of state of dark energy of fiducial LambdaCDM cosmology' + }, { + 'name': "BLINDING", + 'value': blinding, + 'comment': 'String specifying the blinding strategy' + }] len_names = np.array([len(name) for name in names]).max() names = np.array(names, dtype='S' + str(len_names)) results.write( - [ - np.array(names) - ], + [np.array(names)], names=['ABS_IGM'], header=header, comment=['Number of pairs', 'Number of used pairs', 'Absorption name'],