From 4d58c599bcd3e5431fdaebfe88b2cb71d0836940 Mon Sep 17 00:00:00 2001 From: "Johannes U. Lange" Date: Mon, 11 Jul 2022 10:39:59 -0700 Subject: [PATCH] added database capabilities and corrfunc wrappers, switched build system to flit --- .gitignore | 2 + LICENSE | 1 + README.md | 28 ++-- pyproject.toml | 18 +++ scripts/download_snapshot.py | 259 ++++++++++++++++++++++++++++++++ scripts/tabulate_snapshot.py | 235 +++++++++++++++++++++++++++++ setup.py | 12 -- tabcorr/__init__.py | 8 +- tabcorr/aa_cosmos.txt | 41 +++++ tabcorr/aa_test_cosmos.txt | 8 + tabcorr/as_cosmos.csv | 100 +++++++++++++ tabcorr/corrfunc.py | 175 ++++++++++++++++++++++ tabcorr/database.py | 280 +++++++++++++++++++++++++++++++++++ tabcorr/interpolate.py | 38 +++-- tabcorr/tabcorr.py | 5 + 15 files changed, 1163 insertions(+), 47 deletions(-) create mode 100644 pyproject.toml create mode 100644 scripts/download_snapshot.py create mode 100644 scripts/tabulate_snapshot.py delete mode 100644 setup.py create mode 100644 tabcorr/aa_cosmos.txt create mode 100644 tabcorr/aa_test_cosmos.txt create mode 100644 tabcorr/as_cosmos.csv create mode 100644 tabcorr/corrfunc.py create mode 100644 tabcorr/database.py diff --git a/.gitignore b/.gitignore index bee8a64..298aa30 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,3 @@ __pycache__ +dist +.spyproject diff --git a/LICENSE b/LICENSE index 918991a..60c53db 100644 --- a/LICENSE +++ b/LICENSE @@ -1,5 +1,6 @@ MIT License +Copyright (c) 2015-2017 Yao-Yuan Mao Copyright (c) 2022 Johannes U. Lange Permission is hereby granted, free of charge, to any person obtaining a copy diff --git a/README.md b/README.md index 24b41dc..c190043 100644 --- a/README.md +++ b/README.md @@ -1,24 +1,17 @@ -# TabCorr - Tabulated Correlation functions for halotools +# TabCorr: Tabulated Correlation Functions -This Python module provides extremely efficient and precise calculations of galaxy correlation functions in halotools using tabulated values. It is specifically intended for Markov chain monte carlo (MCMC) exploration of the galaxy-halo connection. It implements the method described in [Zheng & Guo (2016)](https://doi.org/10.1093/mnras/stw523) of tabulating correlation functions that only need to be convolved with the mean halo occupation to obtain the full correlation function of galaxies. - -## Prerequisites - -The following python packages (and their prerequisites) are required for running this module. - -* h5py -* numpy -* astropy -* halotools +[![PyPI Version](https://img.shields.io/pypi/v/tabcorr?color=blue)](https://pypi.org/project/tabcorr/) +[![License: MIT](https://img.shields.io/github/license/johannesulf/TabCorr?color=blue)](https://raw.githubusercontent.com/johannesulf/TabCorr/main/LICENSE) +![Language: Python](https://img.shields.io/github/languages/top/johannesulf/TabCorr) -This module has been tested with Python 3.x. +This Python module provides extremely efficient and precise calculations of galaxy correlation functions in halotools using tabulated values. It is specifically intended for Markov chain monte carlo (MCMC) exploration of the galaxy-halo connection. It implements the method described in [Zheng & Guo (2016)](https://doi.org/10.1093/mnras/stw523) of tabulating correlation functions that only need to be convolved with the mean halo occupation to obtain the full correlation function of galaxies. ## Installation -The package can be installed via pip using the following command. +The package can be installed via pip. ``` -pip install git+https://github.com/johannesulf/TabCorr.git +pip install tabcorr ``` ## Usage @@ -93,13 +86,10 @@ plt.ylabel(r'$w_p \ [h^{-1} \ \mathrm{Mpc}]$') plt.tight_layout(pad=0.3) plt.savefig('wp_vs_logm1.png', dpi=300) plt.close() - ``` -The above code will generate the following figures. - -![wp_decomposition](scripts/wp_decomposition.png) -![wp_vs_logm1](scripts/wp_vs_logm1.png) +![](scripts/wp_decomposition.png) +![](scripts/wp_vs_logm1.png) ## Author diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..eb96efe --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,18 @@ +[build-system] +requires = ["flit_core >=3.2,<4"] +build-backend = "flit_core.buildapi" + +[project] +name = "tabcorr" +authors = [{name = "Johannes U. Lange", email = "julange.astro@pm.me"}] +readme = "README.md" +license = {file = "LICENSE"} +classifiers = ["License :: OSI Approved :: MIT License"] +dynamic = ["version", "description"] +dependencies = ["numpy", "scipy", "astropy", "h5py", "tqdm"] + +[project.urls] +Home = "https://johannesulf.github.io/" + +[tool.flit.sdist] +exclude = ["scripts", ".gitignore"] diff --git a/scripts/download_snapshot.py b/scripts/download_snapshot.py new file mode 100644 index 0000000..3f947a0 --- /dev/null +++ b/scripts/download_snapshot.py @@ -0,0 +1,259 @@ +import os +import io +import struct +import requests +import argparse +import numpy as np +from tqdm import tqdm +from astropy.table import Table +from tabcorr import database +from collections import namedtuple + + +def read_gadget_snapshot(bstream, read_pos=False, read_vel=False, + read_id=False, read_mass=False, print_header=False, + single_type=-1, lgadget=False): + """ + This function reads the Gadget-2 snapshot file. + + This is a modified version of the function readGadgetSnapshot by Yao-Yuan + Mao licensed under the MIT License. (https://bitbucket.org/yymao/helpers) + + Parameters + ---------- + bstream : io.BytesIO + Binary stream of the gadget snapshot file. + read_pos : bool, optional + Whether to read the positions or not. Default is false. + read_vel : bool, optional + Whether to read the velocities or not. Default is false. + read_id : bool, optional + Whether to read the particle IDs or not. Default is false. + read_mass : bool, optional + Whether to read the masses or not. Default is false. + print_header : bool, optional + Whether to print out the header or not. Default is false. + single_type : int, optional + Set to -1 (default) to read in all particle types. + Set to 0--5 to read in only the corresponding particle type. + lgadget : bool, optional + Set to True if the particle file comes from l-gadget. + Default is false. + + Returns + ------- + ret : tuple + A tuple of the requested data. + The first item in the returned tuple is always the header. + The header is in the gadget_header namedtuple format. + + """ + gadget_header_fmt = '6I6dddii6Iiiddddii6Ii' + + gadget_header = namedtuple( + 'gadget_header', 'npart mass time redshift flag_sfr flag_feedback ' + + 'npartTotal flag_cooling num_files BoxSize Omega0 OmegaLambda ' + + 'HubbleParam flag_age flag_metals NallHW flag_entr_ics') + + blocks_to_read = (read_pos, read_vel, read_id, read_mass) + ret = [] + + bstream.seek(4, 1) + h = list(struct.unpack(gadget_header_fmt, + bstream.read(struct.calcsize(gadget_header_fmt)))) + if lgadget: + h[30] = 0 + h[31] = h[18] + h[18] = 0 + single_type = 1 + h = tuple(h) + header = gadget_header._make( + (h[0:6],) + (h[6:12],) + h[12:16] + (h[16:22],) + h[22:30] + + (h[30:36],) + h[36:]) + if print_header: + print(header) + if not any(blocks_to_read): + return header + ret.append(header) + bstream.seek(256 - struct.calcsize(gadget_header_fmt), 1) + bstream.seek(4, 1) + + mass_npart = [0 if m else n for m, n in zip(header.mass, header.npart)] + if single_type not in set(range(6)): + single_type = -1 + + for i, b in enumerate(blocks_to_read): + fmt = np.dtype(np.float32) + fmt_64 = np.dtype(np.float64) + item_per_part = 1 + npart = header.npart + + if i < 2: + item_per_part = 3 + elif i == 2: + fmt = np.dtype(np.uint32) + fmt_64 = np.dtype(np.uint64) + elif i == 3: + if sum(mass_npart) == 0: + ret.append(np.array([], fmt)) + break + npart = mass_npart + + size_check = struct.unpack('I', bstream.read(4))[0] + + block_item_size = item_per_part*sum(npart) + if size_check != block_item_size*fmt.itemsize: + fmt = fmt_64 + if size_check != block_item_size*fmt.itemsize: + raise ValueError('Invalid block size in file!') + size_per_part = item_per_part*fmt.itemsize + # + if not b: + bstream.seek(sum(npart)*size_per_part, 1) + else: + if single_type > -1: + bstream.seek(sum(npart[:single_type])*size_per_part, 1) + npart_this = npart[single_type] + else: + npart_this = sum(npart) + data = np.frombuffer(bstream.read(npart_this*size_per_part), fmt) + if item_per_part > 1: + data.shape = (npart_this, item_per_part) + ret.append(data) + if not any(blocks_to_read[i+1:]): + break + if single_type > -1: + bstream.seek(sum(npart[single_type+1:])*size_per_part, 1) + bstream.seek(4, 1) + + return tuple(ret) + + +def download_aemulus_alpha_halos(simulation, redshift): + + try: + username = os.environ['AEMULUS_USERNAME'] + password = os.environ['AEMULUS_PASSWORD'] + except KeyError: + raise RuntimeError("Set the AEMULUS_USERNAME and AEMULUS_PASSWORD " + + "environment variables.") + + scale_factor_snapshots = np.array([0.25, 0.333333, 0.5, 0.540541, 0.588235, + 0.645161, 0.714286, 0.8, 0.909091, 1.0]) + redshift_snapshots = 1 / scale_factor_snapshots - 1 + + try: + assert np.amin(np.abs(redshift_snapshots - redshift)) < 0.005 + except AssertionError: + raise ValueError('No snapshot for redshift {:.2f}.'.format(redshift)) + + snapnum = np.argmin(np.abs(redshift_snapshots - redshift)) + + url = ("https://www.slac.stanford.edu/~jderose/aemulus/phase1/" + + "{}/halos/m200b/outbgc2_{}.list".format(simulation, snapnum)) + + r = requests.get(url, auth=requests.auth.HTTPBasicAuth(username, password)) + halos = Table.read(io.BytesIO(r.content), format='ascii', delimiter=' ') + + url = url.replace('outbgc2', 'out') + r = requests.get(url, auth=requests.auth.HTTPBasicAuth(username, password)) + halos['halo_rs'] = Table.read( + io.BytesIO(r.content), format='ascii', delimiter=' ')['col7'] / 1e3 + halos['R200b'] /= 1e3 + + halos.rename_column('M200b', 'halo_m200m') + halos.rename_column('R200b', 'halo_r200m') + halos.rename_column('Vmax', 'halo_vmax') + for coordinate in ['x', 'y', 'z', 'vx', 'vy', 'vz']: + halos.rename_column(coordinate.upper(), 'halo_{}'.format(coordinate)) + + halos = halos[halos['Parent_ID'] == -1] + + halos.keep_columns([col for col in halos.colnames if col[:5] == 'halo_']) + + return halos + + +def download_aemulus_alpha_particles(simulation, redshift): + + try: + username = os.environ['AEMULUS_USERNAME'] + password = os.environ['AEMULUS_PASSWORD'] + except KeyError: + raise RuntimeError("Set the AEMULUS_USERNAME and AEMULUS_PASSWORD " + + "environment variables.") + + scale_factor_snapshots = np.array([0.25, 0.333333, 0.5, 0.540541, 0.588235, + 0.645161, 0.714286, 0.8, 0.909091, 1.0]) + redshift_snapshots = 1 / scale_factor_snapshots - 1 + + try: + assert np.amin(np.abs(redshift_snapshots - redshift)) < 0.005 + except AssertionError: + raise ValueError('No snapshot for redshift {:.2f}.'.format(redshift)) + + snapnum = np.argmin(np.abs(redshift_snapshots - redshift)) + + ptcls = np.zeros((0, 3)) + + for chunk in tqdm(range(512)): + url = ("https://www.slac.stanford.edu/~jderose/aemulus/phase1/" + + "{}/output/snapdir_{:03}/snapshot_{:03}.{}".format( + simulation, snapnum, snapnum, chunk)) + r = requests.get( + url, auth=requests.auth.HTTPBasicAuth(username, password)) + ptcls_tmp = read_gadget_snapshot( + io.BytesIO(r.content), read_pos=True)[1] + ptcls_tmp = ptcls_tmp[np.random.random(len(ptcls_tmp)) < 0.01] + ptcls = np.vstack([ptcls, ptcls_tmp]) + + return Table([ptcls[:, 0], ptcls[:, 1], ptcls[:, 2]], + names=('x', 'y', 'z')) + + +def main(): + + parser = argparse.ArgumentParser( + description='Download and convert an Aemulus Alpha simulation.') + parser.add_argument('suite', help='simulation suite', + choices=['AemulusAlpha']) + parser.add_argument('redshift', help='simulation redshift', type=float) + parser.add_argument('--cosmo', help='simulation cosmology, default is 0', + type=int, default=0) + parser.add_argument('--phase', help='simulation phase, default is 0', + type=int, default=0) + parser.add_argument('--config', default=None, + help='simulation configuration to assume') + parser.add_argument('--particles', action='store_true', + help='whether to download particles instead of halos') + + args = parser.parse_args() + + name = database.simulation_name( + args.suite, i_cosmo=args.cosmo, i_phase=args.phase, config=args.config) + + print('Downloading data for {} at z={:.2f}...'.format(name, args.redshift)) + + path = database.simulation_snapshot_directory( + args.suite, args.redshift, i_cosmo=args.cosmo, i_phase=args.phase, + config=args.config) + if not os.path.isdir(path): + os.makedirs(path) + + if not args.particles: + subpath = 'halos' + if args.suite == 'AemulusAlpha': + data = download_aemulus_alpha_halos(name, args.redshift) + else: + subpath = 'particles' + if args.suite == 'AemulusAlpha': + data = download_aemulus_alpha_particles(name, args.redshift) + + print('Writing results to {}.'.format(os.path.join(path, 'snapshot.hdf5'))) + data.write(os.path.join(path, 'snapshot.hdf5'), path=subpath, + overwrite=True, append=True) + print('Done!') + + +if __name__ == "__main__": + main() diff --git a/scripts/tabulate_snapshot.py b/scripts/tabulate_snapshot.py new file mode 100644 index 0000000..4147ce4 --- /dev/null +++ b/scripts/tabulate_snapshot.py @@ -0,0 +1,235 @@ +import os +import copy +import argparse +import numpy as np +import multiprocessing +import astropy.units as u +from astropy.table import Table +from astropy.constants import G +from tabcorr import TabCorr, database +from tabcorr.corrfunc import wp, s_mu_tpcf +from halotools.sim_manager import UserSuppliedHaloCatalog +from halotools.mock_observables import tpcf_multipole, mean_delta_sigma +from halotools.empirical_models import TrivialPhaseSpace, BiasedNFWPhaseSpace + + +def read_simulation_snapshot( + suite, redshift, i_cosmo=0, i_phase=None, config=None): + name = database.simulation_name( + suite, i_cosmo=i_cosmo, i_phase=i_phase, config=config) + directory = database.simulation_snapshot_directory( + suite, redshift, i_cosmo=i_cosmo, i_phase=i_phase, config=config) + halos = Table.read(os.path.join(directory, 'snapshot.hdf5'), + path='halos') + try: + ptcls = Table.read(os.path.join(directory, 'snapshot.hdf5'), + path='particles') + except OSError: + ptcls = None + cosmology = database.cosmology(suite, i_cosmo=i_cosmo) + + if suite == 'AbacusSummit': + mdef = '{:.0f}m'.format(halos.meta['SODensityL1']) + lbox = halos.meta['BoxSize'] + particle_mass = halos.meta['ParticleMassHMsun'] + halo_vmax = None + else: + mdef = '200m' + lbox = 1050 + particle_mass = 3.51e10 * cosmology.Om0 / 0.3 + halo_vmax = halos['halo_vmax'] + + return UserSuppliedHaloCatalog( + redshift=redshift, Lbox=lbox, particle_mass=particle_mass, + simname=name, halo_x=halos['halo_x'], halo_y=halos['halo_y'], + halo_z=halos['halo_z'], halo_vx=halos['halo_vx'], + halo_vy=halos['halo_vy'], halo_vz=halos['halo_vz'], + halo_id=np.arange(len(halos)), halo_pid=np.repeat(-1, len(halos)), + halo_upid=np.repeat(-1, len(halos)), + halo_nfw_conc=halos['halo_r{}'.format(mdef)] / halos['halo_rs'], + halo_mvir=halos['halo_m{}'.format(mdef)], + halo_rvir=halos['halo_r{}'.format(mdef)] * 1e-9, + halo_hostid=np.arange(len(halos)), cosmology=cosmology, + halo_vmax=halo_vmax, + **{'halo_m{}'.format(mdef): halos['halo_m{}'.format(mdef)], + 'halo_r{}'.format(mdef): halos['halo_r{}'.format(mdef)]}), ptcls + + +class ScaledBiasedNFWPhaseSpace(BiasedNFWPhaseSpace): + + def __init__(self, profile_integration_tol=1e-5, **kwargs): + BiasedNFWPhaseSpace.__init__( + self, profile_integration_tol=profile_integration_tol, **kwargs) + self.param_dict['alpha_s'] = 1.0 + + def _vrad_disp_from_lookup(self, scaled_radius, *concentration_array, + **kwargs): + return (BiasedNFWPhaseSpace._vrad_disp_from_lookup( + self, scaled_radius, *concentration_array, **kwargs) * + self.param_dict['alpha_s']) + + +class CentralVelocitBiasPhaseSpace(TrivialPhaseSpace): + + def __init__(self, mdef='200m', **kwargs): + TrivialPhaseSpace.__init__(self, **kwargs) + self.param_dict['alpha_c'] = 0.0 + self.mdef = mdef + + def assign_phase_space(self, table, **kwargs): + TrivialPhaseSpace.assign_phase_space(self, table, **kwargs) + vscale = np.sqrt( + G * table['halo_m' + self.mdef] * u.Msun / ( + table['halo_r' + self.mdef] * u.Mpc / (1 + self.redshift))).to( + u.km / u.s).value + for key in ['vx', 'vy', 'vz']: + table[key][:] += (vscale * np.random.normal(size=len(table)) * + self.param_dict['alpha_c'] / np.sqrt(3.0)) + + +def tabcorr_s_mu_to_multipole(halotab_s_mu, mu_bins, order): + halotab_mult = copy.deepcopy(halotab_s_mu) + halotab_mult.tpcf_shape = (halotab_s_mu.tpcf_shape[0], ) + halotab_mult.tpcf_matrix = np.zeros( + (halotab_s_mu.tpcf_shape[0], halotab_s_mu.tpcf_matrix.shape[1])) + + for i in range(halotab_s_mu.tpcf_matrix.shape[1]): + halotab_mult.tpcf_matrix[:, i] = tpcf_multipole( + halotab_s_mu.tpcf_matrix[:, i].reshape( + halotab_s_mu.tpcf_shape), mu_bins, order=order) + + return halotab_mult + + +def main(): + + parser = argparse.ArgumentParser( + description='Tabulate halo correlation functions.') + parser.add_argument('suite', help='simulation suite', + choices=['AemulusAlpha']) + parser.add_argument('redshift', help='simulation redshift', type=float) + parser.add_argument('--cosmo', help='simulation cosmology, default is 0', + type=int, default=0) + parser.add_argument('--phase', help='simulation phase, default is 0', + type=int, default=0) + parser.add_argument('--sim_config', default=None, + help='simulation configuration to assume') + parser.add_argument('--tab_config', default='default', + help='tabulation configuration to assume') + parser.add_argument('--tpcf', default='xi', choices=['xi', 'wp', 'ds'], + help='TPCF to tabulate') + + args = parser.parse_args() + + config = database.configuration(args.tab_config) + + halocat, ptcls = read_simulation_snapshot( + args.suite, args.redshift, i_cosmo=args.cosmo, i_phase=args.phase, + config=args.sim_config) + + for key in halocat.halo_table.colnames: + if key[:6] == 'halo_m' and key[-1] == 'm': + mdef = key[6:] + + if args.tpcf == 'wp' and config['pi_max'] >= 80: + config['alpha_c_bins'] = [0.0] + + if args.tpcf == 'ds': + config['alpha_c_bins'] = [0.0] + config['alpha_s_bins'] = [1.0] + + path = os.path.join(database.simulation_snapshot_directory( + args.suite, args.redshift, i_cosmo=args.cosmo, i_phase=args.phase, + config=args.sim_config), args.tab_config) + if not os.path.isdir(path): + os.makedirs(path) + + phase_space_grid = np.array(np.meshgrid( + config['alpha_c_bins'], config['alpha_s_bins'], + config['conc_gal_bias_bins'])).T.reshape(-1, 3) + table = Table() + table['alpha_c'] = phase_space_grid[:, 0] + table['alpha_s'] = phase_space_grid[:, 1] + table['conc_gal_bias'] = phase_space_grid[:, 2] + table.write(os.path.join(path, args.tpcf + '_grid.csv'), overwrite=True) + + for i, (alpha_c, alpha_s, conc_gal_bias) in enumerate(phase_space_grid): + + cens_prof_model = CentralVelocitBiasPhaseSpace( + redshift=halocat.redshift, mdef=mdef) + cens_prof_model.param_dict['alpha_c'] = alpha_c + + sats_prof_model = ScaledBiasedNFWPhaseSpace( + mdef=mdef, redshift=halocat.redshift, cosmology=halocat.cosmology, + concentration_bins=np.linspace(2.163, 20, 100), + conc_gal_bias_bins=np.array([conc_gal_bias])) + sats_prof_model.param_dict['alpha_s'] = alpha_s + + if args.tpcf == 'ds': + prim_haloprop_bins = 300 + mode = 'cross' + else: + prim_haloprop_bins = 30 + mode = 'auto' + + prim_haloprop_key = 'halo_m' + mdef + + if args.suite == 'AbacusSummit': + sec_haloprop_key = 'halo_nfw_conc' + else: + sec_haloprop_key = 'halo_vmax' + + sec_haloprop_percentile_bins = 0.5 + + if args.suite == 'AbacusSummit': + num_ptcl_requirement = 299 + else: + num_ptcl_requirement = 99 + + kwargs = {'mode': mode, 'cens_prof_model': cens_prof_model, + 'sats_prof_model': sats_prof_model, 'verbose': False, + 'num_threads': multiprocessing.cpu_count(), + 'sats_per_prim_haloprop': config['sats_per_prim_haloprop'], + 'project_xyz': True, + 'prim_haloprop_bins': prim_haloprop_bins, + 'prim_haloprop_key': prim_haloprop_key, + 'sec_haloprop_key': sec_haloprop_key, + 'sec_haloprop_percentile_bins': sec_haloprop_percentile_bins, + 'cosmology_obs': config['cosmo_obs'], + 'Num_ptcl_requirement': num_ptcl_requirement} + + if args.tpcf == 'xi': + halotab_s_mu = TabCorr.tabulate( + halocat, s_mu_tpcf, config['s_bins'], config['mu_bins'], + **kwargs) + for order in [0, 2, 4]: + halotab_multipole = tabcorr_s_mu_to_multipole( + halotab_s_mu, config['mu_bins'], order) + halotab_multipole.write(os.path.join( + path, 'xi{}_{}.hdf5'.format(order, i)), overwrite=True) + + elif args.tpcf == 'wp': + halotab = TabCorr.tabulate( + halocat, wp, config['rp_wp_bins'], config['pi_max'], **kwargs) + halotab.write(os.path.join(path, 'wp_{}.hdf5'.format(i)), + overwrite=True) + + elif args.tpcf == 'ds': + + ptcl_pos = np.vstack([ptcls['x'], ptcls['y'], ptcls['z']]).T + + if args.suite == 'AemulusAlpha': + n_ptcl_tot = 1400**3 + + downsampling_factor = n_ptcl_tot / float(len(ptcl_pos)) + ptcl_mass = halocat.particle_mass * downsampling_factor + + halotab = TabCorr.tabulate( + halocat, mean_delta_sigma, ptcl_pos, ptcl_mass, + config['rp_ds_bins'], **kwargs) + halotab.write(os.path.join(path, 'ds_{}.hdf5'.format(i)), + overwrite=True) + + +if __name__ == "__main__": + main() diff --git a/setup.py b/setup.py deleted file mode 100644 index c1b1b71..0000000 --- a/setup.py +++ /dev/null @@ -1,12 +0,0 @@ -from setuptools import setup -from tabcorr import __version__ - -setup(name='tabcorr', - version=__version__, - description='Tabulated Correlation functions for halotools', - url='https://github.com/johannesulf/TabCorr', - author='Johannes U. Lange', - author_email='julange.astro@pm.me', - packages=['tabcorr'], - install_requires=['numpy', 'scipy', 'astropy', 'h5py', 'tqdm'], - zip_safe=False) diff --git a/tabcorr/__init__.py b/tabcorr/__init__.py index 3f87127..e2f19df 100644 --- a/tabcorr/__init__.py +++ b/tabcorr/__init__.py @@ -1,5 +1,9 @@ +"""Tabulated Correlation Functions""" + from .tabcorr import TabCorr from .interpolate import Interpolator +from . import corrfunc +from . import database -__version__ = '0.9.1' -__all__ = ["TabCorr", "Interpolator"] +__version__ = '1.0.0' +__all__ = ["TabCorr", "Interpolator", "corrfunc", "database"] diff --git a/tabcorr/aa_cosmos.txt b/tabcorr/aa_cosmos.txt new file mode 100644 index 0000000..97d4f3e --- /dev/null +++ b/tabcorr/aa_cosmos.txt @@ -0,0 +1,41 @@ +#boxnum ombh2 omch2 w0 ns ln10As H0 Neff sigma8 +000 2.268325e-02 1.140598e-01 -8.165972e-01 9.755890e-01 3.092918e+00 6.336569e+01 2.918750e+00 7.730715e-01 +001 2.245071e-02 1.172576e-01 -1.134440e+00 9.765227e-01 3.149872e+00 7.309972e+01 3.173750e+00 8.954628e-01 +002 2.299820e-02 1.087487e-01 -6.846595e-01 9.974496e-01 3.093773e+00 6.370848e+01 3.258750e+00 6.929467e-01 +003 2.267590e-02 1.123064e-01 -7.437379e-01 9.480812e-01 3.000900e+00 6.404189e+01 3.556250e+00 6.715524e-01 +004 2.208455e-02 1.062517e-01 -7.666051e-01 9.650928e-01 3.118565e+00 6.504913e+01 2.663750e+00 7.490580e-01 +005 2.066455e-02 1.295035e-01 -1.326359e+00 9.278462e-01 3.023645e+00 7.275418e+01 2.961250e+00 9.339844e-01 +006 2.294146e-02 1.115033e-01 -7.103081e-01 9.706424e-01 3.015966e+00 6.269873e+01 2.706250e+00 7.103632e-01 +007 2.281365e-02 1.195969e-01 -8.666033e-01 9.662983e-01 3.162138e+00 6.437164e+01 3.938750e+00 7.778608e-01 +008 2.067695e-02 1.238351e-01 -1.163731e+00 9.491313e-01 3.147130e+00 6.939996e+01 3.598750e+00 8.919908e-01 +009 2.134390e-02 1.158135e-01 -8.309330e-01 9.475234e-01 3.071769e+00 6.235637e+01 3.896250e+00 7.220269e-01 +010 2.189334e-02 1.289506e-01 -1.241154e+00 9.610168e-01 3.050463e+00 7.208599e+01 4.236250e+00 8.497577e-01 +011 2.262122e-02 1.089506e-01 -8.609277e-01 9.960324e-01 3.158114e+00 6.773404e+01 2.833750e+00 8.055318e-01 +012 2.250605e-02 1.167658e-01 -8.794716e-01 9.539512e-01 3.048009e+00 6.537987e+01 2.876250e+00 7.851298e-01 +013 2.189555e-02 1.171763e-01 -1.119505e+00 9.788153e-01 3.067809e+00 7.108341e+01 3.003750e+00 8.660544e-01 +014 2.264907e-02 1.271080e-01 -1.116826e+00 9.723654e-01 3.093885e+00 6.873003e+01 2.748750e+00 9.232811e-01 +015 2.146760e-02 1.285016e-01 -1.302908e+00 9.335675e-01 3.094304e+00 7.410001e+01 3.726250e+00 9.078274e-01 +016 2.177406e-02 1.206762e-01 -1.131299e+00 9.662156e-01 3.014076e+00 7.007432e+01 3.768750e+00 8.069646e-01 +017 2.228121e-02 1.194239e-01 -1.248070e+00 9.519719e-01 3.035022e+00 7.443912e+01 3.216250e+00 8.695022e-01 +018 2.287409e-02 1.157403e-01 -1.032177e+00 9.533103e-01 3.020205e+00 7.075163e+01 4.278750e+00 7.356345e-01 +019 2.238260e-02 1.132626e-01 -1.091607e+00 9.672909e-01 3.095691e+00 7.242955e+01 3.683750e+00 8.083719e-01 +020 2.225371e-02 1.224880e-01 -9.898374e-01 9.528876e-01 3.119866e+00 6.705524e+01 3.386250e+00 8.378528e-01 +021 2.361513e-02 1.171690e-01 -8.658316e-01 9.757853e-01 3.131585e+00 6.638900e+01 3.853750e+00 7.617592e-01 +022 2.153845e-02 1.210090e-01 -1.031605e+00 9.585567e-01 3.072322e+00 6.806190e+01 2.621250e+00 8.804510e-01 +023 2.272630e-02 1.012181e-01 -5.658486e-01 9.745746e-01 3.019204e+00 6.203350e+01 3.471250e+00 5.746315e-01 +024 2.245361e-02 1.102766e-01 -7.614561e-01 9.589180e-01 3.143937e+00 6.303209e+01 4.151250e+00 6.918467e-01 +025 2.092906e-02 1.170569e-01 -9.478782e-01 9.344636e-01 3.036794e+00 6.571267e+01 3.088750e+00 7.900302e-01 +026 2.236205e-02 1.191811e-01 -1.125103e+00 9.442601e-01 3.127891e+00 7.175536e+01 2.791250e+00 9.043988e-01 +027 2.135801e-02 1.134255e-01 -9.645883e-01 9.663881e-01 3.014713e+00 6.739297e+01 4.023750e+00 7.285257e-01 +028 2.174718e-02 1.317768e-01 -1.399921e+00 9.585702e-01 3.146559e+00 7.476752e+01 3.811250e+00 9.644551e-01 +029 2.232480e-02 1.289151e-01 -1.236071e+00 9.401456e-01 3.159051e+00 7.141356e+01 3.428750e+00 9.344431e-01 +030 2.192792e-02 1.239137e-01 -1.223681e+00 9.552358e-01 3.118282e+00 7.343152e+01 4.066250e+00 8.682105e-01 +031 2.123964e-02 1.275787e-01 -1.381645e+00 9.560725e-01 3.076413e+00 7.375997e+01 3.343750e+00 9.417921e-01 +032 2.249688e-02 1.128450e-01 -9.257578e-01 9.494842e-01 3.042810e+00 6.840342e+01 3.981250e+00 7.208264e-01 +033 2.337319e-02 1.149703e-01 -8.752207e-01 9.891731e-01 3.149163e+00 6.605231e+01 3.641250e+00 7.771575e-01 +034 2.275573e-02 1.222293e-01 -1.032290e+00 9.500281e-01 3.106940e+00 6.907038e+01 3.131250e+00 8.576372e-01 +035 2.341077e-02 1.076406e-01 -6.128197e-01 9.955676e-01 3.140289e+00 6.169472e+01 3.046250e+00 6.845285e-01 +036 2.199894e-02 1.212967e-01 -1.108222e+00 9.673610e-01 3.179424e+00 7.041129e+01 3.301250e+00 9.050606e-01 +037 2.287447e-02 1.097440e-01 -8.487455e-01 9.776227e-01 3.072287e+00 6.672620e+01 3.513750e+00 7.255155e-01 +038 2.371239e-02 1.149747e-01 -9.550075e-01 9.765907e-01 3.053508e+00 6.974680e+01 4.108750e+00 7.368773e-01 +039 2.174072e-02 1.200700e-01 -9.410529e-01 9.602072e-01 3.092704e+00 6.470430e+01 4.193750e+00 7.613921e-01 diff --git a/tabcorr/aa_test_cosmos.txt b/tabcorr/aa_test_cosmos.txt new file mode 100644 index 0000000..8e20043 --- /dev/null +++ b/tabcorr/aa_test_cosmos.txt @@ -0,0 +1,8 @@ +#boxnum ombh2 omch2 w0 ns ln10As H0 Neff sigma8 +000 2.32629e-02 1.07830e-01 -7.26513e-01 9.80515e-01 3.03895e+00 6.32317e+01 2.95000e+00 6.944280e-01 +001 2.27629e-02 1.12830e-01 -8.61513e-01 9.71515e-01 3.06395e+00 6.57317e+01 3.20000e+00 7.542319e-01 +002 2.22629e-02 1.17830e-01 -9.96513e-01 9.62515e-01 3.08895e+00 6.82317e+01 3.45000e+00 8.087561e-01 +003 2.17629e-02 1.22830e-01 -1.13151e+00 9.53515e-01 3.11395e+00 7.07317e+01 3.70000e+00 8.596132e-01 +004 2.12629e-02 1.27830e-01 -1.26651e+00 9.44515e-01 3.13895e+00 7.32317e+01 3.95000e+00 9.077866e-01 +005 2.17629e-02 1.15330e-01 -1.08911e+00 9.51404e-01 3.11895e+00 6.97317e+01 3.70000e+00 8.177121e-01 +006 2.27629e-02 1.20330e-01 -9.03920e-01 9.73626e-01 3.05895e+00 6.67317e+01 3.20000e+00 7.981356e-01 diff --git a/tabcorr/as_cosmos.csv b/tabcorr/as_cosmos.csv new file mode 100644 index 0000000..7338911 --- /dev/null +++ b/tabcorr/as_cosmos.csv @@ -0,0 +1,100 @@ +root , notes , omega_b , omega_cdm , h , A_s , n_s , alpha_s , N_ur , N_ncdm , omega_ncdm , w0_fld , wa_fld , sigma8_m , sigma8_cb +abacus_cosm000 ,"Baseline LCDM, Planck 2018 base_plikHM_TTTEEE_lowl_lowE_lensing mean", 0.02237 , 0.1200 , 0.6736 , 2.0830e-9 , 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.807952 , 0.811355 +abacus_cosm001 ,"WMAP9+ACT+SPT LCDM, Calabrese++ 2017 ", 0.02242 , 0.1134 , 0.7030 , 2.0376e-9 , 0.9638 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.776779 , 0.780222 +abacus_cosm002 ,"wCDM with thawing model w0 = -0.7, wa = -0.5 ", 0.02237 , 0.1200 , 0.6278 , 2.3140e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -0.7 , -0.5 , 0.808189 , 0.811577 +abacus_cosm003 ,"Neff=3.70, from base_nnu_plikHM_TT_lowl_lowE_Riess18_post_BAO ", 0.02260 , 0.1291 , 0.7160 , 2.2438e-9 , 0.9876 , 0.0 , 2.6868 , 1 , 0.00064420 , -1.0 , 0.0 , 0.855190 , 0.858583 +abacus_cosm004 ,"Low sigma8_matter = 0.75, otherwise Baseline LCDM ", 0.02237 , 0.1200 , 0.6736 , 1.7949e-9 , 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.749999 , 0.753159 +abacus_cosm009 ,"Baseline LCDM with massless neutrinos matching omega_cb & sigma8_cb ", 0.02237 , 0.1200 , 0.6736 , 2.0417e-9 , 0.9649 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.811362 ,0.811362 +abacus_cosm010 ,"AbacusCosmos Planck LCDM cosmology ", 0.02222 , 0.1199 , 0.6726 , 2.100e-9 , 0.9652 , 0.0 , 3.04 , 0 , 0.0 , -1.0 , 0.0 , 0.823630 ,0.823630 +abacus_cosm011 ,"AbacusCosmos Planck LCDM cosmology +10% in sigma8 ", 0.02222 , 0.1199 , 0.6726 , 2.541e-9 , 0.9652 , 0.0 , 3.04 , 0 , 0.0 , -1.0 , 0.0 , 0.905993 ,0.905993 +abacus_cosm012 ,"Euclid Flagship1 LCDM, sigma8_m=0.8279 ", 0.02200 , 0.1212 , 0.6700 , 2.1000e-9 , 0.9600 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.827899 ,0.827899 +abacus_cosm013 ,"Euclid Flagship2 LCDM, sigma8_m=0.813715, sigma8_cb=0.817135 ", 0.02200 , 0.1206 , 0.6700 , 2.1000e-9 , 0.9600 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.813715 , 0.817135 +abacus_cosm014 ,"OuterRim LCDM, sigma8=0.80 ", 0.02258 , 0.1109 , 0.7100 , 2.1591e-9 , 0.9630 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.800000 ,0.800000 +abacus_cosm015 ,"DarkSky LCDM, sigma8=0.835 ", 0.02215 , 0.1175 , 0.6880 , 2.1852e-9 , 0.9688 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.835005 ,0.835005 +abacus_cosm016 ,"Horizon Run 4 LCDM, sigma8=0.7937 ", 0.02281 , 0.1120 , 0.7200 , 2.0996e-9 , 0.9600 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.793693 ,0.793693 +abacus_cosm017 ,"IllustrisTNG LCDM, sigma8=0.8159 ", 0.02230 , 0.1194 , 0.6774 , 2.0671e-9 , 0.9667 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.815903 ,0.815903 +abacus_cosm018 ,"MultiDark Planck LCDM, sigma8=0.8228 ", 0.02214 , 0.1189 , 0.6777 , 2.1022e-9 , 0.9600 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.819708 ,0.819708 +abacus_cosm019 ,"Baseline LCDM with two 0.06 eV neutrinos ", 0.02237 , 0.1200 , 0.6683 , 2.1353e-09, 0.9649 , 0.0 , 1.0196 , 2 ,"0.00064420,0.00064420", -1.0 , 0.0 ,0.805050 , 0.811826 +abacus_cosm020 ,"Baseline LCDM w/ massless neutrinos matching theta_star & sigma8_cb ", 0.02237 , 0.1200 , 0.6790 , 2.0342e-09, 0.9649 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 , 0.811350 ,0.811350 +abacus_cosm021 ,"MassiveNUs massless base ", 0.02254 , 0.12446 , 0.700 , 2.1000e-09, 0.9700 , 0.0 , 3.046 , 0 , 0.0 , -1.0 , 0.0 ,0.849842 ,0.849842 +abacus_cosm022 ,"MassiveNUs with one 0.1 eV massive neutrino species ", 0.02254 , 0.12339 , 0.700 , 2.1000e-09, 0.9700 , 0.0 , 2.0328 , 1 , 0.001074 , -1.0 , 0.0 ,0.824961 , 0.830296 +abacus_cosm100 ,"Baseline +2% ln(omega_b) ", 0.02282 , 0.1200 , 0.6777 , 2.0934e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808181 , 0.811575 +abacus_cosm101 ,"Baseline -2% ln(omega_b) ", 0.02193 , 0.1200 , 0.6696 , 2.0751e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808156 , 0.811570 +abacus_cosm102 ,"Baseline +3.3% ln(omega_c) ", 0.02237 , 0.1240 , 0.6597 , 2.0205e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808270 , 0.811574 +abacus_cosm103 ,"Baseline -3.3% ln(omega_c) ", 0.02237 , 0.1161 , 0.6877 , 2.1541e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808075 , 0.811582 +abacus_cosm104 ,"Baseline +0.01 n_s ", 0.02237 , 0.1200 , 0.6736 , 2.0684e-09, 0.9749 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808166 , 0.811572 +abacus_cosm105 ,"Baseline -0.01 n_s ", 0.02237 , 0.1200 , 0.6736 , 2.0999e-09, 0.9549 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808177 , 0.811580 +abacus_cosm106 ,"Baseline +0.02 nrun ", 0.02237 , 0.1200 , 0.6736 , 2.0638e-09, 0.9649 , 0.02 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808181 , 0.811586 +abacus_cosm107 ,"Baseline -0.02 nrun ", 0.02237 , 0.1200 , 0.6736 , 2.1045e-09, 0.9649 , -0.02 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808168 , 0.811571 +abacus_cosm108 ,"Baseline +0.1 w0 ", 0.02237 , 0.1200 , 0.6444 , 2.2390e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -0.9 , 0.0 , 0.808177 , 0.811570 +abacus_cosm109 ,"Baseline -0.1 w0 ", 0.02237 , 0.1200 , 0.7037 , 1.9465e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.1 , 0.0 , 0.808161 , 0.811576 +abacus_cosm110 ,"Baseline +0.4 wa, -0.1 w0 ", 0.02237 , 0.1200 , 0.6698 , 2.1219e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.1 , 0.4 , 0.808170 , 0.811572 +abacus_cosm111 ,"Baseline -0.4 wa, +0.1 w0 ", 0.02237 , 0.1200 , 0.6752 , 2.0639e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -0.9 , -0.4 , 0.808179 , 0.811584 +abacus_cosm112 ,"Baseline +2% sigma8 ", 0.02237 , 0.1200 , 0.6736 , 2.1672e-9 , 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.824120 , 0.827591 +abacus_cosm113 ,"Baseline -2% sigma8 ", 0.02237 , 0.1200 , 0.6736 , 2.0021e-9 , 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.792107 , 0.795443 +abacus_cosm114 ,"Baseline +0.3 nnu, +3.3% ln(omega_c), +0.01 n_s ", 0.02237 , 0.1240 , 0.6947 , 2.0463e-09, 0.9749 , 0.0 , 2.3328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808245 , 0.811563 +abacus_cosm115 ,"Baseline -0.3 nnu, -3.3% ln(omega_c), -0.01 n_s ", 0.02237 , 0.1161 , 0.6517 , 2.1211e-09, 0.9549 , 0.0 , 1.7328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808089 , 0.811582 +abacus_cosm116 ,"Emulator grid around baseline cosmology ", 0.02237 , 0.1200 , 0.6736 , 2.3938e-09, 0.9649 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.866133 , 0.869782 +abacus_cosm117 ,"Baseline +0.83% ln(omega_c) ", 0.02237 , 0.1210 , 0.6701 , 2.0675e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808205 , 0.811584 +abacus_cosm118 ,"Baseline -0.83% ln(omega_c) ", 0.02237 , 0.1190 , 0.6771 , 2.1014e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808152 , 0.811582 +abacus_cosm119 ,"Baseline +0.003 n_s ", 0.02237 , 0.1200 , 0.6736 , 2.0794e-09, 0.9679 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808168 , 0.811573 +abacus_cosm120 ,"Baseline -0.003 n_s ", 0.02237 , 0.1200 , 0.6736 , 2.0889e-09, 0.9619 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.808181 , 0.811585 +abacus_cosm121 ,"Baseline +0.025 w0 ", 0.02237 , 0.1200 , 0.6662 , 2.1211e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -0.975 , 0.0 , 0.808169 , 0.811571 +abacus_cosm122 ,"Baseline -0.025 w0 ", 0.02237 , 0.1200 , 0.6810 , 2.0483e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.025 , 0.0 , 0.808170 , 0.811577 +abacus_cosm123 ,"Baseline +0.1 wa, -0.025 w0 ", 0.02237 , 0.1200 , 0.6729 , 2.0915e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.025 , 0.1 , 0.808176 , 0.811580 +abacus_cosm124 ,"Baseline -0.1 wa, +0.025 w0 ", 0.02237 , 0.1200 , 0.6742 , 2.0778e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -0.975 , -0.1 , 0.808175 , 0.811580 +abacus_cosm125 ,"Baseline +0.5% sigma8 ", 0.02237 , 0.1200 , 0.6736 , 2.1039e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.811995 , 0.815415 +abacus_cosm126 ,"Baseline -0.5% sigma8 ", 0.02237 , 0.1200 , 0.6736 , 2.0622e-09, 0.9649 , 0.0 , 2.0328 , 1 , 0.00064420 , -1.0 , 0.0 , 0.803908 , 0.807294 +abacus_cosm130 ,"Emulator grid around baseline cosmology ", 0.02237 , 0.1200 , 0.6736 , 1.6140e-09, 0.9649 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.711201 , 0.714197 +abacus_cosm131 ,"Emulator grid around baseline cosmology ", 0.02237 , 0.1086 , 0.7165 , 2.3146e-09, 0.9649 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.807866 , 0.811587 +abacus_cosm132 ,"Emulator grid around baseline cosmology ", 0.02237 , 0.1200 , 0.6736 , 2.1791e-09, 0.9049 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.808189 , 0.811584 +abacus_cosm133 ,"Emulator grid around baseline cosmology ", 0.02237 , 0.1200 , 0.6736 , 2.6144e-09, 0.9649 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.905163 , 0.908976 +abacus_cosm134 ,"Emulator grid around baseline cosmology ", 0.02237 , 0.1326 , 0.6319 , 1.9066e-09, 0.9649 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.808458 , 0.811566 +abacus_cosm135 ,"Emulator grid around baseline cosmology ", 0.02237 , 0.1200 , 0.6736 , 1.9904e-09, 1.0249 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.808160 , 0.811573 +abacus_cosm136 ,"Emulator grid around baseline cosmology ", 0.02073 , 0.1192 , 0.6618 , 2.5239e-09, 0.9303 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.882475 , 0.886249 +abacus_cosm137 ,"Emulator grid around baseline cosmology ", 0.02212 , 0.1271 , 0.6472 , 1.6540e-09, 0.9252 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.729693 , 0.732609 +abacus_cosm138 ,"Emulator grid around baseline cosmology ", 0.02108 , 0.1138 , 0.6847 , 2.2282e-09, 0.9723 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.821432 , 0.825095 +abacus_cosm139 ,"Emulator grid around baseline cosmology ", 0.02416 , 0.1128 , 0.7164 , 2.1681e-09, 0.9732 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.793003 , 0.796493 +abacus_cosm140 ,"Emulator grid around baseline cosmology ", 0.02096 , 0.1221 , 0.6536 , 1.8126e-09, 0.9893 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.772033 , 0.775267 +abacus_cosm141 ,"Emulator grid around baseline cosmology ", 0.02381 , 0.1272 , 0.6623 , 1.9945e-09, 0.9384 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.799575 , 0.802737 +abacus_cosm142 ,"Emulator grid around baseline cosmology ", 0.02287 , 0.1130 , 0.7038 , 1.6701e-09, 0.9927 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.707082 , 0.710219 +abacus_cosm143 ,"Emulator grid around baseline cosmology ", 0.02206 , 0.1278 , 0.6443 , 2.1084e-09, 0.9952 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.847522 , 0.850903 +abacus_cosm144 ,"Emulator grid around baseline cosmology ", 0.02210 , 0.1130 , 0.6968 , 2.7653e-09, 0.9279 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.891360 , 0.895324 +abacus_cosm145 ,"Emulator grid around baseline cosmology ", 0.02428 , 0.1186 , 0.6961 , 2.1634e-09, 0.9347 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.801404 , 0.804769 +abacus_cosm146 ,"Emulator grid around baseline cosmology ", 0.02097 , 0.1180 , 0.6682 , 1.9020e-09, 0.9351 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.000 , 0.000 , 0.762696 , 0.765984 +abacus_cosm147 ,"Emulator grid around baseline cosmology ", 0.02113 , 0.1215 , 0.6498 , 1.9081e-09, 0.9587 , 0.000 , 2.0328 , 1 , 0.00064420 , -0.809 , -0.628 , 0.777157 , 0.780416 +abacus_cosm148 ,"Emulator grid around baseline cosmology ", 0.02289 , 0.1201 , 0.6199 , 1.9629e-09, 0.9380 , 0.000 , 2.0328 , 1 , 0.00064420 , -0.925 , 0.393 , 0.714913 , 0.717888 +abacus_cosm149 ,"Emulator grid around baseline cosmology ", 0.02188 , 0.1200 , 0.6688 , 2.1521e-09, 0.9913 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.174 , 0.613 , 0.824854 , 0.828343 +abacus_cosm150 ,"Emulator grid around baseline cosmology ", 0.02248 , 0.1216 , 0.7224 , 2.5099e-09, 0.9355 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.090 , -0.373 , 0.937655 , 0.941570 +abacus_cosm151 ,"Emulator grid around baseline cosmology ", 0.02315 , 0.1254 , 0.6582 , 2.2686e-09, 0.9757 , 0.000 , 2.0328 , 1 , 0.00064420 , -0.826 , -0.615 , 0.860820 , 0.864288 +abacus_cosm152 ,"Emulator grid around baseline cosmology ", 0.02165 , 0.1148 , 0.6204 , 1.7965e-09, 0.9541 , 0.000 , 2.0328 , 1 , 0.00064420 , -0.732 , -0.187 , 0.677885 , 0.680849 +abacus_cosm153 ,"Emulator grid around baseline cosmology ", 0.02192 , 0.1199 , 0.6092 , 2.2946e-09, 0.9917 , 0.000 , 2.0328 , 1 , 0.00064420 , -0.735 , -0.172 , 0.794389 , 0.797729 +abacus_cosm154 ,"Emulator grid around baseline cosmology ", 0.02158 , 0.1148 , 0.7426 , 2.0689e-09, 0.9538 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.076 , -0.395 , 0.838698 , 0.842413 +abacus_cosm155 ,"Emulator grid around baseline cosmology ", 0.02369 , 0.1184 , 0.6247 , 2.0723e-09, 0.9713 , 0.000 , 2.0328 , 1 , 0.00064420 , -0.727 , -0.184 , 0.735302 , 0.738388 +abacus_cosm156 ,"Emulator grid around baseline cosmology ", 0.02202 , 0.1261 , 0.6499 , 1.9738e-09, 0.9678 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.176 , 0.621 , 0.801974 , 0.805210 +abacus_cosm157 ,"Emulator grid around baseline cosmology ", 0.02247 , 0.1214 , 0.6668 , 2.5108e-09, 0.9356 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.167 , 0.616 , 0.872315 , 0.875941 +abacus_cosm158 ,"Emulator grid around baseline cosmology ", 0.02201 , 0.1262 , 0.6958 , 1.8626e-09, 0.9667 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.082 , -0.392 , 0.829816 , 0.833179 +abacus_cosm159 ,"Emulator grid around baseline cosmology ", 0.02206 , 0.1261 , 0.5967 , 1.7977e-09, 0.9673 , 0.000 , 2.0328 , 1 , 0.00064420 , -0.726 , -0.220 , 0.718521 , 0.721402 +abacus_cosm160 ,"Emulator grid around baseline cosmology ", 0.02223 , 0.1185 , 0.7456 , 2.0963e-09, 0.9942 , 0.000 , 2.0328 , 1 , 0.00064420 , -1.271 , 0.217 , 0.876756 , 0.880523 +abacus_cosm161 ,"Emulator grid around baseline cosmology ", 0.02282 , 0.1274 , 0.7330 , 1.9688e-09, 0.9649 , 0.026 , 2.7283 , 1 , 0.00064420 , -1.000 , 0.000 , 0.793066 , 0.796247 +abacus_cosm162 ,"Emulator grid around baseline cosmology ", 0.02136 , 0.1048 , 0.6383 , 2.1592e-09, 0.9297 , -0.026 , 1.3521 , 1 , 0.00064420 , -1.000 , 0.000 , 0.779589 , 0.783285 +abacus_cosm163 ,"Emulator grid around baseline cosmology ", 0.02135 , 0.1290 , 0.7280 , 2.1689e-09, 0.9922 , -0.018 , 2.8575 , 1 , 0.00064420 , -1.000 , 0.000 , 0.838824 , 0.842187 +abacus_cosm164 ,"Emulator grid around baseline cosmology ", 0.02173 , 0.1091 , 0.6042 , 2.0003e-09, 0.9211 , 0.016 , 1.1769 , 1 , 0.00064420 , -1.000 , 0.000 , 0.774159 , 0.777681 +abacus_cosm165 ,"Emulator grid around baseline cosmology ", 0.02225 , 0.1278 , 0.7293 , 2.0676e-09, 1.0241 , 0.026 , 2.7558 , 1 , 0.00064420 , -1.000 , 0.000 , 0.835954 , 0.839319 +abacus_cosm166 ,"Emulator grid around baseline cosmology ", 0.02306 , 0.1151 , 0.6826 , 2.2779e-09, 0.9691 , 0.038 , 1.9089 , 1 , 0.00064420 , -1.000 , 0.000 , 0.837463 , 0.841106 +abacus_cosm167 ,"Emulator grid around baseline cosmology ", 0.02248 , 0.1400 , 0.7050 , 1.7321e-09, 0.9715 , -0.017 , 2.8889 , 1 , 0.00064420 , -1.000 , 0.000 , 0.768419 , 0.771255 +abacus_cosm168 ,"Emulator grid around baseline cosmology ", 0.02157 , 0.1084 , 0.6070 , 2.5378e-09, 0.9275 , 0.016 , 1.1911 , 1 , 0.00064420 , -1.000 , 0.000 , 0.871407 , 0.875401 +abacus_cosm169 ,"Emulator grid around baseline cosmology ", 0.02189 , 0.1222 , 0.6474 , 1.8067e-09, 0.9898 , 0.038 , 1.9109 , 1 , 0.00064420 , -1.000 , 0.000 , 0.777925 , 0.781158 +abacus_cosm170 ,"Emulator grid around baseline cosmology ", 0.02319 , 0.1076 , 0.6418 , 1.7994e-09, 0.9387 , -0.026 , 1.3411 , 1 , 0.00064420 , -1.000 , 0.000 , 0.716059 , 0.719332 +abacus_cosm171 ,"Emulator grid around baseline cosmology ", 0.02348 , 0.1272 , 0.6733 , 2.2585e-09, 0.9762 , -0.038 , 2.1543 , 1 , 0.00064420 , -1.000 , 0.000 , 0.852878 , 0.856267 +abacus_cosm172 ,"Emulator grid around baseline cosmology ", 0.02282 , 0.1032 , 0.6369 , 2.1259e-09, 0.9012 , 0.017 , 1.1954 , 1 , 0.00064420 , -1.000 , 0.000 , 0.765650 , 0.769278 +abacus_cosm173 ,"Emulator grid around baseline cosmology ", 0.02322 , 0.1161 , 0.6800 , 1.8904e-09, 0.9628 , 0.038 , 1.9049 , 1 , 0.00064420 , -1.000 , 0.000 , 0.763962 , 0.767255 +abacus_cosm174 ,"Emulator grid around baseline cosmology ", 0.02302 , 0.1064 , 0.6427 , 2.4889e-09, 0.9437 , -0.026 , 1.3242 , 1 , 0.00064420 , -1.000 , 0.000 , 0.840113 , 0.843996 +abacus_cosm175 ,"Emulator grid around baseline cosmology ", 0.02173 , 0.1336 , 0.6432 , 1.7429e-09, 0.9898 , -0.003 , 2.7456 , 1 , 0.00064420 , -0.911 , 0.350 , 0.708760 , 0.711484 +abacus_cosm176 ,"Emulator grid around baseline cosmology ", 0.02238 , 0.1199 , 0.6972 , 2.3120e-09, 0.9459 , 0.036 , 1.9043 , 1 , 0.00064420 , -1.222 , 0.339 , 0.892483 , 0.896318 +abacus_cosm177 ,"Emulator grid around baseline cosmology ", 0.02239 , 0.1344 , 0.6820 , 2.0641e-09, 1.0002 , 0.007 , 2.8643 , 1 , 0.00064420 , -0.757 , -0.443 , 0.806026 , 0.809110 +abacus_cosm178 ,"Emulator grid around baseline cosmology ", 0.02234 , 0.1204 , 0.6542 , 1.9665e-09, 0.9424 , 0.037 , 1.8987 , 1 , 0.00064420 , -0.874 , -0.455 , 0.791239 , 0.794555 +abacus_cosm179 ,"Emulator grid around baseline cosmology ", 0.02240 , 0.1067 , 0.5881 , 2.2365e-09, 0.9308 , -0.003 , 1.1884 , 1 , 0.00064420 , -0.755 , -0.435 , 0.775969 , 0.779541 +abacus_cosm180 ,"Emulator grid around baseline cosmology ", 0.02258 , 0.1209 , 0.7451 , 2.2790e-09, 0.9796 , -0.034 , 2.1722 , 1 , 0.00064420 , -1.108 , -0.274 , 0.894071 , 0.897833 +abacus_cosm181 ,"Emulator grid around baseline cosmology ", 0.02169 , 0.1112 , 0.5745 , 2.0651e-09, 0.9336 , -0.013 , 1.4059 , 1 , 0.00064420 , -0.910 , 0.350 , 0.730036 , 0.733292 + diff --git a/tabcorr/corrfunc.py b/tabcorr/corrfunc.py new file mode 100644 index 0000000..7c3a261 --- /dev/null +++ b/tabcorr/corrfunc.py @@ -0,0 +1,175 @@ +"""Module providing wrappers around Corrfunc functions.""" + +import numpy as np + + +def wp(sample1, rp_bins, pi_max, sample2=None, period=None, do_auto=True, + do_cross=False): + """Corrfunc implementation of `halotools.mock_observables.wp`. + + This function provides a wrapper around `corrfunc.theory.DDrppi` such that + when used by TabCorr behaves the same as `halotools.mock_observables.wp`. + However, not all functionality of `halotools.mock_observables.wp` is + implemented. + + Parameters + ---------- + sample1 : numpy.ndarray + Numpy array containing the positions of points. + rp_bins : numpy.ndarray + Boundaries defining the radial bins in which pairs are counted. + pi_max : float + Maximum line-of-sight distance. + sample2 : numpy.ndarray, optional + Numpy array containing the positions of a second sample of points. + Default is None. + period : float or numpy.ndarray, optional + If a numpy array the periodic boundary conditions in each + dimension. If a single scalar, period is assumed to be the same in all + dimensions. Default is None. + do_auto : bool, optional + Whether the calculate the auto-correlation function for `sample1`. If + True, returns the auto-correlation function for `sample1`, regardless + of whether `sample2` is given. If False, will return the + cross-correlation between `sample1` and `sample2`. Default is True. + do_cross : bool, optional + Whether the calculate the auto-correlation function for `sample1` and + `sample2`. Default is False. + + Returns + ------- + wp : numpy.ndarray + The two-point correlation function. + + Raises + ------ + RuntimeError + If Corrfunc is not installed. + ValueError + If `do_auto` and `do_cross` have the same value. + + """ + try: + from Corrfunc.theory import DDrppi + except ImportError: + raise RuntimeError('Corrfunc needs to be installed to use wrappers.') + + if (do_auto and do_cross) or (not do_auto and not do_cross): + raise ValueError("'do_auto' and 'do_cross' cannot both be True or " + + "False.") + + if isinstance(period, float) or isinstance(period, int): + period = (period, period, period) + + if isinstance(period, np.ndarray): + period = tuple(period) + + if do_auto: + r = DDrppi(1, 1, pi_max, rp_bins, sample1[:, 0], sample1[:, 1], + sample1[:, 2], periodic=True, boxsize=period, + xbin_refine_factor=1, ybin_refine_factor=1, + zbin_refine_factor=1) + n_exp = (len(sample1) * len(sample1) / np.prod(period) * np.pi * + np.diff(rp_bins**2) * 2 * pi_max) + + else: + r = DDrppi(0, 1, pi_max, rp_bins, sample1[:, 0], sample1[:, 1], + sample1[:, 2], periodic=True, boxsize=period, + X2=sample2[:, 0], Y2=sample2[:, 1], Z2=sample2[:, 2], + xbin_refine_factor=1, ybin_refine_factor=1, + zbin_refine_factor=1) + n_exp = (len(sample1) * len(sample2) / np.prod(period) * np.pi * + np.diff(rp_bins**2) * 2 * pi_max) + + npairs = r['npairs'] + npairs = np.array([np.sum(n) for n in np.split(npairs, len(rp_bins) - 1)]) + + return (npairs / n_exp - 1) * 2 * pi_max + + +def s_mu_tpcf(sample1, s_bins, mu_bins, sample2=None, period=None, + do_auto=True, do_cross=False): + """Corrfunc implementation of `halotools.mock_observables.s_mu_tpcf`. + + This function provides a wrapper around `corrfunc.theory.DDsmu` such that + when used by TabCorr behaves the same as + `halotools.mock_observables.s_mu_tpcf`. However, not all functionality of + `halotools.mock_observables.s_mu_tpcf` is implemented. + + Parameters + ---------- + sample1 : numpy.ndarray + Numpy array containing the positions of points. + s_bins : numpy.ndarray + Boundaries defining the distance bins in which pairs are counted. + mu_bins : numpy.ndarray + Boundaries defining the angular bins in which pairs are counted. + sample2 : numpy.ndarray, optional + Numpy array containing the positions of a second sample of points. + Default is None. + period : float or numpy.ndarray, optional + If a numpy array the periodic boundary conditions in each + dimension. If a single scalar, period is assumed to be the same in all + dimensions. Default is None. + do_auto : bool, optional + Whether the calculate the auto-correlation function for `sample1`. If + True, returns the auto-correlation function for `sample1`, regardless + of whether `sample2` is given. If False, will return the + cross-correlation between `sample1` and `sample2`. Default is True. + do_cross : bool, optional + Whether the calculate the auto-correlation function for `sample1` and + `sample2`. Default is False. + + Returns + ------- + xi : numpy.ndarray + The two-point correlation function. + + Raises + ------ + RuntimeError + If Corrfunc is not installed. + ValueError + If `do_auto` and `do_cross` have the same value or if `mu_bins` are not + uniform bins from 0 to 1. + + """ + try: + from Corrfunc.theory import DDsmu + except ImportError: + raise RuntimeError('Corrfunc needs to be installed to use wrappers.') + + if (do_auto and do_cross) or (not do_auto and not do_cross): + raise ValueError("'do_auto' and 'do_cross' cannot both be True or " + + "False.") + + try: + assert np.all(np.isclose(mu_bins, np.linspace(0, 1, len(mu_bins)))) + except AssertionError: + raise ValueError('Bins in mu must be uniform from 0 to 1.') + + if isinstance(period, float) or isinstance(period, int): + period = (period, period, period) + + if isinstance(period, np.ndarray): + period = tuple(period) + + if do_auto: + r = DDsmu(1, 1, s_bins, 1, len(mu_bins) - 1, sample1[:, 0], + sample1[:, 1], sample1[:, 2], periodic=True, + boxsize=period, xbin_refine_factor=1, + ybin_refine_factor=1, zbin_refine_factor=1) + n_exp = (len(sample1) * len(sample1) / np.prod(period) * 4 * + np.pi / 3 * np.diff(s_bins**3) / (len(mu_bins) - 1)) + + else: + r = DDsmu(0, 1, s_bins, 1, len(mu_bins) - 1, sample1[:, 0], + sample1[:, 1], sample1[:, 2], periodic=True, + boxsize=period, X2=sample2[:, 0], + Y2=sample2[:, 1], Z2=sample2[:, 2], xbin_refine_factor=1, + ybin_refine_factor=1, zbin_refine_factor=1) + n_exp = (len(sample1) * len(sample2) / np.prod(period) * 4 * + np.pi / 3 * np.diff(s_bins**3) / (len(mu_bins) - 1)) + + return (r['npairs'].reshape((len(s_bins) - 1, len(mu_bins) - 1)) / + n_exp[:, np.newaxis] - 1) diff --git a/tabcorr/database.py b/tabcorr/database.py new file mode 100644 index 0000000..730e4d5 --- /dev/null +++ b/tabcorr/database.py @@ -0,0 +1,280 @@ +"""Module providing database capabilities.""" + +import os +import numpy as np +from astropy import units as u +from astropy.table import Table +from . import TabCorr, Interpolator +from astropy.cosmology import Flatw0waCDM, FlatwCDM, Planck15 + + +def configuration(config_str): + """Describe the tabulation configuration used. + + Parameters + ---------- + config_str : str + String describing the configuration. You can specify a mixture of + different configuations by separating them with a `_`. In this case, + the first configurations in the list always take precedence in + case multiple configurations apply to a parameter. + + Returns + ------- + config_dict : dict + Dictionary descibing the configuration, i.e. radial binning, number + of satellites per halo etc. + + Raises + ------ + ValueError + If an unkown configuration was requested. + + """ + config_list = config_str.split('_') + + for config in config_list: + if config not in ['aemulus', 'default']: + raise ValueError('Unkown configuration {}.'.format(config)) + + config_list.append('default') + + config_dict = {} + config_dict['s_bins'] = {'default': np.logspace(-1.0, 1.8, 15), + 'aemulus': np.logspace(-1, 1.78, 10)} + config_dict['rp_wp_bins'] = {'default': np.logspace(-1.0, 1.8, 15), + 'aemulus': np.logspace(-1, 1.78, 10)} + config_dict['pi_max'] = {'default': 80} + config_dict['rp_ds_bins'] = {'default': np.logspace(-1.0, 1.8, 15)} + config_dict['mu_bins'] = {'default': np.linspace(0, 1, 21), + 'aemulus': np.linspace(0, 1, 41)} + config_dict['cosmo_obs'] = {'default': Planck15, 'aemulus': None} + config_dict['alpha_c_bins'] = {'default': np.linspace(0.0, 0.4, 4)} + config_dict['alpha_s_bins'] = {'default': np.linspace(0.8, 1.2, 4)} + config_dict['conc_gal_bias_bins'] = { + 'default': np.geomspace(1.0 / 3.0, 3.0, 4)} + config_dict['sats_per_prim_haloprop'] = {'default': 2e-13} + + for parameter in config_dict.keys(): + for config in config_list: + if config in config_dict[parameter].keys(): + config_dict[parameter] = config_dict[parameter][config] + break + + return config_dict + + +def directory(): + """Return the TabCorr database directory. + + Returns + ------- + dir : str + The TabCorr database directory. + + Raises + ------ + RuntimeError + If the TABCORR_DATABASE environment variable is not set. + + """ + try: + return os.environ['TABCORR_DATABASE'] + except KeyError: + raise RuntimeError( + "You must set the TABCORR_DATABASE environment variable.") + + +def cosmology(suite, i_cosmo=0): + """Return the cosmology of a given simulation. + + Parameters + ---------- + suite : str + The simulation suite. + i_cosmo : int, optional + Number corresponding to the cosmology. Default is 0. + + Returns + ------- + cosmo : str + The TabCorr database directory. + + Raises + ------ + ValueError + If the an unkown simulation or cosmology number is requested. + + """ + if suite == 'AbacusSummit': + table = Table.read(os.path.join( + os.path.dirname(os.path.realpath(__file__)), 'as_cosmos.csv')) + table['i_cosmo'] = np.array([r[-3:] for r in table['root']], dtype=int) + if i_cosmo not in table['i_cosmo']: + raise ValueError('Cosmology number {} not in AbacusSummit.'.format( + i_cosmo)) + cosmo_dict = dict(table[table['i_cosmo'] == i_cosmo][0]) + h = cosmo_dict['h'] + omega_m = cosmo_dict['omega_b'] + cosmo_dict['omega_cdm'] + n_eff = cosmo_dict['N_ur'] + cosmo_dict['N_ncdm'] + m_nu = [float(omega) * 93.04 * u.eV for omega in + cosmo_dict['omega_ncdm'].split(',')] + assert len(m_nu) == cosmo_dict['N_ncdm'] + while len(m_nu) < n_eff - 1: + m_nu.append(0 * u.eV) + return Flatw0waCDM( + H0=h * 100, Om0=omega_m / h**2, Ob0=cosmo_dict['omega_b'] / h**2, + w0=cosmo_dict['w0_fld'], wa=cosmo_dict['wa_fld'], Neff=n_eff, + m_nu=m_nu, Tcmb0=2.7255 * u.K) + + elif suite == 'AemulusAlpha': + path = os.path.dirname(os.path.realpath(__file__)) + if i_cosmo >= 0 and i_cosmo < 40: + cosmo_dict = dict(Table.read( + os.path.join(path, 'aa_cosmos.txt'), format='ascii')[i_cosmo]) + elif i_cosmo >= 0 and i_cosmo < 47: + cosmo_dict = dict(Table.read( + os.path.join(path, 'aa_test_cosmos.txt'), format='ascii')[ + i_cosmo - 40]) + else: + raise ValueError('Unknown cosmology number {}. '.format(i_cosmo) + + 'Must be in the range from 0 to 46.') + cosmo_dict['Ob0'] = cosmo_dict['ombh2'] / (cosmo_dict['H0'] / 100)**2 + cosmo_dict['Oc0'] = cosmo_dict['omch2'] / (cosmo_dict['H0'] / 100)**2 + cosmo_dict['Om0'] = cosmo_dict['Ob0'] + cosmo_dict['Oc0'] + return FlatwCDM(H0=cosmo_dict['H0'], Om0=cosmo_dict['Om0'], + w0=cosmo_dict['w0'], Neff=cosmo_dict['Neff'], + Ob0=cosmo_dict['Ob0'], Tcmb0=2.7255 * u.K) + + else: + raise ValueError('Unkown simulation suite {}.'.format(suite)) + + +def simulation_name(suite, i_cosmo=0, i_phase=0, config=None): + """Return the name of a given simulation. + + Parameters + ---------- + suite : str + The simulation suite. + i_cosmo : int, optional + If applicable, number corresponding to the cosmology. Default is 0. + i_phase : int, optional + If applicable, number corresponding to the simulation phase. Default + is 0. + config : str + Simulation configuration. Only applicable to AbacusSummit. If None, + will default to 'base' for AbacusSummit. Default is None. + + Returns + ------- + name : str + The TabCorr database directory. + + Raises + ------ + ValueError + If the an unkown simulation, simulation number or phase number is + requested. + + """ + if suite == 'AbacusSummit': + + if config is None: + config = 'base' + + return '{}_c{:03d}_ph{:03d}'.format(config, i_cosmo, i_phase) + + elif suite == 'AemulusAlpha': + + if i_cosmo >= 0 and i_cosmo < 40: + return 'Box{:03d}'.format(i_cosmo) + elif i_cosmo >= 0 and i_cosmo < 47: + if i_phase > 6: + raise ValueError( + 'Unknown phase number {}.'.format(i_phase)) + return 'TestBox{:03d}-{:03d}'.format(i_cosmo - 40, i_phase) + else: + raise ValueError('Unknown cosmology number {}. '.format(i_cosmo) + + 'Must be in the range from 0 to 46.') + + else: + raise ValueError('Unkown simulation suite {}.'.format(suite)) + + +def simulation_snapshot_directory( + suite, redshift, i_cosmo=0, i_phase=0, config=None): + """Return the directory where all data for a simulation snapshot is stored. + + Parameters + ---------- + suite : str + The simulation suite. + redshift : float + The redshift of the simulation output. + i_cosmo : int, optional + If applicable, number corresponding to the cosmology. Default is 0. + i_phase : int, optional + If applicable, number corresponding to the simulation phase. Default + is 0. + config : str + Simulation configuration. Only applicable to AbacusSummit. If None, + will default to 'base' for AbacusSummit. Default is None. + + Returns + ------- + name : str + The directory where all data for a simulation snapshot is stored. + + """ + name = simulation_name(suite, i_cosmo=i_cosmo, i_phase=i_phase, + config=config) + return os.path.join( + directory(), suite, name, + '{:.2f}'.format(redshift).replace('.', 'p')) + + +def tabcorr(suite, redshift, tpcf, i_cosmo=0, i_phase=0, sim_config=None, + tab_config='default'): + """Return the TabCorr tabulation for a given simulation, redshift etc. + + Parameters + ---------- + suite : str + The simulation suite. + redshift : float + The redshift of the simulation output. + tpcf : str + String describing the two-point correlation function. + i_cosmo : int, optional + If applicable, number corresponding to the cosmology. Default is 0. + i_phase : int, optional + If applicable, number corresponding to the simulation phase. Default + is 0. + sim_config : str + Simulation configuration. Only applicable to AbacusSummit. If None, + will default to 'base' for AbacusSummit. Default is None. + tab_config : config_str : str + String describing the configuration of the tabulation, i.e. binning, + cosmology etc. + + Returns + ------- + halotab : tabcorr.TabCorr or tabcorr.Interpolator + The tabcorr tabulation. + + """ + directory = os.path.join(simulation_snapshot_directory( + suite, redshift, i_cosmo=i_cosmo, i_phase=i_phase, config=sim_config), + tab_config) + + param_dict_table = Table.read(os.path.join( + directory, tpcf[:2] + '_grid.csv')) + param_dict_table['log_eta'] = np.log10(param_dict_table['conc_gal_bias']) + param_dict_table.remove_column('conc_gal_bias') + for key in ['alpha_c', 'alpha_s', 'log_eta']: + if len(np.unique(param_dict_table[key])) == 1: + param_dict_table.remove_column(key) + tabcorr_list = [TabCorr.read(os.path.join(directory, '{}_{}.hdf5'.format( + tpcf, i))) for i in range(len(param_dict_table))] + return Interpolator(tabcorr_list, param_dict_table) diff --git a/tabcorr/interpolate.py b/tabcorr/interpolate.py index b00782d..f7bafec 100644 --- a/tabcorr/interpolate.py +++ b/tabcorr/interpolate.py @@ -78,8 +78,16 @@ def __init__(self, tabcorr_list, param_dict_table, spline=True): self.xp[:, i] = self.param_dict_table[key].data self.delaunay = Delaunay(self.xp) - def predict(self, model, n_gauss_prim=10, extrapolate=False, - same_halos=False, **occ_kwargs): + # Determine unique halo tables such that we can save computation time + # if halo tables are repeated. + all_gal_type = [np.array(tabcorr.gal_type.as_array().tolist()).ravel() + for tabcorr in tabcorr_list] + unique = np.unique( + all_gal_type, axis=0, return_index=True, return_inverse=True) + self.unique_gal_type_index = unique[1] + self.unique_gal_type_inverse = unique[2] + + def predict(self, model, n_gauss_prim=10, extrapolate=False, **occ_kwargs): """Interpolate the predictions from multiple TabCorr instances. The values of parameters to interpolate should be in the parameter @@ -97,11 +105,6 @@ def predict(self, model, n_gauss_prim=10, extrapolate=False, extrapolate : bool, optional Whether to allow extrapolation beyond points sampled by the input TabCorr instances. Default is False. - same_halos : bool, optional - If True, assume that all TabCorr instances are using the same - simulation and halo bins. This is typically the case if we want - to interpolate between phase-space parameters. Can speed up - calculation considerably. Default is False. **occ_kwargs : dict, optional Keyword arguments passed to the ``mean_occupation`` functions of the model. @@ -130,14 +133,19 @@ def predict(self, model, n_gauss_prim=10, extrapolate=False, 'dictionary of the model.') if self.spline: + + # Calculate the mean occupation numbers, avoiding to calculate + # those repeatedly for identical halo tables. + mean_occupation = [self.tabcorr_list[i].mean_occupation( + model, n_gauss_prim=n_gauss_prim, **occ_kwargs) for i in + self.unique_gal_type_index] + for i in range(len(self.param_dict_table)): - tabcorr = self.tabcorr_list[ - self.param_dict_table['tabcorr_index'][i]] - if i == 0 and same_halos: - model = tabcorr.mean_occupation( - model, n_gauss_prim=n_gauss_prim, **occ_kwargs) + k = self.param_dict_table['tabcorr_index'][i] + tabcorr = self.tabcorr_list[k] ngal_i, xi_i = tabcorr.predict( - model, n_gauss_prim=n_gauss_prim, **occ_kwargs) + mean_occupation[self.unique_gal_type_inverse[k]], + n_gauss_prim=n_gauss_prim, **occ_kwargs) if i == 0: ngal = np.zeros(np.prod([len(xp) for xp in self.xp])) xi = np.zeros([np.prod([len(xp) for xp in self.xp])] + @@ -189,6 +197,7 @@ def spline_interpolation_matrix(xp): Returns ------- + a : numpy.ndarray Matrix `a` such that ``np.einsum('ij,j,i', a[i], y, x0**np.arange(4))`` is the value `i`-th spline evalualated at `x0`. @@ -261,7 +270,8 @@ def spline_interpolate(x, xp, a, yp, extrapolate=False): Returns ------- - Interpolated values with shape ``y.shape[len(x):]``. + yp : float or numpy.ndarray + Interpolated value(s) with shape ``y.shape[len(x):]``. Raises ------ diff --git a/tabcorr/tabcorr.py b/tabcorr/tabcorr.py index 83a8c03..e70077b 100644 --- a/tabcorr/tabcorr.py +++ b/tabcorr/tabcorr.py @@ -124,6 +124,7 @@ def tabulate(cls, halocat, tpcf, *tpcf_args, Returns ------- + halotab : tabcorr.TabCorr `TabCorr` object. Raises @@ -474,6 +475,7 @@ def mean_occupation(self, model, n_gauss_prim=10, **occ_kwargs): Returns ------- + n : numpy.ndarray Array containing the mean occuaption numbers. Has the same length as `self.gal_type`. @@ -686,6 +688,7 @@ def sort_into_bins(log_prim_haloprop, log_prim_haloprop_bins, Returns ------- + x_sorted : list Values of `x` sorted into bins. """ @@ -733,6 +736,7 @@ def distribution_index(x_min, x_max, x_mean): Returns ------- + n : float Index :math:`n` to reproduce the mean, but not smaller than -10 or larger than +10. """ @@ -758,6 +762,7 @@ def symmetric_matrix_to_array(matrix, check_symmetry=True): Returns ------- + m_array : numpy.ndarray Array containing all unique values of `matrix`. Raises