From fea36f2162dcbb87a671ea686a5ea8815fe772b0 Mon Sep 17 00:00:00 2001 From: Santiago Casas Date: Mon, 2 Dec 2024 00:39:41 +0100 Subject: [PATCH] updated photo likelih scripts --- .../likelihod_photometric-Debole-CAMB.py | 441 ++++++++++++++++++ notebooks/likelihod_photometric-Debole.py | 20 +- notebooks/likelihod_photometric.py | 21 +- requirements.txt | 33 +- 4 files changed, 484 insertions(+), 31 deletions(-) create mode 100644 notebooks/likelihod_photometric-Debole-CAMB.py diff --git a/notebooks/likelihod_photometric-Debole-CAMB.py b/notebooks/likelihod_photometric-Debole-CAMB.py new file mode 100644 index 0000000..7d0088b --- /dev/null +++ b/notebooks/likelihod_photometric-Debole-CAMB.py @@ -0,0 +1,441 @@ +#!/usr/bin/env python +# coding: utf-8 + +import numpy as np +import logging +from itertools import product +from copy import deepcopy, copy +from collections.abc import Sequence + +from cosmicfishpie.fishermatrix import cosmicfish +from cosmicfishpie.LSSsurvey import photo_obs as pobs +from cosmicfishpie.LSSsurvey import photo_cov as pcov +from cosmicfishpie.utilities.utils import printing as upr +from nautilus import Prior +from nautilus import Sampler +import re +import time +from pathlib import Path + + +def is_indexable_iterable(var): + return isinstance(var, (list, np.ndarray, Sequence)) and not isinstance(var, (str, bytes)) + + +logger = logging.getLogger("cosmicfishpie.cosmology.nuisance") +logger.setLevel(logging.INFO) + + +upr.debug = False + + +upr.debug_print("test") + + +outfolder = "nautichain_results" +outroot = "cosmicjellyfish_Euclid-DeboleR1-3x2photo_camb_withnuis" + +outpath = Path(outfolder) +outpath.mkdir(parents=True, exist_ok=True) # parents=True creates parent directories if needed +outfilepath = outpath / f"{outroot}" +outfilepath = str(outfilepath) +print(f"Results will be saved to {outfilepath}") + +sampler_settings = { + "n_live": 2000, + "n_networks": 16, + "n_batch": 256, + "pool": 64, +} + + +fiducial = { + "Omegam": 0.3145714273, + "Omegab": 0.0491989, + "h": 0.6737, + "ns": 0.96605, + "sigma8": 0.81, + "w0": -1.0, + "wa": 0.0, + "mnu": 0.0, + "num_massive_neutrinos" : 0, + "num_nu_massive": 0, + "num_nu_massless": 3.044, + "tau": 0.0561, + 'kmax': 10, + 'extrap_kmax': None +} +observables = ['WL', 'GCph'] + + +options = { + "accuracy": 1, + "feedback": 1, + "code": "camb", + "outroot": outroot, + "survey_name": "Euclid", + "survey_name_photo": "Euclid-Photometric-DeboleR1", + "survey_name_spectro": False, + "specs_dir": "./cosmicfishpie/configs/other_survey_specifications/", + "cosmo_model": "LCDM", +} +cosmoFM_fid = cosmicfish.FisherMatrix( + fiducialpars=fiducial, + options=options, + observables=observables, + cosmoModel=options["cosmo_model"], + surveyName=options["survey_name"], +) + + +cosmoFM_fid.IApars + + +photo_fid = pobs.ComputeCls(cosmopars=cosmoFM_fid.fiducialcosmopars, + photopars=cosmoFM_fid.photopars, + IApars=cosmoFM_fid.IApars, + biaspars=cosmoFM_fid.photobiaspars) + +photo_fid.compute_all() + +photo_cov_fid = pcov.PhotoCov(cosmopars=cosmoFM_fid.fiducialcosmopars, + photopars=cosmoFM_fid.photopars, + IApars=cosmoFM_fid.IApars, + biaspars=cosmoFM_fid.photobiaspars, + fiducial_Cls=photo_fid) + + +photo_cov_fid.allparsfid + + +print(photo_fid.binrange_WL) +print(photo_fid.binrange_GCph) +print(photo_cov_fid.ngalbin_WL) +print(photo_cov_fid.ngalbin_GCph) + + +def observable_Cell(photo_th: pobs.ComputeCls): + photo_th.compute_all() + binrange_GCph = photo_th.binrange_GCph + binrange_WL = photo_th.binrange_WL + nbin_GCph = len(binrange_GCph) + nbin_WL = len(binrange_WL) + + ells = photo_th.result["ells"] + output = dict(ells=ells) + + observables = photo_th.observables + if "WL" in observables: + Cell_LL = np.empty((len(ells), nbin_WL, nbin_WL), dtype=np.float64) + if "GCph" in observables: + Cell_GG = np.empty((len(ells), nbin_GCph, nbin_GCph), dtype=np.float64) + if "WL" in observables and "GCph" in observables: + Cell_GL = np.empty((len(ells), nbin_GCph, nbin_WL), dtype=np.float64) + + for i,j in product(binrange_WL, binrange_GCph): + + if "WL" in observables: + Cell_LL[:,i-1,j-1] = (photo_th.result["WL {}xWL {}".format(i,j)] + + np.eye(nbin_WL)[i-1,j-1] + * photo_cov_fid.ellipt_error**2.0 / photo_cov_fid.ngalbin_WL[i-1] + ) + + if "GCph" in observables: + Cell_GG[:,i-1,j-1] = (photo_th.result["GCph {}xGCph {}".format(i,j)] + + np.eye(nbin_GCph)[i-1,j-1] + * 1 / photo_cov_fid.ngalbin_GCph[i-1] + ) + + if "WL" in observables and "GCph" in observables: + Cell_GL[:,i-1,j-1] = photo_th.result["GCph {}xWL {}".format(i,j)] + + if "WL" in observables: + output["Cell_LL"] = Cell_LL + if "GCph" in observables: + output["Cell_GG"] = Cell_GG + if "WL" in observables and "GCph" in observables: + output["Cell_GL"] = Cell_GL + + return output + + +Cells_fid = observable_Cell(photo_fid) +print(Cells_fid.keys()) + + +print(Cells_fid["Cell_LL"].shape) +print(Cells_fid["Cell_GG"].shape) +print(Cells_fid["Cell_GL"].shape) + + +ellmax_WL = cosmoFM_fid.specs["lmax_WL"] +ellmax_GC = cosmoFM_fid.specs["lmax_GCph"] +ellmax_XC = np.minimum(ellmax_GC,ellmax_WL) +nbins_Glob = min(len(list(photo_fid.binrange_WL)), len(list(photo_fid.binrange_GCph))) +print(nbins_Glob) + + +def compute_chi2_per_obs(Cell_fid, Cell_th, ells, dells): + + dfid = np.linalg.det(Cell_fid) + dth = np.linalg.det(Cell_th) + + nells = len(ells) + _, _, nbin = Cell_fid.shape + + dmix = 0 + for i in range(nbin): + Cth_mix = copy(Cell_th) + Cth_mix[:,i,:] = Cell_fid[:,i,:] + dmix += np.linalg.det(Cth_mix) + + ingrd = ( + (2*ells+1) + *( + dmix[:nells]/dth[:nells] + + np.log(dth[:nells]/dfid[:nells]) + - nbin) + ) + ingrd = [*((ingrd[1:]+ingrd[:-1])/2 * dells[:-1]), ingrd[-1]*dells[-1]] + + chi2 = np.sum(ingrd) + return chi2 + + +def compute_chi2(Cells_fid, Cells_th): + """ + Compute χ² for wedges using fully vectorized operations. + Matches the loop implementation exactly. + + Parameters: + ---------- + Cells_fid: Dict + + Cells_th: Dict + + Returns: + ------- + float + χ² value + """ + chi2 = 0 + ells = Cells_fid["ells"] + + if "WL" in observables and not "GCph" in observables: + Cells_WL_th = Cells_th["Cell_LL"] + Cells_WL_fid = Cells_fid["Cell_LL"] + + iWL = np.searchsorted(ells, ellmax_WL) + ells_WL = np.insert(ells, iWL, ellmax_WL) + Dl_WL = np.diff(ells_WL)[:iWL] + ells_WL = ells_WL[:iWL] + + + chi2 += ( + photo_cov_fid.fsky_WL + * compute_chi2_per_obs(Cells_WL_fid, Cells_WL_th, ells_WL, Dl_WL) + ) + + if "GCph" in observables and not "WL" in observables: + Cells_GC_th = Cells_th["Cell_GG"] + Cells_GC_fid = Cells_fid["Cell_GG"] + + iGC = np.searchsorted(ells, ellmax_GC) + ells_GC = np.insert(ells, iGC, ellmax_GC) + Dl_GC = np.diff(ells_GC)[:iGC] + ells_GC = ells_GC[:iGC] + + chi2 += (photo_cov_fid.fsky_GCph + * compute_chi2_per_obs(Cells_GC_fid, Cells_GC_th, ells_GC, Dl_GC)) + + if "GCph" in observables and "WL" in observables: + Cells_XC_th = Cells_th["Cell_GL"] + Cells_XC_fid = Cells_fid["Cell_GL"] + Cells_GC_th = Cells_th["Cell_GG"] + Cells_GC_fid = Cells_fid["Cell_GG"] + Cells_WL_th = Cells_th["Cell_LL"] + Cells_WL_fid = Cells_fid["Cell_LL"] + + iGC = np.searchsorted(ells, ellmax_GC) + ells_GC = np.insert(ells, iGC, ellmax_GC) + Dl_GC = np.diff(ells_GC)[:iGC] + ells_GC = ells_GC[:iGC] + iWL = np.searchsorted(ells, ellmax_WL) + ells_WL = np.insert(ells, iWL, ellmax_WL) + Dl_WL = np.diff(ells_WL)[:iWL] + ells_WL = ells_WL[:iWL] + iXC = np.searchsorted(ells, ellmax_XC) + ells_XC = np.insert(ells, iXC, ellmax_XC) + Dl_XC = np.diff(ells_XC)[:iXC] + ells_XC = ells_GC[:iXC] + + big_th = np.block([[Cells_WL_th[:iXC],np.transpose(Cells_XC_th,(0,2,1))[:iXC]], + [Cells_XC_th[:iXC], Cells_GC_th[:iXC]]]) + big_fid = np.block([[Cells_WL_fid[:iXC],np.transpose(Cells_XC_fid,(0,2,1))[:iXC]], + [Cells_XC_fid[:iXC], Cells_GC_fid[:iXC]]]) + + chi2 += np.sqrt(photo_cov_fid.fsky_WL*photo_cov_fid.fsky_GCph) * compute_chi2_per_obs(big_fid, big_th,ells_XC, Dl_XC) + chi2 += photo_cov_fid.fsky_WL * compute_chi2_per_obs(Cells_WL_fid[:iXC], Cells_WL_th[:iXC],ells_WL[:iXC], Dl_WL[:iXC]) + + return chi2 + + +def loglike(param_vec, prior=None): + + if type(param_vec) == dict: + param_dict = deepcopy(param_vec) + elif is_indexable_iterable(param_vec) and prior is not None: + #print(f'Loading prior with keys: {prior.keys}') + param_dict={key: param_vec[i] for i, key in enumerate(prior.keys)} + + photopars = deepcopy(cosmoFM_fid.photopars) + for ii, pp in enumerate(cosmoFM_fid.photopars.keys()): + photopars[ii] = param_dict.pop(pp, cosmoFM_fid.photopars[pp]) + + photobiaspars = deepcopy(cosmoFM_fid.photobiaspars) + for ii, pp in enumerate(cosmoFM_fid.photobiaspars.keys()): + photobiaspars[pp] = param_dict.pop(pp, cosmoFM_fid.photobiaspars[pp]) + + IApars = deepcopy(cosmoFM_fid.IApars) + for ii, pp in enumerate(cosmoFM_fid.IApars.keys()): + IApars[pp] = param_dict.pop(pp, cosmoFM_fid.IApars[pp]) + + photo_vary = pobs.ComputeCls( + param_dict, + photopars, + IApars, + photobiaspars, + ) + Cells_th = observable_Cell(photo_vary) + + return -0.5 * compute_chi2(Cells_fid,Cells_th) + + +cosmoFM_fid.freeparams + + +cosmoFM_fid.allparams + + +samp1dic = { + 'Omegam': 0.3145714273, + 'Omegab': 0.0491989, + 'h': 0.6737, + 'ns': 0.96605, + 'sigma8': 0.81, + 'b1': 1.0997727037892875, + 'b2': 1.220245876862528, + 'b3': 1.2723993083933989, + 'b4': 1.316624471897739, + 'b5': 1.35812370570578, + 'b6': 1.3998214171814918, + } +print("Sample likelihood", loglike(samp1dic)) + + +loglike(photo_cov_fid.allparsfid) + + +photo_cov_fid.allparsfid + + +fishmat_photo = cosmoFM_fid.compute() + + +prior_nonuis = Prior() +prior_withnuis = Prior() + + +prior_dict ={ + 'Omegam': [0.24, 0.4], + 'Omegab': [0.04, 0.06], + 'h': [0.61, 0.75], + 'ns': [0.92, 1.00], + 'sigma8': [0.79, 0.83], + 'AIA': [1.0, 3.0], + 'etaIA' :[-6.0, 6.0], + 'b1': [1.0, 3.0], + 'b2': [1.0, 3.0], + 'b3': [1.0, 3.0], + 'b4': [1.0, 3.0], + 'b5': [1.0, 3.0], + 'b6': [1.0, 3.0], + 'b7': [1.0, 3.0], + 'b8': [1.0, 3.0], + 'b9': [1.0, 3.0], + 'b10': [1.0, 3.0] + } + + +cosmoFM_fid.freeparams + + +for par in prior_dict.keys(): + if par in cosmoFM_fid.freeparams.keys(): + dist_prior = (prior_dict[par][0], prior_dict[par][1]) + if re.match(r'b\d+', par): + prior_withnuis.add_parameter(par, dist_prior) + elif re.search(r'IA', par): + prior_withnuis.add_parameter(par, dist_prior) + else: + prior_nonuis.add_parameter(par, dist_prior) + prior_withnuis.add_parameter(par, dist_prior) + + +print(prior_nonuis.keys) +print(prior_withnuis.keys) + + +if "withnuis" in options["outroot"]: + prior_chosen = prior_withnuis +elif "nonuis" in options["outroot"]: + prior_chosen = prior_nonuis +else: + raise ValueError("No prior specified in the outroot") +print("Loading prior with keys: ", prior_chosen.keys) + + +tini = time.time() +print("Starting sampler at", time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(tini))) + + +sampler = Sampler(prior_chosen, + loglike, + n_live=sampler_settings["n_live"], + n_networks=sampler_settings["n_networks"], + n_batch=sampler_settings["n_batch"], + pool=sampler_settings["pool"], + pass_dict=False, + filepath=outfilepath+".hdf5", + resume=True, + likelihood_kwargs={'prior': prior_chosen} + ) +sampler.run(verbose=True, discard_exploration=True) +log_z_all = sampler.evidence() +print('Evidence:', log_z_all) +points_all, log_w_all, log_l_all = sampler.posterior() + + +tfin = time.time() +elapsed = tfin - tini +hours = int(elapsed // 3600) +minutes = int((elapsed % 3600) // 60) +seconds = int(elapsed % 60) +print("Sampler finished at", time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(tfin))) +print(f"Total time elapsed: {hours:02d}:{minutes:02d}:{seconds:02d}") + + +sample_wghlkl = (np.vstack((points_all.T, np.exp(log_w_all), log_l_all)).T) + + +outfile_chain = outfilepath+".txt" +print(f"Saving chain to text file {outfile_chain}") + + +headerlist = ['loglike', 'weights'] + list(prior_chosen.keys) +header = " ".join(headerlist) +print("Saving header: ", header) + + +np.savetxt(outfile_chain, sample_wghlkl, header=header) +print("Finished Sampling Run") diff --git a/notebooks/likelihod_photometric-Debole.py b/notebooks/likelihod_photometric-Debole.py index c1224a5..0c8735c 100644 --- a/notebooks/likelihod_photometric-Debole.py +++ b/notebooks/likelihod_photometric-Debole.py @@ -1,14 +1,12 @@ #!/usr/bin/env python # coding: utf-8 -import seaborn as sns import numpy as np import logging from itertools import product from copy import deepcopy, copy from collections.abc import Sequence - from cosmicfishpie.fishermatrix import cosmicfish from cosmicfishpie.LSSsurvey import photo_obs as pobs from cosmicfishpie.LSSsurvey import photo_cov as pcov @@ -17,9 +15,9 @@ from nautilus import Sampler import re import time +from pathlib import Path -snscolors = sns.color_palette("colorblind") def is_indexable_iterable(var): return isinstance(var, (list, np.ndarray, Sequence)) and not isinstance(var, (str, bytes)) @@ -34,14 +32,20 @@ def is_indexable_iterable(var): upr.debug_print("test") +outfolder = "nautichain_results" outroot = "cosmicjellyfish_Euclid-DeboleR1-3x2photo_symb_withnuis" +outpath = Path(outfolder) +outpath.mkdir(parents=True, exist_ok=True) # parents=True creates parent directories if needed +outfilepath = outpath / f"{outroot}" +outfilepath = str(outfilepath) +print(f"Results will be saved to {outfilepath}") sampler_settings = { "n_live": 2000, "n_networks": 16, "n_batch": 256, - "pool": 8, + "pool": 64, } @@ -67,7 +71,7 @@ def is_indexable_iterable(var): "survey_name": "Euclid", "survey_name_photo": "Euclid-Photometric-DeboleR1", "survey_name_spectro": False, - "specs_dir": "../cosmicfishpie/configs/other_survey_specifications/", + "specs_dir": "./cosmicfishpie/configs/other_survey_specifications/", "cosmo_model": "LCDM", } cosmoFM_fid = cosmicfish.FisherMatrix( @@ -415,7 +419,7 @@ def loglike(param_vec, prior=None): n_batch=sampler_settings["n_batch"], pool=sampler_settings["pool"], pass_dict=False, - filepath=options["outroot"]+".hdf5", + filepath=outfilepath+".hdf5", resume=True, likelihood_kwargs={'prior': prior_chosen} ) @@ -437,7 +441,7 @@ def loglike(param_vec, prior=None): sample_wghlkl = (np.vstack((points_all.T, np.exp(log_w_all), log_l_all)).T) -outfile_chain = options["outroot"]+".txt" +outfile_chain = outfilepath+".txt" print(f"Saving chain to text file {outfile_chain}") @@ -447,4 +451,4 @@ def loglike(param_vec, prior=None): np.savetxt(outfile_chain, sample_wghlkl, header=header) - +print("Finished Sampling Run") diff --git a/notebooks/likelihod_photometric.py b/notebooks/likelihod_photometric.py index 4c29bff..56536e9 100644 --- a/notebooks/likelihod_photometric.py +++ b/notebooks/likelihod_photometric.py @@ -1,7 +1,6 @@ #!/usr/bin/env python # coding: utf-8 -import seaborn as sns import numpy as np import logging from itertools import product @@ -17,9 +16,9 @@ from nautilus import Sampler import re import time +from pathlib import Path -snscolors = sns.color_palette("colorblind") def is_indexable_iterable(var): return isinstance(var, (list, np.ndarray, Sequence)) and not isinstance(var, (str, bytes)) @@ -33,15 +32,20 @@ def is_indexable_iterable(var): upr.debug_print("test") - +outfolder = "nautichain_results" outroot = "cosmicjellyfish_Euclid-ISTF-Pess-3x2photo_symb_withnuis" +outpath = Path(outfolder) +outpath.mkdir(parents=True, exist_ok=True) # parents=True creates parent directories if needed +outfilepath = outpath / f"{outroot}" +outfilepath = str(outfilepath) +print(f"Results will be saved to {outfilepath}") sampler_settings = { "n_live": 2000, "n_networks": 16, "n_batch": 256, - "pool": 8, + "pool": 64, } @@ -67,7 +71,8 @@ def is_indexable_iterable(var): "survey_name": "Euclid", "survey_name_photo": "Euclid-Photometric-ISTF-Pessimistic", "survey_name_spectro": False, - "specs_dir": "../cosmicfishpie/configs/default_survey_specifications/", + "specs_dir": "./cosmicfishpie/configs/default_survey_specifications/", + # relative to where script is launched from not where it is located "cosmo_model": "LCDM", } cosmoFM_fid = cosmicfish.FisherMatrix( @@ -415,7 +420,7 @@ def loglike(param_vec, prior=None): n_batch=sampler_settings["n_batch"], pool=sampler_settings["pool"], pass_dict=False, - filepath=options["outroot"]+".hdf5", + filepath=outfilepath+".hdf5", resume=True, likelihood_kwargs={'prior': prior_chosen} ) @@ -437,7 +442,7 @@ def loglike(param_vec, prior=None): sample_wghlkl = (np.vstack((points_all.T, np.exp(log_w_all), log_l_all)).T) -outfile_chain = options["outroot"]+".txt" +outfile_chain = outfilepath+".txt" print(f"Saving chain to text file {outfile_chain}") @@ -447,4 +452,4 @@ def loglike(param_vec, prior=None): np.savetxt(outfile_chain, sample_wghlkl, header=header) - +print("Finished Sampling Run") diff --git a/requirements.txt b/requirements.txt index 2047910..7f96418 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,16 +1,19 @@ -pip>=19.0 -astropy>=4.0 -jupyter>=1.0 -cython==3.0.11 -matplotlib>=3.1 -numpy>=1.17 -scipy>=1.3 -getdist>=1.1.3 -pyyaml -joblib +astropy==6.1.7 +camb==1.5.8 +Cython==3.0.11 +getdist==1.5.4 +h5py==3.12.1 +ipython==8.30.0 +joblib==1.4.2 +jupyter==1.1.1 +matplotlib==3.9.3 +nautilus-sampler==1.0.5 +numpy +pandas +PyYAML==6.0.2 +scikit-learn==1.5.2 +SciPy seaborn -# camb>=1.3 -# - pytest>=6.2 -# - pytest-cov>=2.10 -# - pytest-pycodestyle>=2.2.0 -# - pytest-pydocstyle>=2.2.0 +uv==0.5.5 +colossus +symbolic-pofk@git+https://github.com/DeaglanBartlett/symbolic_pofk.git