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

# Source Code for Module csb.statistics.samplers.mc.propagators

      1  """
      2  Provides various deterministic and stochastic propagators.
      3  """
      4
      5  import numpy
      6
      7  from abc import ABCMeta, abstractmethod
      8
      9  from csb.statistics.samplers.mc import TrajectoryBuilder
     10  from csb.numeric.integrators import FastLeapFrog, VelocityVerlet
     11  from csb.numeric import InvertibleMatrix
     12
     13 -class AbstractPropagator(object):
     14      """
     15      Abstract propagator class. Subclasses serve to propagate
     16      an inital state by some dynamics to a final state.
     17      """
     18
     19      __metaclass__ = ABCMeta
     20
     21      @abstractmethod
     22 -    def generate(self, init_state, length, return_trajectory=False):
     23          """
     24          Generate a trajectory, starting from an initial state with a certain length.
     25
     26          @param init_state: Initial state from which to propagate
     27          @type init_state: L{State}
     28
     29          @param length: Length of the trajectory (in integration steps or stochastic moves)
     30          @type length: int
     31
     32          @param return_trajectory: Return complete L{Trajectory} instead of the initial
     33                                    and final states only (L{PropagationResult})
     34          @type return_trajectory: boolean
     35
     36          @rtype: L{AbstractPropagationResult}
     37          """
     38          pass
     39
     40 -class MDPropagator(AbstractPropagator):
     41      """
     42      Molecular Dynamics propagator. Generates a trajectory
     43      by integration of Hamiltionian equations of motion.
     44
     45      @param gradient: Gradient of potential energy. Guides the dynamics.
     46      @type gradient: L{AbstractGradient}
     47
     48      @param timestep: Timestep to be used for integration
     49      @type timestep: float
     50
     51      @param mass_matrix: Mass matrix
     52      @type mass_matrix: n-dimensional L{InvertibleMatrix} with n being the dimension
     53                                 of the configuration space, that is, the dimension of
     54                                 the position / momentum vectors
     55
     56      @param integrator: Subclass of L{AbstractIntegrator} to be used to integrate
     57                         Hamiltonian equations of motion
     58      @type integrator: type
     59      """
     60
     61 -    def __init__(self, gradient, timestep, mass_matrix=None, integrator=FastLeapFrog):
     62
     63          super(MDPropagator, self).__init__()
     64
     65          self._gradient = None
     66          self.gradient = gradient
     67          self._mass_matrix = mass_matrix
     68          self._timestep = None
     69          self.timestep = timestep
     70          self._integrator = integrator
     71
     72          self._first_run = True
     73
     74      @property
     75 -    def gradient(self):
     76          return self._gradient
     77      @gradient.setter
     78 -    def gradient(self, value):
     79          self._gradient = value
     80
     81      @property
     82 -    def timestep(self):
     83          return self._timestep
     84      @timestep.setter
     85 -    def timestep(self, value):
     86          self._timestep = float(value)
     87
     88      @property
     89 -    def mass_matrix(self):
     90          return self._mass_matrix
     91      @mass_matrix.setter
     92 -    def mass_matrix(self, value):
     93          self._mass_matrix = value
     94
     95 -    def generate(self, init_state, length, return_trajectory=False):
     96
     97          integrator = self._integrator(self.timestep, self.gradient)
     98
     99          result = integrator.integrate(init_state, length,
    100                                        mass_matrix=self.mass_matrix,
    101                                        return_trajectory=return_trajectory)
    102
    103          return result
    104
    105 -class Looper(object):
    106      """
    107      Implements an iterable list with a ring-like topology,
    108      that is, if the iterator points on the last element,
    109      next() returns the first element.
    110      """
    111
    112 -    def __init__(self, items):
    113
    114          self._items = items
    115          self._n_items = len(self._items)
    116          self._current = 0
    117
    118 -    def __iter__(self):
    119
    120          return self
    121
    122 -    def next(self):
    123
    124          if self._current == self._n_items:
    125              self._current = 0
    126
    127          self._current += 1
    128
    129          return self._items[self._current - 1]
    130
    131
    132 -class ThermostattedMDPropagator(MDPropagator):
    133      """
    134      Thermostatted Molecular Dynamics propagator. Employs the Andersen thermostat
    135      which simulates collision with particles of a heat bath at a given temperature.
    136
    137      @param gradient: Gradient of potential energy. Guides the dynamics.
    138      @type gradient: L{AbstractGradient}
    139
    140      @param timestep: Timestep to be used for integration
    141      @type timestep: float
    142
    143      @param mass_matrix: Mass matrix
    144      @type mass_matrix: n-dimensional L{InvertibleMatrix} with n being the dimension
    145                                 of the configuration space, that is, the dimension of
    146                                 the position / momentum vectors
    147
    148      @param temperature: Time-dependent temperature
    149      @type temperature: Real-valued function
    150
    151      @param collision_probability: collision probability within duration of one timestep
    152      @type collision_probability: float
    153
    154      @param update_interval: Interval with which momenta are redrawn
    155      @type update_interval: int
    156
    157      @param integrator: Subclass of L{AbstractIntegrator} to be used to perform
    158                         integration steps between momentum updates
    159      @type integrator: type
    160      """
    161
    162 -    def __init__(self, gradient, timestep, mass_matrix=None, temperature=lambda t: 1.,
    163                   collision_probability=0.1, update_interval=1, integrator=VelocityVerlet):
    164
    165          super(ThermostattedMDPropagator, self).__init__(gradient, timestep,
    166                                                          mass_matrix, integrator)
    167
    168          self._collision_probability = collision_probability
    169          self._update_interval = update_interval
    170          self._temperature = temperature
    171
    172 -    def _update(self, momentum, T, collision_probability):
    173          """
    174          Simulate collision with heat bath particles.
    175
    176          @param momentum: Momentum
    177          @type momentum: one-dimensional numpy array of numbers
    178
    179          @param T: Temperature of the heat bath
    180          @type T: float
    181
    182          @param collision_probability: collision probability within duration of one timestep
    183          @type collision_probability: float
    184
    185          @rtype: tuple (updated momentum, heat induced by the update)
    186          """
    187
    188          d = len(momentum)
    189
    190          heat = 0.
    191          update_list = numpy.where(numpy.random.random(d) < collision_probability)[0]
    192
    193          if len(update_list) > 0:
    194              K = None
    195              if self.mass_matrix.is_unity_multiple:
    196                  K = lambda x: 0.5 * sum(x ** 2) / self.mass_matrix[0][0]
    197              else:
    198                  K = lambda x: 0.5 * numpy.dot(x.T, numpy.dot(self.mass_matrix.inverse, x))
    199
    200              ke_old = K(momentum)
    201
    202              updated_momentum = [numpy.sqrt(T) * self._random_loopers[i].next() for i in update_list]
    203              momentum[update_list] = updated_momentum
    204              heat = (K(momentum) - ke_old) / T
    205
    206          return momentum, heat
    207
    208 -    def _step(self, i, state, heat, integrator):
    209          """
    210          Performs one step consisting of an integration step
    211          and possibly a momentum update
    212
    213          @param i: integration step count
    214          @type i: int
    215
    216          @param state: state to be updated
    217          @type state: L{State}
    218
    219          @param heat: heat produced up to the current integration step
    220          @type heat: float
    221
    222          @param integrator: integration scheme used to evolve the state deterministically
    223          @type integrator: L{AbstractIntegrator}
    224          """
    225
    226          state = integrator.integrate_once(state, i, mass_matrix=self.mass_matrix)
    227
    228          if i % self._update_interval == 0:
    229              state.momentum, stepheat = self._update(state.momentum,
    230                                                      self._temperature(i * self.timestep),
    231                                                      self._collision_probability)
    232
    233              heat += stepheat
    234
    235          return state, heat
    236
    237 -    def generate(self, init_state, length, return_trajectory=False):
    238
    239          if self._first_run == True and self.mass_matrix is None:
    240              d = len(init_state.position)
    241              self.mass_matrix = InvertibleMatrix(numpy.eye(d), numpy.eye(d))
    242
    243          integrator = self._integrator(self.timestep, self.gradient)
    244          builder = TrajectoryBuilder.create(full=return_trajectory)
    245
    246          builder.add_initial_state(init_state)
    247
    248          heat = 0.
    249          state = init_state.clone()
    250
    251          d = len(state.position)
    252
    253          n_randoms = int(1.5 * length * self._collision_probability / float(self._update_interval))
    254
    255          if n_randoms < 5:
    256              n_randoms = 5
    257
    258          if not self.mass_matrix.is_unity_multiple:
    259              randoms = numpy.random.multivariate_normal(mean=numpy.zeros(d),
    260                                                         cov=self.mass_matrix,
    261                                                         size=n_randoms).T
    262          else:
    263              randoms = numpy.random.normal(scale=numpy.sqrt(self.mass_matrix[0][0]),
    264                                            size=(d, n_randoms))
    265          self._random_loopers = [Looper(x) for x in randoms]
    266
    267          for i in range(length - 1):
    268              state, heat = self._step(i, state, heat, integrator)
    269              builder.add_intermediate_state(state)
    270
    271          state, heat = self._step(length - 1, state, heat, integrator)
    272          builder.add_final_state(state)
    273
    274          traj = builder.product
    275          traj.heat = heat
    276
    277          return traj
    278
    279 -class AbstractMCPropagator(AbstractPropagator):
    280      """
    281      Provides the interface for MC trajectory generators. Implementations
    282      generate a sequence of states according to some implementation of
    283      L{AbstractSingleChainMC}.
    284
    285      @param pdf: PDF to sample from
    286      @type pdf: L{AbstractDensity}
    287
    288      @param temperature: See documentation of L{AbstractSingleChainMC}
    289      @type temperature: float
    290      """
    291
    292      __metaclass__ = ABCMeta
    293
    294 -    def __init__(self, pdf, temperature=1.):
    295
    296          self._pdf = pdf
    297          self._temperature = temperature
    298          self._acceptance_rate = 0.0
    299
    300 -    def generate(self, init_state, length, return_trajectory=True):
    301
    302          self._init_sampler(init_state)
    303          self._sampler.state = init_state
    304
    305          builder = TrajectoryBuilder.create(full=return_trajectory)
    306
    307          builder.add_initial_state(init_state)
    308
    309          for i in range(length):
    310              self._sampler.sample()
    311              if i != length - 1:
    312                  builder.add_intermediate_state(self._sampler.state)
    313
    314          builder.add_final_state(self._sampler.state)
    315
    316          self._acceptance_rate = self._sampler.acceptance_rate
    317
    318          return builder.product
    319
    320      @abstractmethod
    321 -    def _init_sampler(self, init_state):
    322          """
    323          Initializes the sampler with which to obtain the MC state
    324          trajectory.
    325          """
    326
    327          pass
    328
    329      @property
    330 -    def acceptance_rate(self):
    331          """
    332          Acceptance rate of the MC sampler that generated the
    333          trajectory.
    334          """
    335          return self._acceptance_rate
    336
    337
    338 -class RWMCPropagator(AbstractMCPropagator):
    339      """
    340      Draws a number of samples from a PDF using the L{RWMCSampler} and
    341      returns them as a L{Trajectory}.
    342
    343      @param pdf: PDF to sample from
    344      @type pdf: L{AbstractDensity}
    345      @param stepsize: Serves to set the step size in
    346                       proposal_density, e.g. for automatic acceptance
    347                       rate adaption
    348      @type stepsize: float
    349      @param proposal_density: The proposal density as a function f(x, s)
    350                               of the current state x and the stepsize s.
    351                               By default, the proposal density is uniform,
    352                               centered around x, and has width s.
    353      @type proposal_density: callable
    354
    355      @param temperature: See documentation of L{AbstractSingleChainMC}
    356      @type temperature: float
    357      """
    358
    359 -    def __init__(self, pdf, stepsize=1., proposal_density=None, temperature=1.):
    360
    361          super(RWMCPropagator, self).__init__(pdf, temperature)
    362
    363          self._stepsize = stepsize
    364          self._proposal_density = proposal_density
    365
    366 -    def _init_sampler(self, init_state):
    367
    368          from csb.statistics.samplers.mc.singlechain import RWMCSampler
    369
    370          self._sampler = RWMCSampler(self._pdf, init_state, self._stepsize,
    371                                      self._proposal_density, self._temperature)
    372
    373 -class HMCPropagator(AbstractMCPropagator):
    374      """
    375      Draws a number of samples from a PDF using the L{HMCSampler} and
    376      returns them as a L{Trajectory}.
    377
    378      @param pdf: PDF to sample from
    379      @type pdf: L{AbstractDensity}
    380      @param gradient: Gradient of the negative log-probability
    381      @type gradient: L{AbstractGradient}
    382
    383      @param timestep: Timestep used for integration
    384      @type timestep: float
    385
    386      @param nsteps: Number of integration steps to be performed in
    387                     each iteration
    388      @type nsteps: int
    389
    390      @param mass_matrix: Mass matrix
    391      @type mass_matrix: n-dimensional L{InvertibleMatrix} with n being the dimension
    392                                 of the configuration space, that is, the dimension of
    393                                 the position / momentum vectors
    394
    395      @param integrator: Subclass of L{AbstractIntegrator} to be used for
    396                         integrating Hamiltionian equations of motion
    397      @type integrator: type
    398
    399      @param temperature: See documentation of L{AbstractSingleChainMC}
    400      @type temperature: float
    401      """
    402
    403 -    def __init__(self, pdf, gradient, timestep, nsteps, mass_matrix=None,
    404                   integrator=FastLeapFrog, temperature=1.):
    405
    406          super(HMCPropagator, self).__init__(pdf, temperature)
    407
    408          self._gradient = gradient
    409          self._timestep = timestep
    410          self._nsteps = nsteps
    411          self._mass_matrix = mass_matrix
    412          self._integrator = integrator
    413
    414 -    def _init_sampler(self, init_state):
    415
    416          from csb.statistics.samplers.mc.singlechain import HMCSampler
    417
    418          self._sampler = HMCSampler(self._pdf, init_state, self._gradient,
    419                                     self._timestep, self._nsteps,
    420                                     mass_matrix=self.mass_matrix,
    421                                     integrator=self._integrator, temperature=self._temperature)
    422
    423      @property
    424 -    def mass_matrix(self):
    425          return self._mass_matrix
    426      @mass_matrix.setter
    427 -    def mass_matrix(self, value):
    428          self._mass_matrix = value
    429
    430
    431 -class AbstractNCMCPropagator(AbstractMCPropagator):
    432      """
    433      Draws a number of samples from a PDF using the L{AbstractNCMCSampler}.
    434
    435      @param protocol: The nonequilibrium protocol specifying a sequence of
    436                       perturbation and propagation steps
    437      @type protocol: L{Protocol}
    438
    439      @param reverse_protocol: The protocol with the order of perturbation and
    440                               propagation reversed in each step.
    441      @type reverse_protocol: L{Protocol}
    442      """
    443
    444      __metaclass__ = ABCMeta
    445
    446 -    def __init__(self, protocol, reverse_protocol):
    447
    448          self._protocol = None
    449          self.protocol = protocol
    450          self._reverse_protocol = None
    451          self.reverse_protocol = reverse_protocol
    452
    453          pdf = self.protocol.steps[0].perturbation.sys_before.hamiltonian
    454          temperature = self.protocol.steps[0].perturbation.sys_before.hamiltonian.temperature
    455
    456          super(AbstractNCMCPropagator, self).__init__(pdf, temperature)
    457
    458      @property
    459 -    def protocol(self):
    460          return self._protocol
    461      @protocol.setter
    462 -    def protocol(self, value):
    463          self._protocol = value
    464
    465      @property
    466 -    def reverse_protocol(self):
    467          return self._reverse_protocol
    468      @reverse_protocol.setter
    469 -    def reverse_protocol(self, value):
    470          self._reverse_protocol = value
    471