Consistent ffi

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

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:

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()

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()

cluster_toolkit/concentration.py

@@ -2,15 +2,15 @@
 
 """
 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"):
     """Concentration of the NFW profile at mass M [Msun/h].
     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)
-
-