Allow access to gsl spline data (#893)

* Allow access to gsl spline data.

* Added docstrings.

* use memcpy

pyccl/ccl_utils.i

@@ -11,6 +11,13 @@
   (double *ys_in, int n_in_y)};
 %apply (int DIM1, double* ARGOUT_ARRAY1) {(int nout, double* output)};
 
+%apply (int DIM1, double* ARGOUT_ARRAY1) {(int x_size, double* xarr)};
+%apply (int DIM1, double* ARGOUT_ARRAY1) {(int y_size, double* yarr)};
+%apply (int DIM1, double* ARGOUT_ARRAY1) {(int z_size, double* zarr)};
+
+%apply (int* OUTPUT) {(int *size)};
+%apply int *OUTPUT { int *x_size, int *y_size };
+
 
 %feature("pythonprepend") spline_integrate %{
     if n_integ * x_in.size != ys_in.size:
@@ -49,4 +56,70 @@ void spline_integrate(int n_integ,
   free(_ys_in);
 }
 
+void get_spline1d_array_size(gsl_spline *spline, int* size, int* status) {
+  if(spline == NULL) {
+    *status = CCL_ERROR_MEMORY;
+    return;
+  }
+  if(spline->interp == NULL) {
+    *status = CCL_ERROR_MEMORY;
+    return;
+  }
+  *size = spline->interp->size;
+}
+
+void get_spline2d_array_sizes(gsl_spline2d *spline2d, int* x_size, int* y_size,
+                              int* status) {
+  if(spline2d == NULL) {
+    *status = CCL_ERROR_MEMORY;
+    return;
+  }
+  *x_size = spline2d->interp_object.xsize;
+  *y_size = spline2d->interp_object.ysize;
+}
+
+void get_spline1d_arrays(gsl_spline *spline,
+                         int x_size, double* xarr,
+                         int y_size, double* yarr,
+                         int *status)
+{
+  if(spline == NULL) {
+    *status = CCL_ERROR_MEMORY;
+    return;
+  }
+  if(x_size != spline->interp->size) {
+    *status = CCL_ERROR_INCONSISTENT;
+    return;
+  }
+  if(y_size != spline->interp->size) {
+    *status = CCL_ERROR_INCONSISTENT;
+    return;
+  }
+  memcpy(xarr, spline->x, sizeof(double)*x_size);
+  memcpy(yarr, spline->y, sizeof(double)*y_size);
+}
+
+void get_spline2d_arrays(gsl_spline2d *spline2d,
+                         int x_size, double* xarr,
+                         int y_size, double* yarr,
+                         int z_size, double* zarr,
+                         int *status)
+{
+  if(spline2d == NULL) {
+    *status = CCL_ERROR_MEMORY;
+    return;
+  }
+  if(x_size != spline2d->interp_object.xsize) {
+    *status = CCL_ERROR_INCONSISTENT;
+    return;
+  }
+  if(y_size != spline2d->interp_object.ysize) {
+    *status = CCL_ERROR_INCONSISTENT;
+    return;
+  }
+  memcpy(xarr, spline2d->xarr, sizeof(double)*x_size);
+  memcpy(yarr, spline2d->yarr, sizeof(double)*y_size);
+  memcpy(zarr, spline2d->zarr, sizeof(double)*x_size*y_size);
+}
+
 %}

pyccl/pyutils.py

@@ -518,3 +518,55 @@ def assert_warns(wtype, f, *args, **kwargs):
         "Warning raised was the wrong type (got %s, expected %s)" % (
             w[0].category, wtype)
     return res
+
+
+def _get_spline1d_arrays(gsl_spline):
+    """Get array data from a 1D GSL spline.
+
+    Args:
+        gsl_spline: `SWIGObject` of gsl_spline
+            The SWIG object of the GSL spline.
+
+    Returns:
+        xarr: array_like
+            The x array of the spline.
+        yarr: array_like
+            The y array of the spline.
+    """
+    status = 0
+    size, status = lib.get_spline1d_array_size(gsl_spline, status)
+    check(status)
+
+    xarr, yarr, status = lib.get_spline1d_arrays(gsl_spline, size, size,
+                                                 status)
+    check(status)
+
+    return xarr, yarr
+
+
+def _get_spline2d_arrays(gsl_spline):
+    """Get array data from a 2D GSL spline.
+
+    Args:
+        gsl_spline: `SWIGObject` of gsl_spline2d
+            The SWIG object of the 2D GSL spline.
+
+    Returns:
+        yarr: array_like
+            The y array of the spline.
+        xarr: array_like
+            The x array of the spline.
+        zarr: array_like
+            The z array of the spline. The shape is (yarr.size, xarr.size).
+    """
+    status = 0
+    x_size, y_size, status = lib.get_spline2d_array_sizes(gsl_spline, status)
+    check(status)
+
+    z_size = x_size*y_size
+    xarr, yarr, zarr, status = lib.get_spline2d_arrays(gsl_spline,
+                                                       x_size, y_size, z_size,
+                                                       status)
+    check(status)
+
+    return yarr, xarr, zarr.reshape(y_size, x_size)

pyccl/tests/test_gsl_splines.py

@@ -0,0 +1,47 @@
+import pytest
+import pyccl as ccl
+
+import numpy as np
+
+
+def test_spline_params():
+    cosmo = ccl.Cosmology(
+        Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.8, n_s=0.96,
+        transfer_function='bbks', matter_power_spectrum='linear')
+
+    assert cosmo.cosmo.spline_params.A_SPLINE_MAX == 1.0
+
+    with pytest.raises(RuntimeError):
+        cosmo.cosmo.spline_params.A_SPLINE_MAX = 0.9
+
+
+def test_spline1d():
+    cosmo = ccl.CosmologyVanillaLCDM()
+    cosmo.compute_distances()
+
+    chi_gsl_spline = cosmo.cosmo.data.chi
+    a_arr, chi_arr = ccl.pyutils._get_spline1d_arrays(chi_gsl_spline)
+    chi = ccl.comoving_radial_distance(cosmo, a_arr)
+
+    assert np.allclose(chi_arr, chi)
+
+
+def test_spline2d():
+    x = np.linspace(0.1, 1, 10)
+    log_y = np.linspace(-3, 1, 20)
+    zarr_in = np.outer(x, np.exp(log_y))
+
+    pk2d = ccl.Pk2D(a_arr=x, lk_arr=log_y, pk_arr=zarr_in,
+                    is_logp=False)
+
+    pk2d_gsl_spline2d = pk2d.psp.fka
+    xarr, yarr, zarr_out_spline = \
+        ccl.pyutils._get_spline2d_arrays(pk2d_gsl_spline2d)
+
+    cosmo = ccl.CosmologyVanillaLCDM()
+    zarr_out_eval = pk2d.eval(k=np.exp(log_y), a=x[-1], cosmo=cosmo)
+
+    assert np.allclose(x, xarr)
+    assert np.allclose(log_y, yarr)
+    assert np.allclose(zarr_in, zarr_out_spline)
+    assert np.allclose(zarr_in[-1], zarr_out_eval)