From 553f384381c6df0f61b13a4db25189959597f96d Mon Sep 17 00:00:00 2001 From: sambit-giri Date: Tue, 4 Apr 2023 15:29:07 +0200 Subject: [PATCH] modelling outside the training range of k and z --- .DS_Store | Bin 0 -> 6148 bytes .gitignore | 8 - BCemu.egg-info/PKG-INFO | 164 +++++++++++++ BCemu.egg-info/SOURCES.txt | 15 ++ BCemu.egg-info/dependency_links.txt | 1 + BCemu.egg-info/not-zip-safe | 1 + BCemu.egg-info/requires.txt | 8 + BCemu.egg-info/top_level.txt | 1 + build/lib/BCemu/BaryonEffectsEmulator.py | 292 +++++++++++++++++++++++ build/lib/BCemu/__init__.py | 19 ++ build/lib/BCemu/download.py | 28 +++ build/lib/BCemu/functions.py | 134 +++++++++++ build/lib/BCemu/input_data/README.rst | 1 + build/lib/BCemu/kpls.py | 61 +++++ dist/BCemu-1.0-py3.10.egg | Bin 0 -> 25789 bytes examples/test.py | 34 ++- src/BaryonEffectsEmulator.py | 57 ++++- tests/sample_BCparam.py | 93 ++++++++ 18 files changed, 896 insertions(+), 21 deletions(-) create mode 100644 .DS_Store delete mode 100644 .gitignore create mode 100644 BCemu.egg-info/PKG-INFO create mode 100644 BCemu.egg-info/SOURCES.txt create mode 100644 BCemu.egg-info/dependency_links.txt create mode 100644 BCemu.egg-info/not-zip-safe create mode 100644 BCemu.egg-info/requires.txt create mode 100644 BCemu.egg-info/top_level.txt create mode 100644 build/lib/BCemu/BaryonEffectsEmulator.py create mode 100644 build/lib/BCemu/__init__.py create mode 100644 build/lib/BCemu/download.py create mode 100644 build/lib/BCemu/functions.py create mode 100644 build/lib/BCemu/input_data/README.rst create mode 100644 build/lib/BCemu/kpls.py create mode 100644 dist/BCemu-1.0-py3.10.egg create mode 100644 tests/sample_BCparam.py diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..3a27c8b1f2a659bf256501a2fc49b013dfe8de69 GIT binary patch literal 6148 zcmeHK%}T>S5Z>*NO(;SS3Oz1(Em+%XEnY&bFJMFuDz&jigE1>fY7eE5v%Zi|;`2DO zyOq}3iw6~*f!S|%b|zuIgxwv+7; zqzwI;-S?!x;>K4}U${B26s4bI>h3{%ac7e}FtZoq@lMV7K5 z)S0s6r|WFC+BMN^w6V@klT0QrdlV&LC0fct|0 zMRW{i8r9YTgK7Z)@fxWFZ0sd4M;LSrW*T7ygzHp5oytv%!F4+Lg^6 + +The package also has a three parameter barynification model. Model A assumes all the three parameters to be independent of redshift while model B assumes the parameters to be redshift dependent via the following form: + +![](https://latex.codecogs.com/svg.latex?\inline&space;X(z)&space;=&space;X_0(1+z)^\nu). + +Below an example fit to the BAHAMAS simulation result is shown. + +```python +import numpy as np +import matplotlib.pyplot as plt +import BCemu +import pickle + +BAH = pickle.load(open('examples/BAHAMAS_data.pkl', 'rb')) + +bfcemu = BCemu.BCM_3param(Ob=0.0463, Om=0.2793) +bcmdict = {'log10Mc': 13.25, + 'thej' : 4.711, + 'deta' : 0.097} + +zs = [0,0.5] +k_eval = 10**np.linspace(-1,1.08,50) +p0_eval1 = bfcemu.get_boost(zs[0], bcmdict, k_eval) +p1_eval1 = bfcemu.get_boost(zs[1], bcmdict, k_eval) + +bfcemu = BCemu.BCM_3param(Ob=0.0463, Om=0.2793) +bcmdict = {'log10Mc': 13.25, + 'thej' : 4.711, + 'deta' : 0.097, + 'nu_Mc' : 0.038, + 'nu_thej': 0.0, + 'nu_deta': 0.060} + +zs = [0,0.5] +k_eval = 10**np.linspace(-1,1.08,50) +p0_eval2 = bfcemu.get_boost(zs[0], bcmdict, k_eval) +p1_eval2 = bfcemu.get_boost(zs[1], bcmdict, k_eval) + +plt.figure(figsize=(10,4.5)) +plt.subplot(121); plt.title('z=0') +plt.semilogx(BAH['z=0']['k'], BAH['z=0']['S'], '-', c='k', lw=5, alpha=0.2, label='BAHAMAS') +plt.semilogx(k_eval, p0_eval1, c='C0', lw=3, label='A', ls='--') +plt.semilogx(k_eval, p0_eval1, c='C2', lw=3, label='B', ls=':') +plt.axis([1e-1,12,0.73,1.04]) +plt.yticks(fontsize=14) +plt.xticks(fontsize=14) +plt.legend() +plt.xlabel(r'$k$ (h/Mpc)', fontsize=14) +plt.ylabel(r'$\frac{P_{\rm DM+baryon}}{P_{\rm DM}}$', fontsize=21) +plt.subplot(122); plt.title('z=0.5') +plt.semilogx(BAH['z=0.5']['k'], BAH['z=0.5']['S'], '-', c='k', lw=5, alpha=0.2, label='BAHAMAS') +plt.semilogx(k_eval, p1_eval1, c='C0', lw=3, label='A', ls='--') +plt.semilogx(k_eval, p1_eval2, c='C2', lw=3, label='B', ls=':') +plt.axis([1e-1,12,0.73,1.04]) +plt.yticks(fontsize=14) +plt.xticks(fontsize=14) +plt.xlabel(r'$k$ (h/Mpc)', fontsize=14) +plt.ylabel(r'$\frac{P_{\rm DM+baryon}}{P_{\rm DM}}$', fontsize=21) +plt.tight_layout() +plt.show() + + + +``` + + + + +## CONTRIBUTING + +If you find any bugs or unexpected behavior in the code, please feel free to open a [Github issue](https://github.com/sambit-giri/BCMemu/issues). The issue page is also good if you seek help or have suggestions for us. + +## References +[1] +Schneider, A., Teyssier, R., Stadel, J., Chisari, N. E., Le Brun, A. M., Amara, A., & Refregier, A. (2019). Quantifying baryon effects on the matter power spectrum and the weak lensing shear correlation. Journal of Cosmology and Astroparticle Physics, 2019(03), 020. [arXiv:1810.08629](https://arxiv.org/abs/1810.08629). + +[2] +Giri, S. K. & Schneider, A. (2021). Emulation of baryonic effects on the matter power spectrum and constraints from galaxy cluster data. Journal of Cosmology and Astroparticle Physics, 2021(12), 046. [arXiv:2108.08863](https://arxiv.org/abs/2108.08863). + diff --git a/BCemu.egg-info/SOURCES.txt b/BCemu.egg-info/SOURCES.txt new file mode 100644 index 0000000..4921f2f --- /dev/null +++ b/BCemu.egg-info/SOURCES.txt @@ -0,0 +1,15 @@ +LICENSE +README.md +setup.py +BCemu.egg-info/PKG-INFO +BCemu.egg-info/SOURCES.txt +BCemu.egg-info/dependency_links.txt +BCemu.egg-info/not-zip-safe +BCemu.egg-info/requires.txt +BCemu.egg-info/top_level.txt +src/BaryonEffectsEmulator.py +src/__init__.py +src/download.py +src/functions.py +src/kpls.py +src/input_data/README.rst \ No newline at end of file diff --git a/BCemu.egg-info/dependency_links.txt b/BCemu.egg-info/dependency_links.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/BCemu.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/BCemu.egg-info/not-zip-safe b/BCemu.egg-info/not-zip-safe new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/BCemu.egg-info/not-zip-safe @@ -0,0 +1 @@ + diff --git a/BCemu.egg-info/requires.txt b/BCemu.egg-info/requires.txt new file mode 100644 index 0000000..3eb46cd --- /dev/null +++ b/BCemu.egg-info/requires.txt @@ -0,0 +1,8 @@ +numpy +scipy +matplotlib +astropy +scikit-learn +smt==1.0.0 +cython +wget diff --git a/BCemu.egg-info/top_level.txt b/BCemu.egg-info/top_level.txt new file mode 100644 index 0000000..621ec32 --- /dev/null +++ b/BCemu.egg-info/top_level.txt @@ -0,0 +1 @@ +BCemu diff --git a/build/lib/BCemu/BaryonEffectsEmulator.py b/build/lib/BCemu/BaryonEffectsEmulator.py new file mode 100644 index 0000000..9c9097e --- /dev/null +++ b/build/lib/BCemu/BaryonEffectsEmulator.py @@ -0,0 +1,292 @@ +""" +Created by Sambit K. Giri +""" + +import numpy as np +from scipy import special +from scipy.interpolate import splev, splrep, bisplrep, bisplev +import os, pickle, pkg_resources +from glob import glob + +path_to_emu0_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z0_nComp3.pkl') +path_to_emu0p5_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z0p5_nComp3.pkl') +path_to_emu1_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z1_nComp3.pkl') +path_to_emu1p5_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z1p5_nComp3.pkl') +path_to_emu2_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z2_nComp3.pkl') + +ks_emulated = np.array([ 0.0341045 , 0.05861015, 0.08348237, 0.10855948, 0.13396836, + 0.15800331, 0.18254227, 0.20724802, 0.23212444, 0.2567891 , + 0.28103674, 0.30569611, 0.33079354, 0.35528021, 0.37965142, + 0.404148 , 0.42867984, 0.45328332, 0.4779432 , 0.50261366, + 0.5271016 , 0.55149765, 0.57618047, 0.60068947, 0.62521509, + 0.64995821, 0.67464717, 0.69909405, 0.72353962, 0.74794006, + 0.77257944, 0.79721068, 0.8218201 , 0.84641959, 0.87094783, + 0.89552692, 0.91995039, 0.94446984, 0.96905847, 0.99352717, + 1.01807451, 1.04265491, 1.06732874, 1.09198076, 1.11652679, + 1.14094596, 1.16538844, 1.18997249, 1.21454137, 1.23910819, + 1.26372463, 1.28832777, 1.31294019, 1.33736488, 1.36188676, + 1.38647099, 1.41088754, 1.43543816, 1.46000179, 1.48461337, + 1.50917878, 1.5336324 , 1.55828985, 1.58282587, 1.60737312, + 1.63194377, 1.65639122, 1.68101056, 1.70559414, 1.73006237, + 1.7545306 , 1.77905407, 1.80374906, 1.82836316, 1.85291022, + 1.87746965, 1.90194743, 1.9264312 , 1.95096646, 1.97544012, + 1.99995779, 2.02464604, 2.0492416 , 2.07372144, 2.09825877, + 2.12287081, 2.14742218, 2.17183929, 2.19641594, 2.22108696, + 2.24561567, 2.28249413, 2.31919863, 2.34372827, 2.36822366, + 2.39275106, 2.41736366, 2.45422376, 2.49100691, 2.51555959, + 2.5401016 , 2.57703286, 2.61373473, 2.63829902, 2.67525014, + 2.71184415, 2.73638728, 2.77337491, 2.81011536, 2.83469182, + 2.8716096 , 2.90829837, 2.94512739, 2.98188263, 3.01883393, + 3.05558917, 3.0923785 , 3.12909011, 3.16603785, 3.21513235, + 3.2519052 , 3.28875968, 3.32549358, 3.36239263, 3.41145867, + 3.44823639, 3.48510286, 3.5341483 , 3.5832604 , 3.62008455, + 3.65703267, 3.7060856 , 3.75515501, 3.80424913, 3.8532721 , + 3.902349 , 3.95147159, 4.00055598, 4.04966883, 4.09877157, + 4.14784735, 4.19693436, 4.24601388, 4.29511175, 4.34416169, + 4.39321611, 4.45468133, 4.5159563 , 4.56506615, 4.61413444, + 4.67561052, 4.73691582, 4.78599272, 4.84740625, 4.90865463, + 4.97012846, 5.03147234, 5.09286565, 5.15410242, 5.21558674, + 5.27687181, 5.33827449, 5.41191905, 5.47320637, 5.53460643, + 5.60825511, 5.68190978, 5.75552213, 5.82909965, 5.89035852, + 5.9518537 , 6.02549713, 6.09910125, 6.17273869, 6.25873037, + 6.34454978, 6.41815464, 6.49181793, 6.5654902 , 6.65144106, + 6.73724886, 6.8108717 , 6.89684204, 6.98270751, 7.06869507, + 7.15448714, 7.24045448, 7.33862776, 7.42446926, 7.51047555, + 7.60865781, 7.70682359, 7.79264862, 7.87861071, 7.97679073, + 8.07497898, 8.17312454, 8.27126411, 8.38183955, 8.49219854, + 8.59034035, 8.68854658, 8.79907933, 8.90941584, 9.01995908, + 9.13035252, 9.24083808, 9.35119033, 9.46172083, 9.58445502, + 9.70714802, 9.829853 , 9.94022996, 10.05076646, 10.17348567, + 10.30854414, 10.44348801, 10.56619074, 10.68889573, 10.82396168, + 10.95882092, 11.09389811, 11.2288637 , 11.36392292, 11.51118508, + 11.64610273, 11.78114023, 11.9283919 , 12.0756593 , 12.22293269, + 12.37017986, 12.51740232]) + + +def check_trained_emulators(): + if not ( len(glob(path_to_emu0_file))*len(glob(path_to_emu0p5_file))*len(glob(path_to_emu1_file))*len(glob(path_to_emu1p5_file))*len(glob(path_to_emu2_file)) ): + from . import download + download.download_emulators() + +def ps_suppression_8param(theta, emul, return_std=False): + log10Mc, mu, thej, gamma, delta, eta, deta, fb = theta + # fb = 0.1612 + mins = np.array([11, 0.0, 2, 1, 3, 0.05, 0.05, 0.10]) + maxs = np.array([15, 2.0, 8, 4, 11, 0.40, 0.40, 0.25]) + assert mins[0]<=log10Mc<=maxs[0] + assert mins[1]<=mu<=maxs[1] + assert mins[2]<=thej<=maxs[2] + assert mins[3]<=gamma<=maxs[3] + assert mins[4]<=delta<=maxs[4] + assert mins[5]<=eta<=maxs[5] + assert mins[6]<=deta<=maxs[6] + assert mins[7]<=fb<=maxs[7] + theta = [log10Mc, mu, thej, gamma, delta, eta, deta, fb] + if type(theta)==list: theta = np.array(theta) + if theta.ndim==1: theta = theta[None,:] + out = emul.predict_values(theta)#, return_std=True) + # print(out.shape) + return out.squeeze() + +def nearest_element_idx(arr, a, both=True): + if both: + dist = np.abs(arr-a) + dist_arg = np.argsort(dist) + return dist_arg[0], dist_arg[1] + else: + return np.abs(arr-a).argmin() + +def clip_range(x, mn, mx, message=None): + if xmx: out = mx + else: out = x + if out!=x: + if message is not None: print('{} set from {:.3f} to {:.3f}'.format(message,x,out)) + return out + +class use_emul: + def __init__(self, emul_names=None, Ob=0.0463, Om=0.2793, verbose=True): + if emul_names is None: + emul_names = {'0' : path_to_emu0_file, + '0.5': path_to_emu0p5_file, + '1' : path_to_emu1_file, + '1.5': path_to_emu1p5_file, + '2' : path_to_emu2_file + } + self.emul_names = emul_names + self.verbose = verbose + check_trained_emulators() + self.load_emulators() + self.update_cosmology(Ob, Om) + self.ks0 = ks_emulated + self.__dict__.update({'nu_Mc': 0, 'nu_mu': 0, 'nu_thej': 0, 'nu_gamma': 0, 'nu_delta': 0, 'nu_eta': 0, 'nu_deta': 0}) + + def load_emulators(self, emul_names=None): + if emul_names is not None: self.emul_names = emul_names + emulators = [] + zs = [] + for ke in self.emul_names: + zs.append(float(ke)) + emu = pickle.load(open(self.emul_names[ke],'rb')) + emu.options['print_prediction'] = False + emulators.append(emu) + print('Emulators loaded.') + self.emulators = np.array(emulators) + self.emul_zs = np.array(zs) + + def update_cosmology(self, Ob, Om): + self.fb = Ob/Om + if self.verbose: print('Baryon fraction is set to {:.3f}'.format(self.fb)) + + def check_range(self): + mins = [11, 0.0, 2, 1, 3, 0.05, 0.05, 0.10] + maxs = [15, 2.0, 8, 4, 11, 0.40, 0.40, 0.25] + self.mins = mins + self.maxs = maxs + knob = False + if mins[0]<=self.log10Mc<=maxs[0] and mins[1]<=self.mu<=maxs[1] and mins[2]<=self.thej<=maxs[2] and mins[3]<=self.gamma<=maxs[3] and mins[4]<=self.delta<=maxs[4] and mins[5]<=self.eta<=maxs[5] and mins[6]<=self.deta<=maxs[6] and mins[7]<=self.fb<=maxs[7]: + knob = True + if not knob: + print('The parameters provided are outside the allowed range.') + print('log10Mc : [{},{}]'.format(self.mins[0],self.maxs[0])) + print('mu : [{},{}]'.format(self.mins[1],self.maxs[1])) + print('thej : [{},{}]'.format(self.mins[2],self.maxs[2])) + print('gamma : [{},{}]'.format(self.mins[3],self.maxs[3])) + print('delta : [{},{}]'.format(self.mins[4],self.maxs[4])) + print('eta : [{},{}]'.format(self.mins[5],self.maxs[5])) + print('deta : [{},{}]'.format(self.mins[6],self.maxs[6])) + print('fb : [{},{}]'.format(self.mins[7],self.maxs[7])) + assert knob + return None + + def z_evolve_param(self, z): + self.log10Mc = clip_range(self.log10Mc*(1+z)**(-self.nu_Mc), self.mins[0], self.maxs[0], message='log10Mc' if self.verbose else None) + self.mu = clip_range(self.mu*(1+z)**(-self.nu_mu), self.mins[1], self.maxs[1], message='mu' if self.verbose else None) + self.thej = clip_range(self.thej*(1+z)**(-self.nu_thej), self.mins[2], self.maxs[2], message='thej' if self.verbose else None) + self.gamma = clip_range(self.gamma*(1+z)**(-self.nu_gamma), self.mins[3], self.maxs[3], message='gamma' if self.verbose else None) + self.delta = clip_range(self.delta*(1+z)**(-self.nu_delta), self.mins[4], self.maxs[4], message='delta' if self.verbose else None) + self.eta = clip_range(self.eta*(1+z)**(-self.nu_eta), self.mins[5], self.maxs[5], message='eta' if self.verbose else None) + self.deta = clip_range(self.deta*(1+z)**(-self.nu_deta), self.mins[6], self.maxs[6], message='deta' if self.verbose else None) + if self.verbose: + if not self.log10Mc or not self.mu or not self.thej or not self.gamma or not self.delta or not self.eta or not self.deta: + if not z: + print('Parameters evolved.') + return None + + def make_theta(self, z): + self.z_evolve_param(z) + theta = [self.log10Mc, self.mu, self.thej, self.gamma, self.delta, self.eta, self.deta, self.fb] + return theta + + def run(self, BCM_dict, z=0, nu_dict=None, return_std=False): + assert 0<=z<=2 + self.BCM_dict = BCM_dict + + self.__dict__.update(BCM_dict) + if nu_dict is not None: self.__dict__.update(nu_dict) + + self.check_range() + + if z in self.emul_zs: + theta = self.make_theta(z) + emu0 = self.emulators[self.emul_zs==z][0] + ps = ps_suppression_8param(theta, emu0, return_std=return_std) + if return_std: return ps[0], self.ks0, ps[1] + else: return ps, self.ks0 + else: + i0, i1 = nearest_element_idx(self.emul_zs, z) + theta0 = self.make_theta(self.emul_zs[i0]) + emu0 = self.emulators[i0] + theta1 = self.make_theta(self.emul_zs[i1]) + emu1 = self.emulators[i1] + ps0 = ps_suppression_8param(theta0, emu0, return_std=return_std) + ps1 = ps_suppression_8param(theta1, emu1, return_std=return_std) + if return_std: return ps0[0] + (ps1[0]-ps0[0])*(z-self.emul_zs[i0])/(self.emul_zs[i1]-self.emul_zs[i0]), self.ks0, ps0[1] + (ps1[1]-ps0[1])*(z-self.emul_zs[i0])/(self.emul_zs[i1]-self.emul_zs[i0]) + else: return ps0 + (ps1-ps0)*(z-self.emul_zs[i0])/(self.emul_zs[i1]-self.emul_zs[i0]), self.ks0 + + +class BCM_7param(use_emul): + def __init__(self, emul_names=None, Ob=0.0463, Om=0.2793, verbose=True, + below_kmin='extrapolate', above_kmax='extrapolate', + above_zmax='extrapolate'): + super().__init__(emul_names=emul_names, Ob=Ob, Om=Om, verbose=verbose) + self.below_kmin = below_kmin + self.above_kmax = above_kmax + self.above_zmax = above_zmax + + def print_param_names(self): + print('\nBaryonification parameters:') + print('--------------------------') + print('log10Mc, mu, thej, gamma, delta, eta, deta\n') + print('Redshift evolution parameters:') + print('-----------------------------') + print('nu_Mc, nu_mu, nu_thej, nu_gamma, nu_delta, nu_eta, nu_deta\n') + + def get_boost(self, z, BCM_params, k_eval, fb=None): + if fb is not None: self.fb = fb + if z in [0,0.5,1,1.5,2]: + pp, kk = self.run(BCM_params, z=z) + else: + p0p0, kk = self.run(BCM_params, z=0) + p0p5, kk = self.run(BCM_params, z=0.5) + p1p0, kk = self.run(BCM_params, z=1) + p1p5, kk = self.run(BCM_params, z=1.5) + p2p0, kk = self.run(BCM_params, z=2) + P = np.array([p0p0,p0p5,p1p5,p1p5,p2p0]) + pp = np.array([splev(z,splrep([0,0.5,1,1.5,2], P[:,ii], s=0, k=1)) + for ii,ki in enumerate(kk)]) + if z>2: + if str(self.above_zmax).lower()=='zmax': + print('Redshift z>2 is outside the training set and therefore set to suppressions at z=2.') + pp = p2p0 + else: + print('Redshift z>2 is outside the training set and therefore extrapolated.') + + pp_tck = splrep(kk, pp, k=3) + p_eval = splev(k_eval, pp_tck, ext=0) + if k_eval.min()self.ks0.max(): + if str(self.above_kmax).lower()=='extrapolate': + print('Suppressions modelled at k>{:.3f} h/Mpc are outside the training set and therefore extrapolated.'.format(self.ks0.max())) + else: + print('Suppressions modelled at k>{:.3f} h/Mpc are outside the training set.'.format(self.ks0.max())) + try: p_eval[k_eval>self.ks0.max()] = self.above_kmax(k_eval[k_eval>self.ks0.max()]) + except: p_eval[k_eval>self.ks0.max()] = self.above_kmax + return p_eval + + +class BCM_3param(use_emul): + def __init__(self, emul_names=None, Ob=0.0463, Om=0.2793, verbose=True, + below_kmin='extrapolate', above_kmax='extrapolate'): + super().__init__(emul_names=emul_names, Ob=Ob, Om=Om, verbose=verbose, + below_kmin=below_kmin, above_kmax=above_kmax) + + def print_param_names(self): + print('\nBaryonification parameters:') + print('--------------------------') + print('log10Mc, thej, deta\n') + print('Redshift evolution parameters:') + print('-----------------------------') + print('nu_Mc, nu_thej, nu_deta\n') + + def get_boost(self, z, BCM_params, k_eval, fb=None): + BCM_params['delta'] = 7.0 + BCM_params['eta'] = 0.2 + BCM_params['mu'] = 1.0 + BCM_params['gamma'] = 2.5 + if k_eval.min()self.ks0.max(): + print('k values above {:.3f} h/Mpc are errornous.'.format(self.ks0.max())) + if fb is not None: self.fb = fb + pp, kk = self.run(BCM_params, z=z) + pp_tck = splrep(kk, pp, k=3) + return splev(k_eval, pp_tck, ext=0) + + diff --git a/build/lib/BCemu/__init__.py b/build/lib/BCemu/__init__.py new file mode 100644 index 0000000..bdd5470 --- /dev/null +++ b/build/lib/BCemu/__init__.py @@ -0,0 +1,19 @@ +''' +Emulator is a Python package for constructing emulators. + +You can also get documentation for all routines directory from +the interpreter using Python's built-in help() function. +For example: +>>> import BCemu +>>> help(BCemu.use_emul) +''' +import sys +from .BaryonEffectsEmulator import * +from . import download +from . import kpls +import smt + + +#Suppress warnings from zero-divisions and nans +import numpy +numpy.seterr(all='ignore') diff --git a/build/lib/BCemu/download.py b/build/lib/BCemu/download.py new file mode 100644 index 0000000..7bfad0f --- /dev/null +++ b/build/lib/BCemu/download.py @@ -0,0 +1,28 @@ +import os, time +import wget, pkg_resources + +def download_emulators(): + print('Downloading emulators...') + + path_to_emul_folder = pkg_resources.resource_filename('BCemu', 'input_data/') + + path_to_emu0_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z0_nComp3.pkl') + path_to_emu0p5_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z0p5_nComp3.pkl') + path_to_emu1_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z1_nComp3.pkl') + path_to_emu1p5_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z1p5_nComp3.pkl') + path_to_emu2_file = pkg_resources.resource_filename('BCemu', 'input_data/kpls_emulator_z2_nComp3.pkl') + + link_to_emu0_file = 'https://github.com/sambit-giri/BCemu/releases/download/v1.0/kpls_emulator_z0_nComp3.pkl' + link_to_emu0p5_file = 'https://github.com/sambit-giri/BCemu/releases/download/v1.0/kpls_emulator_z0p5_nComp3.pkl' + link_to_emu1_file = 'https://github.com/sambit-giri/BCemu/releases/download/v1.0/kpls_emulator_z1_nComp3.pkl' + link_to_emu1p5_file = 'https://github.com/sambit-giri/BCemu/releases/download/v1.0/kpls_emulator_z1p5_nComp3.pkl' + link_to_emu2_file = 'https://github.com/sambit-giri/BCemu/releases/download/v1.0/kpls_emulator_z2_nComp3.pkl' + + print('at z = 0.0'); wget.download(link_to_emu0_file, path_to_emul_folder) + print('\nat z = 0.5'); wget.download(link_to_emu0p5_file, path_to_emul_folder) + print('\nat z = 1.0'); wget.download(link_to_emu1_file, path_to_emul_folder) + print('\nat z = 1.5'); wget.download(link_to_emu1p5_file, path_to_emul_folder) + print('\nat z = 2.0'); wget.download(link_to_emu2_file, path_to_emul_folder) + + print('...done') + return None \ No newline at end of file diff --git a/build/lib/BCemu/functions.py b/build/lib/BCemu/functions.py new file mode 100644 index 0000000..6b93d85 --- /dev/null +++ b/build/lib/BCemu/functions.py @@ -0,0 +1,134 @@ +import numpy as np +import pickle + +def ps_suppression_8param(theta, emul, return_std=False): + log10Mc, mu, thej, gamma, delta, eta, deta, fb = theta + # fb = 0.1612 + mins = [11, 0.0, 2, 1, 3, 0.2, 0.2, 0.10] + maxs = [15, 2.0, 8, 4, 11, 0.4, 0.4, 0.25] + assert mins[0]<=log10Mc<=maxs[0] + assert mins[1]<=mu<=maxs[1] + assert mins[2]<=thej<=maxs[2] + assert mins[3]<=gamma<=maxs[3] + assert mins[4]<=delta<=maxs[4] + assert mins[5]<=eta<=maxs[5] + assert mins[6]<=deta<=maxs[6] + assert mins[7]<=fb<=maxs[7] + theta = [log10Mc, mu, thej, gamma, delta, eta, deta, fb] + if type(theta)==list: theta = np.array(theta) + if theta.ndim==1: theta = theta[None,:] + out = emul.predict_values(theta)#, return_std=True) + # print(out.shape) + return out.squeeze() + +class use_emul: + def __init__(self, emul_names, Ob, Om): + self.emul_names = emul_names + self.fb = Ob/Om + print('Baryon fraction is set to {:.3f}'.format(self.fb)) + self.load_emulators() + self.fix_params() + + def load_emulators(self, emul_names=None): + if emul_names is not None: self.emul_names = emul_names + emulators = [] + zs = [] + for ke in self.emul_names: + zs.append(float(ke)) + emu = pickle.load(open(self.emul_names[ke],'rb')) + emu.options['print_prediction'] = False + emulators.append(emu) + print('Emulators loaded.') + self.emulators = np.array(emulators) + self.emul_zs = np.array(zs) + + def fix_params(self): + mu, gamma, delta, deta = 0.93, 2.25, 6.40, 0.240 + thej, eta = 4.235, 0.22 + self.mu = mu + self.gamma = gamma + self.delta = delta + self.deta = deta + self.thej = thej + self.eta = eta + + def run(self, log10Mc=13.322, nu_Mc=-0.015, z=0): + assert 12.3<=log10Mc<=14.5 + assert -0.1<=nu_Mc<=0.01 + assert 0<=z<=2 + + if z in self.emul_zs: + theta = [log10Mc*(1+z)**nu_Mc, self.mu, self.thej, self.gamma, self.delta, self.eta, self.deta, self.fb] + emu0 = self.emulators[self.emul_zs==z][0] + # print(emu0) + ps = ps_suppression_8param(theta, emu0, return_std=False) + return ps, ks0 + else: + i0, i1 = nearest_element_idx(self.emul_zs, z) + theta0 = [log10Mc*(1+self.emul_zs[i0])**nu_Mc, self.mu, self.thej, self.gamma, self.delta, self.eta, self.deta, self.fb] + emu0 = self.emulators[i0] + theta1 = [log10Mc*(1+self.emul_zs[i1])**nu_Mc, self.mu, self.thej, self.gamma, self.delta, self.eta, self.deta, self.fb] + emu1 = self.emulators[i1] + ps0 = ps_suppression_8param(theta0, emu0, return_std=False) + ps1 = ps_suppression_8param(theta1, emu1, return_std=False) + return ps0 + (ps1-ps0)*(z-self.emul_zs[i0])/(self.emul_zs[i1]-self.emul_zs[i0]), ks0 + + +ks0 = np.array([ 0.0341045 , 0.05861015, 0.08348237, 0.10855948, 0.13396836, + 0.15800331, 0.18254227, 0.20724802, 0.23212444, 0.2567891 , + 0.28103674, 0.30569611, 0.33079354, 0.35528021, 0.37965142, + 0.404148 , 0.42867984, 0.45328332, 0.4779432 , 0.50261366, + 0.5271016 , 0.55149765, 0.57618047, 0.60068947, 0.62521509, + 0.64995821, 0.67464717, 0.69909405, 0.72353962, 0.74794006, + 0.77257944, 0.79721068, 0.8218201 , 0.84641959, 0.87094783, + 0.89552692, 0.91995039, 0.94446984, 0.96905847, 0.99352717, + 1.01807451, 1.04265491, 1.06732874, 1.09198076, 1.11652679, + 1.14094596, 1.16538844, 1.18997249, 1.21454137, 1.23910819, + 1.26372463, 1.28832777, 1.31294019, 1.33736488, 1.36188676, + 1.38647099, 1.41088754, 1.43543816, 1.46000179, 1.48461337, + 1.50917878, 1.5336324 , 1.55828985, 1.58282587, 1.60737312, + 1.63194377, 1.65639122, 1.68101056, 1.70559414, 1.73006237, + 1.7545306 , 1.77905407, 1.80374906, 1.82836316, 1.85291022, + 1.87746965, 1.90194743, 1.9264312 , 1.95096646, 1.97544012, + 1.99995779, 2.02464604, 2.0492416 , 2.07372144, 2.09825877, + 2.12287081, 2.14742218, 2.17183929, 2.19641594, 2.22108696, + 2.24561567, 2.28249413, 2.31919863, 2.34372827, 2.36822366, + 2.39275106, 2.41736366, 2.45422376, 2.49100691, 2.51555959, + 2.5401016 , 2.57703286, 2.61373473, 2.63829902, 2.67525014, + 2.71184415, 2.73638728, 2.77337491, 2.81011536, 2.83469182, + 2.8716096 , 2.90829837, 2.94512739, 2.98188263, 3.01883393, + 3.05558917, 3.0923785 , 3.12909011, 3.16603785, 3.21513235, + 3.2519052 , 3.28875968, 3.32549358, 3.36239263, 3.41145867, + 3.44823639, 3.48510286, 3.5341483 , 3.5832604 , 3.62008455, + 3.65703267, 3.7060856 , 3.75515501, 3.80424913, 3.8532721 , + 3.902349 , 3.95147159, 4.00055598, 4.04966883, 4.09877157, + 4.14784735, 4.19693436, 4.24601388, 4.29511175, 4.34416169, + 4.39321611, 4.45468133, 4.5159563 , 4.56506615, 4.61413444, + 4.67561052, 4.73691582, 4.78599272, 4.84740625, 4.90865463, + 4.97012846, 5.03147234, 5.09286565, 5.15410242, 5.21558674, + 5.27687181, 5.33827449, 5.41191905, 5.47320637, 5.53460643, + 5.60825511, 5.68190978, 5.75552213, 5.82909965, 5.89035852, + 5.9518537 , 6.02549713, 6.09910125, 6.17273869, 6.25873037, + 6.34454978, 6.41815464, 6.49181793, 6.5654902 , 6.65144106, + 6.73724886, 6.8108717 , 6.89684204, 6.98270751, 7.06869507, + 7.15448714, 7.24045448, 7.33862776, 7.42446926, 7.51047555, + 7.60865781, 7.70682359, 7.79264862, 7.87861071, 7.97679073, + 8.07497898, 8.17312454, 8.27126411, 8.38183955, 8.49219854, + 8.59034035, 8.68854658, 8.79907933, 8.90941584, 9.01995908, + 9.13035252, 9.24083808, 9.35119033, 9.46172083, 9.58445502, + 9.70714802, 9.829853 , 9.94022996, 10.05076646, 10.17348567, + 10.30854414, 10.44348801, 10.56619074, 10.68889573, 10.82396168, + 10.95882092, 11.09389811, 11.2288637 , 11.36392292, 11.51118508, + 11.64610273, 11.78114023, 11.9283919 , 12.0756593 , 12.22293269, + 12.37017986, 12.51740232]) + + +def nearest_element_idx(arr, a, both=True): + if both: + dist = np.abs(arr-a) + dist_arg = np.argsort(dist) + return dist_arg[0], dist_arg[1] + else: + return np.abs(arr-a).argmin() + + diff --git a/build/lib/BCemu/input_data/README.rst b/build/lib/BCemu/input_data/README.rst new file mode 100644 index 0000000..34ee8ca --- /dev/null +++ b/build/lib/BCemu/input_data/README.rst @@ -0,0 +1 @@ +The pretrained models will be downloaded later when the package is used for the first time. diff --git a/build/lib/BCemu/kpls.py b/build/lib/BCemu/kpls.py new file mode 100644 index 0000000..9890b24 --- /dev/null +++ b/build/lib/BCemu/kpls.py @@ -0,0 +1,61 @@ +import numpy as np +# from smt.surrogate_models import KPLS + +from packaging import version +from sklearn import __version__ as sklversion + +if version.parse(sklversion) < version.parse("0.22"): + from sklearn.cross_decomposition.pls_ import PLSRegression as pls +else: + from sklearn.cross_decomposition import PLSRegression as pls + +from smt.surrogate_models.krg_based import KrgBased +from smt.utils.kriging_utils import componentwise_distance_PLS + + +class KPLS(KrgBased): + name = "KPLS" + + def _initialize(self): + super(KPLS, self)._initialize() + declare = self.options.declare + declare("n_comp", 1, types=int, desc="Number of principal components") + # KPLS used only with "abs_exp" and "squar_exp" correlations + declare( + "corr", + "squar_exp", + values=("abs_exp", "squar_exp"), + desc="Correlation function type", + types=(str), + ) + self.name = "KPLS" + + def _compute_pls(self, X, y): + self._pls = pls(self.options["n_comp"]) + self.coeff_pls = self._pls.fit(X.copy(), y.copy()).x_rotations_ + return X, y + + def _componentwise_distance(self, dx, opt=0, theta=None, return_derivative=False): + d = componentwise_distance_PLS( + dx, + self.options["corr"], + self.options["n_comp"], + self.coeff_pls, + theta=theta, + return_derivative=return_derivative, + ) + return d + + def transform(self, X): + return self._pls.transform(X) + + def inverse_transform(self, X): + return self._pls.inverse_transform(X) + + def predict_values_from_projections(self, X_projections): + X_orig = self.inverse_transform(X_projections) + return self.predict_values(X_orig) + + + + diff --git a/dist/BCemu-1.0-py3.10.egg b/dist/BCemu-1.0-py3.10.egg new file mode 100644 index 0000000000000000000000000000000000000000..7354b619ffcbc07dfc2653f7e1ac83d1e3d9b02f GIT binary patch literal 25789 zcmZ_0Q;_ap+bsNyZQHhO+qSXBw%6FUXV%!ZZQHZPw)fg6@BZJfzW*SpN>a%|cUSJN zt2#aYbO z)yBxh-ihA9b6Q)=8HX#;FV#=Xb{0fRAtR;;ph`96cN)9xyLRe2I@4u^8=a_{*v^}# z_BiJ9+Iy2Gxdfz~y0o}6h}1%z5Fxsg5%4sKAfWr*P&-lDA!El+;ZL>75k>yiYIR>5 z^>uI$VaR}SvF)}?nlw!^nv z@t(X8=i4lklqL44=4(B>v0~Mrvqi6g=0{KP_oHG5|G8Cd`jbF?Gb@$M*FgVz;_qSC z`b_1M(9{h^CACzdR+sCBpXX<(uEyc~Y_H)UrJQ)w^!41-%|z*^W4q4s5QLNc%R^DR zz2AC`yH|%!jd;>&e$IYT&!UjwOLfuF%d(W;NB&J^l8-pWldr4tolx2&dPH9?k@=rv zKLbB}0}5km7jiN4L{_4{BXwvi%rVzgnk!yHF+16WV^~YQHQmj}`vYPq|2k_YqNY^Y zM4BY3;j9z~>dD{6f;8m~v?hI;X&B5AT;;mZws=XFMm9!djmz3?!dU@K;(@>V+6J&Q z%F|WInKO=|ex;Rfbr7FSXUtQtW|~c8(=b=3PeGcMv`}nd9UEfdNIEF*({7G)(l{s7 zSExf~VbJAACyJvoQDmY09uHodTQHXge=kH7BTuhZ+31_1&7UMB)Ev#uRsp}Lil@m= zJ1f`=E^^lwsu}Im(28eGw2>c6JV zqPO6`tLhm8m0=0vNEb#y9e}Gj z7a2Oa602MzuKqT|gT~}^#>(ndWDkb+0UxNUV>gipiw$Ft=R|uo$>u1@g!VvcnZ_K0 zPH_|m4lqN60#fnF65A-JstVehEulz2lq8bEbYNcneZXQGXcJh%3qVc~Nk~tch)yZf zl}#i;Mhd_Wpg~+^2i!u}$(GvIEJ>;LTNh*}?9@Qtg}jm=BPhNt(E}Lk>%gMyiUTR+ zkeq`mH{L<}xB_(#ry=?Ou#7Tm<cbUxEdDlxIY6VU;t67v(uQcIYq z!RadE+w)#{fFlPI9bJ-~@(WA;!Tc>kQ>;4aa9$P)W-XEHhN>JC^zIerZk8k4{*8ih z!X?tlX5S$rM7S=$Kr<@mav-Lx7_h#qD&I0eHHx^K_McMN20G4_I+RN|;V2#lEN! z4f@~&a5jCj36O9K=pae*D()!87q#L1V8U&d-E>$nIidOyfLC4}|2y1p7#d|T6GfX| z|MnD=xUlZWvmqrkEOo-OpgfQPDA0p2p>a$kS3do1S#+uKs43>h9H!n;C+<~P)lnl&nBM#zvjE0XI0BUV zP9@g}gf;7k2jB{(d%%>50wm5{m<7rS6`W&grNC}*Zv%vtlKA~B8kh}*4>^$kNYzjp z0_V7~gplW$8>1b9wGejH62XwuQThlrfTSWxk%_`XH1w@I|5e1{4Ot~Ww&)3iwMeRi z&oKXC8y65Akcqi0D>cZAdM^SQKno%zKMWTjbS)w(*90Y)q3b&s@kS|7V}d}n=k_x2yJusrGf48cH((WDq?F|fP3{FH;h+Te^()QJZw$Kf9JWJB(4{qyj? z+nPn)ZG62c66^{+MiumA=eRi~kWwk41`+BXdGJ{KA;Jb#(6Y(EK8^YD=V{alBQ=Bl zuJXsVtioiOe5jN<0d~4eTkXPnaSFpd&O|boT=B_#)5mm27Zb7+v@5tmXf!LcY*$!B zKR4|oO#Lei`;DtqZ2%l8F%|kQ^%)R!Q7GJ?_mqs}Z44pMqP$y!-uc$E#jAQK;b9xx z`-q~Y$FIu5Y>N1Hx9#POR@J^Uht>4D_+GE$0}BE1I&VMtA~20Fg$W3>Z7B1#BE8n= zpDJQ7Nb0agHc0Hzm3#VHEnW|;oqzP}G)SA~WRxUOxy|Xr?B+?FoZLh^VC7m$sa5Ap zLd~AV$P#18i@UIPolUWlkF#smpqJa(ol?a5lo|M7<4++7;8**|vek3tCpr75E)KxJ z50C1&Fk0_ z{!o9H@-66MlRU;%(TPwlOa#CPzZ1o5P>q}vOxC66=9|v}TbpMve<9m8} zZpWPo8(gi{OmnksEXp2+mPoW-Z|CcTE^?WID5XJ%A;8$x`*!!0xy=%NQuF$$do5tk zLCIgDk&G)beA)#7OHERToB`D2g*=>uW*n&%K_)!l-dEfeEmPy}`@8&_3iC6xYG2=h z-4a`F%~2P&9wCuYw!41U`i=DrA`_=Z$BwPx)n#{gSGRy-N+-ly@d#Q7H zSxi;HTj>g+`GB*Q-&P~9g3BPKO&9M&dU<5phBFzC$$sR1c62<8{%Y6o?xU6z?64xk<=}9>i#- zBwvj|0`+j5aSVd_b$ky^xvAWN{`Y8h)P5R;6|$K$vusKrDsMuUa5Q_-%0^IPP)G3;AQ-mkI_){`s3IdOH2& z1DnwN1*8tC`Q#qg=05~j))#B`zBq>iDCTt! z6HWAK%UcFvvy}@{%iPFe*Bt^&t<7gdSFW-U8zms!3-U!xt7JYqyFP7gr^82$LHQqp zAvmQ3r@YNZGQ8!%-si-UR(+i@s^~=}&joAzg9&7{yuU3?V zxA$2o$OpgU3SZcDIBWzvJ1RO!f8t&U7gN6JW3d-W$NIk5-^D~f5g5q@)Tn7WaziNM zSjYwpjOt|F2{YAAhLx|?_hNPpC**lYIDYr(GCRLPqzfPk5!X@d?24REJ}-wmZPrI@LKKg? zdkYng66D60@$Ra&aZjaLidlNaP4^t_65*Cec{Yj_=oQlC7!Cx(9mZLB)U1T*H4`Y5 z=CpwbeJC){?tahjGTdbOU6Cp+N%@uk}1zrbY`r;NsL8|d~(Zttlm3^P%YOn6H z`;ZuHG1$rYlib{xY83U;3~t(gV<~rktG(N2?f3kYU7yEAeT=-!4hoQuMo{`$zI;WR zLIqd;Kq5CUX8lxTx|N{C(7h_MJ!DJ&l5PB6?=>$rPT;$vlh514tC_0sB|`pEYFDtgm35ZD zbN*7mfXwFOSIb3&m_U=(stz;f3C`nGt5;XG+w_fj&JaF5!Dv$6yBuh)XPdp1|Jd(DydJU-!QUyoH1*&~CIBw3=S-tGrwAellMa_cGkiSk^557@Tz6aJ2q5(B?ua z`iVBh1bz8=6CF96{kcC~p4QzpsyV+~-6xAVSd<%vs<}@e`gxAz?yp9`2g~N&9E}SZ zQjwd!&{mss&$wmazJ1NjF7a#LwSURFD)4K`>@9T2%{@81aL-*`PE9*m(wHt}fCVqm z&O*RyH!_$|SsBs)~9lB2>Et z6*@4!gd&_DNrtfmCHEIX14MnV=(nbO;9dxmeG4P&eWRU4t{P^q?gSFSlwWUHJaZa! z7=>#qF8V}+o7}B0dpDcZS zX@viV@2Qmr*7P5A9pL|8_)P8H?QHCgO#i{xoFZoz%!v4Vbsu#H2)UU_jibE7s^c1p zT1(y`{C7Q_w`Fm5&kO;Y9Fd%un`ZRx%Z?>$Y785E{2dnXdblM1?54C}sEr;IMVHmK z^eZLxFm{T0Wy(Goid-cp=1GzZ{&ZCmcH-;RqsBS3D;VY>XN)jh7#xUm+xUuj%0NLH zKS*Bhnl^QTCrRR~D-!(>R%A&FU0oswg}~Fj2%qrJa%J_s9&yWLhJN}8zp+kbzmFFf~+`TZBd_Q7(rJ6hewi zBCg~R6j>WaZuV%lk8*ZgtOUX9_NY|gm)~>zBR|#|3u)eyqWNQxg0|m-H730D!{;0FeKqXmeLP6BjFcJLmt)iKdJ_ z&Nvc3ga0tV8}Y29UGK}TzYA2-E`p&Q3Rp;s#9RGAJ0av862BelN}#@ zd-gTAUH<=R|I6kfq;DtGP-|_8g zPfr?Z2y>=Sy!O0qV|M=ECczc8A3b7QK4UTMUv*VA{+bOrI1lO(`Au9aCy&qV!_M&S z)@!ZU@Ysp56$d*N9%Ygh6_Ngx$oUJoXgw%$s5lleuj(X1hx;}B3h%nq{cj&X)@50p zVdk$b_Gl!H$M@wIVHb_#SAO;UxGTHohY#?n=$jvQjk*g_Ej~M+Catz7tJc2mxOFpq zI)LPt?LS6OH&2W8SC#z+47abdgY`LTl2-`ZZt2E9g+C2Dt7gl#rcVZ=ExvYt8pn?o zQJ9x-Kiq~P^c4K|ZE z{Qe?0*vkmgrsR7(0*=o{NTYZhtD?e!B|gt|wuP6?)Wh z<*FY6s>V*Nlwt=9eR^n(<mL`ptLsm{$4Kx^jv6>ht-VajSY@&5V&5+l}zkrOmeTn&$f&Q|*um z?y$7^u?J7Pn={(i9JN4*@I8lSp(v1o;I+661=y@I@>Fg7aa$~u>G-H!thX&YIc!E~oOz8kBb?ubD^CD-xJzd$0NBQn}0`K%e(!^!3d5<-s%V_OA6N_n^b_ zyZ+^s+ndAJ@2AwESHXbdvN|@%lG!ske@51A@P+AQ{B>^3b}9r{N;~uSrB4fegdlmC zqK|riJ91-`bH?CjaC3Z&W2`FC=3e!3^$O`lV?8ye%!+rda5xo)^@y7d`Ju5gn1wKH z!cv(tIU2@!5+be53Pt&l;7;t7$Km*i_~AYKsA_}VQ~&y9_WPj<^^TSr%!s(t9dV5z z3Nvygusju&wV!Zga?C1W5ff8dl3{R(;l!H!V@ZfEp199EMa4!*`auk~oS{FdMve^r zO#khgyfnL^nurvc)<8dTpGQmR(qc47SxBMq01#DYeU5H?zvg4Uu)_@Ja;PJL!Z>w3 z0F|vK=uyY*!zN5{Z=IBYlOA$sWKRaoSh%nx5zzq@wFpua$>I-liwHh{s4|+d1DE%< zbUUG}foR=PFa;19s78QHN1aJgu#-zgd}AY3c?(IDP2z?lefz1iI?!)W!GpEf$+CE= ziWu~9F0iLjkYM^GuX802k@5PXKtb(?W2N^Sst>8G!RBF&O&QPDEXyNXjAZ za!uX^G=D|WI2xW0gY8i7j%HJ|3=*H;EE4yWSE44cSbiDpwH9HxM=gfR=I#ouC*F=fxz@yCfDu}+y z=cV`#j87C&I8zs-j!|6@I*kukuAXk?U=sS3Ong~xCYjX_Mbsv z0vzZz>vOUx3yaScm^j&ZrhzxoY@9Mr66clG6OYqL1l!Dp1T;m;_P&du!Jx*OLvX2~ z_*sSD5sg_tB|+n&FNUUiB#br(aeNdLZK9)?h|*W~p$4)-@<{G;#oF7k`_|lIC|} zZ3@xlFvQr=rSn9DOn(7~3^I;Wdo*D*R0<2JC}IvDx1$N}eoFTPpW=?i5a+MLwpY2) z#u7?b`~@US*8c;$4cF2fpD$N6M!`$?P9`17AV<(DAPv6Lu`fIvoRp&_SXc*D^NNtGuP)wiuLFWH28D>F zO~#cT1H=cMrzMyh03Xqa)gKn%kXVff-4;YjhEyD?4&olsXL>r7enPSu|6;Cm3yBEz zNS$B8RKH}o=V@DoKbl6VaQGpv%XnifFRldOA=FxKj~0JeDN|x}VwS1_f=VLD1!K9V za3tb=_EPHeuC>qtCFL`glt7w8+rzB{&YW|89&sZ!Au{8G{Q6y!(i6}ZOfk7k6!4BKOTA?u_xPK{$VP=W~GCAopD+t;jugZLg-8V=Aio_Vt2oyxH2~}3WA&)i} zA?d7G1mzRW%7gS82~qq21?Cf;0ioEtM%AV=7@pEYnVompFiy;w4;{vdN2K?icj|1L51O7*S0g zq!iE$uC+Aff_y<=6J({%-*&{Jwh6m|5_EB-*3<;W=)J}Ku^_&4jk%NMbwJ6gZ*7eD zp8WK_S@O`x`nAQxWO0w%W?lvLvXIMtT(%Hv)b6vrA@$}r{Bi^(mBX@ ztmda}aGRO564Rd33;7A{<+Wd%iSc)Rt90{T0D|HN|?pNFE9Wg{?Dkx{l^vl za98-)@r;mT%QSs0{;$F}h~dsNJQH9w*YM6BH9#sQFB_lFlkQoBiWmslkjw zj)yg@Z^4_CfWWG#uPdk(!p=5!2p58R+LVeGP4f+-!{8^Zd)|LgmT!c^-> zYn^u@XQ31HIw70cOL*HPF41KmBwdL}m`famM@28Ps+u%x7&`U#3@p`Pd3DwEjS1u% z52FHyc_soKXAmvto~hKGn9A^2z(L|vxFtWZKU!z-<@nt`^HT2cTw)`_`>alnws>2P zvzDnyBcUuPgE$wN3wo~XE6Isc`?zZf4jMhdOlwi*cCZiwzF;TcBV08qXLYL2=HiUSUSaO=?RWVRN z(*>u*7mTjMJZU#}kPkT>o|bBxlIEkCB`r)kZP7ZWpR(u!gqP|E(B9<^oM#KK+#260 zP*i-jg-LvAw&YwKOfbm!0607>1?t|~Vr~wig6pLze$MuZ z%VJ10s92Mq?AY1aW2OCfe@+XIg~tnbY89st22X0eXwIF0Ud|+6rxZ#yoVA5JuhbQ{ zhzPc}jU##Wb#0R^mjpFlj}hfg#(s8>X+uc0byxa|MhHWAovEM6*0@niGph$0(_Rh&pN5>F9VZfGMnrjnu7F7wbw9dCfh^k>aL1 zkNfi5UH!bD^51#!ygBh%6je$P?B9)L0go&SQ(5w>7$kr(An;KU07jw`AX^|tDfxos zrBjlLerjY*tyiKqf6_5w@1-00cuU{5|KqCP&k+A-@P<&~&n}b#y*4V%J;r?Uk&RP? z`uxh&D-TvWgszPt&mJ-eOn?mBB|w*)N$!}LjGJR_x~aEOle8mM<(uAEpYCZf%>*eFptUYFRth#7 z4KFHH`Q4Kd$U0;qc3@6QQG{0v4u^hZ3Cq{OovAilxDXpBXJT|Iup%a01r%9)IsDKl zh4FZS+PEa^;*a}l;^0(v@1n^x)w8Da1Rp%pohN^if&lNZueeF0Mhz;0{Mjh)K}v1< z^EAvhu!6^NFov9u{qcqW!$e`c(-qgDSglX)!{1hwldgl5k84hc!x(4%hlx%7W;U1J zZ`ax>1}v|!mm?hgH~8M zo|%0C?_7^2KaXwz8Fq5*T+J{0TA%F&qID7tpMw?@;7+PVeo+rohKdk^Loj-CK-f$rpz& z)K}Nz-LukoaG&$B84vZ$+XPJO3Aa=YX%2H>|Bw{wi%^DqIi zV$thG1;5u)V`FBC5ug!g@EV3h(i=wR&-fE-uma3E*Y>xCFp3sp_>Be)yd%g z=A;*XuLta|F;4C4G^NNF;*!%3C!m5J*!MvJ5H#<21K;*ZQOkQXE8ye(3I?dumxSbiba^V(FtNY||GTo5d!+!h=FryWB6N&YtLxmX?Ki3n zfOtX3!+PTgplIdO+i{*ys^y($a=96f19(_$3?l66dQH(-sH}3#X#kJC9)!()BKDEF zSM7MMy+mnzkM6nuTf~t1wiRj(R#iB!Q(o}bQL}I$94`J`)Y^K&^!=5v*je-+f~xJh zt#{yf!v~OkSnzdTfCk79)@b)yLLMhE_+J#P?w|m*rK!nrc3}58xB9pkF0lh^L`or9fhv&P?LV1TrFJLfGR?}7uWGeBDY zdH-HHD<7@mKB^u7c285l1& zu13?{kJAi3)*XBRP3^WmMPz}cDSQ*UdkDUO2598kWcIjV1khxDKLWEfBq{JeXx^P& ztN>#^ob=D`N~3)o=^OND{_Gk$S)o3!J{f z{Ir}|flQfk^~ubPA$V#qkq%)8GWYm5vdH}kbncH?-$!GljPLzG$E}InA>h;bX+w^n z62Kdl8oU048L&=bf7$d`2VeoK&2<*sW5)lMJ^C~Mstef{;U25E(gxI*D)#ZY-whyt z;dR?-zjn&|TrWxjR|f^~b{Mxj#kc}y{z!a&l^_6B>E8W$NPRhl&1!#bC!a^yQ#(7y zx3>eW^Bn9qs)PrmGag=UM=W-}9!ZpzGtt2tPX3j!&uv2M!+k#N#NX=(Dr4|99FK?s ze6oGH?N+w}B$w!X%_p+}_@7kIVsY1K{BNas7cCZ~fMD@xxSLJb{N1O)Ro`j%hjD~1 zZj;5Cc0kYO$Nb897y#Gf`dRq<@MszNQko}_qNUxP6=jz-X(d_wMxZmU=uJidSy z#C%nnE}{dJAKI>Ue(MuTzO3L9vI8qyu7Ky_rBVf4)a~F@3dFO;U^q>BAf#-<|i&3Dqio5WSfnh)_^i0B! z4{Tp$3>yNQ4^;kriZgp0|Ah-@f2yBr&R>wWkGhKr-z=$pbGy_No{m>IDRo`Dh;F`GpIl*nOwJHZL2-aFSVJ`#J`okLCMz zmFEvn$ocu^pyeu%_j!Bp_%37xIGBp-oZW&peC2bQSHG6c^D$^!P2ea6V0OD-Y}nfY z%^kM$sWtKg>TDM}d2-KI3AlB8Y*h`w`s%}a@cjlG0e09r*Bcqv5qe4mzoa@>6%8O zB4`fhMHL?*Gs&oj(QuM~k&WkbGvt#q_LRf3h7={Lo6`rCg`;M^S-iS z^~&j?tEiP6u3%#sFKlgHIIp4Rl(53assJqpT!9QbO18*&=#b_~*xZD|f9Y7l(j^P6 zB};J-#E>VD@L@FBuypJ*lW5|sJCoUX{F=>2QKpIlN5mOJ&s-z#G)x zh~xn}&^SJ1GTxv?|7&`q!Dl5G%`{#kWn7u?E)}wfH>xwE#E`=19ka!s;i|k9L=|`_aPg_xW%2 z9NB80YLD0IvLxx%7WEtNlpR}Lo-XK|8=Z2&1vn`Nn{lFgVxi{yHCFQ6rMCP+Pz` zRBB5$^C7w_zCw~F##X0JxvU0aEF*9JP``1j9{W;>DHip|yd+_*=T`)R$=;nBfU_nYV zu3W-dbefq28=@D)1XgMF0#D2XImmwj4Vh~ti0dKlo8l1Ou|BYAWx8%Ix9sApbKzOo zLahwCfuV7wD600N0caNo-{|U=vEe?{^Au0|>HG9QJd1%0&)Oav+E7ExDuJk|hl|^= zA){az2BnJljIA(e>(c4XlAZV1hTJp@*m?(~+Jw^`LzHc_G?ljC?PWvMbhW^D1X~1lr8#S?S3*OOVJh(Di!k0ZG;w=i*mPiL7%62Y)_=RU(5XX;4UPvOM^7TZ=O?37 z<5eY;cydaos^1oGcv5^0IzJJ8pb{b<{}G#*4bKsiZ=(P6g_Sr<`-0_-@-^bbX85m#W^X20V ztRtSA&e%;$cqd$=q)Md9OKOD}7?$?r@hPpOe86Bgyr$YF^ESf$bf~~o;SdvR3AEx3 z;NrlB8acy_U{YI$rP^B#K_@IrC#POj18(s&Gr9ykPb8c&NZxrW`UF)FUo*QG`lDMf4|x2RR&KpOj6fsCbD&8QD=Kv$A`3 zDRNPYsZ1KkA-*1TEkrDnB(XO)UUV8+O6`{qq2QE^VB`5#9Vt7R7*4X)8kw^ZpO47| zST8aC94jt@$Rg#n7a3x%IW)lqQ{}R>j*v^X;I>OCseh*0n;X%HONe^ofNqL_Fg6-E zJP7Qs*g%G|K_6>iwm%q@(Irr^>@Y$I_~F)8#%l1t8uuv0ys$t@gHbYQId-YFt%bG) zGdPcQ_mpjkqo7OfJ&vZ=dIjr#8$xjD*2tSAze{}^7s>$*!Un%hCLfVRB;#=CB3JUt z`?gl~n<_^}VXp$v@&m^mD3Qgy;cc+90=TDBC5QSpIo7^;&iU^U5|5T_rRH&BIwY@IHMFxE;5*J<53#cb5Fe2r7Ev+0p2 z8z71lb!gv}8N`_{Otn<#K^uT6LRe5v^<7GJjOg%6ctx|+2Q!B&GKCOgl(c5TN1dYS z_4D6U1JLOR4Ot_)%HE=piXLZ}y)H|tB7$wh8Gzwq9gaxatoQb0aUTK(r8)2KHZ}L^ zGCApi;BofHBIqws)pmE($Bn%if{qd<^w8!6CMo~lMrhl`2Yg3d38tRyhIMP6#Mm@K(6pK}()vMi(F zxT+V;SXa<6o!f}#Osi^GFWrgb*;ml8oL`FPY^!S6EkTLr7*<@^EOAA1OZ3E6wZMoM zidJ+kQI#o6BM>I}u~eySNxlV**QaQ(~%5LCN;g0weZF zl1gb5yNXP16Dl+vUTJFAd?JT7xoBs3iJ&$&3|*8g9avF`VN;isOS(3 zpZNu=CL5%G7xjwG91F}M%1b6K&zQA)sNI;GY01s1ht@Ag^~jScR0)@)%A)2DRDY$) zH&>JxJMLh4N}s!+l{^fb7o^OCWmHo_+O-BfK2wL+brO#L^Viu_on>-yB_pCuNWOkp zP+_im>1*u4IhOuK8y!Tah{D z8WGEp$r4ca0z_xj;+zRTmR9YKOvE$9b`hp{tdQBJEO4lf5flf^q{gf+ zVPM|yf&(o%s39`I+h+~MCZ%O&Qq5AVw6cuwJp}joaA;Y|XxC*D_cxyKQ<>6OJn}@n zMv))AUE~Q(Zr7X(M&zw=i@_acpiH2DGezx9C%9tU)2H`8tP?YXU=)CAU|#g+(zvLR zeYK=K<#(Jr(-Y_Jic$9eSxK#GpuCts&w1BuBS4J4=+G^=XrQg9HXE8bLKc1lL7rOX%W-M{H3pNHSFFDlL%|SvWL~l{OG~dwGc^;EY~y(Fv;98VwT3<6_ASTH6o+Bd*Ex1 zoJBTRd+6VIQ)>Q82io|I$g&GqJMnE&UEcuHp!goBKtWof(xUL3z`h_C#)?(SG(~YQ z!9q)4XW`$|8LRqd8^U%B@sRR8>6qYw=#g*4xyL*Y^cWJ>3LVg61_8yPDz+LxZYwkE zT=SlM%YujgthB8k{RR}@U+&x>8!E``-AsY&klm&pX(e+_-!3>d;4q8w0&H0F6 zPn9f4NfB-SfqFs~wF;D^#a#_C^aZB1Bt_RvAI0j!-$|}%zLAw8lp@v=OrUY&G3}i< zqy`&ZY&qmyke_G^q;EBxBMHd(wHezowVr)IXbubVvb}ElMl3b8H;$4EQU7dqT@Wq>G<^*jw5^Hh;T z@DysNh~=&pwu^~o{=k{**C%whY!jv#d(%^K2%Ge)FqB2|8Vs3wktk*IS~N0*2TAF; z@bwDVE2n8qhQ5%|`qR>r3y6qnlEyl1_y{O8gN76C)D*f=3TnGy;61qx%aJOQxKvr! zDIr%9$>^|_hcj}2LvXI6M=!+ANI7)=Z+T5Aibqi`nWUnOkfPUU#R<2CUkUTSm>#_t z9#s~OCLT^(VehCZ>8><01G|3(@!#dL?vi*7*d(+Ur*?!yc3MV~zxX^vaj)n@>-Ur= zG0}L4ucO>L|bd4LNuI zG+A7NFnpK^^$`|7YDBew0PZ70vJH%re&LUGNGM#~tMMHqE*7IR2_kql?M`#}Rv`lQZsL){@^F>6Ud zYr8G=48I%6_TahkQzEqtsGmuGxmcMhVl9a51-1@~Izuy92lH##hxEOv5=+z@7*+S; zmA1cm9Irv9r^PrSDtKHhcR$7<6}wrd=|3+Kx$lr_xyWDGK_Q%fD_1`uykeI_8e21Ur!1_a~i>dj*ZhRCD*vWZ!D3&_aO!#-o@m ziZTM*)>9b`JtTs0X9xqr{O5V&)69xHaY_ zh!(n!-j7zH9eC6`76P^@gc-!*{5aq~a29<#(WGgGIz9n{stMWN9k8 ztN#wY6Y4CgvJ=qEVr$QWHT^WN)J8bWj@OWy*HhL(EXq%&E(^glrl;5) zyTpjt!ZA*cvn{nc*BMM`n~D++T!6;yqCpe_3TH)Ze`z}XMsXi^J;2iEEjWHrT*8eX z80iR)J#v4NQE&0-Ya*{81NwZtxx5b{eoFB^i@a8A*vSEaWeYcINo?zWx}kqIG5?#lX{w`5+p z&VX|9~SFM=vA%2Ygq#cHFNauIU8D#w^3Sy85nU;~nFCy~|M z5xNrf`lWgnmkoKrd(|{@9LHD68&#y3!3(8OZBJ~y*~9|(9y@(P7Y9U z*&6CAsCr=3BkT{*f74!wxyPdYk8**3N&f#%`+w@U|AY2sRRQ}&LDXL7>>(*^$qm@( z3~FK|9(gJd;1EdT%YlGfW3KkfU&GH@JuG0AxrRp~O*lSa{18bRpBxVY+Nwpah7*hyT0^^ppdNn1m{30SW;+fS42<5tH`uR6&0mV94vE(v4;=|j1gqCbUU?_lYPg3 zHwL3AN(D^x9rD3tIQ+VX={c*k?IOx{nG5|c!WS#H!%&ETzjYUFB(g=FUHf@sbau9o zOg5qL=mdJumD9Adonla5ihUQB$5#gZyC(axEIKDMpFS$IdV4cDu-96K81Z{><(H(k z+Ew8h(7HQ{G~7$;7DIZo8S`S#t@t4E42ri6t+I%YX`wNa|9r8r#^I&Cr%I^2EtL|h zt%${PkLVjlps3V!Nk62whlGmv3d=j)z;#wta%nX$)^470tjobnxPBg|psH24RL8%& z@_dnZszv5es_q?sxHHD#rdc{4I23rT1l`S=;MKU<{N8E&gWY^p7R@3*0 z7Nz(7-DsL2g78GR{Xa5%@r_X3KkxuR^S_@D(EP72rvKMQ{vQ&JNr%xy(Lf)yAWwn- z5k*zajYK1;KQ5uvoX1j{FC4aOu90hIV%*0jM4*5|lSPHuaP{q^>*x`zrqirMH%M0V zsqbjBZnmy0I*6($d#zxIKhNa>1?~NM-aN`m+MW75*~xI8={}LC-PwT!$^#pYpxXHc zbvZRD^_Pu+&k_f&f7!tOipEw7jqK5op?Ip zR<8U%m0Sf-980@hWJz#$PlCHUBsjq>!F6$mMFRwavv_cKx8SZpf@^Sh*M#7Y{P*1( zZeHrWdQ&yEQ?+wW_cuLts=uDE&vLU1YARGr6jcv8!I8$iwY@=Mu_~AOv`RU0r*_E; z_zo!BQ&ZZj*^a2d zVX3G#9Cmh7)7o58Q&TCn^iXjWf_hChXm)!kw!%vTqm5eK5#!U@xuHKidotloV=w`>fQ@|(k1pLSWMW0_nfF^|z>e;rw4qKW)0#jI*Ztj`jA!|m0G38X54 zE3wMFNN+i`z)3q+a9qFKET+Z`#fQ8nT@hk`t>UANYGz>%3J%dL=Yy>rGTst1$_0u^ zG^f!jQ9h{>oQz=Xi(CYaLF^YJ%%&i$b#mFRmIYN<3|6d(R(jvYtuhBz%j3|5aqFu6 z||I8%t&tGG7j;|2jY+!9Q{EYJ`F40 zq)vnjtIBp3?POJ^ZP{te8U<$0>|muxeH-`NUO5oimo^?Kh|&_#aK+{c5wM#q)CjNV zh5QDl2L93DcoUI0?I4{C1ySUgIo_Jir3+R++4+HJOsDaPzC_NL%>V}Y>^Kx^?L~@j zn5~kWI@*i4uqqK!?wK)t@m`;EFoxJUukcmOnl1?yPE0apVmfXG!{!M4D<+<26Q+pO z`1=`7(%T`<&#!Xo*%2b!7xo6+KDJZBis)v})8Y*}DMlKj%eXb>4^?fu?vQ+Qg6Lr; zu&;qr*7=hBaUKMSP1XgmbX>_$S&vPM#hg{!&Ve>%$dq%|u>o+y&W(?~tU#(H<`&6w zkd!8rsUDO3y86_&qgg67p2qb96;XD&<>}LBT3D9Tn4fJ8>~V1qqOrZ5lffh|{x9|` zp3O4$n4~bvo`6VzBJxAZa<06hWF-6=qP!bwJhE1APuaicXsfQTz!FytOVWRwRz5A@ zSI_p17A*;1f1dIV&CA4LHr({P0<}}Tcld~llqWj2ThXN6rkE>1GN4r>v;aTBMbBF% zo55C2MWa|L6=C7ioO`)hIMX+m1E(CUT%sg#^7D=@yK&*@^$Vh;*gj@cqjZfQnvHYvzIY@lumT8zpi5iV);wGi{B` zq?9EU;~C(IVWEqWjAex-ItQe-*P{rHz7P-)yN$gYTlH})erKJ@eKCR4mQlgzaClR& z4_lNmQq?|CMo99*RuX2FsyY1AGT?^!gAgZoo-dgtuOTM&pq4DkRaz$3?F6O_wLz(yrPU9+=$5t8IE~;>`56^w(ZIRV?CHPxmbVSHa@eUq_#JKu#$sknenB*>Ej zQg_(Cllz&ZfQclGi&6XP-TRir$WplGe$2y>tPPZAz3_HvfZ6ce&W7(`n(67&!@<@9 ztAetro)76uxm4rKUkbc#AAj^Kperm+ zD4e_JvKlMb}tJJtzrvDUjuD)I?dhzwEMC^61nOVrGN=b>a8UhGJZu(=oI zrlhHMz7rpC4Fl#=SGi~pbU+Nh^AY1WRMy?G^AgKy!b$mrTntj2{*rikMn$TxIH_-> zv5YIBo}D;(Bv`x04QWv+(Pr+oj}l9YOHO2KJi&tSd(I~bI)5?4(Z1rk)Sj6$QCEqq z=gV~(W}j?J=!LUuMRPmL_o+@4)gS5z_CJ2S zgE{(SD%m-v=}>S@Zx9|-dX)v}3JW(bh$1Q9-U^wXYpenI{-9eFTD?ZQGGCTENJA+% zV4QS5r|?ov-s(X07c>}6#r;B>#0AZ-*sMmuf4hYZfF78|WR@YLOep3o*`7F^@<(pEL{B$r}RWACV3nhdW2Bom+A zUO;%+D|*eYpM~Bakp-|#T{JvB-eN-R+p8B+J1^_4jCA0HE~ci_!t7`QjZXR-jAOC2 zjV1Xxy!M2fu5r96T>Mm~#_rL`0w{qteecn+I9>*r{g$R4zYw=PlAXpy%F_emyD5hB zUJ;=1>4(x_IxV0~HoVr#1?t#GwRuC1`%A`gaM<%bFd&!scIyQ_OyJXD0xebeV_VQf zK(C~Jwdnhff7q0=&Qp(S}+gGEPS&o#i@n9XuZ0?<>1%4OCGA_^`I6?CzMe%WJ9 zJnOIrW$s`etvd+?T2Gm^E1(<9pjl`wMYu&vcmIG)tN~MA6nFJO*3#2GB3T4^ zw&Swe3;>jHpwjsPwMjP^N7Pv4-El#%UHqY9;+p|V1QfX0S@rcW}=Gg(mY&58u-=TjD5a`Ee zo_cA9-mDQ0p8#NtyTi40;1u|51~&*<{Y9!+X?28f$eh(CqbKR0-cZOQ2G>E2aLe8} zDo0TTIA!a4j{U%X8D8sLX0HvzYq2=KxqjZn{W&)-;gk#NW`*VUs+kHP6rV|1O0+we z;U`|vI?M&wnd)<9G?*Z2SS~!mG=<8aS3TE1H1m>&v@FVZRj(q=i;^B{R?=C0^p)AZM{s?PD z(c5bZ3SP4SEIv29cRhaA*5Rrf-iNyTbSPO-@w3IlcRb@KDbr8Otvs~vDL4R~JzXlA z4bH10u`Es>P2s2)Y8({}c-lZ$r>8a>Dddw5ttsDjVZ2C*O*;B;Vbr1=E}~Nb)8_+V zD(aCm5yVs~(tZ|zZx0HG&K%sx83f9LBwugz$H9tSS6(jQ#`Bw!>x8jaVsNV?jXNCk z_5K3nz#rHVX;DdeT?+luU}!xc9FY2v(9M(mMWBek##H+kAmH}UrY;JMs`m=S%)-Kg zXT6p;Tq;@y#3MZ~1TQOlD%ce~h`Bl7?$QUYK;4+0m zH4kdczSuKq>98+jxj588`V0U!Na|)h!r+zjjhHLYLPt~`9#aJP3#-NUfADu6Q+Wp z`rhMTRj$x_IXk*G#zE;5MctEN*t>9x_0HEJzFLC}9Jk=n#5850fkiTFw9ER+w8PENXF+o!P6#7i_&h(Yy8@o=)2zSX;dy_B0OjF$gxH zU*(3#tjn$2&N{GiT|7%;J0OCa6zOh2Hhm{i$OTQp?vSC!^?rZZ6q zj`()wFv>I3VnVA@Vj@!wJ}SUO#(MwR{~!)D=+c{{C+{QG5B()9y2vlu;o|G|5b^x} z?fF_(+Bd8!bvr*)eK#gf4RQ}2nzKTBu;P6z0fVafZk=S6Xu>)@*n}gA)OG~N>BX@_ zY9;(ztsoQqby|t~VYT%FZ;nKX7V{#ML)gM1H|~JKvZ0lEx3HFw*{4QXcBD~)ExFRm zxW(Yacia=jn~2|4h=QMY&wM!G<~zh}DR;97uWGTiOQVy&q*c&ADEf+sN<5>(j8ErH#mjNl)**uk zw^o8+>cCp$$dRmT#svwulrD}2R|i%)4CxVUk-U?ORm)sta0DiGm)A0KXE$5rYaea> zo3Upl;Ha*2LMn#cZWa!9YxFGwWN_#ZiML6$HsRI_A>{YFhq0-0wTBB+mzEZ_eFOQtUQjxxg09!lC>|uI zP`JX!AE)jci3grw<&7pr236v)WeVIJjx#u04lqAhV@dDRbju-V1&Nrnqb@wJL1QNg zB1VCKIac4H+^&42;9 zcan1502m`bKKNI0HgcX;YG7^!fS$#25gzwHb!_?j!MqUWfJz` zo*~CsMt90G6eK~Yd;O-$JP?tMhFl2et4p1y`DQ^0^XC=L1<@&|%O;xxZtUB)gU2+; zY3f3Z->kBM)%*I4nX7f?zM;H(8)@VQ)TT zH3;s$O;=5vuvgzLvJ7zxHv(>i)f}{m-WjS3^9m3)@7W;q!$)5>%uj*fZrx{_2Dj4G z0Qb*7z|ntnnKl!qDKL_psmd9U=I9SRd>ynvubJBR1cX!#m450)NabCzy-cY-qgfp< zYJ&f|soVBsVtk@^8-uNUa`m+I>9rNpzVA?*I%jieH=VJo<49npxrS=;e)o7)JaoMf z!hm(IA+=XBX|-MTRSt&^wqaf7Aa*CUH%1aMr$nF5M~qGMH>I%lm|El~bocFq+H9@@ ztZZKP1cswoZ`o@)wB!4%uK1v1Ji(Fi#*XIlO03GWlNCxnT*?QI+?_VDa(A(5D!RQ5 za&c+69^j`fH%(`co;i(cS?} z+~1Du#f${oE~KgjR1tykC;Hs+qL=0MHEavtECpsQn^qO`N@38~kx3-FNOPt)6v!;_ zyl%#?+`}#Nh!k+!>WNTzlAB zuymqNwLWW=Tu&4}5A`Ud4UVG2M9BI=kA3=KAUh9oxA=N zo&dWw_d^U3QgF|Ul!3%|!I*gc7BJmw$LrabWR}%C2Hou+qxZG8T_rYml1WGf4`o37 ztxg6X3RDYIK7?#!GQ;MxxJF5sVLKF|@N_VZj=+O1TnHl!dUDH z%IX{St{iakMu~lhXZ|Xk?W^r=XHSJHoyl9|(mi?^41J#llt%NFpL`fYY__UfP{+yH zG#m8in`E5Gcxj$$Q%ZHD#>g-+#acE*iFO(v!;&-B)nY|oXSeNULz9!M7TwDa&#qS2 z76VyrpuI6@QK&f#^aGZYcA!1gDaVkzG`5@A8StX)+_+P9Q7`<8<*tjhQb~FTvo$26 zkINtKN&-LeDe((!s96UtbYw&E7^{y16l#7_9W5^0q^8!vM9|TA%a_-Jc{CVHnc*Gn ziW)Kr+XF_`u@lpY2^vMk$M8=&mC8rgLO9S~7AVwgzrGtrunxhkPcLo@vQvoZ={IfE zhMSEkW8#oi3;>Pvdfleb=tGX7h@Ve+JLI~{&UOZ?b|(enJ?~Z%GGn|U>X*E<9*r{^ z+i;yyHLvx<%CuvQL=oQ(?xMdumO`SVZ;bFC#MY?Y7%!cu)O@1V1xk2f8(KxvN{yyj zfMW5cx-AmE7lerCJYkell@c;fQ;cRLwMZKKRe5VZb0Q#XlunCB;!Qw7nD5>mGZ6!` zC-ClgKTkd_nE>^eMq1W2CPOejvnjocV60d&N3IL@j$PhgC)@2+*~98&=&^h!FQ3x| zidF;KVXFq&oRZ_W1IrgCJHgJ3F)--ZPE|6h#1}lfq z?_Fo5wyAA=d>dPJl>u0A;iOxB;lYZ%_>23#R_ua&2PNqn*bvzv3&BiB-bkgLP1Pjx zqpB}K=B^}72Xv|MKXxm2-pr{?+>e%j!eyV0)vXpGTDm?R3AH{Gv>xtCfpUF5f$ z;89qU;m9*Jjjmhmv>h+eU93&5fR|4)T3eYBF&}DpA~x}%LGyaE--9%d0wUPTcYE_` zV7teMR13xPvO?qPbYsO4y?G>P7%#S<)#Z_@$jEFq7I27gsxK%&tZs4DhZZtumngmKo2VgRqunSVv=Yj7k+g?y zn-$M7NSb0ZQYN}ldkAxh@qQ<(H`4beQnu5imAy>k+nVhsFadj~gqTKw$=j0F_XBj} zSNrzc&-I8}LAj()-e~Hvbnl1hgrkT{)XuSEe1F)N8QPxXNfVJy6JVA#MNNFMnHSY1 z3-8Q9Se;X-Ip=J?iQdXB9cdA2KW18e=v#2-7(Y9I$mG3ic2B1Ol*xE^vMq(T_&IBD z8lB(6;8^(8F851s>lAvuwWSqty;?u#X8{CHsGQR!&^=?Qnjl-wYq^d$gGI`{B|*0- znxb@b7&ei`TGF#9E&UhGx+~Y$4|`miy_*Xp3!}JBOrR?8;SOZq7pWQ=a6E_ul6TT@ zFI1kKYo#3r5<$tgIp`ebnmwwuZgrTeQYNAVKk+<$rbssD05TnV2IB!=97MSBfY63? zYXn^#F-MJ3U2*BdNOgry@w_uDXg4n4{#R(Y(EDy++$pTy+)dh&QxWrbYroN@2BwJg zd88c8tHBTWo;J*T+4rrmXBVR-?m4=4%x=LDWk^wVs!`83`rMSIckf(oJl|h*p4_d_ z$uukTWy+4~FAH~zzM<$tVo@Dp8ss9h_V2{!!&NH^pdLpc!=Rz$w-)r(+hvnQJhX98 z%>)v(R`^geTB@F?KDRq^aRVIh&JGrChznf{`w`#7?3H5Buzni^q?~_Jq)N&oh4|1tPw{8M(3m#SA`X(61I(TQa%06!mm9h;x z>1t~3nTs^uE}A`OOT91_{-jj3Gm;nQ^4>=r*cq#z1=AGH>aqy$3IfOSqs!59ynFQn zYwOCQ%bp?kY|^BDL!$nw5{Ce1oMrKcfq(mw5)(IHR`iw=h6>xzYtc;ldW$IcEA=-< zVMTDaD>}D+Hr5Feuj4Qhyo<$)zFBj=eWql^YRuu|StR?4>DX+bPxtdB**v^m_<`j` zYAJM2pR}Ddl{Tm1!X(Og-=y3@dqMMT!6qSf)7H4IJGbFrS2w<--4ENDdF}hc?}K-N zXeQdM-i5DDzl31KR`m0TZL=XmtcE=OB;3f+#$7-T?{*?}@b;IR_~ve}`737i`A;9F zEq|N0|4cEPPd16lALdrP5A_~DE#3c_-6+eeDv5|GGlE^fRuSzeADGcZk8Ti;k%K4E zu#rAHjE^vwQP0gZaBJoV+2c30EU%xc;(`pKyj;z6p)wA8q<7kl7Lcq9m zrXt-&{&$2Q;igz&9wK`P5J3DmLXeTIku}K3+R#6k?ZXLHh$33&|6oGx7H-7sCSA8%Y`wRY09m@ZN z|87$D3%-r{H~d$T!2iVmZfEfefA{e6{MGz_n_E2Y`OhKv-C^JtG?(me=>OYq;P*y; zug&_^$SB{R8u{l*|F%fi!=EJZG1Z?G@E6s>W2(n`_?Y5Pa`=k^Q}?ETgtR+z9~u7m)S&&;S4c literal 0 HcmV?d00001 diff --git a/examples/test.py b/examples/test.py index 73c0a81..5eb3b9d 100644 --- a/examples/test.py +++ b/examples/test.py @@ -1,6 +1,6 @@ import numpy as np import pickle -from BCMemu import * +from BCemu import * ### Cosmology Ob, Om = 0.0463, 0.2793 @@ -21,7 +21,6 @@ p1p5 = bfcemu.get_boost(1.5, bcmdict, k_eval) p2 = bfcemu.get_boost(2.0, bcmdict, k_eval) - # Read the BAHAMAS data BAH = pickle.load(open('BAHAMAS_data.pkl', 'rb')) @@ -122,4 +121,35 @@ plt.tight_layout() plt.show() +exit() + +import numpy as np +import pickle +from BCemu import * + +### Cosmology +Ob, Om = 0.0463, 0.2793 +bfcemu = BCM_7param(Ob=Ob, Om=Om, verbose=False, + below_kmin=1, #'extrapolate', + above_kmax='extrapolate' + ) + +bcmdict = {'log10Mc': 13.32, + 'mu' : 0.93, + 'thej' : 4.235, + 'gamma' : 2.25, + 'delta' : 6.40, + 'eta' : 0.15, + 'deta' : 0.14, + } +# k0,k0p5,k1,k1p5,k2 = [k_eval for _ in range(5)] +# X = np.array([k0,k0p5,k1,k1p5,k2]) +# Y = np.array([[0,0.5,1,1.5,2] for _ in range(x.shape[1])]).T +# P = np.array([p0,p0p5,p1,p1p5,p2]) +# pnew = np.array([splev(z,splrep([0,0.5,1,1.5,2], P[:,ii], s=0, k=1)) for ii,ki in enumerate(k0)]) + + +## Beyond the trained k and redshift range. +kk = 10**np.linspace(-2,1.08,50) +pp = bfcemu.get_boost(0.0, bcmdict, kk) diff --git a/src/BaryonEffectsEmulator.py b/src/BaryonEffectsEmulator.py index e96aea0..9c9097e 100644 --- a/src/BaryonEffectsEmulator.py +++ b/src/BaryonEffectsEmulator.py @@ -4,7 +4,7 @@ import numpy as np from scipy import special -from scipy.interpolate import splev, splrep +from scipy.interpolate import splev, splrep, bisplrep, bisplev import os, pickle, pkg_resources from glob import glob @@ -145,7 +145,7 @@ def check_range(self): self.maxs = maxs knob = False if mins[0]<=self.log10Mc<=maxs[0] and mins[1]<=self.mu<=maxs[1] and mins[2]<=self.thej<=maxs[2] and mins[3]<=self.gamma<=maxs[3] and mins[4]<=self.delta<=maxs[4] and mins[5]<=self.eta<=maxs[5] and mins[6]<=self.deta<=maxs[6] and mins[7]<=self.fb<=maxs[7]: - knob = True + knob = True if not knob: print('The parameters provided are outside the allowed range.') print('log10Mc : [{},{}]'.format(self.mins[0],self.maxs[0])) @@ -206,8 +206,13 @@ def run(self, BCM_dict, z=0, nu_dict=None, return_std=False): class BCM_7param(use_emul): - def __init__(self, emul_names=None, Ob=0.0463, Om=0.2793, verbose=True): + def __init__(self, emul_names=None, Ob=0.0463, Om=0.2793, verbose=True, + below_kmin='extrapolate', above_kmax='extrapolate', + above_zmax='extrapolate'): super().__init__(emul_names=emul_names, Ob=Ob, Om=Om, verbose=verbose) + self.below_kmin = below_kmin + self.above_kmax = above_kmax + self.above_zmax = above_zmax def print_param_names(self): print('\nBaryonification parameters:') @@ -218,19 +223,49 @@ def print_param_names(self): print('nu_Mc, nu_mu, nu_thej, nu_gamma, nu_delta, nu_eta, nu_deta\n') def get_boost(self, z, BCM_params, k_eval, fb=None): - if k_eval.min()self.ks0.max(): - print('k values above {:.3f} h/Mpc are errornous.'.format(self.ks0.max())) if fb is not None: self.fb = fb - pp, kk = self.run(BCM_params, z=z) + if z in [0,0.5,1,1.5,2]: + pp, kk = self.run(BCM_params, z=z) + else: + p0p0, kk = self.run(BCM_params, z=0) + p0p5, kk = self.run(BCM_params, z=0.5) + p1p0, kk = self.run(BCM_params, z=1) + p1p5, kk = self.run(BCM_params, z=1.5) + p2p0, kk = self.run(BCM_params, z=2) + P = np.array([p0p0,p0p5,p1p5,p1p5,p2p0]) + pp = np.array([splev(z,splrep([0,0.5,1,1.5,2], P[:,ii], s=0, k=1)) + for ii,ki in enumerate(kk)]) + if z>2: + if str(self.above_zmax).lower()=='zmax': + print('Redshift z>2 is outside the training set and therefore set to suppressions at z=2.') + pp = p2p0 + else: + print('Redshift z>2 is outside the training set and therefore extrapolated.') + pp_tck = splrep(kk, pp, k=3) - return splev(k_eval, pp_tck, ext=0) + p_eval = splev(k_eval, pp_tck, ext=0) + if k_eval.min()self.ks0.max(): + if str(self.above_kmax).lower()=='extrapolate': + print('Suppressions modelled at k>{:.3f} h/Mpc are outside the training set and therefore extrapolated.'.format(self.ks0.max())) + else: + print('Suppressions modelled at k>{:.3f} h/Mpc are outside the training set.'.format(self.ks0.max())) + try: p_eval[k_eval>self.ks0.max()] = self.above_kmax(k_eval[k_eval>self.ks0.max()]) + except: p_eval[k_eval>self.ks0.max()] = self.above_kmax + return p_eval class BCM_3param(use_emul): - def __init__(self, emul_names=None, Ob=0.0463, Om=0.2793, verbose=True): - super().__init__(emul_names=emul_names, Ob=Ob, Om=Om, verbose=verbose) + def __init__(self, emul_names=None, Ob=0.0463, Om=0.2793, verbose=True, + below_kmin='extrapolate', above_kmax='extrapolate'): + super().__init__(emul_names=emul_names, Ob=Ob, Om=Om, verbose=verbose, + below_kmin=below_kmin, above_kmax=above_kmax) def print_param_names(self): print('\nBaryonification parameters:') diff --git a/tests/sample_BCparam.py b/tests/sample_BCparam.py new file mode 100644 index 0000000..c96a0ca --- /dev/null +++ b/tests/sample_BCparam.py @@ -0,0 +1,93 @@ +import numpy as np +import matplotlib.pyplot as plt +import BCemu +from tqdm import tqdm + +Ob, Om = 0.0463, 0.2793 + +bfcemu = BCemu.BCM_7param(Ob=Ob, Om=Om) +bcmdict = {'log10Mc': 13.32, + 'mu' : 0.93, + 'thej' : 4.235, + 'gamma' : 2.25, + 'delta' : 6.40, + 'eta' : 0.15, + 'deta' : 0.14, + } + +z = 0 +k_eval = 10**np.linspace(-1,1.08,50) +p_eval = bfcemu.get_boost(z, bcmdict, k_eval) + +plt.semilogx(k_eval, p_eval, c='C0', lw=3) +plt.axis([1e-1,12,0.73,1.04]) +plt.yticks(fontsize=14) +plt.xticks(fontsize=14) +plt.xlabel(r'$k$ (h/Mpc)', fontsize=14) +plt.ylabel(r'$\frac{P_{\rm DM+baryon}}{P_{\rm DM}}$', fontsize=21) +plt.tight_layout() +plt.show() + +n_samples = 1000 +mins = [11, 0.0, 2, 1, 3, 0.05, 0.05, 0.10] +maxs = [15, 2.0, 8, 4, 11, 0.40, 0.40, 0.25] + +pses1 = [] +for i in tqdm(range(n_samples)): + bfcemu = BCemu.BCM_7param(Ob=Ob, Om=Om) + bcmdict = {'log10Mc': np.random.uniform(mins[0], maxs[0], 1)[0], + 'mu' : np.random.uniform(mins[1], maxs[1], 1)[0], + 'thej' : np.random.uniform(mins[2], maxs[2], 1)[0], + 'gamma' : np.random.uniform(mins[3], maxs[3], 1)[0], + 'delta' : np.random.uniform(mins[4], maxs[4], 1)[0], + 'eta' : np.random.uniform(mins[5], maxs[5], 1)[0], + 'deta' : np.random.uniform(mins[6], maxs[6], 1)[0], + } + p_eval = bfcemu.get_boost(z, bcmdict, k_eval) + pses1.append(p_eval) +pses1 = np.array(pses1) + + +# plt.semilogx(k_eval, p_eval, c='C0', lw=3) +plt.fill_between(k_eval, pses1.min(axis=0), pses1.max(axis=0), color='C1', alpha=0.3, label='Planck cosmology') +plt.axis([1e-1,12,0.3,1.6]) +plt.xscale('log') +plt.yticks(fontsize=14) +plt.xticks(fontsize=14) +plt.xlabel(r'$k$ (h/Mpc)', fontsize=14) +plt.ylabel(r'$\frac{P_{\rm DM+baryon}}{P_{\rm DM}}$', fontsize=21) +plt.tight_layout() +plt.show() + +pses2 = [] +for i in tqdm(range(n_samples)): + bfcemu = BCemu.BCM_7param(Ob=Om*np.random.uniform(mins[7], maxs[7], 1)[0], Om=Om) + bcmdict = {'log10Mc': np.random.uniform(mins[0], maxs[0], 1)[0], + 'mu' : np.random.uniform(mins[1], maxs[1], 1)[0], + 'thej' : np.random.uniform(mins[2], maxs[2], 1)[0], + 'gamma' : np.random.uniform(mins[3], maxs[3], 1)[0], + 'delta' : np.random.uniform(mins[4], maxs[4], 1)[0], + 'eta' : np.random.uniform(mins[5], maxs[5], 1)[0], + 'deta' : np.random.uniform(mins[6], maxs[6], 1)[0], + } + p_eval = bfcemu.get_boost(z, bcmdict, k_eval) + pses2.append(p_eval) +pses2 = np.array(pses2) + + +# plt.semilogx(k_eval, p_eval, c='C0', lw=3) +plt.fill_between(k_eval, pses2.min(axis=0), pses2.max(axis=0), color='k', alpha=0.3, label='Cosmology left free') +plt.fill_between(k_eval, pses1.min(axis=0), pses1.max(axis=0), color='C1', alpha=0.3, label='Planck cosmology') +plt.axis([1e-1,12,0.3,1.6]) +plt.xscale('log') +plt.yticks(fontsize=14) +plt.xticks(fontsize=14) +plt.xlabel(r'$k$ (h/Mpc)', fontsize=14) +plt.ylabel(r'$\frac{P_{\rm DM+baryon}}{P_{\rm DM}}$', fontsize=21) +plt.legend() +plt.tight_layout() +plt.show() + + + +