diff --git a/.gitignore b/.gitignore index 8ef9160b..b992e3ca 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,7 @@ -scratch_* -*ipynb_checkpoints -*.DS_Store -*__pycache__ -build/* -epsie.egg-info/* +build/ +test/ +epsie.egg-info/ +.DS_Store reinstall_epsie.sh +*.pyc +scratch_* diff --git a/epsie/chain/chain.py b/epsie/chain/chain.py index 4b9fe626..7b68db1b 100644 --- a/epsie/chain/chain.py +++ b/epsie/chain/chain.py @@ -537,7 +537,12 @@ def step(self): if self.transdimensional: current_pos.update({'_state': self._active_props}) # create a proposal and test it - proposal = self.proposal_dist.jump(current_pos) + # Some proposals may require the chain + if "differential_evolution" in self.proposal_dist.name: + proposal = self.proposal_dist.jump(current_pos, chain=self) + else: + proposal = self.proposal_dist.jump(current_pos) + self.proposed_position = proposal.copy() r = self.model(**proposal) diff --git a/epsie/proposals/__init__.py b/epsie/proposals/__init__.py index 306fb09f..73cba16a 100644 --- a/epsie/proposals/__init__.py +++ b/epsie/proposals/__init__.py @@ -18,6 +18,7 @@ from .joint import JointProposal from .normal import (Normal, AdaptiveNormal, SSAdaptiveNormal, ATAdaptiveNormal) +from .diff_evolution import DifferentialEvolution from .eigenvector import (Eigenvector, AdaptiveEigenvector) from .bounded_eigenvector import (BoundedEigenvector, AdaptiveBoundedEigenvector) diff --git a/epsie/proposals/base.py b/epsie/proposals/base.py index 38d5a35a..59063f8e 100644 --- a/epsie/proposals/base.py +++ b/epsie/proposals/base.py @@ -237,17 +237,17 @@ def _call_jump(self): return False return True - def jump(self, fromx): + def jump(self, fromx, **kwargs): """Depending on the current number of chain iterations and the ``jump_interval`` either calls the ``_jump`` method to provide a random sample or returns ``fromx``. """ if not self._call_jump(): return fromx - return self._jump(fromx) + return self._jump(fromx, **kwargs) @abstractmethod - def _jump(self, fromx): + def _jump(self, fromx, **kwargs): """This should provide random samples from the proposal distribution. Samples should be returned as a dictionary mapping parameters to diff --git a/epsie/proposals/diff_evolution.py b/epsie/proposals/diff_evolution.py new file mode 100644 index 00000000..ed885247 --- /dev/null +++ b/epsie/proposals/diff_evolution.py @@ -0,0 +1,88 @@ +# Copyright (C) 2020 Collin Capano, Richard Stiskalek +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +import numpy +from scipy import stats + +from .base import (BaseProposal, BaseAdaptiveSupport) +from .normal import Normal + + +class DifferentialEvolution(BaseProposal): + """ + Differential evolution proposal with fixed... + + This proposal may handle one or more parameters. + + Parameters + ---------- + parameters : (list of) str + The names of the parameters to produce proposals for. + + jump_interval : int, optional + The jump interval of the proposal, the proposal only gets called every + jump interval-th time. After ``jump_interval_duration`` number of + proposal steps elapses the proposal will again be called on every + chain iteration. By default ``jump_interval`` = 1. + jump_interval_duration : int, optional + The number of proposals steps during which values of ``jump_interval`` + other than 1 are used. After this elapses the proposal is called on + each iteration. + + Attributes + ---------- + + """ + + name = 'differential_evolution' + symmetric = True + + def __init__(self, parameters, jump_interval=1, + jump_interval_duration=None): + self.parameters = parameters + self.ndim = len(self.parameters) + self.set_jump_interval(jump_interval, jump_interval_duration) + self._scale = numpy.ones(self.ndim) * 2.38 * numpy.sqrt(self.ndim) + self._jump_stds = numpy.ones(self.ndim) + + @property + def state(self): + return {'random_state': self.random_state} + + def set_state(self, state): + self.random_state = state['random_state'] + + def _jump(self, fromx, chain): + # Randomly pick two points from the chain + N = chain.iteration + if N > chain.positions.size: + raise NotImplementedError("Account for having saved some points.") + if not N > 1: + raise NotImplementedError("Need to do a few initial steps with a " + "different proposal") + i, j = chain.random_generator.choice(range(N), size=2, replace=False) + + self._jump_stds[:] = 1. + self._jump_stds[:] = self.random_generator.normal(scale=self._scale) + + newpt = [fromx[p] + + self._jump_stds[k] + * (chain.positions[i][p] - chain.positions[j][p]) + for k, p in enumerate(self.parameters)] + + return dict(zip(self.parameters, newpt)) + + def _logpdf(self, xi, givenx): + return stats.norm.logpdf(x=self._jump_stds, scale=self._scale).sum() diff --git a/epsie/proposals/joint.py b/epsie/proposals/joint.py index 144b1c67..4ebb973b 100644 --- a/epsie/proposals/joint.py +++ b/epsie/proposals/joint.py @@ -72,14 +72,14 @@ def _update(self, chain): for prop in self.proposals: prop.update(chain) - def _jump(self, fromx): + def _jump(self, fromx, **kwargs): out = {} for prop in self.proposals: # we'll only pass the parameters that the proposal needs point = {p: fromx[p] for p in prop.parameters} if prop.transdimensional: point.update({'_state': fromx['_state']}) - out.update(prop.jump(point)) + out.update(prop.jump(point, **kwargs)) return out @property