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>
This commit is contained in:
s-kuberski 2023-07-14 15:21:59 +02:00 committed by GitHub
parent 8736d1cd3c
commit 6dcd0c3518
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 141 additions and 2 deletions

View file

@ -485,5 +485,6 @@ from . import input
from . import linalg from . import linalg
from . import mpm from . import mpm
from . import roots from . import roots
from . import integrate
from .version import __version__ from .version import __version__

87
pyerrors/integrate.py Normal file
View file

@ -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:])

51
tests/integrate_test.py Normal file
View file

@ -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

View file

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