diff --git a/cluster_toolkit/__init__.py b/cluster_toolkit/__init__.py index 0f8ab93..04d6f8b 100644 --- a/cluster_toolkit/__init__.py +++ b/cluster_toolkit/__init__.py @@ -26,12 +26,62 @@ lib_file = alt_files[0] _ffi = cffi.FFI() -for file_name in glob.glob(os.path.join(include_dir,'*.h')): +for file_name in glob.glob(os.path.join(include_dir, '*.h')): _ffi.cdef(open(file_name).read()) +_ffi.cdef('const char * gsl_strerror(const int gsl_errno);') + _lib = _ffi.dlopen(lib_file) +_lib.gsl_set_error_handler_off() + + +def _handle_gsl_error(err, fn): + if err != 0: + msg = _ffi.string(_lib.gsl_strerror(err)) + raise Exception('GSL error in function {}: {}'.format(fn.__name__, msg)) + + +class _ArrayWrapper: + def __init__(self, obj, name=None, allow_multidim=False): + self.arr = np.require(obj, dtype=np.float64, + requirements=['C_CONTIGUOUS']) + self.scalar = self.arr.ndim == 0 + self.ndim = self.arr.ndim + self.shape = self.arr.shape + + if (self.ndim > 1) and not allow_multidim: + if name is not None: + raise ValueError('{} cannot be >1 dim'.format(name)) + raise ValueError('array cannot be >1 dim') + + def cast(self): + return _ffi.cast('double*', self.arr.ctypes.data) + + def finish(self): + if self.scalar: + return self.arr[()] + return self.arr + + def __len__(self): + if self.scalar: + return 1 + return self.arr.size + + @classmethod + def zeros_like(cls, obj): + if isinstance(obj, _ArrayWrapper): + return cls(np.zeros_like(obj.arr)) + return cls(np.zeros_like(obj)) + + @classmethod + def zeros(cls, shape): + return cls(np.zeros(shape, dtype=np.double)) + + @classmethod + def ones_like(cls, obj): + return cls(np.ones_like(obj)) -def _dcast(x): - if isinstance(x, list): x = np.asarray(x, dtype=np.float64, order='C') - return _ffi.cast('double*', x.ctypes.data) + @classmethod + def ones(cls, shape): + return cls(np.ones(shape, dtype=np.double)) from . import averaging, bias, boostfactors, concentration, deltasigma, density, exclusion, massfunction, miscentering, peak_height, profile_derivatives, sigma_reconstruction, xi diff --git a/cluster_toolkit/averaging.py b/cluster_toolkit/averaging.py index 45a3807..0f14e03 100644 --- a/cluster_toolkit/averaging.py +++ b/cluster_toolkit/averaging.py @@ -2,13 +2,14 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper, _handle_gsl_error import numpy as np + def average_profile_in_bins(Redges, R, prof): """Average profile in bins. - Calculates the average of some projected profile in a + Calculates the average of some projected profile in a radial bins in Mpc/h comoving. Args: @@ -20,27 +21,34 @@ def average_profile_in_bins(Redges, R, prof): numpy.array: Average profile in bins between the edges provided. """ - Redges = np.asarray(Redges) - if Redges.ndim == 0: #Is scalar - Redges = r[None] #makes r 1D + Redges = _ArrayWrapper(Redges) + R = _ArrayWrapper(R) + prof = _ArrayWrapper(prof) + + if Redges.ndim == 0: raise Exception("Must supply a left and right edge.") if Redges.ndim > 1: raise Exception("Redges cannot be a >1D array.") - if np.min(Redges) < np.min(R): + if np.min(Redges.arr) < np.min(R.arr): raise Exception("Minimum edge must be >= minimum R") - if np.max(Redges) > np.max(R): + if np.max(Redges.arr) > np.max(R.arr): raise Exception("Maximum edge must be <= maximum R") - - ave_prof = np.zeros(len(Redges)-1) - cluster_toolkit._lib.average_profile_in_bins(_dcast(Redges), len(Redges), - _dcast(R), len(R), _dcast(prof), - _dcast(ave_prof)) - return ave_prof + + ave_prof = _ArrayWrapper(np.zeros(len(Redges) - 1)) + r = cluster_toolkit._lib.average_profile_in_bins(Redges.cast(), len(Redges), + R.cast(), len(R), + prof.cast(), + ave_prof.cast()) + + _handle_gsl_error(r, average_profile_in_bins) + + return ave_prof.finish() + def average_profile_in_bin(Rlow, Rhigh, R, prof): """Average profile in a bin. - Calculates the average of some projected profile in a + Calculates the average of some projected profile in a radial bin in Mpc/h comoving. Args: diff --git a/cluster_toolkit/bias.py b/cluster_toolkit/bias.py index d71d380..59deb63 100644 --- a/cluster_toolkit/bias.py +++ b/cluster_toolkit/bias.py @@ -2,9 +2,9 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper import numpy as np -from .peak_height import * +# from .peak_height import * def bias_at_M(M, k, P, Omega_m, delta=200): """Tinker et al. 2010 bais at mass M [Msun/h]. @@ -20,21 +20,17 @@ def bias_at_M(M, k, P, Omega_m, delta=200): float or array like: Halo bias. """ - M = np.asarray(M) - scalar_input = False - if M.ndim == 0: - M = M[None] #makes r 1D - scalar_input = True - if M.ndim > 1: - raise Exception("M cannot be a >1D array.") - M = np.require(M, dtype=np.float64, requirements=["C"]) - k = np.require(k, dtype=np.float64, requirements=["C"]) - P = np.require(P, dtype=np.float64, requirements=["C"]) - bias = np.zeros_like(M) - cluster_toolkit._lib.bias_at_M_arr(_dcast(M), len(M), delta, _dcast(k), _dcast(P), len(k), Omega_m, _dcast(bias)) - if scalar_input: - return np.squeeze(bias) - return bias + M = _ArrayWrapper(M, 'M') + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + if k.shape != P.shape: + raise ValueError('k and P must have the same shape') + + bias = _ArrayWrapper.zeros_like(M) + cluster_toolkit._lib.bias_at_M_arr(M.cast(), len(M), delta, + k.cast(), P.cast(), len(k), + Omega_m, bias.cast()) + return bias.finish() def bias_at_R(R, k, P, delta=200): """Tinker 2010 bais at mass M [Msun/h] corresponding to radius R [Mpc/h comoving]. @@ -49,22 +45,16 @@ def bias_at_R(R, k, P, delta=200): float or array like: Halo bias. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes r 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - R = np.require(R, dtype=np.float64, requirements=["C"]) - k = np.require(k, dtype=np.float64, requirements=["C"]) - P = np.require(P, dtype=np.float64, requirements=["C"]) - bias = np.zeros_like(R) - cluster_toolkit._lib.bias_at_R_arr(_dcast(R), len(R), delta, _dcast(k), _dcast(P), len(k), _dcast(bias)) - if scalar_input: - return np.squeeze(bias) - return bias - + R = _ArrayWrapper(R, 'R') + k = _ArrayWrapper(k) + P = _ArrayWrapper(P) + + bias = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.bias_at_R_arr(R.cast(), len(R), delta, + k.cast(), P.cast(), len(k), + bias.cast()) + return bias.finish() + def bias_at_nu(nu, delta=200): """Tinker 2010 bais at peak height nu. @@ -76,19 +66,12 @@ def bias_at_nu(nu, delta=200): float or array like: Halo bias. """ - nu = np.asarray(nu) - scalar_input = False - if nu.ndim == 0: - nu = nu[None] #makes nu 1D - scalar_input = True - if nu.ndim > 1: - raise Exception("nu cannot be a >1D array.") - - bias = np.zeros_like(nu) - cluster_toolkit._lib.bias_at_nu_arr(_dcast(nu), len(nu), delta, _dcast(bias)) - if scalar_input: - return np.squeeze(bias) - return bias + nu = _ArrayWrapper(nu, 'nu') + + bias = _ArrayWrapper.zeros_like(nu) + cluster_toolkit._lib.bias_at_nu_arr(nu.cast(), len(nu), delta, + bias.cast()) + return bias.finish() def dbiasdM_at_M(M, k, P, Omega_m, delta=200): """d/dM of Tinker et al. 2010 bais at mass M [Msun/h]. @@ -102,31 +85,25 @@ def dbiasdM_at_M(M, k, P, Omega_m, delta=200): Returns: float or array like: Derivative of the halo bias. - + """ - M = np.asarray(M) - scalar_input = False - if M.ndim == 0: - M = M[None] #makes M 1D - scalar_input = True - if M.ndim > 1: - raise Exception("M cannot be a >1D array.") - M = np.require(M, dtype=np.float64, requirements=["C"]) - k = np.require(k, dtype=np.float64, requirements=["C"]) - P = np.require(P, dtype=np.float64, requirements=["C"]) - deriv = np.zeros_like(M) - cluster_toolkit._lib.dbiasdM_at_M_arr(_dcast(M), len(M), delta, _dcast(k), - _dcast(P), len(k), Omega_m, - _dcast(deriv)) - if scalar_input: - return np.squeeze(deriv) - return deriv + M = _ArrayWrapper(M, 'M') + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + + deriv = _ArrayWrapper.zeros_like(M) + cluster_toolkit._lib.dbiasdM_at_M_arr(M.cast(), len(M), delta, k.cast(), + P.cast(), len(k), Omega_m, + deriv.cast()) + return deriv.finish() def _bias_at_nu_FREEPARAMS(nu, A, a, B, b, C, c, delta=200): """A special function used only for quickly computing best fit parameters for the halo bias models. """ - bias = np.zeros_like(nu) - cluster_toolkit._lib.bias_at_nu_arr_FREEPARAMS(_dcast(nu), len(nu), delta, - A, a, B, b, C, c, _dcast(bias)) - return bias + nu = _ArrayWrapper(nu, allow_multidim=True) + bias = _ArrayWrapper.zeros_like(nu) + cluster_toolkit._lib.bias_at_nu_arr_FREEPARAMS(nu.cast(), len(nu), delta, + A, a, B, b, C, c, + bias.cast()) + return bias.finish() diff --git a/cluster_toolkit/boostfactors.py b/cluster_toolkit/boostfactors.py index d74add7..d111aba 100644 --- a/cluster_toolkit/boostfactors.py +++ b/cluster_toolkit/boostfactors.py @@ -2,7 +2,7 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper import numpy as np def boost_nfw_at_R(R, B0, R_scale): @@ -17,20 +17,12 @@ def boost_nfw_at_R(R, B0, R_scale): float or array like: NFW boost factor profile; B = (1-fcl)^-1. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - - boost = np.zeros_like(R) - cluster_toolkit._lib.boost_nfw_at_R_arr(_dcast(R), len(R), B0, R_scale, - _dcast(boost)) - if scalar_input: - return np.squeeze(boost) - return boost + R = _ArrayWrapper(R, 'R') + + boost = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.boost_nfw_at_R_arr(R.cast(), len(R), B0, R_scale, + boost.cast()) + return boost.finish() def boost_powerlaw_at_R(R, B0, R_scale, alpha): """Power law boost factor model. @@ -45,17 +37,9 @@ def boost_powerlaw_at_R(R, B0, R_scale, alpha): float or array like: Power law boost factor profile; B = (1-fcl)^-1. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - - boost = np.zeros_like(R) - cluster_toolkit._lib.boost_powerlaw_at_R_arr(_dcast(R), len(R), B0, - R_scale, alpha, _dcast(boost)) - if scalar_input: - return np.squeeze(boost) - return boost + R = _ArrayWrapper(R, 'R') + + boost = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.boost_powerlaw_at_R_arr(R.cast(), len(R), B0, + R_scale, alpha, boost.cast()) + return boost.finish() diff --git a/cluster_toolkit/concentration.py b/cluster_toolkit/concentration.py index 570b603..c60567a 100644 --- a/cluster_toolkit/concentration.py +++ b/cluster_toolkit/concentration.py @@ -2,7 +2,7 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper import numpy as np def concentration_at_M(Mass, k, P, n_s, Omega_b, Omega_m, h, T_CMB=2.7255, delta=200, Mass_type="crit"): @@ -10,7 +10,7 @@ def concentration_at_M(Mass, k, P, n_s, Omega_b, Omega_m, h, T_CMB=2.7255, delta Only implemented relation at the moment is Diemer & Kravtsov (2015). Note: only single concentrations at a time are allowed at the moment. - + Args: Mass (float): Mass in Msun/h. k (array like): Wavenumbers of power spectrum in h/Mpc comoving. @@ -29,11 +29,13 @@ def concentration_at_M(Mass, k, P, n_s, Omega_b, Omega_m, h, T_CMB=2.7255, delta """ if delta != 200: raise Exception("ConcentrationError: delta=%d. Currently only delta=200 supported"%delta) + + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + if Mass_type is "mean": - return cluster_toolkit._lib.DK15_concentration_at_Mmean(Mass, _dcast(k), _dcast(P), len(k), delta, n_s, Omega_b, Omega_m, h, T_CMB) + return cluster_toolkit._lib.DK15_concentration_at_Mmean(Mass, k.cast(), P.cast(), len(k), delta, n_s, Omega_b, Omega_m, h, T_CMB) elif Mass_type is "crit": - return cluster_toolkit._lib.DK15_concentration_at_Mcrit(Mass, _dcast(k), _dcast(P), len(k), delta, n_s, Omega_b, Omega_m, h, T_CMB) + return cluster_toolkit._lib.DK15_concentration_at_Mcrit(Mass, k.cast(), P.cast(), len(k), delta, n_s, Omega_b, Omega_m, h, T_CMB) else: raise Exception("ConcentrationError: must choose either 'mean' or 'crit', %s is not supported"%Mass_type) - - diff --git a/cluster_toolkit/deltasigma.py b/cluster_toolkit/deltasigma.py index 44befce..0c45044 100644 --- a/cluster_toolkit/deltasigma.py +++ b/cluster_toolkit/deltasigma.py @@ -2,7 +2,7 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper, _handle_gsl_error import numpy as np def Sigma_nfw_at_R(R, mass, concentration, Omega_m, delta=200): @@ -19,21 +19,13 @@ def Sigma_nfw_at_R(R, mass, concentration, Omega_m, delta=200): float or array like: Surface mass density Msun h/pc^2 comoving. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - - Sigma = np.zeros_like(R) - cluster_toolkit._lib.Sigma_nfw_at_R_arr(_dcast(R), len(R), mass, + R = _ArrayWrapper(R, 'R') + + Sigma = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.Sigma_nfw_at_R_arr(R.cast(), len(R), mass, concentration, delta, - Omega_m, _dcast(Sigma)) - if scalar_input: - return np.squeeze(Sigma) - return Sigma + Omega_m, Sigma.cast()) + return Sigma.finish() def Sigma_at_R(R, Rxi, xi, mass, concentration, Omega_m, delta=200): """Surface mass density given some 3d profile [Msun h/pc^2 comoving]. @@ -51,25 +43,23 @@ def Sigma_at_R(R, Rxi, xi, mass, concentration, Omega_m, delta=200): float or array like: Surface mass density Msun h/pc^2 comoving. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - if np.min(R) < np.min(Rxi): + R = _ArrayWrapper(R, 'R') + Rxi = _ArrayWrapper(Rxi, allow_multidim=True) + xi = _ArrayWrapper(xi, allow_multidim=True) + + if np.min(R.arr) < np.min(Rxi.arr): raise Exception("Minimum R for Sigma(R) must be >= than min(r) of xi(r).") - if np.max(R) > np.max(Rxi): + if np.max(R.arr) > np.max(Rxi.arr): raise Exception("Maximum R for Sigma(R) must be <= than max(r) of xi(r).") - Sigma = np.zeros_like(R) - cluster_toolkit._lib.Sigma_at_R_full_arr(_dcast(R), len(R), _dcast(Rxi), - _dcast(xi), len(Rxi), mass, concentration, - delta, Omega_m, _dcast(Sigma)) - if scalar_input: - return np.squeeze(Sigma) - return Sigma + Sigma = _ArrayWrapper.zeros_like(R) + rc = cluster_toolkit._lib.Sigma_at_R_full_arr(R.cast(), len(R), Rxi.cast(), + xi.cast(), len(Rxi), mass, concentration, + delta, Omega_m, Sigma.cast()) + + _handle_gsl_error(rc, Sigma_at_R) + + return Sigma.finish() def DeltaSigma_at_R(R, Rs, Sigma, mass, concentration, Omega_m, delta=200): """Excess surface mass density given Sigma [Msun h/pc^2 comoving]. @@ -87,26 +77,23 @@ def DeltaSigma_at_R(R, Rs, Sigma, mass, concentration, Omega_m, delta=200): float or array like: Excess surface mass density Msun h/pc^2 comoving. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - - if np.min(R) < np.min(Rs): + R = _ArrayWrapper(R, 'R') + Rs = _ArrayWrapper(Rs, allow_multidim=True) + Sigma = _ArrayWrapper(Sigma, allow_multidim=True) + + if np.min(R.arr) < np.min(Rs.arr): raise Exception("Minimum R for DeltaSigma(R) must be "+ ">= than min(R) of Sigma(R).") - if np.max(R) > np.max(Rs): + if np.max(R.arr) > np.max(Rs.arr): raise Exception("Maximum R for DeltaSigma(R) must be "+ "<= than max(R) of Sigma(R).") - DeltaSigma = np.zeros_like(R) - cluster_toolkit._lib.DeltaSigma_at_R_arr(_dcast(R), len(R), _dcast(Rs), - _dcast(Sigma), len(Rs), mass, - concentration, delta, Omega_m, - _dcast(DeltaSigma)) - if scalar_input: - return np.squeeze(DeltaSigma) - return DeltaSigma + DeltaSigma = _ArrayWrapper.zeros_like(R) + rc = cluster_toolkit._lib.DeltaSigma_at_R_arr(R.cast(), len(R), Rs.cast(), + Sigma.cast(), len(Rs), mass, + concentration, delta, Omega_m, + DeltaSigma.cast()) + + _handle_gsl_error(rc, DeltaSigma_at_R) + + return DeltaSigma.finish() diff --git a/cluster_toolkit/density.py b/cluster_toolkit/density.py index 6540548..0f5e229 100644 --- a/cluster_toolkit/density.py +++ b/cluster_toolkit/density.py @@ -2,7 +2,7 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper import numpy as np def rho_nfw_at_r(r, M, c, Omega_m, delta=200): @@ -19,20 +19,12 @@ def rho_nfw_at_r(r, M, c, Omega_m, delta=200): float or array like: NFW halo density profile in Msun h^2/Mpc^3 comoving. """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") + r = _ArrayWrapper(r, 'r') - rho = np.zeros_like(r) - cluster_toolkit._lib.calc_rho_nfw(_dcast(r), len(r), M, c, delta, - Omega_m, _dcast(rho)) - if scalar_input: - return np.squeeze(rho) - return rho + rho = _ArrayWrapper.zeros_like(r) + cluster_toolkit._lib.calc_rho_nfw(r.cast(), len(r), M, c, delta, + Omega_m, rho.cast()) + return rho.finish() def rho_einasto_at_r(r, M, rs, alpha, Omega_m, delta=200, rhos=-1.): @@ -51,18 +43,10 @@ def rho_einasto_at_r(r, M, rs, alpha, Omega_m, delta=200, rhos=-1.): float or array like: Einasto halo density profile in Msun h^2/Mpc^3 comoving. """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") + r = _ArrayWrapper(r, 'r') - rho = np.zeros_like(r) - cluster_toolkit._lib.calc_rho_einasto(_dcast(r), len(r), M, rhos, rs, - alpha, delta, Omega_m, _dcast(rho)) - if scalar_input: - return np.squeeze(rho) - return rho + rho = _ArrayWrapper.zeros_like(r) + cluster_toolkit._lib.calc_rho_einasto(r.cast(), len(r), M, rhos, rs, + alpha, delta, Omega_m, rho.cast()) + return rho.finish() diff --git a/cluster_toolkit/exclusion.py b/cluster_toolkit/exclusion.py index 4cd21c4..a9af6af 100644 --- a/cluster_toolkit/exclusion.py +++ b/cluster_toolkit/exclusion.py @@ -1,7 +1,7 @@ """Correlation functions with halo exclusion. """ import cluster_toolkit -from cluster_toolkit import _dcast as dc +from cluster_toolkit import _ArrayWrapper, _handle_gsl_error import numpy as np def xi_hm_exclusion_at_r(radii, Mass, conc, alpha, @@ -102,7 +102,7 @@ def xi_2h_exclusion_at_r(radii, r_eff, beta_eff, bias, xi_mm): r_eff (float): effective radius for 2-halo subtraction in Mpc/h beta_eff (float): width for effective radius truncation bias (float): halo bias at large scales - xi_mm (float or array-like): matter correlation function. + xi_mm (float or array-like): matter correlation function. Must have same shape as the radii. Returns: @@ -165,7 +165,7 @@ def xi_C_at_r(radii, r_A, r_B, beta_ex, xi_2h): def theta_at_r(radii, rt, beta): """Truncation function. - + Args: radii (float or array-like): Radii of the profile in Mpc/h rt (float): truncation radius in Mpc/h @@ -184,8 +184,9 @@ def theta_at_r(radii, rt, beta): raise Exception("radii cannot be a >1D array.") theta = np.zeros_like(radii) - cluster_toolkit._lib.theta_erfc_at_r_arr(dc(radii), len(radii), - rt, beta, dc(theta)) + rc = cluster_toolkit._lib.theta_erfc_at_r_arr(dc(radii), len(radii), + rt, beta, dc(theta)) + _handle_gsl_error(rc) if scalar_input: return np.squeeze(theta) return theta diff --git a/cluster_toolkit/massfunction.py b/cluster_toolkit/massfunction.py index 05a9c86..aefb6f6 100644 --- a/cluster_toolkit/massfunction.py +++ b/cluster_toolkit/massfunction.py @@ -2,11 +2,11 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper, _handle_gsl_error import numpy as np def dndM_at_M(M, k, P, Omega_m, d=1.97, e=1.0, f=0.51, g=1.228): - """Tinker et al. 2008 appendix C mass function at a given mass. + """Tinker et al. 2008 appendix C mass function at a given mass. Default behavior is for :math:`M_{200m}` mass definition. NOTE: by default, this function is only valid at :math:`z=0`. For use @@ -27,28 +27,19 @@ def dndM_at_M(M, k, P, Omega_m, d=1.97, e=1.0, f=0.51, g=1.228): float or array like: Mass function :math:`dn/dM`. """ - M = np.asarray(M) - scalar_input = False - if M.ndim == 0: - M = M[None] #makes M 1D - scalar_input = True - if M.ndim > 1: - raise Exception("M cannot be a >1D array.") - M = np.require(M, dtype=np.float64, requirements=["C"]) - k = np.require(k, dtype=np.float64, requirements=["C"]) - P = np.require(P, dtype=np.float64, requirements=["C"]) - - dndM = np.zeros_like(M) - cluster_toolkit._lib.dndM_at_M_arr(_dcast(M), len(M), _dcast(k), - _dcast(P), len(k), Omega_m, - d, e, f, g, _dcast(dndM)) - if scalar_input: - return np.squeeze(dndM) - return dndM + M = _ArrayWrapper(M, 'M') + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + + dndM = _ArrayWrapper.zeros_like(M) + cluster_toolkit._lib.dndM_at_M_arr(M.cast(), len(M), k.cast(), + P.cast(), len(k), Omega_m, + d, e, f, g, dndM.cast()) + return dndM.finish() def G_at_M(M, k, P, Omega_m, d=1.97, e=1.0, f=0.51, g=1.228): - """Tinker et al. 2008 appendix C multiplicity funciton G(M) as - a function of mass. Default behavior is for :math:`M_{200m}` mass + """Tinker et al. 2008 appendix C multiplicity funciton G(M) as + a function of mass. Default behavior is for :math:`M_{200m}` mass definition. Args: @@ -64,25 +55,18 @@ def G_at_M(M, k, P, Omega_m, d=1.97, e=1.0, f=0.51, g=1.228): Returns: float or array like: Halo multiplicity :math:`G(M)`. """ - M = np.asarray(M) - scalar_input = False - if M.ndim == 0: - M = M[None] #makes M 1D - scalar_input = True - if M.ndim > 1: - raise Exception("M cannot be a >1D array.") - M = np.require(M, dtype=np.float64, requirements=["C"]) - k = np.require(k, dtype=np.float64, requirements=["C"]) - P = np.require(P, dtype=np.float64, requirements=["C"]) - - G = np.zeros_like(M) - cluster_toolkit._lib.G_at_M_arr(_dcast(M), len(M), _dcast(k), _dcast(P), len(k), Omega_m, d, e, f, g, _dcast(G)) - if scalar_input: - return np.squeeze(G) - return G + M = _ArrayWrapper(M, 'M') + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + + G = _ArrayWrapper.zeros_like(M) + cluster_toolkit._lib.G_at_M_arr(M.cast(), len(M), + k.cast(), P.cast(), len(k), + Omega_m, d, e, f, g, G.cast()) + return G.finish() def G_at_sigma(sigma, d=1.97, e=1.0, f=0.51, g=1.228): - """Tinker et al. 2008 appendix C multiplicity funciton G(sigma) as + """Tinker et al. 2008 appendix C multiplicity funciton G(sigma) as a function of sigma. NOTE: by default, this function is only valid at :math:`z=0`. For use @@ -99,19 +83,12 @@ def G_at_sigma(sigma, d=1.97, e=1.0, f=0.51, g=1.228): Returns: float or array like: Halo multiplicity G(sigma). """ - sigma = np.asarray(sigma) - scalar_input = False - if sigma.ndim == 0: - sigma = sigma[None] - scalar_input = True - if sigma.ndim > 1: - raise Exception("sigma cannot be a >1D array.") - - G = np.zeros_like(sigma) - cluster_toolkit._lib.G_at_sigma_arr(_dcast(sigma), len(sigma), d, e, f, g, _dcast(G)) - if scalar_input: - return np.squeeze(G) - return G + sigma = _ArrayWrapper(sigma, 'sigma') + + G = _ArrayWrapper.zeros_like(sigma) + cluster_toolkit._lib.G_at_sigma_arr(sigma.cast(), len(sigma), + d, e, f, g, G.cast()) + return G.finish() def n_in_bins(edges, Marr, dndM): """Tinker et al. 2008 appendix C binned mass function. @@ -125,15 +102,16 @@ def n_in_bins(edges, Marr, dndM): numpy.ndarray: number density of halos in the mass bins. Length is :code:`len(edges)-1`. """ - edges = np.asarray(edges) - if edges.ndim == 0: - edges = edges[None] - if edges.ndim > 1: - raise Exception("edges cannot be a >1D array.") + edges = _ArrayWrapper(edges, 'edges') - n = np.zeros(len(edges)-1) - cluster_toolkit._lib.n_in_bins(_dcast(edges), len(edges), _dcast(Marr), _dcast(dndM), len(Marr), _dcast(n)) - return n + n = _ArrayWrapper.zeros(len(edges)-1) + Marr = _ArrayWrapper(Marr, 'Marr') + dndM = _ArrayWrapper(dndM, 'dndM') + rc = cluster_toolkit._lib.n_in_bins(edges.cast(), len(edges), + Marr.cast(), dndM.cast(), len(Marr), + n.cast()) + _handle_gsl_error(rc, n_in_bins) + return n.finish() def n_in_bin(Mlo, Mhi, Marr, dndM): """Tinker et al. 2008 appendix C binned mass function. @@ -151,6 +129,12 @@ def n_in_bin(Mlo, Mhi, Marr, dndM): return np.squeeze(n_in_bins([Mlo, Mhi], Marr, dndM)) def _dndM_sigma2_precomputed(M, sigma2, dsigma2dM, Omega_m, d=1.97, e=1.0, f=0.51, g=1.228): - dndM = np.zeros_like(M) - cluster_toolkit._lib.dndM_sigma2_precomputed(_dcast(M), _dcast(sigma2), _dcast(dsigma2dM), len(M), Omega_m, d, e, f, g, _dcast(dndM)) - return dndM + M = _ArrayWrapper(M, allow_multidim=True) + sigma2 = _ArrayWrapper(sigma2, allow_multidim=True) + dsigma2dM = _ArrayWrapper(dsigma2dM, allow_multidim=True) + dndM = _ArrayWrapper.zeros_like(M) + cluster_toolkit._lib.dndM_sigma2_precomputed(M.cast(), sigma2.cast(), + dsigma2dM.cast(), len(M), + Omega_m, d, e, f, g, + dndM.cast()) + return dndM.finish() diff --git a/cluster_toolkit/miscentering.py b/cluster_toolkit/miscentering.py index ec3a7ac..bcb6852 100644 --- a/cluster_toolkit/miscentering.py +++ b/cluster_toolkit/miscentering.py @@ -2,11 +2,11 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper import numpy as np def Sigma_mis_single_at_R(R, Rsigma, Sigma, M, conc, Omega_m, Rmis, delta=200): - """Miscentered surface mass density [Msun h/pc^2 comoving] of a profile miscentered by an + """Miscentered surface mass density [Msun h/pc^2 comoving] of a profile miscentered by an amount Rmis Mpc/h comoving. Units are Msun h/pc^2 comoving. Args: @@ -23,26 +23,29 @@ def Sigma_mis_single_at_R(R, Rsigma, Sigma, M, conc, Omega_m, Rmis, delta=200): float or array like: Miscentered projected surface mass density. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - if np.min(R) < np.min(Rsigma): + R = _ArrayWrapper(R, 'R') + + if np.min(R.arr) < np.min(Rsigma): raise Exception("Minimum R must be >= min(R_Sigma)") - if np.max(R) > np.max(Rsigma): + if np.max(R.arr) > np.max(Rsigma): raise Exception("Maximum R must be <= max(R_Sigma)") - Sigma_mis = np.zeros_like(R) - cluster_toolkit._lib.Sigma_mis_single_at_R_arr(_dcast(R), len(R), _dcast(Rsigma), _dcast(Sigma), - len(Rsigma), M, conc, delta, Omega_m, Rmis, _dcast(Sigma_mis)) - if scalar_input: - return np.squeeze(Sigma_mis) - return Sigma_mis + + Rsigma = _ArrayWrapper(Rsigma, allow_multidim=True) + Sigma = _ArrayWrapper(Sigma, allow_multidim=True) + + if Rsigma.shape != Sigma.shape: + raise ValueError('Rsigma and Sigma must have the same shape') + + Sigma_mis = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.Sigma_mis_single_at_R_arr(R.cast(), len(R), + Rsigma.cast(), Sigma.cast(), + len(Rsigma), M, conc, delta, + Omega_m, Rmis, + Sigma_mis.cast()) + return Sigma_mis.finish() def Sigma_mis_at_R(R, Rsigma, Sigma, M, conc, Omega_m, Rmis, delta=200, kernel="rayleigh"): - """Miscentered surface mass density [Msun h/pc^2 comoving] + """Miscentered surface mass density [Msun h/pc^2 comoving] convolved with a distribution for Rmis. Units are Msun h/pc^2 comoving. Args: @@ -60,18 +63,12 @@ def Sigma_mis_at_R(R, Rsigma, Sigma, M, conc, Omega_m, Rmis, delta=200, kernel=" float or array like: Miscentered projected surface mass density. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - - #Exception checking - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - if np.min(R) < np.min(Rsigma): + R = _ArrayWrapper(R, 'R') + + # Exception checking + if np.min(R.arr) < np.min(Rsigma): raise Exception("Minimum R must be >= min(R_Sigma)") - if np.max(R) > np.max(Rsigma): + if np.max(R.arr) > np.max(Rsigma): raise Exception("Maximum R must be <= max(R_Sigma)") if kernel == "rayleigh": integrand_switch = 0 @@ -81,12 +78,18 @@ def Sigma_mis_at_R(R, Rsigma, Sigma, M, conc, Omega_m, Rmis, delta=200, kernel=" raise Exception("Miscentering kernel must be either "+ "'rayleigh' or 'gamma'") - Sigma_mis = np.zeros_like(R) - cluster_toolkit._lib.Sigma_mis_at_R_arr(_dcast(R), len(R), _dcast(Rsigma), _dcast(Sigma), len(Rsigma), - M, conc, delta, Omega_m, Rmis, integrand_switch, _dcast(Sigma_mis)) - if scalar_input: - return np.squeeze(Sigma_mis) - return Sigma_mis + Rsigma = _ArrayWrapper(Rsigma, allow_multidim=True) + Sigma = _ArrayWrapper(Sigma, allow_multidim=True) + + if Rsigma.shape != Sigma.shape: + raise ValueError('Rsigma and Sigma must have the same shape') + + Sigma_mis = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.Sigma_mis_at_R_arr(R.cast(), len(R), Rsigma.cast(), + Sigma.cast(), len(Rsigma), + M, conc, delta, Omega_m, Rmis, + integrand_switch, Sigma_mis.cast()) + return Sigma_mis.finish() def DeltaSigma_mis_at_R(R, Rsigma, Sigma_mis): """Miscentered excess surface mass density profile at R. Units are Msun h/pc^2 comoving. @@ -100,21 +103,22 @@ def DeltaSigma_mis_at_R(R, Rsigma, Sigma_mis): float array like: Miscentered excess surface mass density profile. """ - R = np.asarray(R) - scalar_input = False - if R.ndim == 0: - R = R[None] #makes R 1D - scalar_input = True - if R.ndim > 1: - raise Exception("R cannot be a >1D array.") - if np.min(R) < np.min(Rsigma): + R = _ArrayWrapper(R, 'R') + if np.min(R.arr) < np.min(Rsigma): raise Exception("Minimum R must be >= min(R_Sigma)") - if np.max(R) > np.max(Rsigma): + if np.max(R.arr) > np.max(Rsigma): raise Exception("Maximum R must be <= max(R_Sigma)") - - DeltaSigma_mis = np.zeros_like(R) - cluster_toolkit._lib.DeltaSigma_mis_at_R_arr(_dcast(R), len(R), _dcast(Rsigma), _dcast(Sigma_mis), - len(Rsigma), _dcast(DeltaSigma_mis)) - if scalar_input: - return np.squeeze(DeltaSigma_mis) - return DeltaSigma_mis + + Rsigma = _ArrayWrapper(Rsigma, allow_multidim=True) + Sigma_mis = _ArrayWrapper(Sigma_mis, allow_multidim=True) + + if Rsigma.shape != Sigma_mis.shape: + raise ValueError('Rsigma and Sigma must have the same shape') + + DeltaSigma_mis = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.DeltaSigma_mis_at_R_arr(R.cast(), len(R), + Rsigma.cast(), + Sigma_mis.cast(), + len(Rsigma), + DeltaSigma_mis.cast()) + return DeltaSigma_mis.finish() diff --git a/cluster_toolkit/peak_height.py b/cluster_toolkit/peak_height.py index 72c3e72..ac72b90 100644 --- a/cluster_toolkit/peak_height.py +++ b/cluster_toolkit/peak_height.py @@ -2,9 +2,8 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper, _handle_gsl_error import numpy as np -from numpy import ascontiguousarray as ACA def sigma2_at_M(M, k, P, Omega_m): """RMS variance in top hat sphere of lagrangian radius R [Mpc/h comoving] corresponding to a mass M [Msun/h] of linear power spectrum. @@ -19,12 +18,16 @@ def sigma2_at_M(M, k, P, Omega_m): float or array like: RMS variance of top hat sphere. """ - if type(M) is list or type(M) is np.ndarray: - s2 = np.zeros_like(M) - cluster_toolkit._lib.sigma2_at_M_arr(_dcast(M), len(M), _dcast(k), _dcast(P), len(k), Omega_m, _dcast(s2)) - return s2 + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + if isinstance(M, list) or isinstance(M, np.ndarray): + M = _ArrayWrapper(M, allow_multidim=True) + s2 = _ArrayWrapper.zeros_like(M) + rc = cluster_toolkit._lib.sigma2_at_M_arr(M.cast(), len(M), k.cast(), P.cast(), len(k), Omega_m, s2.cast()) + _handle_gsl_error(rc, sigma2_at_M) + return s2.finish() else: - return cluster_toolkit._lib.sigma2_at_M(M, _dcast(k), _dcast(P), len(k), Omega_m) + return cluster_toolkit._lib.sigma2_at_M(M, k.cast(), P.cast(), len(k), Omega_m) def sigma2_at_R(R, k, P): """RMS variance in top hat sphere of radius R [Mpc/h comoving] of linear power spectrum. @@ -38,12 +41,16 @@ def sigma2_at_R(R, k, P): float or array like: RMS variance of a top hat sphere. """ - if type(R) is list or type(R) is np.ndarray: - s2 = np.zeros_like(R) - cluster_toolkit._lib.sigma2_at_R_arr(_dcast(R), len(R), _dcast(k), _dcast(P), len(k), _dcast(s2)) - return s2 + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + if isinstance(R, list) or isinstance(R, np.ndarray): + R = _ArrayWrapper(R) + s2 = _ArrayWrapper.zeros_like(R) + rc = cluster_toolkit._lib.sigma2_at_R_arr(R.cast(), len(R), k.cast(), P.cast(), len(k), s2.cast()) + _handle_gsl_error(rc, sigma2_at_R) + return s2.finish() else: - return cluster_toolkit._lib.sigma2_at_R(R, _dcast(k), _dcast(P), len(k)) + return cluster_toolkit._lib.sigma2_at_R(R, k.cast(), P.cast(), len(k)) def nu_at_M(M, k, P, Omega_m): """Peak height of top hat sphere of lagrangian radius R [Mpc/h comoving] corresponding to a mass M [Msun/h] of linear power spectrum. @@ -58,12 +65,15 @@ def nu_at_M(M, k, P, Omega_m): nu (float or array like): Peak height. """ - if type(M) is list or type(M) is np.ndarray: - nu = np.zeros_like(M) - cluster_toolkit._lib.nu_at_M_arr(_dcast(M), len(M), _dcast(k), _dcast(P), len(k), Omega_m, _dcast(nu)) - return nu + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + if isinstance(M, list) or isinstance(M, np.ndarray): + M = _ArrayWrapper(M) + nu = _ArrayWrapper.zeros_like(M) + cluster_toolkit._lib.nu_at_M_arr(M.cast(), len(M), k.cast(), P.cast(), len(k), Omega_m, nu.cast()) + return nu.finish() else: - return cluster_toolkit._lib.nu_at_M(M, _dcast(k), _dcast(P), len(k), Omega_m) + return cluster_toolkit._lib.nu_at_M(M, k.cast(), P.cast(), len(k), Omega_m) def nu_at_R(R, k, P): """Peak height of top hat sphere of radius R [Mpc/h comoving] of linear power spectrum. @@ -77,16 +87,19 @@ def nu_at_R(R, k, P): float or array like: Peak height. """ - if type(R) is list or type(R) is np.ndarray: - nu = np.zeros_like(R) - cluster_toolkit._lib.nu_at_R_arr(_dcast(R), len(R), _dcast(k), _dcast(P), len(k), _dcast(nu)) - return nu + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + if isinstance(R, list) or isinstance(R, np.ndarray): + R = _ArrayWrapper(R) + nu = _ArrayWrapper.zeros_like(R) + cluster_toolkit._lib.nu_at_R_arr(R.cast(), len(R), k.cast(), P.cast(), len(k), nu.cast()) + return nu.finish() else: - return cluster_toolkit._lib.nu_at_R(R, _dcast(k), _dcast(P), len(k)) + return cluster_toolkit._lib.nu_at_R(R, k.cast(), P.cast(), len(k)) def dsigma2dM_at_M(M, k, P, Omega_m): - """Derivative w.r.t. mass of RMS variance in top hat sphere of - lagrangian radius R [Mpc/h comoving] corresponding to a mass + """Derivative w.r.t. mass of RMS variance in top hat sphere of + lagrangian radius R [Mpc/h comoving] corresponding to a mass M [Msun/h] of linear power spectrum. Args: @@ -99,44 +112,61 @@ def dsigma2dM_at_M(M, k, P, Omega_m): float or array like: d/dM of RMS variance of top hat sphere. """ - if type(M) is list or type(M) is np.ndarray: - ds2dM = np.zeros_like(M) - cluster_toolkit._lib.dsigma2dM_at_M_arr(_dcast(M), len(M), _dcast(k), - _dcast(P), len(k), Omega_m, - _dcast(ds2dM)) - return ds2dM + P = _ArrayWrapper(P, allow_multidim=True) + k = _ArrayWrapper(k, allow_multidim=True) + if isinstance(M, list) or isinstance(M, np.ndarray): + M = _ArrayWrapper(M, allow_multidim=True) + ds2dM = _ArrayWrapper.zeros_like(M) + rc = cluster_toolkit._lib.dsigma2dM_at_M_arr(M.cast(), len(M), k.cast(), + P.cast(), len(k), Omega_m, + ds2dM.cast()) + _handle_gsl_error(rc, dsigma2dM_at_M) + return ds2dM.finish() else: - return cluster_toolkit._lib.dsigma2dM_at_M(M, _dcast(k), _dcast(P), + return cluster_toolkit._lib.dsigma2dM_at_M(M, k.cast(), P.cast(), len(k), Omega_m) - return 0 - def _calc_sigma2_at_R(R, k, P, s2): - """Direct call to vectorized version of RMS variance in top hat + """Direct call to vectorized version of RMS variance in top hat sphere of radius R [Mpc/h comoving] of linear power spectrum. """ - cluster_toolkit._lib.sigma2_at_R_arr(_dcast(R), len(R), _dcast(k), _dcast(P), len(k), _dcast(s2)) - return + R = _ArrayWrapper(R, allow_multidim=True) + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + s2 = _ArrayWrapper(s2, allow_multidim=True) + cluster_toolkit._lib.sigma2_at_R_arr(R.cast(), len(R), k.cast(), P.cast(), len(k), s2.cast()) def _calc_sigma2_at_M(M, k, P, Omega_m, s2): """Direct call to vectorized version of RMS variance in top hat sphere of lagrangian radius R [Mpc/h comoving] corresponding to a mass M [Msun/h] of linear power spectrum. """ - cluster_toolkit._lib.sigma2_at_M_arr(_dcast(M), len(M), _dcast(k), _dcast(P), len(k), Omega_m, _dcast(s2)) - return + M = _ArrayWrapper(M, allow_multidim=True) + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + s2 = _ArrayWrapper(s2, allow_multidim=True) + rc = cluster_toolkit._lib.sigma2_at_M_arr(M.cast(), len(M), k.cast(), P.cast(), len(k), Omega_m, s2.cast()) + _handle_gsl_error(rc, _calc_sigma2_at_M) def _calc_nu_at_R(R, k, P, nu): """Direct call to vectorized version of peak height of R. """ - cluster_toolkit._lib.nu_at_R_arr(_dcast(R), len(R), _dcast(k), _dcast(P), len(k), _dcast(nu)) - return + R = _ArrayWrapper(R, allow_multidim=True) + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + nu = _ArrayWrapper(nu, allow_multidim=True) + rc = cluster_toolkit._lib.nu_at_R_arr(R.cast(), len(R), k.cast(), P.cast(), len(k), nu.cast()) + _handle_gsl_error(rc, _calc_nu_at_R) def _calc_nu_at_M(M, k, P, Omega_m, nu): """Direct call to vectorized version of peak height of M. """ - cluster_toolkit._lib.nu_at_M_arr(_dcast(M), len(M), _dcast(k), _dcast(P), len(k), Omega_m, _dcast(nu)) - return + M = _ArrayWrapper(M, allow_multidim=True) + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + nu = _ArrayWrapper(nu, allow_multidim=True) + rc = cluster_toolkit._lib.nu_at_M_arr(M.cast(), len(M), k.cast(), P.cast(), len(k), Omega_m, nu.cast()) + _handle_gsl_error(rc, _calc_nu_at_M) diff --git a/cluster_toolkit/profile_derivatives.py b/cluster_toolkit/profile_derivatives.py index fd3c38e..ed89ae0 100644 --- a/cluster_toolkit/profile_derivatives.py +++ b/cluster_toolkit/profile_derivatives.py @@ -2,9 +2,9 @@ Derivatives of halo profiles. Used to plot splashback results. """ import cluster_toolkit as ct -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper import numpy as np -from numpy import ascontiguousarray as ACA + def drho_nfw_dr_at_R(Radii, Mass, conc, Omega_m, delta=200): """Derivative of the NFW halo density profile. @@ -20,12 +20,11 @@ def drho_nfw_dr_at_R(Radii, Mass, conc, Omega_m, delta=200): float or array like: derivative of the NFW profile. """ - if type(Radii) is list or type(Radii) is np.ndarray: - drhodr = np.zeros_like(Radii) - ct._lib.drho_nfw_dr_at_R_arr(_dcast(R), len(R), Mass, conc, - delta, Omega_m, _dcast(drhodr)) - return xi + Radii = _ArrayWrapper(Radii, allow_multidim=True) + if isinstance(Radii, list) or isinstance(Radii, np.ndarray): + drhodr = _ArrayWrapper.zeros_like(Radii) + ct._lib.drho_nfw_dr_at_R_arr(Radii.cast(), len(Radii), Mass, conc, + delta, Omega_m, drhodr.cast()) + return drhodr.finish() else: - return cluster_toolkit._lib.drho_nfw_dr_at_R(Radii, Mass, conc, - delta, Omega_m) - + return ct._lib.drho_nfw_dr_at_R(Radii, Mass, conc, delta, Omega_m) diff --git a/cluster_toolkit/sigma_reconstruction.py b/cluster_toolkit/sigma_reconstruction.py index 7a85c05..2249b5b 100644 --- a/cluster_toolkit/sigma_reconstruction.py +++ b/cluster_toolkit/sigma_reconstruction.py @@ -2,7 +2,7 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper import numpy as np def Sigma_REC_from_DeltaSigma(R, DeltaSigma): @@ -22,8 +22,8 @@ def Sigma_REC_from_DeltaSigma(R, DeltaSigma): Returns: Reconstructed surface mass density. """ - R = np.asarray(R) - DeltaSigma = np.asarray(DeltaSigma) + R = _ArrayWrapper(R, allow_multidim=True) + DeltaSigma = _ArrayWrapper(DeltaSigma, allow_multidim=True) if R.shape != DeltaSigma.shape: raise Exception("R and DeltaSigma must have the same shape.") @@ -37,8 +37,8 @@ def Sigma_REC_from_DeltaSigma(R, DeltaSigma): #Note the order here, we integrate DOWNWARD dlnR = lnR[0] - lnR[1] - - Sigma = np.zeros(len(R)-1) - cluster_toolkit._lib.Sigma_REC_from_DeltaSigma(dlnR, _dcast(DeltaSigma), - len(R), _dcast(Sigma)) - return Sigma + + Sigma = _ArrayWrapper.zeros(len(R)-1) + cluster_toolkit._lib.Sigma_REC_from_DeltaSigma(dlnR, DeltaSigma.cast(), + len(R), Sigma.cast()) + return Sigma.finish() diff --git a/cluster_toolkit/xi.py b/cluster_toolkit/xi.py index 3756426..12d6fdb 100644 --- a/cluster_toolkit/xi.py +++ b/cluster_toolkit/xi.py @@ -2,7 +2,7 @@ """ import cluster_toolkit -from cluster_toolkit import _dcast +from cluster_toolkit import _ArrayWrapper, _handle_gsl_error import numpy as np def xi_nfw_at_r(r, M, c, Omega_m, delta=200): @@ -19,20 +19,12 @@ def xi_nfw_at_r(r, M, c, Omega_m, delta=200): float or array like: NFW halo profile. """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") + r = _ArrayWrapper(r, 'r') - xi = np.zeros_like(r) - cluster_toolkit._lib.calc_xi_nfw(_dcast(r), len(r), M, c, delta, - Omega_m, _dcast(xi)) - if scalar_input: - return np.squeeze(xi) - return xi + xi = _ArrayWrapper.zeros_like(r) + cluster_toolkit._lib.calc_xi_nfw(r.cast(), len(r), M, c, delta, + Omega_m, xi.cast()) + return xi.finish() def xi_einasto_at_r(r, M, conc, alpha, om, delta=200, rhos=-1.): """Einasto halo profile. @@ -50,21 +42,12 @@ def xi_einasto_at_r(r, M, conc, alpha, om, delta=200, rhos=-1.): float or array like: Einasto halo profile. """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") - - xi = np.zeros_like(r) - cluster_toolkit._lib.calc_xi_einasto(_dcast(r), len(r), M, rhos, - conc, alpha, delta, om, _dcast(xi)) - - if scalar_input: - return np.squeeze(xi) - return xi + r = _ArrayWrapper(r, 'r') + + xi = _ArrayWrapper.zeros_like(r) + cluster_toolkit._lib.calc_xi_einasto(r.cast(), len(r), M, rhos, + conc, alpha, delta, om, xi.cast()) + return xi.finish() def xi_mm_at_r(r, k, P, N=500, step=0.005, exact=False): """Matter-matter correlation function. @@ -81,29 +64,25 @@ def xi_mm_at_r(r, k, P, N=500, step=0.005, exact=False): float or array like: Matter-matter correlation function """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") - - xi = np.zeros_like(r) + r = _ArrayWrapper(r, 'r') + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + + xi = _ArrayWrapper.zeros_like(r) if not exact: - cluster_toolkit._lib.calc_xi_mm(_dcast(r), len(r), _dcast(k), - _dcast(P), len(k), _dcast(xi), - N, step) + rc = cluster_toolkit._lib.calc_xi_mm(r.cast(), len(r), k.cast(), + P.cast(), len(k), xi.cast(), + N, step) + _handle_gsl_error(rc, xi_mm_at_r) else: - if max(r) > 1e3: + if r.arr.max() > 1e3: raise Exception("max(r) cannot be >1e3 for numerical stability.") - cluster_toolkit._lib.calc_xi_mm_exact(_dcast(r), len(r), - _dcast(k), _dcast(P), - len(k), _dcast(xi)) + rc = cluster_toolkit._lib.calc_xi_mm_exact(r.cast(), len(r), + k.cast(), P.cast(), + len(k), xi.cast()) + _handle_gsl_error(rc, xi_mm_at_r) - if scalar_input: - return np.squeeze(xi) - return xi + return xi.finish() def xi_2halo(bias, xi_mm): """2-halo term in halo-matter correlation function @@ -116,10 +95,11 @@ def xi_2halo(bias, xi_mm): float or array like: 2-halo term in halo-matter correlation function """ - xi = np.zeros_like(xi_mm) - cluster_toolkit._lib.calc_xi_2halo(len(xi_mm), bias, _dcast(xi_mm), - _dcast(xi)) - return xi + xi_mm = _ArrayWrapper(xi_mm, allow_multidim=True) + xi = _ArrayWrapper.zeros_like(xi_mm) + cluster_toolkit._lib.calc_xi_2halo(len(xi_mm), bias, xi_mm.cast(), + xi.cast()) + return xi.finish() def xi_hm(xi_1halo, xi_2halo, combination="max"): """Halo-matter correlation function @@ -142,10 +122,12 @@ def xi_hm(xi_1halo, xi_2halo, combination="max"): else: raise Exception("Combinations other than maximum not implemented yet") - xi = np.zeros_like(xi_1halo) - cluster_toolkit._lib.calc_xi_hm(len(xi_1halo), _dcast(xi_1halo), - _dcast(xi_2halo), _dcast(xi), switch) - return xi + xi_1halo = _ArrayWrapper(xi_1halo, allow_multidim=True) + xi_2halo = _ArrayWrapper(xi_2halo, allow_multidim=True) + xi = _ArrayWrapper.zeros_like(xi_1halo) + cluster_toolkit._lib.calc_xi_hm(len(xi_1halo), xi_1halo.cast(), + xi_2halo.cast(), xi.cast(), switch) + return xi.finish() def xi_DK(r, M, conc, be, se, k, P, om, delta=200, rhos=-1., alpha=-1., beta=-1., gamma=-1.): """Diemer-Kravtsov 2014 profile. @@ -169,20 +151,14 @@ def xi_DK(r, M, conc, be, se, k, P, om, delta=200, rhos=-1., alpha=-1., beta=-1. float or array like: DK profile evaluated at the input radii """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") - - xi = np.zeros_like(r) - cluster_toolkit._lib.calc_xi_DK(_dcast(r), len(r), M, rhos, conc, be, se, alpha, beta, gamma, delta, _dcast(k), _dcast(P), len(k), om, _dcast(xi)) + r = _ArrayWrapper(r, 'r') + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + + xi = _ArrayWrapper.zeros_like(r) + cluster_toolkit._lib.calc_xi_DK(r.cast(), len(r), M, rhos, conc, be, se, alpha, beta, gamma, delta, k.cast(), P.cast(), len(k), om, xi.cast()) - if scalar_input: - return np.squeeze(xi) - return xi + return xi.finish() def xi_DK_appendix1(r, M, conc, be, se, k, P, om, bias, xi_mm, delta=200, rhos=-1., alpha=-1., beta=-1., gamma=-1.): """Diemer-Kravtsov 2014 profile, first form from the appendix, eq. A3. @@ -208,20 +184,15 @@ def xi_DK_appendix1(r, M, conc, be, se, k, P, om, bias, xi_mm, delta=200, rhos=- float or array like: DK profile evaluated at the input radii """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") + r = _ArrayWrapper(r, 'r') + k = _ArrayWrapper(k, allow_multidim=True) + P = _ArrayWrapper(P, allow_multidim=True) + xi_mm = _ArrayWrapper(xi_mm, allow_multidim=True) xi = np.zeros_like(r) - cluster_toolkit._lib.calc_xi_DK_app1(_dcast(r), len(r), M, rhos, conc, be, se, alpha, beta, gamma, delta, _dcast(k), _dcast(P), len(k), om, bias, _dcast(xi_mm), _dcast(xi)) + cluster_toolkit._lib.calc_xi_DK_app1(r.cast(), len(r), M, rhos, conc, be, se, alpha, beta, gamma, delta, k.cast(), P.cast(), len(k), om, bias, xi_mm.cast(), xi.cast()) - if scalar_input: - return np.squeeze(xi) - return xi + return xi.finish() def xi_DK_appendix2(r, M, conc, be, se, k, P, om, bias, xi_mm, delta=200, rhos=-1., alpha=-1., beta=-1., gamma=-1.): """Diemer-Kravtsov 2014 profile, second form from the appendix, eq. A4. @@ -246,17 +217,14 @@ def xi_DK_appendix2(r, M, conc, be, se, k, P, om, bias, xi_mm, delta=200, rhos=- Returns: float or array like: DK profile evaluated at the input radii """ - r = np.asarray(r) - scalar_input = False - if r.ndim == 0: - r = r[None] #makes r 1D - scalar_input = True - if r.ndim > 1: - raise Exception("r cannot be a >1D array.") - - xi = np.zeros_like(r) - cluster_toolkit._lib.calc_xi_DK_app2(_dcast(r), len(r), M, rhos, conc, be, se, alpha, beta, gamma, delta, _dcast(k), _dcast(P), len(k), om, bias, _dcast(xi_mm), _dcast(xi)) - - if scalar_input: - return np.squeeze(xi) - return xi + r = _ArrayWrapper(r, 'r') + k = _ArrayWrapper(k) + P = _ArrayWrapper(P) + xi_mm = _ArrayWrapper(xi_mm) + + xi = _ArrayWrapper.zeros_like(r) + cluster_toolkit._lib.calc_xi_DK_app2(r.cast(), len(r), M, rhos, conc, be, + se, alpha, beta, gamma, delta, + k.cast(), P.cast(), len(k), om, bias, + xi_mm.cast(), xi.cast()) + return xi.finish() diff --git a/include/C_bias.h b/include/C_bias.h index 7679599..e541166 100644 --- a/include/C_bias.h +++ b/include/C_bias.h @@ -1,13 +1,13 @@ -int bias_at_nu_arr(double*nu, int Nnu, int delta, double*bias); -int bias_at_R_arr(double*R, int NR, int delta, double*k, double*P, int Nk, +void bias_at_nu_arr(double*nu, int Nnu, int delta, double*bias); +void bias_at_R_arr(double*R, int NR, int delta, double*k, double*P, int Nk, double*bias); -int bias_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, +void bias_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, double Omega_m, double*bias); -int bias_at_nu_arr_FREEPARAMS(double*nu, int Nnu, int delta, double A, +void bias_at_nu_arr_FREEPARAMS(double*nu, int Nnu, int delta, double A, double a, double B, double b, double C, double c, double*bias); -int dbiasdnu_at_nu_arr(double*nu, int Nnu, int delta, double*deriv); -int dbiasdM_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, +void dbiasdnu_at_nu_arr(double*nu, int Nnu, int delta, double*deriv); +void dbiasdM_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, double Omega_m, double*deriv); diff --git a/include/C_boostfactors.h b/include/C_boostfactors.h index 35a199d..8127051 100644 --- a/include/C_boostfactors.h +++ b/include/C_boostfactors.h @@ -1,3 +1,3 @@ -int boost_nfw_at_R_arr(double*R, int NR, double B0, double Rs, double*boost); -int boost_powerlaw_at_R_arr(double*R, int NR, double B0, double Rs, +void boost_nfw_at_R_arr(double*R, int NR, double B0, double Rs, double*boost); +void boost_powerlaw_at_R_arr(double*R, int NR, double B0, double Rs, double alpha, double*boost); diff --git a/include/C_deltasigma.h b/include/C_deltasigma.h index c5fee82..58510dc 100644 --- a/include/C_deltasigma.h +++ b/include/C_deltasigma.h @@ -1,5 +1,5 @@ double Sigma_nfw_at_R(double R, double M, double c, int delta, double om); -int Sigma_nfw_at_R_arr(double*R, int NR, double M,double c, int delta, double om, double*Sigma); +void Sigma_nfw_at_R_arr(double*R, int NR, double M,double c, int delta, double om, double*Sigma); int Sigma_at_R_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double M, double conc, int delta, double om, double*Sigma); diff --git a/include/C_density.h b/include/C_density.h index 4d665e0..be99d9b 100644 --- a/include/C_density.h +++ b/include/C_density.h @@ -1,2 +1,2 @@ -int calc_rho_nfw(double*, int, double, double, int, double, double*); -int calc_rho_einasto(double*, int, double, double, double, double, int, double, double*); +void calc_rho_nfw(double*, int, double, double, int, double, double*); +void calc_rho_einasto(double*, int, double, double, double, double, int, double, double*); diff --git a/include/C_error.h b/include/C_error.h new file mode 100644 index 0000000..69e7878 --- /dev/null +++ b/include/C_error.h @@ -0,0 +1 @@ +void *gsl_set_error_handler_off (void); diff --git a/include/C_exclusion.h b/include/C_exclusion.h index 6dc2bb9..b4bb31b 100644 --- a/include/C_exclusion.h +++ b/include/C_exclusion.h @@ -1,4 +1,4 @@ -int xi_hm_exclusion_at_r_arr(double*r, int Nr, +void xi_hm_exclusion_at_r_arr(double*r, int Nr, double M, double c, double alpha, double rt, double D, double r_eff, double D_eff, @@ -6,14 +6,14 @@ int xi_hm_exclusion_at_r_arr(double*r, int Nr, double bias, double*ximm, int delta, double Omega_m, double*xihm); -int xi_1h_at_r_arr(double*r, int Nr, double M, double c, double alpha, +void xi_1h_at_r_arr(double*r, int Nr, double M, double c, double alpha, double rt, double D, int delta, double Omega_m, double*xi_1h); -int xi_2h_at_r_arr(double*r, int Nr, double r_eff, double D_eff, +void xi_2h_at_r_arr(double*r, int Nr, double r_eff, double D_eff, double bias, double*ximm, double*xi2h); -int xi_C_at_r_arr(double*r, int Nr, double r_A, double r_B, double D, +void xi_C_at_r_arr(double*r, int Nr, double r_A, double r_B, double D, double*xi_2h, double*xi_C); int theta_erfc_at_r_arr(double*r, int Nr, double rt, double D, diff --git a/include/C_massfunction.h b/include/C_massfunction.h index d296da1..f2f42e3 100644 --- a/include/C_massfunction.h +++ b/include/C_massfunction.h @@ -1,16 +1,16 @@ -int G_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, +void G_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double d, double e, double f, double g, double*G); -int G_at_sigma_arr(double*sigma, int Ns, double d, double e, +void G_at_sigma_arr(double*sigma, int Ns, double d, double e, double f, double g, double*G); -int dndM_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, +void dndM_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double d, double e, double f, double g, double*dndM); int n_in_bins(double*edges, int Nedges, double*M, double*dndM, int NM, double*N); //int dndM_sigma2_precomputed(double*M, double*sigma2, double*sigma2_top, double*sigma2_bot, int NM, double Omega_m, double d, double e, double f, double g, double*dndM); -int dndM_sigma2_precomputed(double*M, double*sigma2, double*dsigma2dM, +void dndM_sigma2_precomputed(double*M, double*sigma2, double*dsigma2dM, int NM, double Omega_m, double d, double e, double f, double g, double*dndM); diff --git a/include/C_sigma_reconstruction.h b/include/C_sigma_reconstruction.h index 0fd81ef..cc2a85f 100644 --- a/include/C_sigma_reconstruction.h +++ b/include/C_sigma_reconstruction.h @@ -1,2 +1,2 @@ -int Sigma_REC_from_DeltaSigma(double dlnR, double*DeltaSigma, int N, +void Sigma_REC_from_DeltaSigma(double dlnR, double*DeltaSigma, int N, double*Sigma); diff --git a/include/C_xi.h b/include/C_xi.h index fc71de4..2835627 100644 --- a/include/C_xi.h +++ b/include/C_xi.h @@ -1,18 +1,18 @@ double xi_nfw_at_r(double, double, double, int, double); -int calc_xi_nfw(double*, int, double, double, int, double, double*); +void calc_xi_nfw(double*, int, double, double, int, double, double*); -int calc_xi_einasto(double*, int, double, double, double, double, int, double, double*); +void calc_xi_einasto(double*, int, double, double, double, double, int, double, double*); double rhos_einasto_at_M(double Mass, double conc, double alpha, int delta, double om); double xi_mm_at_r_exact(double r, double*k, double*P, int Nk); int calc_xi_mm_exact(double*r, int Nr, double*k, double*P, int Nk, double*xi); -int calc_xi_2halo(int, double, double*, double*); -int calc_xi_hm(int, double*, double*, double*, int); +void calc_xi_2halo(int, double, double*, double*); +void calc_xi_hm(int, double*, double*, double*, int); int calc_xi_mm(double*, int, double*, double*, int, double*, int, double); -int calc_xi_DK(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double*xi); +void calc_xi_DK(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double*xi); -int calc_xi_DK_app1(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi); +void calc_xi_DK_app1(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi); -int calc_xi_DK_app2(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi); +void calc_xi_DK_app2(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi); diff --git a/src/C_averaging.c b/src/C_averaging.c index 78c9906..febfc43 100644 --- a/src/C_averaging.c +++ b/src/C_averaging.c @@ -11,6 +11,7 @@ #include "C_averaging.h" +#include "gsl/gsl_errno.h" #include "gsl/gsl_integration.h" #include "gsl/gsl_spline.h" #include @@ -25,14 +26,19 @@ typedef struct integrand_params{ gsl_spline*spline; gsl_interp_accel*acc; + int retcode; }integrand_params; double ave_integrand(double lR, void*params){ double R = exp(lR); - integrand_params pars = *(integrand_params *)params; - gsl_spline*spline = pars.spline; - gsl_interp_accel*acc = pars.acc; - return R*R*gsl_spline_eval(spline, R, acc); + integrand_params *pars = (integrand_params *)params; + gsl_spline*spline = pars->spline; + gsl_interp_accel*acc = pars->acc; + double result = 0.0; + int rc = gsl_spline_eval_e(spline, R, acc, &result); + if (pars->retcode == GSL_SUCCESS) + pars->retcode = rc; + return R*R*result; } int average_profile_in_bins(double*Redges, int Nedges, double*R, int NR, @@ -41,27 +47,41 @@ int average_profile_in_bins(double*Redges, int Nedges, double*R, int NR, gsl_interp_accel *acc = gsl_interp_accel_alloc(); gsl_integration_workspace *ws = gsl_integration_workspace_alloc(workspace_size); gsl_function F; + + if (!spline || !acc || !ws) + return GSL_FAILURE; + double result, err; + int rc = GSL_SUCCESS; gsl_spline_init(spline, R, profile, NR); integrand_params params; params.acc = acc; params.spline = spline; + params.retcode = GSL_SUCCESS; F.params = ¶ms; F.function = &ave_integrand; //Loop over bins and compute the average int i; for(i = 0; i < Nedges-1; i++){ - gsl_integration_qag(&F, log(Redges[i]), log(Redges[i+1]), ABSERR, RELERR, + rc = gsl_integration_qag(&F, log(Redges[i]), log(Redges[i+1]), ABSERR, RELERR, workspace_size, 6, ws, &result, &err); ave_profile[i] = 2*result/(Redges[i+1]*Redges[i+1]-Redges[i]*Redges[i]); + + if (rc != GSL_SUCCESS) + break; + if (params.retcode != GSL_SUCCESS) + break; } //Free everything gsl_spline_free(spline); gsl_interp_accel_free(acc); gsl_integration_workspace_free(ws); - return 0; + + if (params.retcode != GSL_SUCCESS) + return params.retcode; + return rc; } diff --git a/src/C_bias.c b/src/C_bias.c index 4280ac2..dbee6c2 100644 --- a/src/C_bias.c +++ b/src/C_bias.c @@ -27,7 +27,7 @@ * * This is the Tinker et al. (2010) bias model. */ -int bias_at_nu_arr(double*nu, int Nnu, int delta, double*bias){ +void bias_at_nu_arr(double*nu, int Nnu, int delta, double*bias){ double y = log10(delta); double xp = exp(-1.0*pow(4./y,4.)); double A = 1.+0.24*y*xp, a = 0.44*y-0.88; @@ -36,7 +36,6 @@ int bias_at_nu_arr(double*nu, int Nnu, int delta, double*bias){ int i; for(i = 0; i < Nnu; i++) bias[i] = 1 - A*pow(nu[i],a)/(pow(nu[i],a)+pow(delta_c,a)) + B*pow(nu[i],b) + C*pow(nu[i],c); - return 0; } /** @@ -45,13 +44,12 @@ int bias_at_nu_arr(double*nu, int Nnu, int delta, double*bias){ * * This is the Tinker et al. (2010) bias model. */ -int bias_at_R_arr(double*R, int NR, int delta, double*k, double*P, int Nk, +void bias_at_R_arr(double*R, int NR, int delta, double*k, double*P, int Nk, double*bias){ double*nu = malloc(sizeof(double)*NR); nu_at_R_arr(R, NR, k, P, Nk, nu); bias_at_nu_arr(nu, NR, delta, bias); free(nu); - return 0; } /** @@ -60,13 +58,12 @@ int bias_at_R_arr(double*R, int NR, int delta, double*k, double*P, int Nk, * * This is the Tinker et al. (2010) bias model. */ -int bias_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, +void bias_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, double Omega_m, double*bias){ double*nu = malloc(sizeof(double)*NM); nu_at_M_arr(M, NM, k, P, Nk, Omega_m, nu); bias_at_nu_arr(nu, NM, delta, bias); free(nu); - return 0; } /** @@ -74,13 +71,12 @@ int bias_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, * of peak heights, with arbitrary free parameters in a Tinker-like model. * */ -int bias_at_nu_arr_FREEPARAMS(double*nu, int Nnu, int delta, +void bias_at_nu_arr_FREEPARAMS(double*nu, int Nnu, int delta, double A, double a, double B, double b, double C, double c, double*bias){ int i; for(i = 0; i < Nnu; i++) bias[i] = 1 - A*pow(nu[i],a)/(pow(nu[i],a)+pow(delta_c,a)) + B*pow(nu[i],b) + C*pow(nu[i],c); - return 0; } /****** @@ -91,7 +87,7 @@ Derivatives of the bias * * This is the Tinker et al. (2010) bias model. */ -int dbiasdnu_at_nu_arr(double*nu, int Nnu, int delta, double*deriv){ +void dbiasdnu_at_nu_arr(double*nu, int Nnu, int delta, double*deriv){ double y = log10(delta); double xp = exp(-1.0*pow(4./y,4.)); double A = 1.+0.24*y*xp, a = 0.44*y-0.88; @@ -102,7 +98,6 @@ int dbiasdnu_at_nu_arr(double*nu, int Nnu, int delta, double*deriv){ deriv[i] = a*A*pow(delta_c, a) * pow(nu[i], a-1)* pow(pow(nu[i], a) + pow(delta_c, a), -2) + B*b*pow(nu[i], b-1) + C*c*pow(nu[i], c-1); - return 0; } /** @@ -110,7 +105,7 @@ int dbiasdnu_at_nu_arr(double*nu, int Nnu, int delta, double*deriv){ * * This is the Tinker et al. (2010) bias model. */ -int dbiasdM_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, +void dbiasdM_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, double Omega_m, double*deriv){ double*nu = malloc(sizeof(double)*NM); double*dbiasdnu = malloc(sizeof(double)*NM); @@ -130,5 +125,4 @@ int dbiasdM_at_M_arr(double*M, int NM, int delta, double*k, double*P, int Nk, free(dbiasdnu); free(sigma2); free(dsigma2dM); - return 0; } diff --git a/src/C_boostfactors.c b/src/C_boostfactors.c index fdbcef7..cb3b067 100644 --- a/src/C_boostfactors.c +++ b/src/C_boostfactors.c @@ -19,7 +19,7 @@ * * Used in McClintock et al. (2018). */ -int boost_nfw_at_R_arr(double*R, int NR, double B0, double Rs, +void boost_nfw_at_R_arr(double*R, int NR, double B0, double Rs, double*boost){ int i; double x2m1; @@ -36,8 +36,6 @@ int boost_nfw_at_R_arr(double*R, int NR, double B0, double Rs, boost[i] = 1; } } - - return 0; } /** @@ -45,11 +43,10 @@ int boost_nfw_at_R_arr(double*R, int NR, double B0, double Rs, * * Used in Melchior et al. (2017). */ -int boost_powerlaw_at_R_arr(double*R, int NR, double B0, double Rs, +void boost_powerlaw_at_R_arr(double*R, int NR, double B0, double Rs, double alpha, double*boost){ int i; for(i = 0; i < NR; i++){ boost[i] = 1+B0*pow(R[i]/Rs,alpha); } - return 0; } diff --git a/src/C_deltasigma.c b/src/C_deltasigma.c index e8147d0..dba4583 100644 --- a/src/C_deltasigma.c +++ b/src/C_deltasigma.c @@ -48,7 +48,7 @@ double Sigma_nfw_at_R(double R, double M, double c, int delta, double om){ * * Note: all distances are comoving. */ -int Sigma_nfw_at_R_arr(double*R, int NR, double M, double c, int delta, double om, double*Sigma){ +void Sigma_nfw_at_R_arr(double*R, int NR, double M, double c, int delta, double om, double*Sigma){ double rhom = om*rhocrit;//SM h^2/Mpc^3 double deltac = delta*0.3333333333*c*c*c/(log(1.+c)-c/(1.+c)); double Rdelta = pow(M/(1.333333333*M_PI*rhom*delta),0.333333333);//Mpc/h @@ -65,7 +65,6 @@ int Sigma_nfw_at_R_arr(double*R, int NR, double M, double c, int delta, double o } Sigma[i] = 2*Rscale*deltac*rhom*gx*1.e-12; //SM h/pc^2 } - return 0; } typedef struct integrand_params{ @@ -135,7 +134,6 @@ double integrand_large_scales(double lRz, void*params){ * Note: all distances are comoving. */ int Sigma_at_R_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double M, double conc, int delta, double om, double*Sigma){ - gsl_set_error_handler_off(); double rhom = om*rhocrit*1e-12; //SM h^2/pc^2/Mpc; integral is over Mpc/h double Rxi0 = Rxi[0]; double Rxi_max = Rxi[Nxi-1]; @@ -146,15 +144,22 @@ int Sigma_at_R_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double M, d gsl_interp_accel*acc= gsl_interp_accel_alloc(); gsl_integration_workspace*workspace = gsl_integration_workspace_alloc(workspace_size); + + // If allocation fails + if (!spline || !acc || !workspace) + return GSL_ENOMEM; + integrand_params params; gsl_function F; double result1, err1, result2, err2; - int i; + int i, rc = GSL_SUCCESS; double*lnRxi = (double*)malloc(Nxi*sizeof(double)); for(i = 0; i < Nxi; i++){ lnRxi[i] = log(Rxi[i]); } - gsl_spline_init(spline, lnRxi, xi, Nxi); + + rc = gsl_spline_init(spline, lnRxi, xi, Nxi); + params.acc = acc; params.spline = spline; params.M = M; @@ -162,32 +167,32 @@ int Sigma_at_R_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double M, d params.delta = delta; params.om = om; F.params = ¶ms; - int status; for(i = 0; i < NR; i++){ ln_z_max = log(sqrt(Rxi_max*Rxi_max - R[i]*R[i])); //Max distance to integrate to params.Rp = R[i]; if(R[i] < Rxi0){ F.function = &integrand_small_scales; - status = gsl_integration_qag(&F, log(Rxi0)-10, log(sqrt(Rxi0*Rxi0-R[i]*R[i])), ABSERR, RELERR, workspace_size, KEY, workspace, &result1, &err1); - if (status) - fprintf(stderr, "Error in C_deltasigma.c in small scales with\n\t%e\n\t%e\n",M,conc); + rc = gsl_integration_qag(&F, log(Rxi0)-10, log(sqrt(Rxi0*Rxi0-R[i]*R[i])), ABSERR, RELERR, workspace_size, KEY, workspace, &result1, &err1); + if (rc != GSL_SUCCESS) + break; F.function = &integrand_medium_scales; - status = gsl_integration_qag(&F, log(sqrt(Rxi0*Rxi0-R[i]*R[i])), ln_z_max, ABSERR, RELERR, workspace_size, KEY, workspace, &result2, &err2); - if (status) - fprintf(stderr, "Error in C_deltasigma.c in medium scales with\n\t%e\n\t%e\n",M,conc); - + rc = gsl_integration_qag(&F, log(sqrt(Rxi0*Rxi0-R[i]*R[i])), ln_z_max, ABSERR, RELERR, workspace_size, KEY, workspace, &result2, &err2); }else{ //R[i] > Rxi0 result1 = 0; F.function = &integrand_medium_scales; - gsl_integration_qag(&F, -10, ln_z_max, ABSERR, RELERR, workspace_size, KEY, workspace, &result2, &err2); + rc = gsl_integration_qag(&F, -10, ln_z_max, ABSERR, RELERR, workspace_size, KEY, workspace, &result2, &err2); } Sigma[i] = (result1+result2)*rhom*2; + if (rc != GSL_SUCCESS) + break; } + gsl_spline_free(spline); gsl_interp_accel_free(acc); gsl_integration_workspace_free(workspace); free(lnRxi); - return 0; + + return rc; } /** @@ -205,10 +210,15 @@ int Sigma_at_R_full_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double double Rxi_max = Rxi[Nxi-1]; double ln_z_max; double result3, err3; + gsl_integration_workspace*workspace = gsl_integration_workspace_alloc(workspace_size); + // Handle allocation failure + if (!workspace) + return GSL_ENOMEM; + integrand_params params; gsl_function F; - int i; + int i, rc = GSL_SUCCESS; Sigma_at_R_arr(R, NR, Rxi, xi, Nxi, M, conc, delta, om, Sigma); params.slope = log(xi[Nxi-1]/xi[Nxi-2])/log(Rxi[Nxi-1]/Rxi[Nxi-2]); params.intercept = xi[Nxi-1]/pow(Rxi[Nxi-1], params.slope); @@ -217,11 +227,14 @@ int Sigma_at_R_full_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double for(i = 0; i < NR; i++){ ln_z_max = log(sqrt(Rxi_max*Rxi_max - R[i]*R[i])); params.Rp = R[i]; - gsl_integration_qag(&F, ln_z_max, ln_z_max+ulim, ABSERR, RELERR, workspace_size, KEY, workspace, &result3, &err3); + rc = gsl_integration_qag(&F, ln_z_max, ln_z_max+ulim, ABSERR, RELERR, workspace_size, KEY, workspace, &result3, &err3); Sigma[i] += (result3*rhom*2); + if (rc != GSL_SUCCESS) + break; } + gsl_integration_workspace_free(workspace); - return 0; + return rc; } ////////////// DELTASIGMA FUNCTIONS BELOW//////////////// @@ -247,15 +260,22 @@ int DeltaSigma_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, doubl gsl_spline*spline = gsl_spline_alloc(gsl_interp_cspline, Ns); gsl_interp_accel*acc = gsl_interp_accel_alloc(); gsl_integration_workspace* workspace = gsl_integration_workspace_alloc(workspace_size); + + // Handle allocation failures + if (!spline || !acc || !workspace) + return GSL_ENOMEM; + integrand_params params; double result1, result2, err1, err2; gsl_function F; - int i; + int i, rc = GSL_SUCCESS; double*lnRs = (double*)malloc(Ns*sizeof(double)); + for(i = 0; i < Ns; i++){ lnRs[i] = log(Rs[i]); } - gsl_spline_init(spline, lnRs, Sigma, Ns); + + rc = gsl_spline_init(spline, lnRs, Sigma, Ns); params.spline = spline; params.acc = acc; params.M = M; @@ -264,15 +284,28 @@ int DeltaSigma_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, doubl params.om = om; F.params = ¶ms; F.function = &DS_integrand_small_scales; - gsl_integration_qag(&F, lrmin-10, lrmin, ABSERR, RELERR, workspace_size, KEY, workspace, &result1, &err1); - F.function = &DS_integrand_medium_scales; - for(i = 0; i < NR; i++){ - gsl_integration_qag(&F, lrmin, log(R[i]), ABSERR, RELERR, workspace_size, KEY, workspace, &result2, &err2); - DeltaSigma[i] = (result1+result2)*2/(R[i]*R[i]) - gsl_spline_eval(spline, log(R[i]), acc); + if (rc == GSL_SUCCESS) { + rc = gsl_integration_qag(&F, lrmin-10, lrmin, ABSERR, RELERR, workspace_size, KEY, workspace, &result1, &err1); + F.function = &DS_integrand_medium_scales; + if (rc == GSL_SUCCESS) { + for(i = 0; i < NR; i++){ + rc = gsl_integration_qag(&F, lrmin, log(R[i]), ABSERR, RELERR, workspace_size, KEY, workspace, &result2, &err2); + if (rc != GSL_SUCCESS) + break; + + double spline_eval = 0.0; + rc = gsl_spline_eval_e(spline, log(R[i]), acc, &spline_eval); + if (rc != GSL_SUCCESS) + break; + DeltaSigma[i] = (result1+result2)*2/(R[i]*R[i]) - spline_eval; + } + } } + gsl_spline_free(spline); gsl_interp_accel_free(acc); gsl_integration_workspace_free(workspace); free(lnRs); - return 0; + + return rc; } diff --git a/src/C_density.c b/src/C_density.c index 52d38cc..7ca0e02 100644 --- a/src/C_density.c +++ b/src/C_density.c @@ -35,17 +35,16 @@ * @Omega_m Matter fraction. * @return NFW halo density. */ -int calc_rho_nfw(double*r, int Nr, double Mass, double conc, int delta, double Omega_m, double*rho_nfw){ +void calc_rho_nfw(double*r, int Nr, double Mass, double conc, int delta, double Omega_m, double*rho_nfw){ int i; double rhom = Omega_m*rhomconst;//Msun h^2/Mpc^3 calc_xi_nfw(r, Nr, Mass, conc, delta, Omega_m, rho_nfw); //rho_nfw actually holds xi_nfw here for(i = 0; i < Nr; i++){ rho_nfw[i] = rhom*(1+rho_nfw[i]); } - return 0; } -int calc_rho_einasto(double*R, int NR, double Mass, double rhos, double conc, double alpha, int delta, double Omega_m, double*rho_einasto){ +void calc_rho_einasto(double*R, int NR, double Mass, double rhos, double conc, double alpha, int delta, double Omega_m, double*rho_einasto){ double rhom = rhomconst*Omega_m; //SM h^2/Mpc^3 double Rdelta = pow(Mass/(1.33333333333*M_PI*rhom*delta), 0.33333333333); double rs = Rdelta / conc; //compute scale radius from concentration @@ -57,6 +56,4 @@ int calc_rho_einasto(double*R, int NR, double Mass, double rhos, double conc, do x = 2./alpha * pow(R[i]/rs, alpha); rho_einasto[i] = rhos * exp(-x); } - - return 0; } diff --git a/src/C_exclusion.c b/src/C_exclusion.c index 99e85dd..4c0aa56 100644 --- a/src/C_exclusion.c +++ b/src/C_exclusion.c @@ -22,7 +22,7 @@ #define rm_min 0.0001 //Mpc/h minimum of the radial splines #define rm_max 10000. //Mpc/h maximum of the radial splines -int xi_hm_exclusion_at_r_arr(double*r, int Nr, +void xi_hm_exclusion_at_r_arr(double*r, int Nr, double M, double c, double alpha, double rt, double D, double r_eff, double D_eff, @@ -43,10 +43,9 @@ int xi_hm_exclusion_at_r_arr(double*r, int Nr, free(xi_1h); free(xi_2h); free(xi_C); - return 0; } -int xi_1h_at_r_arr(double*r, int Nr, double M, double c, double alpha, +void xi_1h_at_r_arr(double*r, int Nr, double M, double c, double alpha, double rt, double D, int delta, double Omega_m, double*xi_1h){ int i; @@ -57,10 +56,9 @@ int xi_1h_at_r_arr(double*r, int Nr, double M, double c, double alpha, xi_1h[i] = (1+xi_1h[i]) * theta[i]; } free(theta); - return 0; //success } -int xi_2h_at_r_arr(double*r, int Nr, double r_eff, double D_eff, +void xi_2h_at_r_arr(double*r, int Nr, double r_eff, double D_eff, double bias, double*ximm, double*xi2h){ int i; double*theta_eff = malloc(sizeof(double)*Nr); @@ -69,10 +67,9 @@ int xi_2h_at_r_arr(double*r, int Nr, double r_eff, double D_eff, xi2h[i] = (1-theta_eff[i]) * bias * ximm[i]; } free(theta_eff); - return 0; } -int xi_C_at_r_arr(double*r, int Nr, double r_A, double r_B, double D, +void xi_C_at_r_arr(double*r, int Nr, double r_A, double r_B, double D, double*xi_2h, double*xi_C){ int i; double*theta_A = malloc(sizeof(double)*Nr); @@ -84,7 +81,6 @@ int xi_C_at_r_arr(double*r, int Nr, double r_A, double r_B, double D, } free(theta_A); free(theta_B); - return 0; } @@ -97,9 +93,14 @@ int theta_erfc_at_r_arr(double*r, int Nr, double rt, double D, int i; double invD_rt = 1./(D*rt); for(i = 0; i < Nr; i++){ - theta[i] = 0.5*gsl_sf_erfc((r[i]-rt) * invD_rt * invsqrt2); + gsl_sf_result erf_result; + int rc = gsl_sf_erfc_e((r[i]-rt) * invD_rt * invsqrt2, + &erf_result); + if (rc != GSL_SUCCESS) + return rc; + theta[i] = 0.5*erf_result.val; } - return 0; + return GSL_SUCCESS; } /////////////////////////////////// diff --git a/src/C_massfunction.c b/src/C_massfunction.c index f4848f0..b372c54 100644 --- a/src/C_massfunction.c +++ b/src/C_massfunction.c @@ -20,7 +20,7 @@ #define del 1e-6 ///////////////// G multiplicity function/////////////////// -int G_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double d, double e, double f, double g, double*G){ +void G_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double d, double e, double f, double g, double*G){ double*sigma = (double*)malloc(sizeof(double)*NM); sigma2_at_M_arr(M, NM, k, P, Nk, om, sigma); int i; @@ -29,10 +29,9 @@ int G_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double d } G_at_sigma_arr(sigma, NM, d, e, f, g, G); free(sigma); - return 0; } -int G_at_sigma_arr(double*sigma, int Ns, double d, double e, double f, double g, double*G){ +void G_at_sigma_arr(double*sigma, int Ns, double d, double e, double f, double g, double*G){ //Compute the prefactor B double d2 = 0.5*d; double gamma_d2 = gsl_sf_gamma(d2); @@ -43,12 +42,11 @@ int G_at_sigma_arr(double*sigma, int Ns, double d, double e, double f, double g, for(i = 0; i < Ns; i++){ G[i] = B*exp(-g/(sigma[i]*sigma[i]))*(pow(sigma[i]/e, -d)+pow(sigma[i], -f)); } - return 0; } ///////////////// dndM functions below /////////////////// -int dndM_sigma2_precomputed(double*M, double*sigma2, double*dsigma2dM, int NM, double Omega_m, double d, double e, double f, double g, double*dndM){ +void dndM_sigma2_precomputed(double*M, double*sigma2, double*dsigma2dM, int NM, double Omega_m, double d, double e, double f, double g, double*dndM){ //This function exists to make emulator tests faster with sigma^2(M) precomputed. double rhom = Omega_m*rhocrit; //normalization coefficient double*sigma = (double*)malloc(sizeof(double)*NM); @@ -63,10 +61,9 @@ int dndM_sigma2_precomputed(double*M, double*sigma2, double*dsigma2dM, int NM, d } free(sigma); free(Gsigma); - return 0; } -int dndM_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double d, double e, double f, double g, double*dndM){ +void dndM_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double d, double e, double f, double g, double*dndM){ double*dsigma2dM = (double*)malloc(sizeof(double)*NM); double*sigma2 = (double*)malloc(sizeof(double)*NM); sigma2_at_M_arr(M, NM, k, P, Nk, om, sigma2); @@ -74,7 +71,6 @@ int dndM_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, doubl dndM_sigma2_precomputed(M, sigma2, dsigma2dM, NM, om, d, e, f, g, dndM); free(sigma2); free(dsigma2dM); - return 0; } ///////////////// derivatives of the MF below /////////////////// @@ -91,11 +87,18 @@ int n_in_bins(double*edges, int Nedges, double*M, double*dndM, int NM, double*N) gsl_spline*spline = gsl_spline_alloc(gsl_interp_cspline, NM); gsl_spline_init(spline, M, dndM, NM); gsl_interp_accel*acc = gsl_interp_accel_alloc(); - int i; + + if (!spline || !acc) + return GSL_ENOMEM; + + int i, rc = GSL_SUCCESS; for(i = 0; i < Nedges-1; i++){ - N[i] = gsl_spline_eval_integ(spline, edges[i], edges[i+1], acc); + rc = gsl_spline_eval_integ_e(spline, edges[i], edges[i+1], acc, &N[i]); + if (rc != GSL_SUCCESS) + break; } + gsl_spline_free(spline); gsl_interp_accel_free(acc); - return 0; + return rc; } diff --git a/src/C_miscentering.c b/src/C_miscentering.c index 9555ef7..233a93a 100644 --- a/src/C_miscentering.c +++ b/src/C_miscentering.c @@ -15,6 +15,7 @@ #include "gsl/gsl_integration.h" #include "gsl/gsl_spline.h" +#include "gsl/gsl_errno.h" #include #include @@ -108,18 +109,21 @@ int Sigma_mis_single_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, //Allocate things static int init_flag = 0; + // FIXME - the static variables are not thread safe static gsl_spline*spline = NULL; static gsl_interp_accel*acc = NULL; static gsl_integration_workspace*workspace = NULL; static integrand_params params; if (init_flag == 0){ - init_flag = 1; spline = gsl_spline_alloc(gsl_interp_cspline, Ns); acc = gsl_interp_accel_alloc(); workspace = gsl_integration_workspace_alloc(workspace_size); + if (!spline || !acc || !workspace) + return GSL_ENOMEM; + init_flag = 1; } - gsl_spline_init(spline, lnRs, Sigma, Ns); + int rc = gsl_spline_init(spline, lnRs, Sigma, Ns); params.acc = acc; params.spline = spline; @@ -140,16 +144,19 @@ int Sigma_mis_single_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, F.params = ¶ms; for(i = 0; i < NR; i++){ + if (rc != GSL_SUCCESS) + break; + params.Rp = R[i]; params.Rp2 = R[i] * R[i]; - gsl_integration_qag(&F, 0, M_PI, ABSERR, RELERR, workspace_size, + rc = gsl_integration_qag(&F, 0, M_PI, ABSERR, RELERR, workspace_size, KEY, workspace, &result, &err); Sigma_mis[i] = result/M_PI; } //Static objects aren't freed free(lnRs); - return 0; + return rc; } /////////////////// SIGMA(R) INTEGRANDS BELOW ////////////////////// @@ -273,14 +280,16 @@ int Sigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, static gsl_integration_workspace*workspace2 = NULL; static integrand_params params; if (init_flag == 0){ - init_flag = 1; spline = gsl_spline_alloc(gsl_interp_cspline, Ns); acc = gsl_interp_accel_alloc(); workspace = gsl_integration_workspace_alloc(workspace_size); workspace2 = gsl_integration_workspace_alloc(workspace_size); + if (!spline || !acc || !workspace || !workspace2) + return GSL_ENOMEM; + init_flag = 1; } - gsl_spline_init(spline, lnRs, Sigma, Ns); + int rc = gsl_spline_init(spline, lnRs, Sigma, Ns); params.spline = spline; params.acc = acc; @@ -296,7 +305,7 @@ int Sigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, params.rmax = Rs[Ns-1]; params.lrmin = log(Rs[0]); params.lrmax = log(Rs[Ns-1]); - + //Angular integral F.function = &angular_integrand; //Radial integral. Choice between rayliegh and gamma distributions @@ -309,24 +318,27 @@ int Sigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, F_radial.function = &Gamma_integrand; break; } - + //Assign the params struct to the GSL functions. F_radial.params = ¶ms; params.F_radial = F_radial; F.params = ¶ms; - + //Angular integral first for(i = 0; i < NR; i++){ + if (rc != GSL_SUCCESS) + break; + params.Rp = R[i]; params.Rp2 = R[i] * R[i]; //Optimization - - gsl_integration_qag(&F, 0, M_PI, ABSERR, RELERR, workspace_size, + + rc = gsl_integration_qag(&F, 0, M_PI, ABSERR, RELERR, workspace_size, KEY, workspace, &result, &err); Sigma_mis[i] = result/(M_PI*Rmis*Rmis); //Normalization } //Static objects aren't freed free(lnRs); - return 0; + return rc; } /////////////////////////////////////////////////// @@ -384,13 +396,15 @@ int DeltaSigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma_mis, int N static gsl_integration_workspace*workspace = NULL; static integrand_params params; if (init_flag == 0){ - init_flag = 1; spline = gsl_spline_alloc(gsl_interp_cspline, Ns); acc = gsl_interp_accel_alloc(); workspace = gsl_integration_workspace_alloc(workspace_size); + if (!spline || !acc || !workspace) + return GSL_ENOMEM; + init_flag = 1; } - gsl_spline_init(spline, Rs, Sigma_mis, Ns); + int rc = gsl_spline_init(spline, Rs, Sigma_mis, Ns); params.spline = spline; params.acc = acc; @@ -398,11 +412,13 @@ int DeltaSigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma_mis, int N F.function = &DS_mis_integrand; for(i = 0; i < NR; i++){ - gsl_integration_qag(&F, lrmin, log(R[i]), ABSERR, RELERR, workspace_size, + if (rc != GSL_SUCCESS) + break; + rc = gsl_integration_qag(&F, lrmin, log(R[i]), ABSERR, RELERR, workspace_size, KEY, workspace, &result, &err); DeltaSigma_mis[i] = (low_part+result)*2/(R[i]*R[i]) - gsl_spline_eval(spline, R[i], acc); } //No free() since we use static variables. - return 0; + return rc; } diff --git a/src/C_peak_height.c b/src/C_peak_height.c index 05e75e1..e96634b 100644 --- a/src/C_peak_height.c +++ b/src/C_peak_height.c @@ -1,6 +1,7 @@ #include "C_peak_height.h" #include "C_power.h" +#include "gsl/gsl_errno.h" #include "gsl/gsl_integration.h" #include "gsl/gsl_spline.h" @@ -86,14 +87,20 @@ int sigma2_at_R_arr(double*R, int NR, double*k, double*P, int Nk, double*s2){ gsl_spline*spline = gsl_spline_alloc(gsl_interp_cspline,Nk); gsl_interp_accel*acc = gsl_interp_accel_alloc(); gsl_integration_workspace*workspace = gsl_integration_workspace_alloc(workspace_size); + + // Handle allocation failure + if (!spline || !acc || !workspace) + return GSL_ENOMEM; + gsl_function F; integrand_params params; double lkmin = log(k[0]); double lkmax = log(k[Nk-1]); double result,abserr; double denom_inv = 1./(2*M_PI*M_PI); - int i; - gsl_spline_init(spline,k,P,Nk); + int i, rc; + + rc = gsl_spline_init(spline,k,P,Nk); params.spline = spline; params.acc = acc; params.kp = k; @@ -101,16 +108,19 @@ int sigma2_at_R_arr(double*R, int NR, double*k, double*P, int Nk, double*s2){ params.Nk = Nk; F.function = &sigma2_integrand; F.params = ¶ms; + for(i = 0; i < NR; i++){ + if (rc != GSL_SUCCESS) + break; params.r = R[i]; - gsl_integration_qag(&F, lkmin, lkmax, ABSERR, RELERR, + rc = gsl_integration_qag(&F, lkmin, lkmax, ABSERR, RELERR, workspace_size, KEY, workspace, &result, &abserr); s2[i] = result * denom_inv; //divide by 2pi^2 } gsl_spline_free(spline); gsl_interp_accel_free(acc); gsl_integration_workspace_free(workspace); - return 0; + return rc; } int sigma2_at_M_arr(double*M, int NM, double*k, double*P, int Nk, @@ -120,9 +130,9 @@ int sigma2_at_M_arr(double*M, int NM, double*k, double*P, int Nk, for(i = 0; i < NM; i++){ R[i] = M_to_R(M[i], Omega_m); } - sigma2_at_R_arr(R, NM, k, P, Nk, s2); + int rc = sigma2_at_R_arr(R, NM, k, P, Nk, s2); free(R); - return 0; + return rc; } /* The derivative with respect to R of sigma^2. This is needed for the mass @@ -134,6 +144,9 @@ int dsigma2dR_at_R_arr(double*R, int NR, double*k, double*P, int Nk, gsl_spline*spline = gsl_spline_alloc(gsl_interp_cspline,Nk); gsl_interp_accel*acc = gsl_interp_accel_alloc(); gsl_integration_workspace*workspace = gsl_integration_workspace_alloc(workspace_size); + if (!spline || !acc || !workspace) + return GSL_ENOMEM; + gsl_function F; integrand_params params; double lkmin = log(k[0]); @@ -143,8 +156,9 @@ int dsigma2dR_at_R_arr(double*R, int NR, double*k, double*P, int Nk, //out of the integrand because we take dw^2/dR //We also pull a factor of 3 out of the integrand. double denom_inv = 1./(M_PI*M_PI); - int i; - gsl_spline_init(spline,k,P,Nk); + int i, rc; + + rc = gsl_spline_init(spline,k,P,Nk); params.spline = spline; params.acc = acc; params.kp = k; @@ -153,15 +167,18 @@ int dsigma2dR_at_R_arr(double*R, int NR, double*k, double*P, int Nk, F.function = &dsigma2dR_integrand; F.params = ¶ms; for(i = 0; i < NR; i++){ + if (rc != GSL_SUCCESS) + break; params.r = R[i]; - gsl_integration_qag(&F, lkmin, lkmax, ABSERR, RELERR, + rc = gsl_integration_qag(&F, lkmin, lkmax, ABSERR, RELERR, workspace_size, KEY, workspace, &result, &abserr); ds2dR[i] = result * denom_inv; //divide by 2pi^2 } gsl_spline_free(spline); gsl_interp_accel_free(acc); gsl_integration_workspace_free(workspace); - return 0; + + return rc; } double dsigma2dR_at_R(double R, double*k, double*P, int Nk){ @@ -184,14 +201,14 @@ int dsigma2dM_at_M_arr(double*M, int NM, double*k, double*P, int Nk, dRdM[i] = dRdM_at_M(M[i], Omega_m); } //Note: ds2dM holds ds2dR for now - dsigma2dR_at_R_arr(R, NM, k, P, Nk, ds2dM); + int rc = dsigma2dR_at_R_arr(R, NM, k, P, Nk, ds2dM); //Chain rule to convert dsigma2dR to dsigma2dM for(i = 0; i < NM; i++){ ds2dM[i] = ds2dM[i] * dRdM[i]; } free(R); free(dRdM); - return 0; + return rc; } double dsigma2dM_at_M(double M, double*k, double*P, int Nk, double Omega_m){ @@ -217,12 +234,12 @@ int nu_at_R_arr(double*R, int NR, double*k, double*P, int Nk, double*nu){ //peak height at an array of R int i; double*s2 = (double*)malloc(sizeof(double)*NR); - sigma2_at_R_arr(R, NR, k, P, Nk, s2); + int rc = sigma2_at_R_arr(R, NR, k, P, Nk, s2); for(i = 0; i < NR; i++){ nu[i] = delta_c/sqrt(s2[i]); } free(s2); - return 0; + return rc; } int nu_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double*nu){ @@ -232,7 +249,7 @@ int nu_at_M_arr(double*M, int NM, double*k, double*P, int Nk, double om, double* for(i = 0; i < NM; i++){ R[i] = M_to_R(M[i], om); } - nu_at_R_arr(R, NM, k, P, Nk, nu); + int rc = nu_at_R_arr(R, NM, k, P, Nk, nu); free(R); - return 0; + return rc; } diff --git a/src/C_sigma_reconstruction.c b/src/C_sigma_reconstruction.c index 0d20004..3af6ad5 100644 --- a/src/C_sigma_reconstruction.c +++ b/src/C_sigma_reconstruction.c @@ -23,7 +23,7 @@ * the input DeltaSigma, also we assume a regular grid * spacing in the ln(R). */ -int Sigma_REC_from_DeltaSigma(double dlnR, double*DeltaSigma, int N, +void Sigma_REC_from_DeltaSigma(double dlnR, double*DeltaSigma, int N, double*Sigma){ /* The transformation is: @@ -46,10 +46,9 @@ int Sigma_REC_from_DeltaSigma(double dlnR, double*DeltaSigma, int N, //Downweight the first and last contributions (midpoint formula) //if ((j == N-2-i) || (j == N-1)){ if ((j == i) || (j == N-1)){ - Sigma[i] += dlnR*DeltaSigma[j]; + Sigma[i] += dlnR*DeltaSigma[j]; } //printf("%e\n",Sigma[i]); } } - return 0; } diff --git a/src/C_xi.c b/src/C_xi.c index 2f8e717..eba6094 100644 --- a/src/C_xi.c +++ b/src/C_xi.c @@ -32,7 +32,7 @@ double xi_nfw_at_r(double r, double Mass, double conc, int delta, double om){ return xi; } -int calc_xi_nfw(double*r, int Nr, double Mass, double conc, int delta, double om, double*xi_nfw){ +void calc_xi_nfw(double*r, int Nr, double Mass, double conc, int delta, double om, double*xi_nfw){ int i; double rhom = om*rhomconst;//SM h^2/Mpc^3 //double rho0_rhom = delta/(3.*(log(1.+conc)-conc/(1.+conc))); @@ -44,9 +44,7 @@ int calc_xi_nfw(double*r, int Nr, double Mass, double conc, int delta, double om r_rs = r[i]/rscale; //xi_nfw[i] = rho0_rhom/(r_rs*(1+r_rs)*(1+r_rs)) - 1.; xi_nfw[i] = Mass/(4.*M_PI*rscale*rscale*rscale*fc)/(r_rs*(1+r_rs)*(1+r_rs))/rhom - 1.0; - } - return 0; } double rhos_einasto_at_M(double Mass, double conc, double alpha, int delta, @@ -63,7 +61,7 @@ double rhos_einasto_at_M(double Mass, double conc, double alpha, int delta, return num/den; } -int calc_xi_einasto(double*r, int Nr, double Mass, double rhos, double conc, +void calc_xi_einasto(double*r, int Nr, double Mass, double rhos, double conc, double alpha, int delta, double Omega_m, double*xi_einasto){ double rhom = Omega_m*rhomconst;//SM h^2/Mpc^3 double rdelta = pow(Mass/(1.3333333333333*M_PI*rhom*delta), 0.333333333333); @@ -76,18 +74,16 @@ int calc_xi_einasto(double*r, int Nr, double Mass, double rhos, double conc, x = 2./alpha * pow(r[i]/rs, alpha); xi_einasto[i] = rhos/rhom * exp(-x) - 1; } - return 0; } -int calc_xi_2halo(int Nr, double bias, double*xi_mm, double*xi_2halo){ +void calc_xi_2halo(int Nr, double bias, double*xi_mm, double*xi_2halo){ int i; for(i = 0; i < Nr; i++){ xi_2halo[i] = bias * xi_mm[i]; } - return 0; } -int calc_xi_hm(int Nr, double*xi_1h, double*xi_2h, double*xi_hm, int flag){ +void calc_xi_hm(int Nr, double*xi_1h, double*xi_2h, double*xi_hm, int flag){ //Flag specifies how to combine the two terms int i; if (flag == 0) { //Take the max @@ -100,7 +96,6 @@ int calc_xi_hm(int Nr, double*xi_1h, double*xi_2h, double*xi_hm, int flag){ xi_hm[i] = 1 + xi_1h[i] + xi_2h[i]; } } - return 0; } int calc_xi_mm(double*r, int Nr, double*k, double*P, int Nk, double*xi, int N, double h){ @@ -109,6 +104,7 @@ int calc_xi_mm(double*r, int Nr, double*k, double*P, int Nk, double*xi, int N, d double PI_2 = M_PI*0.5; double sum; + // FIXME - static allocation is not thread-safe static int init_flag = 0; static gsl_spline*Pspl = NULL; static gsl_interp_accel*acc = NULL; @@ -122,10 +118,12 @@ int calc_xi_mm(double*r, int Nr, double*k, double*P, int Nk, double*xi, int N, d //Create the spline and accelerator if (init_flag == 0){ - Pspl = gsl_spline_alloc(gsl_interp_cspline, Nk); - acc = gsl_interp_accel_alloc(); + Pspl = gsl_spline_alloc(gsl_interp_cspline, Nk); + acc = gsl_interp_accel_alloc(); + if (!Pspl || !acc) + return GSL_ENOMEM; } - gsl_spline_init(Pspl, k, P, Nk); + int rc = gsl_spline_init(Pspl, k, P, Nk); //Compute things if ((init_flag == 0) || (h_old != h) || (N_old < N)){ @@ -163,7 +161,7 @@ int calc_xi_mm(double*r, int Nr, double*k, double*P, int Nk, double*xi, int N, d xi[j] = sum/(r[j]*r[j]*r[j]*M_PI*2); } - return 0; //Note: factor of pi picked up in the quadrature rule + return rc; //Note: factor of pi picked up in the quadrature rule //See Ogata 2005 for details, especially eq. 5.2 } @@ -207,14 +205,17 @@ int calc_xi_mm_exact(double*r, int Nr, double*k, double*P, int Nk, double*xi){ double kmin = 5e-8; double result, err; int i; - int status; gsl_spline*Pspl = gsl_spline_alloc(gsl_interp_cspline, Nk); gsl_interp_accel*acc= gsl_interp_accel_alloc(); - gsl_spline_init(Pspl, k, P, Nk); + if (!Pspl || !acc) + return GSL_ENOMEM; + int rc = gsl_spline_init(Pspl, k, P, Nk); gsl_integration_workspace*workspace = gsl_integration_workspace_alloc(workspace_size); gsl_integration_qawo_table*wf; + if (!Pspl || !acc) + return GSL_ENOMEM; integrand_params_xi_mm_exact params; params.acc = acc; @@ -228,18 +229,15 @@ int calc_xi_mm_exact(double*r, int Nr, double*k, double*P, int Nk, double*xi){ wf = gsl_integration_qawo_table_alloc(r[0], kmax-kmin, GSL_INTEG_SINE, (size_t)workspace_num); for(i = 0; i < Nr; i++){ - status = gsl_integration_qawo_table_set(wf, r[i], kmax-kmin, GSL_INTEG_SINE); - if (status){ - printf("Error in calc_xi_mm_exact, first integral.\n"); - exit(-1); - } + if (rc) + break; + rc = gsl_integration_qawo_table_set(wf, r[i], kmax-kmin, GSL_INTEG_SINE); + if (rc) + break; + params.r=r[i]; - status = gsl_integration_qawo(&F, kmin, ABSERR, RELERR, (size_t)workspace_num, + rc = gsl_integration_qawo(&F, kmin, ABSERR, RELERR, (size_t)workspace_num, workspace, wf, &result, &err); - if (status){ - printf("Error in calc_xi_mm_exact, second integral.\n"); - exit(-1); - } xi[i] = result/(M_PI*M_PI*2); } @@ -249,7 +247,7 @@ int calc_xi_mm_exact(double*r, int Nr, double*k, double*P, int Nk, double*xi){ gsl_integration_workspace_free(workspace); gsl_integration_qawo_table_free(wf); - return 0; + return rc; } double xi_mm_at_r_exact(double r, double*k, double*P, int Nk){ @@ -262,7 +260,7 @@ double xi_mm_at_r_exact(double r, double*k, double*P, int Nk){ * Diemer-Kravtsov 2014 profiles below. */ -int calc_xi_DK(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double*xi){ +void calc_xi_DK(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double*xi){ double rhom = rhomconst*om; //SM h^2/Mpc^3 //Compute rDeltam double rdelta = pow(M/(1.33333333333*M_PI*rhom*delta), 0.33333333333); @@ -297,13 +295,12 @@ int calc_xi_DK(double*r, int Nr, double M, double rhos, double conc, double be, free(rho_ein); free(f_trans); free(rho_outer); - return 0; } ////////////////////////////// //////Appendix version 1////// ////////////////////////////// -int calc_xi_DK_app1(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi){ +void calc_xi_DK_app1(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi){ double rhom = rhomconst*om; //SM h^2/Mpc^3 //Compute r200m double rdelta = pow(M/(1.33333333333*M_PI*rhom*delta), 0.33333333333); @@ -339,13 +336,12 @@ int calc_xi_DK_app1(double*r, int Nr, double M, double rhos, double conc, double free(rho_ein); free(f_trans); free(rho_outer); - return 0; } ////////////////////////////// //////Appendix version 2////// ////////////////////////////// -int calc_xi_DK_app2(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi){ +void calc_xi_DK_app2(double*r, int Nr, double M, double rhos, double conc, double be, double se, double alpha, double beta, double gamma, int delta, double*k, double*P, int Nk, double om, double bias, double*xi_mm, double*xi){ double rhom = rhomconst*om; //SM h^2/Mpc^3 //Compute r200m double rdelta = pow(M/(1.33333333333*M_PI*rhom*delta), 0.33333333333); @@ -381,5 +377,4 @@ int calc_xi_DK_app2(double*r, int Nr, double M, double rhos, double conc, double free(rho_ein); free(f_trans); free(rho_outer); - return 0; } diff --git a/tests/test_bias.py b/tests/test_bias.py index 2eeaff0..1009a45 100644 --- a/tests/test_bias.py +++ b/tests/test_bias.py @@ -44,7 +44,6 @@ def test_derivatives(): deriv = db/dM pd = dbdM[:-1] / deriv npt.assert_array_almost_equal(pd, np.ones_like(pd), 1e-2) - return def test_s2_and_nu_functions(): #Test the mass calls