diff --git a/.travis.yml b/.travis.yml index 2c69c280..c531c385 100644 --- a/.travis.yml +++ b/.travis.yml @@ -25,6 +25,6 @@ install: - pip install -e . script: - - ipcluster start -n 2 --daemon + - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then ipcluster start -n 2 --daemon ; fi #- travis_wait 20 make test - make test diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 2c669ea5..fbf1bfd0 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,19 @@ Changelog ========= + +0.7.4 (2019-03-07) +------------------ +- Add sampler option `algorithm` for bolfi-posterior-sampling +- Add a check whether the option given for `algorithm` is one if the implemented samplers +- Add metropolis sampler `algorithm=metropolis` for bolfi-posterior-sampling +- Add option `warmup` to metropolis-sampler +- Add a small test of metropolis-sampler +- Fix bug in plot_discrepancy for more than 6 parameters +- Implement plot_gp for BayesianOptimization classes for plotting discrepancies + and pair-wise contours in case when we have arbitrary number of parameters +- Fix lint + 0.7.3 (2018-08-30) ------------------ - Fix bug in plot_pairs which crashes in case of 1 parameter diff --git a/Dockerfile b/Dockerfile index 3f3e0b08..0c88d830 100644 --- a/Dockerfile +++ b/Dockerfile @@ -3,8 +3,14 @@ FROM ubuntu:bionic RUN apt-get update && DEBIAN_FRONTEND=noninteractive apt-get install -y graphviz make python3-scipy jupyter python3-matplotlib RUN pip3 install graphviz -WORKDIR /elfi -ADD . /elfi +ENV LANG C.UTF-8 +ENV LC_ALL C.UTF-8 + +WORKDIR /elfi0 +ADD . /elfi0 RUN pip3 install -e . RUN pip3 install -r requirements-dev.txt + +# Note: The created image contains a static version of ELFI. To use the live, up-to-date version +# you should mount the elfi directory when running the container ('make docker' does it for you). diff --git a/Makefile b/Makefile index 5b18667c..85918af6 100644 --- a/Makefile +++ b/Makefile @@ -103,3 +103,11 @@ install: clean ## install the package to the active Python's site-packages dev: install ## install the development requirements to the active Python's site-packages pip install -r requirements-dev.txt + +docker-build: ## build a docker image suitable for running ELFI + docker build --rm -t elfi . + +docker: ## run a docker container with a live elfi directory and publish port 8888 for Jupyter + docker run --rm -v ${PWD}:/elfi -w /elfi -it -p 8888:8888 elfi + # to run Jupyter from within the container: jupyter notebook --ip 0.0.0.0 --no-browser --allow-root + # and then from host open page: http://localhost:8888 diff --git a/README.md b/README.md index c1a04662..33919121 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -**Version 0.7.3 released!** See the CHANGELOG and [notebooks](https://github.com/elfi-dev/notebooks). +**Version 0.7.4 released!** See the CHANGELOG and [notebooks](https://github.com/elfi-dev/notebooks). **NOTE:** For the time being NetworkX 2 is incompatible with ELFI. @@ -77,15 +77,21 @@ pip install elfi ### Docker container -A simple Dockerfile for command-line interface is also provided. This is especially suitable for running tests. Please see [Docker documentation](https://docs.docker.com/) for details. +A simple Dockerfile with Jupyter support is also provided. This is especially suitable for running tests. Please see [Docker documentation](https://docs.docker.com/) for details. ``` git clone --depth 1 https://github.com/elfi-dev/elfi.git cd elfi -docker build -t elfi . -docker run -it elfi +make docker-build # builds the image with requirements for dev +make docker # runs a container with live elfi directory ``` +To open a Jupyter notebook, run +``` +jupyter notebook --ip 0.0.0.0 --no-browser --allow-root +``` +within the container and then on host open the page http://localhost:8888. + ### Potential problems with installation ELFI depends on several other Python packages, which have their own dependencies. @@ -101,13 +107,17 @@ Resolving these may sometimes go wrong: Citation -------- -If you wish to cite ELFI, please use the paper in [arXiv](https://arxiv.org/abs/1708.00707): +If you wish to cite ELFI, please use the paper in [JMLR](http://www.jmlr.org/papers/v19/17-374.html): ``` -@misc{1708.00707, -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 = {2018}, -Eprint = {arXiv:1708.00707}, +@article{JMLR:v19:17-374, + author = {Jarno Lintusaari and Henri Vuollekoski and Antti Kangasr{\"a}{\"a}si{\"o} and Kusti Skyt{\'e}n and Marko J{\"a}rvenp{\"a}{\"a} and Pekka Marttinen and Michael U. Gutmann and Aki Vehtari and Jukka Corander and Samuel Kaski}, + title = {ELFI: Engine for Likelihood-Free Inference}, + journal = {Journal of Machine Learning Research}, + year = {2018}, + volume = {19}, + number = {16}, + pages = {1-7}, + url = {http://jmlr.org/papers/v19/17-374.html} } ``` diff --git a/docs/index.rst b/docs/index.rst index d88f94a0..86a485f9 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -75,15 +75,19 @@ Additionally, ELFI integrates tools for visualization, model comparison, diagnos Citation -------- -If you wish to cite ELFI, please use the paper in arXiv_: +If you wish to cite ELFI, please use the paper in JMLR_: -.. _arXiv: https://arxiv.org/abs/1708.00707 +.. _JMLR: http://www.jmlr.org/papers/v19/17-374.html .. 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 Pekka Marttinen and Michael Gutmann and Aki Vehtari and Jukka Corander and Samuel Kaski}, - Title = {ELFI: Engine for Likelihood Free Inference}, - Year = {2018}, - Eprint = {arXiv:1708.00707}, + @article{JMLR:v19:17-374, + author = {Jarno Lintusaari and Henri Vuollekoski and Antti Kangasr{\"a}{\"a}si{\"o} and Kusti Skyt{\'e}n and Marko J{\"a}rvenp{\"a}{\"a} and Pekka Marttinen and Michael U. Gutmann and Aki Vehtari and Jukka Corander and Samuel Kaski}, + title = {ELFI: Engine for Likelihood-Free Inference}, + journal = {Journal of Machine Learning Research}, + year = {2018}, + volume = {19}, + number = {16}, + pages = {1-7}, + url = {http://jmlr.org/papers/v19/17-374.html} } diff --git a/docs/installation.rst b/docs/installation.rst index cd016a68..27b888f0 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -92,7 +92,7 @@ means the current folder. Docker container ---------------- -A simple Dockerfile for command-line interface is also provided. This is especially suitable +A simple Dockerfile with Jupyter support is also provided. This is especially suitable for running tests. Please see `Docker documentation`_ for details. .. _Docker documentation: https://docs.docker.com/ @@ -101,5 +101,15 @@ for running tests. Please see `Docker documentation`_ for details. git clone --depth 1 https://github.com/elfi-dev/elfi.git cd elfi - docker build -t elfi . - docker run -it elfi + make docker-build # builds the image with requirements for dev + make docker # runs a container with live elfi directory + +To open a Jupyter notebook, run + +.. code-block:: console + + jupyter notebook --ip 0.0.0.0 --no-browser --allow-root + +within the container and then on host open the page `localhost:8888`_. + +.. _localhost:8888: http://localhost:8888 diff --git a/elfi/__init__.py b/elfi/__init__.py index 745d5784..53ce03c9 100644 --- a/elfi/__init__.py +++ b/elfi/__init__.py @@ -26,4 +26,4 @@ __email__ = 'elfi-support@hiit.fi' # make sure __version_ is on the last non-empty line (read by setup.py) -__version__ = '0.7.3' +__version__ = '0.7.4' diff --git a/elfi/examples/daycare.py b/elfi/examples/daycare.py index 066ce2a1..ac3e0374 100644 --- a/elfi/examples/daycare.py +++ b/elfi/examples/daycare.py @@ -165,7 +165,7 @@ def get_model(true_params=None, seed_obs=None, **kwargs): """ logger = logging.getLogger() if true_params is None: - true_params = [3.6, 0.6, 0.1] + true_params = [3.6, 0.6, 0.1] m = elfi.ElfiModel() y_obs = daycare(*true_params, random_state=np.random.RandomState(seed_obs), **kwargs) diff --git a/elfi/loader.py b/elfi/loader.py index 738c7b49..eae8688e 100644 --- a/elfi/loader.py +++ b/elfi/loader.py @@ -157,7 +157,7 @@ def load(cls, context, compiled_net, batch_index): key = 'output' seed = context.seed - if seed is 'global': + if seed == 'global': # Get the random_state of the respective worker by delaying the evaluation random_state = get_np_random key = 'operation' diff --git a/elfi/methods/mcmc.py b/elfi/methods/mcmc.py index 7540d1de..cd87d71e 100644 --- a/elfi/methods/mcmc.py +++ b/elfi/methods/mcmc.py @@ -374,7 +374,7 @@ def _build_tree_nuts(params, momentum, log_slicevar, step, depth, log_joint0, ta mh_ratio, n_steps, is_div, is_out -def metropolis(n_samples, params0, target, sigma_proposals, seed=0): +def metropolis(n_samples, params0, target, sigma_proposals, warmup=0, seed=0): """Sample the target with a Metropolis Markov Chain Monte Carlo using Gaussian proposals. Parameters @@ -387,6 +387,8 @@ def metropolis(n_samples, params0, target, sigma_proposals, seed=0): The target log density to sample (possibly unnormalized). sigma_proposals : np.array Standard deviations for Gaussian proposals of each parameter. + warmup : int + Number of warmup samples. seed : int, optional Seed for pseudo-random number generator. @@ -397,22 +399,29 @@ def metropolis(n_samples, params0, target, sigma_proposals, seed=0): """ random_state = np.random.RandomState(seed) - samples = np.empty((n_samples + 1, ) + params0.shape) + samples = np.empty((n_samples + warmup + 1, ) + params0.shape) samples[0, :] = params0 target_current = target(params0) + if np.isinf(target_current): + raise ValueError( + "Metropolis: Bad initialization point {},logpdf -> -inf.".format(params0)) + n_accepted = 0 - for ii in range(1, n_samples + 1): + for ii in range(1, n_samples + warmup + 1): samples[ii, :] = samples[ii - 1, :] + sigma_proposals * random_state.randn(*params0.shape) target_prev = target_current target_current = target(samples[ii, :]) - - if np.exp(target_current - target_prev) < random_state.rand(): # reject proposal + if ((np.exp(target_current - target_prev) < random_state.rand()) + or np.isinf(target_current) + or np.isnan(target_current)): # reject proposal samples[ii, :] = samples[ii - 1, :] target_current = target_prev else: n_accepted += 1 logger.info( - "{}: Total acceptance ratio: {:.3f}".format(__name__, float(n_accepted) / n_samples)) - return samples[1:, :] + "{}: Total acceptance ratio: {:.3f}".format(__name__, + float(n_accepted) / (n_samples+warmup))) + + return samples[(1+warmup):, :] diff --git a/elfi/methods/parameter_inference.py b/elfi/methods/parameter_inference.py index 85e1e942..b9e0553d 100644 --- a/elfi/methods/parameter_inference.py +++ b/elfi/methods/parameter_inference.py @@ -808,7 +808,7 @@ def __init__(self, exploration_rate=10, batch_size=1, batches_per_acquisition=None, - async=False, + async_acq=False, **kwargs): """Initialize Bayesian optimization. @@ -839,7 +839,7 @@ def __init__(self, batches_per_acquisition : int, optional How many batches will be requested from the acquisition function at one go. Defaults to max_parallel_batches. - async : bool, optional + async_acq : bool, optional Allow acquisitions to be made asynchronously, i.e. do not wait for all the results from the previous acquisition before making the next. This can be more efficient with a large amount of workers (e.g. in cluster environments) but @@ -875,7 +875,7 @@ def __init__(self, self.n_initial_evidence = n_initial self.n_precomputed_evidence = n_precomputed self.update_interval = update_interval - self.async = async + self.async_acq = async_acq self.state['n_evidence'] = self.n_precomputed_evidence self.state['last_GP_update'] = self.n_initial_evidence @@ -1029,7 +1029,7 @@ def _allow_submit(self, batch_index): if not super(BayesianOptimization, self)._allow_submit(batch_index): return False - if self.async: + if self.async_acq: return True # Allow submitting freely as long we are still submitting initial evidence @@ -1115,22 +1115,37 @@ def acq(x): def plot_discrepancy(self, axes=None, **kwargs): """Plot acquired parameters vs. resulting discrepancy. - TODO: refactor + Parameters + ---------- + axes : plt.Axes or arraylike of plt.Axes + + Return + ------ + axes : np.array of plt.Axes + """ - n_plots = self.target_model.input_dim - ncols = kwargs.pop('ncols', 5) - kwargs['sharey'] = kwargs.get('sharey', True) - shape = (max(1, n_plots // ncols), min(n_plots, ncols)) - axes, kwargs = vis._create_axes(axes, shape, **kwargs) - axes = axes.ravel() + return vis.plot_discrepancy(self.target_model, self.parameter_names, axes=axes, **kwargs) - for ii in range(n_plots): - axes[ii].scatter(self.target_model._gp.X[:, ii], self.target_model._gp.Y[:, 0]) - axes[ii].set_xlabel(self.parameter_names[ii]) + def plot_gp(self, axes=None, resol=50, const=None, bounds=None, **kwargs): + """Plot pairwise relationships as a matrix with parameters vs. discrepancy. - axes[0].set_ylabel('Discrepancy') + Parameters + ---------- + axes : matplotlib.axes.Axes, optional + resol : int, optional + Resolution of the plotted grid. + const : np.array, optional + Values for parameters in plots where held constant. Defaults to minimum evidence. + bounds: list of tuples, optional + List of tuples for axis boundaries. - return axes + Returns + ------- + axes : np.array of plt.Axes + + """ + return vis.plot_gp(self.target_model, self.parameter_names, axes, + resol, const, bounds, **kwargs) class BOLFI(BayesianOptimization): @@ -1203,6 +1218,7 @@ def sample(self, threshold=None, initials=None, algorithm='nuts', + sigma_proposals=None, n_evidence=None, **kwargs): r"""Sample the posterior distribution of BOLFI. @@ -1232,7 +1248,10 @@ def sample(self, Initial values for the sampled parameters for each chain. Defaults to best evidence points. algorithm : string, optional - Sampling algorithm to use. Currently only 'nuts' is supported. + Sampling algorithm to use. Currently 'nuts'(default) and 'metropolis' are supported. + sigma_proposals : np.array + Standard deviations for Gaussian proposals of each parameter for Metropolis + Markov Chain sampler. n_evidence : int If the regression model is not fitted yet, specify the amount of evidence @@ -1244,7 +1263,9 @@ def sample(self, if self.state['n_batches'] == 0: self.fit(n_evidence) - # TODO: other MCMC algorithms + # TODO: add more MCMC algorithms + if algorithm not in ['nuts', 'metropolis']: + raise ValueError("Unknown posterior sampler.") posterior = self.extract_posterior(threshold) warmup = warmup or n_samples // 2 @@ -1261,6 +1282,13 @@ 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.") # sampling is embarrassingly parallel, so depending on self.client this may parallelize for ii in range(n_chains): @@ -1272,16 +1300,30 @@ def sample(self, raise ValueError( "BOLFI.sample: Cannot find enough acceptable initialization points!") - tasks_ids.append( - self.client.apply( - mcmc.nuts, - n_samples, - initials[ii_initial], - posterior.logpdf, - posterior.gradient_logpdf, - n_adapt=warmup, - seed=seed, - **kwargs)) + if algorithm == 'nuts': + tasks_ids.append( + self.client.apply( + mcmc.nuts, + n_samples, + initials[ii_initial], + posterior.logpdf, + posterior.gradient_logpdf, + n_adapt=warmup, + seed=seed, + **kwargs)) + + elif algorithm == 'metropolis': + tasks_ids.append( + self.client.apply( + mcmc.metropolis, + n_samples, + initials[ii_initial], + posterior.logpdf, + sigma_proposals, + warmup, + seed=seed, + **kwargs)) + ii_initial += 1 # get results from completed tasks or run sampling (client-specific) @@ -1290,14 +1332,12 @@ def sample(self, chains.append(self.client.get_result(id)) chains = np.asarray(chains) - 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): print(node, mcmc.eff_sample_size(chains[:, :, ii]), mcmc.gelman_rubin(chains[:, :, ii])) - self.target_model.is_sampling = False return BolfiSample( diff --git a/elfi/visualization/visualization.py b/elfi/visualization/visualization.py index 89b94056..1d56fb47 100644 --- a/elfi/visualization/visualization.py +++ b/elfi/visualization/visualization.py @@ -359,3 +359,93 @@ def plot_params_vs_node(node, n_samples=100, func=None, seed=None, axes=None, ** axes[idx].set_axis_off() return axes + + +def plot_discrepancy(gp, parameter_names, axes=None, **kwargs): + """Plot acquired parameters vs. resulting discrepancy. + + Parameters + ---------- + axes : plt.Axes or arraylike of plt.Axes + gp : GPyRegression target model, required + parameter_names : dict, required + Parameter names from model.parameters dict('parameter_name':(lower, upper), ... )` + + Returns + ------- + axes : np.array of plt.Axes + + """ + n_plots = gp.input_dim + ncols = len(gp.bounds) if len(gp.bounds) < 5 else 5 + ncols = kwargs.pop('ncols', ncols) + kwargs['sharey'] = kwargs.get('sharey', True) + if n_plots > 10: + shape = (1 + (1 + n_plots) // (ncols + 1), ncols) + else: + shape = (1 + n_plots // (ncols + 1), ncols) + axes, kwargs = _create_axes(axes, shape, **kwargs) + axes = axes.ravel() + + for ii in range(n_plots): + axes[ii].scatter(gp.X[:, ii], gp.Y[:, 0], **kwargs) + axes[ii].set_xlabel(parameter_names[ii]) + if ii % ncols == 0: + axes[ii].set_ylabel('Discrepancy') + + for idx in range(len(parameter_names), len(axes)): + axes[idx].set_axis_off() + + return axes + + +def plot_gp(gp, parameter_names, axes=None, resol=50, const=None, bounds=None, **kwargs): + """Plot pairwise relationships as a matrix with parameters vs. discrepancy. + + Parameters + ---------- + gp : GPyRegression, required + parameter_names : list, required + Parameter names in format ['mu_0', 'mu_1', ..] + axes : plt.Axes or arraylike of plt.Axes + resol : int, optional + Resolution of the plotted grid. + const : np.array, optional + Values for parameters in plots where held constant. Defaults to minimum evidence. + bounds: list of tuples, optional + List of tuples for axis boundaries. + + Returns + ------- + axes : np.array of plt.Axes + + """ + n_plots = gp.input_dim + shape = (n_plots, n_plots) + axes, kwargs = _create_axes(axes, shape, **kwargs) + + x_evidence = gp.X + y_evidence = gp.Y + if const is None: + const = x_evidence[np.argmin(y_evidence), :] + bounds = bounds or gp.bounds + + for ix in range(n_plots): + for jy in range(n_plots): + if ix == jy: + axes[jy, ix].scatter(x_evidence[:, ix], y_evidence) + axes[jy, ix].set_xlim(bounds[ix]) + axes[jy, ix].set_ylabel('Discrepancy') + else: + x1 = np.linspace(bounds[ix][0], bounds[ix][1], resol) + y1 = np.linspace(bounds[jy][0], bounds[jy][1], resol) + x, y = np.meshgrid(x1, y1) + predictors = np.tile(const, (resol * resol, 1)) + predictors[:, ix] = x.ravel() + predictors[:, jy] = y.ravel() + z = gp.predict_mean(predictors).reshape(resol, resol) + axes[jy, ix].contourf(x, y, z) + axes[jy, ix].set_ylabel(parameter_names[jy]) + axes[jy, ix].set_xlabel(parameter_names[ix]) + + return axes diff --git a/requirements-dev.txt b/requirements-dev.txt index 9dac5a62..49f47a7f 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,5 +1,5 @@ # Testing -pytest>=3.0.3 +pytest>=3.8 tox>=2.4.1 coverage>=4.2 pytest-cov>=2.4.0 diff --git a/setup.py b/setup.py index 2ef202e4..5c983047 100644 --- a/setup.py +++ b/setup.py @@ -27,7 +27,7 @@ url='http://elfi.readthedocs.io', install_requires=requirements, extras_require=optionals, - description='Modular ABC inference framework for python', + description='ELFI - Engine for Likelihood-free Inference', long_description=(open('docs/description.rst').read()), license='BSD', classifiers=[ diff --git a/tests/unit/test_bo.py b/tests/unit/test_bo.py index eb16e116..a8414d7b 100644 --- a/tests/unit/test_bo.py +++ b/tests/unit/test_bo.py @@ -45,7 +45,7 @@ def test_BO(ma2): def test_async(ma2): bounds = {n: (-2, 2) for n in ma2.parameter_names} bo = elfi.BayesianOptimization( - ma2, 'd', initial_evidence=0, update_interval=2, batch_size=2, bounds=bounds, async=True) + ma2, 'd', initial_evidence=0, update_interval=2, batch_size=2, bounds=bounds, async_acq=True) n_samples = 5 bo.infer(n_samples) diff --git a/tests/unit/test_methods.py b/tests/unit/test_methods.py index 1ad77ed0..be7f6877 100644 --- a/tests/unit/test_methods.py +++ b/tests/unit/test_methods.py @@ -73,9 +73,13 @@ def test_BOLFI_short(ma2, distribution_test): n_samples = 10 n_chains = 2 - res_sampling = bolfi.sample(n_samples, n_chains=n_chains) - assert res_sampling.samples_array.shape[1] == 2 - assert len(res_sampling.samples_array) == n_samples // 2 * n_chains + res_sampling_nuts = bolfi.sample(n_samples, n_chains=n_chains) + 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)) + assert res_sampling_metropolis.samples_array.shape[1] == 2 + assert len(res_sampling_metropolis.samples_array) == n_samples // 2 * n_chains # check the cached predictions for RBF x = np.random.random((1, 2))