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

Source Code for Package csb.statistics.samplers.mc

   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 
12 13 -class AbstractMC(AbstractSampler):
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
24 - def __init__(self, state):
25 26 self._state = None 27 self.state = state
28
29 - def _checkstate(self, state):
30 31 if not isinstance(state, AbstractState): 32 raise TypeError(state)
33 34 @abstractproperty
35 - def energy(self):
36 """ 37 Energy of the current state. 38 """ 39 pass
40 41 @property
42 - def state(self):
43 """ 44 Current state. 45 """ 46 return self._state
47 @state.setter
48 - def state(self, value):
49 self._checkstate(value) 50 self._state = value
51 52 @abstractmethod
53 - def sample(self):
54 """ 55 Draw a sample. 56 @rtype: L{AbstractState} 57 """ 58 pass
59
60 -class AbstractPropagationResult(object):
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
69 - def initial(self):
70 """ 71 Initial state 72 """ 73 pass
74 75 @abstractproperty
76 - def final(self):
77 """ 78 Final state 79 """ 80 pass
81 82 @abstractproperty
83 - def heat(self):
84 """ 85 Heat produced during propagation 86 @rtype: float 87 """ 88 pass
89
90 -class PropagationResult(AbstractPropagationResult):
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):
109 110 if not isinstance(initial, AbstractState): 111 raise TypeError(initial) 112 113 if not isinstance(final, AbstractState): 114 raise TypeError(final) 115 116 self._initial = initial 117 self._final = final 118 self._heat = None 119 120 self.heat = heat
121 122 @property
123 - def initial(self):
124 return self._initial
125 126 @property
127 - def final(self):
128 return self._final
129 130 @property
131 - def heat(self):
132 return self._heat
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):
150 151 super(Trajectory, self).__init__(items, type=AbstractState) 152 153 self._heat = heat 154 self._work = work
155 156 @property
157 - def initial(self):
158 return self[0]
159 160 @property
161 - def final(self):
162 return self[self.last_index]
163 164 @property
165 - def heat(self):
166 return self._heat
167 @heat.setter
168 - def heat(self, value):
169 self._heat = float(value)
170 171 @property
172 - def work(self):
173 return self._work
174 @work.setter
175 - def work(self, value):
176 self._work = float(value)
177
178 -class TrajectoryBuilder(object):
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
194 - def create(full=True):
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
213 - def product(self):
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
220 - def add_initial_state(self, state):
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
229 - def add_intermediate_state(self, state):
230 """ 231 Adds a state to the end of the trajectory 232 233 @param state: state to be added 234 @type state: L{State} 235 """ 236 self._states.append(state.clone())
237
238 - def add_final_state(self, state):
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
247 -class ShortTrajectoryBuilder(TrajectoryBuilder):
248
249 - def add_intermediate_state(self, state):
250 pass
251 252 @property
253 - def product(self):
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
265 -class SimpleProposalCommunicator(object):
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
280 - def __init__(self, proposal):
281 282 self._proposal = proposal
283 284 @property
285 - def proposal(self):
286 return self._proposal
287
288 -class AbstractSingleChainMC(AbstractMC):
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.):
309 310 super(AbstractSingleChainMC, self).__init__(state) 311 312 self._pdf = pdf 313 self._temperature = temperature 314 self._nmoves = 0 315 self._accepted = 0 316 self._last_move_accepted = None
317
318 - def _checkstate(self, state):
319 if not isinstance(state, State): 320 raise TypeError(state)
321
322 - def sample(self):
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
335 - def _propose(self):
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
345 - def _calc_pacc(self, proposal_communicator):
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
379 - def energy(self):
380 """ 381 Log-likelihood of the current state. 382 @rtype: float 383 """ 384 return self._pdf.log_prob(self.state.position)
385 386 @property
387 - def acceptance_rate(self):
388 """ 389 Acceptance rate. 390 """ 391 return float(self._accepted) / float(self._nmoves)
392 393 @property
394 - def last_move_accepted(self):
395 """ 396 Information whether the last MC move was accepted or not. 397 """ 398 return self._last_move_accepted
399 400 @property
401 - def temperature(self):
402 return self._temperature
403
404 -class MCCollection(csb.core.BaseCollectionContainer):
405 """ 406 Collection of single-chain samplers. 407 408 @param items: samplers 409 @type items: list of L{AbstractSingleChainMC} 410 """ 411
412 - def __init__(self, items):
415
416 -class AbstractEnsembleMC(AbstractMC):
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
426 - def __init__(self, samplers):
427 428 self._samplers = MCCollection(samplers) 429 state = EnsembleState([x.state for x in self._samplers]) 430 431 super(AbstractEnsembleMC, self).__init__(state)
432
433 - def sample(self):
434 """ 435 Draw an ensemble sample. 436 437 @rtype: L{EnsembleState} 438 """ 439 440 sample = EnsembleState([sampler.sample() for sampler in self._samplers]) 441 self.state = sample 442 443 return sample
444 445 @property
446 - def energy(self):
447 """ 448 Total ensemble energy. 449 """ 450 return sum([x.energy for x in self._samplers])
451
452 -class AbstractSwapParameterInfo(object):
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
473 - def sampler1(self):
474 return self._sampler1
475 476 @property
477 - def sampler2(self):
478 return self._sampler2
479
480 -class RESwapParameterInfo(AbstractSwapParameterInfo):
481 """ 482 Holds parameters for a standard Replica Exchange swap. 483 """ 484 pass
485
486 -class MDRENSSwapParameterInfo(RESwapParameterInfo):
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):
512 513 super(MDRENSSwapParameterInfo, self).__init__(sampler1, sampler2) 514 515 self._mass_matrix = mass_matrix 516 if self.mass_matrix is None: 517 d = len(sampler1.state.position) 518 self.mass_matrix = csb.numeric.InvertibleMatrix(numpy.eye(d), numpy.eye(d)) 519 520 self._traj_length = traj_length 521 self._gradient = gradient 522 self._timestep = timestep
523 524 @property
525 - def timestep(self):
526 """ 527 Integration timestep. 528 """ 529 return self._timestep
530 @timestep.setter
531 - def timestep(self, value):
532 self._timestep = float(value)
533 534 @property
535 - def traj_length(self):
536 """ 537 Trajectory length in number of integration steps. 538 """ 539 return self._traj_length
540 @traj_length.setter
541 - def traj_length(self, value):
542 self._traj_length = int(value)
543 544 @property
545 - def gradient(self):
546 """ 547 Gradient which governs the equations of motion. 548 """ 549 return self._gradient
550 551 @property
552 - def mass_matrix(self):
553 return self._mass_matrix
554 @mass_matrix.setter
555 - def mass_matrix(self, value):
556 self._mass_matrix = value
557
558 -class ThermostattedMDRENSSwapParameterInfo(MDRENSSwapParameterInfo):
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
593 - def __init__(self, sampler1, sampler2, timestep, traj_length, gradient, mass_matrix=None, 594 temperature=lambda l: 1., collision_probability=0.1, collision_interval=1):
595 596 super(ThermostattedMDRENSSwapParameterInfo, self).__init__(sampler1, sampler2, timestep, 597 traj_length, gradient, 598 mass_matrix=mass_matrix) 599 600 self._collision_probability = None 601 self._collision_interval = None 602 self._temperature = temperature 603 self.collision_probability = collision_probability 604 self.collision_interval = collision_interval
605 606 @property
607 - def collision_probability(self):
608 """ 609 Probability for a collision with the heatbath during one timestep. 610 """ 611 return self._collision_probability
612 @collision_probability.setter
613 - def collision_probability(self, value):
614 self._collision_probability = float(value)
615 616 @property
617 - def collision_interval(self):
618 """ 619 Interval during which collision may occur with probability 620 C{collision_probability}. 621 """ 622 return self._collision_interval
623 @collision_interval.setter
624 - def collision_interval(self, value):
625 self._collision_interval = int(value)
626 627 @property
628 - def temperature(self):
629 return self._temperature
630
631 -class AbstractSwapCommunicator(object):
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):
649 650 self._sampler1 = param_info.sampler1 651 self._sampler2 = param_info.sampler2 652 653 self._traj12 = traj12 654 self._traj21 = traj21 655 656 self._param_info = param_info 657 658 self._acceptance_probability = None 659 self._accepted = False
660 661 @property
662 - def sampler1(self):
663 return self._sampler1
664 665 @property
666 - def sampler2(self):
667 return self._sampler2
668 669 @property
670 - def traj12(self):
671 return self._traj12
672 673 @property
674 - def traj21(self):
675 return self._traj21
676 677 @property
678 - def acceptance_probability(self):
679 return self._acceptance_probability
680 @acceptance_probability.setter
681 - def acceptance_probability(self, value):
682 self._acceptance_probability = value
683 684 @property
685 - def accepted(self):
686 return self._accepted
687 @accepted.setter
688 - def accepted(self, value):
689 self._accepted = value
690 691 @property
692 - def param_info(self):
693 return self._param_info
694
695 -class RESwapCommunicator(AbstractSwapCommunicator):
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
704 -class RENSSwapCommunicator(AbstractSwapCommunicator):
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
714 -class SingleSwapStatistics(object):
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
722 - def __init__(self, param_info):
723 self._total_swaps = 0 724 self._accepted_swaps = 0
725 726 @property
727 - def total_swaps(self):
728 return self._total_swaps
729 730 @property
731 - def accepted_swaps(self):
732 return self._accepted_swaps
733 734 @property
735 - def acceptance_rate(self):
736 """ 737 Acceptance rate of the sampler pair. 738 """ 739 if self.total_swaps > 0: 740 return float(self.accepted_swaps) / float(self.total_swaps) 741 else: 742 return 0.
743
744 - def update(self, accepted):
745 """ 746 Updates swap statistics. 747 """ 748 self._total_swaps += 1 749 self._accepted_swaps += int(accepted)
750
751 -class SwapStatistics(object):
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
760 - def __init__(self, param_infos):
761 self._stats = [SingleSwapStatistics(x) for x in param_infos]
762 763 @property
764 - def stats(self):
765 return tuple(self._stats)
766 767 @property
768 - def acceptance_rates(self):
769 """ 770 Returns acceptance rates for all swaps. 771 """ 772 return [x.acceptance_rate for x in self._stats]
773
774 -class AbstractExchangeMC(AbstractEnsembleMC):
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):
789 super(AbstractExchangeMC, self).__init__(samplers) 790 791 self._swaplist1 = [] 792 self._swaplist2 = [] 793 self._currentswaplist = self._swaplist1 794 795 self._param_infos = param_infos 796 self._statistics = SwapStatistics(self._param_infos)
797
798 - def _checkstate(self, state):
799 if not isinstance(state, EnsembleState): 800 raise TypeError(state)
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
824 - def _propose_swap(self, param_info):
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
836 - def _calc_pacc_swap(self, swapcom):
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
848 - def _accept_swap(self, swapcom):
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
868 - def acceptance_rates(self):
869 """ 870 Return swap acceptance rates. 871 872 @rtype: list of floats 873 """ 874 return self.statistics.acceptance_rates
875 876 @property
877 - def param_infos(self):
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
886 - def statistics(self):
887 return self._statistics
888
889 - def _update_statistics(self, index, accepted):
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
903 -class RENSTrajInfo(object):
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):
918 919 self._param_info = param_info 920 self._protocol = protocol 921 self._init_state = init_state
922 923 @property
924 - def param_info(self):
925 return self._param_info
926 927 @property
928 - def protocol(self):
929 return self._protocol
930 931 @property
932 - def init_state(self):
933 return self._init_state
934
935 -class AbstractRENS(AbstractExchangeMC):
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
945 - def _propose_swap(self, param_info):
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
980 - def _calc_pacc_swap(self, swapcom):
981 982 T1 = swapcom.param_info.sampler1.temperature 983 T2 = swapcom.param_info.sampler2.temperature 984 985 heat12 = swapcom.traj12.heat 986 heat21 = swapcom.traj21.heat 987 988 proposal1 = swapcom.traj21.final 989 proposal2 = swapcom.traj12.final 990 991 state1 = swapcom.traj12.initial 992 state2 = swapcom.traj21.initial 993 994 if swapcom.param_info.mass_matrix.is_unity_multiple: 995 inverse_mass_matrix = 1. / swapcom.param_info.mass_matrix[0][0] 996 else: 997 inverse_mass_matrix = swapcom.param_info.mass_matrix.inverse 998 999 E1 = lambda x:-swapcom.sampler1._pdf.log_prob(x) 1000 E2 = lambda x:-swapcom.sampler2._pdf.log_prob(x) 1001 K = lambda x: 0.5 * numpy.dot(x.T, numpy.dot(inverse_mass_matrix, x)) 1002 1003 w12 = (K(proposal2.momentum) + E2(proposal2.position)) / T2 - \ 1004 (K(state1.momentum) + E1(state1.position)) / T1 - heat12 1005 w21 = (K(proposal1.momentum) + E1(proposal1.position)) / T1 - \ 1006 (K(state2.momentum) + E2(state2.position)) / T2 - heat21 1007 1008 swapcom.acceptance_probability = csb.numeric.exp(-w12 - w21) 1009 1010 return swapcom
1011 1012 @abstractmethod
1013 - def _run_traj_generator(self, traj_info):
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
1026 -class AbstractSwapScheme(object):
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
1037 - def __init__(self, algorithm):
1038 1039 self._algorithm = algorithm
1040 1041 @abstractmethod
1042 - def swap_all(self):
1043 """ 1044 Advises the Replica Exchange-like algorithm to perform swaps according to 1045 the some schedule defined here. 1046 """ 1047 1048 pass
1049
1050 -class AlternatingAdjacentSwapScheme(AbstractSwapScheme):
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
1059 - def __init__(self, algorithm):
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
1068 - def _create_swap_lists(self):
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
1086 - def swap_all(self):
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