diff --git a/.dockerignore b/.dockerignore new file mode 100644 index 00000000..d09be028 --- /dev/null +++ b/.dockerignore @@ -0,0 +1,106 @@ +# IDE settings +.idea + +# Byte-compiled / optimized / DLL files +**/__pycache__ +*.py[cod] +*$py.class + +# C/C++ extensions +*.so +*.bin + +# Distribution / packaging +.Python +env +build +develop-eggs +dist +downloads +eggs +.eggs +lib +lib64 +parts +sdist +var +*.egg-info +.installed.cfg +*.egg + +# Images +*.png +*.svg +*.jpg +*.jpeg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*cover +.hypothesis + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py + +# Flask stuff: +instance +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build + +# PyBuilder +target + +# IPython Notebook +.ipynb_checkpoints + +# pyenv +.python-version + +# celery beat schedule file +celerybeat-schedule + +# dotenv +.env + +# virtualenv +venv +ENV +.venv + +# Spyder project settings +.spyderproject + +# Rope project settings +.ropeproject + +# Editor temp files +*~ +*.swo +*.swp + diff --git a/CHANGELOG.rst b/CHANGELOG.rst index fa5096ee..dd2511b2 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,12 @@ Changelog ========= +0.7.2 (2018-06-20) +------------------ +- Added support for kwargs in elfi.set_client +- Added new example: inference of transmission dynamics of bacteria in daycare centers +- Added new example: Lorenz model + 0.7.1 (2018-04-11) ------------------ - Implemented model selection (elfi.compare_models). See API documentation. diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 00000000..7bedd679 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,16 @@ +FROM python:3.6 + +RUN apt-get update && apt-get install -y graphviz + +WORKDIR /elfi +ADD . /elfi + +RUN pip install numpy graphviz +RUN pip install -e . +RUN pip install -r requirements-dev.txt + +# matplotlib unable to show figures +RUN mkdir -p /root/.config/matplotlib +RUN echo "backend : Agg" >> /root/.config/matplotlib/matplotlibrc + +CMD ["/bin/bash"] diff --git a/Makefile b/Makefile index 8d615921..5b18667c 100644 --- a/Makefile +++ b/Makefile @@ -52,7 +52,7 @@ lint: ## check style with flake8 flake8 elfi tests test: ## run tests quickly with the default Python - PYTHONPATH=$$PYTHONPATH:. py.test + PYTHONPATH=$$PYTHONPATH:. py.test --reruns 1 test-all: ## run tests on every Python version with tox tox diff --git a/README.md b/README.md index 2be28d32..d1a7764f 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -**Version 0.7.1 released!** See the CHANGELOG and [notebooks](https://github.com/elfi-dev/notebooks). +**Version 0.7.2 released!** See the CHANGELOG and [notebooks](https://github.com/elfi-dev/notebooks). **NOTE:** For the time being NetworkX 2 is incompatible with ELFI. @@ -75,6 +75,17 @@ source activate elfi pip install elfi ``` +### Docker container + +A simple Dockerfile for command-line interface is also provided. Please see [Docker documentation](https://docs.docker.com/). + +``` +git clone --depth 1 https://github.com/elfi-dev/elfi.git +cd elfi +docker build -t elfi . +docker run -it elfi +``` + ### Potential problems with installation ELFI depends on several other Python packages, which have their own dependencies. @@ -94,9 +105,9 @@ If you wish to cite ELFI, please use the paper in [arXiv](https://arxiv.org/abs/ ``` @misc{1708.00707, -Author = {Jarno Lintusaari and Henri Vuollekoski and Antti Kangasrääsiö and Kusti Skytén and Marko Järvenpää and Michael Gutmann and Aki Vehtari and Jukka Corander and Samuel Kaski}, +Author = {Jarno Lintusaari and Henri Vuollekoski and Antti Kangasrääsiö and Kusti Skytén and Marko Järvenpää and Pekka Marttinen and Michael Gutmann and Aki Vehtari and Jukka Corander and Samuel Kaski}, Title = {ELFI: Engine for Likelihood Free Inference}, -Year = {2017}, +Year = {2018}, Eprint = {arXiv:1708.00707}, } ``` diff --git a/docs/index.rst b/docs/index.rst index d5102660..d88f94a0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -82,8 +82,8 @@ If you wish to cite ELFI, please use the paper in arXiv_: .. code-block:: console @misc{1708.00707, - Author = {Jarno Lintusaari and Henri Vuollekoski and Antti Kangasrääsiö and Kusti Skytén and Marko Järvenpää and Michael Gutmann and Aki Vehtari and Jukka Corander and Samuel Kaski}, + Author = {Jarno Lintusaari and Henri Vuollekoski and Antti Kangasrääsiö and Kusti Skytén and Marko Järvenpää and Pekka Marttinen and Michael Gutmann and Aki Vehtari and Jukka Corander and Samuel Kaski}, Title = {ELFI: Engine for Likelihood Free Inference}, - Year = {2017}, + Year = {2018}, Eprint = {arXiv:1708.00707}, } diff --git a/docs/installation.rst b/docs/installation.rst index bd5e6b73..65d3252a 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -89,3 +89,16 @@ means the current folder. .. _Github repo: https://github.com/elfi-dev/elfi .. _tarball: https://github.com/elfi-dev/elfi/tarball/dev +Docker container +---------------- + +A simple Dockerfile for command-line interface is also provided. Please see `Docker documentation`_. + +.. _Docker documentation: https://docs.docker.com/ + +.. code-block:: console + + git clone --depth 1 https://github.com/elfi-dev/elfi.git + cd elfi + docker build -t elfi . + docker run -it elfi diff --git a/docs/usage/parallelization.rst b/docs/usage/parallelization.rst index 763eb125..42e4409a 100644 --- a/docs/usage/parallelization.rst +++ b/docs/usage/parallelization.rst @@ -56,11 +56,11 @@ in your computer. You can activate it simply by elfi.set_client('multiprocessing') -Any inference instance created after you have set the new client will -automatically use it to perform the computations. Let's try it with our -MA2 example model from the tutorial. When running the next command, take -a look at the system monitor of your operating system; it should show -that all of your cores are doing heavy computation simultaneously. +Any inference instance created **after** you have set the new client +will automatically use it to perform the computations. Let's try it with +our MA2 example model from the tutorial. When running the next command, +take a look at the system monitor of your operating system; it should +show that all of your cores are doing heavy computation simultaneously. .. code:: ipython3 @@ -70,8 +70,8 @@ that all of your cores are doing heavy computation simultaneously. .. parsed-literal:: - CPU times: user 272 ms, sys: 28 ms, total: 300 ms - Wall time: 2.41 s + CPU times: user 298 ms, sys: 25.7 ms, total: 324 ms + Wall time: 3.93 s And that is it. The result object is also just like in the basic case: @@ -91,14 +91,63 @@ And that is it. The result object is also just like in the basic case: Method: Rejection Number of samples: 5000 Number of simulations: 1000000 - Threshold: 0.0817 - Sample means: t1: 0.68, t2: 0.133 + Threshold: 0.0826 + Sample means: t1: 0.694, t2: 0.226 .. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/parallelization_files/parallelization_11_1.png +Note that for reproducibility a reference to the activated client is +saved in the inference instance: + +.. code:: ipython3 + + rej.client + + + + +.. parsed-literal:: + + + + + +If you want to change the client for an existing inference instance, you +have to do something like this: + +.. code:: ipython3 + + elfi.set_client('native') + rej.client = elfi.get_client() + rej.client + + + + +.. parsed-literal:: + + + + + +By default the multiprocessing client will use all cores on your system. +This is not always desirable, as the operating system may prioritize +some other process, leaving ELFI queuing for the promised resources. You +can define some other number of processes like so: + +.. code:: ipython3 + + elfi.set_client(elfi.clients.multiprocessing.Client(num_processes=3)) + +**Note:** The ``multiprocessing`` library may require additional care +under Windows. If you receive a RuntimeError mentioning +``freeze_support``, please include a call to +``multiprocessing.freeze_support()``, see +`documentation `__. + Ipyparallel client ------------------ @@ -136,8 +185,8 @@ take care of the parallelization from now on: .. parsed-literal:: - CPU times: user 3.16 s, sys: 184 ms, total: 3.35 s - Wall time: 13.4 s + CPU times: user 3.47 s, sys: 288 ms, total: 3.76 s + Wall time: 18.1 s To summarize, the only thing that needed to be changed from the basic @@ -230,8 +279,8 @@ The above may look a bit cumbersome, but now this works: Method: Rejection Number of samples: 1000 Number of simulations: 100000 - Threshold: 0.0136 - Sample means: t1: 0.676, t2: 0.129 + Threshold: 0.0146 + Sample means: t1: 0.693, t2: 0.233 @@ -250,5 +299,5 @@ Remember to stop the ipcluster when done .. parsed-literal:: - 2017-07-19 16:20:58.662 [IPClusterStop] Stopping cluster [pid=21020] with [signal=] + 2018-04-24 19:14:56.997 [IPClusterStop] Stopping cluster [pid=39639] with [signal=] diff --git a/docs/usage/tutorial.rst b/docs/usage/tutorial.rst index f9ccadf9..4116c946 100644 --- a/docs/usage/tutorial.rst +++ b/docs/usage/tutorial.rst @@ -26,7 +26,7 @@ settings. %precision 2 # Set an arbitrary seed and a global random state to keep the randomly generated quantities the same between runs - seed = 20170530 + seed = 20170530 # this will be separately given to ELFI np.random.seed(seed) Inference with ELFI: case MA(2) model @@ -79,15 +79,16 @@ Vectorization What is the purpose of the ``batch_size`` argument? In ELFI, operations are vectorized, meaning that instead of simulating a single MA2 sequence -at a time, we simulate a batch of them. A vectorized function takes +at a time, we simulate a *batch* of them. A vectorized function takes vectors as inputs, and computes the output for each element in the vector. Vectorization is a way to make operations efficient in Python. -Above we rely on numpy to carry out the vectorized calculations. -In this case the arguments ``t1`` and ``t2`` are going to be vectors of -length ``batch_size`` and the method returns a 2d array with the -simulations on the rows. Notice that for convenience, the funtion also -works with scalars that are first converted to vectors. +In the MA2 simulator above we rely on numpy to carry out the vectorized +calculations. The arguments ``t1`` and ``t2`` are going to be +**vectors** of length ``batch_size`` and the function returns a 2D array +of shape ``(batch_size, n_obs)`` with each row corresponding to a single +argument pair. Notice that for convenience, the funtion also works with +scalars as they are first converted to vectors. .. note:: There is a built-in tool (`elfi.tools.vectorize`) in ELFI to vectorize operations that are not vectorized. It is basically a for loop wrapper. @@ -117,7 +118,7 @@ observed data alone. -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_10_0.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_11_0.png Approximate Bayesian Computation @@ -205,7 +206,9 @@ between those. Here, we will apply the intuition arising from the definition of the MA(2) process, and use the autocovariances with lags 1 and 2 as the -summary statistics: +summary statistics. Note that since the rows of ``x`` correspond to +independent simulations, we have to tell this numpy function to take +row-wise means by the keyword argument ``axis=1``: .. code:: ipython3 @@ -246,7 +249,7 @@ a DAG. -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_27_0.svg +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_28_0.svg @@ -267,6 +270,8 @@ which :math:`-2<\theta_1<2` with :math:`\theta_1+\theta_2>-1` and :math:`\theta_1-\theta_2<1` i.e. the parameters are sampled from a triangle (see below). +.. note:: By default all created nodes (even independent ones) will belong to the same ELFI model. It's good practice to always check with `elfi.draw` that the result looks as intended. A new default model can be started and set with the call `elfi.new_model()`. One can also create a new model with `my_model = elfi.ElfiModel()` and pass this as a keyword argument `model=my_model` when creating new nodes. Several ELFI models can coexist. + Custom priors ~~~~~~~~~~~~~ @@ -309,7 +314,7 @@ These indeed sample from a triangle: -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_33_0.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_36_0.png Let's change the earlier priors to the new ones in the inference model: @@ -324,7 +329,7 @@ Let's change the earlier priors to the new ones in the inference model: -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_35_0.svg +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_38_0.svg @@ -384,7 +389,7 @@ time is spent in drawing. -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_42_0.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_45_0.png @@ -395,8 +400,8 @@ time is spent in drawing. .. parsed-literal:: - CPU times: user 1.6 s, sys: 166 ms, total: 1.77 s - Wall time: 1.76 s + CPU times: user 1.89 s, sys: 173 ms, total: 2.06 s + Wall time: 2.13 s The ``sample`` method returns a ``Sample`` object, which contains @@ -416,7 +421,7 @@ quantile. .. parsed-literal:: - 0.56 + 0.55600915483879665 @@ -451,8 +456,8 @@ as long as it takes to generate the requested number of samples. .. parsed-literal:: - CPU times: user 198 ms, sys: 35.5 ms, total: 233 ms - Wall time: 231 ms + CPU times: user 221 ms, sys: 40 ms, total: 261 ms + Wall time: 278 ms Method: Rejection Number of samples: 1000 Number of simulations: 40000 @@ -496,9 +501,9 @@ been reached or a maximum of one second of time has been used. Method: Rejection Number of samples: 1000 - Number of simulations: 180000 - Threshold: 0.088 - Sample means: t1: 0.561, t2: 0.221 + Number of simulations: 130000 + Threshold: 0.104 + Sample means: t1: 0.558, t2: 0.219 @@ -538,7 +543,7 @@ in our model: .. code:: ipython3 pool = elfi.OutputPool(['t1', 't2', 'S1', 'S2']) - rej = elfi.Rejection(d, pool=pool) + rej = elfi.Rejection(d, batch_size=10000, pool=pool) %time result3 = rej.sample(N, n_sim=1000000) result3 @@ -546,8 +551,8 @@ in our model: .. parsed-literal:: - CPU times: user 5.01 s, sys: 60.9 ms, total: 5.07 s - Wall time: 5.09 s + CPU times: user 6.13 s, sys: 1.15 s, total: 7.29 s + Wall time: 8.57 s @@ -557,8 +562,8 @@ in our model: Method: Rejection Number of samples: 1000 Number of simulations: 1000000 - Threshold: 0.0363 - Sample means: t1: 0.554, t2: 0.216 + Threshold: 0.0364 + Sample means: t1: 0.557, t2: 0.211 @@ -571,7 +576,7 @@ anything. Let's do that. # Replace the current distance with a cityblock (manhattan) distance and recreate the inference d.become(elfi.Distance('cityblock', S1, S2, p=1)) - rej = elfi.Rejection(d, pool=pool) + rej = elfi.Rejection(d, batch_size=10000, pool=pool) %time result4 = rej.sample(N, n_sim=1000000) result4 @@ -579,8 +584,8 @@ anything. Let's do that. .. parsed-literal:: - CPU times: user 423 ms, sys: 3.35 ms, total: 426 ms - Wall time: 429 ms + CPU times: user 161 ms, sys: 2.84 ms, total: 163 ms + Wall time: 167 ms @@ -590,8 +595,8 @@ anything. Let's do that. Method: Rejection Number of samples: 1000 Number of simulations: 1000000 - Threshold: 0.0457 - Sample means: t1: 0.55, t2: 0.216 + Threshold: 0.0456 + Sample means: t1: 0.56, t2: 0.214 @@ -609,8 +614,8 @@ simulations and only have to simulate the new ones: .. parsed-literal:: - CPU times: user 1.44 s, sys: 17.9 ms, total: 1.46 s - Wall time: 1.47 s + CPU times: user 1.14 s, sys: 185 ms, total: 1.33 s + Wall time: 1.34 s @@ -621,7 +626,7 @@ simulations and only have to simulate the new ones: Number of samples: 1000 Number of simulations: 1200000 Threshold: 0.0415 - Sample means: t1: 0.55, t2: 0.215 + Sample means: t1: 0.56, t2: 0.216 @@ -633,14 +638,14 @@ standard numpy .npy files: .. code:: ipython3 arraypool = elfi.ArrayPool(['t1', 't2', 'Y', 'd']) - rej = elfi.Rejection(d, pool=arraypool) + rej = elfi.Rejection(d, batch_size=10000, pool=arraypool) %time result5 = rej.sample(100, threshold=0.3) .. parsed-literal:: - CPU times: user 28.7 ms, sys: 4.5 ms, total: 33.2 ms - Wall time: 33.4 ms + CPU times: user 66.6 ms, sys: 70.3 ms, total: 137 ms + Wall time: 175 ms This stores the simulated data in binary ``npy`` format under @@ -657,7 +662,7 @@ This stores the simulated data in binary ``npy`` format under .. parsed-literal:: - Files in pools/arraypool_3375867934 are ['d.npy', 't1.npy', 't2.npy', 'Y.npy'] + Files in pools/arraypool_4290044000 are ['d.npy', 't1.npy', 't2.npy', 'Y.npy'] Now lets load all the parameters ``t1`` that were generated with numpy: @@ -671,7 +676,7 @@ Now lets load all the parameters ``t1`` that were generated with numpy: .. parsed-literal:: - array([ 0.36, 0.47, -1.66, ..., 0.09, 0.45, 0.2 ]) + array([ 0.42, -1.15, 1.3 , ..., 0.64, 1.06, -0.47]) @@ -686,7 +691,7 @@ We can also close (or save) the whole pool if we wish to continue later: .. parsed-literal:: - arraypool_3375867934 + arraypool_4290044000 And open it up later to continue where we were left. We can open it @@ -703,7 +708,7 @@ using its name: .. parsed-literal:: - This pool has 3 batches + This pool has 1 batches You can delete the files with: @@ -740,7 +745,7 @@ For example one can plot the marginal distributions: -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_74_0.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_77_0.png Often "pairwise relationships" are more informative: @@ -751,7 +756,7 @@ Often "pairwise relationships" are more informative: -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_76_0.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_79_0.png Note that if working in a non-interactive environment, you can use e.g. @@ -839,8 +844,8 @@ sampling. .. parsed-literal:: - CPU times: user 1.6 s, sys: 156 ms, total: 1.75 s - Wall time: 1.38 s + CPU times: user 1.72 s, sys: 183 ms, total: 1.9 s + Wall time: 1.65 s We can have summaries and plots of the results just like above: @@ -901,15 +906,15 @@ Or just the means: -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_89_0.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_92_0.png -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_89_1.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_92_1.png -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_89_2.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_92_2.png Obviously one still has direct access to the samples as well, which @@ -932,7 +937,7 @@ allows custom plotting: -.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_91_0.png +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_94_0.png It can be seen that the populations iteratively concentrate more and @@ -941,5 +946,71 @@ SMC are weighed, and the weights should be accounted for when interpreting the results. ELFI does this automatically when computing the mean, for example. +What next? +---------- + +If you want to play with different ABC algorithms, no need to repeat the +simulator definitions etc. from this notebook. ELFI provides a +convenient way to quickly get you going: + +.. code:: ipython3 + + from elfi.examples import ma2 + ma2_model = ma2.get_model() + +This constructs the same ELFI graph discussed in this tutorial. The +example is self-contained and includes implementations of all relevant +operations. + +.. code:: ipython3 + + elfi.draw(ma2_model) + + + + +.. image:: http://research.cs.aalto.fi/pml/software/elfi/docs/0.6.2/usage/tutorial_files/tutorial_100_0.svg + + + +.. code:: ipython3 + + elfi.Rejection(ma2_model['d'], batch_size=10000).sample(1000) + + + + +.. parsed-literal:: + + Method: Rejection + Number of samples: 1000 + Number of simulations: 100000 + Threshold: 0.128 + Sample means: t1: 0.719, t2: 0.412 + + + +ELFI ships with many other common example cases from ABC literature, and +they all include the ``get_model`` method mentioned above. The source +codes of these examples are also good for learning more about best +practices with ELFI. + +.. code:: ipython3 + + !ls {elfi.examples.__path__[0] + '/*.py'} | xargs basename + + +.. parsed-literal:: + + __init__.py + bdm.py + bignk.py + gauss.py + gnk.py + lotka_volterra.py + ma2.py + ricker.py + + That's it! See the other documentation for more advanced topics on e.g. BOLFI, external simulators and parallelization. diff --git a/elfi/__init__.py b/elfi/__init__.py index fa178453..dd9da71f 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.1' +__version__ = '0.7.2' diff --git a/elfi/client.py b/elfi/client.py index 3b589cf7..071cc274 100644 --- a/elfi/client.py +++ b/elfi/client.py @@ -28,13 +28,22 @@ def get_client(): return _client -def set_client(client=None): - """Set the current ELFI client instance.""" +def set_client(client=None, **kwargs): + """Set the current ELFI client instance. + + Parameters + ---------- + client : ClientBase or str + Instance of a client from ClientBase, + or a string from ['native', 'multiprocessing', 'ipyparallel']. + If string, the respective constructor is called with `kwargs`. + + """ global _client if isinstance(client, str): m = importlib.import_module('elfi.clients.{}'.format(client)) - client = m.Client() + client = m.Client(**kwargs) _client = client diff --git a/elfi/clients/ipyparallel.py b/elfi/clients/ipyparallel.py index 5c229b45..96f26ad9 100644 --- a/elfi/clients/ipyparallel.py +++ b/elfi/clients/ipyparallel.py @@ -25,7 +25,7 @@ class Client(elfi.client.ClientBase): http://ipyparallel.readthedocs.io """ - def __init__(self, ipp_client=None): + def __init__(self, ipp_client=None, **kwargs): """Create an ipyparallel client for ELFI. Parameters @@ -34,7 +34,7 @@ def __init__(self, ipp_client=None): Use this ipyparallel client with ELFI. """ - self.ipp_client = ipp_client or ipp.Client() + self.ipp_client = ipp_client or ipp.Client(**kwargs) self.view = self.ipp_client.load_balanced_view() self.tasks = {} diff --git a/elfi/clients/multiprocessing.py b/elfi/clients/multiprocessing.py index 34152134..7bcc04a9 100644 --- a/elfi/clients/multiprocessing.py +++ b/elfi/clients/multiprocessing.py @@ -18,7 +18,7 @@ def set_as_default(): class Client(elfi.client.ClientBase): """Client based on Python's built-in multiprocessing module.""" - def __init__(self, num_processes=None): + def __init__(self, num_processes=None, **kwargs): """Create a multiprocessing client. Parameters @@ -27,7 +27,8 @@ def __init__(self, num_processes=None): Number of worker processes to use. Defaults to os.cpu_count(). """ - self.pool = multiprocessing.Pool(processes=num_processes) + num_processes = num_processes or kwargs.pop('processes', None) + self.pool = multiprocessing.Pool(processes=num_processes, **kwargs) self.tasks = {} self._id_counter = itertools.count() diff --git a/elfi/clients/native.py b/elfi/clients/native.py index f72323c8..c67fa234 100644 --- a/elfi/clients/native.py +++ b/elfi/clients/native.py @@ -20,7 +20,7 @@ class Client(elfi.client.ClientBase): Responsible for sending computational graphs to be executed in an Executor """ - def __init__(self): + def __init__(self, **kwargs): """Create a native client.""" self.tasks = {} self._ids = itertools.count() diff --git a/elfi/examples/daycare.py b/elfi/examples/daycare.py new file mode 100644 index 00000000..066ce2a1 --- /dev/null +++ b/elfi/examples/daycare.py @@ -0,0 +1,305 @@ +"""Example of inference of transmission dynamics of bacteria in day care centers. + +Treatment roughly follows: +Numminen, E., Cheng, L., Gyllenberg, M. and Corander, J.: Estimating the transmission dynamics +of Streptococcus pneumoniae from strain prevalence data, Biometrics, 69, 748-757, 2013. +""" + +import logging +from functools import partial + +import numpy as np + +import elfi + + +def daycare(t1, t2, t3, n_dcc=29, n_ind=53, n_strains=33, freq_strains_commun=None, + n_obs=36, time_end=10., batch_size=1, random_state=None): + r"""Generate cross-sectional data from a stochastic variant of the SIS-model. + + This function simulates the transmission dynamics of bacterial infections in daycare centers + (DCC) as described in Nummelin et al. [2013]. The observation model is however simplified to + an equal number of sampled individuals among the daycare centers. + + The model is defined as a continuous-time Markov process with transition probabilities: + + Pr(I_{is}(t+dt)=1 | I_{is}(t)=0) = t1 * E_s(I(t)) + t2 * P_s, if \sum_{j=1}^N_s I_{ij}(t)=0 + Pr(I_{is}(t+dt)=1 | I_{is}(t)=0) = t3 * (t1 * E_s(I(t)) + t2 * P_s), otherwise + Pr(I_{is}(t+dt)=0 | I_{is}(t)=1) = \gamma + + where: + I_{is}(t) is the status of carriage of strain s for individual i. + E_s(I(t)) is the probability of sampling the strain s + t1 is the rate of transmission from other children at the DCC (\beta in paper). + t2 is the rate of transmission from the community outside the DCC (\Lambda in paper). + t3 scales the rate of an infected child being infected with another strain (\theta in paper). + \gamma is the relative probability of healing from a strain. + + As in the paper, \gamma=1, and the other inferred parameters are relative to it. + + The system is solved using the Direct method [Gillespie, 1977]. + + References + ---------- + Numminen, E., Cheng, L., Gyllenberg, M. and Corander, J. (2013) Estimating the transmission + dynamics of Streptococcus pneumoniae from strain prevalence data, Biometrics, 69, 748-757. + Gillespie, D. T. (1977) Exact stochastic simulation of coupled chemical reactions. + The Journal of Physical Chemistry 81 (25), 2340–2361. + + Parameters + ---------- + t1 : float or np.array + Rate of transmission from other individuals at the DCC. + t2 : float or np.array + Rate of transmission from the community outside the DCC. + t3 : float or np.array + Scaling of co-infection for individuals infected with another strain. + n_dcc : int, optional + Number of daycare centers. + n_ind : int, optional + Number of individuals in a DCC (same for all). + n_strains : int, optional + Number of bacterial strains considered. + freq_strains_commun : np.array of shape (n_strains,), optional + Prevalence of each strain in the community outside the DCC. Defaults to 0.1. + n_obs : int, optional + Number of individuals sampled from each DCC (same for all). + time_end : float, optional + The system is solved using the Direct method until all cases within the batch exceed this. + batch_size : int, optional + random_state : np.random.RandomState, optional + + Returns + ------- + state_obs : np.array + Observations in shape (batch_size, n_dcc, n_obs, n_strains). + + """ + random_state = random_state or np.random + + t1 = np.asanyarray(t1).reshape((-1, 1, 1, 1)) + t2 = np.asanyarray(t2).reshape((-1, 1, 1, 1)) + t3 = np.asanyarray(t3).reshape((-1, 1, 1, 1)) + + if freq_strains_commun is None: + freq_strains_commun = np.full(n_strains, 0.1) + + prob_commun = t2 * freq_strains_commun + + # the state (infection status) is a 4D tensor for computational performance + state = np.zeros((batch_size, n_dcc, n_ind, n_strains), dtype=np.bool) + + # time for each DCC in the batch + time = np.zeros((batch_size, n_dcc)) + + n_factor = 1. / (n_ind - 1) + gamma = 1. # relative, see paper + ind_b_dcc = [np.repeat(np.arange(batch_size), n_dcc), np.tile(np.arange(n_dcc), batch_size)] + + while np.any(time < time_end): + with np.errstate(divide='ignore', invalid='ignore'): + # probability of sampling a strain; in paper: E_s(I(t)) + prob_strain = np.sum(state / np.sum(state, axis=3, keepdims=True), + axis=2, keepdims=True) * n_factor + prob_strain = np.where(np.isfinite(prob_strain), prob_strain, 0) + + # init prob to get infected, same for all + hazards = t1 * prob_strain + prob_commun # shape (batch_size, n_dcc, 1, n_strains) + + # co-infection, depends on the individual's state + hazards = np.tile(hazards, (1, 1, n_ind, 1)) + any_infection = np.any(state, axis=3, keepdims=True) + hazards = np.where(any_infection, t3 * hazards, hazards) + + # (relative) probability to be cured + hazards[state] = gamma + + # normalize to probabilities + inv_sum_hazards = 1. / np.sum(hazards, axis=(2, 3), keepdims=True) + probs = hazards * inv_sum_hazards + + # times until next transition (for each DCC in the batch) + delta_t = random_state.exponential(inv_sum_hazards[:, :, 0, 0]) + time = time + delta_t + + # choose transition + probs = probs.reshape((batch_size, n_dcc, -1)) + cumprobs = np.cumsum(probs[:, :, :-1], axis=2) + x = random_state.uniform(size=(batch_size, n_dcc, 1)) + ind_transit = np.sum(x >= cumprobs, axis=2) + + # update state, need to find the correct indices first + ind_transit = ind_b_dcc + list(np.unravel_index(ind_transit.ravel(), (n_ind, n_strains))) + state[ind_transit] = np.logical_not(state[ind_transit]) + + # observation model: simply take the first n_obs individuals + state_obs = state[:, :, :n_obs, :] + + return state_obs + + +def get_model(true_params=None, seed_obs=None, **kwargs): + """Return a complete ELFI graph ready for inference. + + Selection of true values, priors etc. follows the approach in + + Numminen, E., Cheng, L., Gyllenberg, M. and Corander, J.: Estimating the transmission dynamics + of Streptococcus pneumoniae from strain prevalence data, Biometrics, 69, 748-757, 2013. + + and + + Gutmann M U, Corander J (2016). Bayesian Optimization for Likelihood-Free Inference + of Simulator-Based Statistical Models. JMLR 17(125):1−47, 2016. + + Parameters + ---------- + true_params : list, optional + Parameters with which the observed data is generated. + seed_obs : int, optional + Seed for the observed data generation. + + Returns + ------- + m : elfi.ElfiModel + + """ + logger = logging.getLogger() + if true_params is None: + true_params = [3.6, 0.6, 0.1] + + m = elfi.ElfiModel() + y_obs = daycare(*true_params, random_state=np.random.RandomState(seed_obs), **kwargs) + sim_fn = partial(daycare, **kwargs) + priors = [] + sumstats = [] + + priors.append(elfi.Prior('uniform', 0, 11, model=m, name='t1')) + priors.append(elfi.Prior('uniform', 0, 2, model=m, name='t2')) + priors.append(elfi.Prior('uniform', 0, 1, model=m, name='t3')) + + elfi.Simulator(sim_fn, *priors, observed=y_obs, name='DCC') + + sumstats.append(elfi.Summary(ss_shannon, m['DCC'], name='Shannon')) + sumstats.append(elfi.Summary(ss_strains, m['DCC'], name='n_strains')) + sumstats.append(elfi.Summary(ss_prevalence, m['DCC'], name='prevalence')) + sumstats.append(elfi.Summary(ss_prevalence_multi, m['DCC'], name='multi')) + + elfi.Discrepancy(distance, *sumstats, name='d') + + logger.info("Generated observations with true parameters " + "t1: %.1f, t2: %.3f, t3: %.1f, ", *true_params) + + return m + + +def ss_shannon(data): + r"""Calculate the Shannon index of diversity of the distribution of observed strains. + + H = -\sum p \log(p) + + https://en.wikipedia.org/wiki/Diversity_index#Shannon_index + + Parameters + ---------- + data : np.array of shape (batch_size, n_dcc, n_obs, n_strains) + + Returns + ------- + np.array of shape (batch_size, n_dcc) + + """ + proportions = np.sum(data, axis=2) / data.shape[2] + shannon = -np.sum(proportions * np.log(proportions + 1e-9), axis=2) # axis 3 is now 2 + + return shannon + + +def ss_strains(data): + """Calculate the number of different strains observed. + + Parameters + ---------- + data : np.array of shape (batch_size, n_dcc, n_obs, n_strains) + + Returns + ------- + np.array of shape (batch_size, n_dcc) + + """ + strain_active = np.any(data, axis=2) + n_strain_obs = np.sum(strain_active, axis=2) # axis 3 is now 2 + + return n_strain_obs + + +def ss_prevalence(data): + """Calculate the prevalence of carriage among the observed individuals. + + Parameters + ---------- + data : np.array of shape (batch_size, n_dcc, n_obs, n_strains) + + Returns + ------- + np.array of shape (batch_size, n_dcc) + + """ + any_infection = np.any(data, axis=3) + n_infected = np.sum(any_infection, axis=2) + + return n_infected / data.shape[2] + + +def ss_prevalence_multi(data): + """Calculate the prevalence of multiple infections among the observed individuals. + + Parameters + ---------- + data : np.array of shape (batch_size, n_dcc, n_obs, n_strains) + + Returns + ------- + np.array of shape (batch_size, n_dcc) + + """ + n_infections = np.sum(data, axis=3) + n_multi_infections = np.sum(n_infections > 1, axis=2) + + return n_multi_infections / data.shape[2] + + +def distance(*summaries, observed): + """Calculate an L1-based distance between the simulated and observed summaries. + + Follows the simplified single-distance approach in: + Gutmann M U, Corander J (2016). Bayesian Optimization for Likelihood-Free Inference + of Simulator-Based Statistical Models. JMLR 17(125):1−47, 2016. + + Parameters + ---------- + *summaries : k np.arrays of shape (m, n) + observed : list of k np.arrays of shape (1, n) + + Returns + ------- + np.array of shape (m,) + + """ + summaries = np.stack(summaries) + observed = np.stack(observed) + n_ss, _, n_dcc = summaries.shape + + # scale summaries with max observed + obs_max = np.max(observed, axis=2, keepdims=True) + obs_max = np.where(obs_max == 0, 1, obs_max) + y = observed / obs_max + x = summaries / obs_max + + # sort to make comparison more robust + y = np.sort(y, axis=2) + x = np.sort(x, axis=2) + + # L1 norm divided by the dimension + dist = np.sum(np.abs(x - y), axis=(0, 2)) / (n_ss * n_dcc) + + return dist diff --git a/elfi/examples/lorenz.py b/elfi/examples/lorenz.py new file mode 100644 index 00000000..bdb8aab5 --- /dev/null +++ b/elfi/examples/lorenz.py @@ -0,0 +1,320 @@ +"""Example implementation of the Lorenz forecast model. + +References +---------- +- Dutta, R., Corander, J., Kaski, S. and Gutmann, M.U., 2016. + Likelihood-free inference by ratio estimation. arXiv preprint arXiv:1611.10242. + https://arxiv.org/abs/1611.10242 + +""" + +from functools import partial + +import numpy as np + +import elfi + + +def _lorenz_ode(y, params): + """Parametrized Lorenz 96 system defined by a coupled stochastic differential equations (SDE). + + Parameters + ---------- + y : numpy.ndarray of dimension (batch_size, n_obs) + Current state of the SDE. + params : list + The list of parameters needed to evaluate function. In this case it is + list of four elements - eta, theta1, theta2 and f. + + Returns + ------- + dy_dt : np.array + Rate of change of the SDE. + + """ + dy_dt = np.empty_like(y) + + eta = params[0] + theta1 = params[1] + theta2 = params[2] + + f = params[3] + + g = theta1 + y * theta2 + + dy_dt[:, 0] = -y[:, -2] * y[:, -1] + y[:, -1] * y[:, 1] - y[:, 0] + f - g[:, 0] + eta[:, 0] + + dy_dt[:, 1] = -y[:, -1] * y[:, 0] + y[:, 0] * y[:, 2] - y[:, 1] + f - g[:, 1] + eta[:, 1] + + dy_dt[:, 2:-1] = (-y[:, :-3] * y[:, 1:-2] + y[:, 1:-2] * y[:, 3:] - y[:, 2:-1] + f - g[:, 2:-1] + + eta[:, 2:-1]) + + dy_dt[:, -1] = (-y[:, -3] * y[:, -2] + y[:, -2] * y[:, 0] - y[:, -1] + f - g[:, -1] + + eta[:, -1]) + + return dy_dt + + +def runge_kutta_ode_solver(ode, time_step, y, params): + """4th order Runge-Kutta ODE solver. + + Carnahan, B., Luther, H. A., and Wilkes, J. O. (1969). + Applied Numerical Methods. Wiley, New York. + + Parameters + ---------- + ode : function + Ordinary differential equation function. In the Lorenz model it is SDE. + time_step : float + y : np.ndarray of dimension (batch_size, n_obs) + Current state of the time-series. + params : list of parameters + The parameters needed to evaluate the ode. In this case it is + list of four elements - eta, theta1, theta2 and f. + + Returns + ------- + np.ndarray + Resulting state initiated at y and satisfying ode solved by this solver. + + """ + k1 = time_step * ode(y, params) + + k2 = time_step * ode(y + k1 / 2, params) + + k3 = time_step * ode(y + k2 / 2, params) + + k4 = time_step * ode(y + k3, params) + + y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6 + + return y + + +def forecast_lorenz(theta1=None, theta2=None, f=10., phi=0.984, n_obs=40, n_timestep=160, + batch_size=1, initial_state=None, random_state=None, total_duration=4): + """Forecast Lorenz model. + + Wilks, D. S. (2005). Effects of stochastic parametrizations in the + Lorenz ’96 system. Quarterly Journal of the Royal Meteorological Society, + 131(606), 389–407. + + Parameters + ---------- + theta1, theta2: list or numpy.ndarray + Closure parameters. + phi : float, optional + This value is used to express stochastic forcing term. It should be configured according + to force term and eventually impacts to the result of eta. + More details in Wilks (2005) et al. + initial_state: numpy.ndarray, optional + Initial state value of the time-series. + f : float, optional + Force term + n_obs : int, optional + Size of the observed 1D grid + n_timestep : int, optional + Number of the time step intervals + batch_size : int, optional + random_state : np.random.RandomState, optional + total_duration : float, optional + + Returns + ------- + np.ndarray of size (b, n, m) which is (batch_size, time, n_obs) + The computed SDE with time series. + + """ + if not initial_state: + initial_state = np.tile([2.40711741e-01, 4.75597337e+00, 1.19145654e+01, 1.31324866e+00, + 2.82675744e+00, 3.96016971e+00, 2.10479504e+00, 5.47742826e+00, + 5.42519447e+00, -1.45166074e+00, 2.01991521e+00, 3.93873313e+00, + 8.22837848e+00, 4.89401702e+00, -5.66278973e+00, 1.58617220e+00, + -1.23849251e+00, -6.04649288e-01, 6.04132264e+00, 7.47588536e+00, + 1.82761402e+00, 3.19209639e+00, -7.58539653e-02, -6.00928508e-03, + 4.52902964e-01, 3.22063602e+00, 7.18613523e+00, 2.39210634e+00, + -2.65743666e+00, 2.32046235e-01, 1.28079141e+00, 4.23344286e+00, + 6.94213238e+00, -1.15939497e+00, -5.23037351e-01, 1.54618811e+00, + 1.77863869e+00, 3.30139201e+00, 7.47769309e+00, -3.91312909e-01], + (batch_size, 1)) + + y = initial_state + eta = 0 + + theta1 = np.asarray(theta1).reshape(-1, 1) + + theta2 = np.asarray(theta2).reshape(-1, 1) + + time_step = total_duration / n_timestep + + random_state = random_state or np.random + + time_series = np.empty(shape=(batch_size, n_timestep, n_obs)) + time_series[:, 0, :] = y + + for i in range(1, n_timestep): + e = random_state.normal(0, 1, y.shape) + eta = phi * eta + e * np.sqrt(1 - pow(phi, 2)) + params = (eta, theta1, theta2, f) + + y = runge_kutta_ode_solver(ode=_lorenz_ode, time_step=time_step, y=y, params=params) + time_series[:, i, :] = y + + return time_series + + +def get_model(true_params=None, seed_obs=None, initial_state=None, n_obs=40, f=10., phi=0.984, + total_duration=4): + """Return a complete Lorenz model in inference task. + + This is a simplified example that achieves reasonable predictions. + + Hakkarainen, J., Ilin, A., Solonen, A., Laine, M., Haario, H., Tamminen, + J., Oja, E., and Järvinen, H. (2012). On closure parameter estimation in + chaotic systems. Nonlinear Processes in Geophysics, 19(1), 127–143. + + Parameters + ---------- + true_params : list, optional + Parameters with which the observed data is generated. + seed_obs : int, optional + Seed for the observed data generation. + initial_state : ndarray + Initial state value of the time-series. + n_obs : int, optional + Number of observed variables + f : float, optional + Force term + phi : float, optional + This value is used to express stochastic forcing term. It should be configured according + to force term and eventually impacts to the result of eta. + More details in Wilks (2005) et al. + total_duration : float, optional + + Returns + ------- + m : elfi.ElfiModel + + """ + simulator = partial(forecast_lorenz, initial_state=initial_state, f=f, n_obs=n_obs, phi=phi, + total_duration=total_duration) + + if not true_params: + true_params = [2.0, 0.1] + + m = elfi.ElfiModel() + + y_obs = simulator(*true_params, random_state=np.random.RandomState(seed_obs)) + sumstats = [] + + elfi.Prior('uniform', 0.5, 3., model=m, name='theta1') + elfi.Prior('uniform', 0, 0.3, model=m, name='theta2') + elfi.Simulator(simulator, m['theta1'], m['theta2'], observed=y_obs, name='Lorenz') + + sumstats.append(elfi.Summary(mean, m['Lorenz'], name='Mean')) + + sumstats.append(elfi.Summary(var, m['Lorenz'], name='Var')) + + sumstats.append(elfi.Summary(autocov, m['Lorenz'], name='Autocov')) + + sumstats.append(elfi.Summary(cov, m['Lorenz'], name='Cov')) + + sumstats.append(elfi.Summary(xcov, m['Lorenz'], True, name='CrosscovPrev')) + + sumstats.append(elfi.Summary(xcov, m['Lorenz'], False, name='CrosscovNext')) + + elfi.Distance('euclidean', *sumstats, name='d') + + return m + + +def mean(x): + """Return the mean of Y_{k}. + + Parameters + ---------- + x : np.array of size (b, n, m) which is (batch_size, time, n_obs) + + Returns + ------- + np.array of size (b,) + The computed mean of statistics over time and space. + + """ + return np.mean(x, axis=(1, 2)) + + +def var(x): + """Return the variance of Y_{k}. + + Parameters + ---------- + x : np.array of size (b, n, m) which is (batch_size, time, n_obs) + + Returns + ------- + np.array of size (b,) + The average over space of computed variance with respect to time. + + """ + return np.mean(np.var(x, axis=1), axis=1) + + +def cov(x): + """Return the covariance of Y_{k} with its neighbour Y_{k+1}. + + Parameters + ---------- + x : np.array of size (b, n, m) which is (batch_size, time, n_obs) + + Returns + ------- + np.array of size (b,) + The average over space of computed covariance with respect to time. + + """ + x_next = np.roll(x, -1, axis=2) + return np.mean(np.mean((x - np.mean(x, keepdims=True, axis=1)) * + (x_next - np.mean(x_next, keepdims=True, axis=1)), + axis=1), axis=1) + + +def xcov(x, prev=True): + """Return the cross-covariance of Y_{k} with its neighbours from previous time step. + + Parameters + ---------- + x : np.array of size (b, n, m) which is (batch_size, time, n_obs) + prev : bool + The side of previous neighbour. True for previous neighbour, False for next. + + Returns + ------- + np.array of size (b,) + The average over space of computed cross-covariance with respect to time. + + """ + x_lag = np.roll(x, 1, axis=2) if prev else np.roll(x, -1, axis=2) + return np.mean((x[:, :-1, :] - np.mean(x[:, :-1, :], keepdims=True, axis=1)) * + (x_lag[:, 1:, :] - np.mean(x_lag[:, 1:, :], keepdims=True, axis=1)), + axis=(1, 2)) + + +def autocov(x): + """Return the auto-covariance with time lag 1. + + Parameters + ---------- + x : np.array of size (b, n, m) which is (batch_size, time, n_obs) + + Returns + ------- + C : np.array of size (b,) + The average over space of computed auto-covariance with respect to time. + + """ + c = np.mean((x[:, :-1, :] - np.mean(x[:, :-1, :], keepdims=True, axis=1)) * + (x[:, 1:, :] - np.mean(x[:, 1:, :], keepdims=True, axis=1)), + axis=(1, 2)) + + return c diff --git a/elfi/examples/lotka_volterra.py b/elfi/examples/lotka_volterra.py index f9cc6482..2180d922 100644 --- a/elfi/examples/lotka_volterra.py +++ b/elfi/examples/lotka_volterra.py @@ -141,7 +141,7 @@ 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, stochastic=True): +def get_model(n_obs=50, true_params=None, seed_obs=None, **kwargs): """Return a complete Lotka-Volterra model in inference task. Parameters @@ -159,13 +159,14 @@ def get_model(n_obs=50, true_params=None, seed_obs=None, stochastic=True): """ logger = logging.getLogger() - simulator = partial(lotka_volterra, n_obs=n_obs) if true_params is None: - true_params = [1.0, 0.005, 0.6, 50, 100, 10.] + true_params = [1.0, 0.005, 0.6, 50, 100, 10.] + + kwargs['n_obs'] = n_obs + y_obs = lotka_volterra(*true_params, random_state=np.random.RandomState(seed_obs), **kwargs) m = elfi.ElfiModel() - y_obs = simulator(*true_params, n_obs=n_obs, random_state=np.random.RandomState(seed_obs)) - sim_fn = partial(simulator, n_obs=n_obs) + sim_fn = partial(lotka_volterra, **kwargs) priors = [] sumstats = [] diff --git a/requirements-dev.txt b/requirements-dev.txt index 74ef217c..9dac5a62 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -3,6 +3,7 @@ pytest>=3.0.3 tox>=2.4.1 coverage>=4.2 pytest-cov>=2.4.0 +pytest-rerunfailures>=4.1 # Linting flake8>=3.0.4 diff --git a/tests/unit/test_client.py b/tests/unit/test_client.py index fb5cda04..20762df8 100644 --- a/tests/unit/test_client.py +++ b/tests/unit/test_client.py @@ -29,4 +29,15 @@ def test_batch_handler(simple_model): assert not np.array_equal(out0['k2'], out1['k2']) +def test_multiprocessing_kwargs(simple_model): + m = simple_model + num_proc = 2 + elfi.set_client('multiprocessing', num_processes=num_proc) + rej = elfi.Rejection(m['k1']) + assert rej.client.num_cores == num_proc + + elfi.set_client('multiprocessing', processes=num_proc) + rej = elfi.Rejection(m['k1']) + assert rej.client.num_cores == num_proc + # TODO: add testing that client is cleared from tasks after they are retrieved diff --git a/tests/unit/test_examples.py b/tests/unit/test_examples.py index f12c8a67..dff8d8f9 100644 --- a/tests/unit/test_examples.py +++ b/tests/unit/test_examples.py @@ -5,7 +5,7 @@ import pytest import elfi -from elfi.examples import bdm, bignk, gauss, gnk, lotka_volterra, ricker +from elfi.examples import bdm, bignk, gauss, gnk, lotka_volterra, ricker, daycare, lorenz def test_bdm(): @@ -72,6 +72,12 @@ def test_Ricker(): rej.sample(20) +def test_Lorenz(): + m = lorenz.get_model() + rej = elfi.Rejection(m, m['d'], batch_size=10) + rej.sample(20) + + def test_gnk(): m = gnk.get_model() rej = elfi.Rejection(m, m['d'], batch_size=10) @@ -85,6 +91,11 @@ def test_bignk(stats_summary=['ss_octile']): def test_Lotka_Volterra(): - m = lotka_volterra.get_model() + m = lotka_volterra.get_model(time_end=0.05) rej = elfi.Rejection(m, m['d'], batch_size=10) - rej.sample(20) + rej.sample(10, quantile=0.5) + +def test_daycare(): + m = daycare.get_model(time_end=0.05) + rej = elfi.Rejection(m['d'], batch_size=10) + rej.sample(10, quantile=0.5)