1 """
2 Implements several extended-ensemble Monte Carlo sampling algorithms.
3
4 Here is a short example which shows how to sample from a PDF using the replica
5 exchange with non-equilibrium switches (RENS) method. It draws 5000 samples from
6 a 1D normal distribution using the RENS algorithm working on three Markov chains
7 being generated by the HMC algorithm:
8
9
10 >>> import numpy
11 >>> from numpy import sqrt
12 >>> from csb.io.plots import Chart
13 >>> from csb.statistics.pdf import Normal
14 >>> from csb.statistics.samplers import State
15 >>> from csb.statistics.samplers.mc.multichain import ThermostattedMDRENSSwapParameterInfo
16 >>> from csb.statistics.samplers.mc.multichain import ThermostattedMDRENS, AlternatingAdjacentSwapScheme
17 >>> from csb.statistics.samplers.mc.singlechain import HMCSampler
18
19 >>> # Pick some initial state for the different Markov chains:
20 >>> initial_state = State(numpy.array([1.]))
21
22 >>> # Set standard deviations:
23 >>> std_devs = [1./sqrt(5), 1. / sqrt(3), 1.]
24
25 >>> # Set HMC timesteps and trajectory length:
26 >>> hmc_timesteps = [0.6, 0.7, 0.6]
27 >>> hmc_trajectory_length = 20
28 >>> hmc_gradients = [lambda q, t: 1 / (std_dev ** 2) * q for std_dev in std_devs]
29
30 >>> # Set parameters for the thermostatted RENS algorithm:
31 >>> rens_trajectory_length = 30
32 >>> rens_timesteps = [0.3, 0.5]
33
34 >>> # Set interpolation gradients as a function of the work parameter l:
35 >>> rens_gradients = [lambda q, l, i=i: (l / (std_devs[i + 1] ** 2) + (1 - l) / (std_devs[i] ** 2)) * q
36 for i in range(len(std_devs)-1)]
37
38 >>> # Initialize HMC samplers:
39 >>> samplers = [HMCSampler(Normal(sigma=std_devs[i]), initial_state, hmc_gradients[i], hmc_timesteps[i],
40 hmc_trajectory_length) for i in range(len(std_devs))]
41
42 >>> # Create swap parameter objects:
43 params = [ThermostattedMDRENSSwapParameterInfo(samplers[0], samplers[1], rens_timesteps[0],
44 rens_trajectory_length, rens_gradients[0]),
45 ThermostattedMDRENSSwapParameterInfo(samplers[1], samplers[2], rens_timesteps[1],
46 rens_trajectory_length, rens_gradients[1])]
47
48 >>> # Initialize thermostatted RENS algorithm:
49 >>> algorithm = ThermostattedMDRENS(samplers, params)
50
51 >>> # Initialize swapping scheme:
52 >>> swapper = AlternatingAdjacentSwapScheme(algorithm)
53
54 >>> # Initialize empty list which will store the samples:
55 >>> states = []
56 >>> for i in range(5000):
57 if i % 5 == 0:
58 swapper.swap_all()
59 states.append(algorithm.sample())
60
61 >>> # Print acceptance rates:
62 >>> print('HMC acceptance rates:', [s.acceptance_rate for s in samplers])
63 >>> print('swap acceptance rates:', algorithm.acceptance_rates)
64
65 >>> # Create and plot histogram for first sampler and numpy.random.normal reference:
66 >>> chart = Chart()
67 >>> rawstates = [state[0].position[0] for state in states]
68 >>> chart.plot.hist([numpy.random.normal(size=5000, scale=std_devs[0]), rawstates], bins=30, normed=True)
69 >>> chart.plot.legend(['numpy.random.normal', 'RENS + HMC'])
70 >>> chart.show()
71
72
73 For L{ReplicaExchangeMC} (RE), the procedure is easier because apart from the
74 two sampler instances the corresponding L{RESwapParameterInfo} objects take
75 no arguments.
76
77 Every replica exchange algorithm in this module (L{ReplicaExchangeMC}, L{MDRENS},
78 L{ThermostattedMDRENS}) is used in a similar way. A simulation is always
79 initialized with a list of samplers (instances of classes derived from
80 L{AbstractSingleChainMC}) and a list of L{AbstractSwapParameterInfo} objects
81 suited for the algorithm under consideration. Every L{AbstractSwapParameterInfo}
82 object holds all the information needed to perform a swap between two samplers.
83 The usual scheme is to swap only adjacent replicae in a scheme::
84
85 1 <--> 2, 3 <--> 4, ...
86 2 <--> 3, 4 <--> 5, ...
87 1 <--> 2, 3 <--> 4, ...
88
89 This swapping scheme is implemented in the L{AlternatingAdjacentSwapScheme} class,
90 but different schemes can be easily implemented by deriving from L{AbstractSwapScheme}.
91 Then the simulation is run by looping over the number of samples to be drawn
92 and calling the L{AbstractExchangeMC.sample} method of the algorithm. By calling
93 the L{AbstractSwapScheme.swap_all} method of the specific L{AbstractSwapScheme}
94 implementation, all swaps defined in the list of L{AbstractSwapParameterInfo}
95 objects are performed according to the swapping scheme. The
96 L{AbstractSwapScheme.swap_all} method may be called for example after sampling
97 intervals of a fixed length or randomly.
98 """
99
100 import numpy
101
102 import csb.numeric
103
104 from abc import ABCMeta, abstractmethod
105
106 from csb.statistics.samplers import EnsembleState
107 from csb.statistics.samplers.mc import AbstractMC, Trajectory, MCCollection, augment_state
108 from csb.statistics.samplers.mc.propagators import MDPropagator, ThermostattedMDPropagator
109 from csb.statistics.samplers.mc.neqsteppropagator import NonequilibriumStepPropagator
110 from csb.statistics.samplers.mc.neqsteppropagator import Protocol, Step, ReducedHamiltonian
111 from csb.statistics.samplers.mc.neqsteppropagator import ReducedHamiltonianPerturbation
112 from csb.statistics.samplers.mc.neqsteppropagator import HMCPropagation, HMCPropagationParam
113 from csb.statistics.samplers.mc.neqsteppropagator import HamiltonianSysInfo, NonequilibriumTrajectory
114 from csb.numeric.integrators import AbstractGradient, FastLeapFrog
118 """
119 Abstract class for Monte Carlo sampling algorithms simulating several ensembles.
120
121 @param samplers: samplers which sample from their respective equilibrium distributions
122 @type samplers: list of L{AbstractSingleChainMC}
123 """
124
125 __metaclass__ = ABCMeta
126
133
145
146 @property
148 """
149 Total ensemble energy.
150 """
151 return sum([x.energy for x in self._samplers])
152
155 """
156 Abstract class for Monte Carlo sampling algorithms employing some replica exchange method.
157
158 @param samplers: samplers which sample from their respective equilibrium distributions
159 @type samplers: list of L{AbstractSingleChainMC}
160
161 @param param_infos: list of ParameterInfo instances providing information needed
162 for performing swaps
163 @type param_infos: list of L{AbstractSwapParameterInfo}
164 """
165
166 __metaclass__ = ABCMeta
167
168 - def __init__(self, samplers, param_infos):
177
181
182 - def swap(self, index):
183 """
184 Perform swap between sampler pair described by param_infos[index]
185 and return outcome (true = accepted, false = rejected).
186
187 @param index: index of swap pair in param_infos
188 @type index: int
189
190 @rtype: boolean
191 """
192 param_info = self._param_infos[index]
193 swapcom = self._propose_swap(param_info)
194 swapcom = self._calc_pacc_swap(swapcom)
195 result = self._accept_swap(swapcom)
196
197 self.state = EnsembleState([x.state for x in self._samplers])
198
199 self.statistics.stats[index].update(result)
200
201 return result
202
203 @abstractmethod
205 """
206 Calculate proposal states for a swap between two samplers.
207
208 @param param_info: ParameterInfo instance holding swap parameters
209 @type param_info: L{AbstractSwapParameterInfo}
210
211 @rtype: L{AbstractSwapCommunicator}
212 """
213 pass
214
215 @abstractmethod
217 """
218 Calculate probability to accept a swap given initial and proposal states.
219
220 @param swapcom: SwapCommunicator instance holding information to be communicated
221 between distinct swap substeps
222 @type swapcom: L{AbstractSwapCommunicator}
223
224 @rtype: L{AbstractSwapCommunicator}
225 """
226 pass
227
249
250 @property
258
259 @property
261 """
262 List of SwapParameterInfo instances holding all necessary parameters.
263
264 @rtype: list of L{AbstractSwapParameterInfo}
265 """
266 return self._param_infos
267
268 @property
270 return self._statistics
271
273 """
274 Update statistics of a given swap process.
275
276 @param index: position of swap statistics to be updated
277 @type index: int
278
279 @param accepted: outcome of the swap
280 @type accepted: boolean
281 """
282
283 self._stats[index][0] += 1
284 self._stats[index][1] += int(accepted)
285
288 """
289 Subclass instances hold all parameters necessary for performing a swap
290 between two given samplers.
291 """
292
293 __metaclass__ = ABCMeta
294
295 - def __init__(self, sampler1, sampler2):
296 """
297 @param sampler1: First sampler
298 @type sampler1: L{AbstractSingleChainMC}
299
300 @param sampler2: Second sampler
301 @type sampler2: L{AbstractSingleChainMC}
302 """
303
304 self._sampler1 = sampler1
305 self._sampler2 = sampler2
306
307 @property
309 return self._sampler1
310
311 @property
313 return self._sampler2
314
317 """
318 Holds all the information which needs to be communicated between
319 distinct swap substeps.
320
321 @param param_info: ParameterInfo instance holding swap parameters
322 @type param_info: L{AbstractSwapParameterInfo}
323
324 @param traj12: Forward trajectory
325 @type traj12: L{Trajectory}
326
327 @param traj21: Reverse trajectory
328 @type traj21: L{Trajectory}
329 """
330
331 __metaclass__ = ABCMeta
332
333 - def __init__(self, param_info, traj12, traj21):
345
346 @property
348 return self._sampler1
349
350 @property
352 return self._sampler2
353
354 @property
357
358 @property
361
362 @property
364 return self._acceptance_probability
365 @acceptance_probability.setter
367 self._acceptance_probability = value
368
369 @property
371 return self._accepted
372 @accepted.setter
374 self._accepted = value
375
376 @property
378 return self._param_info
379
382 """
383 Replica Exchange (RE, Swendsen & Yang 1986) implementation.
384 """
385
392
413
416 """
417 Holds parameters for a standard Replica Exchange swap.
418 """
419 pass
420
423 """
424 Holds all the information which needs to be communicated between distinct
425 RE swap substeps.
426
427 See L{AbstractSwapCommunicator} for constructor signature.
428 """
429 pass
430
433 """
434 Abstract Replica Exchange with Nonequilibrium Switches
435 (RENS, Ballard & Jarzynski 2009) class.
436 Subclasses implement various ways of generating trajectories
437 (deterministic or stochastic).
438 """
439
440 __metaclass__ = ABCMeta
441
454
456 """
457 Sets the protocol lambda(t) to either the forward or the reverse protocol.
458
459 @param traj_info: TrajectoryInfo object holding information neccessary to
460 calculate the rens trajectories.
461 @type traj_info: L{RENSTrajInfo}
462 """
463
464 if traj_info.direction == "fw":
465 return traj_info.param_info.protocol
466 else:
467 return lambda t, tau: traj_info.param_info.protocol(tau - t, tau)
468
469 return protocol
470
472 """
473 Determine the initial temperature of a RENS trajectory.
474
475 @param traj_info: TrajectoryInfo object holding information neccessary to
476 calculate the RENS trajectory.
477 @type traj_info: L{RENSTrajInfo}
478 """
479
480 if traj_info.direction == "fw":
481 return traj_info.param_info.sampler1.temperature
482 else:
483 return traj_info.param_info.sampler2.temperature
484
485 @abstractmethod
487 """
488 Calculates the works expended during the nonequilibrium
489 trajectories.
490
491 @param swapcom: Swap communicator object holding all the
492 neccessary information.
493 @type swapcom: L{RENSSwapCommunicator}
494
495 @return: The expended during the forward and the backward
496 trajectory.
497 @rtype: 2-tuple of floats
498 """
499
500 pass
501
508
509 @abstractmethod
511 """
512 Factory method which produces the propagator object used to calculate
513 the RENS trajectories.
514
515 @param traj_info: TrajectoryInfo object holding information neccessary to
516 calculate the rens trajectories.
517 @type traj_info: L{RENSTrajInfo}
518 @rtype: L{AbstractPropagator}
519 """
520 pass
521
523 """
524 Run the trajectory generator which generates a trajectory
525 of a given length between the states of two samplers.
526
527 @param traj_info: TrajectoryInfo instance holding information
528 needed to generate a nonequilibrium trajectory
529 @type traj_info: L{RENSTrajInfo}
530
531 @rtype: L{Trajectory}
532 """
533
534 init_temperature = self._get_init_temperature(traj_info)
535
536 init_state = traj_info.init_state.clone()
537
538 if init_state.momentum is None:
539 init_state = augment_state(init_state,
540 init_temperature,
541 traj_info.param_info.mass_matrix)
542
543 gen = self._propagator_factory(traj_info)
544
545 traj = gen.generate(init_state, int(traj_info.param_info.traj_length))
546
547 return traj
548
551 """
552 Holds parameters for a RENS swap.
553 """
554
555 __metaclass__ = ABCMeta
556
557 - def __init__(self, sampler1, sampler2, protocol):
568
569 @property
571 """
572 Switching protocol determining the time dependence
573 of the switching parameter.
574 """
575 return self._protocol
576 @protocol.setter
578 self._protocol = value
579
582 """
583 Holds all the information which needs to be communicated between distinct
584 RENS swap substeps.
585
586 See L{AbstractSwapCommunicator} for constructor signature.
587 """
588
589 pass
590
593 """
594 Holds information necessary for calculating a RENS trajectory.
595
596 @param param_info: ParameterInfo instance holding swap parameters
597 @type param_info: L{AbstractSwapParameterInfo}
598
599 @param init_state: state from which the trajectory is supposed to start
600 @type init_state: L{State}
601
602 @param direction: Either "fw" or "bw", indicating a forward or backward
603 trajectory. This is neccessary to pick the protocol or
604 the reversed protocol, respectively.
605 @type direction: string, either "fw" or "bw"
606 """
607
608 - def __init__(self, param_info, init_state, direction):
613
614 @property
616 return self._param_info
617
618 @property
620 return self._init_state
621
622 @property
624 return self._direction
625
626
627 -class MDRENS(AbstractRENS):
628 """
629 Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski 2009)
630 with Molecular Dynamics (MD) trajectories.
631
632 @param samplers: Samplers which sample their
633 respective equilibrium distributions
634 @type samplers: list of L{AbstractSingleChainMC}
635
636 @param param_infos: ParameterInfo instance holding
637 information required to perform a MDRENS swap
638 @type param_infos: list of L{MDRENSSwapParameterInfo}
639
640 @param integrator: Subclass of L{AbstractIntegrator} to be used to
641 calculate the non-equilibrium trajectories
642 @type integrator: type
643 """
644
651
663
693
696 """
697 Holds parameters for a MDRENS swap.
698
699 @param sampler1: First sampler
700 @type sampler1: L{AbstractSingleChainMC}
701
702 @param sampler2: Second sampler
703 @type sampler2: L{AbstractSingleChainMC}
704
705 @param timestep: Integration timestep
706 @type timestep: float
707
708 @param traj_length: Trajectory length in number of timesteps
709 @type traj_length: int
710
711 @param gradient: Gradient which determines the dynamics during a trajectory
712 @type gradient: L{AbstractGradient}
713
714 @param protocol: Switching protocol determining the time dependence of the
715 switching parameter. It is a function M{f} taking the running
716 time t and the switching time tau to yield a value in M{[0, 1]}
717 with M{f(0, tau) = 0} and M{f(tau, tau) = 1}. Default is a linear
718 protocol, which is being set manually due to an epydoc bug
719 @type protocol: callable
720
721 @param mass_matrix: Mass matrix
722 @type mass_matrix: n-dimensional matrix of type L{InvertibleMatrix} with n being the dimension
723 of the configuration space, that is, the dimension of
724 the position / momentum vectors
725 """
726
727 - def __init__(self, sampler1, sampler2, timestep, traj_length, gradient,
728 protocol=None, mass_matrix=None):
748
749 @property
751 """
752 Integration timestep.
753 """
754 return self._timestep
755 @timestep.setter
757 self._timestep = float(value)
758
759 @property
761 """
762 Trajectory length in number of integration steps.
763 """
764 return self._traj_length
765 @traj_length.setter
767 self._traj_length = int(value)
768
769 @property
771 """
772 Gradient which governs the equations of motion.
773 """
774 return self._gradient
775
776 @property
778 return self._mass_matrix
779 @mass_matrix.setter
781 self._mass_matrix = value
782
783 @property
785 """
786 Switching protocol determining the time dependence
787 of the switching parameter.
788 """
789 return self._protocol
790 @protocol.setter
792 self._protocol = value
793
796 """
797 Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski, 2009)
798 with Andersen-thermostatted Molecular Dynamics (MD) trajectories.
799
800 @param samplers: Samplers which sample their
801 respective equilibrium distributions
802 @type samplers: list of L{AbstractSingleChainMC}
803
804 @param param_infos: ParameterInfo instance holding
805 information required to perform a MDRENS swap
806 @type param_infos: list of L{ThermostattedMDRENSSwapParameterInfo}
807
808 @param integrator: Subclass of L{AbstractIntegrator} to be used to
809 calculate the non-equilibrium trajectories
810 @type integrator: type
811 """
812
816
834
836 """
837 @param sampler1: First sampler
838 @type sampler1: subclass instance of L{AbstractSingleChainMC}
839
840 @param sampler2: Second sampler
841 @type sampler2: subclass instance of L{AbstractSingleChainMC}
842
843 @param timestep: Integration timestep
844 @type timestep: float
845
846 @param traj_length: Trajectory length in number of timesteps
847 @type traj_length: int
848
849 @param gradient: Gradient which determines the dynamics during a trajectory
850 @type gradient: subclass instance of L{AbstractGradient}
851
852 @param mass_matrix: Mass matrix
853 @type mass_matrix: n-dimensional L{InvertibleMatrix} with n being the dimension
854 of the configuration space, that is, the dimension of
855 the position / momentum vectors
856
857 @param protocol: Switching protocol determining the time dependence of the
858 switching parameter. It is a function f taking the running
859 time t and the switching time tau to yield a value in [0, 1]
860 with f(0, tau) = 0 and f(tau, tau) = 1
861 @type protocol: callable
862
863 @param temperature: Temperature interpolation function.
864 @type temperature: Real-valued function mapping from [0,1] to R.
865 T(0) = temperature of the ensemble sampler1 samples from, T(1) = temperature
866 of the ensemble sampler2 samples from
867
868 @param collision_probability: Probability for a collision with the heatbath during one timestep
869 @type collision_probability: float
870
871 @param collision_interval: Interval during which collision may occur with probability
872 collision_probability
873 @type collision_interval: int
874 """
875
876 - def __init__(self, sampler1, sampler2, timestep, traj_length, gradient, mass_matrix=None,
877 protocol=None, temperature=lambda l: 1.0,
878 collision_probability=0.1, collision_interval=1):
890
891 @property
893 """
894 Probability for a collision with the heatbath during one timestep.
895 """
896 return self._collision_probability
897 @collision_probability.setter
899 self._collision_probability = float(value)
900
901 @property
903 """
904 Interval during which collision may occur with probability
905 C{collision_probability}.
906 """
907 return self._collision_interval
908 @collision_interval.setter
910 self._collision_interval = int(value)
911
912 @property
914 return self._temperature
915
918 """
919 Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski 2009)
920 with stepwise trajectories as described in Nilmeier et al., "Nonequilibrium candidate
921 Monte Carlo is an efficient tool for equilibrium simulation", PNAS 2011.
922 The switching parameter dependence of the Hamiltonian is a linear interpolation
923 between the PDFs of the sampler objects,
924 M{H(S{lambda}) = H_2 * S{lambda} + (1 - S{lambda}) * H_1}.
925 The perturbation kernel is a thermodynamic perturbation and the propagation is subclass
926 responsibility.
927 Note that due to the linear interpolations between the two Hamiltonians, the
928 log-probability has to be evaluated four times per perturbation step which can be
929 costly. In this case it is advisable to define the intermediate log probabilities
930 in _run_traj_generator differently.
931
932 @param samplers: Samplers which sample their respective equilibrium distributions
933 @type samplers: list of L{AbstractSingleChainMC}
934
935 @param param_infos: ParameterInfo instances holding
936 information required to perform a HMCStepRENS swaps
937 @type param_infos: list of L{AbstractSwapParameterInfo}
938 """
939
940 __metaclass__ = ABCMeta
941
942 - def __init__(self, samplers, param_infos):
947
948 @abstractmethod
950 """
951 Set up the propagation steps using the information about the current system
952 setup and parameters from the SwapParameterInfo object.
953
954 @param im_sys_infos: Information about the intermediate system setups
955 @type im_sys_infos: List of L{AbstractSystemInfo}
956
957 @param param_info: SwapParameterInfo object containing parameters for the
958 propagations like timesteps, trajectory lengths etc.
959 @type param_info: L{AbstractSwapParameterInfo}
960 """
961
962 pass
963
964 @abstractmethod
966 """
967 If needed, set im_sys_infos.hamiltonian.gradient.
968
969 @param im_sys_infos: Information about the intermediate system setups
970 @type im_sys_infos: List of L{AbstractSystemInfo}
971
972 @param param_info: SwapParameterInfo object containing parameters for the
973 propagations like timesteps, trajectory lengths etc.
974 @type param_info: L{AbstractSwapParameterInfo}
975
976 @param t_prot: Switching protocol defining the time dependence of the switching
977 parameter.
978 @type t_prot: callable
979 """
980
981 pass
982
1024
1031
1042
1045 """
1046 Replica Exchange with Nonequilibrium Switches (RENS, Ballard & Jarzynski 2009)
1047 with stepwise trajectories as described in Nilmeier et al., "Nonequilibrium candidate
1048 Monte Carlo is an efficient tool for equilibrium simulation", PNAS 2011.
1049 The switching parameter dependence of the Hamiltonian is a linear interpolation
1050 between the PDFs of the sampler objects,
1051 M{H(S{lambda}) = H_2 * S{lambda} + (1 - S{lambda}) * H_1}.
1052 The perturbation kernel is a thermodynamic perturbation and the propagation is done using HMC.
1053
1054 Note that due to the linear interpolations between the two Hamiltonians, the
1055 log-probability and its gradient has to be evaluated four times per perturbation step which
1056 can be costly. In this case it is advisable to define the intermediate log probabilities
1057 in _run_traj_generator differently.
1058
1059 @param samplers: Samplers which sample their respective equilibrium distributions
1060 @type samplers: list of L{AbstractSingleChainMC}
1061
1062 @param param_infos: ParameterInfo instances holding
1063 information required to perform a HMCStepRENS swaps
1064 @type param_infos: list of L{HMCStepRENSSwapParameterInfo}
1065 """
1066
1067 - def __init__(self, samplers, param_infos):
1070
1080
1095
1099
1102 """
1103 Holds all required information for performing HMCStepRENS swaps.
1104
1105 @param sampler1: First sampler
1106 @type sampler1: subclass instance of L{AbstractSingleChainMC}
1107
1108 @param sampler2: Second sampler
1109 @type sampler2: subclass instance of L{AbstractSingleChainMC}
1110
1111 @param timestep: integration timestep
1112 @type timestep: float
1113
1114 @param hmc_traj_length: HMC trajectory length
1115 @type hmc_traj_length: int
1116
1117 @param hmc_iterations: number of HMC iterations in the propagation step
1118 @type hmc_iterations: int
1119
1120 @param gradient: gradient governing the equations of motion, function of
1121 position array and switching protocol
1122 @type gradient: callable
1123
1124 @param intermediate_steps: number of steps in the protocol; this is a discrete version
1125 of the switching time in "continuous" RENS implementations
1126 @type intermediate_steps: int
1127
1128 @param protocol: Switching protocol determining the time dependence of the
1129 switching parameter. It is a function f taking the running
1130 time t and the switching time tau to yield a value in [0, 1]
1131 with f(0, tau) = 0 and f(tau, tau) = 1
1132 @type protocol: callable
1133
1134 @param mass_matrix: mass matrix for kinetic energy definition
1135 @type mass_matrix: L{InvertibleMatrix}
1136
1137 @param integrator: Integration scheme to be utilized
1138 @type integrator: l{AbstractIntegrator}
1139 """
1140
1141 - def __init__(self, sampler1, sampler2, timestep, hmc_traj_length, hmc_iterations,
1142 gradient, intermediate_steps, parametrization=None, protocol=None,
1143 mass_matrix=None, integrator=FastLeapFrog):
1165
1166 @property
1168 """
1169 Integration timestep.
1170 """
1171 return self._timestep
1172 @timestep.setter
1174 self._timestep = float(value)
1175
1176 @property
1178 """
1179 HMC trajectory length in number of integration steps.
1180 """
1181 return self._hmc_traj_length
1182 @hmc_traj_length.setter
1184 self._hmc_traj_length = int(value)
1185
1186 @property
1188 """
1189 Gradient which governs the equations of motion.
1190 """
1191 return self._gradient
1192 @gradient.setter
1194 self._gradient = value
1195
1196 @property
1198 return self._mass_matrix
1199 @mass_matrix.setter
1201 self._mass_matrix = value
1202
1203 @property
1205 return self._hmc_iterations
1206 @hmc_iterations.setter
1208 self._hmc_iterations = value
1209
1210 @property
1213 @intermediate_steps.setter
1216
1217 @property
1219 return self._integrator
1220 @integrator.setter
1222 self._integrator = value
1223
1226 """
1227 Provides the interface for classes defining schemes according to which swaps in
1228 Replica Exchange-like simulations are performed.
1229
1230 @param algorithm: Exchange algorithm that performs the swaps
1231 @type algorithm: L{AbstractExchangeMC}
1232 """
1233
1234 __metaclass__ = ABCMeta
1235
1237
1238 self._algorithm = algorithm
1239
1240 @abstractmethod
1242 """
1243 Advises the Replica Exchange-like algorithm to perform swaps according to
1244 the schedule defined here.
1245 """
1246
1247 pass
1248
1251 """
1252 Provides a swapping scheme in which tries exchanges between neighbours only
1253 following the scheme 1 <-> 2, 3 <-> 4, ... and after a sampling period 2 <-> 3, 4 <-> 5, ...
1254
1255 @param algorithm: Exchange algorithm that performs the swaps
1256 @type algorithm: L{AbstractExchangeMC}
1257 """
1258
1260
1261 super(AlternatingAdjacentSwapScheme, self).__init__(algorithm)
1262
1263 self._current_swap_list = None
1264 self._swap_list1 = []
1265 self._swap_list2 = []
1266 self._create_swap_lists()
1267
1269
1270 if len(self._algorithm.param_infos) == 1:
1271 self._swap_list1.append(0)
1272 self._swap_list2.append(0)
1273 else:
1274 i = 0
1275 while i < len(self._algorithm.param_infos):
1276 self._swap_list1.append(i)
1277 i += 2
1278
1279 i = 1
1280 while i < len(self._algorithm.param_infos):
1281 self._swap_list2.append(i)
1282 i += 2
1283
1284 self._current_swap_list = self._swap_list1
1285
1287
1288 for x in self._current_swap_list:
1289 self._algorithm.swap(x)
1290
1291 if self._current_swap_list == self._swap_list1:
1292 self._current_swap_list = self._swap_list2
1293 else:
1294 self._current_swap_list = self._swap_list1
1295
1298 """
1299 Tracks swap statistics of a single sampler pair.
1300
1301 @param param_info: ParameterInfo instance holding swap parameters
1302 @type param_info: L{AbstractSwapParameterInfo}
1303 """
1304
1306 self._total_swaps = 0
1307 self._accepted_swaps = 0
1308
1309 @property
1311 return self._total_swaps
1312
1313 @property
1315 return self._accepted_swaps
1316
1317 @property
1319 """
1320 Acceptance rate of the sampler pair.
1321 """
1322 if self.total_swaps > 0:
1323 return float(self.accepted_swaps) / float(self.total_swaps)
1324 else:
1325 return 0.
1326
1328 """
1329 Updates swap statistics.
1330 """
1331 self._total_swaps += 1
1332 self._accepted_swaps += int(accepted)
1333
1336 """
1337 Tracks swap statistics for an AbstractExchangeMC subclass instance.
1338
1339 @param param_infos: list of ParameterInfo instances providing information
1340 needed for performing swaps
1341 @type param_infos: list of L{AbstractSwapParameterInfo}
1342 """
1343
1345 self._stats = [SingleSwapStatistics(x) for x in param_infos]
1346
1347 @property
1349 return tuple(self._stats)
1350
1351 @property
1353 """
1354 Returns acceptance rates for all swaps.
1355 """
1356 return [x.acceptance_rate for x in self._stats]
1357
1360 """
1361 Produces interpolations for functions changed during non-equilibrium
1362 trajectories.
1363
1364 @param protocol: protocol to be used to generate non-equilibrium trajectories
1365 @type protocol: function mapping t to [0...1] for fixed tau
1366 @param tau: switching time
1367 @type tau: float
1368 """
1369
1377
1378 @property
1380 return self._protocol
1381 @protocol.setter
1383 if not hasattr(value, '__call__'):
1384 raise TypeError(value)
1385 self._protocol = value
1386
1387 @property
1390 @tau.setter
1391 - def tau(self, value):
1392 self._tau = float(value)
1393
1395 """
1396 Create a gradient instance with according to given protocol
1397 and switching time.
1398
1399 @param gradient: gradient with G(0) = G_1 and G(1) = G_2
1400 @type gradient: callable
1401 """
1402 return Gradient(gradient, self._protocol, self._tau)
1403
1405 """
1406 Create a temperature function according to given protocol and
1407 switching time.
1408
1409 @param temperature: temperature with T(0) = T_1 and T(1) = T_2
1410 @type temperature: callable
1411 """
1412 return lambda t: temperature(self.protocol(t, self.tau))
1413
1416
1417 - def __init__(self, gradient, protocol, tau):
1422
1424 return self._gradient(q, self._protocol(t, self._tau))
1425
1426
1427 -class ReplicaHistory(object):
1428 '''
1429 Replica history object, works with both RE and RENS for
1430 the AlternatingAdjacentSwapScheme.
1431
1432 @param samples: list holding ensemble states
1433 @type samples: list
1434
1435 @param swap_interval: interval with which swaps were attempted, e.g.,
1436 5 means that every 5th regular MC step is replaced
1437 by a swap
1438 @type swap_interval: int
1439
1440 @param first_swap: sample index of the first sample generated by a swap attempt.
1441 If None, the first RE sampled is assumed to have sample index
1442 swap_interval. If specified, it has to be greater than zero
1443 @type first_swap: int
1444 '''
1445
1446 - def __init__(self, samples, swap_interval, first_swap=None):
1447 self.samples = samples
1448 self.swap_interval = swap_interval
1449 if first_swap == None:
1450 self.first_swap = swap_interval - 1
1451 elif first_swap > 0:
1452 self.first_swap = first_swap - 1
1453 else:
1454 raise(ValueError("Sample index of first swap has to be greater than zero!"))
1455 self.n_replicas = len(samples[0])
1456
1457 @staticmethod
1459 if x == 1:
1460 return -1
1461 if x == -1:
1462 return 1
1463
1464 - def calculate_history(self, start_ensemble):
1465 '''
1466 Calculates the replica history of the first state of ensemble #start_ensemble.
1467
1468 @param start_ensemble: index of the ensemble to start at, zero-indexed
1469 @type start_ensemble: int
1470
1471 @return: replica history as a list of ensemble indices
1472 @rtype: list of ints
1473 '''
1474
1475 sample_counter = 0
1476
1477
1478
1479 if start_ensemble % 2 == 0:
1480 direction = +1
1481 else:
1482 direction = -1
1483
1484
1485
1486 if start_ensemble % 2 == 0 and start_ensemble == self.n_replicas - 1:
1487 direction = -1
1488
1489
1490 history = []
1491
1492
1493 ens = start_ensemble
1494
1495 while sample_counter < len(self.samples):
1496 if self.n_replicas == 2:
1497 if (sample_counter - self.first_swap - 1) % self.swap_interval == 0 and \
1498 sample_counter >= self.first_swap:
1499
1500
1501 current_sample = self.samples[sample_counter][ens]
1502
1503
1504 previous_sample = self.samples[sample_counter - 1][history[-1]]
1505
1506
1507
1508
1509 swap_accepted = not numpy.all(current_sample.position ==
1510 previous_sample.position)
1511
1512 if swap_accepted:
1513 if ens == 0:
1514 ens = 1
1515 else:
1516 ens = 0
1517 history.append(ens)
1518 else:
1519 history.append(ens)
1520
1521 else:
1522
1523 if (sample_counter - self.first_swap - 1) % self.swap_interval == 0 and \
1524 sample_counter >= self.first_swap:
1525
1526 current_sample = self.samples[sample_counter][ens]
1527
1528
1529 previous_sample = self.samples[sample_counter - 1][ens]
1530
1531
1532
1533
1534 swap_accepted = not numpy.all(current_sample.position == previous_sample.position)
1535
1536 if swap_accepted:
1537 ens += direction
1538 else:
1539 if ens == self.n_replicas - 1:
1540
1541 direction = -1
1542 elif ens == 0:
1543
1544 direction = +1
1545 else:
1546
1547
1548 direction = self._change_direction(direction)
1549
1550 history.append(ens)
1551
1552 sample_counter += 1
1553
1554 return history
1555
1557 '''
1558 Calculates sequentially correlated trajectories projected on a specific ensemble.
1559
1560 @param ensemble: ensemble index of ensemble of interest, zero-indexed
1561 @type ensemble: int
1562
1563 @return: list of Trajectory objects containg sequentially correlated trajectories
1564 @rtype: list of L{Trajectory} objects.
1565 '''
1566
1567 trajectories = []
1568
1569 for i in range(self.n_replicas):
1570 history = self.calculate_history(i)
1571 traj = [x[ensemble] for k, x in enumerate(self.samples) if history[k] == ensemble]
1572 trajectories.append(Trajectory(traj))
1573
1574 return trajectories
1575
1577 '''
1578 Calculates sequentially correlated trajectories.
1579
1580 @return: list of Trajectory objects containg sequentially correlated trajectories
1581 @rtype: list of L{Trajectory} objects.
1582 '''
1583
1584 trajectories = []
1585
1586 for i in range(self.n_replicas):
1587 history = self.calculate_history(i)
1588 traj = [x[history[k]] for k, x in enumerate(self.samples)]
1589 trajectories.append(Trajectory(traj))
1590
1591 return trajectories
1592