From 2cd6a25f335f669b296c051fe73eae91b2a02d2d Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Fri, 9 Sep 2022 10:42:35 -0400 Subject: [PATCH 01/15] Split up test_gp.py --- .github/workflows/tests.yml | 7 +- pymc/tests/gp/__init__.py | 0 pymc/tests/{test_gp.py => gp/test_cov.py} | 593 ---------------------- pymc/tests/gp/test_gp.py | 457 +++++++++++++++++ pymc/tests/gp/test_mean.py | 90 ++++ pymc/tests/gp/test_util.py | 102 ++++ 6 files changed, 654 insertions(+), 595 deletions(-) create mode 100644 pymc/tests/gp/__init__.py rename pymc/tests/{test_gp.py => gp/test_cov.py} (51%) create mode 100644 pymc/tests/gp/test_gp.py create mode 100644 pymc/tests/gp/test_mean.py create mode 100644 pymc/tests/gp/test_util.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 218bce856e1..1e8b823284f 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -74,7 +74,10 @@ jobs: - | pymc/tests/distributions/test_timeseries.py - pymc/tests/test_gp.py + pymc/tests/gp/test_cov.py + pymc/tests/gp/test_gp.py + pymc/tests/gp/test_mean.py + pymc/tests/gp/test_util.py pymc/tests/test_model.py pymc/tests/test_model_graph.py pymc/tests/test_ode.py @@ -148,7 +151,7 @@ jobs: test-subset: - pymc/tests/test_variational_inference.py pymc/tests/test_initial_point.py - pymc/tests/test_pickling.py pymc/tests/test_profile.py pymc/tests/test_step.py - - pymc/tests/test_gp.py pymc/tests/test_ode.py pymc/tests/test_smc.py pymc/tests/test_parallel_sampling.py + - pymc/tests/gp/test_cov.py pymc/tests/gp/test_gp.py pymc/tests/gp/test_mean.py pymc/tests/gp/test_util.py pymc/tests/test_ode.py pymc/tests/test_smc.py pymc/tests/test_parallel_sampling.py - pymc/tests/test_sampling.py pymc/tests/test_posteriors.py fail-fast: false diff --git a/pymc/tests/gp/__init__.py b/pymc/tests/gp/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pymc/tests/test_gp.py b/pymc/tests/gp/test_cov.py similarity index 51% rename from pymc/tests/test_gp.py rename to pymc/tests/gp/test_cov.py index 999578b17cf..550f68727d5 100644 --- a/pymc/tests/test_gp.py +++ b/pymc/tests/gp/test_cov.py @@ -12,10 +12,6 @@ # See the License for the specific language governing permissions and # limitations under the License. -# pylint:disable=unused-variable -from functools import reduce -from operator import add - import aesara import aesara.tensor as at import numpy as np @@ -26,80 +22,6 @@ from pymc.math import cartesian, kronecker -np.random.seed(101) - - -class TestZeroMean: - def test_value(self): - X = np.linspace(0, 1, 10)[:, None] - with pm.Model() as model: - zero_mean = pm.gp.mean.Zero() - M = zero_mean(X).eval() - assert np.all(M == 0) - assert M.shape == (10,) - - -class TestConstantMean: - def test_value(self): - X = np.linspace(0, 1, 10)[:, None] - with pm.Model() as model: - const_mean = pm.gp.mean.Constant(6) - M = const_mean(X).eval() - assert np.all(M == 6) - assert M.shape == (10,) - - -class TestLinearMean: - def test_value(self): - X = np.linspace(0, 1, 10)[:, None] - with pm.Model() as model: - linear_mean = pm.gp.mean.Linear(2, 0.5) - M = linear_mean(X).eval() - npt.assert_allclose(M[1], 0.7222, atol=1e-3) - assert M.shape == (10,) - - -class TestAddProdMean: - def test_add(self): - X = np.linspace(0, 1, 10)[:, None] - with pm.Model() as model: - mean1 = pm.gp.mean.Linear(coeffs=2, intercept=0.5) - mean2 = pm.gp.mean.Constant(2) - mean = mean1 + mean2 + mean2 - M = mean(X).eval() - npt.assert_allclose(M[1], 0.7222 + 2 + 2, atol=1e-3) - - def test_prod(self): - X = np.linspace(0, 1, 10)[:, None] - with pm.Model() as model: - mean1 = pm.gp.mean.Linear(coeffs=2, intercept=0.5) - mean2 = pm.gp.mean.Constant(2) - mean = mean1 * mean2 * mean2 - M = mean(X).eval() - npt.assert_allclose(M[1], 0.7222 * 2 * 2, atol=1e-3) - - def test_add_multid(self): - X = np.linspace(0, 1, 30).reshape(10, 3) - A = np.array([1, 2, 3]) - b = 10 - with pm.Model() as model: - mean1 = pm.gp.mean.Linear(coeffs=A, intercept=b) - mean2 = pm.gp.mean.Constant(2) - mean = mean1 + mean2 + mean2 - M = mean(X).eval() - npt.assert_allclose(M[1], 10.8965 + 2 + 2, atol=1e-3) - - def test_prod_multid(self): - X = np.linspace(0, 1, 30).reshape(10, 3) - A = np.array([1, 2, 3]) - b = 10 - with pm.Model() as model: - mean1 = pm.gp.mean.Linear(coeffs=A, intercept=b) - mean2 = pm.gp.mean.Constant(2) - mean = mean1 * mean2 * mean2 - M = mean(X).eval() - npt.assert_allclose(M[1], 10.8965 * 2 * 2, atol=1e-3) - class TestCovAdd: def test_symadd_cov(self): @@ -793,518 +715,3 @@ def test_raises3(self): with pm.Model() as model: with pytest.raises(ValueError): B = pm.gp.cov.Coregion(1) - - -class TestMarginalVsLatent: - R""" - Compare the logp of models Marginal, noise=0 and Latent. - """ - - def setup_method(self): - X = np.random.randn(20, 3) - y = np.random.randn(20) - Xnew = np.random.randn(30, 3) - pnew = np.random.randn(30) - with pm.Model() as model: - cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - mean_func = pm.gp.mean.Constant(0.5) - gp = pm.gp.Marginal(mean_func=mean_func, cov_func=cov_func) - f = gp.marginal_likelihood("f", X, y, noise=0.0) - p = gp.conditional("p", Xnew) - self.logp = model.compile_logp()({"p": pnew}) - self.X = X - self.Xnew = Xnew - self.y = y - self.pnew = pnew - - def testLatent1(self): - with pm.Model() as model: - cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - mean_func = pm.gp.mean.Constant(0.5) - gp = pm.gp.Latent(mean_func=mean_func, cov_func=cov_func) - f = gp.prior("f", self.X, reparameterize=False) - p = gp.conditional("p", self.Xnew) - assert tuple(f.shape.eval()) == (self.X.shape[0],) - assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) - latent_logp = model.compile_logp()({"f": self.y, "p": self.pnew}) - npt.assert_allclose(latent_logp, self.logp, atol=0, rtol=1e-2) - - def testLatent2(self): - with pm.Model() as model: - cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - mean_func = pm.gp.mean.Constant(0.5) - gp = pm.gp.Latent(mean_func=mean_func, cov_func=cov_func) - f = gp.prior("f", self.X, reparameterize=True) - p = gp.conditional("p", self.Xnew) - assert tuple(f.shape.eval()) == (self.X.shape[0],) - assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) - chol = np.linalg.cholesky(cov_func(self.X).eval()) - y_rotated = np.linalg.solve(chol, self.y - 0.5) - latent_logp = model.compile_logp()({"f_rotated_": y_rotated, "p": self.pnew}) - npt.assert_allclose(latent_logp, self.logp, atol=5) - - -class TestMarginalVsMarginalApprox: - R""" - Compare test fits of models Marginal and MarginalApprox. - """ - - def setup_method(self): - self.sigma = 0.1 - self.x = np.linspace(-5, 5, 30) - self.y = np.random.normal(0.25 * self.x, self.sigma) - with pm.Model() as model: - cov_func = pm.gp.cov.Linear(1, c=0.0) - c = pm.Normal("c", mu=20.0, sigma=100.0) # far from true value - mean_func = pm.gp.mean.Constant(c) - self.gp = pm.gp.Marginal(mean_func=mean_func, cov_func=cov_func) - sigma = pm.HalfNormal("sigma", sigma=100) - self.gp.marginal_likelihood("lik", self.x[:, None], self.y, sigma) - self.map_full = pm.find_MAP(method="bfgs") # bfgs seems to work much better than lbfgsb - - self.x_new = np.linspace(-6, 6, 20) - - # Include additive Gaussian noise, return diagonal of predicted covariance matrix - with model: - self.pred_mu, self.pred_var = self.gp.predict( - self.x_new[:, None], point=self.map_full, pred_noise=True, diag=True - ) - - # Dont include additive Gaussian noise, return full predicted covariance matrix - with model: - self.pred_mu, self.pred_covar = self.gp.predict( - self.x_new[:, None], point=self.map_full, pred_noise=False, diag=False - ) - - @pytest.mark.parametrize("approx", ["FITC", "VFE", "DTC"]) - def test_fits_and_preds(self, approx): - """Get MAP estimate for GP approximation, compare results and predictions to what's returned - by an unapproximated GP. The tolerances are fairly wide, but narrow relative to initial - values of the unknown parameters. - """ - - with pm.Model() as model: - cov_func = pm.gp.cov.Linear(1, c=0.0) - c = pm.Normal("c", mu=20.0, sigma=100.0, initval=-500.0) - mean_func = pm.gp.mean.Constant(c) - gp = pm.gp.MarginalApprox(mean_func=mean_func, cov_func=cov_func, approx=approx) - sigma = pm.HalfNormal("sigma", sigma=100, initval=50.0) - gp.marginal_likelihood("lik", self.x[:, None], self.x[:, None], self.y, sigma) - map_approx = pm.find_MAP(method="bfgs") - - # Check MAP gets approximately correct result - npt.assert_allclose(self.map_full["c"], map_approx["c"], atol=0.01, rtol=0.1) - npt.assert_allclose(self.map_full["sigma"], map_approx["sigma"], atol=0.01, rtol=0.1) - - # Check that predict (and conditional) work, include noise, with diagonal non-full pred var. - with model: - pred_mu_approx, pred_var_approx = gp.predict( - self.x_new[:, None], point=map_approx, pred_noise=True, diag=True - ) - npt.assert_allclose(self.pred_mu, pred_mu_approx, atol=0.0, rtol=0.1) - npt.assert_allclose(self.pred_var, pred_var_approx, atol=0.0, rtol=0.1) - - # Check that predict (and conditional) work, no noise, full pred covariance. - with model: - pred_mu_approx, pred_var_approx = gp.predict( - self.x_new[:, None], point=map_approx, pred_noise=True, diag=True - ) - npt.assert_allclose(self.pred_mu, pred_mu_approx, atol=0.0, rtol=0.1) - npt.assert_allclose(self.pred_var, pred_var_approx, atol=0.0, rtol=0.1) - - -class TestGPAdditive: - def setup_method(self): - self.X = np.random.randn(50, 3) - self.y = np.random.randn(50) - self.Xnew = np.random.randn(60, 3) - self.noise = pm.gp.cov.WhiteNoise(0.1) - self.covs = ( - pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]), - pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]), - pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]), - ) - self.means = (pm.gp.mean.Constant(0.5), pm.gp.mean.Constant(0.5), pm.gp.mean.Constant(0.5)) - - def testAdditiveMarginal(self): - with pm.Model() as model1: - gp1 = pm.gp.Marginal(mean_func=self.means[0], cov_func=self.covs[0]) - gp2 = pm.gp.Marginal(mean_func=self.means[1], cov_func=self.covs[1]) - gp3 = pm.gp.Marginal(mean_func=self.means[2], cov_func=self.covs[2]) - - gpsum = gp1 + gp2 + gp3 - fsum = gpsum.marginal_likelihood("f", self.X, self.y, noise=self.noise) - model1_logp = model1.compile_logp()({}) - - with pm.Model() as model2: - gptot = pm.gp.Marginal( - mean_func=reduce(add, self.means), cov_func=reduce(add, self.covs) - ) - fsum = gptot.marginal_likelihood("f", self.X, self.y, noise=self.noise) - model2_logp = model2.compile_logp()({}) - npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) - - with model1: - fp1 = gpsum.conditional( - "fp1", self.Xnew, given={"X": self.X, "y": self.y, "noise": self.noise, "gp": gpsum} - ) - with model2: - fp2 = gptot.conditional("fp2", self.Xnew) - - fp = np.random.randn(self.Xnew.shape[0]) - logp1 = model1.compile_logp()({"fp1": fp}) - logp2 = model2.compile_logp()({"fp2": fp}) - npt.assert_allclose(logp1, logp2, atol=0, rtol=1e-2) - - @pytest.mark.parametrize("approx", ["FITC", "VFE", "DTC"]) - def testAdditiveMarginalApprox(self, approx): - Xu = np.random.randn(10, 3) - sigma = 0.1 - with pm.Model() as model1: - gp1 = pm.gp.MarginalApprox( - mean_func=self.means[0], cov_func=self.covs[0], approx=approx - ) - gp2 = pm.gp.MarginalApprox( - mean_func=self.means[1], cov_func=self.covs[1], approx=approx - ) - gp3 = pm.gp.MarginalApprox( - mean_func=self.means[2], cov_func=self.covs[2], approx=approx - ) - - gpsum = gp1 + gp2 + gp3 - fsum = gpsum.marginal_likelihood("f", self.X, Xu, self.y, noise=sigma) - model1_logp = model1.compile_logp()({}) - - with pm.Model() as model2: - gptot = pm.gp.MarginalApprox( - mean_func=reduce(add, self.means), cov_func=reduce(add, self.covs), approx=approx - ) - fsum = gptot.marginal_likelihood("f", self.X, Xu, self.y, noise=sigma) - model2_logp = model2.compile_logp()({}) - npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) - - with model1: - fp1 = gpsum.conditional( - "fp1", - self.Xnew, - given={"X": self.X, "Xu": Xu, "y": self.y, "sigma": sigma, "gp": gpsum}, - ) - with model2: - fp2 = gptot.conditional("fp2", self.Xnew) - - fp = np.random.randn(self.Xnew.shape[0]) - - model1_logp = model1.compile_logp()({"fp1": fp}) - model2_logp = model2.compile_logp()({"fp2": fp}) - - npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) - - def testAdditiveLatent(self): - with pm.Model() as model1: - gp1 = pm.gp.Latent(mean_func=self.means[0], cov_func=self.covs[0]) - gp2 = pm.gp.Latent(mean_func=self.means[1], cov_func=self.covs[1]) - gp3 = pm.gp.Latent(mean_func=self.means[2], cov_func=self.covs[2]) - - gpsum = gp1 + gp2 + gp3 - fsum = gpsum.prior("fsum", self.X, reparameterize=False) - model1_logp = model1.compile_logp()({"fsum": self.y}) - - with pm.Model() as model2: - gptot = pm.gp.Latent(mean_func=reduce(add, self.means), cov_func=reduce(add, self.covs)) - fsum = gptot.prior("fsum", self.X, reparameterize=False) - model2_logp = model2.compile_logp()({"fsum": self.y}) - npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) - - with model1: - fp1 = gpsum.conditional("fp1", self.Xnew, given={"X": self.X, "f": self.y, "gp": gpsum}) - with model2: - fp2 = gptot.conditional("fp2", self.Xnew) - - fp = np.random.randn(self.Xnew.shape[0]) - logp1 = model1.compile_logp()({"fsum": self.y, "fp1": fp}) - logp2 = model2.compile_logp()({"fsum": self.y, "fp2": fp}) - npt.assert_allclose(logp1, logp2, atol=0, rtol=1e-2) - - def testAdditiveSparseRaises(self): - # cant add different approximations - with pm.Model() as model: - cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - gp1 = pm.gp.MarginalApprox(cov_func=cov_func, approx="DTC") - gp2 = pm.gp.MarginalApprox(cov_func=cov_func, approx="FITC") - with pytest.raises(Exception) as e_info: - gp1 + gp2 - - def testAdditiveTypeRaises1(self): - with pm.Model() as model: - cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - gp1 = pm.gp.MarginalApprox(cov_func=cov_func, approx="DTC") - gp2 = pm.gp.Marginal(cov_func=cov_func) - with pytest.raises(Exception) as e_info: - gp1 + gp2 - - def testAdditiveTypeRaises2(self): - with pm.Model() as model: - cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - gp1 = pm.gp.Latent(cov_func=cov_func) - gp2 = pm.gp.Marginal(cov_func=cov_func) - with pytest.raises(Exception) as e_info: - gp1 + gp2 - - -class TestTP: - R""" - Compare TP with high degress of freedom to GP - """ - - def setup_method(self): - X = np.random.randn(20, 3) - y = np.random.randn(20) - Xnew = np.random.randn(30, 3) - pnew = np.random.randn(30) - - with pm.Model() as model1: - cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - gp = pm.gp.Latent(cov_func=cov_func) - f = gp.prior("f", X, reparameterize=False) - p = gp.conditional("p", Xnew) - self.gp_latent_logp = model1.compile_logp()({"f": y, "p": pnew}) - self.X = X - self.y = y - self.Xnew = Xnew - self.pnew = pnew - self.nu = 10000 - - def testTPvsLatent(self): - with pm.Model() as model: - scale_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - tp = pm.gp.TP(scale_func=scale_func, nu=self.nu) - f = tp.prior("f", self.X, reparameterize=False) - p = tp.conditional("p", self.Xnew) - assert tuple(f.shape.eval()) == (self.X.shape[0],) - assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) - tp_logp = model.compile_logp()({"f": self.y, "p": self.pnew}) - npt.assert_allclose(self.gp_latent_logp, tp_logp, atol=0, rtol=1e-2) - - def testTPvsLatentReparameterized(self): - with pm.Model() as model: - scale_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - tp = pm.gp.TP(scale_func=scale_func, nu=self.nu) - f = tp.prior("f", self.X, reparameterize=True) - p = tp.conditional("p", self.Xnew) - assert tuple(f.shape.eval()) == (self.X.shape[0],) - assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) - chol = np.linalg.cholesky(scale_func(self.X).eval()) - f_rotated = np.linalg.solve(chol, self.y) - tp_logp = model.compile_logp()({"f_rotated_": f_rotated, "p": self.pnew}) - npt.assert_allclose(self.gp_latent_logp, tp_logp, atol=0, rtol=1e-2) - - def testAdditiveTPRaises(self): - with pm.Model() as model: - scale_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) - gp1 = pm.gp.TP(scale_func=scale_func, nu=10) - gp2 = pm.gp.TP(scale_func=scale_func, nu=10) - with pytest.raises(Exception) as e_info: - gp1 + gp2 - - -class TestLatentKron: - """ - Compare gp.LatentKron to gp.Latent, both with Gaussian noise. - """ - - def setup_method(self): - self.Xs = [ - np.linspace(0, 1, 7)[:, None], - np.linspace(0, 1, 5)[:, None], - np.linspace(0, 1, 6)[:, None], - ] - self.X = cartesian(*self.Xs) - self.N = np.prod([len(X) for X in self.Xs]) - self.y = np.random.randn(self.N) * 0.1 - self.Xnews = (np.random.randn(5, 1), np.random.randn(5, 1), np.random.randn(5, 1)) - self.Xnew = np.concatenate(self.Xnews, axis=1) - self.pnew = np.random.randn(len(self.Xnew)) - ls = 0.2 - with pm.Model() as latent_model: - self.cov_funcs = ( - pm.gp.cov.ExpQuad(1, ls), - pm.gp.cov.ExpQuad(1, ls), - pm.gp.cov.ExpQuad(1, ls), - ) - cov_func = pm.gp.cov.Kron(self.cov_funcs) - self.mean = pm.gp.mean.Constant(0.5) - gp = pm.gp.Latent(mean_func=self.mean, cov_func=cov_func) - f = gp.prior("f", self.X) - p = gp.conditional("p", self.Xnew) - chol = np.linalg.cholesky(cov_func(self.X).eval()) - self.y_rotated = np.linalg.solve(chol, self.y - 0.5) - self.logp = latent_model.compile_logp()({"f_rotated_": self.y_rotated, "p": self.pnew}) - - def testLatentKronvsLatent(self): - with pm.Model() as kron_model: - kron_gp = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - f = kron_gp.prior("f", self.Xs) - p = kron_gp.conditional("p", self.Xnew) - assert tuple(f.shape.eval()) == (self.X.shape[0],) - assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) - kronlatent_logp = kron_model.compile_logp()({"f_rotated_": self.y_rotated, "p": self.pnew}) - npt.assert_allclose(kronlatent_logp, self.logp, atol=0, rtol=1e-3) - - def testLatentKronRaisesAdditive(self): - with pm.Model() as kron_model: - gp1 = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - gp2 = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - with pytest.raises(TypeError): - gp1 + gp2 - - def testLatentKronRaisesSizes(self): - with pm.Model() as kron_model: - gp = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - with pytest.raises(ValueError): - gp.prior("f", Xs=[np.linspace(0, 1, 7)[:, None], np.linspace(0, 1, 5)[:, None]]) - - -class TestMarginalKron: - """ - Compare gp.MarginalKron to gp.Marginal. - """ - - def setup_method(self): - self.Xs = [ - np.linspace(0, 1, 7)[:, None], - np.linspace(0, 1, 5)[:, None], - np.linspace(0, 1, 6)[:, None], - ] - self.X = cartesian(*self.Xs) - self.N = np.prod([len(X) for X in self.Xs]) - self.y = np.random.randn(self.N) * 0.1 - self.Xnews = (np.random.randn(5, 1), np.random.randn(5, 1), np.random.randn(5, 1)) - self.Xnew = np.concatenate(self.Xnews, axis=1) - self.sigma = 0.2 - self.pnew = np.random.randn(len(self.Xnew)) - - ls = 0.2 - with pm.Model() as model: - self.cov_funcs = [ - pm.gp.cov.ExpQuad(1, ls), - pm.gp.cov.ExpQuad(1, ls), - pm.gp.cov.ExpQuad(1, ls), - ] - cov_func = pm.gp.cov.Kron(self.cov_funcs) - self.mean = pm.gp.mean.Constant(0.5) - gp = pm.gp.Marginal(mean_func=self.mean, cov_func=cov_func) - f = gp.marginal_likelihood("f", self.X, self.y, noise=self.sigma) - p = gp.conditional("p", self.Xnew) - self.mu, self.cov = gp.predict(self.Xnew) - self.logp = model.compile_logp()({"p": self.pnew}) - - def testMarginalKronvsMarginalpredict(self): - with pm.Model() as kron_model: - kron_gp = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - f = kron_gp.marginal_likelihood("f", self.Xs, self.y, sigma=self.sigma) - p = kron_gp.conditional("p", self.Xnew) - mu, cov = kron_gp.predict(self.Xnew) - assert tuple(f.shape.eval()) == (self.X.shape[0],) - assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) - npt.assert_allclose(mu, self.mu, atol=1e-5, rtol=1e-2) - npt.assert_allclose(cov, self.cov, atol=1e-5, rtol=1e-2) - with kron_model: - _, var = kron_gp.predict(self.Xnew, diag=True) - npt.assert_allclose(np.diag(cov), var, atol=1e-5, rtol=1e-2) - - def testMarginalKronvsMarginal(self): - with pm.Model() as kron_model: - kron_gp = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - f = kron_gp.marginal_likelihood("f", self.Xs, self.y, sigma=self.sigma) - p = kron_gp.conditional("p", self.Xnew) - kron_logp = kron_model.compile_logp()({"p": self.pnew}) - npt.assert_allclose(kron_logp, self.logp, atol=0, rtol=1e-2) - - def testMarginalKronRaises(self): - with pm.Model() as kron_model: - gp1 = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - gp2 = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) - with pytest.raises(TypeError): - gp1 + gp2 - - -class TestPlotGP: - def test_plot_gp_dist(self): - """Test that the plotting helper works with the stated input shapes.""" - import matplotlib.pyplot as plt - - X = 100 - S = 500 - fig, ax = plt.subplots() - pm.gp.util.plot_gp_dist( - ax, x=np.linspace(0, 50, X), samples=np.random.normal(np.arange(X), size=(S, X)) - ) - plt.close() - pass - - def test_plot_gp_dist_warn_nan(self): - """Test that the plotting helper works with the stated input shapes.""" - import matplotlib.pyplot as plt - - X = 100 - S = 500 - samples = np.random.normal(np.arange(X), size=(S, X)) - samples[15, 3] = np.nan - fig, ax = plt.subplots() - with pytest.warns(UserWarning): - pm.gp.util.plot_gp_dist(ax, x=np.linspace(0, 50, X), samples=samples) - plt.close() - pass - - -class TestKmeansInducing: - def setup_method(self): - self.centers = (-5, 5) - self.x = np.concatenate( - (self.centers[0] + np.random.randn(500), self.centers[1] + np.random.randn(500)) - ) - - def test_kmeans(self): - X = self.x[:, None] - Xu = pm.gp.util.kmeans_inducing_points(2, X).flatten() - npt.assert_allclose(np.asarray(self.centers), np.sort(Xu), atol=0.1) - - X = at.as_tensor_variable(self.x[:, None]) - Xu = pm.gp.util.kmeans_inducing_points(2, X).flatten() - npt.assert_allclose(np.asarray(self.centers), np.sort(Xu), atol=0.1) - - def test_kmeans_raises(self): - with pytest.raises(TypeError): - Xu = pm.gp.util.kmeans_inducing_points(2, "str is the wrong type").flatten() - - -class TestReplaceWithValues: - def test_basic_replace(self): - with pm.Model() as model: - a = pm.Normal("a") - b = pm.Normal("b", mu=a) - c = a * b - - (c_val,) = pm.gp.util.replace_with_values( - [c], replacements={"a": 2, "b": 3, "x": 100}, model=model - ) - assert c_val == np.array(6.0) - - def test_replace_no_inputs_needed(self): - with pm.Model() as model: - a = at.as_tensor_variable(2.0) - b = 1.0 + a - c = a * b - (c_val,) = pm.gp.util.replace_with_values([c], replacements={"x": 100}) - assert c_val == np.array(6.0) - - def test_missing_input(self): - with pm.Model() as model: - a = pm.Normal("a") - b = pm.Normal("b", mu=a) - c = a * b - - with pytest.raises(ValueError): - (c_val,) = pm.gp.util.replace_with_values( - [c], replacements={"a": 2, "x": 100}, model=model - ) diff --git a/pymc/tests/gp/test_gp.py b/pymc/tests/gp/test_gp.py new file mode 100644 index 00000000000..c56d24336db --- /dev/null +++ b/pymc/tests/gp/test_gp.py @@ -0,0 +1,457 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from functools import reduce +from operator import add + +import numpy as np +import numpy.testing as npt +import pytest + +import pymc as pm + +from pymc.math import cartesian + + +class TestMarginalVsMarginalApprox: + R""" + Compare test fits of models Marginal and MarginalApprox. + """ + + def setup_method(self): + self.sigma = 0.1 + self.x = np.linspace(-5, 5, 30) + self.y = np.random.normal(0.25 * self.x, self.sigma) + with pm.Model() as model: + cov_func = pm.gp.cov.Linear(1, c=0.0) + c = pm.Normal("c", mu=20.0, sigma=100.0) # far from true value + mean_func = pm.gp.mean.Constant(c) + self.gp = pm.gp.Marginal(mean_func=mean_func, cov_func=cov_func) + sigma = pm.HalfNormal("sigma", sigma=100) + self.gp.marginal_likelihood("lik", self.x[:, None], self.y, sigma) + self.map_full = pm.find_MAP(method="bfgs") # bfgs seems to work much better than lbfgsb + + self.x_new = np.linspace(-6, 6, 20) + + # Include additive Gaussian noise, return diagonal of predicted covariance matrix + with model: + self.pred_mu, self.pred_var = self.gp.predict( + self.x_new[:, None], point=self.map_full, pred_noise=True, diag=True + ) + + # Dont include additive Gaussian noise, return full predicted covariance matrix + with model: + self.pred_mu, self.pred_covar = self.gp.predict( + self.x_new[:, None], point=self.map_full, pred_noise=False, diag=False + ) + + @pytest.mark.parametrize("approx", ["FITC", "VFE", "DTC"]) + def test_fits_and_preds(self, approx): + """Get MAP estimate for GP approximation, compare results and predictions to what's returned + by an unapproximated GP. The tolerances are fairly wide, but narrow relative to initial + values of the unknown parameters. + """ + + with pm.Model() as model: + cov_func = pm.gp.cov.Linear(1, c=0.0) + c = pm.Normal("c", mu=20.0, sigma=100.0, initval=-500.0) + mean_func = pm.gp.mean.Constant(c) + gp = pm.gp.MarginalApprox(mean_func=mean_func, cov_func=cov_func, approx=approx) + sigma = pm.HalfNormal("sigma", sigma=100, initval=50.0) + gp.marginal_likelihood("lik", self.x[:, None], self.x[:, None], self.y, sigma) + map_approx = pm.find_MAP(method="bfgs") + + # Check MAP gets approximately correct result + npt.assert_allclose(self.map_full["c"], map_approx["c"], atol=0.01, rtol=0.1) + npt.assert_allclose(self.map_full["sigma"], map_approx["sigma"], atol=0.01, rtol=0.1) + + # Check that predict (and conditional) work, include noise, with diagonal non-full pred var. + with model: + pred_mu_approx, pred_var_approx = gp.predict( + self.x_new[:, None], point=map_approx, pred_noise=True, diag=True + ) + npt.assert_allclose(self.pred_mu, pred_mu_approx, atol=0.0, rtol=0.1) + npt.assert_allclose(self.pred_var, pred_var_approx, atol=0.0, rtol=0.1) + + # Check that predict (and conditional) work, no noise, full pred covariance. + with model: + pred_mu_approx, pred_var_approx = gp.predict( + self.x_new[:, None], point=map_approx, pred_noise=True, diag=True + ) + npt.assert_allclose(self.pred_mu, pred_mu_approx, atol=0.0, rtol=0.1) + npt.assert_allclose(self.pred_var, pred_var_approx, atol=0.0, rtol=0.1) + + +class TestGPAdditive: + def setup_method(self): + self.X = np.random.randn(50, 3) + self.y = np.random.randn(50) + self.Xnew = np.random.randn(60, 3) + self.noise = pm.gp.cov.WhiteNoise(0.1) + self.covs = ( + pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]), + pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]), + pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]), + ) + self.means = (pm.gp.mean.Constant(0.5), pm.gp.mean.Constant(0.5), pm.gp.mean.Constant(0.5)) + + def testAdditiveMarginal(self): + with pm.Model() as model1: + gp1 = pm.gp.Marginal(mean_func=self.means[0], cov_func=self.covs[0]) + gp2 = pm.gp.Marginal(mean_func=self.means[1], cov_func=self.covs[1]) + gp3 = pm.gp.Marginal(mean_func=self.means[2], cov_func=self.covs[2]) + + gpsum = gp1 + gp2 + gp3 + fsum = gpsum.marginal_likelihood("f", self.X, self.y, noise=self.noise) + model1_logp = model1.compile_logp()({}) + + with pm.Model() as model2: + gptot = pm.gp.Marginal( + mean_func=reduce(add, self.means), cov_func=reduce(add, self.covs) + ) + fsum = gptot.marginal_likelihood("f", self.X, self.y, noise=self.noise) + model2_logp = model2.compile_logp()({}) + npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) + + with model1: + fp1 = gpsum.conditional( + "fp1", self.Xnew, given={"X": self.X, "y": self.y, "noise": self.noise, "gp": gpsum} + ) + with model2: + fp2 = gptot.conditional("fp2", self.Xnew) + + fp = np.random.randn(self.Xnew.shape[0]) + logp1 = model1.compile_logp()({"fp1": fp}) + logp2 = model2.compile_logp()({"fp2": fp}) + npt.assert_allclose(logp1, logp2, atol=0, rtol=1e-2) + + @pytest.mark.parametrize("approx", ["FITC", "VFE", "DTC"]) + def testAdditiveMarginalApprox(self, approx): + Xu = np.random.randn(10, 3) + sigma = 0.1 + with pm.Model() as model1: + gp1 = pm.gp.MarginalApprox( + mean_func=self.means[0], cov_func=self.covs[0], approx=approx + ) + gp2 = pm.gp.MarginalApprox( + mean_func=self.means[1], cov_func=self.covs[1], approx=approx + ) + gp3 = pm.gp.MarginalApprox( + mean_func=self.means[2], cov_func=self.covs[2], approx=approx + ) + + gpsum = gp1 + gp2 + gp3 + fsum = gpsum.marginal_likelihood("f", self.X, Xu, self.y, noise=sigma) + model1_logp = model1.compile_logp()({}) + + with pm.Model() as model2: + gptot = pm.gp.MarginalApprox( + mean_func=reduce(add, self.means), cov_func=reduce(add, self.covs), approx=approx + ) + fsum = gptot.marginal_likelihood("f", self.X, Xu, self.y, noise=sigma) + model2_logp = model2.compile_logp()({}) + npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) + + with model1: + fp1 = gpsum.conditional( + "fp1", + self.Xnew, + given={"X": self.X, "Xu": Xu, "y": self.y, "sigma": sigma, "gp": gpsum}, + ) + with model2: + fp2 = gptot.conditional("fp2", self.Xnew) + + fp = np.random.randn(self.Xnew.shape[0]) + + model1_logp = model1.compile_logp()({"fp1": fp}) + model2_logp = model2.compile_logp()({"fp2": fp}) + + npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) + + def testAdditiveLatent(self): + with pm.Model() as model1: + gp1 = pm.gp.Latent(mean_func=self.means[0], cov_func=self.covs[0]) + gp2 = pm.gp.Latent(mean_func=self.means[1], cov_func=self.covs[1]) + gp3 = pm.gp.Latent(mean_func=self.means[2], cov_func=self.covs[2]) + + gpsum = gp1 + gp2 + gp3 + fsum = gpsum.prior("fsum", self.X, reparameterize=False) + model1_logp = model1.compile_logp()({"fsum": self.y}) + + with pm.Model() as model2: + gptot = pm.gp.Latent(mean_func=reduce(add, self.means), cov_func=reduce(add, self.covs)) + fsum = gptot.prior("fsum", self.X, reparameterize=False) + model2_logp = model2.compile_logp()({"fsum": self.y}) + npt.assert_allclose(model1_logp, model2_logp, atol=0, rtol=1e-2) + + with model1: + fp1 = gpsum.conditional("fp1", self.Xnew, given={"X": self.X, "f": self.y, "gp": gpsum}) + with model2: + fp2 = gptot.conditional("fp2", self.Xnew) + + fp = np.random.randn(self.Xnew.shape[0]) + logp1 = model1.compile_logp()({"fsum": self.y, "fp1": fp}) + logp2 = model2.compile_logp()({"fsum": self.y, "fp2": fp}) + npt.assert_allclose(logp1, logp2, atol=0, rtol=1e-2) + + def testAdditiveSparseRaises(self): + # cant add different approximations + with pm.Model() as model: + cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + gp1 = pm.gp.MarginalApprox(cov_func=cov_func, approx="DTC") + gp2 = pm.gp.MarginalApprox(cov_func=cov_func, approx="FITC") + with pytest.raises(Exception) as e_info: + gp1 + gp2 + + def testAdditiveTypeRaises1(self): + with pm.Model() as model: + cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + gp1 = pm.gp.MarginalApprox(cov_func=cov_func, approx="DTC") + gp2 = pm.gp.Marginal(cov_func=cov_func) + with pytest.raises(Exception) as e_info: + gp1 + gp2 + + def testAdditiveTypeRaises2(self): + with pm.Model() as model: + cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + gp1 = pm.gp.Latent(cov_func=cov_func) + gp2 = pm.gp.Marginal(cov_func=cov_func) + with pytest.raises(Exception) as e_info: + gp1 + gp2 + + +class TestMarginalVsLatent: + R""" + Compare the logp of models Marginal, noise=0 and Latent. + """ + + def setup_method(self): + X = np.random.randn(20, 3) + y = np.random.randn(20) + Xnew = np.random.randn(30, 3) + pnew = np.random.randn(30) + with pm.Model() as model: + cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + mean_func = pm.gp.mean.Constant(0.5) + gp = pm.gp.Marginal(mean_func=mean_func, cov_func=cov_func) + f = gp.marginal_likelihood("f", X, y, noise=0.0) + p = gp.conditional("p", Xnew) + self.logp = model.compile_logp()({"p": pnew}) + self.X = X + self.Xnew = Xnew + self.y = y + self.pnew = pnew + + def testLatent1(self): + with pm.Model() as model: + cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + mean_func = pm.gp.mean.Constant(0.5) + gp = pm.gp.Latent(mean_func=mean_func, cov_func=cov_func) + f = gp.prior("f", self.X, reparameterize=False) + p = gp.conditional("p", self.Xnew) + assert tuple(f.shape.eval()) == (self.X.shape[0],) + assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) + latent_logp = model.compile_logp()({"f": self.y, "p": self.pnew}) + npt.assert_allclose(latent_logp, self.logp, atol=0, rtol=1e-2) + + def testLatent2(self): + with pm.Model() as model: + cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + mean_func = pm.gp.mean.Constant(0.5) + gp = pm.gp.Latent(mean_func=mean_func, cov_func=cov_func) + f = gp.prior("f", self.X, reparameterize=True) + p = gp.conditional("p", self.Xnew) + assert tuple(f.shape.eval()) == (self.X.shape[0],) + assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) + chol = np.linalg.cholesky(cov_func(self.X).eval()) + y_rotated = np.linalg.solve(chol, self.y - 0.5) + latent_logp = model.compile_logp()({"f_rotated_": y_rotated, "p": self.pnew}) + npt.assert_allclose(latent_logp, self.logp, atol=5) + + +class TestTP: + R""" + Compare TP with high degress of freedom to GP + """ + + def setup_method(self): + X = np.random.randn(20, 3) + y = np.random.randn(20) + Xnew = np.random.randn(30, 3) + pnew = np.random.randn(30) + + with pm.Model() as model1: + cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + gp = pm.gp.Latent(cov_func=cov_func) + f = gp.prior("f", X, reparameterize=False) + p = gp.conditional("p", Xnew) + self.gp_latent_logp = model1.compile_logp()({"f": y, "p": pnew}) + self.X = X + self.y = y + self.Xnew = Xnew + self.pnew = pnew + self.nu = 10000 + + def testTPvsLatent(self): + with pm.Model() as model: + scale_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + tp = pm.gp.TP(scale_func=scale_func, nu=self.nu) + f = tp.prior("f", self.X, reparameterize=False) + p = tp.conditional("p", self.Xnew) + assert tuple(f.shape.eval()) == (self.X.shape[0],) + assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) + tp_logp = model.compile_logp()({"f": self.y, "p": self.pnew}) + npt.assert_allclose(self.gp_latent_logp, tp_logp, atol=0, rtol=1e-2) + + def testTPvsLatentReparameterized(self): + with pm.Model() as model: + scale_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + tp = pm.gp.TP(scale_func=scale_func, nu=self.nu) + f = tp.prior("f", self.X, reparameterize=True) + p = tp.conditional("p", self.Xnew) + assert tuple(f.shape.eval()) == (self.X.shape[0],) + assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) + chol = np.linalg.cholesky(scale_func(self.X).eval()) + f_rotated = np.linalg.solve(chol, self.y) + tp_logp = model.compile_logp()({"f_rotated_": f_rotated, "p": self.pnew}) + npt.assert_allclose(self.gp_latent_logp, tp_logp, atol=0, rtol=1e-2) + + def testAdditiveTPRaises(self): + with pm.Model() as model: + scale_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3]) + gp1 = pm.gp.TP(scale_func=scale_func, nu=10) + gp2 = pm.gp.TP(scale_func=scale_func, nu=10) + with pytest.raises(Exception) as e_info: + gp1 + gp2 + + +class TestLatentKron: + """ + Compare gp.LatentKron to gp.Latent, both with Gaussian noise. + """ + + def setup_method(self): + self.Xs = [ + np.linspace(0, 1, 7)[:, None], + np.linspace(0, 1, 5)[:, None], + np.linspace(0, 1, 6)[:, None], + ] + self.X = cartesian(*self.Xs) + self.N = np.prod([len(X) for X in self.Xs]) + self.y = np.random.randn(self.N) * 0.1 + self.Xnews = (np.random.randn(5, 1), np.random.randn(5, 1), np.random.randn(5, 1)) + self.Xnew = np.concatenate(self.Xnews, axis=1) + self.pnew = np.random.randn(len(self.Xnew)) + ls = 0.2 + with pm.Model() as latent_model: + self.cov_funcs = ( + pm.gp.cov.ExpQuad(1, ls), + pm.gp.cov.ExpQuad(1, ls), + pm.gp.cov.ExpQuad(1, ls), + ) + cov_func = pm.gp.cov.Kron(self.cov_funcs) + self.mean = pm.gp.mean.Constant(0.5) + gp = pm.gp.Latent(mean_func=self.mean, cov_func=cov_func) + f = gp.prior("f", self.X) + p = gp.conditional("p", self.Xnew) + chol = np.linalg.cholesky(cov_func(self.X).eval()) + self.y_rotated = np.linalg.solve(chol, self.y - 0.5) + self.logp = latent_model.compile_logp()({"f_rotated_": self.y_rotated, "p": self.pnew}) + + def testLatentKronvsLatent(self): + with pm.Model() as kron_model: + kron_gp = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + f = kron_gp.prior("f", self.Xs) + p = kron_gp.conditional("p", self.Xnew) + assert tuple(f.shape.eval()) == (self.X.shape[0],) + assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) + kronlatent_logp = kron_model.compile_logp()({"f_rotated_": self.y_rotated, "p": self.pnew}) + npt.assert_allclose(kronlatent_logp, self.logp, atol=0, rtol=1e-3) + + def testLatentKronRaisesAdditive(self): + with pm.Model() as kron_model: + gp1 = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + gp2 = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + with pytest.raises(TypeError): + gp1 + gp2 + + def testLatentKronRaisesSizes(self): + with pm.Model() as kron_model: + gp = pm.gp.LatentKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + with pytest.raises(ValueError): + gp.prior("f", Xs=[np.linspace(0, 1, 7)[:, None], np.linspace(0, 1, 5)[:, None]]) + + +class TestMarginalKron: + """ + Compare gp.MarginalKron to gp.Marginal. + """ + + def setup_method(self): + self.Xs = [ + np.linspace(0, 1, 7)[:, None], + np.linspace(0, 1, 5)[:, None], + np.linspace(0, 1, 6)[:, None], + ] + self.X = cartesian(*self.Xs) + self.N = np.prod([len(X) for X in self.Xs]) + self.y = np.random.randn(self.N) * 0.1 + self.Xnews = (np.random.randn(5, 1), np.random.randn(5, 1), np.random.randn(5, 1)) + self.Xnew = np.concatenate(self.Xnews, axis=1) + self.sigma = 0.2 + self.pnew = np.random.randn(len(self.Xnew)) + + ls = 0.2 + with pm.Model() as model: + self.cov_funcs = [ + pm.gp.cov.ExpQuad(1, ls), + pm.gp.cov.ExpQuad(1, ls), + pm.gp.cov.ExpQuad(1, ls), + ] + cov_func = pm.gp.cov.Kron(self.cov_funcs) + self.mean = pm.gp.mean.Constant(0.5) + gp = pm.gp.Marginal(mean_func=self.mean, cov_func=cov_func) + f = gp.marginal_likelihood("f", self.X, self.y, noise=self.sigma) + p = gp.conditional("p", self.Xnew) + self.mu, self.cov = gp.predict(self.Xnew) + self.logp = model.compile_logp()({"p": self.pnew}) + + def testMarginalKronvsMarginalpredict(self): + with pm.Model() as kron_model: + kron_gp = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + f = kron_gp.marginal_likelihood("f", self.Xs, self.y, sigma=self.sigma) + p = kron_gp.conditional("p", self.Xnew) + mu, cov = kron_gp.predict(self.Xnew) + assert tuple(f.shape.eval()) == (self.X.shape[0],) + assert tuple(p.shape.eval()) == (self.Xnew.shape[0],) + npt.assert_allclose(mu, self.mu, atol=1e-5, rtol=1e-2) + npt.assert_allclose(cov, self.cov, atol=1e-5, rtol=1e-2) + with kron_model: + _, var = kron_gp.predict(self.Xnew, diag=True) + npt.assert_allclose(np.diag(cov), var, atol=1e-5, rtol=1e-2) + + def testMarginalKronvsMarginal(self): + with pm.Model() as kron_model: + kron_gp = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + f = kron_gp.marginal_likelihood("f", self.Xs, self.y, sigma=self.sigma) + p = kron_gp.conditional("p", self.Xnew) + kron_logp = kron_model.compile_logp()({"p": self.pnew}) + npt.assert_allclose(kron_logp, self.logp, atol=0, rtol=1e-2) + + def testMarginalKronRaises(self): + with pm.Model() as kron_model: + gp1 = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + gp2 = pm.gp.MarginalKron(mean_func=self.mean, cov_funcs=self.cov_funcs) + with pytest.raises(TypeError): + gp1 + gp2 diff --git a/pymc/tests/gp/test_mean.py b/pymc/tests/gp/test_mean.py new file mode 100644 index 00000000000..eaceba6aca5 --- /dev/null +++ b/pymc/tests/gp/test_mean.py @@ -0,0 +1,90 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import numpy as np +import numpy.testing as npt + +import pymc as pm + + +class TestZeroMean: + def test_value(self): + X = np.linspace(0, 1, 10)[:, None] + with pm.Model() as model: + zero_mean = pm.gp.mean.Zero() + M = zero_mean(X).eval() + assert np.all(M == 0) + assert M.shape == (10,) + + +class TestConstantMean: + def test_value(self): + X = np.linspace(0, 1, 10)[:, None] + with pm.Model() as model: + const_mean = pm.gp.mean.Constant(6) + M = const_mean(X).eval() + assert np.all(M == 6) + assert M.shape == (10,) + + +class TestLinearMean: + def test_value(self): + X = np.linspace(0, 1, 10)[:, None] + with pm.Model() as model: + linear_mean = pm.gp.mean.Linear(2, 0.5) + M = linear_mean(X).eval() + npt.assert_allclose(M[1], 0.7222, atol=1e-3) + assert M.shape == (10,) + + +class TestAddProdMean: + def test_add(self): + X = np.linspace(0, 1, 10)[:, None] + with pm.Model() as model: + mean1 = pm.gp.mean.Linear(coeffs=2, intercept=0.5) + mean2 = pm.gp.mean.Constant(2) + mean = mean1 + mean2 + mean2 + M = mean(X).eval() + npt.assert_allclose(M[1], 0.7222 + 2 + 2, atol=1e-3) + + def test_prod(self): + X = np.linspace(0, 1, 10)[:, None] + with pm.Model() as model: + mean1 = pm.gp.mean.Linear(coeffs=2, intercept=0.5) + mean2 = pm.gp.mean.Constant(2) + mean = mean1 * mean2 * mean2 + M = mean(X).eval() + npt.assert_allclose(M[1], 0.7222 * 2 * 2, atol=1e-3) + + def test_add_multid(self): + X = np.linspace(0, 1, 30).reshape(10, 3) + A = np.array([1, 2, 3]) + b = 10 + with pm.Model() as model: + mean1 = pm.gp.mean.Linear(coeffs=A, intercept=b) + mean2 = pm.gp.mean.Constant(2) + mean = mean1 + mean2 + mean2 + M = mean(X).eval() + npt.assert_allclose(M[1], 10.8965 + 2 + 2, atol=1e-3) + + def test_prod_multid(self): + X = np.linspace(0, 1, 30).reshape(10, 3) + A = np.array([1, 2, 3]) + b = 10 + with pm.Model() as model: + mean1 = pm.gp.mean.Linear(coeffs=A, intercept=b) + mean2 = pm.gp.mean.Constant(2) + mean = mean1 * mean2 * mean2 + M = mean(X).eval() + npt.assert_allclose(M[1], 10.8965 * 2 * 2, atol=1e-3) diff --git a/pymc/tests/gp/test_util.py b/pymc/tests/gp/test_util.py new file mode 100644 index 00000000000..52568b6a09b --- /dev/null +++ b/pymc/tests/gp/test_util.py @@ -0,0 +1,102 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import aesara.tensor as at +import numpy as np +import numpy.testing as npt +import pytest + +import pymc as pm + + +class TestPlotGP: + def test_plot_gp_dist(self): + """Test that the plotting helper works with the stated input shapes.""" + import matplotlib.pyplot as plt + + X = 100 + S = 500 + fig, ax = plt.subplots() + pm.gp.util.plot_gp_dist( + ax, x=np.linspace(0, 50, X), samples=np.random.normal(np.arange(X), size=(S, X)) + ) + plt.close() + pass + + def test_plot_gp_dist_warn_nan(self): + """Test that the plotting helper works with the stated input shapes.""" + import matplotlib.pyplot as plt + + X = 100 + S = 500 + samples = np.random.normal(np.arange(X), size=(S, X)) + samples[15, 3] = np.nan + fig, ax = plt.subplots() + with pytest.warns(UserWarning): + pm.gp.util.plot_gp_dist(ax, x=np.linspace(0, 50, X), samples=samples) + plt.close() + pass + + +class TestKmeansInducing: + def setup_method(self): + self.centers = (-5, 5) + self.x = np.concatenate( + (self.centers[0] + np.random.randn(500), self.centers[1] + np.random.randn(500)) + ) + + def test_kmeans(self): + X = self.x[:, None] + Xu = pm.gp.util.kmeans_inducing_points(2, X).flatten() + npt.assert_allclose(np.asarray(self.centers), np.sort(Xu), atol=0.1) + + X = at.as_tensor_variable(self.x[:, None]) + Xu = pm.gp.util.kmeans_inducing_points(2, X).flatten() + npt.assert_allclose(np.asarray(self.centers), np.sort(Xu), atol=0.1) + + def test_kmeans_raises(self): + with pytest.raises(TypeError): + Xu = pm.gp.util.kmeans_inducing_points(2, "str is the wrong type").flatten() + + +class TestReplaceWithValues: + def test_basic_replace(self): + with pm.Model() as model: + a = pm.Normal("a") + b = pm.Normal("b", mu=a) + c = a * b + + (c_val,) = pm.gp.util.replace_with_values( + [c], replacements={"a": 2, "b": 3, "x": 100}, model=model + ) + assert c_val == np.array(6.0) + + def test_replace_no_inputs_needed(self): + with pm.Model() as model: + a = at.as_tensor_variable(2.0) + b = 1.0 + a + c = a * b + (c_val,) = pm.gp.util.replace_with_values([c], replacements={"x": 100}) + assert c_val == np.array(6.0) + + def test_missing_input(self): + with pm.Model() as model: + a = pm.Normal("a") + b = pm.Normal("b", mu=a) + c = a * b + + with pytest.raises(ValueError): + (c_val,) = pm.gp.util.replace_with_values( + [c], replacements={"a": 2, "x": 100}, model=model + ) From 9c680596aa80efc38461d0c808ebafca4a24d118 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Fri, 9 Sep 2022 11:12:01 -0400 Subject: [PATCH 02/15] Split up test_ode.py --- .github/workflows/tests.yml | 7 ++-- pymc/tests/ode/__init__.py | 0 pymc/tests/{ => ode}/test_ode.py | 40 ----------------------- pymc/tests/ode/test_utils.py | 56 ++++++++++++++++++++++++++++++++ 4 files changed, 60 insertions(+), 43 deletions(-) create mode 100644 pymc/tests/ode/__init__.py rename pymc/tests/{ => ode}/test_ode.py (92%) create mode 100644 pymc/tests/ode/test_utils.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 1e8b823284f..d6a49cc1945 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -80,7 +80,8 @@ jobs: pymc/tests/gp/test_util.py pymc/tests/test_model.py pymc/tests/test_model_graph.py - pymc/tests/test_ode.py + pymc/tests/ode/test_ode.py + pymc/tests/ode/test_utils.py pymc/tests/test_profile.py pymc/tests/test_quadpotential.py @@ -151,7 +152,7 @@ jobs: test-subset: - pymc/tests/test_variational_inference.py pymc/tests/test_initial_point.py - pymc/tests/test_pickling.py pymc/tests/test_profile.py pymc/tests/test_step.py - - pymc/tests/gp/test_cov.py pymc/tests/gp/test_gp.py pymc/tests/gp/test_mean.py pymc/tests/gp/test_util.py pymc/tests/test_ode.py pymc/tests/test_smc.py pymc/tests/test_parallel_sampling.py + - pymc/tests/gp/test_cov.py pymc/tests/gp/test_gp.py pymc/tests/gp/test_mean.py pymc/tests/gp/test_util.py pymc/tests/ode/test_ode.py pymc/tests/ode/test_utils.py pymc/tests/test_smc.py pymc/tests/test_parallel_sampling.py - pymc/tests/test_sampling.py pymc/tests/test_posteriors.py fail-fast: false @@ -364,7 +365,7 @@ jobs: floatx: [float32] python-version: ["3.10"] test-subset: - - pymc/tests/test_sampling.py pymc/tests/test_ode.py + - pymc/tests/test_sampling.py pymc/tests/ode/test_ode.py pymc/tests/ode/test_utils.py fail-fast: false runs-on: ${{ matrix.os }} env: diff --git a/pymc/tests/ode/__init__.py b/pymc/tests/ode/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pymc/tests/test_ode.py b/pymc/tests/ode/test_ode.py similarity index 92% rename from pymc/tests/test_ode.py rename to pymc/tests/ode/test_ode.py index 96bbf22ee7b..8f90c677e9f 100644 --- a/pymc/tests/test_ode.py +++ b/pymc/tests/ode/test_ode.py @@ -18,54 +18,14 @@ import numpy as np import pytest -from scipy.integrate import odeint from scipy.stats import norm import pymc as pm from pymc.ode import DifferentialEquation -from pymc.ode.utils import augment_system from pymc.tests.helpers import fast_unstable_sampling_mode -def test_gradients(): - """Tests the computation of the sensitivities from the Aesara computation graph""" - - # ODE system for which to compute gradients - def ode_func(y, t, p): - return np.exp(-t) - p[0] * y[0] - - # Computation of graidients with Aesara - augmented_ode_func = augment_system(ode_func, 1, 1 + 1) - - # This is the new system, ODE + Sensitivities, which will be integrated - def augmented_system(Y, t, p): - dydt, ddt_dydp = augmented_ode_func(Y[:1], t, p, Y[1:]) - derivatives = np.concatenate([dydt, ddt_dydp]) - return derivatives - - # Create real sensitivities - y0 = 0.0 - t = np.arange(0, 12, 0.25).reshape(-1, 1) - a = 0.472 - p = np.array([y0, a]) - - # Derivatives of the analytic solution with respect to y0 and alpha - # Treat y0 like a parameter and solve analytically. Then differentiate. - # I used CAS to get these derivatives - y0_sensitivity = np.exp(-a * t) - a_sensitivity = ( - -(np.exp(t * (a - 1)) - 1 + (a - 1) * (y0 * a - y0 - 1) * t) * np.exp(-a * t) / (a - 1) ** 2 - ) - - sensitivity = np.c_[y0_sensitivity, a_sensitivity] - - integrated_solutions = odeint(func=augmented_system, y0=[y0, 1, 0], t=t.ravel(), args=(p,)) - simulated_sensitivity = integrated_solutions[:, 1:] - - np.testing.assert_allclose(sensitivity, simulated_sensitivity, rtol=1e-5) - - def test_simulate(): """Tests the integration in DifferentialEquation""" diff --git a/pymc/tests/ode/test_utils.py b/pymc/tests/ode/test_utils.py new file mode 100644 index 00000000000..9faf60e3b60 --- /dev/null +++ b/pymc/tests/ode/test_utils.py @@ -0,0 +1,56 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import numpy as np +import scipy.integrate as ode + +from pymc.ode.utils import augment_system + + +def test_gradients(): + """Tests the computation of the sensitivities from the Aesara computation graph""" + + # ODE system for which to compute gradients + def ode_func(y, t, p): + return np.exp(-t) - p[0] * y[0] + + # Computation of graidients with Aesara + augmented_ode_func = augment_system(ode_func, 1, 1 + 1) + + # This is the new system, ODE + Sensitivities, which will be integrated + def augmented_system(Y, t, p): + dydt, ddt_dydp = augmented_ode_func(Y[:1], t, p, Y[1:]) + derivatives = np.concatenate([dydt, ddt_dydp]) + return derivatives + + # Create real sensitivities + y0 = 0.0 + t = np.arange(0, 12, 0.25).reshape(-1, 1) + a = 0.472 + p = np.array([y0, a]) + + # Derivatives of the analytic solution with respect to y0 and alpha + # Treat y0 like a parameter and solve analytically. Then differentiate. + # I used CAS to get these derivatives + y0_sensitivity = np.exp(-a * t) + a_sensitivity = ( + -(np.exp(t * (a - 1)) - 1 + (a - 1) * (y0 * a - y0 - 1) * t) * np.exp(-a * t) / (a - 1) ** 2 + ) + + sensitivity = np.c_[y0_sensitivity, a_sensitivity] + + integrated_solutions = ode.odeint(func=augmented_system, y0=[y0, 1, 0], t=t.ravel(), args=(p,)) + simulated_sensitivity = integrated_solutions[:, 1:] + + np.testing.assert_allclose(sensitivity, simulated_sensitivity, rtol=1e-5) From 2e2aa6c16d5415227c489b049c3957e2b9f7e560 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 20:35:57 -0400 Subject: [PATCH 03/15] Split up test_tuning.py --- .github/workflows/tests.yml | 3 +- pymc/tests/tuning/__init__.py | 0 pymc/tests/tuning/test_scaling.py | 33 +++++++++++++++++++ .../test_starting.py} | 21 +----------- 4 files changed, 36 insertions(+), 21 deletions(-) create mode 100644 pymc/tests/tuning/__init__.py create mode 100644 pymc/tests/tuning/test_scaling.py rename pymc/tests/{test_tuning.py => tuning/test_starting.py} (72%) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index d6a49cc1945..7743491b897 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -59,7 +59,8 @@ jobs: pymc/tests/distributions/test_simulator.py - | - pymc/tests/test_tuning.py + pymc/tests/tuning/test_scaling.py + pymc/tests/tuning/test_starting.py pymc/tests/test_shared.py pymc/tests/test_types.py diff --git a/pymc/tests/tuning/__init__.py b/pymc/tests/tuning/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pymc/tests/tuning/test_scaling.py b/pymc/tests/tuning/test_scaling.py new file mode 100644 index 00000000000..3306ab08165 --- /dev/null +++ b/pymc/tests/tuning/test_scaling.py @@ -0,0 +1,33 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import warnings + +import numpy as np + +from pymc.tests import models +from pymc.tuning import scaling + + +def test_adjust_precision(): + a = np.array([-10, -0.01, 0, 10, 1e300, -np.inf, np.inf]) + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", "divide by zero encountered in log", RuntimeWarning) + a1 = scaling.adjust_precision(a) + assert all((a1 > 0) & (a1 < 1e200)) + + +def test_guess_scaling(): + start, model, _ = models.non_normal(n=5) + a1 = scaling.guess_scaling(start, model=model) + assert all((a1 > 0) & (a1 < 1e200)) diff --git a/pymc/tests/test_tuning.py b/pymc/tests/tuning/test_starting.py similarity index 72% rename from pymc/tests/test_tuning.py rename to pymc/tests/tuning/test_starting.py index 7ed406ee8c5..4c98a940199 100644 --- a/pymc/tests/test_tuning.py +++ b/pymc/tests/tuning/test_starting.py @@ -11,30 +11,12 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -import warnings - import numpy as np import pytest -from numpy import inf - from pymc.step_methods.metropolis import tune from pymc.tests import models -from pymc.tuning import find_MAP, scaling - - -def test_adjust_precision(): - a = np.array([-10, -0.01, 0, 10, 1e300, -inf, inf]) - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", "divide by zero encountered in log", RuntimeWarning) - a1 = scaling.adjust_precision(a) - assert all((a1 > 0) & (a1 < 1e200)) - - -def test_guess_scaling(): - start, model, _ = models.non_normal(n=5) - a1 = scaling.guess_scaling(start, model=model) - assert all((a1 > 0) & (a1 < 1e200)) +from pymc.tuning import find_MAP @pytest.mark.parametrize("bounded", [False, True]) @@ -54,4 +36,3 @@ def test_tune_not_inplace(): returned_scaling = tune(orig_scaling, acc_rate=0.6) assert not returned_scaling is orig_scaling assert np.all(orig_scaling == np.array([0.001, 0.1])) - pass From 5ef94da772836afa7bd9014fe3e8ba18926bbda1 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 21:48:35 -0400 Subject: [PATCH 04/15] Move test_starting.py into tuning/ --- .github/workflows/tests.yml | 1 - pymc/tests/test_starting.py | 163 ----------------------------- pymc/tests/tuning/test_starting.py | 138 +++++++++++++++++++++++- 3 files changed, 137 insertions(+), 165 deletions(-) delete mode 100644 pymc/tests/test_starting.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 7743491b897..1c7c73740ba 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -46,7 +46,6 @@ jobs: pymc/tests/test_hmc.py pymc/tests/test_func_utils.py pymc/tests/distributions/test_shape_utils.py - pymc/tests/test_starting.py pymc/tests/distributions/test_mixture.py - | diff --git a/pymc/tests/test_starting.py b/pymc/tests/test_starting.py deleted file mode 100644 index 3c9ce07f98a..00000000000 --- a/pymc/tests/test_starting.py +++ /dev/null @@ -1,163 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import numpy as np -import pytest - -from pymc import ( - Beta, - Binomial, - Deterministic, - Gamma, - Model, - Normal, - Point, - Uniform, - find_MAP, -) -from pymc.exceptions import ImputationWarning -from pymc.tests.checks import close_to -from pymc.tests.helpers import select_by_precision -from pymc.tests.models import non_normal, simple_arbitrary_det, simple_model -from pymc.tuning import starting - - -def test_accuracy_normal(): - _, model, (mu, _) = simple_model() - with model: - newstart = find_MAP(Point(x=[-10.5, 100.5])) - close_to(newstart["x"], [mu, mu], select_by_precision(float64=1e-5, float32=1e-4)) - - -def test_accuracy_non_normal(): - _, model, (mu, _) = non_normal(4) - with model: - newstart = find_MAP(Point(x=[0.5, 0.01, 0.95, 0.99])) - close_to(newstart["x"], mu, select_by_precision(float64=1e-5, float32=1e-4)) - - -def test_find_MAP_discrete(): - tol1 = 2.0**-11 - tol2 = 2.0**-6 - alpha = 4 - beta = 4 - n = 20 - yes = 15 - - with Model() as model: - p = Beta("p", alpha, beta) - Binomial("ss", n=n, p=p) - Binomial("s", n=n, p=p, observed=yes) - - map_est1 = starting.find_MAP() - map_est2 = starting.find_MAP(vars=model.value_vars) - - close_to(map_est1["p"], 0.6086956533498806, tol1) - - close_to(map_est2["p"], 0.695642178810167, tol2) - assert map_est2["ss"] == 14 - - -def test_find_MAP_no_gradient(): - _, model = simple_arbitrary_det() - with model: - find_MAP() - - -def test_find_MAP(): - tol = 2.0**-11 # 16 bit machine epsilon, a low bar - data = np.random.randn(100) - # data should be roughly mean 0, std 1, but let's - # normalize anyway to get it really close - data = (data - np.mean(data)) / np.std(data) - - with Model(): - mu = Uniform("mu", -1, 1) - sigma = Uniform("sigma", 0.5, 1.5) - Normal("y", mu=mu, tau=sigma**-2, observed=data) - - # Test gradient minimization - map_est1 = starting.find_MAP(progressbar=False) - # Test non-gradient minimization - map_est2 = starting.find_MAP(progressbar=False, method="Powell") - - close_to(map_est1["mu"], 0, tol) - close_to(map_est1["sigma"], 1, tol) - - close_to(map_est2["mu"], 0, tol) - close_to(map_est2["sigma"], 1, tol) - - -def test_find_MAP_issue_5923(): - # Test that gradient-based minimization works well regardless of the order - # of variables in `vars`, and even when starting a reasonable distance from - # the MAP. - tol = 2.0**-11 # 16 bit machine epsilon, a low bar - data = np.random.randn(100) - # data should be roughly mean 0, std 1, but let's - # normalize anyway to get it really close - data = (data - np.mean(data)) / np.std(data) - - with Model(): - mu = Uniform("mu", -1, 1) - sigma = Uniform("sigma", 0.5, 1.5) - Normal("y", mu=mu, tau=sigma**-2, observed=data) - - start = {"mu": -0.5, "sigma": 1.25} - map_est1 = starting.find_MAP(progressbar=False, vars=[mu, sigma], start=start) - map_est2 = starting.find_MAP(progressbar=False, vars=[sigma, mu], start=start) - - close_to(map_est1["mu"], 0, tol) - close_to(map_est1["sigma"], 1, tol) - - close_to(map_est2["mu"], 0, tol) - close_to(map_est2["sigma"], 1, tol) - - -def test_find_MAP_issue_4488(): - # Test for https://github.com/pymc-devs/pymc/issues/4488 - with Model() as m: - with pytest.warns(ImputationWarning): - x = Gamma("x", alpha=3, beta=10, observed=np.array([1, np.nan])) - y = Deterministic("y", x + 1) - map_estimate = find_MAP() - - assert not set.difference({"x_missing", "x_missing_log__", "y"}, set(map_estimate.keys())) - np.testing.assert_allclose(map_estimate["x_missing"], 0.2, rtol=1e-4, atol=1e-4) - np.testing.assert_allclose(map_estimate["y"], [2.0, map_estimate["x_missing"][0] + 1]) - - -def test_allinmodel(): - model1 = Model() - model2 = Model() - with model1: - x1 = Normal("x1", mu=0, sigma=1) - y1 = Normal("y1", mu=0, sigma=1) - with model2: - x2 = Normal("x2", mu=0, sigma=1) - y2 = Normal("y2", mu=0, sigma=1) - - x1 = model1.rvs_to_values[x1] - y1 = model1.rvs_to_values[y1] - x2 = model2.rvs_to_values[x2] - y2 = model2.rvs_to_values[y2] - - starting.allinmodel([x1, y1], model1) - starting.allinmodel([x1], model1) - with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2', 'y2'\]"): - starting.allinmodel([x2, y2], model1) - with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2'\]"): - starting.allinmodel([x2, y1], model1) - with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2'\]"): - starting.allinmodel([x2], model1) diff --git a/pymc/tests/tuning/test_starting.py b/pymc/tests/tuning/test_starting.py index 4c98a940199..528e1fc3864 100644 --- a/pymc/tests/tuning/test_starting.py +++ b/pymc/tests/tuning/test_starting.py @@ -14,9 +14,15 @@ import numpy as np import pytest +import pymc as pm + +from pymc.exceptions import ImputationWarning from pymc.step_methods.metropolis import tune from pymc.tests import models -from pymc.tuning import find_MAP +from pymc.tests.checks import close_to +from pymc.tests.helpers import select_by_precision +from pymc.tests.models import non_normal, simple_arbitrary_det, simple_model +from pymc.tuning import find_MAP, starting @pytest.mark.parametrize("bounded", [False, True]) @@ -36,3 +42,133 @@ def test_tune_not_inplace(): returned_scaling = tune(orig_scaling, acc_rate=0.6) assert not returned_scaling is orig_scaling assert np.all(orig_scaling == np.array([0.001, 0.1])) + + +def test_accuracy_normal(): + _, model, (mu, _) = simple_model() + with model: + newstart = find_MAP(pm.Point(x=[-10.5, 100.5])) + close_to(newstart["x"], [mu, mu], select_by_precision(float64=1e-5, float32=1e-4)) + + +def test_accuracy_non_normal(): + _, model, (mu, _) = non_normal(4) + with model: + newstart = find_MAP(pm.Point(x=[0.5, 0.01, 0.95, 0.99])) + close_to(newstart["x"], mu, select_by_precision(float64=1e-5, float32=1e-4)) + + +def test_find_MAP_discrete(): + tol1 = 2.0**-11 + tol2 = 2.0**-6 + alpha = 4 + beta = 4 + n = 20 + yes = 15 + + with pm.Model() as model: + p = pm.Beta("p", alpha, beta) + pm.Binomial("ss", n=n, p=p) + pm.Binomial("s", n=n, p=p, observed=yes) + + map_est1 = find_MAP() + map_est2 = find_MAP(vars=model.value_vars) + + close_to(map_est1["p"], 0.6086956533498806, tol1) + + close_to(map_est2["p"], 0.695642178810167, tol2) + assert map_est2["ss"] == 14 + + +def test_find_MAP_no_gradient(): + _, model = simple_arbitrary_det() + with model: + find_MAP() + + +def test_find_MAP(): + tol = 2.0**-11 # 16 bit machine epsilon, a low bar + data = np.random.randn(100) + # data should be roughly mean 0, std 1, but let's + # normalize anyway to get it really close + data = (data - np.mean(data)) / np.std(data) + + with pm.Model(): + mu = pm.Uniform("mu", -1, 1) + sigma = pm.Uniform("sigma", 0.5, 1.5) + pm.Normal("y", mu=mu, tau=sigma**-2, observed=data) + + # Test gradient minimization + map_est1 = find_MAP(progressbar=False) + # Test non-gradient minimization + map_est2 = find_MAP(progressbar=False, method="Powell") + + close_to(map_est1["mu"], 0, tol) + close_to(map_est1["sigma"], 1, tol) + + close_to(map_est2["mu"], 0, tol) + close_to(map_est2["sigma"], 1, tol) + + +def test_find_MAP_issue_5923(): + # Test that gradient-based minimization works well regardless of the order + # of variables in `vars`, and even when starting a reasonable distance from + # the MAP. + tol = 2.0**-11 # 16 bit machine epsilon, a low bar + data = np.random.randn(100) + # data should be roughly mean 0, std 1, but let's + # normalize anyway to get it really close + data = (data - np.mean(data)) / np.std(data) + + with pm.Model(): + mu = pm.Uniform("mu", -1, 1) + sigma = pm.Uniform("sigma", 0.5, 1.5) + pm.Normal("y", mu=mu, tau=sigma**-2, observed=data) + + start = {"mu": -0.5, "sigma": 1.25} + map_est1 = find_MAP(progressbar=False, vars=[mu, sigma], start=start) + map_est2 = find_MAP(progressbar=False, vars=[sigma, mu], start=start) + + close_to(map_est1["mu"], 0, tol) + close_to(map_est1["sigma"], 1, tol) + + close_to(map_est2["mu"], 0, tol) + close_to(map_est2["sigma"], 1, tol) + + +def test_find_MAP_issue_4488(): + # Test for https://github.com/pymc-devs/pymc/issues/4488 + with pm.Model() as m: + with pytest.warns(ImputationWarning): + x = pm.Gamma("x", alpha=3, beta=10, observed=np.array([1, np.nan])) + y = pm.Deterministic("y", x + 1) + map_estimate = find_MAP() + + assert not set.difference({"x_missing", "x_missing_log__", "y"}, set(map_estimate.keys())) + np.testing.assert_allclose(map_estimate["x_missing"], 0.2, rtol=1e-4, atol=1e-4) + np.testing.assert_allclose(map_estimate["y"], [2.0, map_estimate["x_missing"][0] + 1]) + + +def test_allinmodel(): + model1 = pm.Model() + model2 = pm.Model() + with model1: + x1 = pm.Normal("x1", mu=0, sigma=1) + y1 = pm.Normal("y1", mu=0, sigma=1) + with model2: + x2 = pm.Normal("x2", mu=0, sigma=1) + y2 = pm.Normal("y2", mu=0, sigma=1) + + x1 = model1.rvs_to_values[x1] + y1 = model1.rvs_to_values[y1] + x2 = model2.rvs_to_values[x2] + y2 = model2.rvs_to_values[y2] + + starting.allinmodel([x1, y1], model1) + starting.allinmodel([x1], model1) + with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2', 'y2'\]"): + starting.allinmodel([x2, y2], model1) + with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2'\]"): + starting.allinmodel([x2, y1], model1) + with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2'\]"): + starting.allinmodel([x2], model1) From 36912a399d39f25b577790910f7509ebee2b6614 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 20:39:44 -0400 Subject: [PATCH 05/15] Replace test_ndarray_backend.py --- .github/workflows/tests.yml | 2 +- pymc/tests/backends/__init__.py | 0 pymc/tests/{backend_fixtures.py => backends/fixtures.py} | 0 .../tests/{test_ndarray_backend.py => backends/test_ndarray.py} | 2 +- 4 files changed, 2 insertions(+), 2 deletions(-) create mode 100644 pymc/tests/backends/__init__.py rename pymc/tests/{backend_fixtures.py => backends/fixtures.py} (100%) rename pymc/tests/{test_ndarray_backend.py => backends/test_ndarray.py} (99%) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 1c7c73740ba..843135b10cd 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -42,7 +42,7 @@ jobs: pymc/tests/test_aesaraf.py pymc/tests/test_math.py pymc/tests/test_posdef_sym.py - pymc/tests/test_ndarray_backend.py + pymc/tests/backends/test_ndarray.py pymc/tests/test_hmc.py pymc/tests/test_func_utils.py pymc/tests/distributions/test_shape_utils.py diff --git a/pymc/tests/backends/__init__.py b/pymc/tests/backends/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pymc/tests/backend_fixtures.py b/pymc/tests/backends/fixtures.py similarity index 100% rename from pymc/tests/backend_fixtures.py rename to pymc/tests/backends/fixtures.py diff --git a/pymc/tests/test_ndarray_backend.py b/pymc/tests/backends/test_ndarray.py similarity index 99% rename from pymc/tests/test_ndarray_backend.py rename to pymc/tests/backends/test_ndarray.py index a2aa231a437..9abbfe5fbe3 100644 --- a/pymc/tests/test_ndarray_backend.py +++ b/pymc/tests/backends/test_ndarray.py @@ -17,7 +17,7 @@ import pytest from pymc.backends import base, ndarray -from pymc.tests import backend_fixtures as bf +from pymc.tests.backends import fixtures as bf STATS1 = [{"a": np.float64, "b": bool}] From dbd67ac6bacda19c11b3f014d161c207af6491e2 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 20:43:04 -0400 Subject: [PATCH 06/15] Rename test_data_container.py into test_data.py --- .github/workflows/tests.yml | 2 +- pymc/tests/{test_data_container.py => test_data.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename pymc/tests/{test_data_container.py => test_data.py} (100%) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 843135b10cd..00e0ad36e78 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -224,7 +224,7 @@ jobs: test-subset: - | pymc/tests/test_parallel_sampling.py - pymc/tests/test_data_container.py + pymc/tests/test_data.py pymc/tests/test_missing.py - | diff --git a/pymc/tests/test_data_container.py b/pymc/tests/test_data.py similarity index 100% rename from pymc/tests/test_data_container.py rename to pymc/tests/test_data.py From ab6a27058c0e93db7f4737d15a1c9f2ea6f7c15d Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 20:47:07 -0400 Subject: [PATCH 07/15] Move test_idata_conversion.py into backends --- .github/workflows/tests.yml | 2 +- pymc/tests/{test_idata_conversion.py => backends/test_arviz.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename pymc/tests/{test_idata_conversion.py => backends/test_arviz.py} (100%) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 00e0ad36e78..e264475c103 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -231,7 +231,7 @@ jobs: pymc/tests/test_sampling.py - | - pymc/tests/test_idata_conversion.py + pymc/tests/backends/test_arviz.py pymc/tests/test_updates.py fail-fast: false runs-on: ${{ matrix.os }} diff --git a/pymc/tests/test_idata_conversion.py b/pymc/tests/backends/test_arviz.py similarity index 100% rename from pymc/tests/test_idata_conversion.py rename to pymc/tests/backends/test_arviz.py From 83efffcf195989ac972c0f5312969323eed08f36 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 20:56:32 -0400 Subject: [PATCH 08/15] Merge test_minibatches.py into test_data.py --- .github/workflows/tests.yml | 1 - pymc/tests/test_data.py | 344 ++++++++++++++++++++++++++++++- pymc/tests/test_minibatches.py | 357 --------------------------------- 3 files changed, 338 insertions(+), 364 deletions(-) delete mode 100644 pymc/tests/test_minibatches.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index e264475c103..9a981c97ea4 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -66,7 +66,6 @@ jobs: - | pymc/tests/test_modelcontext.py pymc/tests/distributions/test_dist_math.py - pymc/tests/test_minibatches.py pymc/tests/test_pickling.py pymc/tests/distributions/test_transform.py pymc/tests/test_parallel_sampling.py diff --git a/pymc/tests/test_data.py b/pymc/tests/test_data.py index 4ccd86426b2..a6f1dd36dfe 100644 --- a/pymc/tests/test_data.py +++ b/pymc/tests/test_data.py @@ -13,19 +13,23 @@ # limitations under the License. import io +import itertools as it +import aesara +import aesara.tensor as at +import cloudpickle import numpy as np import pytest +import scipy.stats as st from aesara import shared -from aesara.tensor import TensorConstant from aesara.tensor.var import TensorVariable import pymc as pm -from pymc.aesaraf import floatX +from pymc.aesaraf import GeneratorOp, floatX from pymc.exceptions import ShapeError -from pymc.tests.helpers import SeededTest +from pymc.tests.helpers import SeededTest, select_by_precision class TestData(SeededTest): @@ -315,8 +319,8 @@ def test_explicit_coords(self): # pass coordinates explicitly, use numpy array in Data container with pm.Model(coords=coords) as pmodel: # Dims created from coords are constant by default - assert isinstance(pmodel.dim_lengths["rows"], TensorConstant) - assert isinstance(pmodel.dim_lengths["columns"], TensorConstant) + assert isinstance(pmodel.dim_lengths["rows"], at.TensorConstant) + assert isinstance(pmodel.dim_lengths["columns"], at.TensorConstant) pm.MutableData("observations", data, dims=("rows", "columns")) # new data with same (!) shape pm.set_data({"observations": data + 1}) @@ -432,7 +436,7 @@ def test_data_mutable_default_warning(self): with pm.Model(): with pytest.warns(UserWarning, match="`mutable` kwarg was not specified"): data = pm.Data("x", [1, 2, 3]) - assert isinstance(data, TensorConstant) + assert isinstance(data, at.TensorConstant) pass @@ -451,3 +455,331 @@ def test_data_naming(): def test_get_data(): data = pm.get_data("radon.csv") assert type(data) == io.BytesIO + + +class _DataSampler: + """ + Not for users + """ + + def __init__(self, data, batchsize=50, random_seed=42, dtype="floatX"): + self.dtype = aesara.config.floatX if dtype == "floatX" else dtype + self.rng = np.random.RandomState(random_seed) + self.data = data + self.n = batchsize + + def __iter__(self): + return self + + def __next__(self): + idx = self.rng.uniform(size=self.n, low=0.0, high=self.data.shape[0] - 1e-16).astype( + "int64" + ) + return np.asarray(self.data[idx], self.dtype) + + next = __next__ + + +@pytest.fixture(scope="module") +def datagen(): + return _DataSampler(np.random.uniform(size=(1000, 10))) + + +def integers(): + i = 0 + while True: + yield pm.floatX(i) + i += 1 + + +def integers_ndim(ndim): + i = 0 + while True: + yield np.ones((2,) * ndim) * i + i += 1 + + +@pytest.mark.usefixtures("strict_float32") +class TestGenerator: + def test_basic(self): + generator = pm.GeneratorAdapter(integers()) + gop = GeneratorOp(generator)() + assert gop.tag.test_value == np.float32(0) + f = aesara.function([], gop) + assert f() == np.float32(0) + assert f() == np.float32(1) + for _ in range(2, 100): + f() + assert f() == np.float32(100) + + def test_ndim(self): + for ndim in range(10): + res = list(it.islice(integers_ndim(ndim), 0, 2)) + generator = pm.GeneratorAdapter(integers_ndim(ndim)) + gop = GeneratorOp(generator)() + f = aesara.function([], gop) + assert ndim == res[0].ndim + np.testing.assert_equal(f(), res[0]) + np.testing.assert_equal(f(), res[1]) + + def test_cloning_available(self): + gop = pm.generator(integers()) + res = gop**2 + shared = aesara.shared(pm.floatX(10)) + res1 = aesara.clone_replace(res, {gop: shared}) + f = aesara.function([], res1) + assert f() == np.float32(100) + + def test_default_value(self): + def gen(): + for i in range(2): + yield pm.floatX(np.ones((10, 10)) * i) + + gop = pm.generator(gen(), np.ones((10, 10)) * 10) + f = aesara.function([], gop) + np.testing.assert_equal(np.ones((10, 10)) * 0, f()) + np.testing.assert_equal(np.ones((10, 10)) * 1, f()) + np.testing.assert_equal(np.ones((10, 10)) * 10, f()) + with pytest.raises(ValueError): + gop.set_default(1) + + def test_set_gen_and_exc(self): + def gen(): + for i in range(2): + yield pm.floatX(np.ones((10, 10)) * i) + + gop = pm.generator(gen()) + f = aesara.function([], gop) + np.testing.assert_equal(np.ones((10, 10)) * 0, f()) + np.testing.assert_equal(np.ones((10, 10)) * 1, f()) + with pytest.raises(StopIteration): + f() + gop.set_gen(gen()) + np.testing.assert_equal(np.ones((10, 10)) * 0, f()) + np.testing.assert_equal(np.ones((10, 10)) * 1, f()) + + def test_pickling(self, datagen): + gen = pm.generator(datagen) + cloudpickle.loads(cloudpickle.dumps(gen)) + bad_gen = pm.generator(integers()) + with pytest.raises(TypeError): + cloudpickle.dumps(bad_gen) + + def test_gen_cloning_with_shape_change(self, datagen): + gen = pm.generator(datagen) + gen_r = pm.at_rng().normal(size=gen.shape).T + X = gen.dot(gen_r) + res, _ = aesara.scan(lambda x: x.sum(), X, n_steps=X.shape[0]) + assert res.eval().shape == (50,) + shared = aesara.shared(datagen.data.astype(gen.dtype)) + res2 = aesara.clone_replace(res, {gen: shared**2}) + assert res2.eval().shape == (1000,) + + +def gen1(): + i = 0 + while True: + yield np.ones((10, 100)) * i + i += 1 + + +def gen2(): + i = 0 + while True: + yield np.ones((20, 100)) * i + i += 1 + + +class TestScaling: + """ + Related to minibatch training + """ + + def test_density_scaling(self): + with pm.Model() as model1: + pm.Normal("n", observed=[[1]], total_size=1) + p1 = aesara.function([], model1.logp()) + + with pm.Model() as model2: + pm.Normal("n", observed=[[1]], total_size=2) + p2 = aesara.function([], model2.logp()) + assert p1() * 2 == p2() + + def test_density_scaling_with_generator(self): + # We have different size generators + + def true_dens(): + g = gen1() + for i, point in enumerate(g): + yield st.norm.logpdf(point).sum() * 10 + + t = true_dens() + # We have same size models + with pm.Model() as model1: + pm.Normal("n", observed=gen1(), total_size=100) + p1 = aesara.function([], model1.logp()) + + with pm.Model() as model2: + gen_var = pm.generator(gen2()) + pm.Normal("n", observed=gen_var, total_size=100) + p2 = aesara.function([], model2.logp()) + + for i in range(10): + _1, _2, _t = p1(), p2(), next(t) + decimals = select_by_precision(float64=7, float32=1) + np.testing.assert_almost_equal(_1, _t, decimal=decimals) # Value O(-50,000) + np.testing.assert_almost_equal(_1, _2) + # Done + + def test_gradient_with_scaling(self): + with pm.Model() as model1: + genvar = pm.generator(gen1()) + m = pm.Normal("m") + pm.Normal("n", observed=genvar, total_size=1000) + grad1 = aesara.function([m.tag.value_var], at.grad(model1.logp(), m.tag.value_var)) + with pm.Model() as model2: + m = pm.Normal("m") + shavar = aesara.shared(np.ones((1000, 100))) + pm.Normal("n", observed=shavar) + grad2 = aesara.function([m.tag.value_var], at.grad(model2.logp(), m.tag.value_var)) + + for i in range(10): + shavar.set_value(np.ones((100, 100)) * i) + g1 = grad1(1) + g2 = grad2(1) + np.testing.assert_almost_equal(g1, g2) + + def test_multidim_scaling(self): + with pm.Model() as model0: + pm.Normal("n", observed=[[1, 1], [1, 1]], total_size=[]) + p0 = aesara.function([], model0.logp()) + + with pm.Model() as model1: + pm.Normal("n", observed=[[1, 1], [1, 1]], total_size=[2, 2]) + p1 = aesara.function([], model1.logp()) + + with pm.Model() as model2: + pm.Normal("n", observed=[[1], [1]], total_size=[2, 2]) + p2 = aesara.function([], model2.logp()) + + with pm.Model() as model3: + pm.Normal("n", observed=[[1, 1]], total_size=[2, 2]) + p3 = aesara.function([], model3.logp()) + + with pm.Model() as model4: + pm.Normal("n", observed=[[1]], total_size=[2, 2]) + p4 = aesara.function([], model4.logp()) + + with pm.Model() as model5: + pm.Normal("n", observed=[[1]], total_size=[2, Ellipsis, 2]) + p5 = aesara.function([], model5.logp()) + _p0 = p0() + assert ( + np.allclose(_p0, p1()) + and np.allclose(_p0, p2()) + and np.allclose(_p0, p3()) + and np.allclose(_p0, p4()) + and np.allclose(_p0, p5()) + ) + + def test_common_errors(self): + with pytest.raises(ValueError) as e: + with pm.Model() as m: + pm.Normal("n", observed=[[1]], total_size=[2, Ellipsis, 2, 2]) + m.logp() + assert "Length of" in str(e.value) + with pytest.raises(ValueError) as e: + with pm.Model() as m: + pm.Normal("n", observed=[[1]], total_size=[2, 2, 2]) + m.logp() + assert "Length of" in str(e.value) + with pytest.raises(TypeError) as e: + with pm.Model() as m: + pm.Normal("n", observed=[[1]], total_size="foo") + m.logp() + assert "Unrecognized" in str(e.value) + with pytest.raises(TypeError) as e: + with pm.Model() as m: + pm.Normal("n", observed=[[1]], total_size=["foo"]) + m.logp() + assert "Unrecognized" in str(e.value) + with pytest.raises(ValueError) as e: + with pm.Model() as m: + pm.Normal("n", observed=[[1]], total_size=[Ellipsis, Ellipsis]) + m.logp() + assert "Double Ellipsis" in str(e.value) + + def test_mixed1(self): + with pm.Model(): + data = np.random.rand(10, 20, 30, 40, 50) + mb = pm.Minibatch(data, [2, None, 20, Ellipsis, 10]) + pm.Normal("n", observed=mb, total_size=(10, None, 30, Ellipsis, 50)) + + def test_mixed2(self): + with pm.Model(): + data = np.random.rand(10, 20, 30, 40, 50) + mb = pm.Minibatch(data, [2, None, 20]) + pm.Normal("n", observed=mb, total_size=(10, None, 30)) + + def test_free_rv(self): + with pm.Model() as model4: + pm.Normal("n", observed=[[1, 1], [1, 1]], total_size=[2, 2]) + p4 = aesara.function([], model4.logp()) + + with pm.Model() as model5: + n = pm.Normal("n", total_size=[2, Ellipsis, 2], size=(2, 2)) + p5 = aesara.function([n.tag.value_var], model5.logp()) + assert p4() == p5(pm.floatX([[1]])) + assert p4() == p5(pm.floatX([[1, 1], [1, 1]])) + + +@pytest.mark.usefixtures("strict_float32") +class TestMinibatch: + data = np.random.rand(30, 10, 40, 10, 50) + + def test_1d(self): + mb = pm.Minibatch(self.data, 20) + assert mb.eval().shape == (20, 10, 40, 10, 50) + + def test_2d(self): + mb = pm.Minibatch(self.data, [(10, 42), (4, 42)]) + assert mb.eval().shape == (10, 4, 40, 10, 50) + + @pytest.mark.parametrize( + "batch_size, expected", + [ + ([(10, 42), None, (4, 42)], (10, 10, 4, 10, 50)), + ([(10, 42), Ellipsis, (4, 42)], (10, 10, 40, 10, 4)), + ([(10, 42), None, Ellipsis, (4, 42)], (10, 10, 40, 10, 4)), + ([10, None, Ellipsis, (4, 42)], (10, 10, 40, 10, 4)), + ], + ) + def test_special_batch_size(self, batch_size, expected): + mb = pm.Minibatch(self.data, batch_size) + assert mb.eval().shape == expected + + def test_cloning_available(self): + gop = pm.Minibatch(np.arange(100), 1) + res = gop**2 + shared = aesara.shared(np.array([10])) + res1 = aesara.clone_replace(res, {gop: shared}) + f = aesara.function([], res1) + assert f() == np.array([100]) + + def test_align(self): + m = pm.Minibatch(np.arange(1000), 1, random_seed=1) + n = pm.Minibatch(np.arange(1000), 1, random_seed=1) + f = aesara.function([], [m, n]) + n.eval() # not aligned + a, b = zip(*(f() for _ in range(1000))) + assert a != b + pm.align_minibatches() + a, b = zip(*(f() for _ in range(1000))) + assert a == b + n.eval() # not aligned + pm.align_minibatches([m]) + a, b = zip(*(f() for _ in range(1000))) + assert a != b + pm.align_minibatches([m, n]) + a, b = zip(*(f() for _ in range(1000))) + assert a == b diff --git a/pymc/tests/test_minibatches.py b/pymc/tests/test_minibatches.py deleted file mode 100644 index 58bf66c1386..00000000000 --- a/pymc/tests/test_minibatches.py +++ /dev/null @@ -1,357 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import itertools - -import aesara -import cloudpickle -import numpy as np -import pytest - -from aesara import tensor as at -from scipy import stats as stats - -import pymc as pm - -from pymc import GeneratorAdapter, Normal, at_rng, floatX, generator -from pymc.aesaraf import GeneratorOp -from pymc.tests.helpers import select_by_precision - - -class _DataSampler: - """ - Not for users - """ - - def __init__(self, data, batchsize=50, random_seed=42, dtype="floatX"): - self.dtype = aesara.config.floatX if dtype == "floatX" else dtype - self.rng = np.random.RandomState(random_seed) - self.data = data - self.n = batchsize - - def __iter__(self): - return self - - def __next__(self): - idx = self.rng.uniform(size=self.n, low=0.0, high=self.data.shape[0] - 1e-16).astype( - "int64" - ) - return np.asarray(self.data[idx], self.dtype) - - next = __next__ - - -@pytest.fixture(scope="module") -def datagen(): - return _DataSampler(np.random.uniform(size=(1000, 10))) - - -def integers(): - i = 0 - while True: - yield pm.floatX(i) - i += 1 - - -def integers_ndim(ndim): - i = 0 - while True: - yield np.ones((2,) * ndim) * i - i += 1 - - -@pytest.mark.usefixtures("strict_float32") -class TestGenerator: - def test_basic(self): - generator = GeneratorAdapter(integers()) - gop = GeneratorOp(generator)() - assert gop.tag.test_value == np.float32(0) - f = aesara.function([], gop) - assert f() == np.float32(0) - assert f() == np.float32(1) - for _ in range(2, 100): - f() - assert f() == np.float32(100) - - def test_ndim(self): - for ndim in range(10): - res = list(itertools.islice(integers_ndim(ndim), 0, 2)) - generator = GeneratorAdapter(integers_ndim(ndim)) - gop = GeneratorOp(generator)() - f = aesara.function([], gop) - assert ndim == res[0].ndim - np.testing.assert_equal(f(), res[0]) - np.testing.assert_equal(f(), res[1]) - - def test_cloning_available(self): - gop = generator(integers()) - res = gop**2 - shared = aesara.shared(floatX(10)) - res1 = aesara.clone_replace(res, {gop: shared}) - f = aesara.function([], res1) - assert f() == np.float32(100) - - def test_default_value(self): - def gen(): - for i in range(2): - yield floatX(np.ones((10, 10)) * i) - - gop = generator(gen(), np.ones((10, 10)) * 10) - f = aesara.function([], gop) - np.testing.assert_equal(np.ones((10, 10)) * 0, f()) - np.testing.assert_equal(np.ones((10, 10)) * 1, f()) - np.testing.assert_equal(np.ones((10, 10)) * 10, f()) - with pytest.raises(ValueError): - gop.set_default(1) - - def test_set_gen_and_exc(self): - def gen(): - for i in range(2): - yield floatX(np.ones((10, 10)) * i) - - gop = generator(gen()) - f = aesara.function([], gop) - np.testing.assert_equal(np.ones((10, 10)) * 0, f()) - np.testing.assert_equal(np.ones((10, 10)) * 1, f()) - with pytest.raises(StopIteration): - f() - gop.set_gen(gen()) - np.testing.assert_equal(np.ones((10, 10)) * 0, f()) - np.testing.assert_equal(np.ones((10, 10)) * 1, f()) - - def test_pickling(self, datagen): - gen = generator(datagen) - cloudpickle.loads(cloudpickle.dumps(gen)) - bad_gen = generator(integers()) - with pytest.raises(TypeError): - cloudpickle.dumps(bad_gen) - - def test_gen_cloning_with_shape_change(self, datagen): - gen = generator(datagen) - gen_r = at_rng().normal(size=gen.shape).T - X = gen.dot(gen_r) - res, _ = aesara.scan(lambda x: x.sum(), X, n_steps=X.shape[0]) - assert res.eval().shape == (50,) - shared = aesara.shared(datagen.data.astype(gen.dtype)) - res2 = aesara.clone_replace(res, {gen: shared**2}) - assert res2.eval().shape == (1000,) - - -def gen1(): - i = 0 - while True: - yield np.ones((10, 100)) * i - i += 1 - - -def gen2(): - i = 0 - while True: - yield np.ones((20, 100)) * i - i += 1 - - -class TestScaling: - """ - Related to minibatch training - """ - - def test_density_scaling(self): - with pm.Model() as model1: - Normal("n", observed=[[1]], total_size=1) - p1 = aesara.function([], model1.logp()) - - with pm.Model() as model2: - Normal("n", observed=[[1]], total_size=2) - p2 = aesara.function([], model2.logp()) - assert p1() * 2 == p2() - - def test_density_scaling_with_generator(self): - # We have different size generators - - def true_dens(): - g = gen1() - for i, point in enumerate(g): - yield stats.norm.logpdf(point).sum() * 10 - - t = true_dens() - # We have same size models - with pm.Model() as model1: - Normal("n", observed=gen1(), total_size=100) - p1 = aesara.function([], model1.logp()) - - with pm.Model() as model2: - gen_var = generator(gen2()) - Normal("n", observed=gen_var, total_size=100) - p2 = aesara.function([], model2.logp()) - - for i in range(10): - _1, _2, _t = p1(), p2(), next(t) - decimals = select_by_precision(float64=7, float32=1) - np.testing.assert_almost_equal(_1, _t, decimal=decimals) # Value O(-50,000) - np.testing.assert_almost_equal(_1, _2) - # Done - - def test_gradient_with_scaling(self): - with pm.Model() as model1: - genvar = generator(gen1()) - m = Normal("m") - Normal("n", observed=genvar, total_size=1000) - grad1 = aesara.function([m.tag.value_var], at.grad(model1.logp(), m.tag.value_var)) - with pm.Model() as model2: - m = Normal("m") - shavar = aesara.shared(np.ones((1000, 100))) - Normal("n", observed=shavar) - grad2 = aesara.function([m.tag.value_var], at.grad(model2.logp(), m.tag.value_var)) - - for i in range(10): - shavar.set_value(np.ones((100, 100)) * i) - g1 = grad1(1) - g2 = grad2(1) - np.testing.assert_almost_equal(g1, g2) - - def test_multidim_scaling(self): - with pm.Model() as model0: - Normal("n", observed=[[1, 1], [1, 1]], total_size=[]) - p0 = aesara.function([], model0.logp()) - - with pm.Model() as model1: - Normal("n", observed=[[1, 1], [1, 1]], total_size=[2, 2]) - p1 = aesara.function([], model1.logp()) - - with pm.Model() as model2: - Normal("n", observed=[[1], [1]], total_size=[2, 2]) - p2 = aesara.function([], model2.logp()) - - with pm.Model() as model3: - Normal("n", observed=[[1, 1]], total_size=[2, 2]) - p3 = aesara.function([], model3.logp()) - - with pm.Model() as model4: - Normal("n", observed=[[1]], total_size=[2, 2]) - p4 = aesara.function([], model4.logp()) - - with pm.Model() as model5: - Normal("n", observed=[[1]], total_size=[2, Ellipsis, 2]) - p5 = aesara.function([], model5.logp()) - _p0 = p0() - assert ( - np.allclose(_p0, p1()) - and np.allclose(_p0, p2()) - and np.allclose(_p0, p3()) - and np.allclose(_p0, p4()) - and np.allclose(_p0, p5()) - ) - - def test_common_errors(self): - with pytest.raises(ValueError) as e: - with pm.Model() as m: - Normal("n", observed=[[1]], total_size=[2, Ellipsis, 2, 2]) - m.logp() - assert "Length of" in str(e.value) - with pytest.raises(ValueError) as e: - with pm.Model() as m: - Normal("n", observed=[[1]], total_size=[2, 2, 2]) - m.logp() - assert "Length of" in str(e.value) - with pytest.raises(TypeError) as e: - with pm.Model() as m: - Normal("n", observed=[[1]], total_size="foo") - m.logp() - assert "Unrecognized" in str(e.value) - with pytest.raises(TypeError) as e: - with pm.Model() as m: - Normal("n", observed=[[1]], total_size=["foo"]) - m.logp() - assert "Unrecognized" in str(e.value) - with pytest.raises(ValueError) as e: - with pm.Model() as m: - Normal("n", observed=[[1]], total_size=[Ellipsis, Ellipsis]) - m.logp() - assert "Double Ellipsis" in str(e.value) - - def test_mixed1(self): - with pm.Model(): - data = np.random.rand(10, 20, 30, 40, 50) - mb = pm.Minibatch(data, [2, None, 20, Ellipsis, 10]) - Normal("n", observed=mb, total_size=(10, None, 30, Ellipsis, 50)) - - def test_mixed2(self): - with pm.Model(): - data = np.random.rand(10, 20, 30, 40, 50) - mb = pm.Minibatch(data, [2, None, 20]) - Normal("n", observed=mb, total_size=(10, None, 30)) - - def test_free_rv(self): - with pm.Model() as model4: - Normal("n", observed=[[1, 1], [1, 1]], total_size=[2, 2]) - p4 = aesara.function([], model4.logp()) - - with pm.Model() as model5: - n = Normal("n", total_size=[2, Ellipsis, 2], size=(2, 2)) - p5 = aesara.function([n.tag.value_var], model5.logp()) - assert p4() == p5(pm.floatX([[1]])) - assert p4() == p5(pm.floatX([[1, 1], [1, 1]])) - - -@pytest.mark.usefixtures("strict_float32") -class TestMinibatch: - data = np.random.rand(30, 10, 40, 10, 50) - - def test_1d(self): - mb = pm.Minibatch(self.data, 20) - assert mb.eval().shape == (20, 10, 40, 10, 50) - - def test_2d(self): - mb = pm.Minibatch(self.data, [(10, 42), (4, 42)]) - assert mb.eval().shape == (10, 4, 40, 10, 50) - - @pytest.mark.parametrize( - "batch_size, expected", - [ - ([(10, 42), None, (4, 42)], (10, 10, 4, 10, 50)), - ([(10, 42), Ellipsis, (4, 42)], (10, 10, 40, 10, 4)), - ([(10, 42), None, Ellipsis, (4, 42)], (10, 10, 40, 10, 4)), - ([10, None, Ellipsis, (4, 42)], (10, 10, 40, 10, 4)), - ], - ) - def test_special_batch_size(self, batch_size, expected): - mb = pm.Minibatch(self.data, batch_size) - assert mb.eval().shape == expected - - def test_cloning_available(self): - gop = pm.Minibatch(np.arange(100), 1) - res = gop**2 - shared = aesara.shared(np.array([10])) - res1 = aesara.clone_replace(res, {gop: shared}) - f = aesara.function([], res1) - assert f() == np.array([100]) - - def test_align(self): - m = pm.Minibatch(np.arange(1000), 1, random_seed=1) - n = pm.Minibatch(np.arange(1000), 1, random_seed=1) - f = aesara.function([], [m, n]) - n.eval() # not aligned - a, b = zip(*(f() for _ in range(1000))) - assert a != b - pm.align_minibatches() - a, b = zip(*(f() for _ in range(1000))) - assert a == b - n.eval() # not aligned - pm.align_minibatches([m]) - a, b = zip(*(f() for _ in range(1000))) - assert a != b - pm.align_minibatches([m, n]) - a, b = zip(*(f() for _ in range(1000))) - assert a == b From 690ae428fb34d1b40015486f4592880df9f3da5f Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 21:26:22 -0400 Subject: [PATCH 09/15] Merge test_modelcontext.py into test_model.py --- .github/workflows/tests.yml | 1 - pymc/tests/test_model.py | 66 +++++++++++++++++++++++++- pymc/tests/test_modelcontext.py | 83 --------------------------------- 3 files changed, 65 insertions(+), 85 deletions(-) delete mode 100644 pymc/tests/test_modelcontext.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 9a981c97ea4..d2250f9d42b 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -64,7 +64,6 @@ jobs: pymc/tests/test_types.py - | - pymc/tests/test_modelcontext.py pymc/tests/distributions/test_dist_math.py pymc/tests/test_pickling.py pymc/tests/distributions/test_transform.py diff --git a/pymc/tests/test_model.py b/pymc/tests/test_model.py index 6a74b88d8dc..ff2e22b6e6c 100644 --- a/pymc/tests/test_model.py +++ b/pymc/tests/test_model.py @@ -11,6 +11,7 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. +import threading import unittest import warnings @@ -37,7 +38,7 @@ from pymc.blocking import DictToArrayBijection, RaveledVars from pymc.distributions import Normal, transforms from pymc.exceptions import ImputationWarning, ShapeError, ShapeWarning -from pymc.model import Point, ValueGradFunction +from pymc.model import Point, ValueGradFunction, modelcontext from pymc.tests.helpers import SeededTest @@ -1037,3 +1038,66 @@ def test_model_parent_set_programmatically(): y = pm.Normal("y") assert "y" in model.named_vars + + +class TestModelContext: + def test_thread_safety(self): + """Regression test for issue #1552: Thread safety of model context manager + + This test creates two threads that attempt to construct two + unrelated models at the same time. + For repeatable testing, the two threads are syncronised such + that thread A enters the context manager first, then B, + then A attempts to declare a variable while B is still in the context manager. + """ + aInCtxt, bInCtxt, aDone = (threading.Event() for _ in range(3)) + modelA = pm.Model() + modelB = pm.Model() + + def make_model_a(): + with modelA: + aInCtxt.set() + bInCtxt.wait() + Normal("a", 0, 1) + aDone.set() + + def make_model_b(): + aInCtxt.wait() + with modelB: + bInCtxt.set() + aDone.wait() + Normal("b", 0, 1) + + threadA = threading.Thread(target=make_model_a) + threadB = threading.Thread(target=make_model_b) + threadA.start() + threadB.start() + threadA.join() + threadB.join() + # now let's see which model got which variable + # previous to #1555, the variables would be swapped: + # - B enters it's model context after A, but before a is declared -> a goes into B + # - A leaves it's model context before B attempts to declare b. A's context manager + # takes B from the stack, such that b ends up in model A + assert ( + list(modelA.named_vars), + list(modelB.named_vars), + ) == (["a"], ["b"]) + + +def test_mixed_contexts(): + modelA = pm.Model() + modelB = pm.Model() + with pytest.raises((ValueError, TypeError)): + modelcontext(None) + with modelA: + with modelB: + assert pm.Model.get_context() == modelB + assert modelcontext(None) == modelB + assert pm.Model.get_context() == modelA + assert modelcontext(None) == modelA + assert pm.Model.get_context(error_if_none=False) is None + with pytest.raises(TypeError): + pm.Model.get_context(error_if_none=True) + with pytest.raises((ValueError, TypeError)): + modelcontext(None) diff --git a/pymc/tests/test_modelcontext.py b/pymc/tests/test_modelcontext.py deleted file mode 100644 index 772db77b954..00000000000 --- a/pymc/tests/test_modelcontext.py +++ /dev/null @@ -1,83 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import threading - -from pytest import raises - -from pymc import Model, Normal -from pymc.model import modelcontext - - -class TestModelContext: - def test_thread_safety(self): - """Regression test for issue #1552: Thread safety of model context manager - - This test creates two threads that attempt to construct two - unrelated models at the same time. - For repeatable testing, the two threads are syncronised such - that thread A enters the context manager first, then B, - then A attempts to declare a variable while B is still in the context manager. - """ - aInCtxt, bInCtxt, aDone = (threading.Event() for _ in range(3)) - modelA = Model() - modelB = Model() - - def make_model_a(): - with modelA: - aInCtxt.set() - bInCtxt.wait() - Normal("a", 0, 1) - aDone.set() - - def make_model_b(): - aInCtxt.wait() - with modelB: - bInCtxt.set() - aDone.wait() - Normal("b", 0, 1) - - threadA = threading.Thread(target=make_model_a) - threadB = threading.Thread(target=make_model_b) - threadA.start() - threadB.start() - threadA.join() - threadB.join() - # now let's see which model got which variable - # previous to #1555, the variables would be swapped: - # - B enters it's model context after A, but before a is declared -> a goes into B - # - A leaves it's model context before B attempts to declare b. A's context manager - # takes B from the stack, such that b ends up in model A - assert ( - list(modelA.named_vars), - list(modelB.named_vars), - ) == (["a"], ["b"]) - - -def test_mixed_contexts(): - modelA = Model() - modelB = Model() - with raises((ValueError, TypeError)): - modelcontext(None) - with modelA: - with modelB: - assert Model.get_context() == modelB - assert modelcontext(None) == modelB - assert Model.get_context() == modelA - assert modelcontext(None) == modelA - assert Model.get_context(error_if_none=False) is None - with raises(TypeError): - Model.get_context(error_if_none=True) - with raises((ValueError, TypeError)): - modelcontext(None) From 4eaf00712150e34b47889fb7cd60a9cf3a3283d9 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 21:35:04 -0400 Subject: [PATCH 10/15] Merge test_pickling.py into test_model.py --- .github/workflows/tests.yml | 3 +-- pymc/tests/test_model.py | 45 +++++++++++++++++++++++++------------ pymc/tests/test_pickling.py | 37 ------------------------------ 3 files changed, 32 insertions(+), 53 deletions(-) delete mode 100644 pymc/tests/test_pickling.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index d2250f9d42b..e6476cb9c2a 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -65,7 +65,6 @@ jobs: - | pymc/tests/distributions/test_dist_math.py - pymc/tests/test_pickling.py pymc/tests/distributions/test_transform.py pymc/tests/test_parallel_sampling.py pymc/tests/test_printing.py @@ -149,7 +148,7 @@ jobs: python-version: ["3.8"] test-subset: - pymc/tests/test_variational_inference.py pymc/tests/test_initial_point.py - - pymc/tests/test_pickling.py pymc/tests/test_profile.py pymc/tests/test_step.py + - pymc/tests/test_model.py pymc/tests/test_profile.py pymc/tests/test_step.py - pymc/tests/gp/test_cov.py pymc/tests/gp/test_gp.py pymc/tests/gp/test_mean.py pymc/tests/gp/test_util.py pymc/tests/ode/test_ode.py pymc/tests/ode/test_utils.py pymc/tests/test_smc.py pymc/tests/test_parallel_sampling.py - pymc/tests/test_sampling.py pymc/tests/test_posteriors.py diff --git a/pymc/tests/test_model.py b/pymc/tests/test_model.py index ff2e22b6e6c..734bf81f6cd 100644 --- a/pymc/tests/test_model.py +++ b/pymc/tests/test_model.py @@ -11,6 +11,7 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. +import pickle import threading import unittest import warnings @@ -40,6 +41,7 @@ from pymc.exceptions import ImputationWarning, ShapeError, ShapeWarning from pymc.model import Point, ValueGradFunction, modelcontext from pymc.tests.helpers import SeededTest +from pymc.tests.models import simple_model class NewModel(pm.Model): @@ -444,24 +446,39 @@ def test_tempered_logp_dlogp(): npt.assert_allclose(func_temp_nograd(x), func_temp(x)[0]) -def test_model_pickle(tmpdir): - """Tests that PyMC models are pickleable""" - with pm.Model() as model: - x = pm.Normal("x") - pm.Normal("y", observed=1) +class TestPickling: + def test_model_pickle(self, tmpdir): + """Tests that PyMC models are pickleable""" + with pm.Model() as model: + x = pm.Normal("x") + pm.Normal("y", observed=1) - cloudpickle.loads(cloudpickle.dumps(model)) + cloudpickle.loads(cloudpickle.dumps(model)) + def test_model_pickle_deterministic(self, tmpdir): + """Tests that PyMC models are pickleable""" + with pm.Model() as model: + x = pm.Normal("x") + z = pm.Normal("z") + pm.Deterministic("w", x / z) + pm.Normal("y", observed=1) -def test_model_pickle_deterministic(tmpdir): - """Tests that PyMC models are pickleable""" - with pm.Model() as model: - x = pm.Normal("x") - z = pm.Normal("z") - pm.Deterministic("w", x / z) - pm.Normal("y", observed=1) + cloudpickle.loads(cloudpickle.dumps(model)) - cloudpickle.loads(cloudpickle.dumps(model)) + def setup_method(self): + _, self.model, _ = simple_model() + + def test_model_roundtrip(self): + m = self.model + for proto in range(pickle.HIGHEST_PROTOCOL + 1): + try: + s = cloudpickle.dumps(m, proto) + cloudpickle.loads(s) + except Exception: + raise AssertionError( + "Exception while trying roundtrip with pickle protocol %d:\n" % proto + + "".join(traceback.format_exc()) + ) def test_model_vars(): diff --git a/pymc/tests/test_pickling.py b/pymc/tests/test_pickling.py deleted file mode 100644 index d0a551a6e49..00000000000 --- a/pymc/tests/test_pickling.py +++ /dev/null @@ -1,37 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import pickle -import traceback - -import cloudpickle - -from pymc.tests.models import simple_model - - -class TestPickling: - def setup_method(self): - _, self.model, _ = simple_model() - - def test_model_roundtrip(self): - m = self.model - for proto in range(pickle.HIGHEST_PROTOCOL + 1): - try: - s = cloudpickle.dumps(m, proto) - cloudpickle.loads(s) - except Exception: - raise AssertionError( - "Exception while trying roundtrip with pickle protocol %d:\n" % proto - + "".join(traceback.format_exc()) - ) From 0eb481ad0e27d6bb11aa2980ab2f2e3fd2e2ca48 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sat, 10 Sep 2022 21:40:13 -0400 Subject: [PATCH 11/15] Move test_quadpotential.py to step_methods/hmc/ --- .github/workflows/tests.yml | 2 +- pymc/tests/step_methods/__init__.py | 0 pymc/tests/step_methods/hmc/__init__.py | 0 pymc/tests/{ => step_methods/hmc}/test_quadpotential.py | 0 4 files changed, 1 insertion(+), 1 deletion(-) create mode 100644 pymc/tests/step_methods/__init__.py create mode 100644 pymc/tests/step_methods/hmc/__init__.py rename pymc/tests/{ => step_methods/hmc}/test_quadpotential.py (100%) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index e6476cb9c2a..3f13b6aacb9 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -80,7 +80,7 @@ jobs: pymc/tests/ode/test_ode.py pymc/tests/ode/test_utils.py pymc/tests/test_profile.py - pymc/tests/test_quadpotential.py + pymc/tests/step_methods/hmc/test_quadpotential.py fail-fast: false runs-on: ${{ matrix.os }} diff --git a/pymc/tests/step_methods/__init__.py b/pymc/tests/step_methods/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pymc/tests/step_methods/hmc/__init__.py b/pymc/tests/step_methods/hmc/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pymc/tests/test_quadpotential.py b/pymc/tests/step_methods/hmc/test_quadpotential.py similarity index 100% rename from pymc/tests/test_quadpotential.py rename to pymc/tests/step_methods/hmc/test_quadpotential.py From ccde6105fc023233ef3bc8ac46bbce60bd48d822 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sun, 11 Sep 2022 00:22:18 -0400 Subject: [PATCH 12/15] Move test_posdef_sym into test_multivariate --- .github/workflows/tests.yml | 1 - pymc/tests/distributions/test_multivariate.py | 26 +++++++++++ pymc/tests/test_posdef_sym.py | 43 ------------------- 3 files changed, 26 insertions(+), 44 deletions(-) delete mode 100644 pymc/tests/test_posdef_sym.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 3f13b6aacb9..c6d58c9341c 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -41,7 +41,6 @@ jobs: pymc/tests/distributions/test_logprob.py pymc/tests/test_aesaraf.py pymc/tests/test_math.py - pymc/tests/test_posdef_sym.py pymc/tests/backends/test_ndarray.py pymc/tests/test_hmc.py pymc/tests/test_func_utils.py diff --git a/pymc/tests/distributions/test_multivariate.py b/pymc/tests/distributions/test_multivariate.py index 76052513df0..6800ced10a8 100644 --- a/pymc/tests/distributions/test_multivariate.py +++ b/pymc/tests/distributions/test_multivariate.py @@ -36,6 +36,7 @@ from pymc.distributions.multivariate import ( _LKJCholeskyCov, _OrderedMultinomial, + posdef, quaddist_matrix, ) from pymc.distributions.shape_utils import change_dist_size, to_tuple @@ -1874,3 +1875,28 @@ def test_car_rng_fn(sparse): ) f -= 1 assert p > delta + + +def test_posdef_symmetric1(): + data = np.array([[1.0, 0], [0, 1]], dtype=aesara.config.floatX) + assert posdef(data) == 1 + + +def test_posdef_symmetric2(): + data = np.array([[1.0, 2], [2, 1]], dtype=aesara.config.floatX) + assert posdef(data) == 0 + + +def test_posdef_symmetric3(): + """The test return 0 if the matrix has 0 eigenvalue. + + Is this correct? + """ + data = np.array([[1.0, 1], [1, 1]], dtype=aesara.config.floatX) + assert posdef(data) == 0 + + +def test_posdef_symmetric4(): + d = np.array([[1, 0.99, 1], [0.99, 1, 0.999], [1, 0.999, 1]], aesara.config.floatX) + + assert posdef(d) == 0 diff --git a/pymc/tests/test_posdef_sym.py b/pymc/tests/test_posdef_sym.py deleted file mode 100644 index fd429b43369..00000000000 --- a/pymc/tests/test_posdef_sym.py +++ /dev/null @@ -1,43 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import aesara -import numpy as np - -from pymc.distributions import multivariate as mv - - -def test_posdef_symmetric1(): - data = np.array([[1.0, 0], [0, 1]], dtype=aesara.config.floatX) - assert mv.posdef(data) == 1 - - -def test_posdef_symmetric2(): - data = np.array([[1.0, 2], [2, 1]], dtype=aesara.config.floatX) - assert mv.posdef(data) == 0 - - -def test_posdef_symmetric3(): - """The test return 0 if the matrix has 0 eigenvalue. - - Is this correct? - """ - data = np.array([[1.0, 1], [1, 1]], dtype=aesara.config.floatX) - assert mv.posdef(data) == 0 - - -def test_posdef_symmetric4(): - d = np.array([[1, 0.99, 1], [0.99, 1, 0.999], [1, 0.999, 1]], aesara.config.floatX) - - assert mv.posdef(d) == 0 From 79ba2c5db575f6107463bff245324ac3c21a92d4 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sun, 11 Sep 2022 00:24:50 -0400 Subject: [PATCH 13/15] Move test_profile into test_model --- .github/workflows/tests.yml | 3 +-- pymc/tests/test_model.py | 16 ++++++++++++++++ pymc/tests/test_profile.py | 31 ------------------------------- 3 files changed, 17 insertions(+), 33 deletions(-) delete mode 100644 pymc/tests/test_profile.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index c6d58c9341c..73b64cbf766 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -78,7 +78,6 @@ jobs: pymc/tests/test_model_graph.py pymc/tests/ode/test_ode.py pymc/tests/ode/test_utils.py - pymc/tests/test_profile.py pymc/tests/step_methods/hmc/test_quadpotential.py fail-fast: false @@ -147,7 +146,7 @@ jobs: python-version: ["3.8"] test-subset: - pymc/tests/test_variational_inference.py pymc/tests/test_initial_point.py - - pymc/tests/test_model.py pymc/tests/test_profile.py pymc/tests/test_step.py + - pymc/tests/test_model.py pymc/tests/test_step.py - pymc/tests/gp/test_cov.py pymc/tests/gp/test_gp.py pymc/tests/gp/test_mean.py pymc/tests/gp/test_util.py pymc/tests/ode/test_ode.py pymc/tests/ode/test_utils.py pymc/tests/test_smc.py pymc/tests/test_parallel_sampling.py - pymc/tests/test_sampling.py pymc/tests/test_posteriors.py diff --git a/pymc/tests/test_model.py b/pymc/tests/test_model.py index 734bf81f6cd..c7bea809fa1 100644 --- a/pymc/tests/test_model.py +++ b/pymc/tests/test_model.py @@ -1118,3 +1118,19 @@ def test_mixed_contexts(): pm.Model.get_context(error_if_none=True) with pytest.raises((ValueError, TypeError)): modelcontext(None) + + +class TestProfile: + def setup_method(self): + _, self.model, _ = simple_model() + + def test_profile_model(self): + assert self.model.profile(self.model.logp()).fct_call_time > 0 + + def test_profile_variable(self): + rv = self.model.basic_RVs[0] + assert self.model.profile(self.model.logp(vars=[rv], sum=False)).fct_call_time + + def test_profile_count(self): + count = 1005 + assert self.model.profile(self.model.logp(), n=count).fct_callcount == count diff --git a/pymc/tests/test_profile.py b/pymc/tests/test_profile.py deleted file mode 100644 index 55d54cbaee9..00000000000 --- a/pymc/tests/test_profile.py +++ /dev/null @@ -1,31 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -from pymc.tests.models import simple_model - - -class TestProfile: - def setup_method(self): - _, self.model, _ = simple_model() - - def test_profile_model(self): - assert self.model.profile(self.model.logp()).fct_call_time > 0 - - def test_profile_variable(self): - rv = self.model.basic_RVs[0] - assert self.model.profile(self.model.logp(vars=[rv], sum=False)).fct_call_time - - def test_profile_count(self): - count = 1005 - assert self.model.profile(self.model.logp(), n=count).fct_callcount == count From 905bc422dfdc7d81737813d968da47964aca8de6 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sun, 11 Sep 2022 00:28:59 -0400 Subject: [PATCH 14/15] Move test_updates.py into variational/ --- .github/workflows/tests.yml | 2 +- pymc/tests/variational/__init__.py | 0 pymc/tests/{ => variational}/test_updates.py | 0 3 files changed, 1 insertion(+), 1 deletion(-) create mode 100644 pymc/tests/variational/__init__.py rename pymc/tests/{ => variational}/test_updates.py (100%) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 73b64cbf766..31f22eac2a9 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -227,7 +227,7 @@ jobs: - | pymc/tests/backends/test_arviz.py - pymc/tests/test_updates.py + pymc/tests/variational/test_updates.py fail-fast: false runs-on: ${{ matrix.os }} env: diff --git a/pymc/tests/variational/__init__.py b/pymc/tests/variational/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pymc/tests/test_updates.py b/pymc/tests/variational/test_updates.py similarity index 100% rename from pymc/tests/test_updates.py rename to pymc/tests/variational/test_updates.py From 1558472e8c6f200571c31e08dd196e44e0eba280 Mon Sep 17 00:00:00 2001 From: Virgile Andreani Date: Sun, 11 Sep 2022 15:09:52 -0400 Subject: [PATCH 15/15] Parametrize the posdef_symmetric test --- pymc/tests/distributions/test_multivariate.py | 33 ++++++++----------- 1 file changed, 13 insertions(+), 20 deletions(-) diff --git a/pymc/tests/distributions/test_multivariate.py b/pymc/tests/distributions/test_multivariate.py index 6800ced10a8..2dde22c915f 100644 --- a/pymc/tests/distributions/test_multivariate.py +++ b/pymc/tests/distributions/test_multivariate.py @@ -1877,26 +1877,19 @@ def test_car_rng_fn(sparse): assert p > delta -def test_posdef_symmetric1(): - data = np.array([[1.0, 0], [0, 1]], dtype=aesara.config.floatX) - assert posdef(data) == 1 - - -def test_posdef_symmetric2(): - data = np.array([[1.0, 2], [2, 1]], dtype=aesara.config.floatX) - assert posdef(data) == 0 - - -def test_posdef_symmetric3(): - """The test return 0 if the matrix has 0 eigenvalue. +@pytest.mark.parametrize( + "matrix, result", + [ + ([[1.0, 0], [0, 1]], 1), + ([[1.0, 2], [2, 1]], 0), + ([[1.0, 1], [1, 1]], 0), + ([[1, 0.99, 1], [0.99, 1, 0.999], [1, 0.999, 1]], 0), + ], +) +def test_posdef_symmetric(matrix, result): + """The test returns 0 if the matrix has 0 eigenvalue. Is this correct? """ - data = np.array([[1.0, 1], [1, 1]], dtype=aesara.config.floatX) - assert posdef(data) == 0 - - -def test_posdef_symmetric4(): - d = np.array([[1, 0.99, 1], [0.99, 1, 0.999], [1, 0.999, 1]], aesara.config.floatX) - - assert posdef(d) == 0 + data = np.array(matrix, dtype=aesara.config.floatX) + assert posdef(data) == result