pyerrors.integrate

 1import numpy as np
 2from autograd import jacobian
 3from scipy.integrate import quad as squad
 4
 5from .obs import Obs, derived_observable
 6
 7
 8def quad(func, p, a, b, **kwargs):
 9    '''Performs a (one-dimensional) numeric integration of f(p, x) from a to b.
10
11    The integration is performed using scipy.integrate.quad().
12    All parameters that can be passed to scipy.integrate.quad may also be passed to this function.
13    The output is the same as for scipy.integrate.quad, the first element being an Obs.
14
15    Parameters
16    ----------
17    func : object
18        function to integrate, has to be of the form
19
20        ```python
21        import autograd.numpy as anp
22
23        def func(p, x):
24            return p[0] + p[1] * x + p[2] * anp.sinh(x)
25        ```
26        where x is the integration variable.
27    p : list of floats or Obs
28        parameters of the function func.
29    a: float or Obs
30        Lower limit of integration (use -numpy.inf for -infinity).
31    b: float or Obs
32        Upper limit of integration (use -numpy.inf for -infinity).
33    All parameters of scipy.integrate.quad
34
35    Returns
36    -------
37    y : Obs
38        The integral of func from `a` to `b`.
39    abserr : float
40        An estimate of the absolute error in the result.
41    infodict : dict
42        A dictionary containing additional information.
43        Run scipy.integrate.quad_explain() for more information.
44    message
45        A convergence message.
46    explain
47        Appended only with 'cos' or 'sin' weighting and infinite
48        integration limits, it contains an explanation of the codes in
49        infodict['ierlst']
50    '''
51
52    Np = len(p)
53    isobs = [True if isinstance(pi, Obs) else False for pi in p]
54    pval = np.array([p[i].value if isobs[i] else p[i] for i in range(Np)],)
55    pobs = [p[i] for i in range(Np) if isobs[i]]
56
57    bounds = [a, b]
58    isobs_b = [True if isinstance(bi, Obs) else False for bi in bounds]
59    bval = np.array([bounds[i].value if isobs_b[i] else bounds[i] for i in range(2)])
60    bobs = [bounds[i] for i in range(2) if isobs_b[i]]
61    bsign = [-1, 1]
62
63    ifunc = np.vectorize(lambda x: func(pval, x))
64
65    intpars = squad.__code__.co_varnames[3:3 + len(squad.__defaults__)]
66    ikwargs = {k: kwargs[k] for k in intpars if k in kwargs}
67
68    integration_result = squad(ifunc, bval[0], bval[1], **ikwargs)
69    val = integration_result[0]
70
71    jac = jacobian(func)
72
73    derivint = []
74    for i in range(Np):
75        if isobs[i]:
76            ifunc = np.vectorize(lambda x, i=i: jac(pval, x)[i])
77            derivint.append(squad(ifunc, bounds[0], bounds[1], **ikwargs)[0])
78
79    for i in range(2):
80        if isobs_b[i]:
81            derivint.append(bsign[i] * func(pval, bval[i]))
82
83    if len(derivint) == 0:
84        return integration_result
85
86    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)
87
88    return (res, *integration_result[1:])