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

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