.. _algorithms:

Algorithms
==========

Summary of Available Algorithms
-------------------------------

+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| Algorithm                   | Class            | Parallelization | Applications                                                              |
+=============================+==================+=================+===========================================================================+
| `Differential Evolution`_   | Population-based | Synchronous or  | General-purpose parameter fitting                                         |
|                             |                  | Asynchronous    |                                                                           |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Scatter Search`_           | Population-based | Synchronous     | General-purpose parameter fitting, especially difficult problems with high| 
|                             |                  |                 | dimensions or many local minima                                           |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Particle Swarm`_           | Population-based | Asynchronous    | Fitting models with high variability in runtime                           | 
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Metropolis-Hastings MCMC`_ | Metropolis       | Independent     | Finding probability distributions of parameters (deprecated)              |
|                             | sampling         | Markov Chains   |                                                                           |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Simulated Annealing`_      | Metropolis       | Independent     | Problem-specific applications (deprecated)                                |
|                             | sampling         | Markov Chains   |                                                                           |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Parallel Tempering`_       | Metropolis       | Synchronized    | Finding probability distributions in challenging probablity landscapes    |
|                             | sampling         | Markov Chains   | (deprecated)                                                              |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Simplex`_                  | Local search     | Synchronous     | Local optimization, or refinement of a result from another algorithm.     |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Powell`_                   | Local search     | Synchronous     | Local optimization, or refinement of a result from another algorithm.     |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `CMA-ES`_                   | Population-based | Synchronous     | Global or local optimization, or refinement; robust on ill-conditioned    |
|                             |                  |                 | objectives.                                                               |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `Adaptive MCMC`_            | Metropolis       | Independent     | Finding probability distributions in challenging probablity landscapes    |
|                             | sampling         | Markov Chains   |                                                                           |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+
| `DREAM`_                    | Hybrid           | Synchronous     | Finding probability distributions with accelerated convergence            |
|                             | Population /     |                 |                                                                           |
|                             | Metropolis       |                 |                                                                           |
+-----------------------------+------------------+-----------------+---------------------------------------------------------------------------+


In addition to the population-based and Metropolis samplers in the table above,
PyBNF provides gradient-based optimizers (``trf``, ``lbfgs``) and
profile-likelihood analysis — see :ref:`Gradient-based optimization
<alg-gradient>` — the :ref:`Hamiltonian Monte Carlo (NUTS) <alg-hmc>` and
:ref:`Preconditioned DREAM <alg-p_dream>` samplers, and :ref:`model checking
<alg-check>`. Information criteria and posterior export are covered under
:ref:`Model selection and posterior analysis <model_selection>`.

General implementation features for all algorithms
--------------------------------------------------

All algorithms in PyBNF keep track of a list of parameter sets (a "population"), and over the course of the simulation, submit new parameter sets to run on the simulator. Algorithms periodically output the file ``sorted_params.txt`` containing the best parameter sets found so far, and the corresponding objective function values. 

.. _param-init:

Initialization
^^^^^^^^^^^^^^

The initial population of parameter sets is generated based on the keys specified for each free parameter: ``uniform_var``, ``loguniform_var``, ``normal_var`` or ``lognormal_var``. The value of the parameter in each new random parameter set is drawn from the specified probability distribution. PyBNF ships a broad catalog of prior distribution families beyond these four — see :ref:`priors`.

The ``latin_hypercube`` option for initialization is enabled by default. This option only affects initialization of ``uniform_var``\ s and ``loguniform_var``\ s. When enabled, instead of drawing an independent random value for each starting parameter set, the starting parameter sets are generated with Latin hypercube sampling, which ensures a roughly even distribution of the parameter sets throughout the search space. 

.. _objective:

Objective functions
^^^^^^^^^^^^^^^^^^^

All algorithms use an objective function to evaluate the quality of fit for each parameter set. The objective function is set with the ``objfunc`` key. The following options are available. Note that :math:`y_i` are the experimental data points and :math:`a_i` are the simulated data points. The summation is over all experimental data points.

    * Chi squared (``obj_func = chi_sq``): :math:`f(y, a) =  \sum_i \frac{(y_i - a_i)^2}{2 \sigma_i^2}` , where :math:`\sigma_i` is the standard deviation of point :math:`y_i`, which must be specified in the :ref:`exp file <exp-file>`.
    * Sum of squares (``obj_func = sos``): :math:`f(y, a) =  \sum_i (y_i - a_i)^2`
    * Sum of differences (``obj_func = sod``): :math:`f(y, a) =  \sum_i |y_i - a_i|`
    * Normalized sum of squares (``obj_func = norm_sos``): :math:`f(y, a) =  \sum_i \frac{(y_i - a_i)^2}{y_i^2}`
    * Average-normalized sum of squares (``obj_func = ave_norm_sos``): :math:`f(y, a) =  \sum_i \frac{(y_i - a_i)^2}{\bar{y}^2}`, where :math:`\bar{y}` is the average of the entire data column :math:`y`.
    
If you include any :ref:`constraints <con-file>` in your fit, the constraints add extra terms to the objective function. 

Changing parameter values
^^^^^^^^^^^^^^^^^^^^^^^^^

All algorithms perform changes to parameter values as the fitting proceeds. The way these changes are calculated depends on the type of parameter. 

``loguniform_var``\ s and ``lognormal_var``\ s are moved in logarithmic space (base 10) throughout the entire fitting run. 

``uniform_var``\ s and ``loguniform_var``\ s avoid moving outside the defined initialization range. If a move is attempted that would take the parameter outside the bounds, the parameter value is reflected over the boundary, back within bounds. This feature can be disabled by appending ``U`` to the end of the variable definition (e.g. ``uniform_var = x__FREE 10 30 U``)


.. _alg-de:

Differential Evolution
----------------------

Algorithm
^^^^^^^^^
A population of individuals (points in parameter space) are iteratively evaluated with an objective function.  Parent individuals from the current iteration are selected to form new individuals in the next iteration.  The new individual's parameters are derived by combining parameters from the parents. New individuals are accepted into the population if they have an objective value lower than that of a member of the current population.

Parallelization
^^^^^^^^^^^^^^^
Three versions of differential evolution are available: All run in parallel, but they differ in their level of synchronicity.

Asynchronous differential evolution (``fit_type = ade``) never allows processors to sit idle. One new simulation is started every time a simulation completes. This version is the best choice when a large number of processors are available.

Synchronous differential evolution (``fit_type = de``) consists of discrete iterations. In each iteration, n simulations are run in parallel, but all must complete before moving on to the next iteration. 

Island-based differential evolution [Penas2015]_ is partially asynchronous algorithm. To use this version, set ``fit_type = de`` and set a value greater than 1 for the ``islands`` key. In this version, the current population consists of m islands. Each island is able to move on to the next iteration even if other islands are still in progress. If m is set to the number of available processors, then processors will never sit idle. Note however that this might still underperform compared to the synchronous algorithm run on the same number of processors. 

Implementation details
^^^^^^^^^^^^^^^^^^^^^^

We maintain a list of ``population_size`` current parameter sets, and in each iteration, ``population_size`` new parameter sets are proposed. The method to propose a new parameter set is specified by the config key ``de_strategy``. The default setting ``rand1`` works best for most problems, and runs as follows: We choose 3 random parameter sets p1, p2, and p3 in the current population. For each free parameter P, the new parameter set is assigned the value p1[P] + ``mutation_factor`` * (p2[P]-p3[P]) with probability ``mutation_rate``, or p1[P] with probability 1 - ``mutation_rate``. The new parameter set replaces the parameter set with the same index in the current population if it has a lower objective value. 

With ``de_strategy`` of ``best1`` or ``best2``, we force the above p1 to be the parameter set with the lowest objective value. With ``de_strategy`` of ``all1`` or ``all2``, we force p1 to be the parameter set at the same index we are proposing to replace. The ``best`` strategy results in fast convergence to what is likely only a local optimum. The ``all`` strategy converges more slowly, and prevents the entire population from converging to the same value. However, there is still a risk of each member of the population becoming stuck in its own local minimum. For the ``de_strategy``\ s ending in ``2``, we instead choose a total of 5 parameter sets, p1 through p5, and set the new parameter value as p1[P] + ``mutation_factor`` * (p2[P]-p3[P] + p4[P]-p5[P])

Asynchronous version
""""""""""""""""""""

The asynchronous version of the algorithm is identical to the sychronous algorithm, except that whenever a simulation completes, a new parameter set is immediately proposed based on the current population. Therefore, the random parameter sets p1, p2, and p3 might come from different iteration numbers.

.. _alg-island:

Island-based version
""""""""""""""""""""

In the island-based version of the algorithm [Penas2015]_, the population is divided into ``num_islands`` islands, which each follow the above update procedure independently. Every ``migrate_every`` iterations, a migration step occurs in which ``num_to_migrate`` individuals from each island are transferred randomly to others (according to a random permutation of the islands, keeping the number of individuals on each island constant). The migration step does not require synchronization of the islands; it is performed when the last island reaches the appropriate iteration number, regardless of whether other islands are already further along. 

Applications
^^^^^^^^^^^^
In our experience, differential evolution tends to be a good general-purpose algorithm. 

The asynchronous version has similar advantages to `Particle Swarm`_. 

.. _alg-ss:

Scatter Search
--------------

Algorithm
^^^^^^^^^
Scatter Search [Glover2000]_ functions similarly to differential evolution, but maintains a smaller current population than the number of available processors. In each iteration, every possible pair of individuals are combined to propose a new individual.

Parallelization
^^^^^^^^^^^^^^^
In a scatter search run of population size n, each iteration requires n\*(n-1) independent simulations that can all be run in parallel. Scatter search requires synchronization at the end of each iteration, waiting for all simulations to complete before moving to the next iteration. 

Implementation details
^^^^^^^^^^^^^^^^^^^^^^
The PyBNF implementation follows the outline presented in the introduction of [Penas2017]_ and uses the recombination method described in [Egea2009]_.

We maintain a reference set of ``population_size`` individuals, recommended to be a small number (~ 9-18). Each newly proposed parameter set is based on a "parent" parameter set and a "helper" parameter set, both from the current reference set. In each iteration, we consider all possible parent-helper combinations, for a total of n\*(n-1) parameter sets. The new parameter set depends on the rank of the parent and helper (call them :math:`p_i` and :math:`h_i`) when the reference set is sorted from best to worst. 

Then we apply a series of formulas to choose the next parameter value.

Let :math:`\alpha` = -1 if :math:`h_i>p_i` or 1 if :math:`p_i<h_i`, let :math:`\beta = (|h_i-p_i|-1) / (n-2)`, let :math:`d = \textrm{helper}[P] - \textrm{parent}[P]` for some parameter P. 

Then the in the new parameter set, :math:`P = \textrm{parent}[P] + \textrm{rand\_uniform}(-d * (1 + \alpha * \beta), d * (1 - \alpha * \beta))`

Intuitively what we do here is perturb P on the order of d (which acts as a measure of the variability of P in the population). If the parent is better than the helper, we keep P closer to the parent, and if the helper is better, we shift it closer to the helper. 

The proposed new parameter set is accepted if it achieves a lower objective value than its parent.

If a parent goes ``local_min_limit`` iterations without being replaced by a new parameter set, it is assumed to be stuck in a local minimum, and is replaced with a new random parameter set. The random parameter set is drawn from a "reserve queue", which is initialized at the start of the fitting run to contain ``reserve_size`` Latin hypercube distributed samples. The reserve queue ensures that each time we take a new random parameter set, we are sampling a part of parameter space that we have not sampled previously. 


Applications
^^^^^^^^^^^^
We find scatter search is also a good general-purpose fitting algorithm. It performs especially well on fitting problems that are difficult due to a search space that is high dimensional or contains many local minima. 

.. _alg-pso: 

Particle Swarm
--------------

Algorithm
^^^^^^^^^
In particle swarm optimization, each parameter set is represented by a particle moving through parameter space at some velocity. The acceleration of each particle is set in a way that moves it toward more favorable areas of parameter space: the acceleration has contributions pointing toward both the best parameter set seen so far by the individual particle, and the global best parameter set seen by any particle in the population. 

Parallelization
^^^^^^^^^^^^^^^
Particle swarm optimization in PyBNF is an asynchronous, parallel algorithm. As soon as one simulation completes, that particle can calculate its next parameter set and begin a new simulation. Processors will never remain idle.

Implementation details
^^^^^^^^^^^^^^^^^^^^^^
The PyBNF implementation is based on the description in [Moraes2015]_. Each particle keeps track of its current position, velocity, and the best parameter set it has seen during the run. 

After each simulation completes, the velocity of the particle is updated according to the formula :math:`v_{i+1} = w*v_i + c_1*u_1*(x_i-x_{\textrm{min}}) + c_2*u_2*(x_i-x_{\textrm{globalmin}})`. The constants in the above formula may be set with config keys: *w* is ``particle_weight``, :math:`c_1` is ``cognitive``, and :math:`c_2` is ``social``. :math:`x_i` is the current particle position, :math:`v_i` is the current velocity, :math:`v_{i+1}` is the updated velocity, :math:`x_{\textrm{min}}` is the best parameter set this particle has seen, and :math:`x_{\textrm{globalmin}}` is the best parameter set any particle has seen. :math:`u_1` and :math:`u_2` are uniform random numbers in the range [0,1]. Following the velocity update, the position of the particle is updated by adding its current velocity. 

We apply a special treatment if a ``uniform_var`` or ``loguniform_var`` moves outside of the specified box constraints. As with other algorithms, the particle position is reflected back inside the boundaries. In addition, the component of the velocity corresponding to the parameter that moved out of bounds is set to zero, to prevent the particle from immediately crossing the same boundary again. 

.. _pso-adaptive:

An optional feature (discussed in [Moraes2015]_) allows the particle weight *w* to vary over the course of the simulation. In the original algorithm descirption, *w* was called "inertia weight", but when *w* takes a value less than 1, it can be thought of as friction - a force that decelerates particles regardless of the objective function evaluations. The idea is to reduce *w* (increase friction) over the course of the fitting run, to make the particles come to a stop at a good objective value by the end of the run. 

When using the adaptive friction feature, *w* starts at ``particle_weight``, and approaches ``particle_weight_final`` by the end of the simulation. The value of *w* changes based on how many iterations we deem "unproductive" according to the following criterion: An iteration is unproductive if the global best objective function obj_min changes by less than ``adaptive_abs_tol`` + ``adaptive_rel_tol`` \* obj_min, where ``adaptive_abs_tol`` and ``adaptive_rel_tol`` can be set in the config. Then, we keep track of N, the total number of unproductive iterations so far. At each iteration we set *w* = ``particle_weight`` + (``particle_weight_final`` - ``particle_weight``) \* N / (N + ``adaptive_n_max``). As can be seen in the above formula, the config key ``adaptive_n_max`` sets the number of unproductive iterations it takes to reach halfway between ``particle_weight`` and ``particle_weight_final``.

Applications
^^^^^^^^^^^^
Particle swarm optimization tends to be the best-performing algorithm for problems that benefit from asynchronicity, namely, models for which the runtime per simulation can vary greatly depending on the parameter set. Models simulated by SSA or NFsim often fall into this category. In these cases, synchronous algorithms would cause some cores to remain idle while slow-running simulations complete on other cores, whereas asynchronous algorithms like particle swarm  allow you to use all cores at all times. 

.. _alg-mcmc:

Metropolis-Hastings MCMC
------------------------

Algorithm
^^^^^^^^^
Metropolis-Hastings Markov chain Monte Carlo (MCMC) is a Bayesian sampling method. The free parameters are taken to be random variables with 
unknown probability distributions, and the algorithm samples points in parameter space are sampled with a frequency
proportional to the probability of the parameter set given the data. The result is a probability
distribution over parameter space that expresses the probability of each possible parameter set. With this algorithm, we
obtain not just a point estimate of the best fit, but a means to quantify the uncertainty in each parameter value.

When running this algorithm in PyBNF, it is assumed that the objective function is a *likelihood* function, that is, the objective function value gives the negative log probability of the data given the parameter set. There is a solid mathematical basis for making this assumption when the **chi-squared objective function** is used. It is not recommended to use this algorithm with different objective function, or an objective function that includes qualitative data.

PyBNF outputs additional files containing this probability distribution information. The files in ``Results/Histograms/`` give histograms of the marginal probability distributions for each free parameter. The files ``credible##.txt`` (e.g., ``credible95.txt``) use the marginal histogram for each parameter to calculate a *credible interval* - an interval in which the parameter value is expected to fall with the specified probability (e.g. 95%).  Finally, ``samples.txt`` contains all parameter sets sampled over the course of the fitting run, allowing the user to perform further custom analysis on the sampled probability distribution. 

Parallelization
^^^^^^^^^^^^^^^
Metropolis-Hastings is not an inherently parallel algorithm. In the Markov chain, we need to know the current state before proposing the next one. However, PyBNF supports running several independent Markov chains by specifying the number of chains with the ``population_size`` key. All samples from all parallel chains are pooled to obtain a better estimate of the final posterior probability distribution. 

Note that each chain must independently go through the burn-in period, but after the burn-in, your rate of sampling will be improved proportional to the number of parallel chains in your run. 

Implementation details
^^^^^^^^^^^^^^^^^^^^^^
Our implementation is described in [Kozer2013]_. We start at a random point in parameter space, and make a step of size ``step_size`` to move to a new location in parameter space. We take the value of the objective function to be the negative log probability of the data given the parameter set (the *likelihood* in Bayesian statistics).  We assume a prior distribution based on the parameter definitions in the config file -- a uniform, loguniform, normal, or lognormal distribution, depending on the config key used. Priors are evaluated with ``scipy.stats``; normal and lognormal priors therefore include the distribution normalization constant in reported log-probability values. This constant cancels in Metropolis acceptance ratios. Note: If a uniform or loguniform prior is used, the prior does not affect the result other than to confine the distribution within the specified range. If a normal or lognormal prior is used, the prior does affect the probability of accepting each proposed move, and therefore the choice of prior affects the final sampled probability distribution.

The Bayesian *posterior* distribution -- the probability of the parameters given the data -- is given by the product of the above likelihood and prior. We use the value of the posterior to determine whether to accept the proposed move. 

Moves are accepted according to the Metropolis criterion. If a move increases the value of the posterior, it is always accepted. If it decreases the value of the posterior, it is accepted with probability :math:`e^{- \beta \Delta F}`, where :math:`\Delta F` is the change in the posterior, and :math:`\beta` represents the inverse "temperature" at which the Metropolis sampling occurs. To generate the true posterior distribution, :math:`\beta` should be set to 1. The sampled distribution becomes more broad with smaller :math:`\beta` and more narrow with a larger :math:`\beta`. 


Applications
^^^^^^^^^^^^
Metropolis-Hastings is the simplest method available in PyBNF to generate a probability distribution in parameter space. 

.. _alg-sa:

Simulated Annealing
-------------------

Algorithm
^^^^^^^^^
Simulated annealing is another Markov chain-based algorithm, but our goal is not to find a full probability distribution, just find the optimal parameter set. To do so, we start the Markov chain at a high temperature, where unfavorable moves are accepted frequently, and gradually reduce the temperature over the course of the simulation. The idea is that we will explore parameter space broadly at the start of the fitting run, and become more confined to the optimal region of parameter space as the run proceeds. 

Parallelization
^^^^^^^^^^^^^^^
Simulated annealing is not an inherently parallel algorithm. The trajectory is a Markov chain in which we need to know the current state before proposing the next one. However, PyBNF supports running several independent simulated annealing chains in parallel. By running many chains simulatenously, we have a better chance that one of the chains achieves a good final fit. 

Implementation details
^^^^^^^^^^^^^^^^^^^^^^
The Markov chain is implemented in the same way as described above for the Markov chain Monte Carlo algorithm, incorporating both the objective function value and the prior distribution to calculate the posterior probability density. 

The difference is in the Metropolis criterion for acceptance of a proposed move. Here, a move that decreases the value of the posterior is accepted with probability :math:`e^{- \beta \Delta F}`, where :math:`\beta` decreases over the course of the fitting run. 

Applications
^^^^^^^^^^^^
We have not found any problems for which simulated annealing is better than the other available algorithms, but provide the functionality with the hope that it proves useful for some specific problems. 

.. _alg-pt:

Parallel Tempering
------------------

Algorithm
^^^^^^^^^
Parallel tempering is a more sophisticated version of Markov chain Monte Carlo (MCMC). We run several Markov chains in parallel at different temperatures. At specified iterations during the run, there is an opportunity to exchange replicates between the different temperatures. Only the samples recorded at the lowest temperature count towards our final probability distribution, but the presence of the higher temperature replicates makes it easier to escape local minima and explore the full parameter space. 

When running parallel tempering, PyBNF outputs files containing probability distribution information, the same as with Metropolis-Hastings MCMC.

Like Metropolis-Hastings, it is recommended to only use parallel tempering with the **chi-squared objective function**. 

Parallelization
^^^^^^^^^^^^^^^
The replicates are run in parallel. Synchronization is required at every iteration in which we attempt replica exchange. 

Implementation details
^^^^^^^^^^^^^^^^^^^^^^
The PyBNF implementation is based on the description in [Gupta2018a]_. Markov chains are run by the same method as in Metropolis-Hastings, except that the value of :math:`\beta` in the acceptance probability :math:`e^{- \beta \Delta F}` varies between replicas. 

Every ``exchange_every`` iterations, we attempt replica exchange. We propose moves that consist of swapping two replicas between adjacent temperatures. Moves are accepted with probability :math:`\min (1, e^{\Delta \beta \Delta F})` where :math:`\Delta \beta` is the change in :math:`\beta` = 1/Temperature, and :math:`\Delta F` is the difference in the objective values of the replicas. In other words, moves that transfer a lower-objective replica to a lower temperature (higher :math:`\beta`) are always accepted, and those that transfer a higher-objective replica to a lower temperature are accepted with a Metropolis-like probability based on the extent of objective difference. 

The list of :math:`\beta`\ s used is customizable with the ``beta`` or ``beta_range`` key. The number of replicas per temperature is also customizable. To maintain detailed balance, it is required that each temperature contains the same number of replicas. 


Applications
^^^^^^^^^^^^
Like conventional Metropolis-Hastings MCMC, the goal of parallel tempering is to provide a distribution of possible parameter values rather than a single point estimate. 

Compared to Metropolis-Hastings MCMC, parallel tempering offers a trade-off: Parallel tempering generates fewer samples per unit CPU time (because most of the processors run higher temperature simulations that don't sample the distribution of interest), but traverses parameter space more efficiently, making each sample more valuable. The decision between parallel tempering and Metropolis-Hastings therefore depends on the nature of your parameter space: parallel tempering is expected to perform better when the space is complex, with many local minima that make it challenging to explore. 



.. _alg-am:

Adaptive MCMC
------------------

Algorithm
^^^^^^^^^
Adaptive Markov chain Monte Carlo (MCMC) is a Bayesian sampling method. The Bayesian method is described in further detail under the Metropolis-Hasting definition.
When running this algorithm in PyBNF, it is assumed that the objective function is a likelihood function, that is, the objective function value gives the negative log probability of the data given the parameter set. It is recommended that the Chi-squared or negative binomial or their dynamic versions are used with this algorithm.
The output when using the adaptive MCMC in PyBNF consists of the posterior from each chain and a file containing parameters from all chains saved in the ‘/Results/AMCMC/Run’ folder.

Parallelization
^^^^^^^^^^^^^^^
Like Metropolis-Hastings, Adaptive MCMC is not an inherently parallel algorithm. Like the Metropolis-Hastings algorithm PyBNF supports running several independent Markov chains using the populations_size key.
Note as stated in the Metropolis-Hastings description each chain must independently go through the burn-in period, but after the burn-in, your rate of sampling will be improved proportional to the number of parallel chains in your run.

Implementation details
^^^^^^^^^^^^^^^^^^^^^^
The implementation algorithm can be found in Andrieu and Thoms, Stat Comput 18: 343–373 (2008). The algorithm uses a random walk like MCMC during the training period. While in the training period data is collected to determine the covariance of the posterior. Once the training phase is completed real time on the fly calculation of the diffusivity and covariance is performed for the remaining iterations.
Note: as stated in the Metropolis-Hastings description if a uniform or loguniform prior is used, the prior does not affect the result other than to confine the distribution within the specified range. If a normal or lognormal prior is used, the prior does affect the probability of accepting each proposed move, and therefore the choice of prior affects the final sampled probability distribution.
 

.. _alg-dream:


DREAM
-----

Algorithm
^^^^^^^^^
**D**\ iffe\ **R**\ ential **E**\ volution **A**\ daptive **M**\ etropolis (DREAM), described in [Vrugt2016]_, is an
MCMC approach for estimating the joint probability distribution of a model's free parameters.  DREAM combines features
from traditional Bayesian MCMC (e.g. the Metropolis-Hastings acceptance criterion) and differential evolution (parameter
recombination).  DREAM accelerates convergence of the MCMC and facilitates sampling of multimodal distributions.

PyBNF implements the DREAM(ZS) variant, which draws proposal donors from a growing external archive of past states
rather than from the current chain population.  This allows efficient sampling with as few as 3--5 chains regardless
of the number of parameters.

Parallelization
^^^^^^^^^^^^^^^
DREAM uses parallel MCMC chains evaluated synchronously in generations.  After all chains in a generation are evaluated,
the Metropolis-Hastings acceptance criterion is applied, and new proposals are generated for the next generation.

Proposal mechanisms
^^^^^^^^^^^^^^^^^^^
DREAM supports two proposal mechanisms, selected randomly each generation with probability controlled by the
``snooker_prob`` configuration key (default 0.1 for snooker, remainder for parallel direction).

**Parallel direction update.** The proposal for chain :math:`X` on iteration :math:`i` is:

.. math::

   X_{i+1} = X_i + \zeta + (1 + \lambda)\,\gamma \sum_{j=1}^{\delta}\left(Z_{a_j} - Z_{b_j}\right)

where :math:`Z_{a_j}` and :math:`Z_{b_j}` are pairs of states drawn from the external archive, :math:`\delta` is the
number of chain pairs (``delta``, default 1), :math:`\zeta \sim \mathcal{N}(0, \sigma_\zeta)` is a small normal
perturbation, :math:`\lambda \sim \mathcal{U}(-\lambda_0, \lambda_0)` is a uniform perturbation, and :math:`\gamma` is
the jump rate.  By default, :math:`\gamma` is set adaptively to :math:`2.38/\sqrt{2\delta d'}` where :math:`d'` is the
number of dimensions selected by the crossover mask (this can be overridden by setting ``step_size`` explicitly).
With probability ``gamma_prob`` (default 0.1), a mode jump is performed by setting :math:`\gamma = 1` and updating all
dimensions.

**Snooker update** ([terBraak2008]_). Projects archive points onto the line through the current chain state and a
reference archive point, then jumps along that axis with :math:`\gamma_s \sim \mathcal{U}(1.2, 2.2)`.  The acceptance
criterion includes a Hastings correction factor :math:`\left(\|X_p - Z_c\| / \|X - Z_c\|\right)^{d-1}` to maintain
detailed balance.

Subspace sampling
^^^^^^^^^^^^^^^^^
DREAM incorporates subspace sampling in parameter space via a crossover mechanism.  A set of crossover probabilities
is defined by ``crossover_number`` (default 3), giving values :math:`\{1/n_\mathrm{CR}, 2/n_\mathrm{CR}, \ldots, 1\}`.
One value is drawn each iteration and determines which parameters are updated.  During the first half of burn-in, the
crossover selection probabilities are adapted based on the standardized Euclidean distance traveled by each crossover
value, then frozen for the remainder of the run.

ZS archive
^^^^^^^^^^
The external archive is initialized with ``archive_size`` random draws from the prior (default :math:`10d` where
:math:`d` is the number of parameters).  Every ``archive_thin_rate`` generations (default 10), the current chain states
are appended to the archive.  This constitutes diminishing adaptation, ensuring ergodicity.

Convergence diagnostics
^^^^^^^^^^^^^^^^^^^^^^^
Every 10 iterations, the rank-normalized split-:math:`\hat{R}` statistic is computed and reported
([Vehtari2021]_).  This diagnostic splits each chain in half, rank-normalizes the values to normal scores, and
computes :math:`\hat{R}` on both the ranked and folded (absolute deviation from median) values.  The maximum across
location and scale :math:`\hat{R}` is reported for each parameter.  Values near 1.0 indicate convergence; values above
1.01 suggest the chains have not yet mixed.

If ``rhat_threshold`` is set to a positive value (e.g. 1.05), the algorithm will automatically stop sampling once all
parameters have :math:`\hat{R}` below the threshold (checked after burn-in).

Outlier chain detection
^^^^^^^^^^^^^^^^^^^^^^^
During burn-in, outlier chains are detected every 10 iterations and reset to copies of randomly selected non-outlier
chains.  Two detection methods are available via the ``outlier_method`` key:

- ``iqr`` (default): Flags chains whose mean log-posterior (over the last 50% of history) falls below
  :math:`Q_{25} - 2 \cdot \mathrm{IQR}`.
- ``grubbs``: Applies the Grubbs test at significance level :math:`\alpha = 0.01` to detect a single minimum outlier.


.. _alg-p_dream:


Preconditioned DREAM
--------------------

Algorithm
^^^^^^^^^
Preconditioned DREAM (``fit_type = p_dream``) is :ref:`DREAM(ZS) <alg-dream>` with
its proposals computed in a covariance-whitened parameter space. The
preconditioning transform is estimated from the sampled history (and adapted as
sampling proceeds), so the differential-evolution donor moves are scaled and
rotated to the geometry of a correlated posterior rather than the raw parameter
axes. On posteriors with strong parameter correlations this markedly improves
mixing over plain DREAM; on well-conditioned posteriors it reduces to DREAM(ZS).

Parallelization
^^^^^^^^^^^^^^^
As with :ref:`DREAM <alg-dream>`, the chains advance synchronously: one generation
of all chains is proposed and scored together, and the shared ZS archive and the
preconditioning estimate are updated between generations.

Applications
^^^^^^^^^^^^
Sampling posteriors with strongly correlated parameters, where an axis-aligned or
isotropic proposal mixes slowly.


.. _alg-hmc:


Hamiltonian Monte Carlo (NUTS)
------------------------------

Algorithm
^^^^^^^^^
Hamiltonian Monte Carlo (HMC), via the No-U-Turn Sampler (NUTS) [Hoffman2014]_, is a
**gradient-based** Bayesian sampler. Where the other samplers propose moves blindly and accept or
reject them, HMC follows the gradient of the log-posterior, so it mixes far more efficiently on
correlated, curved, or high-dimensional posteriors. PyBNF drives NUTS through the
`blackjax <https://blackjax-devs.github.io/blackjax/>`_ library and uses HMC as a **reference
sampler** — a gradient-quality yardstick against which the gradient-free samplers (``am`` /
``dream`` / ``p_dream``) are benchmarked on hard geometries.

HMC needs a differentiable log-density, so it runs **only on a closed-form analytical objective**:
a built-in target (``objective = banana, …``) or an inline ``objective = expression`` (see
:ref:`Analytical and user-defined objectives <analytical_objectives>`), including a data-bound
curve fit. It does **not** run on a simulator (BNGL / SBML) model, nor on ``objective = callable``
(a general Python callable is not differentiable); those cases raise a clear error pointing back to
the gradient-free samplers. Select it with ``job_type = hmc`` (requires :ref:`edition <edition>`
``>= 2``).

The same closed-form objective is automatically differentiated with `JAX <https://jax.dev>`_:
the expression compiled for the score path is also lambdified to a JAX function whose gradient
``jax.grad`` produces exactly, so the sampler and the score path can never drift. The full prior
catalog composes — each family contributes a JAX log-density — and constrained parameters
(positive, bounded, log-scaled) are sampled through an unconstraining bijection, so NUTS explores an
unbounded space and never hits a support wall.

Each chain runs an independent NUTS run with window adaptation (dual-averaging step size + mass
matrix) over ``num_warmup`` steps, then keeps ``num_samples`` near-independent draws; the draws are
written in the standard samples format, so the convergence diagnostics, the credible-interval
output, and the :ref:`ArviZ bridge <output_inference_data>` all apply unchanged.

Parallelization
^^^^^^^^^^^^^^^
``population_size`` sets the number of independent chains. Each analytical NUTS chain is a tight
in-process numeric loop, so the chains run as independent blackjax runs rather than through the
per-evaluation dask dispatch the simulator samplers use.

Reliability
^^^^^^^^^^^
Alongside the rank-normalized split-:math:`\hat{R}` and bulk/tail ESS shared with the other samplers
([Vehtari2021]_), HMC reports the **number of divergent transitions** per chain — a NUTS-specific
signal that the integrator could not traverse a region (too sharp a curvature for the tuned step
size). A nonzero divergence count, like a high :math:`\hat{R}`, means the draws are not yet
trustworthy: tighten ``target_accept`` (e.g. to 0.95, which shrinks the step size) and re-run.

Installation
^^^^^^^^^^^^
HMC needs the optional gradient stack (JAX + blackjax)::

    pip install pybnf[jax]

The gradient-free samplers need no extra. Requires Python with JAX support.


.. _alg-sim:

Simplex
-------

Algorithm
^^^^^^^^^
Simplex is a local search algorithm that operates solely on objective evaluations at single points (i.e. it does not require calculation of gradients). The algorithm maintains a set on N+1 points in N-dimensional parameter space, which are thought of as defining an N-dimensional solid called a *simplex*. Individual points may be reflected through the lower-dimensional solid defined by the other N points, to obtain a local improvement in objective function value. The simplex algorithm has been nicknamed the "amoeba" algorithm because the simplex crawls through parameter space similar to an amoeba, extending protrusions in favorable directions.

Parallelization
^^^^^^^^^^^^^^^
The PyBNF Simplex implementation is parallel and synchronous. Synchronization is required at the end of every iteration. Parallelization is achieved by simultaneously evaluating a subset of the N+1 points in the simplex. Therefore, this parallelization can take advantage of at most N+1 processors, where N is the number of free parameters. 


Implementation details
^^^^^^^^^^^^^^^^^^^^^^
PyBNF implements the parallelized Simplex algorithm described in [Lee2007]_. 

The initial simplex consists of N+1 points chosen deterministically based on the specified step size (set with the ``simplex_step`` and ``simplex_log_step`` keys, or for individual parameters with the ``var`` and ``log_var`` keys). One point of the simplex is the specified starting point for the search. The other N points are obtained by adding the step size to one parameter, and leaving the other N-1 parameters at the starting values.


.. figure:: simplex.png
   :width: 200px
   :align: center
   :figclass: align-center

   Illustration of the simplex algorithm, modifying point P on a 3-point simplex in 2 dimensions

Each iteration, we operate on the k worst points in the simplex, where k is the number of available processors (``parallel_count``). For each point P, we  consider the hyperplane defined by the other N points in the simplex (blue line). Let d be the distance from P to the hyperplane. We evaluate point P\ :sub:`1` obtained by reflecting P through the hyperplane, to a distance of d \* ``simplex_reflect`` on the other side. Depending on the resulting objective value, we try another point in the second phase of the iteration. Three cases are possible.

1) The new point is better than the current global minimum: We try a second point continuing in the same direction for a distance of d \* ``simplex_expansion`` away from the hyperplane (P\ :sub:`2,1`).
2) The new point is worse than the global minimum, but better than the next worst point in the simplex: We don't try a second point.
3) The new point is worse than the next worst point in the simplex: We try a second point moving closer to the hyperplane. If P was better than P\ :sub:`1`, we try a point a distance of d \* ``simplex_contraction`` from the hyperplane in the direction of P (P\ :sub:`2,3a`). If P\ :sub:`1` was better than P, we instead try the same distance from the hyperplane in the direction of P\ :sub:`1` (P\ :sub:`2,3b`).

In all cases, P in the simplex is set to the best choice among P, P\ :sub:`1`, or whichever second point we tried.

If in a given iteration, all k points resulted in Case 3 and did not update to P\ :sub:`2,3a` or P\ :sub:`2,3b`, the iteration did not effectively change the state of the simplex. Then, we contract the simplex towards the best point: We set each point P to ``simplex_contract`` \* P0 + (1 - ``simplex_contract``) \* P, where P0 is the best point in the simplex.

Applications
^^^^^^^^^^^^
Local optimization with the simplex algorithm is useful for improving on an already known good solution. In PyBNF, the most common application is to apply the simplex algorithm to the best-fit result obtained from one of the other algorithms. You can automatically refine your final result with the simplex algorithm by setting the ``refine`` key to 1, and setting simplex config keys in addition to the config for your main algorithm.

It is also possible to run the Simplex algorithm on its own, using a custom starting point. In this case, you should use the ``var`` and ``log_var`` keys to specify your known starting point.


.. _refinement:

Refinement (choosing a local optimizer)
---------------------------------------
PyBNF offers three derivative-free, black-box local optimizers for the post-fit
polish step (``refine = 1``): `Simplex`_ (Nelder–Mead), `Powell`_, and `CMA-ES`_.
Set ``refine_method`` to ``sim`` (the default, backward-compatible), ``powell``,
or ``cmaes`` to choose; only the chosen optimizer's config keys are read. Each
also runs standalone as its own ``fit_type`` (``sim`` / ``powell`` / ``cmaes``),
started from a single point given with the ``var`` / ``logvar`` keys. CMA-ES
additionally runs standalone as a *global* optimizer over a bounded ``uniform_var``
/ ``loguniform_var`` box (see `CMA-ES`_).

All three need only objective values (no gradients), which suits PyBNF's
black-box, often-noisy simulator objectives. As rough guidance:

* **Simplex** — the long-standing default; a robust, low-overhead amoeba search.
* **Powell** — conjugate-direction line minimization; often converges faster than
  Simplex on smooth, well-behaved objectives.
* **CMA-ES** — adapts a full covariance model of the search distribution, so it is
  the most robust of the three on ill-conditioned or rotated (correlated)
  objectives, at the cost of more evaluations per step. Being population-based, it
  also parallelizes across the whole generation.


.. _alg-powell:

Powell
------

Algorithm
^^^^^^^^^
Powell's method is a derivative-free local optimizer that minimizes the objective
by a sequence of one-dimensional line minimizations along a set of search
directions (initially the coordinate axes). After each cycle of line searches it
may replace one direction with the net direction of the cycle's progress, building
up *conjugate* directions that give quadratic convergence on a quadratic bowl
([Powell1964]_; see also [NumericalRecipes]_). PyBNF performs each line
minimization by **bracketing** the minimum (geometric expansion from
``±powell_step``) and refining it with **Brent's method** (parabolic interpolation
with a golden-section fallback) to ``powell_line_tol``. This follows long, curved,
non-quadratic valleys and adapts its step length, where a single fixed-step
parabola would stall. The line search is confined to the parameter box, so when
refining a bounded fit a minimum that lies past a bound is found on the boundary.

Parallelization
^^^^^^^^^^^^^^^
Powell is fundamentally serial — each line-search evaluation depends on the
previous, so the bracketing + Brent search runs one objective evaluation at a
time. (CMA-ES is the derivative-free optimizer that evaluates a whole generation
in parallel.) The search runs in the parameter sampling space (log-scaled for log
parameters), so a step is geometric for a log parameter.

Applications
^^^^^^^^^^^^
Use Powell to refine an already-good solution (``refine = 1`` with
``refine_method = powell``), or standalone (``fit_type = powell``) from a starting
point given with the ``var`` / ``logvar`` keys. Tune ``powell_step`` (initial
bracketing step), ``powell_line_tol`` (1-D line-minimum precision),
``powell_stop_tol`` (per-cycle convergence), and ``powell_max_iterations`` (cycle
budget; defaults to ``max_iterations``).


.. _alg-cmaes:

CMA-ES
------

Algorithm
^^^^^^^^^
The Covariance Matrix Adaptation Evolution Strategy ([Hansen2001]_) is a
derivative-free, population-based optimizer. Each generation it samples
``population_size`` candidate parameter sets from a multivariate normal
distribution and, from the best of them, updates the distribution's mean, overall
step size, and full covariance matrix. Adapting the covariance lets the search
stretch and rotate to match the local objective geometry, which makes CMA-ES
notably robust on ill-conditioned and correlated problems.

Parallelization
^^^^^^^^^^^^^^^
Synchronous and population-based: every candidate in a generation is evaluated in
parallel, then the distribution is updated. Like the other start-point optimizers
it searches in the parameter sampling space (log-scaled for log parameters).

Applications
^^^^^^^^^^^^
Use CMA-ES to refine a result (``refine = 1`` with ``refine_method = cmaes``) or
standalone (``fit_type = cmaes``). Standalone, it accepts either of two starts:

* a single starting **point** given with the ``var`` / ``logvar`` keys (local
  search, like Simplex and Powell); or
* a bounded **box** given with the ``uniform_var`` / ``loguniform_var`` keys —
  its *global-start* mode. CMA-ES begins at the box center and seeds its
  covariance with the per-coordinate box widths, so the first generation spans the
  whole box; candidates are reflected back into the box. Because covariance
  adaptation makes CMA-ES one of the strongest general-purpose black-box global
  optimizers, this is the recommended way to run it as a primary search (the
  ergonomics of ``de`` / ``pso``, with covariance adaptation). Use ``var`` /
  ``logvar`` only with a single value per parameter, and ``uniform_var`` /
  ``loguniform_var`` only with bounds — do not mix the two.

The population size is ``population_size`` (at least 4) and the generation budget
is ``max_iterations``. ``cmaes_sigma0`` sets the initial step size — in point-start
mode an absolute step in the sampling space, in box mode a fraction of each box
width — and ``cmaes_stop_tol`` the convergence threshold on the search
distribution's spread. As a refiner the start is always the injected best fit.


.. _alg-gradient:


Gradient-based optimization
---------------------------

For edition-2 fits of ODE-network models, PyBNF offers two gradient-based
optimizers driven by exact forward parameter sensitivities:

* ``fit_type = trf`` — a trust-region least-squares method (Trust-Region
  Reflective), for objectives that are a sum of squared residuals; and
* ``fit_type = lbfgs`` — a quasi-Newton method (L-BFGS-B) that minimizes a scalar
  objective from its analytic gradient.

Both converge far faster than the metaheuristics near a good fit, and the same
sensitivity machinery drives profile-likelihood identifiability analysis
(``fit_type = profile_likelihood``). These methods, the noise families and
constraints they support, and the capability gate that decides when a gradient fit
is possible are documented in full on the :ref:`gradient-based fitting
<gradient_fitting>` page.


.. _alg-check:


Model checking
--------------

``fit_type = check`` (equivalently ``job_type = check``) is a first-class checking
method rather than a search: it evaluates the objective value and constraint
satisfaction for a given set of parameter values without exploring parameter
space. Use it to verify that a model satisfies a set of qualitative properties
(:ref:`constraints <con-file>` / BPSL) or to score a specific parameter set — a
quick way to test a candidate fit or validate a model against its specification.
See the `model-checking tutorial lesson
<https://github.com/lanl/PyBNF/tree/main/examples/tutorial/46_model_checking>`__
for a worked example.


.. _model_selection:


Model selection and posterior analysis
--------------------------------------

Beyond a single best fit, PyBNF supports comparing models and characterizing
posteriors.

**Information criteria.** From a fit's maximized log-likelihood PyBNF computes the
Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC), and
the small-sample-corrected AIC (AICc), using the *full* normalized log-likelihood
(including the noise-model normalizer):

.. math::

   \mathrm{AIC} = 2k - 2\ln L, \qquad
   \mathrm{BIC} = k\ln n - 2\ln L, \qquad
   \mathrm{AICc} = \mathrm{AIC} + \frac{2k(k+1)}{n-k-1},

where :math:`k` is the number of free parameters and :math:`n` the number of data
points. Ranking competing models by AIC is the basis of the model-selection
tutorial lesson.

**ArviZ / LOO / WAIC.** A completed MCMC run can be exported as an ArviZ
``InferenceData`` object by setting
:ref:`output_inference_data <output_inference_data>` ``= 1`` (needs the optional
``arviz`` extra, ``pip install pybnf[arviz]``). The posterior then plugs into the
ArviZ / bayesplot ecosystem for trace, rank, forest, and pair plots and for
convergence diagnostics (R-hat, bulk/tail ESS). When the run also records the
per-observation log-likelihood, the exported object carries a ``log_likelihood``
group, so leave-one-out cross-validation (``az.loo``) and the widely-applicable
information criterion (``az.waic``) can be computed directly.


.. [Egea2009] Egea, J. A.; Balsa-Canto, E.; García, M.-S. G.; Banga, J. R. Dynamic Optimization of Nonlinear Processes with an Enhanced Scatter Search Method. Ind. Eng. Chem. Res. 2009, 48 (9), 4388–4401.
.. [Glover2000] Glover, F.; Laguna, M.; Martí, R. Fundamentals of Scatter Search and Path Relinking. Control Cybern. 2000, 29 (3), 652–684.
.. [Hoffman2014] Hoffman, M. D.; Gelman, A. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res. 2014, 15 (1), 1593–1623.
.. [Hansen2001] Hansen, N.; Ostermeier, A. Completely Derandomized Self-Adaptation in Evolution Strategies. Evol. Comput. 2001, 9 (2), 159–195.
.. [Gupta2018a] Gupta, S.; Hainsworth, L.; Hogg, J. S.; Lee, R. E. C.; Faeder, J. R. Evaluation of Parallel Tempering to Accelerate Bayesian Parameter Estimation in Systems Biology. 2018 26th Euromicro International Conference on Parallel, Distributed and Network-based Processing (PDP) 2018, 690–697.
.. [Kozer2013] Kozer, N.; Barua, D.; Orchard, S.; Nice, E. C.; Burgess, A. W.; Hlavacek, W. S.; Clayton, A. H. A. Exploring Higher-Order EGFR Oligomerisation and Phosphorylation—a Combined Experimental and Theoretical Approach. Mol. BioSyst. Mol. BioSyst 2013, 9 (9), 1849–1863.
.. [Lee2007] Lee, D.; Wiswall, M. A Parallel Implementation of the Simplex Function Minimization Routine. Comput. Econ. 2007, 30 (2), 171–187.
.. [Moraes2015] Moraes, A. O. S.; Mitre, J. F.; Lage, P. L. C.; Secchi, A. R. A Robust Parallel Algorithm of the Particle Swarm Optimization Method for Large Dimensional Engineering Problems. Appl. Math. Model. 2015, 39 (14), 4223–4241.
.. [Penas2015] Penas, D. R.; González, P.; Egea, J. A.; Banga, J. R.; Doallo, R. Parallel Metaheuristics in Computational Biology: An Asynchronous Cooperative Enhanced Scatter Search Method. Procedia Comput. Sci. 2015, 51 (1), 630–639.
.. [Penas2017] Penas, D. R.; González, P.; Egea, J. A.; Doallo, R.; Banga, J. R. Parameter Estimation in Large-Scale Systems Biology Models: A Parallel and Self-Adaptive Cooperative Strategy. BMC Bioinformatics 2017, 18 (1), 52.
.. [Powell1964] Powell, M. J. D. An Efficient Method for Finding the Minimum of a Function of Several Variables without Calculating Derivatives. Comput. J. 1964, 7 (2), 155–162.
.. [NumericalRecipes] Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press, 2007 (§10.7, Powell's method).
.. [terBraak2008] ter Braak, C. J. F.; Vrugt, J. A. Differential Evolution Markov Chain with Snooker Updater and Fewer Chains. Stat. Comput. 2008, 18 (4), 435–446.
.. [Vehtari2021] Vehtari, A.; Gelman, A.; Simpson, D.; Carpenter, B.; Bürkner, P.-C. Rank-Normalization, Folding, and Localization: An Improved R-hat for Assessing Convergence of MCMC. Bayesian Anal. 2021, 16 (2), 667–718.
.. [Vrugt2016] Vrugt, J. A. Markov Chain Monte Carlo Simulation Using the DREAM Software Package: Theory, Concepts, and MATLAB Implementation. Environ. Model. Softw. 2016, 75, 273–316.