From 8c61b438ea7175c749554dba53947850a96713da Mon Sep 17 00:00:00 2001
From: Henri Vuollekoski <henri.vuollekoski@gmail.com>
Date: Tue, 3 Apr 2018 12:02:48 +0300
Subject: [PATCH 1/6] Fix Travis (#258)

* Fix Travis, change in Homebrew

* Relax np.array_equal->np.allclose
---
 .travis.yml                          | 7 ++++---
 tests/functional/test_consistency.py | 9 +++++----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index e879401f..e639c446 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -12,9 +12,10 @@ matrix:
     - mkdir -p /Users/travis/.matplotlib
     - "echo 'backend: TkAgg' > /Users/travis/.matplotlib/matplotlibrc"
     - brew update
-    - brew install python3
-    - virtualenv env -p python3
-    - source env/bin/activate
+    - brew upgrade python
+    - pip3 install virtualenv
+    - virtualenv py3env -p python3
+    - source py3env/bin/activate
 
 cache: pip
 
diff --git a/tests/functional/test_consistency.py b/tests/functional/test_consistency.py
index 1908667a..0973ff71 100644
--- a/tests/functional/test_consistency.py
+++ b/tests/functional/test_consistency.py
@@ -11,12 +11,12 @@
 def check_consistent_sample(sample, sample_diff, sample_same):
     assert not np.array_equal(sample.outputs['t1'], sample_diff.outputs['t1'])
 
-    assert np.array_equal(sample.outputs['t1'], sample_same.outputs['t1'])
-    assert np.array_equal(sample.outputs['t2'], sample_same.outputs['t2'])
+    assert np.allclose(sample.outputs['t1'], sample_same.outputs['t1'])
+    assert np.allclose(sample.outputs['t2'], sample_same.outputs['t2'])
 
     # BOLFI does not have d in its outputs
     if 'd' in sample.outputs:
-        assert np.array_equal(sample.outputs['d'], sample_same.outputs['d'])
+        assert np.allclose(sample.outputs['d'], sample_same.outputs['d'])
 
 
 @pytest.mark.usefixtures('with_all_clients')
@@ -101,7 +101,8 @@ def test_bo(ma2):
     check_consistent_sample(res, res_diff, res_same)
 
     assert not np.array_equal(res.x_min, res_diff.x_min)
-    assert np.array_equal(res.x_min, res_same.x_min)
+    assert np.allclose(res.x_min['t1'], res_same.x_min['t1'], atol=1e-07)
+    assert np.allclose(res.x_min['t2'], res_same.x_min['t2'], atol=1e-07)
 
 
 # TODO: skipped in travis due to NUTS initialization failing too often. Should be fixed.

From a35290ae022132e9c5b02be6d84ead457324f904 Mon Sep 17 00:00:00 2001
From: Henri Vuollekoski <henri.vuollekoski@gmail.com>
Date: Thu, 5 Apr 2018 13:55:31 +0300
Subject: [PATCH 2/6] Some updates in Docs (#259)

* Clarifications in tutorial notebook

* Explain vectorization

* Address comments to PR
---
 docs/faq.rst            | 27 +++++++++++++++++
 docs/usage/tutorial.rst | 66 ++++++++++++++++++++---------------------
 2 files changed, 60 insertions(+), 33 deletions(-)

diff --git a/docs/faq.rst b/docs/faq.rst
index 464f2b87..ed3a9a5c 100644
--- a/docs/faq.rst
+++ b/docs/faq.rst
@@ -10,3 +10,30 @@ produces outputs from the interval (1, 3).*
 their definitions. There the uniform distribution uses the location/scale definition, so
 the first argument defines the starting point of the interval and the second its length.
 
+.. _vectorization:
+
+*Q: What is vectorization in ELFI?*
+
+**A**: Looping is relatively inefficient in Python, and so whenever possible, you should *vectorize*
+your operations_. This means that repetitive computations are performed on a batch of data using
+precompiled libraries (typically NumPy_), which effectively runs the loops in faster, compiled C-code.
+Due to the potentially huge saving in CPU-time, operations including user-code are by default assumed to
+be vectorized in ELFI. This must be accounted for.
+
+.. _operations: good-to-know.html#operations
+.. _NumPy: http://www.numpy.org/
+
+For example, imagine you have a simulator that depends on a scalar parameter and produces a vector of 5
+values. When this is used in ELFI with ``batch_size`` set to 1000, ELFI draws 1000 values from the
+parameter's prior distribution and gives this *vector* to the simulator. Ideally, the simulator should
+efficiently process all 1000 parameter cases in one go and output a numpy array of shape (1000, 5). In
+ELFI, the length (i.e. the first dimension) of all output arrays should equal ``batch_size`` **even if
+it is 1**. Note that because of this the evaluation of summary statistics, distances etc. should
+bypass the first dimension (e.g. with NumPy functions using ``axis=1`` in this case).
+
+See ``elfi.examples`` for tips on how to vectorize simulators and work with ELFI. In case you are
+unable to vectorize your simulator, you can use `elfi.tools.vectorize`_ to mimic
+vectorized behaviour, though without the performance benefits.
+
+.. _`elfi.tools.vectorize`: api.html#elfi.tools.vectorize
+
diff --git a/docs/usage/tutorial.rst b/docs/usage/tutorial.rst
index 444a97c5..f9ccadf9 100644
--- a/docs/usage/tutorial.rst
+++ b/docs/usage/tutorial.rst
@@ -18,10 +18,9 @@ settings.
     
     import numpy as np
     import scipy.stats
-    import matplotlib
     import matplotlib.pyplot as plt
     import logging
-    logging.basicConfig(level=logging.INFO)
+    logging.basicConfig(level=logging.INFO)  # sometimes this is required to enable logging inside Jupyter
     
     %matplotlib inline
     %precision 2
@@ -251,7 +250,7 @@ a DAG.
 
 
 
-.. note:: You will need the Graphviz_ software as well as the graphviz `Python package`_ (https://pypi.python.org/pypi/graphviz) for drawing this. The software is already installed in many unix-like OS.
+.. note:: You will need the Graphviz_ software as well as the graphviz `Python package`_ (https://pypi.python.org/pypi/graphviz) for drawing this.
 
 .. _Graphviz: http://www.graphviz.org
 .. _`Python package`: https://pypi.python.org/pypi/graphviz
@@ -396,8 +395,8 @@ time is spent in drawing.
 
 .. parsed-literal::
 
-    CPU times: user 2.28 s, sys: 165 ms, total: 2.45 s
-    Wall time: 2.45 s
+    CPU times: user 1.6 s, sys: 166 ms, total: 1.77 s
+    Wall time: 1.76 s
 
 
 The ``sample`` method returns a ``Sample`` object, which contains
@@ -452,8 +451,8 @@ as long as it takes to generate the requested number of samples.
 
 .. parsed-literal::
 
-    CPU times: user 222 ms, sys: 40.3 ms, total: 263 ms
-    Wall time: 261 ms
+    CPU times: user 198 ms, sys: 35.5 ms, total: 233 ms
+    Wall time: 231 ms
     Method: Rejection
     Number of samples: 1000
     Number of simulations: 40000
@@ -497,9 +496,9 @@ been reached or a maximum of one second of time has been used.
 
     Method: Rejection
     Number of samples: 1000
-    Number of simulations: 190000
-    Threshold: 0.0855
-    Sample means: t1: 0.561, t2: 0.218
+    Number of simulations: 180000
+    Threshold: 0.088
+    Sample means: t1: 0.561, t2: 0.221
     
 
 
@@ -547,8 +546,8 @@ in our model:
 
 .. parsed-literal::
 
-    CPU times: user 5.26 s, sys: 37.1 ms, total: 5.3 s
-    Wall time: 5.3 s
+    CPU times: user 5.01 s, sys: 60.9 ms, total: 5.07 s
+    Wall time: 5.09 s
 
 
 
@@ -558,8 +557,8 @@ in our model:
     Method: Rejection
     Number of samples: 1000
     Number of simulations: 1000000
-    Threshold: 0.036
-    Sample means: t1: 0.561, t2: 0.227
+    Threshold: 0.0363
+    Sample means: t1: 0.554, t2: 0.216
 
 
 
@@ -580,8 +579,8 @@ anything. Let's do that.
 
 .. parsed-literal::
 
-    CPU times: user 636 ms, sys: 1.35 ms, total: 638 ms
-    Wall time: 638 ms
+    CPU times: user 423 ms, sys: 3.35 ms, total: 426 ms
+    Wall time: 429 ms
 
 
 
@@ -591,8 +590,8 @@ anything. Let's do that.
     Method: Rejection
     Number of samples: 1000
     Number of simulations: 1000000
-    Threshold: 0.0452
-    Sample means: t1: 0.56, t2: 0.228
+    Threshold: 0.0457
+    Sample means: t1: 0.55, t2: 0.216
 
 
 
@@ -610,8 +609,8 @@ simulations and only have to simulate the new ones:
 
 .. parsed-literal::
 
-    CPU times: user 1.72 s, sys: 10.6 ms, total: 1.73 s
-    Wall time: 1.73 s
+    CPU times: user 1.44 s, sys: 17.9 ms, total: 1.46 s
+    Wall time: 1.47 s
 
 
 
@@ -621,8 +620,8 @@ simulations and only have to simulate the new ones:
     Method: Rejection
     Number of samples: 1000
     Number of simulations: 1200000
-    Threshold: 0.0417
-    Sample means: t1: 0.561, t2: 0.225
+    Threshold: 0.0415
+    Sample means: t1: 0.55, t2: 0.215
 
 
 
@@ -640,8 +639,8 @@ standard numpy .npy files:
 
 .. parsed-literal::
 
-    CPU times: user 25.8 ms, sys: 3.27 ms, total: 29 ms
-    Wall time: 28.5 ms
+    CPU times: user 28.7 ms, sys: 4.5 ms, total: 33.2 ms
+    Wall time: 33.4 ms
 
 
 This stores the simulated data in binary ``npy`` format under
@@ -658,7 +657,7 @@ This stores the simulated data in binary ``npy`` format under
 
 .. parsed-literal::
 
-    Files in pools/arraypool_3521077242 are ['d.npy', 't1.npy', 't2.npy', 'Y.npy']
+    Files in pools/arraypool_3375867934 are ['d.npy', 't1.npy', 't2.npy', 'Y.npy']
 
 
 Now lets load all the parameters ``t1`` that were generated with numpy:
@@ -672,7 +671,7 @@ Now lets load all the parameters ``t1`` that were generated with numpy:
 
 .. parsed-literal::
 
-    array([ 0.79, -0.01, -1.47, ...,  0.98,  0.18,  0.5 ])
+    array([ 0.36,  0.47, -1.66, ...,  0.09,  0.45,  0.2 ])
 
 
 
@@ -687,7 +686,7 @@ We can also close (or save) the whole pool if we wish to continue later:
 
 .. parsed-literal::
 
-    arraypool_3521077242
+    arraypool_3375867934
 
 
 And open it up later to continue where we were left. We can open it
@@ -718,12 +717,12 @@ You can delete the files with:
         os.listdir(arraypool.path)
         
     except FileNotFoundError:
-        print("The directry is removed")
+        print("The directory is removed")
 
 
 .. parsed-literal::
 
-    The directry is removed
+    The directory is removed
 
 
 Visualizing the results
@@ -820,8 +819,9 @@ sampler:
     smc = elfi.SMC(d, batch_size=10000, seed=seed)
 
 For sampling, one has to define the number of output samples, the number
-of populations and a *schedule* i.e. a list of quantiles to use for each
-population. In essence, a population is just refined rejection sampling.
+of populations and a *schedule* i.e. a list of thresholds to use for
+each population. In essence, a population is just refined rejection
+sampling.
 
 .. code:: ipython3
 
@@ -839,8 +839,8 @@ population. In essence, a population is just refined rejection sampling.
 
 .. parsed-literal::
 
-    CPU times: user 1.72 s, sys: 154 ms, total: 1.87 s
-    Wall time: 1.56 s
+    CPU times: user 1.6 s, sys: 156 ms, total: 1.75 s
+    Wall time: 1.38 s
 
 
 We can have summaries and plots of the results just like above:

From c2b15b5708a0da6e47132c27190dec08f7e0f6d2 Mon Sep 17 00:00:00 2001
From: Henri Vuollekoski <henri.vuollekoski@gmail.com>
Date: Thu, 5 Apr 2018 17:08:11 +0300
Subject: [PATCH 3/6] Fix setting threshold=0 (#261)

---
 elfi/methods/parameter_inference.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/elfi/methods/parameter_inference.py b/elfi/methods/parameter_inference.py
index 8c70a233..61e5b397 100644
--- a/elfi/methods/parameter_inference.py
+++ b/elfi/methods/parameter_inference.py
@@ -553,7 +553,7 @@ def _update_state_meta(self):
 
     def _update_objective_n_batches(self):
         # Only in the case that the threshold is used
-        if not self.objective.get('threshold'):
+        if self.objective.get('threshold') is None:
             return
 
         s = self.state

From 6b46971575dab9814f14fa641b895d05db6b4f84 Mon Sep 17 00:00:00 2001
From: Henri Vuollekoski <henri.vuollekoski@gmail.com>
Date: Fri, 6 Apr 2018 16:11:48 +0300
Subject: [PATCH 4/6] Simpler default: batch_size=1 (#262)

---
 docs/faq.rst                             | 13 +++++++------
 elfi/methods/parameter_inference.py      |  5 ++++-
 tests/functional/test_post_processing.py | 10 ++++++----
 3 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/docs/faq.rst b/docs/faq.rst
index ed3a9a5c..8b842c47 100644
--- a/docs/faq.rst
+++ b/docs/faq.rst
@@ -17,8 +17,8 @@ the first argument defines the starting point of the interval and the second its
 **A**: Looping is relatively inefficient in Python, and so whenever possible, you should *vectorize*
 your operations_. This means that repetitive computations are performed on a batch of data using
 precompiled libraries (typically NumPy_), which effectively runs the loops in faster, compiled C-code.
-Due to the potentially huge saving in CPU-time, operations including user-code are by default assumed to
-be vectorized in ELFI. This must be accounted for.
+ELFI supports vectorized operations, and due to the potentially huge saving in CPU-time it is
+recommended to vectorize all user-code whenever possible.
 
 .. _operations: good-to-know.html#operations
 .. _NumPy: http://www.numpy.org/
@@ -26,14 +26,15 @@ be vectorized in ELFI. This must be accounted for.
 For example, imagine you have a simulator that depends on a scalar parameter and produces a vector of 5
 values. When this is used in ELFI with ``batch_size`` set to 1000, ELFI draws 1000 values from the
 parameter's prior distribution and gives this *vector* to the simulator. Ideally, the simulator should
-efficiently process all 1000 parameter cases in one go and output a numpy array of shape (1000, 5). In
-ELFI, the length (i.e. the first dimension) of all output arrays should equal ``batch_size`` **even if
-it is 1**. Note that because of this the evaluation of summary statistics, distances etc. should
+efficiently process all 1000 parameter cases in one go and output an array of shape (1000, 5). When using
+vectorized operations in ELFI, the length (i.e. the first dimension) of all output arrays should equal
+``batch_size``. Note that because of this the evaluation of summary statistics, distances etc. should
 bypass the first dimension (e.g. with NumPy functions using ``axis=1`` in this case).
 
 See ``elfi.examples`` for tips on how to vectorize simulators and work with ELFI. In case you are
 unable to vectorize your simulator, you can use `elfi.tools.vectorize`_ to mimic
-vectorized behaviour, though without the performance benefits.
+vectorized behaviour, though without the performance benefits. Finally, for simplicity vectorization
+is not assumed (``batch_size=1`` by default).
 
 .. _`elfi.tools.vectorize`: api.html#elfi.tools.vectorize
 
diff --git a/elfi/methods/parameter_inference.py b/elfi/methods/parameter_inference.py
index 61e5b397..07028949 100644
--- a/elfi/methods/parameter_inference.py
+++ b/elfi/methods/parameter_inference.py
@@ -61,7 +61,7 @@ class ParameterInference:
     def __init__(self,
                  model,
                  output_names,
-                 batch_size=1000,
+                 batch_size=1,
                  seed=None,
                  pool=None,
                  max_parallel_batches=None):
@@ -76,6 +76,9 @@ def __init__(self,
         output_names : list
             Names of the nodes whose outputs will be requested from the ELFI graph.
         batch_size : int, optional
+            The number of parameter evaluations in each pass through the ELFI graph.
+            When using a vectorized simulator, using a suitably large batch_size can provide
+            a significant performance boost.
         seed : int, optional
             Seed for the data generation from the ElfiModel
         pool : OutputPool, optional
diff --git a/tests/functional/test_post_processing.py b/tests/functional/test_post_processing.py
index a7e2e0da..0f4ad12f 100644
--- a/tests/functional/test_post_processing.py
+++ b/tests/functional/test_post_processing.py
@@ -22,7 +22,7 @@ def test_single_parameter_linear_adjustment():
     """A regression test against values obtained in the notebook."""
     seed = 20170616
     n_obs = 50
-    batch_size = 100
+    batch_size = 1000
     mu, sigma = (5, 1)
 
     # Hyperparameters
@@ -44,7 +44,8 @@ def test_single_parameter_linear_adjustment():
     elfi.Summary(lambda x: x.mean(axis=1), m['gauss'], name='ss_mean')
     elfi.Distance('euclidean', m['ss_mean'], name='d')
 
-    res = elfi.Rejection(m['d'], output_names=['ss_mean'], seed=seed).sample(1000, threshold=1)
+    res = elfi.Rejection(m['d'], output_names=['ss_mean'], batch_size=batch_size,
+                         seed=seed).sample(1000, threshold=1)
     adj = elfi.adjust_posterior(model=m, sample=res, parameter_names=['mu'], summary_names=['ss_mean'])
 
     assert np.allclose(_statistics(adj.outputs['mu']), (4.9772879640569778, 0.02058680115402544))
@@ -55,7 +56,7 @@ def test_nonfinite_values():
     """A regression test against values obtained in the notebook."""
     seed = 20170616
     n_obs = 50
-    batch_size = 100
+    batch_size = 1000
     mu, sigma = (5, 1)
 
     # Hyperparameters
@@ -77,7 +78,8 @@ def test_nonfinite_values():
     elfi.Summary(lambda x: x.mean(axis=1), m['gauss'], name='ss_mean')
     elfi.Distance('euclidean', m['ss_mean'], name='d')
 
-    res = elfi.Rejection(m['d'], output_names=['ss_mean'], seed=seed).sample(1000, threshold=1)
+    res = elfi.Rejection(m['d'], output_names=['ss_mean'], batch_size=batch_size,
+                         seed=seed).sample(1000, threshold=1)
 
     # Add some invalid values
     res.outputs['mu'] = np.append(res.outputs['mu'], np.array([np.inf]))

From 3bf2b53286c9b6bd88dce54316925ac788b38b94 Mon Sep 17 00:00:00 2001
From: Henri Vuollekoski <henri.vuollekoski@gmail.com>
Date: Tue, 10 Apr 2018 21:13:50 +0300
Subject: [PATCH 5/6] Implemented model comparison. (#247)

* Implemented model comparison.

* Improve documentation

* Address comment

* Relax the requirement for equal n_samples
---
 CHANGELOG.rst                      |  8 +++-
 README.md                          |  2 +
 docs/api.rst                       |  7 ++++
 docs/index.rst                     |  2 +
 elfi/__init__.py                   |  3 +-
 elfi/methods/model_selection.py    | 59 ++++++++++++++++++++++++++++++
 tests/unit/test_model_selection.py | 19 ++++++++++
 7 files changed, 98 insertions(+), 2 deletions(-)
 create mode 100644 elfi/methods/model_selection.py
 create mode 100644 tests/unit/test_model_selection.py

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 71ca6094..fa65a4eb 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,9 +1,15 @@
 Changelog
 =========
 
+dev
+---
+- Implemented model selection (elfi.compare_models). See API documentation.
+
 0.7 (2017-11-30)
 ----------------
-
+- Added new example: the stochastic Lotka-Volterra model
+- Fix methods.bo.utils.minimize to be strictly within bounds
+- Implemented the Two Stage Procedure, a method of summary-statistics diagnostics
 - Added the MaxVar acquisition method
 - Added the RandMaxVar acquisition method
 - Added the ExpIntVar acquisition method
diff --git a/README.md b/README.md
index fc6c944f..dc0b8d81 100644
--- a/README.md
+++ b/README.md
@@ -28,6 +28,8 @@ Other notable included algorithms and methods:
 - Bayesian Optimization
 - [No-U-Turn-Sampler](http://jmlr.org/papers/volume15/hoffman14a/hoffman14a.pdf), a Hamiltonian Monte Carlo MCMC sampler
 
+ELFI also integrates tools for visualization, model comparison, diagnostics and post-processing.
+
 See examples under [notebooks](https://github.com/elfi-dev/notebooks) to get started. Full
 documentation can be found at http://elfi.readthedocs.io/. Limited user-support may be
 asked from elfi-support.at.hiit.fi, but the 
diff --git a/docs/api.rst b/docs/api.rst
index 913f4139..03791f1b 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -265,6 +265,13 @@ Inference API classes
    :members:
    :inherited-members:
 
+**Model selection**
+
+.. currentmodule:: .
+
+.. autofunction:: elfi.compare_models
+
+
 Other
 .....
 
diff --git a/docs/index.rst b/docs/index.rst
index c6a78876..d5102660 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -34,6 +34,8 @@ ELFI also has the following non LFI methods:
 
 .. _No-U-Turn-Sampler: http://jmlr.org/papers/volume15/hoffman14a/hoffman14a.pdf
 
+Additionally, ELFI integrates tools for visualization, model comparison, diagnostics and post-processing.
+
 
 .. toctree::
     :maxdepth: 1
diff --git a/elfi/__init__.py b/elfi/__init__.py
index 5187b668..508bc038 100644
--- a/elfi/__init__.py
+++ b/elfi/__init__.py
@@ -12,6 +12,7 @@
 import elfi.model.tools as tools
 from elfi.client import get_client, set_client
 from elfi.methods.diagnostics import TwoStageSelection
+from elfi.methods.model_selection import *
 from elfi.methods.parameter_inference import *
 from elfi.methods.post_processing import adjust_posterior
 from elfi.model.elfi_model import *
@@ -24,4 +25,4 @@
 __email__ = 'elfi-support@hiit.fi'
 
 # make sure __version_ is on the last non-empty line (read by setup.py)
-__version__ = '0.7'
+__version__ = '0.7_dev'
diff --git a/elfi/methods/model_selection.py b/elfi/methods/model_selection.py
new file mode 100644
index 00000000..794f6e4f
--- /dev/null
+++ b/elfi/methods/model_selection.py
@@ -0,0 +1,59 @@
+"""This module contains methods for model comparison and selection."""
+
+import numpy as np
+
+
+def compare_models(sample_objs, model_priors=None):
+    """Find posterior probabilities for different models.
+
+    The algorithm requires elfi.Sample objects from prerun inference methods. For example the
+    output from elfi.Rejection.sample is valid. The portion of samples for each model in the top
+    discrepancies are adjusted by each models acceptance ratio and prior probability.
+
+    The discrepancies (including summary statistics) must be comparable so that it is
+    meaningful to sort them!
+
+    Parameters
+    ----------
+    sample_objs : list of elfi.Sample
+        Resulting Sample objects from prerun inference models. The objects must include
+        a valid `discrepancies` attribute.
+    model_priors : array_like, optional
+        Prior probability of each model. Defaults to 1 / n_models.
+
+    Returns
+    -------
+    np.array
+        Posterior probabilities for the considered models.
+
+    """
+    n_models = len(sample_objs)
+    n_min = min([s.n_samples for s in sample_objs])
+
+    # concatenate discrepancy vectors
+    try:
+        discrepancies = np.concatenate([s.discrepancies for s in sample_objs])
+    except ValueError:
+        raise ValueError("All Sample objects must include valid discrepancies.")
+
+    # sort and take the smallest n_min
+    inds = np.argsort(discrepancies)[:n_min]
+
+    # calculate the portions of accepted samples for each model in the top discrepancies
+    p_models = np.empty(n_models)
+    up_bound = 0
+    for i in range(n_models):
+        low_bound = up_bound
+        up_bound += sample_objs[i].n_samples
+        p_models[i] = np.logical_and(inds >= low_bound, inds < up_bound).sum()
+
+        # adjust by the number of simulations run
+        p_models[i] /= sample_objs[i].n_sim
+
+        # adjust by the prior model probability
+        if model_priors is not None:
+            p_models[i] *= model_priors[i]
+
+    p_models = p_models / p_models.sum()
+
+    return p_models
diff --git a/tests/unit/test_model_selection.py b/tests/unit/test_model_selection.py
new file mode 100644
index 00000000..39f6b368
--- /dev/null
+++ b/tests/unit/test_model_selection.py
@@ -0,0 +1,19 @@
+import elfi
+from elfi.examples import gauss, ma2
+
+
+def test_compare_models():
+    m = gauss.get_model()
+    res1 = elfi.Rejection(m['d']).sample(100)
+
+    # use less informative prior
+    m['mu'].become(elfi.Prior('uniform', -10, 50))
+    res2 = elfi.Rejection(m['d']).sample(100)
+
+    # use different simulator
+    m['gauss'].become(elfi.Simulator(ma2.MA2, m['mu'], m['sigma'], observed=m.observed['gauss']))
+    res3 = elfi.Rejection(m['d']).sample(100)
+
+    p = elfi.compare_models([res1, res2, res3])
+    assert p[0] > p[1]
+    assert p[1] > p[2]

From 79c7e9c693bc5fd248de8e5e816225e5ac2197c1 Mon Sep 17 00:00:00 2001
From: Henri Vuollekoski <henri.vuollekoski@gmail.com>
Date: Thu, 12 Apr 2018 09:32:36 +0300
Subject: [PATCH 6/6] Preparing version 0.7.1 to dev (#264)

* Release 0.7.1

* Fix flake8 errors
---
 CHANGELOG.rst       | 6 ++++--
 README.md           | 2 +-
 elfi/__init__.py    | 2 +-
 elfi/model/tools.py | 6 +++---
 elfi/store.py       | 6 +++---
 5 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index fa65a4eb..fa5096ee 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,9 +1,11 @@
 Changelog
 =========
 
-dev
----
+0.7.1 (2018-04-11)
+------------------
 - Implemented model selection (elfi.compare_models). See API documentation.
+- Fix threshold=0 in rejection sampling
+- Set default batch_size to 1 in ParameterInference base class
 
 0.7 (2017-11-30)
 ----------------
diff --git a/README.md b/README.md
index dc0b8d81..2be28d32 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-**Version 0.7 released!** See the CHANGELOG and [notebooks](https://github.com/elfi-dev/notebooks).
+**Version 0.7.1 released!** See the CHANGELOG and [notebooks](https://github.com/elfi-dev/notebooks).
 
 **NOTE:** For the time being NetworkX 2 is incompatible with ELFI.
 
diff --git a/elfi/__init__.py b/elfi/__init__.py
index 508bc038..fa178453 100644
--- a/elfi/__init__.py
+++ b/elfi/__init__.py
@@ -25,4 +25,4 @@
 __email__ = 'elfi-support@hiit.fi'
 
 # make sure __version_ is on the last non-empty line (read by setup.py)
-__version__ = '0.7_dev'
+__version__ = '0.7.1'
diff --git a/elfi/model/tools.py b/elfi/model/tools.py
index c0aab4b7..77f3d65c 100644
--- a/elfi/model/tools.py
+++ b/elfi/model/tools.py
@@ -44,10 +44,10 @@ def run_vectorized(operation, *inputs, constants=None, dtype=None, batch_size=No
 
         # Test if a numpy array
         if is_array(inpt):
-            l = len(inpt)
+            length = len(inpt)
             if batch_size is None:
-                batch_size = l
-            elif batch_size != l:
+                batch_size = length
+            elif batch_size != length:
                 raise ValueError('Batch size {} does not match with input {} length of '
                                  '{}. Please check `constants` argument for the '
                                  'vectorize decorator for marking constant inputs.')
diff --git a/elfi/store.py b/elfi/store.py
index 02a3749d..c5aead02 100644
--- a/elfi/store.py
+++ b/elfi/store.py
@@ -208,12 +208,12 @@ def _make_store_for(self, node):
 
     def __len__(self):
         """Return the largest batch index in any of the stores."""
-        l = 0
+        largest = 0
         for output, store in self.stores.items():
             if store is None:
                 continue
-            l = max(l, len(store))
-        return l
+            largest = max(largest, len(store))
+        return largest
 
     def __getitem__(self, batch_index):
         """Return the batch."""