Package csb :: Package statistics :: Package samplers :: Package mc :: Module multichain |
[frames] | no frames]

# 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.

Classes
| AbstractEnsembleMC Abstract class for Monte Carlo sampling algorithms simulating several ensembles.
| AbstractExchangeMC Abstract class for Monte Carlo sampling algorithms employing some replica exchange method.
| AbstractRENS Abstract Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski 2009) class.
| AbstractRENSSwapParameterInfo Holds parameters for a RENS swap.
| AbstractStepRENS Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski 2009) with stepwise trajectories as described in Nilmeier et al., "Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation", PNAS 2011.
| AbstractSwapCommunicator Holds all the information which needs to be communicated between distinct swap substeps.
| AbstractSwapParameterInfo Subclass instances hold all parameters necessary for performing a swap between two given samplers.
| AbstractSwapScheme Provides the interface for classes defining schemes according to which swaps in Replica Exchange-like simulations are performed.
| AlternatingAdjacentSwapScheme Provides a swapping scheme in which tries exchanges between neighbours only following the scheme 1 <-> 2, 3 <-> 4, ...
| Gradient
| HMCStepRENS Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski 2009) with stepwise trajectories as described in Nilmeier et al., "Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation", PNAS 2011.
| HMCStepRENSSwapParameterInfo Holds all required information for performing HMCStepRENS swaps.
| InterpolationFactory Produces interpolations for functions changed during non-equilibrium trajectories.
| MDRENS Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski 2009) with Molecular Dynamics (MD) trajectories.
| MDRENSSwapParameterInfo Holds parameters for a MDRENS swap.
| RENSSwapCommunicator Holds all the information which needs to be communicated between distinct RENS swap substeps.
| RENSTrajInfo Holds information necessary for calculating a RENS trajectory.
| RESwapCommunicator Holds all the information which needs to be communicated between distinct RE swap substeps.
| RESwapParameterInfo Holds parameters for a standard Replica Exchange swap.
| ReplicaExchangeMC Replica Exchange (RE, Swendsen & Yang 1986) implementation.
| ReplicaHistory Replica history object, works with both RE and RENS for the AlternatingAdjacentSwapScheme.
| SingleSwapStatistics Tracks swap statistics of a single sampler pair.
| SwapStatistics Tracks swap statistics for an AbstractExchangeMC subclass instance.
| ThermostattedMDRENS Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski, 2009) with Andersen-thermostatted Molecular Dynamics (MD) trajectories.
| ThermostattedMDRENSSwapParameterInfo