#288 Expose c++ integrator to python and remove scipy dependence

INSTALL.md

@@ -63,17 +63,12 @@ The GalSim package also requires
   functionality, and can be omitted if users are happy to use JSON-style
   config parsing or prefer to write python scripts directly.
 
-* Optional dependency: the scientific Python module SciPy
-  (http://www.scipy.org).  SciPy is only required for generation of
-  shears due to an NFW halo.
-
 These should installed onto your Python system so that they can be imported by:
 
     >>> import numpy
     >>> import pyfits
     >>> import yaml [ if using the galsim_yaml executable or otherwise plan to
                       parse .yaml configuration files ]
-    >>> import scipy [ if generating shears according to an NFW halo ]
 
 within Python.  You can test this by loading up the Python interpreter for the
 version of Python you will be using with the GalSim toolkit. This is usually

galsim/__init__.py

@@ -16,3 +16,4 @@
 from psfcorr import *
 from io import *
 from . import lensing
+from . import integ

galsim/integ.py

@@ -0,0 +1,36 @@
+"""@file integ.py
+A Python layer version of the C++ int1d function in galim::integ
+"""
+
+from . import _galsim
+
+def int1d(func, min, max, rel_err=1.e-6, abs_err=1.e-12):
+    """Integrate a 1-dimensional function from min to max.
+
+    Example usage:
+
+        >>> def func(x): return x**2
+        >>> galsim.integ.int1d(func, 0, 1)
+        0.33333333333333337
+        >>> galsim.integ.int1d(func, 0, 2)
+        2.666666666666667
+        >>> galsim.integ.int1d(func, -1, 1)
+        0.66666666666666674
+
+    @param func     The function to be integrated.  y = func(x) should be valid.
+    @param min      The lower end of the integration bounds (anything < -1.e10 is treated as
+                    negative infinity).
+    @param max      The upper end of the integration bounds (anything > 1.e10 is treated as positive
+                    infinity).
+    @param rel_err  The desired relative error (default `rel_err = 1.e-6`)
+    @param abs_err  The desired absolute error (default `abs_err = 1.e-12`)
+    """
+    min = float(min)
+    max = float(max)
+    rel_err = float(rel_err)
+    abs_err = float(abs_err)
+    success, result = _galsim.PyInt1d(func,min,max,rel_err,abs_err)
+    if success:
+        return result
+    else:
+        raise RuntimeError(result)

galsim/lensing.py

@@ -414,9 +414,8 @@ def Da(self, z, z_ref=0):
                 raise ValueError("Redshift z must not be negative")
             if z < z_ref:
                 raise ValueError("Redshift z must not be smaller than the reference redshift")
-            # TODO: We don't want to depend on scipy, so need to move this down to c++.
-            from scipy.integrate import quad
-            d = quad(self.__angKernel, z_ref+1, z+1)[0]
+
+            d = galsim.integ.int1d(self.__angKernel, z_ref+1, z+1)
             # check for curvature
             rk = (abs(self.omega_c))**0.5
             if (rk*d > 0.01):