pyerrors.integrate

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