Merge pull request #39 from eth-cscs/release-0.5.2

Release 0.5.2
eth-cscs · Jul 31, 2018 · ad12808 · ad12808
2 parents 1beaa48 + 3318d6f
commit ad12808
Show file tree

Hide file tree

Showing 15 changed files with 187 additions and 162 deletions.
diff --git 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.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)
+* [Documentation](http://abcpy.readthedocs.io/en/v0.5.2) 
+* [Examples](https://github.com/eth-cscs/abcpy/tree/v0.5.2/examples) directory and
+* [Reference](http://abcpy.readthedocs.io/en/v0.5.2/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
@@ -1 +1 @@
-0.5.1
+0.5.2
diff --git a/abcpy/inferences.py b/abcpy/inferences.py
@@ -6,7 +6,7 @@
 from abcpy.perturbationkernel import DefaultKernel
 from abcpy.jointdistances import LinearCombination
 from abcpy.jointapprox_lhd import ProductCombination
-
+import copy
 
 import numpy as np
 from abcpy.output import Journal
@@ -1008,6 +1008,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_
         abcpy.output.Journal
             A journal containing simulation results, metadata and optionally intermediate results.
         """
+        global broken_preemptively
         self.sample_from_prior(rng=self.rng)
         self.accepted_parameters_manager.broadcast(self.backend, observations)
         self.epsilon = epsilon
@@ -1051,19 +1052,24 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_
         accept = 0
         samples_until = 0
 
+        ## Counter whether broken preemptively
+        broken_preemptively = False
+
         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]
 
+                #Broadcast Accepted parameters and Accedpted weights
                 self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters, accepted_weights=accepted_weights)
 
                 kernel_parameters = []
                 for kernel in self.kernel.kernels:
                     kernel_parameters.append(
                         self.accepted_parameters_manager.get_accepted_parameters_bds_values(kernel.models))
 
+                #Broadcast Accepted Kernel parameters
                 self.accepted_parameters_manager.update_kernel_values(self.backend, kernel_parameters=kernel_parameters)
 
                 new_cov_mats = self.kernel.calculate_cov(self.accepted_parameters_manager)
@@ -1074,7 +1080,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_
                 else:
                     accepted_cov_mats = [beta*new_cov_mat + 0.0001*(new_cov_mat)*np.eye(accepted_parameters.shape[1]) for new_cov_mat in new_cov_mats]
 
-
+                # Broadcast Accepted Covariance Matrix
                 self.accepted_parameters_manager.update_broadcast(self.backend, accepted_cov_mats=accepted_cov_mats)
 
             # main SABC algorithm
@@ -1100,6 +1106,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_
                                                                                                        zip(
                                                                                                            *params_and_dists)]
 
+            # Keeping counter of number of simulations
             for count in counter:
                 self.simulation_counter+=count
 
@@ -1130,37 +1137,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_
                 U = np.mean(smooth_distances)
             epsilon = self._schedule(U, v)
 
-            self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters)
-
-            kernel_parameters = []
-            for kernel in self.kernel.kernels:
-                kernel_parameters.append(
-                    self.accepted_parameters_manager.get_accepted_parameters_bds_values(kernel.models))
-
-            self.accepted_parameters_manager.update_kernel_values(self.backend, kernel_parameters=kernel_parameters)
-
-            new_cov_mats = self.kernel.calculate_cov(self.accepted_parameters_manager)
-
-            if accepted_parameters.shape[1] > 1:
-                accepted_cov_mats = [beta*new_cov_mat+0.0001*np.trace(new_cov_mat)*np.eye(len(new_cov_mat)) for new_cov_mat in new_cov_mats]
-            else:
-                accepted_cov_mats = [beta * new_cov_mat + 0.0001 * (new_cov_mat) * np.eye(accepted_parameters.shape[1])
-                                     for new_cov_mat in new_cov_mats]
-
             # 4: Show progress and if acceptance rate smaller than a value break the iteration
-
-            self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters,accepted_cov_mats=accepted_cov_mats)
-
-            # print("INFO: Saving intermediate configuration to output journal.")
-            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)
-
             if aStep > 0:
                 accept = accept + np.sum(acceptance)
                 samples_until = samples_until + sample_array[aStep]
@@ -1169,6 +1146,7 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_
                 'updates: ', np.sum(sample_array[1:aStep + 1]) / np.sum(sample_array[1:]) * 100, ' epsilon: ', epsilon, \
                 'u.mean: ', U, 'acceptance rate: ', acceptance_rate)
                 if acceptance_rate < ar_cutoff:
+                    broken_preemptively = True
                     break
 
             # 5: Resampling if number of accepted particles greater than resample
@@ -1190,11 +1168,71 @@ def sample(self, observations, steps, epsilon, n_samples = 10000, n_samples_per_
                 accept = 0
                 samples_until = 0
 
+                ## Compute and broadcast accepted parameters, accepted kernel parameters and accepted Covariance matrix
+                # Broadcast Accepted parameters and add to journal
+                self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters)
+                # Compute Accepetd Kernel parameters and broadcast them
+                kernel_parameters = []
+                for kernel in self.kernel.kernels:
+                    kernel_parameters.append(
+                        self.accepted_parameters_manager.get_accepted_parameters_bds_values(kernel.models))
+                self.accepted_parameters_manager.update_kernel_values(self.backend, kernel_parameters=kernel_parameters)
+                # Compute Kernel Covariance Matrix and broadcast it
+                new_cov_mats = self.kernel.calculate_cov(self.accepted_parameters_manager)
+                if accepted_parameters.shape[1] > 1:
+                    accepted_cov_mats = [beta * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(len(new_cov_mat))
+                                         for new_cov_mat in new_cov_mats]
+                else:
+                    accepted_cov_mats = [
+                        beta * new_cov_mat + 0.0001 * (new_cov_mat) * np.eye(accepted_parameters.shape[1])
+                        for new_cov_mat in new_cov_mats]
+                self.accepted_parameters_manager.update_broadcast(self.backend, accepted_cov_mats=accepted_cov_mats)
+
+                if (full_output == 1 and aStep<= steps-1):
+                    ## Saving intermediate configuration to output journal.
+                    print('Saving after resampling')
+                    journal.add_parameters(copy.deepcopy(accepted_parameters))
+                    journal.add_weights(copy.deepcopy(accepted_weights))
+                    journal.add_distances(copy.deepcopy(distances))
+                    names_and_parameters = self._get_names_and_parameters()
+                    journal.add_user_parameters(names_and_parameters)
+                    journal.number_of_simulations.append(self.simulation_counter)
+            else:
+                ## Compute and broadcast accepted parameters, accepted kernel parameters and accepted Covariance matrix
+                # Broadcast Accepted parameters
+                self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters)
+                # Compute Accepetd Kernel parameters and broadcast them
+                kernel_parameters = []
+                for kernel in self.kernel.kernels:
+                    kernel_parameters.append(
+                        self.accepted_parameters_manager.get_accepted_parameters_bds_values(kernel.models))
+                self.accepted_parameters_manager.update_kernel_values(self.backend, kernel_parameters=kernel_parameters)
+                # Compute Kernel Covariance Matrix and broadcast it
+                new_cov_mats = self.kernel.calculate_cov(self.accepted_parameters_manager)
+                if accepted_parameters.shape[1] > 1:
+                    accepted_cov_mats = [beta * new_cov_mat + 0.0001 * np.trace(new_cov_mat) * np.eye(len(new_cov_mat))
+                                         for new_cov_mat in new_cov_mats]
+                else:
+                    accepted_cov_mats = [
+                        beta * new_cov_mat + 0.0001 * (new_cov_mat) * np.eye(accepted_parameters.shape[1])
+                        for new_cov_mat in new_cov_mats]
+                self.accepted_parameters_manager.update_broadcast(self.backend, accepted_cov_mats=accepted_cov_mats)
+
+                if (full_output == 1 and aStep <= steps-1):
+                    ## Saving intermediate configuration to output journal.
+                    journal.add_parameters(copy.deepcopy(accepted_parameters))
+                    journal.add_weights(copy.deepcopy(accepted_weights))
+                    journal.add_distances(copy.deepcopy(distances))
+                    names_and_parameters = self._get_names_and_parameters()
+                    journal.add_user_parameters(names_and_parameters)
+                    journal.number_of_simulations.append(self.simulation_counter)
+
         # Add epsilon_arr, number of final steps and final output to the journal
         # print("INFO: Saving final configuration to output journal.")
-        if full_output == 0:
-            journal.add_parameters(accepted_parameters)
-            journal.add_weights(accepted_weights)
+        if (full_output == 0) or (full_output ==1 and broken_preemptively and aStep<= steps-1):
+            journal.add_parameters(copy.deepcopy(accepted_parameters))
+            journal.add_weights(copy.deepcopy(accepted_weights))
+            journal.add_distances(copy.deepcopy(distances))
             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)
@@ -1545,8 +1583,8 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1
 
             # print("INFO: Saving intermediate configuration to output journal.")
             if full_output == 1:
-                journal.add_parameters(accepted_parameters)
-                journal.add_weights(accepted_weights)
+                journal.add_parameters(copy.deepcopy(accepted_parameters))
+                journal.add_weights(copy.deepcopy(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)
@@ -1564,8 +1602,8 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1
         # Add anneal_parameter, number of final steps and final output to the journal
         # print("INFO: Saving final configuration to output journal.")
         if full_output == 0:
-            journal.add_parameters(accepted_parameters)
-            journal.add_weights(accepted_weights)
+            journal.add_parameters(copy.deepcopy(accepted_parameters))
+            journal.add_weights(copy.deepcopy(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)
@@ -1912,9 +1950,8 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1
 
             # print("INFO: Saving configuration to output journal.")
             if (full_output == 1 and aStep <= steps - 1) or (full_output == 0 and aStep == steps - 1):
-                journal.add_parameters(accepted_parameters)
+                journal.add_parameters(copy.deepcopy(accepted_parameters))
                 journal.add_weights(np.ones(shape=(len(accepted_parameters), 1)) * (1 / len(accepted_parameters)))
-
                 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)
@@ -2226,9 +2263,8 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1
 
             # print("INFO: Saving configuration to output journal.")
             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)
-
+                journal.add_parameters(copy.deepcopy(accepted_parameters))
+                journal.add_weights(copy.deepcopy(accepted_weights))
                 self.accepted_parameters_manager.update_broadcast(self.backend,
                                                                   accepted_parameters=accepted_parameters,
                                                                   accepted_weights=accepted_weights)
@@ -2544,9 +2580,9 @@ def sample(self, observations, steps, n_samples = 10000, n_samples_per_param = 1
             # print("INFO: Saving configuration to output journal.")
             if (full_output == 1 and aStep <= steps - 1) or (full_output == 0 and aStep == steps - 1):
                 self.accepted_parameters_manager.update_broadcast(self.backend, accepted_parameters=accepted_parameters)
-                journal.add_parameters(accepted_parameters)
-                journal.add_weights(accepted_weights)
-                journal.add_opt_values(accepted_y_sim)
+                journal.add_parameters(copy.deepcopy(accepted_parameters))
+                journal.add_weights(copy.deepcopy(accepted_weights))
+                journal.add_opt_values(copy.deepcopy(accepted_y_sim))
 
                 names_and_parameters = self._get_names_and_parameters()
                 journal.add_user_parameters(names_and_parameters)

diff --git a/abcpy/output.py b/abcpy/output.py
@@ -34,6 +34,7 @@ def __init__(self, type):
 
         self.parameters = []
         self.weights = []
+        self.distances = []
         self.opt_values = []
         self.configuration = {}
 
@@ -183,6 +184,39 @@ def add_weights(self, weights):
         if self._type == 1:
             self.weights.append(weights)
 
+    def get_distances(self, iteration=None):
+        """
+        Returns the distances from a sampling scheme.
+
+        For intermediate results, pass the iteration.
+
+        Parameters
+        ----------
+        iteration: int
+            specify the iteration for which to return distances
+        """
+
+        if iteration is None:
+            end = len(self.distances) - 1
+            return self.distances[end]
+        else:
+            return self.distances[iteration]
+
+    def add_distances(self, distances):
+        """
+        Saves provided distances by appending them to the journal. If type==0, old weights get overwritten.
+
+        Parameters
+        ----------
+        distances: numpy.array
+            vector containing n distances
+        """
+
+        if self._type == 0:
+            self.distances = [distances]
+
+        if self._type == 1:
+            self.distances.append(distances)
 
 
     def add_opt_values(self, opt_values):
@@ -202,7 +236,6 @@ def add_opt_values(self, opt_values):
             self.opt_values.append(opt_values)
 
 
-
     def save(self, filename):
         """
         Stores the journal to disk.

diff --git a/doc/source/conf.py b/doc/source/conf.py
@@ -25,13 +25,13 @@ class Mock(MagicMock):
     def __getattr__(cls, name):
             return MagicMock()
 
-MOCK_MODULES = ['numpy', 'pandas', 'glmnet', 'scipy', 'scipy.stats', 'scipy.special', 'scipy.optimize', 'sklearn', 'sklearn.covariance', 'findspark', 'coverage', 'numpy.random']
+MOCK_MODULES = ['numpy', 'pandas', 'glmnet', 'mpi4py', 'scipy', 'scipy.stats', 'scipy.special', 'scipy.optimize', 'sklearn', 'sklearn.covariance', 'findspark', 'coverage', 'numpy.random']
 sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES)
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
-sys.path.insert(0, os.path.abspath('.'))
+sys.path.insert(0, os.path.abspath('../..'))
 # -- General configuration ------------------------------------------------
 
 # If your documentation needs a minimal Sphinx version, state it here.
@@ -50,6 +50,7 @@ def __getattr__(cls, name):
 ]
 
 autodoc_member_order = 'bysource'
+autodoc_mock_imports = MOCK_MODULES
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']