1 """
2 Abstract Monte Carlo samplers.
3 """
4
5 import numpy.random
6
7 import csb.numeric
8 import csb.core
9
10 from abc import ABCMeta, abstractmethod, abstractproperty
11 from csb.statistics.samplers import AbstractSampler, AbstractState, State, EnsembleState
14 """
15 Abstract Monte Carlo sampler class. Subclasses implement various
16 Monte carlo equilibrium sampling schemes.
17
18 @param state: Initial state
19 @type state: L{AbstractState}
20 """
21
22 __metaclass__ = ABCMeta
23
28
33
34 @abstractproperty
36 """
37 Energy of the current state.
38 """
39 pass
40
41 @property
43 """
44 Current state.
45 """
46 return self._state
47 @state.setter
51
52 @abstractmethod
54 """
55 Draw a sample.
56 @rtype: L{AbstractState}
57 """
58 pass
59
61 """
62 Abstract class providing the interface for the result
63 of a deterministic or stochastic propagation of a state.
64 """
65
66 __metaclass__ = ABCMeta
67
68 @abstractproperty
70 """
71 Initial state
72 """
73 pass
74
75 @abstractproperty
77 """
78 Final state
79 """
80 pass
81
82 @abstractproperty
84 """
85 Heat produced during propagation
86 @rtype: float
87 """
88 pass
89
91 """
92 Describes the result of a deterministic or stochastic
93 propagation of a state.
94
95 @param initial: Initial state from which the
96 propagation started
97 @type initial: L{State}
98
99 @param final: Final state in which the propagation
100 resulted
101 @type final: L{State}
102
103 @param heat: Heat produced during propagation
104 @type heat: float
105 """
106
107
108 - def __init__(self, initial, final, heat=0.0):
121
122 @property
125
126 @property
129
130 @property
133 @heat.setter
134 - def heat(self, value):
135 self._heat = float(value)
136
137 -class Trajectory(csb.core.CollectionContainer, AbstractPropagationResult):
138 """
139 Ordered collection of states, representing a phase-space trajectory.
140
141 @param items: list of states defining a phase-space trajectory
142 @type items: list of L{AbstractState}
143 @param heat: heat produced during the trajectory
144 @type heat: float
145 @param work: work produced during the trajectory
146 @type work: float
147 """
148
149 - def __init__(self, items, heat=0.0, work=0.0):
155
156 @property
159
160 @property
163
164 @property
167 @heat.setter
168 - def heat(self, value):
169 self._heat = float(value)
170
171 @property
174 @work.setter
175 - def work(self, value):
176 self._work = float(value)
177
179 """
180 Allows to build a Trajectory object step by step.
181
182 @param heat: heat produced over the trajectory
183 @type heat: float
184 @param work: work produced during the trajectory
185 @type work: float
186 """
187
188 - def __init__(self, heat=0.0, work=0.0):
189 self._heat = heat
190 self._work = work
191 self._states = []
192
193 @staticmethod
195 """
196 Trajectory builder factory.
197
198 @param full: if True, a TrajectoryBuilder instance designed
199 to build a full trajectory with initial state,
200 intermediate states and a final state. If False,
201 a ShortTrajectoryBuilder instance designed to
202 hold only the initial and the final state is
203 returned
204 @type full: boolean
205 """
206
207 if full:
208 return TrajectoryBuilder()
209 else:
210 return ShortTrajectoryBuilder()
211
212 @property
214 """
215 The L{Trajectory} instance build by a specific instance of
216 this class
217 """
218 return Trajectory(self._states, heat=self._heat, work=self._work)
219
221 """
222 Inserts a state at the beginning of the trajectory
223
224 @param state: state to be added
225 @type state: L{State}
226 """
227 self._states.insert(0, state.clone())
228
237
239 """
240 Adds a state to the end of the trajectory
241
242 @param state: state to be added
243 @type state: L{State}
244 """
245 self._states.append(state.clone())
246
248
251
252 @property
254 """
255 The L{PropagationResult} instance built by a specific instance of
256 this class
257 """
258
259 if len(self._states) != 2:
260 raise ValueError("Can't create a product, two states required")
261
262 initial, final = self._states
263 return PropagationResult(initial, final, heat=self._heat)
264
266 """
267 With the exception of the current state of the Markov chain, this
268 holds all the information needed to calculate the acceptance
269 probability in both the L{RWMCSampler} and L{HMCSampler} classes,
270 that is, only the proposal state.
271 For more advanced algorithms, one may derive classes capable of
272 holding the neccessary additional information from this class.
273
274 @param proposal: Proposal state
275 @type proposal: L{State}
276 """
277
278 __metaclass__ = ABCMeta
279
283
284 @property
286 return self._proposal
287
289 """
290 Abstract class for Monte Carlo sampling algorithms simulating
291 only one ensemble.
292
293 @param pdf: probability density function to sample from
294 @type pdf: subclass of L{csb.statistics.pdf.AbstractDensity}
295
296 @param state: Initial state
297 @type state: L{State}
298
299 @param temperature: Pseudo-temperature of the Boltzmann ensemble
300 M{p(x) = 1/N * exp(-1/T * E(x))} with the
301 pseudo-energy defined as M{E(x) = -log(p(x))}
302 where M{p(x)} is the PDF under consideration
303 @type temperature: float
304 """
305
306 __metaclass__ = ABCMeta
307
308 - def __init__(self, pdf, state, temperature=1.):
317
321
323 """
324 Draw a sample.
325 @rtype: L{State}
326 """
327
328 proposal_communicator = self._propose()
329 pacc = self._calc_pacc(proposal_communicator)
330 self._accept(proposal_communicator.proposal, pacc)
331
332 return self.state
333
334 @abstractmethod
336 """
337 Calculate a new proposal state and gather additional information
338 needed to calculate the acceptance probability.
339
340 @rtype: L{SimpleProposalCommunicator}
341 """
342 pass
343
344 @abstractmethod
346 """
347 Calculate probability with which to accept the proposal.
348
349 @param proposal_communicator: Contains information about the proposal
350 and additional information needed to
351 calculate the acceptance probability
352 @type proposal_communicator: L{SimpleProposalCommunicator}
353 """
354 pass
355
356 - def _accept(self, proposal, pacc):
357 """
358 Accept / reject proposal with given acceptance probability pacc.
359
360 @param proposal: proposal state
361 @type proposal: L{State}
362
363 @param pacc: acceptance probability
364 @type pacc: float
365 """
366
367 self._nmoves += 1
368
369 if numpy.random.random() < pacc:
370 self.state = proposal
371 self._accepted += 1
372 self._last_move_accepted = True
373 return True
374 else:
375 self._last_move_accepted = False
376 return False
377
378 @property
380 """
381 Log-likelihood of the current state.
382 @rtype: float
383 """
384 return self._pdf.log_prob(self.state.position)
385
386 @property
388 """
389 Acceptance rate.
390 """
391 return float(self._accepted) / float(self._nmoves)
392
393 @property
395 """
396 Information whether the last MC move was accepted or not.
397 """
398 return self._last_move_accepted
399
400 @property
402 return self._temperature
403
405 """
406 Collection of single-chain samplers.
407
408 @param items: samplers
409 @type items: list of L{AbstractSingleChainMC}
410 """
411
415
417 """
418 Abstract class for Monte Carlo sampling algorithms simulating several ensembles.
419
420 @param samplers: samplers which sample from their respective equilibrium distributions
421 @type samplers: list of L{AbstractSingleChainMC}
422 """
423
424 __metaclass__ = ABCMeta
425
432
444
445 @property
447 """
448 Total ensemble energy.
449 """
450 return sum([x.energy for x in self._samplers])
451
453 """
454 Subclass instances hold all parameters necessary for performing a swap
455 between two given samplers.
456 """
457
458 __metaclass__ = ABCMeta
459
460 - def __init__(self, sampler1, sampler2):
461 """
462 @param sampler1: First sampler
463 @type sampler1: L{AbstractSingleChainMC}
464
465 @param sampler2: Second sampler
466 @type sampler2: L{AbstractSingleChainMC}
467 """
468
469 self._sampler1 = sampler1
470 self._sampler2 = sampler2
471
472 @property
474 return self._sampler1
475
476 @property
478 return self._sampler2
479
481 """
482 Holds parameters for a standard Replica Exchange swap.
483 """
484 pass
485
487 """
488 Holds parameters for a MDRENS swap.
489
490 @param sampler1: First sampler
491 @type sampler1: L{AbstractSingleChainMC}
492
493 @param sampler2: Second sampler
494 @type sampler2: L{AbstractSingleChainMC}
495
496 @param timestep: Integration timestep
497 @type timestep: float
498
499 @param mass_matrix: Mass matrix
500 @type mass_matrix: n-dimensional matrix of type L{InvertibleMatrix} with n being the dimension
501 of the configuration space, that is, the dimension of
502 the position / momentum vectors
503
504 @param traj_length: Trajectory length in number of timesteps
505 @type traj_length: int
506
507 @param gradient: Gradient which determines the dynamics during a trajectory
508 @type gradient: L{AbstractGradient}
509 """
510
511 - def __init__(self, sampler1, sampler2, timestep, traj_length, gradient, mass_matrix=None):
523
524 @property
526 """
527 Integration timestep.
528 """
529 return self._timestep
530 @timestep.setter
532 self._timestep = float(value)
533
534 @property
536 """
537 Trajectory length in number of integration steps.
538 """
539 return self._traj_length
540 @traj_length.setter
542 self._traj_length = int(value)
543
544 @property
546 """
547 Gradient which governs the equations of motion.
548 """
549 return self._gradient
550
551 @property
553 return self._mass_matrix
554 @mass_matrix.setter
556 self._mass_matrix = value
557
559 """
560 @param sampler1: First sampler
561 @type sampler1: subclass instance of L{AbstractSingleChainMC}
562
563 @param sampler2: Second sampler
564 @type sampler2: subclass instance of L{AbstractSingleChainMC}
565
566 @param timestep: Integration timestep
567 @type timestep: float
568
569 @param mass_matrix: Mass matrix
570 @type mass_matrix: n-dimensional L{InvertibleMatrix} with n being the dimension
571 of the configuration space, that is, the dimension of
572 the position / momentum vectors
573
574 @param traj_length: Trajectory length in number of timesteps
575 @type traj_length: int
576
577 @param gradient: Gradient which determines the dynamics during a trajectory
578 @type gradient: subclass instance of L{AbstractGradient}
579
580 @param temperature: Temperature interpolation function.
581 @type temperature: Real-valued function mapping from [0,1] to R.
582 T(0) = temperature of the ensemble sampler1 samples from, T(1) = temperature
583 of the ensemble sampler2 samples from
584
585 @param collision_probability: Probability for a collision with the heatbath during one timestep
586 @type collision_probability: float
587
588 @param collision_interval: Interval during which collision may occur with probability
589 collision_probability
590 @type collision_interval: int
591 """
592
605
606 @property
608 """
609 Probability for a collision with the heatbath during one timestep.
610 """
611 return self._collision_probability
612 @collision_probability.setter
614 self._collision_probability = float(value)
615
616 @property
618 """
619 Interval during which collision may occur with probability
620 C{collision_probability}.
621 """
622 return self._collision_interval
623 @collision_interval.setter
625 self._collision_interval = int(value)
626
627 @property
629 return self._temperature
630
632 """
633 Holds all the information which needs to be communicated between
634 distinct swap substeps.
635
636 @param param_info: ParameterInfo instance holding swap parameters
637 @type param_info: L{AbstractSwapParameterInfo}
638
639 @param traj12: Forward trajectory
640 @type traj12: L{Trajectory}
641
642 @param traj21: Reverse trajectory
643 @type traj21: L{Trajectory}
644 """
645
646 __metaclass__ = ABCMeta
647
648 - def __init__(self, param_info, traj12, traj21):
660
661 @property
663 return self._sampler1
664
665 @property
667 return self._sampler2
668
669 @property
672
673 @property
676
677 @property
679 return self._acceptance_probability
680 @acceptance_probability.setter
682 self._acceptance_probability = value
683
684 @property
686 return self._accepted
687 @accepted.setter
689 self._accepted = value
690
691 @property
693 return self._param_info
694
696 """
697 Holds all the information which needs to be communicated between distinct
698 RE swap substeps.
699
700 See L{AbstractSwapCommunicator} for constructor signature.
701 """
702 pass
703
705 """
706 Holds all the information which needs to be communicated between distinct
707 RENS swap substeps.
708
709 See L{AbstractSwapCommunicator} for constructor signature.
710 """
711
712 pass
713
715 """
716 Tracks swap statistics of a single sampler pair.
717
718 @param param_info: ParameterInfo instance holding swap parameters
719 @type param_info: L{AbstractSwapParameterInfo}
720 """
721
723 self._total_swaps = 0
724 self._accepted_swaps = 0
725
726 @property
728 return self._total_swaps
729
730 @property
732 return self._accepted_swaps
733
734 @property
743
745 """
746 Updates swap statistics.
747 """
748 self._total_swaps += 1
749 self._accepted_swaps += int(accepted)
750
752 """
753 Tracks swap statistics for an AbstractExchangeMC subclass instance.
754
755 @param param_infos: list of ParameterInfo instances providing information
756 needed for performing swaps
757 @type param_infos: list of L{AbstractSwapParameterInfo}
758 """
759
761 self._stats = [SingleSwapStatistics(x) for x in param_infos]
762
763 @property
765 return tuple(self._stats)
766
767 @property
769 """
770 Returns acceptance rates for all swaps.
771 """
772 return [x.acceptance_rate for x in self._stats]
773
775 """
776 Abstract class for Monte Carlo sampling algorithms employing some replica exchange method.
777
778 @param samplers: samplers which sample from their respective equilibrium distributions
779 @type samplers: list of L{AbstractSingleChainMC}
780
781 @param param_infos: list of ParameterInfo instances providing information needed
782 for performing swaps
783 @type param_infos: list of L{AbstractSwapParameterInfo}
784 """
785
786 __metaclass__ = ABCMeta
787
788 - def __init__(self, samplers, param_infos):
797
801
802 - def swap(self, index):
803 """
804 Perform swap between sampler pair described by param_infos[index]
805 and return outcome (true = accepted, false = rejected).
806
807 @param index: index of swap pair in param_infos
808 @type index: int
809
810 @rtype: boolean
811 """
812 param_info = self._param_infos[index]
813 swapcom = self._propose_swap(param_info)
814 swapcom = self._calc_pacc_swap(swapcom)
815 result = self._accept_swap(swapcom)
816
817 self.state = EnsembleState([x.state for x in self._samplers])
818
819 self.statistics.stats[index].update(result)
820
821 return result
822
823 @abstractmethod
825 """
826 Calculate proposal states for a swap between two samplers.
827
828 @param param_info: ParameterInfo instance holding swap parameters
829 @type param_info: L{AbstractSwapParameterInfo}
830
831 @rtype: L{AbstractSwapCommunicator}
832 """
833 pass
834
835 @abstractmethod
837 """
838 Calculate probability to accept a swap given initial and proposal states.
839
840 @param swapcom: SwapCommunicator instance holding information to be communicated
841 between distinct swap substeps
842 @type swapcom: L{AbstractSwapCommunicator}
843
844 @rtype: L{AbstractSwapCommunicator}
845 """
846 pass
847
849 """
850 Accept / reject an exchange between two samplers given proposal states and
851 the acceptance probability and returns the outcome (true = accepted, false = rejected).
852
853 @param swapcom: SwapCommunicator instance holding information to be communicated
854 between distinct swap substeps
855 @type swapcom: L{AbstractSwapCommunicator}
856
857 @rtype: boolean
858 """
859
860 if numpy.random.random() < swapcom.acceptance_probability:
861 swapcom.sampler1.state = swapcom.traj21.final
862 swapcom.sampler2.state = swapcom.traj12.final
863 return True
864 else:
865 return False
866
867 @property
875
876 @property
878 """
879 List of SwapParameterInfo instances holding all necessary parameters.
880
881 @rtype: list of L{AbstractSwapParameterInfo}
882 """
883 return self._param_infos
884
885 @property
887 return self._statistics
888
890 """
891 Update statistics of a given swap process.
892
893 @param index: position of swap statistics to be updated
894 @type index: int
895
896 @param accepted: outcome of the swap
897 @type accepted: boolean
898 """
899
900 self._stats[index][0] += 1
901 self._stats[index][1] += int(accepted)
902
904 """
905 Holds information necessary for calculating a RENS trajectory.
906
907 @param param_info: ParameterInfo instance holding swap parameters
908 @type param_info: L{AbstractSwapParameterInfo}
909
910 @param init_state: state from which the trajectory is supposed to start
911 @type init_state: L{State}
912
913 @param protocol: Protocol to be used to generate nonequilibrium trajectories
914 @type protocol: Real-valued function that maps [0, switching time] to [0, 1]
915 """
916
917 - def __init__(self, param_info, init_state, protocol):
922
923 @property
925 return self._param_info
926
927 @property
929 return self._protocol
930
931 @property
933 return self._init_state
934
936 """
937 Abstract Replica Exchange with Nonequilibrium Switches
938 (RENS, Ballard & Jarzynski 2009) class.
939 Subclasses implement various ways of generating trajectories
940 (deterministic or stochastic).
941 """
942
943 __metaclass__ = ABCMeta
944
946
947 T1 = param_info.sampler1.temperature
948 T2 = param_info.sampler2.temperature
949
950 momentum_covariance_matrix1 = T1 * param_info.mass_matrix
951 momentum_covariance_matrix2 = T2 * param_info.mass_matrix
952
953 d = len(param_info.sampler1.state.position)
954
955 if param_info.mass_matrix.is_unity_multiple:
956 momentum1 = numpy.random.normal(scale=numpy.sqrt(T1 * param_info.mass_matrix[0][0]),
957 size=d)
958 momentum2 = numpy.random.normal(scale=numpy.sqrt(T2 * param_info.mass_matrix[0][0]),
959 size=d)
960 else:
961 momentum1 = numpy.random.multivariate_normal(mean=numpy.zeros(d),
962 cov=momentum_covariance_matrix1)
963 momentum2 = numpy.random.multivariate_normal(mean=numpy.zeros(d),
964 cov=momentum_covariance_matrix2)
965
966 init_state1 = State(param_info.sampler1.state.position, momentum1)
967 init_state2 = State(param_info.sampler2.state.position, momentum2)
968
969 param_info.sampler1.state = init_state1
970 param_info.sampler2.state = init_state2
971
972 trajinfo12 = RENSTrajInfo(param_info, init_state1, protocol=lambda t, tau: t / tau)
973 trajinfo21 = RENSTrajInfo(param_info, init_state2, protocol=lambda t, tau: (tau - t) / tau)
974
975 traj12 = self._run_traj_generator(trajinfo12)
976 traj21 = self._run_traj_generator(trajinfo21)
977
978 return RENSSwapCommunicator(param_info, traj12, traj21)
979
1011
1012 @abstractmethod
1014 """
1015 Run the trajectory generator which generates a trajectory
1016 of a given length between the states of two samplers.
1017
1018 @param traj_info: TrajectoryInfo instance holding information
1019 needed to generate a nonequilibrium trajectory
1020 @type traj_info: L{RENSTrajInfo}
1021
1022 @rtype: L{Trajectory}
1023 """
1024 pass
1025
1027 """
1028 Provides the interface for classes defining schemes according to which swaps in
1029 Replica Exchange-like simulations are performed.
1030
1031 @param algorithm: Exchange algorithm that performs the swaps
1032 @type algorithm: L{AbstractExchangeMC}
1033 """
1034
1035 __metaclass__ = ABCMeta
1036
1038
1039 self._algorithm = algorithm
1040
1041 @abstractmethod
1043 """
1044 Advises the Replica Exchange-like algorithm to perform swaps according to
1045 the some schedule defined here.
1046 """
1047
1048 pass
1049
1051 """
1052 Provides a swapping scheme in which tries exchanges between neighbours only
1053 following the scheme 1 <-> 2, 3 <-> 4, ... and after a sampling period 2 <-> 3, 4 <-> 5, ...
1054
1055 @param algorithm: Exchange algorithm that performs the swaps
1056 @type algorithm: L{AbstractExchangeMC}
1057 """
1058
1060
1061 super(AlternatingAdjacentSwapScheme, self).__init__(algorithm)
1062
1063 self._current_swap_list = None
1064 self._swap_list1 = []
1065 self._swap_list2 = []
1066 self._create_swap_lists()
1067
1069
1070 if len(self._algorithm.param_infos) == 1:
1071 self._swap_list1.append(0)
1072 self._swap_list2.append(0)
1073 else:
1074 i = 0
1075 while i < len(self._algorithm.param_infos):
1076 self._swap_list1.append(i)
1077 i += 2
1078
1079 i = 1
1080 while i < len(self._algorithm.param_infos):
1081 self._swap_list2.append(i)
1082 i += 2
1083
1084 self._current_swap_list = self._swap_list1
1085
1087
1088 for x in self._current_swap_list:
1089 self._algorithm.swap(x)
1090
1091 if self._current_swap_list == self._swap_list1:
1092 self._current_swap_list = self._swap_list2
1093 else:
1094 self._current_swap_list = self._swap_list1
1095