mirror of
https://github.com/fjosw/pyerrors.git
synced 2026-08-04 11:31:21 +02:00
72 lines
2.3 KiB
Python
72 lines
2.3 KiB
Python
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
|
|
|
|
|
|
def test_integrate_per_parameter_derivatives():
|
|
# \int_0^1 p0*x + p1*x^2 dx = p0/2 + p1/3
|
|
# If the lambda closure in integrate.quad failed to bind `i` per-iteration,
|
|
# all per-parameter derivatives would collapse to a single value.
|
|
def f(p, x):
|
|
return p[0] * x + p[1] * x ** 2
|
|
|
|
p = [pe.cov_Obs(1.0, 0.1 ** 2, "p0"), pe.cov_Obs(2.0, 0.2 ** 2, "p1")]
|
|
res, _ = pe.integrate.quad(f, p, 0.0, 1.0)
|
|
res.gm()
|
|
|
|
ana_val = p[0].value / 2 + p[1].value / 3
|
|
assert np.isclose(res.value, ana_val)
|
|
|
|
grad0 = res.covobs["p0"].grad.item()
|
|
grad1 = res.covobs["p1"].grad.item()
|
|
assert np.isclose(grad0, 0.5)
|
|
assert np.isclose(grad1, 1.0 / 3.0)
|
|
assert not np.isclose(grad0, grad1)
|