Skip to content

Commit

Permalink
Fix tests on FreeBSD, remove util._jacobi (#903)
Browse files Browse the repository at this point in the history 
Closes #901

The private function `iminuit.util._jacobi` is removed without
deprecation. It was internally replaced by
`scipy.optimize.approx_fprime`. Those who used `_jacobi` for some reason
should switch to `jacobi.jacobi` from the `jacobi` package if a
high-accuracy derivative is required or use
`scipy.optimize_approx_fprime` for a low-accuracy derivative.
  • Loading branch information
@HDembinski
HDembinski committed Jun 20, 2023
1 parent 3cbdf06 commit 70036f7
Show file tree
Hide file tree
Showing 3 changed files with 25 additions and 159 deletions.
6 changes: 3 additions & 3 deletions src/iminuit/minuit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -996,6 +996,7 @@ def scipy(
Bounds,
NonlinearConstraint,
LinearConstraint,
approx_fprime,
)
except ModuleNotFoundError as exc:
exc.msg += "\n\nPlease install scipy to use scipy minimizers in iminuit."
Expand Down Expand Up @@ -1271,9 +1272,8 @@ def __call__(self, par, v):
elif "jac" in r:
jac = r.jac
else:
tol = 1e-2
dx = np.sqrt(np.diag(matrix) * tol)
jac = mutil._jacobi(fcn, r.x, dx, tol)[1][0]
dx = np.sqrt(np.diag(matrix) * 1e-10)
jac = approx_fprime(r.x, fcn, epsilon=dx)

edm_goal = self._edm_goal(migrad_factor=True)
fm = FunctionMinimum(
Expand Down
61 changes: 5 additions & 56 deletions src/iminuit/util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -949,60 +949,6 @@ def __getitem__(self, key):
return OrderedDict.__getitem__(self, key)


def _jacobi(
fn: Callable, x: NDArray, dx: NDArray, tol: float, debug: bool = False
) -> Tuple[NDArray, NDArray]:
assert x.ndim == 1
assert dx.ndim == 1
assert np.all(dx >= 0)
assert tol > 0
y = fn(x)
yrank = np.ndim(y)
jac = np.zeros((1 if yrank == 0 else len(y), len(x)))
h = np.zeros(len(x))
divergence = True
for i, hi in enumerate(dx):
if i > 0:
h[i - 1] = 0
if hi == 0:
continue
h[i] = hi
prev_esq = np.inf
for iter in range(20):
assert h[i] > 0
yu = fn(x + h)
yd = fn(x - h)
du = (yu - y) / h[i]
dd = (y - yd) / h[i]
d = 0.5 * (du + dd)
delta = du - dd
if np.all(np.abs(delta) <= tol * np.abs(d)):
if debug:
print(
f"jacobi: iter={iter} converged; delta={delta} "
f"threshold={tol * np.abs(d)}"
)
jac[:, i] = d
break
esq = np.dot(delta, delta)
if debug:
print(f"jacobi: iter={iter} d={d} esq={esq} h={h}")
if iter > 0 and esq < prev_esq:
divergence = False
if esq >= prev_esq:
if divergence:
print(f"jacobi: iter={iter} divergence detected")
jac[:, i] = np.nan
else:
print(f"jacobi: iter={iter} no convergence")
# no convergence, use previous more accurate jac[:, i]
break
jac[:, i] = d
prev_esq = esq
h[i] *= 0.1
return y, jac


@_deprecated.deprecated("use jacobi.propagate instead from jacobi library")
def propagate(
fn: Callable,
Expand Down Expand Up @@ -1031,17 +977,20 @@ def propagate(
y is the result of fn(x)
ycov is the propagated covariance matrix.
"""
from scipy.optimize import approx_fprime

vx = np.atleast_1d(x) # type:ignore
if np.ndim(cov) != 2: # type:ignore
raise ValueError("cov must be 2D array-like")
vcov = np.atleast_2d(cov) # type:ignore
if vcov.shape[0] != vcov.shape[1]:
raise ValueError("cov must have shape (N, N)")
tol = 1e-2
tol = 1e-10
dx = (np.diag(vcov) * tol) ** 0.5
if not np.all(dx >= 0):
raise ValueError("diagonal elements of covariance matrix must be non-negative")
y, jac = _jacobi(fn, vx, dx, tol)
y = fn(vx)
jac = np.atleast_2d(approx_fprime(vx, fn, dx))
ycov = np.einsum("ij,kl,jl", jac, jac, vcov)
return y, np.squeeze(ycov) if np.ndim(y) == 0 else ycov

Expand Down
117 changes: 17 additions & 100 deletions tests/test_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,13 @@
import pickle
from iminuit._hide_modules import hide_modules

try:
import scipy # noqa

scipy_available = True
except ModuleNotFoundError:
scipy_available = False


def test_ndim():
ndim = util._ndim
Expand Down Expand Up @@ -511,6 +518,7 @@ def g(x, a, b):
assert pg == (0, 3, 4)


@pytest.mark.skipif(not scipy_available, reason="needs scipy")
def test_propagate_1():
cov = [
[1.0, 0.1, 0.2],
Expand All @@ -525,14 +533,17 @@ def fn(x):
with pytest.warns(np.VisibleDeprecationWarning):
y, ycov = util.propagate(fn, x, cov)
np.testing.assert_allclose(y, [3, 5, 7])
np.testing.assert_allclose(ycov, [[4, 0.4, 0.8], [0.4, 8, 1.2], [0.8, 1.2, 12]])
np.testing.assert_allclose(
ycov, [[4, 0.4, 0.8], [0.4, 8, 1.2], [0.8, 1.2, 12]], rtol=1e-3
)

with pytest.warns(np.VisibleDeprecationWarning):
y, ycov = util.propagate(fn, [1], [[2]])
np.testing.assert_allclose(y, 3)
np.testing.assert_allclose(ycov, 8)
np.testing.assert_allclose(ycov, 8, rtol=1e-3)


@pytest.mark.skipif(not scipy_available, reason="needs scipy")
def test_propagate_2():
cov = [
[1.0, 0.1, 0.2],
Expand All @@ -549,7 +560,7 @@ def fn(x):
with pytest.warns(np.VisibleDeprecationWarning):
y, ycov = util.propagate(fn, x, cov)
np.testing.assert_equal(y, fn(x))
np.testing.assert_allclose(ycov, np.einsum("ij,kl,jl", a, a, cov))
np.testing.assert_allclose(ycov, np.einsum("ij,kl,jl", a, a, cov), rtol=1e-3)

def fn(x):
return np.linalg.multi_dot([x.T, cov, x])
Expand All @@ -558,9 +569,10 @@ def fn(x):
y, ycov = util.propagate(fn, x, cov)
np.testing.assert_equal(y, fn(np.array(x)))
jac = 2 * np.dot(cov, x)
np.testing.assert_allclose(ycov, np.einsum("i,k,ik", jac, jac, cov))
np.testing.assert_allclose(ycov, np.einsum("i,k,ik", jac, jac, cov), rtol=1e-3)


@pytest.mark.skipif(not scipy_available, reason="needs scipy")
def test_propagate_3():
# matrices with full zero rows and columns are supported
cov = [
Expand All @@ -579,6 +591,7 @@ def fn(x):
np.testing.assert_allclose(ycov, [[4, 0.0, 0.8], [0.0, 0.0, 0.0], [0.8, 0.0, 12]])


@pytest.mark.skipif(not scipy_available, reason="needs scipy")
def test_propagate_on_bad_input():
cov = [[np.nan, 0.0], [0.0, 1.0]]
x = [1.0, 2.0]
Expand All @@ -600,102 +613,6 @@ def fn(x):
util.propagate(fn, x, cov)


def test_jacobi():
n = 100
x = np.linspace(-5, 5, n)
dx = np.abs(x) * 1e-3 + 1e-3
tol = 1e-3

y, jac = util._jacobi(np.exp, x, dx, tol)
np.testing.assert_equal(y, np.exp(x))
jac_ref = np.zeros((n, n))
for i in range(n):
jac_ref[i, i] = np.exp(x[i])
np.testing.assert_allclose(jac, jac_ref)

x = np.linspace(1e-10, 5, n)
dx = np.abs(x) * 1e-3
y, jac = util._jacobi(np.log, x, dx, tol)
jac_ref = np.zeros((n, n))
for i in range(n):
jac_ref[i, i] = 1 / x[i]
np.testing.assert_allclose(jac, jac_ref)


def test_jacobi_on_bad_input():
x = np.array([1])
dx = np.array([0.1])
y, jac = util._jacobi(lambda x: np.nan, x, dx, 1e-3)

np.testing.assert_equal(y, np.nan)
np.testing.assert_equal(jac, np.nan)

x = np.array([1])
dx = np.array([0])
y, jac = util._jacobi(lambda x: x**2, x, dx, 1e-3)

np.testing.assert_equal(y, 1)
np.testing.assert_equal(jac, 0)

x = np.array([np.nan])
dx = np.array([0.1])
y, jac = util._jacobi(lambda x: x**2, x, dx, 1e-3)

np.testing.assert_equal(y, np.nan)
np.testing.assert_equal(jac, np.nan)

x = np.array([1])
dx = np.array([np.nan])
with pytest.raises(AssertionError):
# dx must be >= 0
util._jacobi(lambda x: x**2, x, dx, 1e-3)


@pytest.mark.parametrize("fail", (False, True))
def test_jacobi_low_resolution(fail, capsys):
x = np.array([1])
dx = np.array([1])

y, jac = util._jacobi(
lambda x: np.exp(x.astype(np.float32)),
x,
dx,
1e-10 if fail else 1e-3,
debug=True,
)

assert fail == ("no convergence" in capsys.readouterr()[0])

np.testing.assert_allclose(y, np.exp(1), rtol=1e-3)
np.testing.assert_allclose(jac, np.exp(1), rtol=1e-3)


def test_jacobi_divergence_1(capsys):
rng = np.random.default_rng(1)

x = np.array([1])
dx = np.array([1])

y, jac = util._jacobi(lambda x: rng.normal(), x, dx, 0.1, debug=True)

assert "divergence" in capsys.readouterr()[0]

np.testing.assert_allclose(y, 0, atol=1)
np.testing.assert_equal(jac, np.nan)


def test_jacobi_divergence_2(capsys):
x = np.array([1e-10])
dx = np.array([0.1])

y, jac = util._jacobi(lambda x: 1 / x, x, dx, 1e-3, debug=True)

assert "divergence" in capsys.readouterr()[0]

np.testing.assert_allclose(y, 1e10)
np.testing.assert_equal(jac, np.nan)


def test_iterate():
assert list(util._iterate(1)) == [1]
assert list(util._iterate([1, 2])) == [1, 2]
Expand Down

0 comments on commit 70036f7

Please sign in to comment.