Merge pull request #38 from eth-cscs/release-0.5.1

Release 0.5.1

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

VERSION

@@ -1 +1 @@
-0.5.0
+0.5.1

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:

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
 

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)

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)  
 

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.
 

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,33 +230,37 @@ 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]
 
         statistics_calculator = Identity(degree=2, cross=False)
         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__':

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)