From a9eed9fbdc5c9925c398dc35e66bd6e9823ea4d7 Mon Sep 17 00:00:00 2001 From: ben18785 Date: Fri, 12 Mar 2021 21:56:34 +0000 Subject: [PATCH] Added tests to Multinest addresses #282 --- examples/sampling/nested-multinest.ipynb | 33 +++++ pints/_nested/_multinest.py | 1 + pints/tests/test_nested_multinest.py | 156 +++++++++++++++-------- 3 files changed, 137 insertions(+), 53 deletions(-) diff --git a/examples/sampling/nested-multinest.ipynb b/examples/sampling/nested-multinest.ipynb index 3e66de8f4..ad3483cf0 100644 --- a/examples/sampling/nested-multinest.ipynb +++ b/examples/sampling/nested-multinest.ipynb @@ -1306,6 +1306,39 @@ "pints.plot.pairwise(samples, kde=True, ref_parameters=real_parameters)\n", "plt.show()" ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "ename": "LinAlgError", + "evalue": "Singular matrix", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mpoints\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mcov\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcov\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtranspose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mcov_inv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcov\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36minv\u001b[0;34m(*args, **kwargs)\u001b[0m\n", + "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m 543\u001b[0m \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'D->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 544\u001b[0m \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 545\u001b[0;31m \u001b[0mainv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 546\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 547\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m 86\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 90\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix" + ] + } + ], + "source": [ + "points = np.array([[0, 0], [0.5, 0.5], [1, 1]])\n", + "cov = np.cov(np.transpose(points))\n", + "cov_inv = np.linalg.inv(cov)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/pints/_nested/_multinest.py b/pints/_nested/_multinest.py index 42ee07723..d1e864ca7 100644 --- a/pints/_nested/_multinest.py +++ b/pints/_nested/_multinest.py @@ -178,6 +178,7 @@ def __init__(self, log_prior): self._convert_to_unit_cube = log_prior.convert_to_unit_cube self._convert_from_unit_cube = log_prior.convert_from_unit_cube + self._ellipsoid_tree = None def ask(self, n_points): """ diff --git a/pints/tests/test_nested_multinest.py b/pints/tests/test_nested_multinest.py index cee175523..bc5f5ef06 100644 --- a/pints/tests/test_nested_multinest.py +++ b/pints/tests/test_nested_multinest.py @@ -21,58 +21,63 @@ unittest.TestCase.assertRaisesRegex = unittest.TestCase.assertRaisesRegexp -# class TestMultiNestSampler(unittest.TestCase): -# """ -# Unit (not functional!) tests for :class:`MultinestSampler`. -# """ -# -# @classmethod -# def setUpClass(cls): -# """ Prepare for the test. """ -# # Create toy model -# model = pints.toy.LogisticModel() -# cls.real_parameters = [0.015, 500] -# times = np.linspace(0, 1000, 1000) -# values = model.simulate(cls.real_parameters, times) -# -# # Add noise -# np.random.seed(1) -# cls.noise = 10 -# values += np.random.normal(0, cls.noise, values.shape) -# cls.real_parameters.append(cls.noise) -# -# # Create an object with links to the model and time series -# problem = pints.SingleOutputProblem(model, times, values) -# -# # Create a uniform prior over both the parameters and the new noise -# # variable -# cls.log_prior = pints.UniformLogPrior( -# [0.01, 400], -# [0.02, 600] -# ) -# -# # Create a log-likelihood -# cls.log_likelihood = pints.GaussianKnownSigmaLogLikelihood( -# problem, cls.noise) -# -# def test_getters_and_setters(self): -# # tests various get() and set() methods. -# controller = pints.NestedController(self.log_likelihood, -# self.log_prior, -# method=pints.MultinestSampler) -# self.assertEqual(controller.sampler().f_s_threshold(), 1.1) -# controller.sampler().set_f_s_threshold(4) -# self.assertEqual(controller.sampler().f_s_threshold(), 4) -# self.assertRaises(ValueError, controller.sampler().set_f_s_threshold, -# 0.5) -# -# def test_runs(self): -# # tests that sampler runs -# sampler = pints.NestedController(self.log_likelihood, self.log_prior, -# method=pints.MultinestSampler) -# sampler.set_iterations(100) -# sampler.set_log_to_screen(False) -# sampler.run() +class TestMultiNestSampler(unittest.TestCase): + """ + Unit (not functional!) tests for :class:`MultinestSampler`. + """ + + @classmethod + def setUpClass(cls): + """ Prepare for the test. """ + # Create toy model + model = pints.toy.LogisticModel() + cls.real_parameters = [0.015, 500] + times = np.linspace(0, 1000, 1000) + values = model.simulate(cls.real_parameters, times) + + # Add noise + np.random.seed(1) + cls.noise = 10 + values += np.random.normal(0, cls.noise, values.shape) + cls.real_parameters.append(cls.noise) + + # Create an object with links to the model and time series + problem = pints.SingleOutputProblem(model, times, values) + + # Create a uniform prior over both the parameters and the new noise + # variable + cls.log_prior = pints.UniformLogPrior( + [0.01, 400], + [0.02, 600] + ) + + # Create a log-likelihood + cls.log_likelihood = pints.GaussianKnownSigmaLogLikelihood( + problem, cls.noise) + + def test_getters_and_setters(self): + # tests various get() and set() methods. + controller = pints.NestedController(self.log_likelihood, + self.log_prior, + method=pints.MultinestSampler) + self.assertEqual(controller.sampler().f_s_threshold(), 1.1) + controller.sampler().set_f_s_threshold(4) + self.assertEqual(controller.sampler().f_s_threshold(), 4) + self.assertRaises(ValueError, controller.sampler().set_f_s_threshold, + 0.5) + + def test_runs(self): + # tests that sampler runs + controller = pints.NestedController(self.log_likelihood, + self.log_prior, + method=pints.MultinestSampler) + controller.set_iterations(450) + controller.set_log_to_screen(False) + controller.run() + + # test getting ellipsoid tree post-run + et = controller.sampler().ellipsoid_tree() + self.assertTrue(et.n_leaf_ellipsoids() >= 1) class TestEllipsoidTree(unittest.TestCase): @@ -121,7 +126,7 @@ def test_getters(self): # leaf nodes self.assertTrue(tree.n_leaf_ellipsoids() >= 1) - leaves = tree.leaf_ellipsoids() + leaves = np.copy(tree.leaf_ellipsoids()) self.assertEqual(tree.n_leaf_ellipsoids(), len(leaves)) [self.assertTrue(isinstance(x, Ellipsoid)) for x in leaves] @@ -129,6 +134,20 @@ def test_getters(self): ellipsoid = tree.ellipsoid() self.assertTrue(isinstance(ellipsoid, Ellipsoid)) + # Check a given draw is within bounding ellipsoids + self.assertTrue(tree.count_within_leaf_ellipsoids( + self.draws[0]) >= 1) + self.assertTrue(tree.count_within_leaf_ellipsoids( + [-100, -100]) == 0) + + # check f_s + self.assertTrue(tree.f_s() >= 1) + + # updating triggers changes? + f_s_old = np.copy(tree.f_s()) + tree.update_leaf_ellipsoids(1) + self.assertFalse(tree.f_s() == f_s_old) + def test_calculations(self): # tests the calculations done @@ -152,7 +171,38 @@ def test_calculations(self): ellipsoid = Ellipsoid.minimum_volume_ellipsoid(draws) self.assertRaises(ValueError, ellipsoidTree.vsk, ellipsoid) + # hk + ellipsoid = ellipsoidTree.ellipsoid() + d = Ellipsoid.mahalanobis_distance(self.draws[0], + ellipsoid.weight_matrix(), + ellipsoid.centroid()) + V_S_k = 1 + self.assertEqual(ellipsoid.volume() * d / V_S_k, + ellipsoidTree.h_k(self.draws[0], ellipsoid, V_S_k)) + + def test_sampling(self): + # tests sampling from leaf ellipsoids + # check all samples in leaf ellipsoids + tree = EllipsoidTree(self.draws, 100) + samples = tree.sample_leaf_ellipsoids(100) + for sample in samples: + self.assertTrue(tree.count_within_leaf_ellipsoids( + sample) >= 1) + + def test_splitting(self): + # tests splitting and enlarging ellipsoids + + # singular matrix error? + tree = EllipsoidTree(self.draws, 100) + points = np.array([[0, 0], [0.5, 0.5], [1, 1]]) + assignments = np.array([0, 0, 0]) + _, _, bool = tree.split_ellipsoids(points, assignments, 0) + self.assertFalse(bool) + + # enlarging + ellipsoid = tree.ellipsoid() + tree.compare_enlarge(ellipsoid, 20) if __name__ == '__main__':