feat: Added numerical integration of generic functions (#201)

* feat: Added numerical integration of generic functions

* refactored integration routines

* tests: two trivial tests for integration added.

* docs: quad docstring corrected.

* Small bugfix for integration without obs

---------

Co-authored-by: Fabian Joswig <fabian.joswig@ed.ac.uk>

pyerrors/__init__.py

@@ -485,5 +485,6 @@ def func(a, x):
 from . import linalg
 from . import mpm
 from . import roots
+from . import integrate
 
 from .version import __version__

pyerrors/integrate.py

@@ -0,0 +1,87 @@
+import numpy as np
+from .obs import derived_observable, Obs
+from autograd import jacobian
+from scipy.integrate import quad as squad
+
+
+def quad(func, p, a, b, **kwargs):
+    '''Performs a (one-dimensional) numeric integration of f(p, x) from a to b.
+
+    The integration is performed using scipy.integrate.quad().
+    All parameters that can be passed to scipy.integrate.quad may also be passed to this function.
+    The output is the same as for scipy.integrate.quad, the first element being an Obs.
+
+    Parameters
+    ----------
+    func : object
+        function to integrate, has to be of the form
+
+        ```python
+        import autograd.numpy as anp
+
+        def func(p, x):
+            return p[0] + p[1] * x + p[2] * anp.sinh(x)
+        ```
+        where x is the integration variable.
+    p : list of floats or Obs
+        parameters of the function func.
+    a: float or Obs
+        Lower limit of integration (use -numpy.inf for -infinity).
+    b: float or Obs
+        Upper limit of integration (use -numpy.inf for -infinity).
+    All parameters of scipy.integrate.quad
+
+    Returns
+    -------
+    y : Obs
+        The integral of func from `a` to `b`.
+    abserr : float
+        An estimate of the absolute error in the result.
+    infodict : dict
+        A dictionary containing additional information.
+        Run scipy.integrate.quad_explain() for more information.
+    message
+        A convergence message.
+    explain
+        Appended only with 'cos' or 'sin' weighting and infinite
+        integration limits, it contains an explanation of the codes in
+        infodict['ierlst']
+    '''
+
+    Np = len(p)
+    isobs = [True if isinstance(pi, Obs) else False for pi in p]
+    pval = np.array([p[i].value if isobs[i] else p[i] for i in range(Np)],)
+    pobs = [p[i] for i in range(Np) if isobs[i]]
+
+    bounds = [a, b]
+    isobs_b = [True if isinstance(bi, Obs) else False for bi in bounds]
+    bval = np.array([bounds[i].value if isobs_b[i] else bounds[i] for i in range(2)])
+    bobs = [bounds[i] for i in range(2) if isobs_b[i]]
+    bsign = [-1, 1]
+
+    ifunc = np.vectorize(lambda x: func(pval, x))
+
+    intpars = squad.__code__.co_varnames[3:3 + len(squad.__defaults__)]
+    ikwargs = {k: kwargs[k] for k in intpars if k in kwargs}
+
+    integration_result = squad(ifunc, bval[0], bval[1], **ikwargs)
+    val = integration_result[0]
+
+    jac = jacobian(func)
+
+    derivint = []
+    for i in range(Np):
+        if isobs[i]:
+            ifunc = np.vectorize(lambda x: jac(pval, x)[i])
+            derivint.append(squad(ifunc, bounds[0], bounds[1], **ikwargs)[0])
+
+    for i in range(2):
+        if isobs_b[i]:
+            derivint.append(bsign[i] * func(pval, bval[i]))
+
+    if len(derivint) == 0:
+        return integration_result
+
+    res = derived_observable(lambda x, **kwargs: 0 * (x[0] + np.finfo(np.float64).eps) * (pval[0] + np.finfo(np.float64).eps) + val, pobs + bobs, man_grad=derivint)
+
+    return (res, *integration_result[1:])

tests/integrate_test.py

@@ -0,0 +1,51 @@
+import numpy as np
+import autograd.numpy as anp
+import pyerrors as pe
+
+
+def test_integration():
+    def f(p, x):
+        return p[0] * x + p[1] * x**2 - p[2] / x
+
+    def F(p, x):
+        return p[0] * x**2 / 2. + p[1] * x**3 / 3. - anp.log(x) * p[2]
+
+    def check_ana_vs_int(p, l, u, **kwargs):
+        numint_full = pe.integrate.quad(f, p, l, u, **kwargs)
+        numint = numint_full[0]
+
+        anaint = F(p, u) - F(p, l)
+        diff = (numint - anaint)
+
+        if isinstance(numint, pe.Obs):
+            numint.gm()
+            anaint.gm()
+
+            assert(diff.is_zero())
+        else:
+            assert(np.isclose(0, diff))
+
+    pobs = np.array([pe.cov_Obs(1., .1**2, '0'), pe.cov_Obs(2., .2**2, '1'), pe.cov_Obs(2.2, .17**2, '2')])
+    lobs = pe.cov_Obs(.123, .012**2, 'l')
+    uobs = pe.cov_Obs(1., .05**2, 'u')
+
+    check_ana_vs_int(pobs, lobs, uobs)
+    check_ana_vs_int(pobs, lobs.value, uobs)
+    check_ana_vs_int(pobs, lobs, uobs.value)
+    check_ana_vs_int(pobs, lobs.value, uobs.value)
+    for i in range(len(pobs)):
+        p = [pi for pi in pobs]
+        p[i] = pobs[i].value
+        check_ana_vs_int(p, lobs, uobs)
+
+    check_ana_vs_int([pi.value for pi in pobs], lobs, uobs)
+    check_ana_vs_int([pi.value for pi in pobs], lobs.value, uobs.value)
+
+    check_ana_vs_int(pobs, lobs, uobs, epsabs=1.e-9, epsrel=1.236e-10, limit=100)
+    assert(len(pe.integrate.quad(f, pobs, lobs, uobs, full_output=True)) > 2)
+
+    r1, _ = pe.integrate.quad(F, pobs, 1, 0.1)
+    r2, _ = pe.integrate.quad(F, pobs, 0.1, 1)
+    assert r1 == -r2
+    iamzero, _ = pe.integrate.quad(F, pobs, 1, 1)
+    assert iamzero == 0

tests/linalg_test.py

@@ -14,9 +14,9 @@ def get_real_matrix(dimension):
         exponent_imag = np.random.normal(0, 1)
         base_matrix[n, m] = pe.Obs([np.random.normal(1.0, 0.1, 100)], ['t'])
 
-
     return base_matrix
 
+
 def get_complex_matrix(dimension):
     base_matrix = np.empty((dimension, dimension), dtype=object)
     for (n, m), entry in np.ndenumerate(base_matrix):
@@ -109,7 +109,6 @@ def _perform_complex_check(arr):
         assert np.all([o.imag.is_zero_within_error(0.001) for o in arr])
         assert np.all([o.imag.dvalue < 0.001 for o in arr])
 
-
     tt = [get_real_matrix(4), get_real_matrix(3)]
     q = np.tensordot(tt[0], tt[1], 0)
     c1 = tt[1] @ q
@@ -355,3 +354,4 @@ def test_complex_matrix_real_entries():
     my_mat[0, 1] = 4
     my_mat[2, 0] = pe.Obs([np.random.normal(1.0, 0.1, 100)], ['t'])
     assert np.all((my_mat @ pe.linalg.inv(my_mat) - np.identity(4)) == 0)
+