From 40ec92c902673402fdaef0c6c9f06f91dcc1c263 Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Sat, 30 Sep 2023 21:42:31 +0200 Subject: [PATCH 01/14] document hierarchical optimization sources; add aesara+jax to objective docs --- doc/api.rst | 2 ++ pypesto/hierarchical/__init__.py | 27 ++++++++++++++++++++++++++- tox.ini | 2 +- 3 files changed, 29 insertions(+), 2 deletions(-) diff --git a/doc/api.rst b/doc/api.rst index f23a6a0cb..470318cd0 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -12,6 +12,8 @@ API reference pypesto.history pypesto.logging pypesto.objective + pypesto.objective.aesara + pypesto.objective.jax pypesto.objective.julia pypesto.optimize pypesto.petab diff --git a/pypesto/hierarchical/__init__.py b/pypesto/hierarchical/__init__.py index fc326f25e..74caae7d8 100644 --- a/pypesto/hierarchical/__init__.py +++ b/pypesto/hierarchical/__init__.py @@ -2,7 +2,32 @@ Hierarchical ============ -Hierarchical optimization sub-package. +Contains an implementation of the hierarchical optimization approach, which +decomposes the parameter estimation problem into an outer and an inner problem. +In the outer problem, only dynamic parameters are optimized. +In the inner problem, conditional on the outer solution, static parameters are +optimized. +Static parameters can be parameters affecting directly the model observables, +such as scaling factors, offsets, and noise parameters. + +Hierarchical optimization has the advantage that the outer problem is typically +less complex than the full problem, and thus can be solved more efficiently. +Further, the inner problem can often be solved analytically, which is more +efficient. +Thus, hierarchical optimization can be used to speed up parameter estimation, +finding optimal values more efficiently and reliably. + +The implementation in this package is based on: + +* Loos et al. 2018 (https://doi.org/10.1093/bioinformatics/bty514), + who give an analytic solution for the inner problem for scaling factors and + noise standard deviations, for Gaussian and Laplace noise, using forward + sensitivity analysis (FSA). +* Schmiester et al. 2020 (https://doi.org/10.1093/bioinformatics/btz581), + who give an analytic solution for the inner problem for scaling factors, + offsets and noise standard deviations, for Gaussian and Laplace noise, + using adjoint sensitivity analysis (ASA). ASA allows to calculate gradients + substantially more efficiently in high dimension. """ from .calculator import HierarchicalAmiciCalculator diff --git a/tox.ini b/tox.ini index 9866525c1..8773bafa0 100644 --- a/tox.ini +++ b/tox.ini @@ -171,7 +171,7 @@ description = [testenv:doc] extras = - doc,amici,petab + doc,amici,petab,aesara,jax commands = sphinx-build -W -b html doc/ doc/_build/html description = From 9db727262042799f411ae7e9b55d5a3d7ff28d8f Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Sun, 1 Oct 2023 16:05:42 +0200 Subject: [PATCH 02/14] document origin of the mh + pt samplers --- pypesto/sample/adaptive_metropolis.py | 58 +++++++++++++++++-- pypesto/sample/adaptive_parallel_tempering.py | 19 +++++- pypesto/sample/metropolis.py | 25 +++++++- pypesto/sample/parallel_tempering.py | 21 ++++++- 4 files changed, 116 insertions(+), 7 deletions(-) diff --git a/pypesto/sample/adaptive_metropolis.py b/pypesto/sample/adaptive_metropolis.py index cbe9b1825..7e822b4ad 100644 --- a/pypesto/sample/adaptive_metropolis.py +++ b/pypesto/sample/adaptive_metropolis.py @@ -8,7 +8,54 @@ class AdaptiveMetropolisSampler(MetropolisSampler): - """Metropolis-Hastings sampler with adaptive proposal covariance.""" + """Metropolis-Hastings sampler with adaptive proposal covariance. + + A core problem of the standard Metropolis-Hastings sampler is + its fixed proposal distribution, which must be manually tuned. + This sampler adapts the proposal distribution during the sampling + process based on previous samples. + It adapts the correlation structure and the scaling factor of the + proposal distribution. + For both parts, there exist a variety of methods, see + + * Ballnus et al. 2017. + Comprehensive benchmarking of Markov chain Monte Carlo methods for + dynamical systems + (https://doi.org/10.1186/s12918-017-0433-1) + * Andrieu et al. 2008. + A tutorial on adaptive MCMC + (https://doi.org/10.1007/s11222-008-9110-y) + + for a review. + + Here, we approximate the covariance matrix via a weighted average of + current and earlier samples, + with a decay factor determining the relative contribution of the + current sample and earlier ones to the weighted average of mean and + covariance. + The scaling factor we aim to converge to a fixed target acceptance rate + of 0.234, as suggested by theoretical results. + The implementation is based on: + + * Lacki et al. 2015. + State-dependent swap strategies and automatic reduction of number of + temperatures in adaptive parallel tempering algorithm + (https://doi.org/10.1007/s11222-015-9579-0) + * Miasojedow et al. 2013. + An adaptive parallel tempering algorithm + (https://doi.org/10.1080/10618600.2013.778779) + + In turn, these are based on adaptive MCMC as discussed in: + + * Haario et al. 2001. + An adaptive Metropolis algorithm + (https://doi.org/10.2307/3318737) + + For reference matlab implementations see: + + * https://github.com/ICB-DCM/PESTO/blob/master/private/performPT.m + * https://github.com/ICB-DCM/PESTO/blob/master/private/updateStatistics.m + """ def __init__(self, options: Dict = None): super().__init__(options) @@ -79,7 +126,7 @@ def _update_proposal( decay_constant=decay_constant, ) - # compute covariance scaling factor + # compute covariance scaling factor based on the target acceptance rate self._cov_scale *= np.exp( (np.exp(log_p_acc) - target_acceptance_rate) / np.power(n_sample_cur + 1, decay_constant) @@ -100,8 +147,11 @@ def update_history_statistics( n_cur_sample: int, decay_constant: float, ) -> Tuple[np.ndarray, np.ndarray]: - """ - Update sampling statistics. + """Updatae sampling mean and covariante matrix via weighted average. + + Update sampling mean and covariance matrix based on the previous + estimate and the most recent sample via a weighted average, + with gradually decaying update rate. Parameters ---------- diff --git a/pypesto/sample/adaptive_parallel_tempering.py b/pypesto/sample/adaptive_parallel_tempering.py index 788f752b4..2475db3e3 100644 --- a/pypesto/sample/adaptive_parallel_tempering.py +++ b/pypesto/sample/adaptive_parallel_tempering.py @@ -7,7 +7,24 @@ class AdaptiveParallelTemperingSampler(ParallelTemperingSampler): - """Parallel tempering sampler with adaptive temperature adaptation.""" + """Parallel tempering sampler with adaptive temperature adaptation. + + Compared to the base class, this sampler adapts the temperatures + during the sampling process. + This both simplifies the setup as it avoids manual tuning, + and improves the performance as the temperatures are adapted to the + current state of the chains. + + This implementation is based on: + + * Vousden et al. 2016. + Dynamic temperature selection for parallel tempering in Markov chain + Monte Carlo simulations + (https://doi.org/10.1093/mnras/stv2422), + + via a matlab reference implementation + (https://github.com/ICB-DCM/PESTO/blob/master/private/performPT.m). + """ @classmethod def default_options(cls) -> Dict: diff --git a/pypesto/sample/metropolis.py b/pypesto/sample/metropolis.py index 88a8ca41c..09495f534 100644 --- a/pypesto/sample/metropolis.py +++ b/pypesto/sample/metropolis.py @@ -11,7 +11,30 @@ class MetropolisSampler(InternalSampler): - """Simple Metropolis-Hastings sampler with fixed proposal variance.""" + """Simple Metropolis-Hastings sampler with fixed proposal variance. + + The Metropolis-Hastings sampler is a Markov chain Monte Carlo (MCMC) + method generating a sequence of samples from a probability + distribution. + + This class implements a simple Metropolis algorithm with fixed + symmetric Gaussian proposal distribution. + + For the underlying original publication, see: + + * Metropolis et al. 1953. + Equation of State Calculations by Fast Computing Machines + (https://doi.org/10.1063/1.1699114) + * Hastings 1970. + Monte Carlo sampling methods using Markov chains and their + applications + (https://doi.org/10.1093/biomet/57.1.97) + + For reference matlab implementations see: + + * https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm + * https://github.com/ICB-DCM/PESTO/blob/master/private/performPT.m + """ def __init__(self, options: Dict = None): super().__init__(options) diff --git a/pypesto/sample/parallel_tempering.py b/pypesto/sample/parallel_tempering.py index e82ab9e34..31fd95089 100644 --- a/pypesto/sample/parallel_tempering.py +++ b/pypesto/sample/parallel_tempering.py @@ -10,7 +10,26 @@ class ParallelTemperingSampler(Sampler): - """Simple parallel tempering sampler.""" + """Simple parallel tempering sampler. + + Parallem tempering is a Markov chain Monte Carlo (MCMC) method that + uses multiple chains with different temperatures to sample from a + probability distribution. + The chains are coupled by swapping samples between them. + This allows to sample from distributions with multiple modes more + efficiently, as high-temperature chains can jump between modes, while + low-temperature chains can sample the modes more precisely. + + This implementation is based on: + + * Vousden et al. 2016. + Dynamic temperature selection for parallel tempering in Markov chain + Monte Carlo simulations + (https://doi.org/10.1093/mnras/stv2422), + + via a matlab-based reference implementation + (https://github.com/ICB-DCM/PESTO/blob/master/private/performPT.m). + """ def __init__( self, From aac11dff21f2f924407f74b5c80e27dc8254e224 Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Sun, 1 Oct 2023 16:07:13 +0200 Subject: [PATCH 03/14] fix typo --- pypesto/sample/adaptive_metropolis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypesto/sample/adaptive_metropolis.py b/pypesto/sample/adaptive_metropolis.py index 7e822b4ad..37846e757 100644 --- a/pypesto/sample/adaptive_metropolis.py +++ b/pypesto/sample/adaptive_metropolis.py @@ -147,7 +147,7 @@ def update_history_statistics( n_cur_sample: int, decay_constant: float, ) -> Tuple[np.ndarray, np.ndarray]: - """Updatae sampling mean and covariante matrix via weighted average. + """Update sampling mean and covariante matrix via weighted average. Update sampling mean and covariance matrix based on the previous estimate and the most recent sample via a weighted average, From 9bc34e24affff4859e9511caf3e101e04eb19ee0 Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Sun, 1 Oct 2023 16:07:28 +0200 Subject: [PATCH 04/14] fix typo --- pypesto/sample/adaptive_metropolis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypesto/sample/adaptive_metropolis.py b/pypesto/sample/adaptive_metropolis.py index 37846e757..6b29bd804 100644 --- a/pypesto/sample/adaptive_metropolis.py +++ b/pypesto/sample/adaptive_metropolis.py @@ -147,7 +147,7 @@ def update_history_statistics( n_cur_sample: int, decay_constant: float, ) -> Tuple[np.ndarray, np.ndarray]: - """Update sampling mean and covariante matrix via weighted average. + """Update sampling mean and covariance matrix via weighted average. Update sampling mean and covariance matrix based on the previous estimate and the most recent sample via a weighted average, From d01ec35191c2afa1da5ef9e0ee9acb5e5940cc23 Mon Sep 17 00:00:00 2001 From: Paul Jonas Jost <70631928+PaulJonasJost@users.noreply.github.com> Date: Mon, 2 Oct 2023 09:23:14 +0200 Subject: [PATCH 05/14] Update pypesto/sample/parallel_tempering.py Typo --- pypesto/sample/parallel_tempering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypesto/sample/parallel_tempering.py b/pypesto/sample/parallel_tempering.py index 31fd95089..44a9f4453 100644 --- a/pypesto/sample/parallel_tempering.py +++ b/pypesto/sample/parallel_tempering.py @@ -12,7 +12,7 @@ class ParallelTemperingSampler(Sampler): """Simple parallel tempering sampler. - Parallem tempering is a Markov chain Monte Carlo (MCMC) method that + Parallel tempering is a Markov chain Monte Carlo (MCMC) method that uses multiple chains with different temperatures to sample from a probability distribution. The chains are coupled by swapping samples between them. From 69ba85bff388ea3fb33a9b3d8b2b48c5bee9642e Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 10:22:11 +0200 Subject: [PATCH 06/14] streamline doc deps --- .readthedocs.yml | 1 - setup.cfg | 4 ++++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index 4a2a4d3c1..b75a9c82b 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -17,7 +17,6 @@ python: path: . extra_requirements: - doc - - amici build: os: "ubuntu-20.04" diff --git a/setup.cfg b/setup.cfg index 6194ed8bd..a98b85b3a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -143,6 +143,10 @@ doc = ipykernel >= 6.15.1 %(select)s %(fides)s + %(amici)s + %(petab)s + %(aesara)s + %(jax)s example = notebook >= 6.1.4 select = From a030296c68153b5b548edac880f2d97ef505769f Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 10:33:58 +0200 Subject: [PATCH 07/14] update setup req versions --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ed59d2789..06bda6a3f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,8 +4,8 @@ [build-system] requires = [ - "setuptools >= 52.0.0", - "wheel >= 0.36.2", + "setuptools >= 68.2.0", + "wheel >= 0.41.2", ] build-backend = "setuptools.build_meta" From b8ee596d4ece956be0f680aa9ca803b4d5b9611f Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 10:57:22 +0200 Subject: [PATCH 08/14] try cyipopt fix --- tox.ini | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tox.ini b/tox.ini index 8773bafa0..f6a365cd3 100644 --- a/tox.ini +++ b/tox.ini @@ -124,11 +124,12 @@ description = Test model selection functionality [testenv:notebooks1] -allowlist_externals = bash +allowlist_externals = bash,conda extras = example,amici,petab,pyswarm,pymc3,cmaes,nlopt,fides commands = # workaround as ipopt has incomplete build - pip install git+https://github.com/mechmotum/cyipopt.git@master + # pip install git+https://github.com/mechmotum/cyipopt.git@master + conda install -c conda-forge cyipopt bash test/run_notebook.sh 1 description = Run notebooks 1 From ddccd9ab2af9144c56fc6932afa0a83eb2049a3d Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 11:10:15 +0200 Subject: [PATCH 09/14] test --- tox.ini | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tox.ini b/tox.ini index f6a365cd3..ca527a6c0 100644 --- a/tox.ini +++ b/tox.ini @@ -95,10 +95,12 @@ description = Test Julia interface [testenv:optimize] +allowlist_externals = conda extras = test,dlib,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab commands = # workaround as ipopt has incomplete build - pip install git+https://github.com/mechmotum/cyipopt.git@master + # pip install git+https://github.com/mechmotum/cyipopt.git@master + conda install -c conda-forge cyipopt pytest --cov=pypesto --cov-report=xml --cov-append \ test/optimize -s description = From 009ef5a8beead8454c69eb2be45f6c8570ef5c60 Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 12:02:11 +0200 Subject: [PATCH 10/14] tesT --- setup.cfg | 2 +- tox.ini | 7 ++----- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/setup.cfg b/setup.cfg index a98b85b3a..573406d37 100644 --- a/setup.cfg +++ b/setup.cfg @@ -96,7 +96,7 @@ amici = petab = petab >= 0.2.0 ipopt = - ipopt + cyipopt >= 1.3.0 dlib = dlib >= 19.19.0 nlopt = diff --git a/tox.ini b/tox.ini index ca527a6c0..f09ffc206 100644 --- a/tox.ini +++ b/tox.ini @@ -95,12 +95,10 @@ description = Test Julia interface [testenv:optimize] -allowlist_externals = conda -extras = test,dlib,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab +extras = test,dlib,ipopt,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab commands = # workaround as ipopt has incomplete build # pip install git+https://github.com/mechmotum/cyipopt.git@master - conda install -c conda-forge cyipopt pytest --cov=pypesto --cov-report=xml --cov-append \ test/optimize -s description = @@ -126,12 +124,11 @@ description = Test model selection functionality [testenv:notebooks1] -allowlist_externals = bash,conda +allowlist_externals = bash extras = example,amici,petab,pyswarm,pymc3,cmaes,nlopt,fides commands = # workaround as ipopt has incomplete build # pip install git+https://github.com/mechmotum/cyipopt.git@master - conda install -c conda-forge cyipopt bash test/run_notebook.sh 1 description = Run notebooks 1 From 841e25cd575baa94e75609413d0e7870e8c11fdc Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 12:05:59 +0200 Subject: [PATCH 11/14] test --- tox.ini | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tox.ini b/tox.ini index f09ffc206..65f5562e0 100644 --- a/tox.ini +++ b/tox.ini @@ -95,10 +95,11 @@ description = Test Julia interface [testenv:optimize] -extras = test,dlib,ipopt,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab +extras = test,dlib,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab commands = # workaround as ipopt has incomplete build # pip install git+https://github.com/mechmotum/cyipopt.git@master + pip install --use-pep517 cyipopt pytest --cov=pypesto --cov-report=xml --cov-append \ test/optimize -s description = From 9dd23c6401577e3be4eac6f26a97c5c736100f10 Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 14:00:18 +0200 Subject: [PATCH 12/14] fix test_visualize::test_optimization_scatter_with_x_None --- pypesto/visualize/parameters.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pypesto/visualize/parameters.py b/pypesto/visualize/parameters.py index 52509116f..d175c25b6 100644 --- a/pypesto/visualize/parameters.py +++ b/pypesto/visualize/parameters.py @@ -561,7 +561,7 @@ def optimization_scatter( parameter_indices = process_parameter_indices( parameter_indices=parameter_indices, result=result ) - # remove all start indices, that encounter an inf value at the start + # remove all start indices that encounter an inf value at the start # resulting in optimize_result[start]["x"] being None start_indices_finite = start_indices[ [ @@ -570,7 +570,7 @@ def optimization_scatter( ] ] # compare start_indices with start_indices_finite and log a warning - if not np.all(start_indices == start_indices_finite): + if len(start_indices) != len(start_indices_finite): logger.warning( 'Some start indices were removed due to inf values at the start.' ) From 4f1e49ea61fe6f8bf6c40c97e51179aaa7790b2e Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 14:09:29 +0200 Subject: [PATCH 13/14] test --- .github/workflows/ci.yml | 3 ++- tox.ini | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ec4710d23..1847c1a2e 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -189,7 +189,8 @@ jobs: strategy: matrix: # update to 3.11 when nlopt allows - python-version: ['3.9', '3.10'] + # ipopt does not work on 3.9 (https://github.com/mechmotum/cyipopt/issues/225) + python-version: ['3.10', '3.11'] steps: - name: Check out repository diff --git a/tox.ini b/tox.ini index 65f5562e0..3c40cc937 100644 --- a/tox.ini +++ b/tox.ini @@ -95,11 +95,11 @@ description = Test Julia interface [testenv:optimize] -extras = test,dlib,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab +extras = test,dlib,ipopt,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab commands = # workaround as ipopt has incomplete build # pip install git+https://github.com/mechmotum/cyipopt.git@master - pip install --use-pep517 cyipopt + # pip install --use-pep517 cyipopt pytest --cov=pypesto --cov-report=xml --cov-append \ test/optimize -s description = From b0745c22ee740f127f921281e78dc8e894c6b112 Mon Sep 17 00:00:00 2001 From: Yannik Schaelte Date: Mon, 2 Oct 2023 14:25:41 +0200 Subject: [PATCH 14/14] fixup --- .github/workflows/ci.yml | 3 +-- tox.ini | 5 ----- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1847c1a2e..f3fb62908 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -188,9 +188,8 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - # update to 3.11 when nlopt allows # ipopt does not work on 3.9 (https://github.com/mechmotum/cyipopt/issues/225) - python-version: ['3.10', '3.11'] + python-version: ['3.11'] steps: - name: Check out repository diff --git a/tox.ini b/tox.ini index 3c40cc937..9477f90f0 100644 --- a/tox.ini +++ b/tox.ini @@ -97,9 +97,6 @@ description = [testenv:optimize] extras = test,dlib,ipopt,pyswarm,cmaes,nlopt,fides,mpi,pyswarms,petab commands = - # workaround as ipopt has incomplete build - # pip install git+https://github.com/mechmotum/cyipopt.git@master - # pip install --use-pep517 cyipopt pytest --cov=pypesto --cov-report=xml --cov-append \ test/optimize -s description = @@ -128,8 +125,6 @@ description = allowlist_externals = bash extras = example,amici,petab,pyswarm,pymc3,cmaes,nlopt,fides commands = - # workaround as ipopt has incomplete build - # pip install git+https://github.com/mechmotum/cyipopt.git@master bash test/run_notebook.sh 1 description = Run notebooks 1