Two parallelized functions (#727)

* parallelized integ_spline

* missing omp for

* parallelized fftlog

* visible from python

include/ccl_utils.h

@@ -49,16 +49,19 @@ double ccl_j_bessel(int l,double x);
 
 /**
  * Compute spline integral.
- * @param n number of elements in input array.
+ * @param nx number of elements in input array.
+ * @param ny number of y arrays.
  * @param x input x-values.
- * @param y input y-values.
+ * @param y input y-values (ny arrays with nx elements).
  * @param a lower end of integration range.
  * @param b upper end of integration range (use b<a if you just want to integrate over all of y).
- * @return spline integral
+ * @param result array of output spline integral values.
+ * @param T spline type.
+ * @param status status flag.
  */
-double ccl_integ_spline(int n,double *x,double *y,
-                        double a, double b,
-                        const gsl_interp_type *T, int *status);
+void ccl_integ_spline(int ny, int nx,double *x,double **y,
+                      double a, double b, double *result,
+                      const gsl_interp_type *T, int *status);
 
 CCL_END_DECLS
 

pyccl/ccl.i

@@ -54,8 +54,8 @@ from .errors import CCLError
 %include "ccl_sigM.i"
 %include "ccl_f1d.i"
 %include "ccl_fftlog.i"
+%include "ccl_utils.i"
 
 /* list header files not yet having a .i file here */
 %include "../include/ccl_config.h"
 %include "../include/ccl_error.h"
-%include "../include/ccl_utils.h"

pyccl/ccl_utils.i

@@ -0,0 +1,52 @@
+%module ccl_utils
+
+%{
+/* put additional #include here */
+%}
+
+%include "../include/ccl_utils.h"
+
+%apply (double* IN_ARRAY1, int DIM1) {
+  (double *x_in, int n_in_x),
+  (double *ys_in, int n_in_y)};
+%apply (int DIM1, double* ARGOUT_ARRAY1) {(int nout, double* output)};
+
+
+%feature("pythonprepend") spline_integrate %{
+    if n_integ * x_in.size != ys_in.size:
+        raise CCLError("Input size for `ys_in` must match `n_integ * x_in.size`")
+
+    if nout != n_integ:
+        raise CCLError("Input shape for `output` must match n_integ")
+%}
+
+%inline %{
+
+
+void spline_integrate(int n_integ,
+                      double *x_in, int n_in_x,
+                      double *ys_in, int n_in_y,
+                      double a, double b,
+                      int nout, double *output,
+                      int *status)
+{
+  int ii;
+
+  double **_ys_in=NULL;
+  _ys_in = malloc(n_integ*sizeof(double *));
+  if(_ys_in==NULL)
+    *status = CCL_ERROR_MEMORY;
+
+  if(*status==0) {
+    for(ii=0;ii<n_integ;ii++)
+      _ys_in[ii]=&(ys_in[ii*n_in_x]);
+
+    ccl_integ_spline(n_integ, n_in_x, x_in, _ys_in,
+                     a, b, output, gsl_interp_akima,
+                     status);
+  }
+
+  free(_ys_in);
+}
+
+%}

pyccl/pyutils.py

@@ -432,3 +432,33 @@ def _fftlog_transform(rs, frs,
         fks = fks.squeeze()
 
     return ks, fks
+
+
+def _spline_integrate(x, ys, a, b):
+    if np.ndim(x) != 1:
+        raise ValueError("x should be a 1D array")
+    if np.ndim(ys) < 1 or np.ndim(ys) > 2:
+        raise ValueError("ys should be 1D or a 2D array")
+    if np.ndim(ys) == 1:
+        n_integ = 1
+        n_x = len(ys)
+    else:
+        n_integ, n_x = ys.shape
+
+    if len(x) != n_x:
+        raise ValueError("x should have %d elements" % n_x)
+
+    if np.ndim(a) > 0 or np.ndim(b) > 0:
+        raise TypeError("Integration limits should be scalar")
+
+    status = 0
+    result, status = lib.spline_integrate(n_integ,
+                                          x, ys.flatten(),
+                                          a, b, n_integ,
+                                          status)
+    check(status)
+
+    if np.ndim(ys) == 1:
+        result = result[0]
+
+    return result

pyccl/tests/test_spline_integ.py

@@ -0,0 +1,67 @@
+import numpy as np
+import pytest
+from pyccl.pyutils import _spline_integrate
+from pyccl.errors import CCLError
+
+
+def yf(x, pw):
+    return x**pw
+
+
+def yf_int(x, pw):
+    return x**(pw + 1) / (pw + 1)
+
+
+def test_intspline_accuracy():
+    nx = 1024
+    x_arr = np.linspace(-2, 2, nx)
+    y_arr = np.array([yf(x_arr, p)
+                      for p in [2, 4]])
+    int_spline = _spline_integrate(x_arr, y_arr, -1, 1)
+    int_exact = np.array([yf_int(1, p) - yf_int(-1, p)
+                          for p in [2, 4]])
+    assert np.all(np.fabs(int_spline / int_exact - 1) < 1E5-5)
+
+
+def test_intspline_raises_shapes():
+    x = np.linspace(0, 1, 10)
+    y = np.zeros([4, 10])
+
+    # x is 0- or 2-D
+    with pytest.raises(ValueError):
+        _spline_integrate(x[0], y, 0, 1)
+
+    with pytest.raises(ValueError):
+        _spline_integrate(y, y, 0, 1)
+
+    # y is 0-D
+    with pytest.raises(ValueError):
+        _spline_integrate(x, y[0, 0], 0, 1)
+
+    # Mismatching x and y
+    with pytest.raises(ValueError):
+        _spline_integrate(x[1:], y, 0, 1)
+
+    # Limits as arrays
+    with pytest.raises(TypeError):
+        _spline_integrate(x, y, np.zeros(10), 1)
+
+
+def test_intspline_raises_internal():
+    x = np.linspace(0, 1, 10)
+    y = np.zeros([4, 10])
+
+    # Function is not splineable
+    with pytest.raises(CCLError):
+        _spline_integrate(0 * x, y, 0, 1)
+
+
+def test_intspline_smoke():
+    x = np.linspace(0, 1, 10)
+    y = yf(x, 2)
+
+    r = _spline_integrate(x, y, 0, 1)
+    assert np.all(~np.isnan(r))
+
+    r = _spline_integrate(x, np.array([y, y]), 0, 1)
+    assert np.all(~np.isnan(r))

src/ccl_cls.c

@@ -167,9 +167,9 @@ static void integ_cls_limber_spline(ccl_cosmology *cosmo,
   }
 
   if(*status == 0) {
-    *result = ccl_integ_spline(nk, lk_arr, fk_arr,
-			       1, -1, gsl_interp_akima,
-			       status);
+    ccl_integ_spline(1, nk, lk_arr, &fk_arr,
+                     1, -1, result, gsl_interp_akima,
+                     status);
   }
   free(fk_arr);
   free(lk_arr);