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