diff --git a/README.md b/README.md index a7e0db1d..e673006d 100644 --- a/README.md +++ b/README.md @@ -26,9 +26,9 @@ scientists by providing # Documentation For more information, check out the -* [Documentation](http://abcpy.readthedocs.io/en/v0.5.0) -* [Examples](https://github.com/eth-cscs/abcpy/tree/v0.5.0/examples) directory and -* [Reference](http://abcpy.readthedocs.io/en/v0.5.0/abcpy.html) +* [Documentation](http://abcpy.readthedocs.io/en/v0.5.1) +* [Examples](https://github.com/eth-cscs/abcpy/tree/v0.5.1/examples) directory and +* [Reference](http://abcpy.readthedocs.io/en/v0.5.1/abcpy.html) Further, we provide a [collection of models](https://github.com/eth-cscs/abcpy-models) for which ABCpy diff --git a/VERSION b/VERSION index 8f0916f7..4b9fcbec 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.5.0 +0.5.1 diff --git a/abcpy/continuousmodels.py b/abcpy/continuousmodels.py index 45926ea0..a06be373 100644 --- a/abcpy/continuousmodels.py +++ b/abcpy/continuousmodels.py @@ -114,6 +114,7 @@ def pdf(self, input_values, x): lower_bound = input_values[:self.get_output_dimension()] upper_bound = input_values[self.get_output_dimension():] + if (np.product(np.greater_equal(x, np.array(lower_bound)) * np.less_equal(x, np.array(upper_bound)))): pdf_value = 1. / np.product(np.array(upper_bound) - np.array(lower_bound)) else: diff --git a/abcpy/graphtools.py b/abcpy/graphtools.py index aecd56c6..d0f706ef 100644 --- a/abcpy/graphtools.py +++ b/abcpy/graphtools.py @@ -90,6 +90,31 @@ def pdf_of_prior(self, models, parameters, mapping=None, is_root=True): Calculates the joint probability density function of the prior of the specified models at the given parameter values. Commonly used to check whether new parameters are valid given the prior, as well as to calculate acceptance probabilities. + Parameters + ---------- + models: list of abcpy.ProbabilisticModel objects + Defines the models for which the pdf of their prior should be evaluated + parameters: python list + The parameters at which the pdf should be evaluated + mapping: list of tupels + Defines the mapping of probabilistic models and index in a parameter list. + is_root: boolean + A flag specifying whether the provided models are the root models. This is to ensure that the pdf is calculated correctly. + + Returns + ------- + list + The resulting pdf,as well as the next index to be considered in the parameters list. + """ + self.set_parameters(parameters) + result = self._recursion_pdf_of_prior(models, parameters, mapping, is_root) + return result + + def _recursion_pdf_of_prior(self, models, parameters, mapping=None, is_root=True): + """ + Calculates the joint probability density function of the prior of the specified models at the given parameter values. + Commonly used to check whether new parameters are valid given the prior, as well as to calculate acceptance probabilities. + Parameters ---------- models: list of abcpy.ProbabilisticModel objects @@ -138,23 +163,23 @@ def pdf_of_prior(self, models, parameters, mapping=None, is_root=True): for parent_index, parent in enumerate(model.get_input_models()): # Only calculate the pdf if the parent has never been visited for this model if(not(visited_parents[parent_index])): - pdf = self.pdf_of_prior([parent], parameters, mapping=mapping, is_root=False) + pdf = self._recursion_pdf_of_prior([parent], parameters, mapping=mapping, is_root=False) input_models = model.get_input_models() for j in range(len(input_models)): if input_models[j][0] == parent: visited_parents[j]=True - result[i]*=pdf + result[i] *= pdf if(not(is_root)): if(model.calculated_pdf is None): result[i] *= model.pdf(model.get_input_values(),relevant_parameters) else: - result[i]*=model.calculated_pdf + result[i] *= 1 # Multiply the pdfs of all roots together to give an overall pdf. temporary_result = result result = 1. for individual_result in temporary_result: - result*=individual_result + result *= individual_result return result diff --git a/abcpy/inferences.py b/abcpy/inferences.py index 9e9cc993..cf7e57f6 100644 --- a/abcpy/inferences.py +++ b/abcpy/inferences.py @@ -234,7 +234,7 @@ def sample(self, observations, n_samples, n_samples_per_param, epsilon, full_out journal.add_parameters(accepted_parameters) journal.add_weights(np.ones((n_samples, 1))) - + self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters) names_and_parameters = self._get_names_and_parameters() journal.add_user_parameters(names_and_parameters) @@ -482,7 +482,8 @@ def sample(self, observations, steps, epsilon_init, n_samples = 10000, n_samples if (full_output == 1 and aStep <= steps - 1) or (full_output == 0 and aStep == steps - 1): journal.add_parameters(accepted_parameters) journal.add_weights(accepted_weights) - + self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, + accepted_weights=accepted_weights) names_and_parameters = self._get_names_and_parameters() journal.add_user_parameters(names_and_parameters) @@ -830,7 +831,8 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1 journal.add_parameters(accepted_parameters) journal.add_weights(accepted_weights) journal.add_opt_values(approx_likelihood_new_parameters) - + self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, + accepted_weights=accepted_weights) names_and_parameters = self._get_names_and_parameters() journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) @@ -1050,6 +1052,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_ samples_until = 0 for aStep in range(0, steps): + print(aStep) if(aStep==0 and journal_file is not None): accepted_parameters=journal.parameters[-1] accepted_weights=journal.weights[-1] @@ -1103,7 +1106,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_ new_parameters = np.array(new_parameters) new_distances = np.array(new_distances) new_all_distances = np.concatenate(new_all_distances) - index = index_arr + index = np.array(index) acceptance = np.array(acceptance) # Reading all_distances at Initial step @@ -1152,7 +1155,8 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_ if full_output == 1: journal.add_parameters(accepted_parameters) journal.add_weights(accepted_weights) - + self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, + accepted_weights=accepted_weights) names_and_parameters = self._get_names_and_parameters() journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) @@ -1191,7 +1195,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_ if full_output == 0: journal.add_parameters(accepted_parameters) journal.add_weights(accepted_weights) - + self.accepted_parameters_manager.update_broadcast(self.backend,accepted_parameters=accepted_parameters,accepted_weights=accepted_weights) names_and_parameters = self._get_names_and_parameters() journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) @@ -1311,7 +1315,6 @@ def _accept_parameter(self, data): distance = self.distance.distance(self.accepted_parameters_manager.observations_bds.value(), y_sim) all_distances.append(distance) acceptance = rng.binomial(1, np.exp(-distance / self.epsilon), 1) - acceptance = 1 else: ## Select one arbitrary particle: @@ -1545,6 +1548,8 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1 journal.add_parameters(accepted_parameters) journal.add_weights(accepted_weights) journal.add_opt_values(accepted_cov_mats) + self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, + accepted_weights=accepted_weights) names_and_parameters = self._get_names_and_parameters() journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) @@ -1562,7 +1567,8 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1 journal.add_parameters(accepted_parameters) journal.add_weights(accepted_weights) journal.add_opt_values(accepted_cov_mats) - + self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, + accepted_weights=accepted_weights) names_and_parameters = self._get_names_and_parameters() journal.add_user_parameters(names_and_parameters) journal.number_of_simulations.append(self.simulation_counter) diff --git a/abcpy/statistics.py b/abcpy/statistics.py index 00f6af75..a49ec07b 100644 --- a/abcpy/statistics.py +++ b/abcpy/statistics.py @@ -71,10 +71,10 @@ def _polynomial_expansion(self, summary_statistics): if not isinstance(summary_statistics, (np.ndarray)): raise TypeError('Summary statisticss is not of allowed types') # Include the polynomial expansion - result = summary_statistics + result = summary_statistics for ind in range(2,self.degree+1): result = np.column_stack((result,np.power(summary_statistics,ind))) - + # Include the cross-product term if self.cross == True and summary_statistics.shape[1]>1: # Convert to a matrix @@ -93,7 +93,6 @@ class Identity(Statistics): def __init__(self, degree = 2, cross = True): self.degree = degree self.cross = cross - def statistics(self, data): if isinstance(data, list): @@ -105,7 +104,6 @@ def statistics(self, data): data = np.concatenate(data).reshape(len(data),-1) else: raise TypeError('Input data should be of type list, but found type {}'.format(type(data))) - # Expand the data with polynomial expansion result = self._polynomial_expansion(data) diff --git a/doc/source/installation.rst b/doc/source/installation.rst index 83520cb5..7b65a1b5 100644 --- a/doc/source/installation.rst +++ b/doc/source/installation.rst @@ -34,7 +34,7 @@ To create a package and install it do :: make package - pip3 install build/dist/abcpy-0.4.0-py3-none-any.whl + pip3 install build/dist/abcpy-0.5.1-py3-none-any.whl Note that ABCpy requires Python3. diff --git a/tests/graphtools_tests.py b/tests/graphtools_tests.py index 23d92524..631f04e8 100644 --- a/tests/graphtools_tests.py +++ b/tests/graphtools_tests.py @@ -219,6 +219,7 @@ def test(self): mapping, index = sampler._get_mapping() self.assertTrue(mapping==[(B1, 0),(N2, 1),(N1,2)]) +from abcpy.continuousmodels import Uniform class PdfOfPriorTests(unittest.TestCase): """Tests the implemetation of pdf_of_prior""" @@ -229,11 +230,11 @@ def __init__(self, parameters): def pdf(self, input_values, x): return x - self.N1 = Mockobject([1, 0.1]) - self.N2 = Mockobject([self.N1, 0.1]) - self.N3 = Mockobject([0.1, 0.01]) - self.graph1 = Mockobject([self.N2, self.N3]) - self.graph2 = Mockobject([2, self.N3]) + self.N1 = Uniform([[1.3], [1.55]], name='n1') + self.N2 = Uniform([[self.N1], [1.60]], name='n2') + self.N3 = Uniform([[2], [4]], name='n3') + self.graph1 = Uniform([[self.N1], [self.N2]], name='graph1') + self.graph2 = Uniform([[.5], [self.N3]]) self.graph = [self.graph1, self.graph2] @@ -241,21 +242,25 @@ def pdf(self, input_values, x): distance_calculator = LogReg(statistics_calculator) backend = Backend() - self.sampler = RejectionABC(self.graph, [distance_calculator, distance_calculator], backend) + self.sampler1 = RejectionABC([self.graph1], [distance_calculator], backend) + self.sampler2 = RejectionABC([self.graph2], [distance_calculator], backend) + self.sampler3 = RejectionABC(self.graph, [distance_calculator, distance_calculator], backend) - rng = np.random.RandomState(1) - - self.sampler.sample_from_prior(rng=rng) - - self.pdf = self.sampler.pdf_of_prior(self.sampler.model, [1, 2, 4]) + self.pdf1 = self.sampler1.pdf_of_prior(self.sampler1.model, [1.32088846, 1.42945274]) + self.pdf2 = self.sampler2.pdf_of_prior(self.sampler2.model, [3]) + self.pdf3 = self.sampler3.pdf_of_prior(self.sampler3.model, [1.32088846, 1.42945274, 3]) def test_return_value(self): """Tests whether the return value is float.""" - self.assertTrue(isinstance(self.pdf, float)) + self.assertTrue(isinstance(self.pdf1, float)) + self.assertTrue(isinstance(self.pdf2, float)) + self.assertTrue(isinstance(self.pdf3, float)) def test_result(self): """Test whether pdf calculation works as intended""" - self.assertTrue(self.pdf==32) + self.assertTrue(self.pdf1 == 14.331188169432183) + self.assertTrue(self.pdf2 == 0.5) + self.assertTrue(self.pdf3 == 7.1655940847160915) if __name__ == '__main__': diff --git a/tests/inferences_tests.py b/tests/inferences_tests.py index ab401209..3516e1f8 100644 --- a/tests/inferences_tests.py +++ b/tests/inferences_tests.py @@ -91,8 +91,8 @@ def test_sample(self): self.assertEqual(mu_sample_shape, (10,1)) self.assertEqual(sigma_sample_shape, (10,1)) self.assertEqual(weights_sample_shape, (10,1)) - self.assertLess(abs(mu_post_mean - (-3.64971)), 1e-3) - self.assertLess(abs(sigma_post_mean - 4.1925), 1e-3) + self.assertLess(abs(mu_post_mean - (-3.55548377)), 1e-3) + self.assertLess(abs(sigma_post_mean - 5.87147492), 1e-3) self.assertFalse(journal.number_of_simulations == 0) @@ -111,8 +111,8 @@ def test_sample(self): self.assertEqual(mu_sample_shape, (10,1)) self.assertEqual(sigma_sample_shape, (10,1)) self.assertEqual(weights_sample_shape, (10,1)) - self.assertLess(abs(mu_post_mean - (-3.09297) ), 1e-3) - self.assertLess(abs(sigma_post_mean - 5.78645), 1e-3) + self.assertLess(abs(mu_post_mean - (-3.88512394) ), 1e-3) + self.assertLess(abs(sigma_post_mean - 2.90352432), 1e-3) self.assertFalse(journal.number_of_simulations == 0) @@ -254,8 +254,8 @@ def test_sample(self): self.assertEqual(mu_sample_shape, (10,1)) self.assertEqual(sigma_sample_shape, (10,1)) self.assertEqual(weights_sample_shape, (10,1)) - self.assertLess(mu_post_mean - 2.120856674879079, 10e-2) - self.assertLess(sigma_post_mean - 6.711723792285109, 10e-2) + self.assertLess(mu_post_mean - 0.55859197, 10e-2) + self.assertLess(sigma_post_mean - 7.03987723, 10e-2) self.assertFalse(journal.number_of_simulations == 0) @@ -314,7 +314,7 @@ def test_sample(self): self.assertEqual(sigma_sample_shape, (10,1)) self.assertEqual(weights_sample_shape, (10,1)) self.assertLess(mu_post_mean - (-0.81410299), 10e-2) - self.assertLess(sigma_post_mean - 5.51827908, 10e-2) + self.assertLess(sigma_post_mean - 9.25442675, 10e-2) self.assertFalse(journal.number_of_simulations == 0)