Module multichain

source code

Implements several extended-ensemble Monte Carlo sampling algorithms.

Here is a short example which shows how to sample from a PDF using the replica exchange with non-equilibrium switches (RENS) method. It draws 5000 samples from a 1D normal distribution using the RENS algorithm working on three Markov chains being generated by the HMC algorithm:

>>> import numpy
>>> from numpy import sqrt
>>> from csb.io.plots import Chart
>>> from csb.statistics.pdf import Normal
>>> from csb.statistics.samplers import State
>>> from csb.statistics.samplers.mc.multichain import ThermostattedMDRENSSwapParameterInfo
>>> from csb.statistics.samplers.mc.multichain import ThermostattedMDRENS, AlternatingAdjacentSwapScheme
>>> from csb.statistics.samplers.mc.singlechain import HMCSampler

>>> # Pick some initial state for the different Markov chains:
>>> initial_state = State(numpy.array([1.]))

>>> # Set standard deviations:
>>> std_devs = [1./sqrt(5), 1. / sqrt(3), 1.]

>>> # Set HMC timesteps and trajectory length:
>>> hmc_timesteps = [0.6, 0.7, 0.6]
>>> hmc_trajectory_length = 20
>>> hmc_gradients = [lambda q, t: 1 / (std_dev ** 2) * q for std_dev in std_devs]

>>> # Set parameters for the thermostatted RENS algorithm:
>>> rens_trajectory_length = 30
>>> rens_timesteps = [0.3, 0.5]

>>> # Set interpolation gradients as a function of the work parameter l:
>>> rens_gradients = [lambda q, l, i=i: (l / (std_devs[i + 1] ** 2) + (1 - l) / (std_devs[i] ** 2)) * q 
                      for i in range(len(std_devs)-1)]

>>> # Initialize HMC samplers:
>>> samplers = [HMCSampler(Normal(sigma=std_devs[i]), initial_state, hmc_gradients[i], hmc_timesteps[i],
                hmc_trajectory_length) for i in range(len(std_devs))]

>>> # Create swap parameter objects:
params = [ThermostattedMDRENSSwapParameterInfo(samplers[0], samplers[1], rens_timesteps[0],
          rens_trajectory_length, rens_gradients[0]),
          ThermostattedMDRENSSwapParameterInfo(samplers[1], samplers[2], rens_timesteps[1],
          rens_trajectory_length, rens_gradients[1])]

>>> # Initialize thermostatted RENS algorithm:
>>> algorithm = ThermostattedMDRENS(samplers, params)

>>> # Initialize swapping scheme:
>>> swapper = AlternatingAdjacentSwapScheme(algorithm)

>>> # Initialize empty list which will store the samples:
>>> states = []
>>> for i in range(5000):
        if i % 5 == 0:
            swapper.swap_all()
        states.append(algorithm.sample())

>>> # Print acceptance rates:
>>> print('HMC acceptance rates:', [s.acceptance_rate for s in samplers])
>>> print('swap acceptance rates:', algorithm.acceptance_rates)

>>> # Create and plot histogram for first sampler and numpy.random.normal reference:
>>> chart = Chart()
>>> rawstates = [state[0].position[0] for state in states]
>>> chart.plot.hist([numpy.random.normal(size=5000, scale=std_devs[0]), rawstates], bins=30, normed=True)
>>> chart.plot.legend(['numpy.random.normal', 'RENS + HMC'])
>>> chart.show()

For ReplicaExchangeMC (RE), the procedure is easier because apart from the two sampler instances the corresponding RESwapParameterInfo objects take no arguments.

Every replica exchange algorithm in this module (ReplicaExchangeMC, MDRENS, ThermostattedMDRENS) is used in a similar way. A simulation is always initialized with a list of samplers (instances of classes derived from AbstractSingleChainMC) and a list of AbstractSwapParameterInfo objects suited for the algorithm under consideration. Every AbstractSwapParameterInfo object holds all the information needed to perform a swap between two samplers. The usual scheme is to swap only adjacent replicae in a scheme:

   1 <--> 2, 3 <--> 4, ...
   2 <--> 3, 4 <--> 5, ...
   1 <--> 2, 3 <--> 4, ...

This swapping scheme is implemented in the AlternatingAdjacentSwapScheme class, but different schemes can be easily implemented by deriving from AbstractSwapScheme. Then the simulation is run by looping over the number of samples to be drawn and calling the AbstractExchangeMC.sample method of the algorithm. By calling the AbstractSwapScheme.swap_all method of the specific AbstractSwapScheme implementation, all swaps defined in the list of AbstractSwapParameterInfo objects are performed according to the swapping scheme. The AbstractSwapScheme.swap_all method may be called for example after sampling intervals of a fixed length or randomly.