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
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
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
74
75 @property
78 @gradient.setter
80 self._gradient = value
81
82 @property
85 @timestep.setter
87 self._timestep = float(value)
88
89 @property
91 return self._mass_matrix
92 @mass_matrix.setter
94 self._mass_matrix = value
95
96 - def generate(self, init_state, length, return_trajectory=False):
105
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
114
115 self._items = items
116 self._n_items = len(self._items)
117 self._current = 0
118
122
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
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
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
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):
317
318 @abstractmethod
320 """
321 Initializes the sampler with which to obtain the MC state
322 trajectory.
323 """
324
325 pass
326
327 @property
329 """
330 Acceptance rate of the MC sampler that generated the
331 trajectory.
332 """
333 return self._acceptance_rate
334
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.):
362
369
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.):
410
419
420 @property
422 return self._mass_matrix
423 @mass_matrix.setter
425 self._mass_matrix = value
426