From 41f7486a80bc8cf7f608a2846dc9ea1a41b73e12 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Fri, 25 Feb 2022 10:55:33 +0100
Subject: [PATCH 01/22] Fix acq_noise_var API (#397)

* Change ac_noise_var into float, int or dict

* Add checks for acq_noise_var

* Move acq_noise_var tests to acquisition.py.

* Fix lint.

* Add and modify tests for new acquisition API.

* Add CHANGELOG.

* Fix test_bo.py

* Fix doctrings that refer to AcquisitionBase.

* Remove check for None.

* Store parameter_names in GPyRegression

* Fix a bug on acquisition-tests.

* Modify documents to follow new API

* Fix a bug in docstring.

* Fix a bug in tests.
---
 CHANGELOG.rst                     |  2 ++
 docs/usage/BOLFI.rst              |  2 +-
 elfi/methods/bo/acquisition.py    | 40 +++++++++++++++++++++++++++----
 elfi/methods/bo/gpy_regression.py |  1 +
 elfi/methods/inference/bolfi.py   |  5 ++--
 elfi/methods/inference/bolfire.py |  4 ++--
 elfi/methods/inference/romc.py    |  4 ++--
 tests/unit/test_bo.py             | 13 ++++++++--
 8 files changed, 57 insertions(+), 14 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 00075782..0a435275 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,8 @@
 Changelog
 =========
 
+- Fix acq_noise_var-bug in acquisition.py. Influenced BOLFI.
+
 0.8.3 (2021-02-17)
 ------------------
 - Add a new inference method: BOLFIRE
diff --git a/docs/usage/BOLFI.rst b/docs/usage/BOLFI.rst
index df881533..dc805409 100644
--- a/docs/usage/BOLFI.rst
+++ b/docs/usage/BOLFI.rst
@@ -92,7 +92,7 @@ surface is updated after each batch (especially so if the noise is 0!).
 .. code:: ipython3
 
     bolfi = elfi.BOLFI(log_d, batch_size=1, initial_evidence=20, update_interval=10, 
-                       bounds={'t1':(-2, 2), 't2':(-1, 1)}, acq_noise_var=[0.1, 0.1], seed=seed)
+                       bounds={'t1':(-2, 2), 't2':(-1, 1)}, acq_noise_var=0.1, seed=seed)
 
 Sometimes you may have some samples readily available. You could then
 initialize the GP model with a dictionary of previous results by giving
diff --git a/elfi/methods/bo/acquisition.py b/elfi/methods/bo/acquisition.py
index a5c9d06c..5e959e95 100644
--- a/elfi/methods/bo/acquisition.py
+++ b/elfi/methods/bo/acquisition.py
@@ -62,15 +62,45 @@ def __init__(self,
         self.n_inits = int(n_inits)
         self.max_opt_iters = int(max_opt_iters)
         self.constraints = constraints
-
-        if noise_var is not None and np.asanyarray(noise_var).ndim > 1:
-            raise ValueError("Noise variance must be a float or 1d vector of variances "
-                             "for the different input dimensions.")
-        self.noise_var = noise_var
+        if noise_var is not None:
+            self._check_noise_var(noise_var)
+            self.noise_var = self._transform_noise_var(noise_var)
+        else:
+            self.noise_var = noise_var
         self.exploration_rate = exploration_rate
         self.random_state = np.random if seed is None else np.random.RandomState(seed)
         self.seed = 0 if seed is None else seed
 
+    def _check_noise_var(self, noise_var):
+        if isinstance(noise_var, dict):
+            if not set(noise_var) == set(self.model.parameter_names):
+                raise ValueError("Acquisition noise dictionary should contain all parameters.")
+
+            if not all(isinstance(x, (int, float)) for x in noise_var.values()):
+                raise ValueError("Acquisition noise dictionary values "
+                                 "should all be int or float.")
+
+            if any([x < 0 for x in noise_var.values()]):
+                raise ValueError("Acquisition noises values should all be "
+                                 "non-negative int or float.")
+
+        elif isinstance(noise_var, (int, float)):
+            if noise_var < 0:
+                raise ValueError("Acquisition noise should be non-negative int or float.")
+        else:
+            raise ValueError("Either acquisition noise is a float or "
+                             "it is a dictionary of floats defining "
+                             "variance for each parameter dimension.")
+
+    def _transform_noise_var(self, noise_var):
+        if isinstance(noise_var, (float, int)):
+            return noise_var
+
+        # return a sorted list of noise variances in the same order than
+        # parameter_names of the model
+        if isinstance(noise_var, dict):
+            return list(map(noise_var.get, self.model.parameter_names))
+
     def evaluate(self, x, t=None):
         """Evaluate the acquisition function at 'x'.
 
diff --git a/elfi/methods/bo/gpy_regression.py b/elfi/methods/bo/gpy_regression.py
index 8999a45f..8894f167 100644
--- a/elfi/methods/bo/gpy_regression.py
+++ b/elfi/methods/bo/gpy_regression.py
@@ -73,6 +73,7 @@ def __init__(self,
             raise ValueError("Keyword `bounds` must be a dictionary "
                              "`{'parameter_name': (lower, upper), ... }`")
 
+        self.parameter_names = parameter_names
         self.input_dim = input_dim
         self.bounds = bounds
 
diff --git a/elfi/methods/inference/bolfi.py b/elfi/methods/inference/bolfi.py
index 628444e3..6b814155 100644
--- a/elfi/methods/inference/bolfi.py
+++ b/elfi/methods/inference/bolfi.py
@@ -59,9 +59,9 @@ def __init__(self,
         target_model : GPyRegression, optional
         acquisition_method : Acquisition, optional
             Method of acquiring evidence points. Defaults to LCBSC.
-        acq_noise_var : float or np.array, optional
+        acq_noise_var : float or dict, optional
             Variance(s) of the noise added in the default LCBSC acquisition method.
-            If an array, should be 1d specifying the variance for each dimension.
+            If a dictionary, values should be float specifying the variance for each dimension.
         exploration_rate : float, optional
             Exploration rate of the acquisition method
         batch_size : int, optional
@@ -98,6 +98,7 @@ def __init__(self,
             self.target_model.update(params, precomputed[target_name])
 
         self.batches_per_acquisition = batches_per_acquisition or self.max_parallel_batches
+
         self.acquisition_method = acquisition_method or LCBSC(self.target_model,
                                                               prior=ModelPrior(
                                                                   self.model),
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index 9f41d286..65e052ab 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -58,9 +58,9 @@ def __init__(self,
             custom target_model is given.
         n_initial_evidence: int, optional
             Number of initial evidence.
-        acq_noise_var: float or np.ndarray, optional
+        acq_noise_var : float or dict, optional
             Variance(s) of the noise added in the default LCBSC acquisition method.
-            If an array, should be 1d specifying the variance for each dimension.
+            If a dictionary, values should be float specifying the variance for each dimension.
         exploration_rate: float, optional
             Exploration rate of the acquisition method.
         update_interval : int, optional
diff --git a/elfi/methods/inference/romc.py b/elfi/methods/inference/romc.py
index eb10ee90..92deef21 100644
--- a/elfi/methods/inference/romc.py
+++ b/elfi/methods/inference/romc.py
@@ -86,9 +86,9 @@ def __init__(self,
         target_model : GPyRegression, optional
         acquisition_method : Acquisition, optional
             Method of acquiring evidence points. Defaults to LCBSC.
-        acq_noise_var : float or np.array, optional
+        acq_noise_var : float or dict, optional
             Variance(s) of the noise added in the default LCBSC acquisition method.
-            If an array, should be 1d specifying the variance for each dimension.
+            If a dictionary, values should be float specifying the variance for each dimension.
         exploration_rate : float, optional
             Exploration rate of the acquisition method
         batch_size : int, optional
diff --git a/tests/unit/test_bo.py b/tests/unit/test_bo.py
index a8414d7b..14ead168 100644
--- a/tests/unit/test_bo.py
+++ b/tests/unit/test_bo.py
@@ -99,7 +99,7 @@ def test_acquisition():
     assert np.all((new[:, 1] >= bounds['b'][0]) & (new[:, 1] <= bounds['b'][1]))
 
     # check acquisition with separate variance for dimensions
-    acq_noise_var = np.random.uniform(0, 5, size=2)
+    acq_noise_var = {'a': 0.1, 'b': 0.5}
     t = 1
     acquisition_method = acquisition.LCBSC(target_model, noise_var=acq_noise_var)
     new = acquisition_method.acquire(n2, t=t)
@@ -111,12 +111,21 @@ def test_acquisition():
     acq_noise_cov = np.random.rand(n_params, n_params) * 0.5
     acq_noise_cov += acq_noise_cov.T
     acq_noise_cov += n_params * np.eye(n_params)
-    t = 1
     with pytest.raises(ValueError):
         acquisition.LCBSC(target_model, noise_var=acq_noise_cov)
 
+    # check acquisition with negative variances
+    acq_noise_var = -0.1
+    with pytest.raises(ValueError):
+        acquisition.LCBSC(target_model, noise_var=acq_noise_var)
+
+    acq_noise_var = {'a': 0.1, 'b': -0.1}
+    with pytest.raises(ValueError):
+        acquisition.LCBSC(target_model, noise_var=acq_noise_var)
+
     # test Uniform Acquisition
     t = 1
+    acq_noise_var = 0.1
     acquisition_method = acquisition.UniformAcquisition(target_model, noise_var=acq_noise_var)
     new = acquisition_method.acquire(n2, t=t)
     assert new.shape == (n2, n_params)

From a1aa456da6d5e46337e4fa2a4afe068bf2a6d8c6 Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Fri, 25 Feb 2022 11:24:56 +0100
Subject: [PATCH 02/22] Fix array check (#398)

* update array check

* remove obsolete float conversion

* update CHANGELOG

* Update CHANGELOG.rst

Co-authored-by: Henri Pesonen <henri.e.pesonen@gmail.com>
---
 CHANGELOG.rst                     | 1 +
 elfi/methods/inference/bolfire.py | 2 +-
 elfi/utils.py                     | 2 +-
 3 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 0a435275..9d833430 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Fix is_array in utils
 - Fix acq_noise_var-bug in acquisition.py. Influenced BOLFI.
 
 0.8.3 (2021-02-17)
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index 65e052ab..0cab6b16 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -412,7 +412,7 @@ def _get_observed_summary_values(self, model, summary_names):
 
     def _get_parameter_values(self, batch):
         """Return parameter values from a given batch."""
-        return {parameter_name: float(batch[parameter_name]) for parameter_name
+        return {parameter_name: batch[parameter_name] for parameter_name
                 in self.model.parameter_names}
 
     def _resolve_n_initial_evidence(self, n_initial_evidence):
diff --git a/elfi/utils.py b/elfi/utils.py
index 5b8b7263..5c8989a1 100644
--- a/elfi/utils.py
+++ b/elfi/utils.py
@@ -54,7 +54,7 @@ def args_to_tuple(*args):
 
 def is_array(output):
     """Check if `output` behaves as np.array (simple)."""
-    return hasattr(output, 'shape')
+    return hasattr(output, 'shape') and output.ndim > 0
 
 
 # NetworkX utils

From 173cebb584746ec8d86974b4fc05d90876dccccf Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Mon, 28 Feb 2022 13:49:55 +0100
Subject: [PATCH 03/22] Require a dictionary for sigma_proposals to be fed in
 metropolis. (#399)

* Require a dictionary for sigma_proposals to be fed in metropolis.

* Fix a bug in type annotation
---
 CHANGELOG.rst                     |  1 +
 elfi/methods/bo/acquisition.py    | 34 +++++++++++------------
 elfi/methods/inference/bolfi.py   | 15 +++++------
 elfi/methods/inference/bolfire.py | 11 +++-----
 elfi/methods/utils.py             | 45 ++++++++++++++++++++++++++++++-
 tests/unit/test_methods.py        |  5 +++-
 6 files changed, 74 insertions(+), 37 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9d833430..29825387 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Change sigma_proposals-input in metropolis from list to dict
 - Fix is_array in utils
 - Fix acq_noise_var-bug in acquisition.py. Influenced BOLFI.
 
diff --git a/elfi/methods/bo/acquisition.py b/elfi/methods/bo/acquisition.py
index 5e959e95..ed6b931f 100644
--- a/elfi/methods/bo/acquisition.py
+++ b/elfi/methods/bo/acquisition.py
@@ -8,6 +8,7 @@
 
 import elfi.methods.mcmc as mcmc
 from elfi.methods.bo.utils import minimize
+from elfi.methods.utils import resolve_sigmas
 
 logger = logging.getLogger(__name__)
 
@@ -493,7 +494,7 @@ class RandMaxVar(MaxVar):
     """
 
     def __init__(self, quantile_eps=.01, sampler='nuts', n_samples=50,
-                 limit_faulty_init=10, sigma_proposals_metropolis=None, *args, **opts):
+                 limit_faulty_init=10, sigma_proposals=None, *args, **opts):
         """Initialise RandMaxVar.
 
         Parameters
@@ -506,10 +507,9 @@ def __init__(self, quantile_eps=.01, sampler='nuts', n_samples=50,
             Length of the sampler's chain for obtaining the acquisitions.
         limit_faulty_init : int, optional
             Limit for the iterations used to obtain the sampler's initial points.
-        sigma_proposals_metropolis : array_like, optional
-            Standard deviation proposals for tuning the metropolis sampler.
-            For the default settings, the sigmas are set to the 1/10
-            of the parameter intervals' length.
+        sigma_proposals : dict, optional
+            Standard deviations for Gaussian proposals of each parameter for Metropolis
+            Markov Chain sampler. Defaults to 1/10 of surrogate model bound lengths.
 
         """
         super(RandMaxVar, self).__init__(quantile_eps, *args, **opts)
@@ -517,7 +517,10 @@ def __init__(self, quantile_eps=.01, sampler='nuts', n_samples=50,
         self.name_sampler = sampler
         self._n_samples = n_samples
         self._limit_faulty_init = limit_faulty_init
-        self._sigma_proposals_metropolis = sigma_proposals_metropolis
+        if self.name_sampler == 'metropolis':
+            self._sigma_proposals = resolve_sigmas(self.model.parameter_names,
+                                                   sigma_proposals,
+                                                   self.model.bounds)
 
     def acquire(self, n, t=None):
         """Acquire a batch of acquisition points.
@@ -576,18 +579,10 @@ def _evaluate_logpdf(theta):
 
             # Sampling the acquisition using the chosen sampler.
             if self.name_sampler == 'metropolis':
-                if self._sigma_proposals_metropolis is None:
-                    # Setting the default values of the sigma proposals to 1/10
-                    # of each parameters interval's length.
-                    sigma_proposals = []
-                    for bound in self.model.bounds:
-                        length_interval = bound[1] - bound[0]
-                        sigma_proposals.append(length_interval / 10)
-                    self._sigma_proposals_metropolis = sigma_proposals
                 samples = mcmc.metropolis(self._n_samples,
                                           theta_init,
                                           _evaluate_logpdf,
-                                          sigma_proposals=self._sigma_proposals_metropolis,
+                                          sigma_proposals=self._sigma_proposals,
                                           seed=self.seed)
             elif self.name_sampler == 'nuts':
                 samples = mcmc.nuts(self._n_samples,
@@ -637,7 +632,7 @@ class ExpIntVar(MaxVar):
 
     def __init__(self, quantile_eps=.01, integration='grid', d_grid=.2,
                  n_samples_imp=100, iter_imp=2, sampler='nuts', n_samples=2000,
-                 sigma_proposals_metropolis=None, *args, **opts):
+                 sigma_proposals=None, *args, **opts):
         """Initialise ExpIntVar.
 
         Parameters
@@ -662,8 +657,9 @@ def __init__(self, quantile_eps=.01, integration='grid', d_grid=.2,
         n_samples : int, optional
             Chain length for the sampler that generates the random numbers
             from the proposal distribution for IS.
-        sigma_proposals_metropolis : array_like, optional
-            Standard deviation proposals for tuning the metropolis sampler.
+        sigma_proposals : dict, optional
+            Standard deviations for Gaussian proposals of each parameter for Metropolis
+            Markov Chain sampler. Defaults to 1/10 of surrogate model bound lengths.
 
         """
         super(ExpIntVar, self).__init__(quantile_eps, *args, **opts)
@@ -681,7 +677,7 @@ def __init__(self, quantile_eps=.01, integration='grid', d_grid=.2,
                                          quantile_eps=self.quantile_eps,
                                          sampler=sampler,
                                          n_samples=n_samples,
-                                         sigma_proposals_metropolis=sigma_proposals_metropolis)
+                                         sigma_proposals=sigma_proposals)
         elif self._integration == 'grid':
             grid_param = [slice(b[0], b[1], d_grid) for b in self.model.bounds]
             self.points_int = np.mgrid[grid_param].reshape(len(self.model.bounds), -1).T
diff --git a/elfi/methods/inference/bolfi.py b/elfi/methods/inference/bolfi.py
index 6b814155..2b8351b8 100644
--- a/elfi/methods/inference/bolfi.py
+++ b/elfi/methods/inference/bolfi.py
@@ -17,7 +17,7 @@
 from elfi.methods.inference.parameter_inference import ParameterInference
 from elfi.methods.posteriors import BolfiPosterior
 from elfi.methods.results import BolfiSample, OptimizationResult
-from elfi.methods.utils import arr2d_to_batch, batch_to_arr2d, ceil_to_batch_size
+from elfi.methods.utils import arr2d_to_batch, batch_to_arr2d, ceil_to_batch_size, resolve_sigmas
 from elfi.model.extensions import ModelPrior
 
 logger = logging.getLogger(__name__)
@@ -490,9 +490,9 @@ def sample(self,
             Defaults to best evidence points.
         algorithm : string, optional
             Sampling algorithm to use. Currently 'nuts'(default) and 'metropolis' are supported.
-        sigma_proposals : np.array
+        sigma_proposals : dict, optional
             Standard deviations for Gaussian proposals of each parameter for Metropolis
-            Markov Chain sampler.
+            Markov Chain sampler. Defaults to 1/10 of surrogate model bound lengths.
         n_evidence : int
             If the regression model is not fitted yet, specify the amount of evidence
 
@@ -525,12 +525,9 @@ def sample(self,
         tasks_ids = []
         ii_initial = 0
         if algorithm == 'metropolis':
-            if sigma_proposals is None:
-                raise ValueError("Gaussian proposal standard deviations "
-                                 "have to be provided for Metropolis-sampling.")
-            elif sigma_proposals.shape[0] != self.target_model.input_dim:
-                raise ValueError("The length of Gaussian proposal standard "
-                                 "deviations must be n_params.")
+            sigma_proposals = resolve_sigmas(self.target_model.parameter_names,
+                                             sigma_proposals,
+                                             self.target_model.bounds)
 
         # sampling is embarrassingly parallel, so depending on self.client this may parallelize
         for ii in range(n_chains):
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index 0cab6b16..53b7feb9 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -14,7 +14,7 @@
 from elfi.methods.inference.parameter_inference import ParameterInference
 from elfi.methods.posteriors import BOLFIREPosterior
 from elfi.methods.results import BOLFIRESample
-from elfi.methods.utils import arr2d_to_batch, batch_to_arr2d
+from elfi.methods.utils import arr2d_to_batch, batch_to_arr2d, resolve_sigmas
 from elfi.model.elfi_model import ElfiModel, Summary
 from elfi.model.extensions import ModelPrior
 
@@ -269,12 +269,9 @@ def sample(self,
 
         # Check standard deviations of Gaussian proposals when using Metropolis-Hastings
         if algorithm == 'metropolis':
-            if sigma_proposals is None:
-                raise ValueError('Gaussian proposal standard deviations have '
-                                 'to be provided for Metropolis-sampling.')
-            elif sigma_proposals.shape[0] != self.target_model.input_dim:
-                raise ValueError('The length of Gaussian proposal standard '
-                                 'deviations must be n_params.')
+            sigma_proposals = resolve_sigmas(self.target_model.parameter_names,
+                                             sigma_proposals,
+                                             self.target_model.bounds)
 
         posterior = self.extract_result()
         warmup = warmup or n_samples // 2
diff --git a/elfi/methods/utils.py b/elfi/methods/utils.py
index 5d733c67..7d8d7299 100644
--- a/elfi/methods/utils.py
+++ b/elfi/methods/utils.py
@@ -2,7 +2,7 @@
 
 import logging
 from math import ceil
-from typing import Union
+from typing import Dict, List, Optional, Union
 
 import numpy as np
 import scipy.stats as ss
@@ -455,3 +455,46 @@ def flat_array_to_dict(names, arr):
     for ii, param_name in enumerate(names):
         param_dict[param_name] = np.expand_dims(arr[ii:ii + 1], 0)
     return param_dict
+
+
+def resolve_sigmas(parameter_names: List[str],
+                   sigma_proposals: Optional[Dict] = None,
+                   bounds: Optional[Dict] = None) -> List:
+    """Map dictionary of sigma_proposals into a list order as parameter_names.
+
+    Parameters
+    ----------
+    parameter_names: List[str]
+        names of the parameters
+    sigma_proposals: Dict
+        non-negative standard deviations for each dimension
+        {'parameter_name': float}
+    bounds : Dict, optional
+        the region where to estimate the posterior for each parameter in
+        model.parameters
+        `{'parameter_name':(lower, upper), ... }
+
+    Returns
+    -------
+    List
+       list of sigma_proposals in the same order than in parameter_names
+
+    """
+    if sigma_proposals is None:
+        sigma_proposals = []
+        for bound in bounds:
+            length_interval = bound[1] - bound[0]
+            sigma_proposals.append(length_interval / 10)
+    elif isinstance(sigma_proposals, dict):
+        errmsg = "sigma_proposals' keys have to be identical to " \
+                    "target_model.parameter_names."
+        if len(sigma_proposals) is not len(parameter_names):
+            raise ValueError(errmsg)
+        try:
+            sigma_proposals = [sigma_proposals[x] for x in parameter_names]
+        except ValueError:
+            print(parameter_names)
+    else:
+        raise ValueError("If provided, sigma_proposals need to be input as a dict.")
+
+    return sigma_proposals
diff --git a/tests/unit/test_methods.py b/tests/unit/test_methods.py
index 15baa626..b9aba870 100644
--- a/tests/unit/test_methods.py
+++ b/tests/unit/test_methods.py
@@ -100,7 +100,10 @@ def test_BOLFI_short(ma2, distribution_test):
     assert res_sampling_nuts.samples_array.shape[1] == 2
     assert len(res_sampling_nuts.samples_array) == n_samples // 2 * n_chains
 
-    res_sampling_metropolis = bolfi.sample(n_samples, n_chains=n_chains, algorithm='metropolis',sigma_proposals=np.ones(2))
+    res_sampling_metropolis = bolfi.sample(n_samples,
+                                           n_chains=n_chains,
+                                           algorithm='metropolis',
+                                           sigma_proposals={'t1': 0.2, 't2': 0.1})
     assert res_sampling_metropolis.samples_array.shape[1] == 2
     assert len(res_sampling_metropolis.samples_array) == n_samples // 2 * n_chains
 

From b7bb84cb751286437d2f1aead8c7be9d09e1763a Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Wed, 2 Mar 2022 14:27:25 +0100
Subject: [PATCH 04/22] Add additive acquisition cost (#400)

* initial commit

* update cost function

-standardise input and output shape
-include documentation

* add tests

* add type check

* fix style issues

* update CHANGELOG

* Minor style changes

Co-authored-by: Henri Pesonen <henri.e.pesonen@gmail.com>
---
 CHANGELOG.rst                     |  1 +
 elfi/methods/bo/acquisition.py    | 25 +++++++--------
 elfi/methods/bo/utils.py          | 53 +++++++++++++++++++++++++++++++
 elfi/methods/inference/bolfire.py |  6 +++-
 tests/unit/test_utils.py          | 20 +++++++++++-
 5 files changed, 90 insertions(+), 15 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 29825387..7ffd4e4e 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Add option to use additive acquisition cost in LCBSC
 - Change sigma_proposals-input in metropolis from list to dict
 - Fix is_array in utils
 - Fix acq_noise_var-bug in acquisition.py. Influenced BOLFI.
diff --git a/elfi/methods/bo/acquisition.py b/elfi/methods/bo/acquisition.py
index ed6b931f..f4062465 100644
--- a/elfi/methods/bo/acquisition.py
+++ b/elfi/methods/bo/acquisition.py
@@ -7,7 +7,7 @@
 import scipy.stats as ss
 
 import elfi.methods.mcmc as mcmc
-from elfi.methods.bo.utils import minimize
+from elfi.methods.bo.utils import CostFunction, minimize
 from elfi.methods.utils import resolve_sigmas
 
 logger = logging.getLogger(__name__)
@@ -223,7 +223,7 @@ class LCBSC(AcquisitionBase):
 
     """
 
-    def __init__(self, *args, delta=None, include_prior=False, **kwargs):
+    def __init__(self, *args, delta=None, additive_cost=None, **kwargs):
         """Initialize LCBSC.
 
         Parameters
@@ -231,8 +231,8 @@ def __init__(self, *args, delta=None, include_prior=False, **kwargs):
         delta: float, optional
             In between (0, 1). Default is 1/exploration_rate. If given, overrides the
             exploration_rate.
-        include_prior: bool, optional
-            If true, add negative log prior to model evaluations.
+        additive_cost: CostFunction, optional
+            Cost function output is added to the base acquisition value.
 
         """
         if delta is not None:
@@ -243,7 +243,10 @@ def __init__(self, *args, delta=None, include_prior=False, **kwargs):
         super(LCBSC, self).__init__(*args, **kwargs)
         self.name = 'lcbsc'
         self.label_fn = 'Confidence Bound'
-        self.include_prior = include_prior
+
+        if additive_cost is not None and not isinstance(additive_cost, CostFunction):
+            raise TypeError("Additive cost must be type CostFunction.")
+        self.additive_cost = additive_cost
 
     @property
     def delta(self):
@@ -272,10 +275,8 @@ def evaluate(self, x, t=None):
         """
         mean, var = self.model.predict(x, noiseless=True)
         value = mean - np.sqrt(self._beta(t) * var)
-        if self.include_prior:
-            # we use negative prior, since we minimize
-            negative_log_prior = -1 * self.prior.logpdf(x).reshape(-1, 1)
-            value += negative_log_prior
+        if self.additive_cost is not None:
+            value += self.additive_cost.evaluate(x)
         return value
 
     def evaluate_gradient(self, x, t=None):
@@ -295,10 +296,8 @@ def evaluate_gradient(self, x, t=None):
         mean, var = self.model.predict(x, noiseless=True)
         grad_mean, grad_var = self.model.predictive_gradients(x)
         value = grad_mean - 0.5 * grad_var * np.sqrt(self._beta(t) / var)
-        if self.include_prior:
-            # we use negative prior, since we minimize
-            grad_negative_log_prior = -1 * self.prior.gradient_logpdf(x).reshape(1, -1)
-            value += grad_negative_log_prior
+        if self.additive_cost is not None:
+            value += self.additive_cost.evaluate_gradient(x)
         return value
 
 
diff --git a/elfi/methods/bo/utils.py b/elfi/methods/bo/utils.py
index d695f137..f9442dbd 100644
--- a/elfi/methods/bo/utils.py
+++ b/elfi/methods/bo/utils.py
@@ -108,3 +108,56 @@ def minimize(fun,
         locs_out[i] = np.clip(locs_out[i], *bounds[i])
 
     return locs[ind_min], vals[ind_min]
+
+
+class CostFunction:
+    """Convenience class for modelling acquisition costs."""
+
+    def __init__(self, function, gradient, scale=1):
+        """Initialise CostFunction.
+
+        Parameters
+        ----------
+        function : callable
+            Function that returns cost function value.
+        gradient : callable
+            Function that returns cost function gradient.
+        scale : float, optional
+            Cost function is multiplied with scale.
+
+        """
+        self.function = function
+        self.gradient = gradient
+        self.scale = scale
+
+    def evaluate(self, x):
+        """Return cost function value evaluated at x.
+
+        Parameters
+        ----------
+        x : np.ndarray, shape: (input_dim,) or (n, input_dim)
+
+        Returns
+        -------
+        np.ndarray, shape: (n, 1)
+
+        """
+        x = np.atleast_2d(x)
+        n, input_dim = x.shape
+        return self.scale * self.function(x).reshape(n, 1)
+
+    def evaluate_gradient(self, x):
+        """Return cost function gradient evaluated at x.
+
+        Parameters
+        ----------
+        x : np.ndarray, shape: (input_dim,) or (n, input_dim)
+
+        Returns
+        -------
+        np.ndarray, shape: (n, input_dim)
+
+        """
+        x = np.atleast_2d(x)
+        n, input_dim = x.shape
+        return self.scale * self.gradient(x).reshape(n, input_dim)
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index 53b7feb9..17aa4881 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -11,6 +11,7 @@
 from elfi.loader import get_sub_seed
 from elfi.methods.bo.acquisition import LCBSC, AcquisitionBase
 from elfi.methods.bo.gpy_regression import GPyRegression
+from elfi.methods.bo.utils import CostFunction
 from elfi.methods.inference.parameter_inference import ParameterInference
 from elfi.methods.posteriors import BOLFIREPosterior
 from elfi.methods.results import BOLFIRESample
@@ -92,6 +93,9 @@ def __init__(self,
         # Initialize GP regression
         self.target_model = self._resolve_target_model(target_model)
 
+        # Define acquisition cost
+        self.cost = CostFunction(self.prior.logpdf, self.prior.gradient_logpdf, scale=-1)
+
         # Initialize BO
         self.n_initial_evidence = self._resolve_n_initial_evidence(n_initial_evidence)
         self.acquisition_method = self._resolve_acquisition_method(acquisition_method)
@@ -434,7 +438,7 @@ def _resolve_acquisition_method(self, acquisition_method):
                          noise_var=self.acq_noise_var,
                          exploration_rate=self.exploration_rate,
                          seed=self.seed,
-                         include_prior=True)
+                         additive_cost=self.cost)
         if isinstance(acquisition_method, AcquisitionBase):
             return acquisition_method
         raise TypeError('acquisition_method must be an instance of AcquisitionBase.')
diff --git a/tests/unit/test_utils.py b/tests/unit/test_utils.py
index 2a0eb8fe..92a531c2 100644
--- a/tests/unit/test_utils.py
+++ b/tests/unit/test_utils.py
@@ -6,7 +6,7 @@
 
 import elfi
 from elfi.examples.ma2 import get_model
-from elfi.methods.bo.utils import minimize, stochastic_optimization
+from elfi.methods.bo.utils import CostFunction, minimize, stochastic_optimization
 from elfi.methods.density_ratio_estimation import DensityRatioEstimation
 from elfi.methods.utils import (GMDistribution, normalize_weights, numgrad, numpy_to_python_type,
                                 sample_object_to_dict, weighted_sample_quantile, weighted_var)
@@ -282,3 +282,21 @@ def test_ratio_estimation(self):
 
         assert np.max(np.abs(test_w - test_w_estim)) < 0.1
         assert np.abs(np.max(test_w) - densratio.max_ratio()) < 0.1
+
+class TestCostFunction:
+
+    def test_evaluate(self):
+        def fun(x):
+            return x[0]**2 + (x[1] - 1)**4
+        
+        cost = CostFunction(elfi.tools.vectorize(fun), None, scale=10)
+        x = np.array([0.5, 0.5])
+        assert np.isclose(10 * fun(x), cost.evaluate(x))
+
+    def test_evaluate_gradient(self):
+        def grad(x):
+            return np.array([2 * x[0], 4 * (x[1] - 1)**3])
+
+        cost = CostFunction(None, elfi.tools.vectorize(grad), scale=10)
+        x = np.array([0.5, 0.5])
+        assert np.allclose(10 * grad(x), cost.evaluate_gradient(x))

From 576226e0b7b7c1cb416ed7b4ccc98405b2b78abb Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Mon, 14 Mar 2022 11:28:08 +0100
Subject: [PATCH 05/22] Update tox (#401)

* Update and fix tox.ini

* Add CHANGELOG

* Update MANIFEST.in
---
 CHANGELOG.rst | 1 +
 MANIFEST.in   | 2 ++
 tox.ini       | 3 +--
 3 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 7ffd4e4e..bb9c2206 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Update tox.ini
 - Add option to use additive acquisition cost in LCBSC
 - Change sigma_proposals-input in metropolis from list to dict
 - Fix is_array in utils
diff --git a/MANIFEST.in b/MANIFEST.in
index 0e007070..400460e1 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1 +1,3 @@
 recursive-include elfi/examples/cpp *.txt *.cpp Makefile
+include docs/description.rst
+include requirements.txt
diff --git a/tox.ini b/tox.ini
index d5bae418..2a206b08 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,6 +1,5 @@
 [tox]
-envlist = py27, py35, flake8
-; envlist = py26, py27, py33, py34, py35, flake8
+envlist = py36, py37, py38, flake8
 
 [testenv:flake8]
 basepython=python

From 189e94efb5f091baf5f04c7e2432179c0723991d Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Wed, 16 Mar 2022 13:39:34 +0100
Subject: [PATCH 06/22] Fix/kernel name bug (#402)

* Update CHANGELOG

* Use always surrogate parameter_names in BOLFI to avoid bugs.
---
 CHANGELOG.rst                   |  1 +
 elfi/methods/inference/bolfi.py | 28 +++++++++++++++++-----------
 2 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index bb9c2206..f6b3c78b 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Extract BO results using `target_model.parameter_names` from instead of `model.parameter_names`
 - Update tox.ini
 - Add option to use additive acquisition cost in LCBSC
 - Change sigma_proposals-input in metropolis from list to dict
diff --git a/elfi/methods/inference/bolfi.py b/elfi/methods/inference/bolfi.py
index 2b8351b8..29a200a9 100644
--- a/elfi/methods/inference/bolfi.py
+++ b/elfi/methods/inference/bolfi.py
@@ -93,7 +93,7 @@ def __init__(self,
         n_initial, precomputed = self._resolve_initial_evidence(
             initial_evidence)
         if precomputed is not None:
-            params = batch_to_arr2d(precomputed, self.parameter_names)
+            params = batch_to_arr2d(precomputed, self.target_model.parameter_names)
             n_precomputed = len(params)
             self.target_model.update(params, precomputed[target_name])
 
@@ -188,8 +188,11 @@ def extract_result(self):
         x_min, _ = stochastic_optimization(
             self.target_model.predict_mean, self.target_model.bounds, seed=self.seed)
 
-        batch_min = arr2d_to_batch(x_min, self.parameter_names)
-        outputs = arr2d_to_batch(self.target_model.X, self.parameter_names)
+        batch_min = arr2d_to_batch(x_min, self.target_model.parameter_names)
+        outputs = arr2d_to_batch(self.target_model.X, self.target_model.parameter_names)
+
+        # batch_min = arr2d_to_batch(x_min, self.parameter_names)
+        # outputs = arr2d_to_batch(self.target_model.X, self.parameter_names)
         outputs[self.target_name] = self.target_model.Y
 
         return OptimizationResult(
@@ -209,7 +212,7 @@ def update(self, batch, batch_index):
         super(BayesianOptimization, self).update(batch, batch_index)
         self.state['n_evidence'] += self.batch_size
 
-        params = batch_to_arr2d(batch, self.parameter_names)
+        params = batch_to_arr2d(batch, self.target_model.parameter_names)
         self._report_batch(batch_index, params, batch[self.target_name])
 
         optimize = self._should_optimize()
@@ -245,7 +248,7 @@ def prepare_new_batch(self, batch_index):
                 self.acq_batch_size, t=t)
 
         batch = arr2d_to_batch(
-            acquisition[:self.batch_size], self.parameter_names)
+            acquisition[:self.batch_size], self.target_model.parameter_names)
         self.state['acquisition'] = acquisition[self.batch_size:]
 
         return batch
@@ -311,7 +314,7 @@ def plot_state(self, **options):
         visin.draw_contour(
             gp.predict_mean,
             gp.bounds,
-            self.parameter_names,
+            self.target_model.parameter_names,
             title='GP target surface',
             points=gp.X,
             axes=f.axes[0],
@@ -342,7 +345,7 @@ def acq(x):
         visin.draw_contour(
             acq,
             gp.bounds,
-            self.parameter_names,
+            self.target_model.parameter_names,
             title='Acquisition surface',
             points=None,
             axes=f.axes[1],
@@ -363,7 +366,10 @@ def plot_discrepancy(self, axes=None, **kwargs):
         axes : np.array of plt.Axes
 
         """
-        return vis.plot_discrepancy(self.target_model, self.parameter_names, axes=axes, **kwargs)
+        return vis.plot_discrepancy(self.target_model,
+                                    self.target_model.parameter_names,
+                                    axes=axes,
+                                    **kwargs)
 
     def plot_gp(self, axes=None, resol=50, const=None, bounds=None, true_params=None, **kwargs):
         """Plot pairwise relationships as a matrix with parameters vs. discrepancy.
@@ -385,7 +391,7 @@ def plot_gp(self, axes=None, resol=50, const=None, bounds=None, true_params=None
         axes : np.array of plt.Axes
 
         """
-        return vis.plot_gp(self.target_model, self.parameter_names, axes,
+        return vis.plot_gp(self.target_model, self.target_model.parameter_names, axes,
                            resol, const, bounds, true_params, **kwargs)
 
 
@@ -574,7 +580,7 @@ def sample(self,
         print(
             "{} chains of {} iterations acquired. Effective sample size and Rhat for each "
             "parameter:".format(n_chains, n_samples))
-        for ii, node in enumerate(self.parameter_names):
+        for ii, node in enumerate(self.target_model.parameter_names):
             print(node, mcmc.eff_sample_size(chains[:, :, ii]),
                   mcmc.gelman_rubin_statistic(chains[:, :, ii]))
         self.target_model.is_sampling = False
@@ -582,7 +588,7 @@ def sample(self,
         return BolfiSample(
             method_name='BOLFI',
             chains=chains,
-            parameter_names=self.parameter_names,
+            parameter_names=self.target_model.parameter_names,
             warmup=warmup,
             threshold=float(posterior.threshold),
             n_sim=self.state['n_evidence'],

From 90a1e74fb927501baf48d7d42d62bcf8a484738b Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Thu, 17 Mar 2022 09:32:31 +0100
Subject: [PATCH 07/22] Fix/source of parameter names (#403)

* Add CHANGELOG.rst

* Use target_model.parameters_names instead of model.parameter_names.
---
 CHANGELOG.rst                     |  1 +
 elfi/methods/inference/bolfire.py | 12 ++++++------
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index f6b3c78b..81145faf 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Use `target_model.parameter_names` from instead of `model.parameter_names` in `BOLFIRE`
 - Extract BO results using `target_model.parameter_names` from instead of `model.parameter_names`
 - Update tox.ini
 - Add option to use additive acquisition cost in LCBSC
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index 17aa4881..d49eb9f3 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -127,7 +127,7 @@ def set_objective(self, n_evidence):
 
     def extract_result(self):
         """Extract the results from the current state."""
-        return BOLFIREPosterior(self.parameter_names,
+        return BOLFIREPosterior(self.target_model.parameter_names,
                                 self.target_model,
                                 self.prior,
                                 self.classifier_attributes)
@@ -158,7 +158,7 @@ def update(self, batch, batch_index):
 
         # BO part
         self.state['n_evidence'] += self.batch_size
-        parameter_values = batch_to_arr2d(batch, self.parameter_names)
+        parameter_values = batch_to_arr2d(batch, self.target_model.parameter_names)
         optimize = self._should_optimize()
         self.target_model.update(parameter_values, negative_log_ratio_value, optimize)
         if optimize:
@@ -182,7 +182,7 @@ def prepare_new_batch(self, batch_index):
 
         # Acquire parameter values from the acquisition function
         acquisition = self.acquisition_method.acquire(self.batch_size, t)
-        return arr2d_to_batch(acquisition, self.parameter_names)
+        return arr2d_to_batch(acquisition, self.target_model.parameter_names)
 
     def predict_log_ratio(self, X, y, X_obs):
         """Predict the log-ratio, i.e, logarithm of likelihood / marginal.
@@ -335,7 +335,7 @@ def sample(self,
 
         logger.info(f'{n_chains} chains of {n_samples} iterations acquired. '
                     'Effective sample size and Rhat for each parameter:')
-        for ii, node in enumerate(self.parameter_names):
+        for ii, node in enumerate(self.target_model.parameter_names):
             logger.info(f'{node} {mcmc.eff_sample_size(chains[:, :, ii])} '
                         f'{mcmc.gelman_rubin_statistic(chains[:, :, ii])}')
 
@@ -343,7 +343,7 @@ def sample(self,
 
         return BOLFIRESample(method_name='BOLFIRE',
                              chains=chains,
-                             parameter_names=self.parameter_names,
+                             parameter_names=self.target_model.parameter_names,
                              warmup=warmup,
                              n_sim=self.state['n_sim'],
                              seed=self.seed,
@@ -414,7 +414,7 @@ def _get_observed_summary_values(self, model, summary_names):
     def _get_parameter_values(self, batch):
         """Return parameter values from a given batch."""
         return {parameter_name: batch[parameter_name] for parameter_name
-                in self.model.parameter_names}
+                in self.target_model.parameter_names}
 
     def _resolve_n_initial_evidence(self, n_initial_evidence):
         """Resolve number of initial evidence."""

From 155a127cbc3c3575df0fd6cb86a0bd85e33aeebd Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Thu, 24 Mar 2022 09:46:21 +0100
Subject: [PATCH 08/22] Update BOLFIRE to run simulations in batches (#404)

* use batch system to run simulations

* clean

-rearranged new functions and added docstrings
-removed unused functions and updated unit tests

* update CHANGELOG

* fix typo

* fix random state setting

* update likelihood construction
---
 CHANGELOG.rst                     |   1 +
 elfi/methods/inference/bolfire.py | 119 ++++++++++++++++++------------
 tests/unit/test_bolfire_unit.py   |   6 +-
 3 files changed, 74 insertions(+), 52 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 81145faf..711cdfe9 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Use batch system to run simulations in BOLFIRE
 - Use `target_model.parameter_names` from instead of `model.parameter_names` in `BOLFIRE`
 - Extract BO results using `target_model.parameter_names` from instead of `model.parameter_names`
 - Update tox.ini
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index d49eb9f3..2dcb196c 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -107,6 +107,11 @@ def __init__(self,
         # Initialize classifier attributes list
         self.classifier_attributes = []
 
+        # Initialize data collection
+        self._likelihood = np.zeros((self.n_training_data, self.marginal.shape[1]))
+        self._random_state = np.random.RandomState(self.seed)
+        self._init_round()
+
     @property
     def n_evidence(self):
         """Return the number of acquired evidence points."""
@@ -123,7 +128,7 @@ def set_objective(self, n_evidence):
         """
         if n_evidence < self.n_evidence:
             logger.warning('Requesting less evidence than there already exists.')
-        self.objective['n_sim'] = n_evidence
+        self.objective['n_batches'] = n_evidence * int(self.n_training_data / self.batch_size)
 
     def extract_result(self):
         """Extract the results from the current state."""
@@ -144,25 +149,12 @@ def update(self, batch, batch_index):
             Index of batch.
 
         """
-        # Update the inference state
-        self.state['n_batches'] += 1
-        self.state['n_sim'] += self.batch_size * self.n_training_data
-
-        # Predict log-ratio
-        likelihood = self._generate_likelihood(self._get_parameter_values(batch))
-        X, y = self._generate_training_data(likelihood, self.marginal)
-        negative_log_ratio_value = -1 * self.predict_log_ratio(X, y, self.observed)
-
-        # Update classifier attributes list
-        self.classifier_attributes += [self.classifier.attributes]
+        super(BOLFIRE, self).update(batch, batch_index)
 
-        # BO part
-        self.state['n_evidence'] += self.batch_size
-        parameter_values = batch_to_arr2d(batch, self.target_model.parameter_names)
-        optimize = self._should_optimize()
-        self.target_model.update(parameter_values, negative_log_ratio_value, optimize)
-        if optimize:
-            self.state['last_GP_update'] = self.target_model.n_evidence
+        self._merge_batch(batch)
+        if self._round_sim == self.n_training_data:
+            self._update_logratio_model()
+            self._init_round()
 
     def prepare_new_batch(self, batch_index):
         """Prepare values for a new batch.
@@ -176,13 +168,8 @@ def prepare_new_batch(self, batch_index):
         batch: dict
 
         """
-        t = batch_index - self.n_initial_evidence
-        if t < 0:  # Sample parameter values from the model priors
-            return
-
-        # Acquire parameter values from the acquisition function
-        acquisition = self.acquisition_method.acquire(self.batch_size, t)
-        return arr2d_to_batch(acquisition, self.target_model.parameter_names)
+        batch_parameters = np.repeat(self._params, self.batch_size, axis=0)
+        return arr2d_to_batch(batch_parameters, self.target_model.parameter_names)
 
     def predict_log_ratio(self, X, y, X_obs):
         """Predict the log-ratio, i.e, logarithm of likelihood / marginal.
@@ -360,7 +347,9 @@ def _resolve_model(self, model):
     def _resolve_n_training_data(self, n_training_data):
         """Resolve the size of training data to be used."""
         if isinstance(n_training_data, int) and n_training_data > 0:
-            return n_training_data
+            if n_training_data % self.batch_size == 0:
+                return n_training_data
+            raise ValueError('n_training_data must be a multiple of batch_size.')
         raise TypeError('n_training_data must be a positive int.')
 
     def _get_summary_names(self, model):
@@ -384,20 +373,7 @@ def _generate_marginal(self, seed_marginal=None):
         batch = self.model.generate(self.n_training_data,
                                     outputs=self.summary_names,
                                     seed=seed_marginal)
-        return np.column_stack([batch[summary_name] for summary_name in self.summary_names])
-
-    def _generate_likelihood(self, parameter_values):
-        """Generate likelihood data."""
-        batch = self.model.generate(self.n_training_data,
-                                    outputs=self.summary_names,
-                                    with_values=parameter_values)
-        return np.column_stack([batch[summary_name] for summary_name in self.summary_names])
-
-    def _generate_training_data(self, likelihood, marginal):
-        """Generate training data."""
-        X = np.vstack((likelihood, marginal))
-        y = np.concatenate((np.ones(likelihood.shape[0]), -1 * np.ones(marginal.shape[0])))
-        return X, y
+        return batch_to_arr2d(batch, self.summary_names)
 
     def _resolve_classifier(self, classifier):
         """Resolve classifier."""
@@ -411,11 +387,6 @@ def _get_observed_summary_values(self, model, summary_names):
         """Return observed values for summary statistics."""
         return np.column_stack([model[summary_name].observed for summary_name in summary_names])
 
-    def _get_parameter_values(self, batch):
-        """Return parameter values from a given batch."""
-        return {parameter_name: batch[parameter_name] for parameter_name
-                in self.target_model.parameter_names}
-
     def _resolve_n_initial_evidence(self, n_initial_evidence):
         """Resolve number of initial evidence."""
         if isinstance(n_initial_evidence, int) and n_initial_evidence >= 0:
@@ -443,8 +414,62 @@ def _resolve_acquisition_method(self, acquisition_method):
             return acquisition_method
         raise TypeError('acquisition_method must be an instance of AcquisitionBase.')
 
+    def _init_round(self):
+        """Initialize data collection round."""
+        self._round_sim = 0
+
+        # Set new parameter values
+        if self.n_evidence < self.n_initial_evidence:
+            # Sample parameter values from the model priors
+            self._params = self.prior.rvs(1, random_state=self._random_state)
+        else:
+            # Acquire parameter values from the acquisition function
+            t = self.n_evidence - self.n_initial_evidence
+            self._params = self.acquisition_method.acquire(1, t)
+
+    def _new_round(self, batch_index):
+        """Check whether batch_index starts a new data collection round."""
+        return (batch_index * self.batch_size) % self.n_training_data == 0
+
+    def _allow_submit(self, batch_index):
+        """Check whether batch_index can be prepared."""
+        # Do not prepare batches with new parameter values until the current round is finished
+        if self._new_round(batch_index) and self.batches.has_pending:
+            return False
+        else:
+            return super(BOLFIRE, self)._allow_submit(batch_index)
+
+    def _merge_batch(self, batch):
+        """Add batch to collected data."""
+        data = batch_to_arr2d(batch, self.summary_names)
+        self._likelihood[self._round_sim:self._round_sim + self.batch_size] = data
+        self._round_sim += self.batch_size
+
+    def _update_logratio_model(self):
+        """Calculate log-ratio based on collected data and update surrogate model."""
+        # Predict log-ratio
+        X, y = self._generate_training_data(self._likelihood, self.marginal)
+        negative_log_ratio_value = -1 * self.predict_log_ratio(X, y, self.observed)
+
+        # Update classifier attributes list
+        self.classifier_attributes += [self.classifier.attributes]
+
+        # BO part
+        self.state['n_evidence'] += 1
+        parameter_values = self._params
+        optimize = self._should_optimize()
+        self.target_model.update(parameter_values, negative_log_ratio_value, optimize)
+        if optimize:
+            self.state['last_GP_update'] = self.target_model.n_evidence
+
+    def _generate_training_data(self, likelihood, marginal):
+        """Generate training data."""
+        X = np.vstack((likelihood, marginal))
+        y = np.concatenate((np.ones(likelihood.shape[0]), -1 * np.ones(marginal.shape[0])))
+        return X, y
+
     def _should_optimize(self):
         """Check whether GP hyperparameters should be optimized."""
-        current = self.target_model.n_evidence + self.batch_size
+        current = self.target_model.n_evidence + 1
         next_update = self.state['last_GP_update'] + self.update_interval
         return current >= self.n_initial_evidence and current >= next_update
diff --git a/tests/unit/test_bolfire_unit.py b/tests/unit/test_bolfire_unit.py
index 854f65b8..98ef9980 100644
--- a/tests/unit/test_bolfire_unit.py
+++ b/tests/unit/test_bolfire_unit.py
@@ -52,12 +52,8 @@ def test_generate_marginal(bolfire_method):
     assert bolfire_method._generate_marginal().shape == (10, 10)
 
 
-def test_generate_likelihood(bolfire_method, parameter_values):
-    assert bolfire_method._generate_likelihood(parameter_values).shape == (10, 10)
-
-
 def test_generate_training_data(bolfire_method, parameter_values):
-    likelihood = bolfire_method._generate_likelihood(parameter_values)
+    likelihood = np.random.rand(10, 10)
     X, y = bolfire_method._generate_training_data(likelihood, bolfire_method.marginal)
     assert X.shape == (20, 10)
     assert y.shape == (20,)

From 2d3eb32df9249243023958b02ef554eeed9555f6 Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Thu, 24 Mar 2022 10:12:02 +0100
Subject: [PATCH 09/22] Update BOLFIRE posterior (#407)

* use negative log-posterior in minimisation

* make minimisation optional

* update CHANGELOG

* update tests

* fix style issue

Co-authored-by: Henri Pesonen <henri.e.pesonen@gmail.com>
---
 CHANGELOG.rst                    |  1 +
 elfi/methods/posteriors.py       | 71 +++++++++++++-------------------
 tests/functional/test_bolfire.py |  2 +-
 3 files changed, 31 insertions(+), 43 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 711cdfe9..df0d259a 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Make MAP estimates calculation in BOLFIRE optional and based on log-posterior
 - Use batch system to run simulations in BOLFIRE
 - Use `target_model.parameter_names` from instead of `model.parameter_names` in `BOLFIRE`
 - Extract BO results using `target_model.parameter_names` from instead of `model.parameter_names`
diff --git a/elfi/methods/posteriors.py b/elfi/methods/posteriors.py
index adb76dea..b035dc76 100644
--- a/elfi/methods/posteriors.py
+++ b/elfi/methods/posteriors.py
@@ -49,10 +49,6 @@ def __init__(self, model, threshold=None, prior=None, n_inits=10, max_opt_iters=
             for details. By default, the minimum value of discrepancy estimate mean is used.
         prior : ScipyLikeDistribution, optional
             By default uniform distribution within model bounds.
-        n_inits : int, optional
-            Number of initialization points in internal optimization.
-        max_opt_iters : int, optional
-            Maximum number of iterations performed in internal optimization.
         seed : int, optional
 
         """
@@ -267,8 +263,6 @@ def __init__(self,
                  model,
                  prior,
                  classifier_attributes,
-                 n_opt_inits=10,
-                 max_opt_iters=1000,
                  *args, **kwargs):
         """Initialize BOLFIREPosterior.
 
@@ -282,34 +276,13 @@ def __init__(self,
             Joint prior of the elfi model.
         classifier_attributes: list
             Classifier's attributes on each inference round.
-        n_opt_inits: int, optional
-            Number of initialization points in internal optimization.
-        max_opt_iters: int, optional
-            Maximum number of iterations performed in internal optimization.
 
         """
         self._parameter_names = parameter_names
         self._model = model
         self._prior = prior
-        self._n_opt_inits = n_opt_inits
-        self._max_opt_iters = max_opt_iters
         self._classifier_attributes = classifier_attributes
 
-        # compute map estimates
-        self._map_estimates = self._compute_map_estimates()
-
-    @property
-    def map_estimates(self):
-        """Return the maximum a posterior estimate for each parameter.
-
-        Returns
-        -------
-        OrderedDict
-
-        """
-        return OrderedDict([(parameter_name, self._map_estimates[i]) for i, parameter_name
-                            in enumerate(self._parameter_names)])
-
     @property
     def classifier_attributes(self):
         """Return the classifier's attributes."""
@@ -338,10 +311,6 @@ def pdf(self, x):
         """
         return np.exp(self.logpdf(x))
 
-    def _negative_pdf(self, x):
-        """Return the negative unnormalized posterior at x."""
-        return -1 * self.pdf(x)
-
     def logpdf(self, x):
         """Return the unnormalized log-posterior at x.
 
@@ -356,6 +325,10 @@ def logpdf(self, x):
         """
         return self._prior.logpdf(x).reshape(-1, 1) - self._model.predict_mean(x)
 
+    def _negative_logpdf(self, x):
+        """Return the negative unnormalized log-posterior at x."""
+        return -1 * self.logpdf(x)
+
     def gradient_pdf(self, x):
         """Return the gradient of the unnormalized posterior pdf at x.
 
@@ -370,10 +343,6 @@ def gradient_pdf(self, x):
         """
         return np.exp(self.logpdf(x)) * self.gradient_logpdf(x)
 
-    def _negative_gradient_pdf(self, x):
-        """Return the negative gradient of the unnormalized posterior pdf at x."""
-        return -1 * self.gradient_pdf(x)
-
     def gradient_logpdf(self, x):
         """Return the gradient of unnormalized log-posterior pdf at x.
 
@@ -389,17 +358,35 @@ def gradient_logpdf(self, x):
         return self._prior.gradient_logpdf(x).reshape(1, -1) \
             - self._model.predictive_gradient_mean(x)
 
-    def _compute_map_estimates(self):
-        """Return the maximum a posterior estimate for each parameter."""
+    def _negative_gradient_logpdf(self, x):
+        """Return the negative gradient of the unnormalized log-posterior pdf at x."""
+        return -1 * self.gradient_logpdf(x)
+
+    def compute_map_estimates(self, n_opt_inits=10, max_opt_iters=1000):
+        """Return the maximum a posterior estimate for each parameter.
+
+        Parameters
+        ----------
+        n_inits : int, optional
+            Number of initialization points in optimization.
+        max_opt_iters : int, optional
+            Maximum number of iterations performed in optimization.
+
+        Returns
+        -------
+        OrderedDict
+
+        """
         minimum_location, _ = minimize(
-            fun=self._negative_pdf,
+            fun=self._negative_logpdf,
             bounds=self._model.bounds,
-            grad=self._negative_gradient_pdf,
+            grad=self._negative_gradient_logpdf,
             prior=self._prior,
-            n_start_points=self._n_opt_inits,
-            maxiter=self._max_opt_iters
+            n_start_points=n_opt_inits,
+            maxiter=max_opt_iters
         )
-        return minimum_location
+        return OrderedDict([(param, minimum_location[i]) for i, param in
+                            enumerate(self._model.parameter_names)])
 
 
 class RomcPosterior:
diff --git a/tests/functional/test_bolfire.py b/tests/functional/test_bolfire.py
index e59871ee..b46cf611 100644
--- a/tests/functional/test_bolfire.py
+++ b/tests/functional/test_bolfire.py
@@ -90,7 +90,7 @@ def test_bolfire(true_param, seed):
     assert bolfire_method.n_evidence == n_evidence
 
     # check the map estimates
-    map_estimates = bolfire_posterior.map_estimates
+    map_estimates = bolfire_posterior.compute_map_estimates()
     assert np.abs(map_estimates['mu'] - true_param) <= 0.5
 
     # run sampling

From 3d5e595e4da646913143c5851e4937a02d0ec816 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Mon, 28 Mar 2022 15:11:04 +0200
Subject: [PATCH 10/22] Add arch description (#408)

* Add model details to ARCH-docstring

* Add CHANGELOG
---
 CHANGELOG.rst         |  1 +
 elfi/examples/arch.py | 17 ++++++++++++++++-
 2 files changed, 17 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index df0d259a..1b68dcd2 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Add docstring description to ARCH-model
 - Make MAP estimates calculation in BOLFIRE optional and based on log-posterior
 - Use batch system to run simulations in BOLFIRE
 - Use `target_model.parameter_names` from instead of `model.parameter_names` in `BOLFIRE`
diff --git a/elfi/examples/arch.py b/elfi/examples/arch.py
index ce86cbc0..6d5a2ca6 100644
--- a/elfi/examples/arch.py
+++ b/elfi/examples/arch.py
@@ -63,7 +63,22 @@ def get_model(n_obs=100, true_params=None, seed_obs=None, n_lags=5):
 
 
 def arch(t1, t2, n_obs=100, batch_size=1, random_state=None):
-    """Generate a sequence of samples from the ARCH(1) model.
+    r"""Generate a sequence of samples from the ARCH(1) regression model.
+
+    Autoregressive conditional heteroskedasticity (ARCH) sequence describes the variance
+    of the error term as a function of previous error terms.
+
+        x_i = t_1 x_{i-1} + \epsilon_i
+
+        \epsilon_i = w_i \sqrt{0.2 + t_2 \epsilon_{i-1}^2}
+
+    where w_i is white noise ~ N(0,1) independent of \epsilon_0 ~  N(0,1)
+
+    References
+    ----------
+
+    Engle, R.F. (1982). Autoregressive Conditional Heteroscedasticity with
+        Estimates of the Variance of United Kingdom Inflation. Econometrica, 50(4): 987-1007
 
     Parameters
     ----------

From 6ad69e5566c3628a799626f586d4feb8b4ae58cc Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Tue, 29 Mar 2022 10:04:42 +0200
Subject: [PATCH 11/22] Fix generate outputs (#409)

* fix nodes type

* add new test

* update CHANGELOG
---
 CHANGELOG.rst                 |  1 +
 elfi/model/elfi_model.py      |  2 +-
 tests/unit/test_elfi_model.py | 11 +++++++++++
 3 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 1b68dcd2..9d6af898 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Fix default outputs in generate
 - Add docstring description to ARCH-model
 - Make MAP estimates calculation in BOLFIRE optional and based on log-posterior
 - Use batch system to run simulations in BOLFIRE
diff --git a/elfi/model/elfi_model.py b/elfi/model/elfi_model.py
index 4d4259d1..846ea525 100644
--- a/elfi/model/elfi_model.py
+++ b/elfi/model/elfi_model.py
@@ -277,7 +277,7 @@ def generate(self, batch_size=1, outputs=None, with_values=None, seed=None):
 
         """
         if outputs is None:
-            outputs = self.source_net.nodes()
+            outputs = list(self.source_net.nodes())
         elif isinstance(outputs, str):
             outputs = [outputs]
         if not isinstance(outputs, list):
diff --git a/tests/unit/test_elfi_model.py b/tests/unit/test_elfi_model.py
index 3556c415..28068af5 100644
--- a/tests/unit/test_elfi_model.py
+++ b/tests/unit/test_elfi_model.py
@@ -19,6 +19,17 @@ def test_generate(ma2):
     assert res.ndim == 1
 
 
+@pytest.mark.usefixtures('with_all_clients')
+def test_generate_outputs(ma2):
+    n_gen = 10
+
+    res = ma2.generate(n_gen)
+
+    assert 'd' in res
+    assert res['d'].shape[0] == n_gen
+    assert res['d'].ndim == 1
+
+
 @pytest.mark.usefixtures('with_all_clients')
 def test_observed():
     true_params = [.6, .2]

From 3b92623c40db414eff38f56dd77f4aef7a43dc12 Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Tue, 29 Mar 2022 11:32:48 +0200
Subject: [PATCH 12/22] Fix simulator node observed (#410)

* add len check

* update CHANGELOG

Co-authored-by: Henri Pesonen <henri.e.pesonen@gmail.com>
---
 CHANGELOG.rst    | 1 +
 elfi/executor.py | 3 +++
 2 files changed, 4 insertions(+)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9d6af898..6c8be963 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Fix the observed property in simulator nodes
 - Fix default outputs in generate
 - Add docstring description to ARCH-model
 - Make MAP estimates calculation in BOLFIRE optional and based on log-posterior
diff --git a/elfi/executor.py b/elfi/executor.py
index 4f20fb6f..d281b9a1 100644
--- a/elfi/executor.py
+++ b/elfi/executor.py
@@ -107,6 +107,9 @@ def get_execution_order(cls, G):
         # Filter those output nodes who have an operation to run
         needed = tuple(sorted(node for node in output_nodes if 'operation' in G.nodes[node]))
 
+        if len(needed) == 0:
+            return []
+
         if needed not in cache:
             # Resolve the nodes that need to be executed in the graph
             nodes_to_execute = set(needed)

From f24d44ddf20ef1e3d5e792121a5b88e8527157a7 Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Wed, 30 Mar 2022 14:48:27 +0200
Subject: [PATCH 13/22] Update BOLFIRE inputs (#411)

* add optional input variable feature_names

* make training data size not optional

* update tests

* fix style issue

* update CHANGELOG
---
 CHANGELOG.rst                     |  1 +
 elfi/methods/inference/bolfire.py | 48 ++++++++++++++++++++++---------
 tests/functional/test_bolfire.py  | 10 +++----
 tests/unit/test_bolfire_unit.py   |  2 +-
 4 files changed, 41 insertions(+), 20 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 6c8be963..8f2e8831 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Add feature names as an optional input and make training data size a required input in BOLFIRE
 - Fix the observed property in simulator nodes
 - Fix default outputs in generate
 - Add docstring description to ARCH-model
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index 2dcb196c..6955867b 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -16,7 +16,7 @@
 from elfi.methods.posteriors import BOLFIREPosterior
 from elfi.methods.results import BOLFIRESample
 from elfi.methods.utils import arr2d_to_batch, batch_to_arr2d, resolve_sigmas
-from elfi.model.elfi_model import ElfiModel, Summary
+from elfi.model.elfi_model import ElfiModel, ObservableMixin, Summary
 from elfi.model.extensions import ModelPrior
 
 logger = logging.getLogger(__name__)
@@ -27,7 +27,8 @@ class BOLFIRE(ParameterInference):
 
     def __init__(self,
                  model,
-                 n_training_data=10,
+                 n_training_data,
+                 feature_names=None,
                  marginal=None,
                  seed_marginal=None,
                  classifier=None,
@@ -45,8 +46,10 @@ def __init__(self,
         ----------
         model: ElfiModel
             Elfi graph used by the algorithm.
-        n_training_data: int, optional
+        n_training_data: int
             Size of training data.
+        feature_names: str or list, optional
+            ElfiModel nodes used as features in classification. Default all Summary nodes.
         marginal: np.ndnarray, optional
             Marginal data.
         seed_marginal: int, optional
@@ -78,10 +81,10 @@ def __init__(self,
 
         # Initialize attributes
         self.n_training_data = self._resolve_n_training_data(n_training_data)
-        self.summary_names = self._get_summary_names(self.model)
+        self.feature_names = self._resolve_feature_names(self.model, feature_names)
         self.marginal = self._resolve_marginal(marginal, seed_marginal)
         self.classifier = self._resolve_classifier(classifier)
-        self.observed = self._get_observed_summary_values(self.model, self.summary_names)
+        self.observed = self._get_observed_feature_values(self.model, self.feature_names)
         self.prior = ModelPrior(self.model)
 
         # TODO: write resolvers for the attributes below
@@ -340,8 +343,6 @@ def _resolve_model(self, model):
         """Resolve a given elfi model."""
         if not isinstance(model, ElfiModel):
             raise ValueError('model must be an ElfiModel.')
-        if len(self._get_summary_names(model)) == 0:
-            raise NotImplementedError('model must have at least one Summary node.')
         return model
 
     def _resolve_n_training_data(self, n_training_data):
@@ -352,6 +353,27 @@ def _resolve_n_training_data(self, n_training_data):
             raise ValueError('n_training_data must be a multiple of batch_size.')
         raise TypeError('n_training_data must be a positive int.')
 
+    def _resolve_feature_names(self, model, feature_names):
+        """Resolve feature names to be used."""
+        if feature_names is None:
+            feature_names = self._get_summary_names(model)
+            if len(feature_names) == 0:
+                raise NotImplementedError('Could not resolve feature_names based on the model.')
+            logger.info('Using all summary statistics as features in classification.')
+            return feature_names
+        if isinstance(feature_names, str):
+            feature_names = [feature_names]
+        if isinstance(feature_names, list):
+            if len(feature_names) == 0:
+                raise ValueError('feature_names must include at least one item.')
+            for feature_name in feature_names:
+                if feature_name not in model.nodes:
+                    raise ValueError(f'Node \'{feature_name}\' not found in the model.')
+                if not isinstance(model[feature_name], ObservableMixin):
+                    raise TypeError(f'Node \'{feature_name}\' is not observable.')
+            return feature_names
+        raise TypeError('feature_names must be a string or a list of strings.')
+
     def _get_summary_names(self, model):
         """Return the names of summary statistics."""
         return [node for node in model.nodes if isinstance(model[node], Summary)
@@ -371,9 +393,9 @@ def _resolve_marginal(self, marginal, seed_marginal=None):
     def _generate_marginal(self, seed_marginal=None):
         """Generate marginal data."""
         batch = self.model.generate(self.n_training_data,
-                                    outputs=self.summary_names,
+                                    outputs=self.feature_names,
                                     seed=seed_marginal)
-        return batch_to_arr2d(batch, self.summary_names)
+        return batch_to_arr2d(batch, self.feature_names)
 
     def _resolve_classifier(self, classifier):
         """Resolve classifier."""
@@ -383,9 +405,9 @@ def _resolve_classifier(self, classifier):
             return classifier
         raise ValueError('classifier must be an instance of Classifier.')
 
-    def _get_observed_summary_values(self, model, summary_names):
-        """Return observed values for summary statistics."""
-        return np.column_stack([model[summary_name].observed for summary_name in summary_names])
+    def _get_observed_feature_values(self, model, feature_names):
+        """Return observed feature values."""
+        return np.column_stack([model[feature_name].observed for feature_name in feature_names])
 
     def _resolve_n_initial_evidence(self, n_initial_evidence):
         """Resolve number of initial evidence."""
@@ -441,7 +463,7 @@ def _allow_submit(self, batch_index):
 
     def _merge_batch(self, batch):
         """Add batch to collected data."""
-        data = batch_to_arr2d(batch, self.summary_names)
+        data = batch_to_arr2d(batch, self.feature_names)
         self._likelihood[self._round_sim:self._round_sim + self.batch_size] = data
         self._round_sim += self.batch_size
 
diff --git a/tests/functional/test_bolfire.py b/tests/functional/test_bolfire.py
index b46cf611..1d74ca00 100644
--- a/tests/functional/test_bolfire.py
+++ b/tests/functional/test_bolfire.py
@@ -45,17 +45,15 @@ def test_bolfire_init(true_param, seed):
     m = simple_gaussian_model(true_param, seed)
 
     # define the bolfire method
-    bolfire_method = elfi.BOLFIRE(model=m)
+    bolfire_method = elfi.BOLFIRE(model=m, n_training_data=10)
 
-    # check the default size of training data
-    assert bolfire_method.n_training_data == 10
     # check the size of mariginal data (should be the size of training data x number of summaries)
     assert bolfire_method.marginal.shape == (10, 10)
-    # check the summary names
-    assert bolfire_method.summary_names == [f'power_{i}' for i in range(10)]
+    # check the feature names
+    assert bolfire_method.feature_names == [f'power_{i}' for i in range(10)]
     # check the type of a default classifier
     assert isinstance(bolfire_method.classifier, LogisticRegression)
-    # check the lenght of observed summary values
+    # check the length of observed feature values
     assert len(bolfire_method.observed[0]) == 10
     # check the type of the prior
     assert isinstance(bolfire_method.prior, ModelPrior)
diff --git a/tests/unit/test_bolfire_unit.py b/tests/unit/test_bolfire_unit.py
index 98ef9980..8c75f79a 100644
--- a/tests/unit/test_bolfire_unit.py
+++ b/tests/unit/test_bolfire_unit.py
@@ -45,7 +45,7 @@ def parameter_values():
 @pytest.fixture
 def bolfire_method(true_param, seed):
     m = simple_gaussian_model(true_param, seed)
-    return elfi.BOLFIRE(m)
+    return elfi.BOLFIRE(m, 10)
 
 
 def test_generate_marginal(bolfire_method):

From 7f5cfcffb501bd414d7e1169db757b1125243fb9 Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Tue, 10 May 2022 10:35:17 +0200
Subject: [PATCH 14/22] Fix parameter order in model prior (#412)

* add control over parameter order in model prior

add parameter names as optional input variable in model prior and use that to control parameter order in the priors used in BOLFI and BOLFIRE

* update CHANGELOG
---
 CHANGELOG.rst                     |  1 +
 elfi/methods/inference/bolfi.py   |  7 ++++---
 elfi/methods/inference/bolfire.py |  2 +-
 elfi/model/extensions.py          | 18 +++++++++++++++---
 4 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 8f2e8831..e06517ad 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Add parameter names as an optional input in model prior and fix the parameter order in priors used in BOLFI and BOLFIRE
 - Add feature names as an optional input and make training data size a required input in BOLFIRE
 - Fix the observed property in simulator nodes
 - Fix default outputs in generate
diff --git a/elfi/methods/inference/bolfi.py b/elfi/methods/inference/bolfi.py
index 29a200a9..61855d6b 100644
--- a/elfi/methods/inference/bolfi.py
+++ b/elfi/methods/inference/bolfi.py
@@ -99,9 +99,9 @@ def __init__(self,
 
         self.batches_per_acquisition = batches_per_acquisition or self.max_parallel_batches
 
+        prior = ModelPrior(self.model, parameter_names=self.target_model.parameter_names)
         self.acquisition_method = acquisition_method or LCBSC(self.target_model,
-                                                              prior=ModelPrior(
-                                                                  self.model),
+                                                              prior=prior,
                                                               noise_var=acq_noise_var,
                                                               exploration_rate=exploration_rate,
                                                               seed=self.seed)
@@ -456,7 +456,8 @@ def extract_posterior(self, threshold=None):
             raise ValueError(
                 'Model is not fitted yet, please see the `fit` method.')
 
-        return BolfiPosterior(self.target_model, threshold=threshold, prior=ModelPrior(self.model))
+        prior = ModelPrior(self.model, parameter_names=self.target_model.parameter_names)
+        return BolfiPosterior(self.target_model, threshold=threshold, prior=prior)
 
     def sample(self,
                n_samples,
diff --git a/elfi/methods/inference/bolfire.py b/elfi/methods/inference/bolfire.py
index 6955867b..845a1898 100644
--- a/elfi/methods/inference/bolfire.py
+++ b/elfi/methods/inference/bolfire.py
@@ -85,7 +85,6 @@ def __init__(self,
         self.marginal = self._resolve_marginal(marginal, seed_marginal)
         self.classifier = self._resolve_classifier(classifier)
         self.observed = self._get_observed_feature_values(self.model, self.feature_names)
-        self.prior = ModelPrior(self.model)
 
         # TODO: write resolvers for the attributes below
         self.bounds = bounds
@@ -95,6 +94,7 @@ def __init__(self,
 
         # Initialize GP regression
         self.target_model = self._resolve_target_model(target_model)
+        self.prior = ModelPrior(self.model, parameter_names=self.target_model.parameter_names)
 
         # Define acquisition cost
         self.cost = CostFunction(self.prior.logpdf, self.prior.gradient_logpdf, scale=-1)
diff --git a/elfi/model/extensions.py b/elfi/model/extensions.py
index cd7ef21a..28b293fb 100644
--- a/elfi/model/extensions.py
+++ b/elfi/model/extensions.py
@@ -116,18 +116,30 @@ def name(this):
 # TODO: could use some optimization
 # TODO: support the case where some priors are multidimensional
 class ModelPrior:
-    """Construct a joint prior distribution over all the parameter nodes in `ElfiModel`."""
+    """Construct a joint prior distribution over all or selected parameter nodes in `ElfiModel`."""
 
-    def __init__(self, model):
+    def __init__(self, model, parameter_names=None):
         """Initialize a ModelPrior.
 
         Parameters
         ----------
         model : ElfiModel
+        parameter_names : list, optional
+            Parameters included in the prior and their order. Default model.parameter_names.
 
         """
         model = model.copy()
-        self.parameter_names = model.parameter_names
+
+        if parameter_names is None:
+            self.parameter_names = model.parameter_names
+        elif isinstance(parameter_names, list):
+            for param in parameter_names:
+                if param not in model.parameter_names:
+                    raise ValueError(f"Parameter \'{param}\' not found in model parameters.")
+            self.parameter_names = parameter_names
+        else:
+            raise ValueError("parameter_names must be a list of strings.")
+
         self.dim = len(self.parameter_names)
         self.client = Client()
 

From f67c2ebdcaf3c06452ad8c9d66937191340ce401 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Thu, 9 Jun 2022 10:57:41 +0200
Subject: [PATCH 15/22] Add option to disable observation noise from
 lotka-volterra. (#414)

* Add option to disable observation noise from lotka-volterra.

* Modify observation_noise-option logic.

* Simplify allowed true_parameters -option.

* Change logger.errors to ValueErrors.
---
 CHANGELOG.rst                   |  2 ++
 elfi/examples/lotka_volterra.py | 29 ++++++++++++++++++++++++++---
 2 files changed, 28 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index e06517ad..7d30f700 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,8 @@
 Changelog
 =========
 
+
+- Add boolean `observation_noise` option to `lotka_volterra`
 - Add parameter names as an optional input in model prior and fix the parameter order in priors used in BOLFI and BOLFIRE
 - Add feature names as an optional input and make training data size a required input in BOLFIRE
 - Fix the observed property in simulator nodes
diff --git a/elfi/examples/lotka_volterra.py b/elfi/examples/lotka_volterra.py
index 0d17ff44..384f5a17 100644
--- a/elfi/examples/lotka_volterra.py
+++ b/elfi/examples/lotka_volterra.py
@@ -141,9 +141,11 @@ def lotka_volterra(r1, r2, r3, prey_init=50, predator_init=100, sigma=0., n_obs=
     return stock_out
 
 
-def get_model(n_obs=50, true_params=None, seed_obs=None, **kwargs):
+def get_model(n_obs=50, true_params=None, observation_noise=False, seed_obs=None, **kwargs):
     """Return a complete Lotka-Volterra model in inference task.
 
+    Including observation noise to system is optional.
+
     Parameters
     ----------
     n_obs : int, optional
@@ -152,6 +154,8 @@ def get_model(n_obs=50, true_params=None, seed_obs=None, **kwargs):
         Parameters with which the observed data is generated.
     seed_obs : int, optional
         Seed for the observed data generation.
+    observation_noise : bool, optional
+        Whether or not add normal noise to observations.
 
     Returns
     -------
@@ -160,7 +164,24 @@ def get_model(n_obs=50, true_params=None, seed_obs=None, **kwargs):
     """
     logger = logging.getLogger()
     if true_params is None:
-        true_params = [1.0, 0.005, 0.6, 50, 100, 10.]
+        if observation_noise:
+            true_params = [1.0, 0.005, 0.6, 50, 100, 10.]
+        else:
+            true_params = [1.0, 0.005, 0.6, 50, 100, 0.]
+    else:
+        if observation_noise:
+            if len(true_params) != 6:
+                raise ValueError(
+                        "Option observation_noise = True."
+                        " Provide six input parameters."
+                        )
+        else:
+            if len(true_params) != 5:
+                raise ValueError(
+                        "Option observation_noise = False."
+                        " Provide five input parameters."
+                        )
+            true_params = true_params + [0]
 
     kwargs['n_obs'] = n_obs
     y_obs = lotka_volterra(*true_params, random_state=np.random.RandomState(seed_obs), **kwargs)
@@ -175,7 +196,9 @@ def get_model(n_obs=50, true_params=None, seed_obs=None, **kwargs):
     priors.append(elfi.Prior(ExpUniform, -6., 2., model=m, name='r3'))
     priors.append(elfi.Prior('poisson', 50, model=m, name='prey0'))
     priors.append(elfi.Prior('poisson', 100, model=m, name='predator0'))
-    priors.append(elfi.Prior(ExpUniform, np.log(0.5), np.log(50), model=m, name='sigma'))
+
+    if observation_noise:
+        priors.append(elfi.Prior(ExpUniform, np.log(0.5), np.log(50), model=m, name='sigma'))
 
     elfi.Simulator(sim_fn, *priors, observed=y_obs, name='LV')
     sumstats.append(elfi.Summary(partial(pick_stock, species=0), m['LV'], name='prey'))

From 5e3cddf721711bcae50b53872058e78886f3241a Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Thu, 9 Jun 2022 13:25:07 +0200
Subject: [PATCH 16/22] Add sumstats to lv (#416)

* Add option to disable observation noise from lotka-volterra.

* Modify observation_noise-option logic.

* Simplify allowed true_parameters -option.

* Add summary statistics to LV

* Add sumstats to LV

* Add CHANGELOG.rst

* Enable batch evaluation.
---
 CHANGELOG.rst                   |   2 +-
 elfi/examples/lotka_volterra.py | 100 +++++++++++++++++++++++---------
 2 files changed, 73 insertions(+), 29 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 7d30f700..64ed8e4b 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,7 +1,7 @@
 Changelog
 =========
 
-
+- Add summary statistics to Lotka-Volterra model
 - Add boolean `observation_noise` option to `lotka_volterra`
 - Add parameter names as an optional input in model prior and fix the parameter order in priors used in BOLFI and BOLFIRE
 - Add feature names as an optional input and make training data size a required input in BOLFIRE
diff --git a/elfi/examples/lotka_volterra.py b/elfi/examples/lotka_volterra.py
index 384f5a17..dffc59c7 100644
--- a/elfi/examples/lotka_volterra.py
+++ b/elfi/examples/lotka_volterra.py
@@ -188,22 +188,32 @@ def get_model(n_obs=50, true_params=None, observation_noise=False, seed_obs=None
 
     m = elfi.ElfiModel()
     sim_fn = partial(lotka_volterra, **kwargs)
-    priors = []
-    sumstats = []
-
-    priors.append(elfi.Prior(ExpUniform, -6., 2., model=m, name='r1'))
-    priors.append(elfi.Prior(ExpUniform, -6., 2., model=m, name='r2'))  # easily kills populations
-    priors.append(elfi.Prior(ExpUniform, -6., 2., model=m, name='r3'))
-    priors.append(elfi.Prior('poisson', 50, model=m, name='prey0'))
-    priors.append(elfi.Prior('poisson', 100, model=m, name='predator0'))
+    priors = [
+        elfi.Prior(ExpUniform, -6., 2., model=m, name='r1'),
+        elfi.Prior(ExpUniform, -6., 2., model=m, name='r2'),  # easily kills populations
+        elfi.Prior(ExpUniform, -6., 2., model=m, name='r3'),
+        elfi.Prior('poisson', 50, model=m, name='prey0'),
+        elfi.Prior('poisson', 100, model=m, name='predator0')
+    ]
 
     if observation_noise:
         priors.append(elfi.Prior(ExpUniform, np.log(0.5), np.log(50), model=m, name='sigma'))
 
     elfi.Simulator(sim_fn, *priors, observed=y_obs, name='LV')
-    sumstats.append(elfi.Summary(partial(pick_stock, species=0), m['LV'], name='prey'))
-    sumstats.append(elfi.Summary(partial(pick_stock, species=1), m['LV'], name='predator'))
-    elfi.Distance('sqeuclidean', *sumstats, name='d')
+
+    sumstats = [
+        elfi.Summary(partial(stock_mean, species=0), m['LV'], name='prey_mean'),
+        elfi.Summary(partial(stock_mean, species=1), m['LV'], name='pred_mean'),
+        elfi.Summary(partial(stock_log_variance, species=0), m['LV'], name='prey_log_var'),
+        elfi.Summary(partial(stock_log_variance, species=1), m['LV'], name='pred_log_var'),
+        elfi.Summary(partial(stock_autocorr, species=0, lag=1), m['LV'], name='prey_autocorr_1'),
+        elfi.Summary(partial(stock_autocorr, species=1, lag=1), m['LV'], name='pred_autocorr_1'),
+        elfi.Summary(partial(stock_autocorr, species=0, lag=2), m['LV'], name='prey_autocorr_2'),
+        elfi.Summary(partial(stock_autocorr, species=1, lag=2), m['LV'], name='pred_autocorr_2'),
+        elfi.Summary(stock_crosscorr, m['LV'], name='crosscorr')
+    ]
+
+    elfi.Distance('euclidean', *sumstats, name='d')
 
     logger.info("Generated %i observations with true parameters r1: %.1f, r2: %.3f, r3: %.1f, "
                 "prey0: %i, predator0: %i, sigma: %.1f.", n_obs, *true_params)
@@ -211,6 +221,57 @@ def get_model(n_obs=50, true_params=None, observation_noise=False, seed_obs=None
     return m
 
 
+def stock_mean(stock, species=0, mu=0, std=1):
+    stock = np.atleast_2d(stock[:, :, species])
+    mu_x = np.mean(stock, axis=1)
+
+    return (mu_x - mu) / std
+
+
+def stock_log_variance(stock, species=0, mu=0, std=1):
+    stock = np.atleast_2d(stock[:, :, species])
+    var_x = np.var(stock, axis=1, ddof=1)
+    log_x = np.log(var_x + 1)
+
+    return (log_x - mu) / std
+
+
+def stock_autocorr(stock, species=0, lag=1, mu=0, std=1):
+    stock = np.atleast_2d(stock[:, :, species])
+    n_obs = stock.shape[1]
+
+    mu_x = np.mean(stock, axis=1, keepdims=True)
+    std_x = np.std(stock, axis=1, ddof=1, keepdims=True)
+    sx = ((stock - np.repeat(mu_x, n_obs, axis=1)) / np.repeat(std_x, n_obs, axis=1))
+
+    sx_t = sx[:, lag:]
+    sx_s = sx[:, :-lag]
+
+    C = np.sum(sx_t * sx_s, axis=1) / (n_obs - 1)
+
+    return (C - mu) / std
+
+
+def stock_crosscorr(stock, mu=0, std=1):
+    n_obs = stock.shape[1]
+
+    x_preys = stock[:, :, 0]  # preys
+    x_preds = stock[:, :, 1]  # predators
+
+    mu_preys = np.mean(x_preys, axis=1, keepdims=True)
+    mu_preds = np.mean(x_preds, axis=1, keepdims=True)
+    std_preys = np.std(x_preys, axis=1, keepdims=True)
+    std_preds = np.std(x_preds, axis=1, keepdims=True)
+    s_preys = ((x_preys - np.repeat(mu_preys, n_obs, axis=1))
+               / np.repeat(std_preys, n_obs, axis=1))
+    s_preds = ((x_preds - np.repeat(mu_preds, n_obs, axis=1))
+               / np.repeat(std_preds, n_obs, axis=1))
+
+    C = np.sum(s_preys * s_preds, axis=1) / (n_obs - 1)
+
+    return (C - mu) / std
+
+
 class ExpUniform(elfi.Distribution):
     r"""Prior distribution for parameter.
 
@@ -258,20 +319,3 @@ def pdf(cls, x, a, b):
             p = np.where((x < np.exp(a)) | (x > np.exp(b)), 0, np.reciprocal(x))
             p /= (b - a)  # normalize
         return p
-
-
-def pick_stock(stock, species):
-    """Return the stock for single species.
-
-    Parameters
-    ----------
-    stock : np.array
-    species : int
-        0 for prey, 1 for predator.
-
-    Returns
-    -------
-    np.array
-
-    """
-    return stock[:, :, species]

From fda4c78bfb0afe3dacca7b52354c9abecca3e7d0 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Thu, 9 Jun 2022 13:59:59 +0200
Subject: [PATCH 17/22] Fix lint in arch (#417)

* Add option to disable observation noise from lotka-volterra.

* Modify observation_noise-option logic.

* Simplify allowed true_parameters -option.

* Add summary statistics to LV

* Add sumstats to LV

* Add CHANGELOG.rst

* Fix linting in arch.py
---
 CHANGELOG.rst                   | 1 +
 elfi/examples/arch.py           | 2 +-
 elfi/examples/lotka_volterra.py | 5 ++++-
 3 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 64ed8e4b..09d3eb62 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Fix linting in `arch.py`
 - Add summary statistics to Lotka-Volterra model
 - Add boolean `observation_noise` option to `lotka_volterra`
 - Add parameter names as an optional input in model prior and fix the parameter order in priors used in BOLFI and BOLFIRE
diff --git a/elfi/examples/arch.py b/elfi/examples/arch.py
index 6d5a2ca6..68240e5a 100644
--- a/elfi/examples/arch.py
+++ b/elfi/examples/arch.py
@@ -76,7 +76,6 @@ def arch(t1, t2, n_obs=100, batch_size=1, random_state=None):
 
     References
     ----------
-
     Engle, R.F. (1982). Autoregressive Conditional Heteroscedasticity with
         Estimates of the Variance of United Kingdom Inflation. Econometrica, 50(4): 987-1007
 
@@ -102,6 +101,7 @@ def arch(t1, t2, n_obs=100, batch_size=1, random_state=None):
     e = E(t2, n_obs, batch_size, random_state)
     for i in range(1, n_obs + 1):
         y[:, i] = t1 * y[:, i - 1] + e[:, i]
+
     return y[:, 1:]
 
 
diff --git a/elfi/examples/lotka_volterra.py b/elfi/examples/lotka_volterra.py
index dffc59c7..f5c929ee 100644
--- a/elfi/examples/lotka_volterra.py
+++ b/elfi/examples/lotka_volterra.py
@@ -222,6 +222,7 @@ def get_model(n_obs=50, true_params=None, observation_noise=False, seed_obs=None
 
 
 def stock_mean(stock, species=0, mu=0, std=1):
+    """Calculate the mean of the trajectory by species."""
     stock = np.atleast_2d(stock[:, :, species])
     mu_x = np.mean(stock, axis=1)
 
@@ -229,6 +230,7 @@ def stock_mean(stock, species=0, mu=0, std=1):
 
 
 def stock_log_variance(stock, species=0, mu=0, std=1):
+    """Calculate the log variance of the trajectory by species."""
     stock = np.atleast_2d(stock[:, :, species])
     var_x = np.var(stock, axis=1, ddof=1)
     log_x = np.log(var_x + 1)
@@ -237,13 +239,13 @@ def stock_log_variance(stock, species=0, mu=0, std=1):
 
 
 def stock_autocorr(stock, species=0, lag=1, mu=0, std=1):
+    """Calculate the autocorrelation of lag n of the trajectory by species."""
     stock = np.atleast_2d(stock[:, :, species])
     n_obs = stock.shape[1]
 
     mu_x = np.mean(stock, axis=1, keepdims=True)
     std_x = np.std(stock, axis=1, ddof=1, keepdims=True)
     sx = ((stock - np.repeat(mu_x, n_obs, axis=1)) / np.repeat(std_x, n_obs, axis=1))
-
     sx_t = sx[:, lag:]
     sx_s = sx[:, :-lag]
 
@@ -253,6 +255,7 @@ def stock_autocorr(stock, species=0, lag=1, mu=0, std=1):
 
 
 def stock_crosscorr(stock, mu=0, std=1):
+    """Calculate the cross correlation of the species trajectories."""
     n_obs = stock.shape[1]
 
     x_preys = stock[:, :, 0]  # preys

From e5d3191e49144a17b44e065c90e3b0ba61f9d390 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Fri, 10 Jun 2022 14:57:28 +0200
Subject: [PATCH 18/22] Reformat summary-output. (#418)

* Reformat summary-output.

* Add CHANGELOG.rst

* Bring back sample_means_summary

* Remove double definition of weighted sample quantiles

* Revert to older definition of sample quantiles.
---
 CHANGELOG.rst           |  1 +
 elfi/methods/results.py | 35 +++++++++++++++++++++++++++++++----
 2 files changed, 32 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 09d3eb62..35c166a8 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Reformat `summary()` for `Sample(ParameterInferenceResult)`
 - Fix linting in `arch.py`
 - Add summary statistics to Lotka-Volterra model
 - Add boolean `observation_noise` option to `lotka_volterra`
diff --git a/elfi/methods/results.py b/elfi/methods/results.py
index a4806157..632cea75 100644
--- a/elfi/methods/results.py
+++ b/elfi/methods/results.py
@@ -8,11 +8,12 @@
 import sys
 from collections import OrderedDict
 
+import matplotlib.pyplot as plt
 import numpy as np
-from matplotlib import pyplot as plt
 
 import elfi.visualization.visualization as vis
-from elfi.methods.utils import numpy_to_python_type, sample_object_to_dict
+from elfi.methods.utils import (numpy_to_python_type, sample_object_to_dict,
+                                weighted_sample_quantile)
 
 logger = logging.getLogger(__name__)
 
@@ -174,7 +175,7 @@ def summary(self):
             desc += "Threshold: {:.3g}\n".format(self.threshold)
         print(desc, end='')
         try:
-            self.sample_means_summary()
+            self.sample_summary()
         except TypeError:
             pass
 
@@ -182,7 +183,28 @@ def sample_means_summary(self):
         """Print a representation of sample means."""
         s = "Sample means: "
         s += ', '.join(["{}: {:.3g}".format(k, v) for k, v in self.sample_means.items()])
-        print(s)
+        print(s)        
+
+    def sample_summary(self):
+        """Print sample mean and 95% credible interval."""
+        print("{0:24} {1:18} {2:17} {3:5}".format("Parameter", "Mean", "2.5%", "97.5%"))
+        print(''.join([
+            "{0:10} "
+            "{1:18.3f} "
+            "{2:18.3f} "
+            "{3:18.3f}\n"
+            .format(k[:10] + ":", v[0], v[1], v[2])
+            for k, v in self.sample_means_and_95CIs.items()]))
+
+    @property
+    def sample_means_and_95CIs(self):
+        """Construct OrderedDict for mean and 95% credible interval."""
+        return OrderedDict(
+            [(k, (np.average(v, axis=0, weights=self.weights),
+                  weighted_sample_quantile(v, alpha=0.025, weights=self.weights),
+                  weighted_sample_quantile(v, alpha=0.975, weights=self.weights)))
+             for k, v in self.samples.items()]
+                            )
 
     @property
     def sample_means(self):
@@ -196,6 +218,11 @@ def sample_means(self):
         return OrderedDict([(k, np.average(v, axis=0, weights=self.weights))
                             for k, v in self.samples.items()])
 
+    def sample_quantiles(self, alpha=0.5):
+        """Evaluate weighted sample quantiles of sampled parameters."""
+        return OrderedDict([(k, weighted_sample_quantile(v, alpha=alpha, weights=self.weights))
+                            for k, v in self.samples.items()])
+
     @property
     def sample_means_array(self):
         """Evaluate weighted averages of sampled parameters.

From d5a3828583f2d05ace9693080fb3077253e8c28f Mon Sep 17 00:00:00 2001
From: uremes <uremes@users.noreply.github.com>
Date: Sat, 11 Jun 2022 14:46:17 +0200
Subject: [PATCH 19/22] Fix acquisition index in state plot (#420)

* fix plot_state

* update CHANGELOG
---
 CHANGELOG.rst                   | 1 +
 elfi/methods/inference/bolfi.py | 5 +++--
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 35c166a8..e5b61850 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Fix acquisition index in state plot
 - Reformat `summary()` for `Sample(ParameterInferenceResult)`
 - Fix linting in `arch.py`
 - Add summary statistics to Lotka-Volterra model
diff --git a/elfi/methods/inference/bolfi.py b/elfi/methods/inference/bolfi.py
index 61855d6b..ad7f9064 100644
--- a/elfi/methods/inference/bolfi.py
+++ b/elfi/methods/inference/bolfi.py
@@ -326,7 +326,7 @@ def plot_state(self, **options):
             if len(gp.X) > 1:
                 f.axes[1].scatter(*point, color='red')
 
-        displays = [gp._gp]
+        displays = [gp.instance]
 
         if options.get('interactive'):
             from IPython import display
@@ -338,8 +338,9 @@ def plot_state(self, **options):
         # Update
         visin._update_interactive(displays, options)
 
+        acq_index = self._get_acquisition_index(self.state['n_batches'])
         def acq(x):
-            return self.acquisition_method.evaluate(x, len(gp.X))
+            return self.acquisition_method.evaluate(x, acq_index)
 
         # Draw the acquisition surface
         visin.draw_contour(

From d7ebc33e70fc3c8e7a3d238b0cab9c0ff44dfee7 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Mon, 13 Jun 2022 10:19:52 +0200
Subject: [PATCH 20/22] Modify lv priors (#419)

* Change LV priors from poisson to normal approximations.

* Approximate poisson distributions with gaussians.

* Add CHANGELOG.rst
---
 CHANGELOG.rst                   |  1 +
 elfi/examples/lotka_volterra.py | 10 ++++++----
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index e5b61850..b51e9c6c 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Modify Lotka-Volterra model's priors as many methods do not support discrete random variables.
 - Fix acquisition index in state plot
 - Reformat `summary()` for `Sample(ParameterInferenceResult)`
 - Fix linting in `arch.py`
diff --git a/elfi/examples/lotka_volterra.py b/elfi/examples/lotka_volterra.py
index f5c929ee..3fe0e06e 100644
--- a/elfi/examples/lotka_volterra.py
+++ b/elfi/examples/lotka_volterra.py
@@ -74,8 +74,10 @@ def lotka_volterra(r1, r2, r3, prey_init=50, predator_init=100, sigma=0., n_obs=
 
     n_full = 20000
     stock = np.empty((batch_size, n_full, 2), dtype=np.int32)
-    stock[:, 0, 0] = prey_init
-    stock[:, 0, 1] = predator_init
+    # As we use approximate continuous priors for prey_init and
+    # predator_init, we'll round them down to closest integers
+    stock[:, 0, 0] = np.floor(prey_init)
+    stock[:, 0, 1] = np.floor(predator_init)
     stoichiometry = np.array([[1, 0], [-1, 1], [0, -1], [0, 0]], dtype=np.int32)
     times = np.empty((batch_size, n_full))
     times[:, 0] = 0
@@ -192,8 +194,8 @@ def get_model(n_obs=50, true_params=None, observation_noise=False, seed_obs=None
         elfi.Prior(ExpUniform, -6., 2., model=m, name='r1'),
         elfi.Prior(ExpUniform, -6., 2., model=m, name='r2'),  # easily kills populations
         elfi.Prior(ExpUniform, -6., 2., model=m, name='r3'),
-        elfi.Prior('poisson', 50, model=m, name='prey0'),
-        elfi.Prior('poisson', 100, model=m, name='predator0')
+        elfi.Prior('normal', 50, np.sqrt(50), model=m, name='prey0'),
+        elfi.Prior('normal', 100, np.sqrt(100), model=m, name='predator0')
     ]
 
     if observation_noise:

From f99c9159cbdc38a70175a79e9512cfffd710d7c4 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Mon, 13 Jun 2022 14:13:29 +0200
Subject: [PATCH 21/22] Fix linting (#422)

---
 elfi/methods/inference/bolfi.py | 1 +
 elfi/methods/results.py         | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/elfi/methods/inference/bolfi.py b/elfi/methods/inference/bolfi.py
index ad7f9064..5b7c067f 100644
--- a/elfi/methods/inference/bolfi.py
+++ b/elfi/methods/inference/bolfi.py
@@ -339,6 +339,7 @@ def plot_state(self, **options):
         visin._update_interactive(displays, options)
 
         acq_index = self._get_acquisition_index(self.state['n_batches'])
+
         def acq(x):
             return self.acquisition_method.evaluate(x, acq_index)
 
diff --git a/elfi/methods/results.py b/elfi/methods/results.py
index 632cea75..8f3c5284 100644
--- a/elfi/methods/results.py
+++ b/elfi/methods/results.py
@@ -183,7 +183,7 @@ def sample_means_summary(self):
         """Print a representation of sample means."""
         s = "Sample means: "
         s += ', '.join(["{}: {:.3g}".format(k, v) for k, v in self.sample_means.items()])
-        print(s)        
+        print(s)
 
     def sample_summary(self):
         """Print sample mean and 95% credible interval."""

From 6b1b17aab3e26120e9d9c156aacfb04bb4811421 Mon Sep 17 00:00:00 2001
From: Henri Pesonen <henri.e.pesonen@gmail.com>
Date: Mon, 13 Jun 2022 14:59:41 +0200
Subject: [PATCH 22/22] Bump to version 0.8.4 (#423)

---
 CHANGELOG.rst    | 2 ++
 README.md        | 2 +-
 elfi/__init__.py | 2 +-
 3 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index b51e9c6c..6813da7a 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,8 @@
 Changelog
 =========
 
+0.8.4 (2021-06-13)
+------------------
 - Modify Lotka-Volterra model's priors as many methods do not support discrete random variables.
 - Fix acquisition index in state plot
 - Reformat `summary()` for `Sample(ParameterInferenceResult)`
diff --git a/README.md b/README.md
index 6a3b4f42..5b1c2019 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-**Version 0.8.3 released!** See the [CHANGELOG](CHANGELOG.rst) and [notebooks](https://github.com/elfi-dev/notebooks).
+**Version 0.8.4 released!** See the [CHANGELOG](CHANGELOG.rst) and [notebooks](https://github.com/elfi-dev/notebooks).
 
 <img src="https://raw.githubusercontent.com/elfi-dev/elfi/dev/docs/logos/elfi_logo_text_nobg.png" width="200" />
 
diff --git a/elfi/__init__.py b/elfi/__init__.py
index 293569e2..9bd32190 100644
--- a/elfi/__init__.py
+++ b/elfi/__init__.py
@@ -31,4 +31,4 @@
 __email__ = 'elfi-support@hiit.fi'
 
 # make sure __version_ is on the last non-empty line (read by setup.py)
-__version__ = '0.8.3'
+__version__ = '0.8.4'