Package csb :: Package numeric :: Module integrators

[frames] | no frames]

# Source Code for Module csb.numeric.integrators

      1  """
      2  provides various integration schemes and an abstract gradient class.
      3  """
      4
      5  import numpy
      6
      7  from abc import ABCMeta, abstractmethod
      8  from csb.statistics.samplers.mc import State, TrajectoryBuilder
      9  from csb.numeric import InvertibleMatrix
     10
     11
     12 -class AbstractIntegrator(object):
     13      """
     14      Abstract integrator class. Subclasses implement different integration
     15      schemes for solving deterministic equations of motion.
     16
     17      @param timestep: Integration timestep
     18      @type timestep: float
     19
     20      @param gradient: Gradient of potential energy
     21      @type gradient: L{AbstractGradient}
     22      """
     23
     24      __metaclass__ = ABCMeta
     25
     26 -    def __init__(self, timestep, gradient):
     27
     28          self._timestep = timestep
     29          self._gradient = gradient
     30
     31 -    def integrate(self, init_state, length, mass_matrix=None, return_trajectory=False):
     32          """
     33          Integrates equations of motion starting from an initial state a certain
     34          number of steps.
     35
     36          @param init_state: Initial state from which to start integration
     37          @type init_state: L{State}
     38
     39          @param length: Nubmer of integration steps to be performed
     40          @type length: int
     41
     42          @param mass_matrix: Mass matrix
     43          @type mass_matrix:  n-dimensional L{InvertibleMatrix} with n being the dimension
     44                                      of the configuration space, that is, the dimension of
     45                                      the position / momentum vectors
     46
     47          @param return_trajectory: Return complete L{Trajectory} instead of the initial
     48                                    and final states only (L{PropagationResult}). This reduces
     49                                    performance.
     50          @type return_trajectory: boolean
     51
     52          @rtype: L{AbstractPropagationResult}
     53          """
     54
     55          builder = TrajectoryBuilder.create(full=return_trajectory)
     56
     57          builder.add_initial_state(init_state)
     58          state = init_state.clone()
     59
     60          for i in range(length - 1):
     61              state = self.integrate_once(state, i, mass_matrix=mass_matrix)
     62              builder.add_intermediate_state(state)
     63
     64          state = self.integrate_once(state, length - 1, mass_matrix=mass_matrix)
     65          builder.add_final_state(state)
     66
     67          return builder.product
     68
     69      @abstractmethod
     70 -    def integrate_once(self, state, current_step, mass_matrix=None):
     71          """
     72          Integrates one step starting from an initial state and an initial time
     73          given by the product of the timestep and the current_step parameter.
     74          The input C{state} is changed in place.
     75
     76          @param state: State which to evolve one integration step
     77          @type state: L{State}
     78
     79          @param current_step: Current integration step
     80          @type current_step: int
     81
     82          @param mass_matrix: mass matrix
     83          @type mass_matrix:  n-dimensional numpy array with n being the dimension
     84                              of the configuration space, that is, the dimension of
     85                              the position / momentum vectors
     86          @return: the altered state
     87          @rtype: L{State}
     88          """
     89          pass
     90
     91 -    def _get_inverse(self, mass_matrix):
     92
     93          inverse_mass_matrix = None
     94          if mass_matrix is None:
     95              inverse_mass_matrix = 1.0
     96          else:
     97              if mass_matrix.is_unity_multiple:
     98                  inverse_mass_matrix = mass_matrix.inverse[0][0]
     99              else:
    100                  inverse_mass_matrix = mass_matrix.inverse
    101
    102          return inverse_mass_matrix
    103
    104 -class LeapFrog(AbstractIntegrator):
    105      """
    106      Leap Frog integration scheme implementation that calculates position and
    107      momenta at equal times. Slower than FastLeapFrog, but intermediate points
    108      in trajectories obtained using
    109      LeapFrog.integrate(init_state, length, return_trajectoy=True) are physical.
    110      """
    111
    112 -    def integrate_once(self, state, current_step, mass_matrix=None):
    113
    114          inverse_mass_matrix = self._get_inverse(mass_matrix)
    115
    116          i = current_step
    117
    118          if i == 0:
    119              self._oldgrad = self._gradient(state.position, 0.)
    120
    121          momentumhalf = state.momentum - 0.5 * self._timestep * self._oldgrad
    122          state.position = state.position + self._timestep * numpy.dot(inverse_mass_matrix, momentumhalf)
    123          self._oldgrad = self._gradient(state.position, (i + 1) * self._timestep)
    124          state.momentum = momentumhalf - 0.5 * self._timestep * self._oldgrad
    125
    126          return state
    127
    128 -class FastLeapFrog(LeapFrog):
    129      """
    130      Leap Frog integration scheme implementation that calculates position and
    131      momenta at unequal times by concatenating the momentum updates of two
    132      successive integration steps.
    133      WARNING: intermediate points in trajectories obtained by
    134      FastLeapFrog.integrate(init_state, length, return_trajectories=True)
    135      are NOT to be interpreted as phase-space trajectories, because
    136      position and momenta are not given at equal times! In the initial and the
    137      final state, positions and momenta are given at equal times.
    138      """
    139
    140 -    def integrate(self, init_state, length, mass_matrix=None, return_trajectory=False):
    141
    142          inverse_mass_matrix = self._get_inverse(mass_matrix)
    143
    144          builder = TrajectoryBuilder.create(full=return_trajectory)
    145
    146          builder.add_initial_state(init_state)
    147          state = init_state.clone()
    148
    149          state.momentum = state.momentum - 0.5 * self._timestep * self._gradient(state.position, 0.)
    150
    151          for i in range(length-1):
    152              state.position = state.position + self._timestep * numpy.dot(inverse_mass_matrix, state.momentum)
    153              state.momentum = state.momentum - self._timestep * \
    154                               self._gradient(state.position, (i + 1) * self._timestep)
    155              builder.add_intermediate_state(state)
    156
    157          state.position = state.position + self._timestep * numpy.dot(inverse_mass_matrix, state.momentum)
    158          state.momentum = state.momentum - 0.5 * self._timestep * \
    159                           self._gradient(state.position, length * self._timestep)
    160          builder.add_final_state(state)
    161
    162          return builder.product
    163
    164 -class VelocityVerlet(AbstractIntegrator):
    165      """
    166      Velocity Verlet integration scheme implementation.
    167      """
    168
    169 -    def integrate_once(self, state, current_step, mass_matrix=None):
    170
    171          inverse_mass_matrix = self._get_inverse(mass_matrix)
    172
    173          i = current_step
    174
    175          if i == 0:
    176              self._oldgrad = self._gradient(state.position, 0.)
    177
    178          state.position = state.position + self._timestep * numpy.dot(inverse_mass_matrix, state.momentum) \
    179                           - 0.5 * self._timestep ** 2 * numpy.dot(inverse_mass_matrix, self._oldgrad)
    180          newgrad = self._gradient(state.position, (i + 1) * self._timestep)
    181          state.momentum = state.momentum - 0.5 * self._timestep * (self._oldgrad + newgrad)
    182          self._oldgrad = newgrad
    183
    184          return state
    185
    186 -class AbstractGradient(object):
    187      """
    188      Abstract gradient class. Implementations evaluate the gradient of an energy
    189      function.
    190      """
    191
    192      __metaclass__ = ABCMeta
    193
    194      @abstractmethod
    195 -    def evaluate(self, q, t):
    196          """
    197          Evaluates the gradient at position q and time t.
    198
    199          @param q: Position array
    200          @type q:  One-dimensional numpy array
    201
    202          @param t: Time
    203          @type t: float
    204
    205          @rtype: numpy array
    206          """
    207          pass
    208
    209 -    def __call__(self, q, t):
    210          """
    211          Evaluates the gradient at position q and time t.
    212
    213          @param q: Position array
    214          @type q:  One-dimensional numpy array
    215
    216          @param t: Time
    217          @type t: float
    218
    219          @rtype: numpy array
    220          """
    221          State.check_flat_array(q)
    222          return self.evaluate(q, t)
    223