diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index df2c983a..ee259b24 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -485,5 +485,6 @@ from . import input from . import linalg from . import mpm from . import roots +from . import integrate from .version import __version__ diff --git a/pyerrors/integrate.py b/pyerrors/integrate.py new file mode 100644 index 00000000..3b48f3fb --- /dev/null +++ b/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:]) diff --git a/tests/integrate_test.py b/tests/integrate_test.py new file mode 100644 index 00000000..07a4313b --- /dev/null +++ b/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 diff --git a/tests/linalg_test.py b/tests/linalg_test.py index babf3797..ec0591ba 100644 --- a/tests/linalg_test.py +++ b/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 test_einsum(): 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) +