diff --git a/tests/test_approximations.py b/tests/test_approximations.py index 28b09873..0581a6cd 100644 --- a/tests/test_approximations.py +++ b/tests/test_approximations.py @@ -119,7 +119,7 @@ def test_sufficient_statistics(self, pars): assert np.isclose(ln_t_j, ck_ln_t_j, rtol=1e-3) def test_mean_and_variance(self, pars): - logconst, t_i, var_t_i, t_j, var_t_j = approx.mean_and_variance(*pars) + logconst, t_i, _, var_t_i, t_j, _, var_t_j = approx.mean_and_variance(*pars) ck_normconst = scipy.integrate.dblquad( lambda ti, tj: self.pdf(ti, tj, *pars), 0, @@ -250,80 +250,11 @@ def test_approximate_gamma(self, k): assert kl_1 < kl_0 -@pytest.mark.parametrize( - "pars", - [ - [1.62, 0.00074, 25603.8, 0.6653, 0.0, 0.0011], # "Cancellation error" - ], -) -class Test2F1Failsafe: - """ - Test approximation of marginal pairwise joint distributions by a gamma via - arbitrary precision mean/variance matching, when sufficient statistics - calculation fails - """ - - def test_sufficient_statistics_throws_exception(self, pars): - with pytest.raises(Exception, match="Cancellation error"): - approx.sufficient_statistics(*pars) - - def test_exception_uses_mean_and_variance(self, pars): - _, t_i, va_t_i, t_j, va_t_j = approx.mean_and_variance(*pars) - ai1, bi1 = approx.approximate_gamma_mom(t_i, va_t_i) - aj1, bj1 = approx.approximate_gamma_mom(t_j, va_t_j) - _, par_i, par_j = approx.gamma_projection(*pars) - ai2, bi2 = par_i - aj2, bj2 = par_j - assert np.isclose(ai1, ai2) - assert np.isclose(bi1, bi2) - assert np.isclose(aj1, aj2) - assert np.isclose(bj1, bj2) - - class TestGammaFactorization: """ Test various functions for manipulating factorizations of gamma distributions """ - def test_rescale_gamma(self): - # posterior_shape = prior_shape + sum(in_shape - 1) + sum(out_shape - 1) - # posterior_rate = prior_rate + sum(in_rate) + sum(out_rate) - in_message = np.array([[1.5, 0.25], [1.5, 0.25]]) - out_message = np.array([[1.5, 0.25], [1.5, 0.25]]) - posterior = np.array([4, 1.5]) # prior is implicitly [2, 0.5] - prior = np.array( - [ - posterior[0] - - np.sum(in_message[:, 0] - 1) - - np.sum(out_message[:, 0] - 1), - posterior[1] - np.sum(in_message[:, 1]) - np.sum(out_message[:, 1]), - ] - ) - # rescale - target_shape = 12 - new_post, new_in, new_out = approx.rescale_gamma( - posterior, in_message, out_message, target_shape - ) - new_prior = np.array( - [ - new_post[0] - np.sum(new_in[:, 0] - 1) - np.sum(new_out[:, 0] - 1), - new_post[1] - np.sum(new_in[:, 1]) - np.sum(new_out[:, 1]), - ] - ) - print(prior, new_prior) - assert new_post[0] == target_shape - # mean is conserved - assert np.isclose(new_post[0] / new_post[1], posterior[0] / posterior[1]) - # magnitude of messages (in natural parameterization) is conserved - assert np.isclose( - (new_prior[0] - 1) / np.sum(new_in[:, 0] - 1), - (prior[0] - 1) / np.sum(in_message[:, 0] - 1), - ) - assert np.isclose( - new_prior[1] / np.sum(new_in[:, 1]), - prior[1] / np.sum(in_message[:, 1]), - ) - def test_average_gammas(self): # E[x] = shape/rate # E[log x] = digamma(shape) - log(rate) diff --git a/tests/test_hypergeo.py b/tests/test_hypergeo.py index 39285133..5f6d5875 100644 --- a/tests/test_hypergeo.py +++ b/tests/test_hypergeo.py @@ -56,100 +56,97 @@ def test_betaln(self, x): ) -@pytest.mark.parametrize( - "pars", - list( - itertools.product( - [0.8, 20.0, 200.0], - [1.9, 90.3, 900.3], - [1.6, 30.7, 300.7], - [0.0, 0.1, 0.45], - ) - ), -) -class TestTaylorSeries: +@pytest.mark.parametrize("a_i", [1.0, 10.0, 100.0, 1000.0]) +@pytest.mark.parametrize("b_i", [0.001, 0.01, 0.1, 1.0]) +@pytest.mark.parametrize("a_j", [1.0, 10.0, 100.0, 1000.0]) +@pytest.mark.parametrize("b_j", [0.001, 0.01, 0.1, 1.0]) +@pytest.mark.parametrize("y", [0.0, 1.0, 10.0, 1000.0]) +@pytest.mark.parametrize("mu", [0.005, 0.05, 0.5, 5.0]) +class TestLaplaceApprox: """ - Test Taylor series expansions of 2F1 + Test that Laplace approximation to 2F1 returns reasonable answers """ @staticmethod - def _2f1_validate(a, b, c, z, offset=1.0): - val = mpmath.re(mpmath.hyp2f1(a, b, c, z)) + def _2f1_validate(a_i, b_i, a_j, b_j, y, mu, offset=1.0): + A = a_j + B = a_i + a_j + y + C = a_j + y + 1 + z = (mu - b_j) / (mu + b_i) + val = mpmath.re(mpmath.hyp2f1(A, B, C, z, maxterms=1e6)) return val / offset - @staticmethod - def _2f1_grad_validate(a, b, c, z, offset=1.0): - p = [a, b, c, z] - grad = nd.Gradient( - lambda x: float(TestTaylorSeries._2f1_validate(*x, offset=offset)), - step=1e-7, - richardson_terms=4, + def test_2f1(self, a_i, b_i, a_j, b_j, y, mu): + pars = [a_i, b_i, a_j, b_j, y, mu] + f, *_ = hypergeo._hyp2f1(*pars) + check = float(mpmath.log(self._2f1_validate(*pars))) + assert np.isclose(f, check, rtol=2e-2) + + def test_grad(self, a_i, b_i, a_j, b_j, y, mu): + pars = [a_i, b_i, a_j, b_j, y, mu] + _, *grad = hypergeo._hyp2f1(*pars) + da_i = nd.Derivative( + lambda a_i: hypergeo._hyp2f1(a_i, b_i, a_j, b_j, y, mu)[0], step=1e-3 ) - return grad(p) + db_i = nd.Derivative( + lambda b_i: hypergeo._hyp2f1(a_i, b_i, a_j, b_j, y, mu)[0], step=1e-5 + ) + da_j = nd.Derivative( + lambda a_j: hypergeo._hyp2f1(a_i, b_i, a_j, b_j, y, mu)[0], step=1e-3 + ) + db_j = nd.Derivative( + lambda b_j: hypergeo._hyp2f1(a_i, b_i, a_j, b_j, y, mu)[0], step=1e-5 + ) + check = [da_i(a_i), db_i(b_i), da_j(a_j), db_j(b_j)] + assert np.allclose(grad, check, rtol=1e-3) - def test_2f1(self, pars): - f, s, *_ = hypergeo._hyp2f1_taylor_series(*pars) - check = self._2f1_validate(*pars) - assert s == mpmath.sign(check) - assert np.isclose(f, float(mpmath.log(mpmath.fabs(check)))) - def test_2f1_grad(self, pars): - _, _, *grad = hypergeo._hyp2f1_taylor_series(*pars) - grad = grad[:-1] - offset = self._2f1_validate(*pars) - check = self._2f1_grad_validate(*pars, offset=offset) - assert np.allclose(grad, check) +# ------------------------------------------------- # +# The routines below aren't used in tsdate anymore, # +# but may be useful in the future # +# ------------------------------------------------- # @pytest.mark.parametrize( "pars", list( itertools.product( - [0.8, 20.3, 200.2], - [0.0, 1.0, 10.0, 31.0], - [1.6, 30.5, 300.7], - [1.1, 1.5, 1.9, 4.2], + [0.8, 20.0, 200.0], + [1.9, 90.3, 900.3], + [1.6, 30.7, 300.7], + [0.0, 0.1, 0.45], ) ), ) -class TestRecurrence: +class TestTaylorSeries: """ - Test recurrence for 2F1 when one parameter is a negative integer + Test Taylor series expansions of 2F1 """ - @staticmethod - def _transform_pars(a, b, c, z): - return a, b, c + a, z - @staticmethod def _2f1_validate(a, b, c, z, offset=1.0): - val = mpmath.re(mpmath.hyp2f1(a, -b, c, z)) + val = mpmath.re(mpmath.hyp2f1(a, b, c, z)) return val / offset @staticmethod def _2f1_grad_validate(a, b, c, z, offset=1.0): p = [a, b, c, z] grad = nd.Gradient( - lambda x: float(TestRecurrence._2f1_validate(*x, offset=offset)), - step=1e-6, + lambda x: float(TestTaylorSeries._2f1_validate(*x, offset=offset)), + step=1e-7, richardson_terms=4, ) return grad(p) def test_2f1(self, pars): - pars = self._transform_pars(*pars) - f, s, *_ = hypergeo._hyp2f1_recurrence(*pars) + f, *_ = hypergeo._hyp2f1_taylor_series(*pars) check = self._2f1_validate(*pars) - assert s == mpmath.sign(check) assert np.isclose(f, float(mpmath.log(mpmath.fabs(check)))) def test_2f1_grad(self, pars): - pars = self._transform_pars(*pars) - _, _, *grad = hypergeo._hyp2f1_recurrence(*pars) - grad = grad[:-1] + _, *grad = hypergeo._hyp2f1_taylor_series(*pars) offset = self._2f1_validate(*pars) check = self._2f1_grad_validate(*pars, offset=offset) - check[1] = 0.0 # integer parameter has no gradient assert np.allclose(grad, check) @@ -192,93 +189,3 @@ def test_is_valid_2f1(self, pars): dz *= 1 + 1e-3 d2z *= 1 - 1e-3 assert not hypergeo._is_valid_2f1(dz, d2z, *pars, 1e-10) - - -@pytest.mark.parametrize("muts", [0.0, 1.0, 5.0, 10.0]) -@pytest.mark.parametrize( - "hyp2f1_func, pars", - [ - (hypergeo._hyp2f1_dlmf1521, [1.4, 0.018, 2.34, 2.3e-05, 0.0, 0.0395]), - (hypergeo._hyp2f1_dlmf1581, [1.4, 0.018, 20.3, 0.04, 0.0, 2.3e-05]), - (hypergeo._hyp2f1_dlmf1583, [5.4, 0.018, 10.34, 0.04, 0.0, 2.3e-05]), - ], -) -class TestTransforms: - """ - Test numerically stable transformations of hypergeometric functions - """ - - @staticmethod - def _2f1_validate(a_i, b_i, a_j, b_j, y, mu, offset=1.0): - A = a_j - B = a_i + a_j + y - C = a_j + y + 1 - z = (mu - b_j) / (mu + b_i) - val = mpmath.re(mpmath.hyp2f1(A, B, C, z, maxterms=1e6)) - return val / offset - - @staticmethod - def _2f1_grad_validate(a_i, b_i, a_j, b_j, y, mu, offset=1.0): - p = [a_i, b_i, a_j, b_j] - grad = nd.Gradient( - lambda x: float(TestTransforms._2f1_validate(*x, y, mu, offset=offset)), - step=1e-6, - richardson_terms=4, - ) - return grad(p) - - def test_2f1(self, muts, hyp2f1_func, pars): - pars[4] = muts - f, s, *_ = hyp2f1_func(*pars) - assert s > 0 - check = float(mpmath.log(self._2f1_validate(*pars))) - assert np.isclose(f, check) - - def test_2f1_grad(self, muts, hyp2f1_func, pars): - pars[4] = muts - _, s, *grad = hyp2f1_func(*pars) - assert s > 0 - offset = self._2f1_validate(*pars) - check = self._2f1_grad_validate(*pars, offset=offset) - assert np.allclose(grad, check) - - -@pytest.mark.parametrize( - "func, pars, err", - [ - [ - hypergeo._hyp2f1_dlmf1583, - [-21.62, 0.00074, 1003.8, 0.7653, 100.0, 0.0011], - "Cancellation error", - ], - [ - hypergeo._hyp2f1_dlmf1583, - [1.62, 0.00074, 25603.8, 0.6653, 0.0, 0.0011], - "Cancellation error", - ], - # TODO: gives zero function value, then reroutes through dlmf1581 - # [ - # hypergeo._hyp2f1_dlmf1583, - # [9007.39, 0.241, 10000, 0.2673, 2.0, 0.01019], - # "Cancellation error", - # ], - [ - hypergeo._hyp2f1_dlmf1581, - [1.62, 0.00074, 25603.8, 0.7653, 100.0, 0.0011], - "Maximum terms", - ], - [ - hypergeo._hyp2f1_dlmf1583, - [1.0, 1.0, 1.0, 1.0, 3.0, 0.0], - "Zero division", - ], - ], -) -class TestInvalid2F1: - """ - Test cases where homegrown 2F1 fails to converge - """ - - def test_hyp2f1_error(self, func, pars, err): - with pytest.raises(hypergeo.Invalid2F1, match=err): - func(*pars) diff --git a/tsdate/core.py b/tsdate/core.py index 238d65c1..abe0b8a5 100644 --- a/tsdate/core.py +++ b/tsdate/core.py @@ -1632,7 +1632,7 @@ def variational_dates( # convert priors to natural parameterization and average for n in priors.nonfixed_nodes: priors[n][0] -= 1.0 - assert priors[n][0] > -1.0 # TODO: throw error + assert priors[n][0] > -1.0 assert priors[n][1] > 0.0 if global_prior: logging.info("Pooling node-specific priors into global prior") diff --git a/tsdate/hypergeo.py b/tsdate/hypergeo.py index 1db02e41..6b9a6590 100644 --- a/tsdate/hypergeo.py +++ b/tsdate/hypergeo.py @@ -269,40 +269,40 @@ def _hyp2f1_laplace_approx(a, b, c, x): return f, df_da, df_db, df_dc, df_dx -@numba.njit("UniTuple(float64, 5)(float64, float64, float64, float64)") -def _hyp2f1_laplace_recurrence(a, b, c, x): - """ - Use contiguous relations to stabilixe the calculation of 2F1 - """ - - if x == 0.0: - return 0.0, 0.0, 0.0, 0.0, a * b / c - - A = x * (1 - x) * (a + 1) * (b + 1) / (c * (c + 1)) - dA_da = A / (a + 1) - dA_db = A / (b + 1) - dA_dc = -A / (c * (c + 1)) * (2 * c + 1) - - B = 1 - (a + b + 1) * x / c - dB_da = -x / c - dB_db = -x / c - dB_dc = (a + b + 1) * x / c**2 - - Af, dAf_da, dAf_db, dAf_dc, _ = _hyp2f1_laplace_approx(a + 2, b + 2, c + 2, x) - Bf, dBf_da, dBf_db, dBf_dc, _ = _hyp2f1_laplace_approx(a + 1, b + 1, c + 1, x) - Mx = max(Af, Bf) - Cf = np.log(A * np.exp(Af - Mx) + B * np.exp(Bf - Mx)) + Mx - - sA = np.exp(Af - Cf) - sB = np.exp(Bf - Cf) - - v = Cf - da = sA * (dA_da + A * dAf_da) + sB * (dB_da + B * dBf_da) - db = sA * (dA_db + A * dAf_db) + sB * (dB_db + B * dBf_db) - dc = sA * (dA_dc + A * dAf_dc) + sB * (dB_dc + B * dBf_dc) - dx = sB * a * b / c - - return v, da, db, dc, dx +# @numba.njit("UniTuple(float64, 5)(float64, float64, float64, float64)") +# def _hyp2f1_laplace_recurrence(a, b, c, x): +# """ +# Use contiguous relations to stabilize the calculation of 2F1 +# """ +# +# if x == 0.0: +# return 0.0, 0.0, 0.0, 0.0, a * b / c +# +# A = x * (1 - x) * (a + 1) * (b + 1) / (c * (c + 1)) +# dA_da = A / (a + 1) +# dA_db = A / (b + 1) +# dA_dc = -A / (c * (c + 1)) * (2 * c + 1) +# +# B = 1 - (a + b + 1) * x / c +# dB_da = -x / c +# dB_db = -x / c +# dB_dc = (a + b + 1) * x / c**2 +# +# Af, dAf_da, dAf_db, dAf_dc, _ = _hyp2f1_laplace_approx(a + 2, b + 2, c + 2, x) +# Bf, dBf_da, dBf_db, dBf_dc, _ = _hyp2f1_laplace_approx(a + 1, b + 1, c + 1, x) +# Mx = max(Af, Bf) +# Cf = np.log(A * np.exp(Af - Mx) + B * np.exp(Bf - Mx)) + Mx +# +# sA = np.exp(Af - Cf) +# sB = np.exp(Bf - Cf) +# +# v = Cf +# da = sA * (dA_da + A * dAf_da) + sB * (dB_da + B * dBf_da) +# db = sA * (dA_db + A * dAf_db) + sB * (dB_db + B * dBf_db) +# dc = sA * (dA_dc + A * dAf_dc) + sB * (dB_dc + B * dBf_dc) +# dx = sB * a * b / c +# +# return v, da, db, dc, dx @numba.njit(