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