nebm fs: polishing integrator; removed base class dep

fidimag/common/chain_method_integrators.py

@@ -42,6 +42,8 @@ def set_step(self, step):
         self._step = step_choices[step]
 
 
+# TODO: these integrators are not general anymore, as they rely on the
+# structure of the chain method classes
 class VerletIntegrator(BaseIntegrator):
     """A quick Verlet integration in Cartesian coordinates
 
@@ -148,19 +150,24 @@ def _step(self, t, y, h, f):
 # NEBM integrator for F-S algorithm
 
 
-class FSIntegrator(BaseIntegrator):
+# TODO: these integrators are not general anymore, as they rely on the
+# structure of the chain method classes
+# TODO: type checks
+class FSIntegrator(object):
     """A step integrator considering the action of the band
     """
-    def __init__(self, band, forces, rhs_fun, action_fun, n_images, n_dofs_image,
+    def __init__(self, ChainObj,
+                 # band, forces, distances, rhs_fun, action_fun,
+                 # n_images, n_dofs_image,
                  maxSteps=1000,
                  maxCreep=5, actionTol=1e-10, forcesTol=1e-6,
                  etaScale=1.0, dEta=2, minEta=0.001,
                  # perturbSeed=42, perturbFactor=0.1,
                  nTrail=10, resetMax=20
                  ):
-        super(FSIntegrator, self).__init__(band, rhs_fun)
+        # super(FSIntegrator, self).__init__(band, rhs_fun)
 
-        self.action_fun = action_fun
+        self.ChainObj = ChainObj
 
         # Integration parameters
         # TODO: move to run function?
@@ -172,24 +179,22 @@ def __init__(self, band, forces, rhs_fun, action_fun, n_images, n_dofs_image,
         self.forcesTol = forcesTol
         self.actionTol = actionTol
         self.maxSteps = maxSteps
-
+        self.step = 0
+        self.nTrail = nTrail
         self.i_step = 0
-        self.n_images = n_images
-        self.n_dofs_image = n_dofs_image
+
+        # Chain objects:
+        self.n_images = self.ChainObj.n_images
+        self.n_dofs_image = self.ChainObj.n_dofs_image
         # self.forces_prev = np.zeros_like(band).reshape(n_images, -1)
         # self.G :
-        self.forces = forces
-        self.forces_old = np.zeros_like(forces)
+        self.forces = self.ChainObj.forces
+        self.distances = self.ChainObj.distances
+        self.forces_old = np.zeros_like(self.ChainObj.forces)
 
-        # Rename y to band
-        self.band = self.y
-        self.band_old = np.zeros_like(self.band)  # y -> band
-
-        self.step = 0
-        self.nTrail = nTrail
-
-    def run_until(self, t):
-        pass
+        # self.band should be just a reference to the band in the ChainObj
+        self.band = self.ChainObj.band
+        self.band_old = np.zeros_like(self.band)
 
     def run_for(self, n_steps):
 
@@ -202,6 +207,7 @@ def run_for(self, n_steps):
         self.trailAction = np.zeros(self.nTrail)
         trailPool = cycle(range(self.nTrail))  # cycle through 0,1,...,(nTrail-1),0,1,...
         eta = 1.0
+        self.i_step = 0
 
         # In __init__:
         # self.band_last[:] = self.band
@@ -216,7 +222,7 @@ def run_for(self, n_steps):
                 self.band[:] = self.band_old
 
                 # Compute from self.band. Do not update the step at this stage:
-                # This step updates the forces in the G array of the nebm module,
+                # This step updates the forces and distances in the G array of the nebm module,
                 # using the current band state self.y
                 # TODO: remove time from chain method rhs
                 #       make a specific function to update G??
@@ -243,12 +249,14 @@ def run_for(self, n_steps):
 
                 self.trailAction
 
-                self.rhs(t, self.y)
-                self.action = self.action_fun()
+                self.ChainObj.nebm_step(self.band, ensure_zero_extrema=True)
+                self.action = self.ChainObj.action_fun()
 
                 self.trailAction[nStart] = self.action
                 nStart = next(trailPool)
 
+                self.i_step += 1
+
                 # Getting averages of forces from the INNER images in the band (no extrema)
                 # (forces are given by vector G in the chain method code)
                 # TODO: we might use all band images, not only inner ones
@@ -266,7 +274,7 @@ def run_for(self, n_steps):
                     exitFlag = True
                     break  # creep loop
 
-                if (n_steps >= self.maxSteps):
+                if (self.i_step >= self.maxSteps):
                     print('Number of steps reached maximum')
                     exitFlag = True
                     break  # creep loop
@@ -282,7 +290,7 @@ def run_for(self, n_steps):
                         print('')
                         resetCount += 1
                         bestAction = self.action_old
-                        RefinePath(self.band)  # Resets the path to equidistant structures (smoothing  kinks?)
+                        self.refine_path(self.ChainObj.distances, self.band)  # Resets the path to equidistant structures (smoothing  kinks?)
                         # PathChanged[:] = True
 
                         if resetCount > self.resetMax:
@@ -309,13 +317,13 @@ def run_for(self, n_steps):
             resetCount = 0
 
     # Taken from the string method class
-    def refine_path(self, distances):
+    def refine_path(self, distances, band):
         """
         """
         new_dist = np.linspace(distances[0], distances[-1], distances.shape[0])
         # Restructure the string by interpolating every spin component
         # print(self.integrator.y[self.n_dofs_image:self.n_dofs_image + 10])
-        bandrs = self.band.reshape(self.n_images, self.n_dofs_image)
+        bandrs = band.reshape(self.n_images, self.n_dofs_image)
         for i in range(self.n_dofs_image):
 
             cs = si.CubicSpline(distances, bandrs[:, i])

fidimag/common/nebm_FS.py

@@ -324,7 +324,7 @@ def compute_min_action(self):
         minAction = np.sum(np.abs(self.energies[1:] - self.energies[:-1]))
         return 2 * (dE + minAction)
 
-    def nebm_step(self, y):
+    def nebm_step(self, y, ensure_zero_extrema=False):
 
         self.compute_effective_field_and_energy(y)
         nebm_clib.project_images(self.gradientE, y,
@@ -342,6 +342,12 @@ def nebm_step(self, y):
                                           self.n_dofs_image
                                           )
 
+        # The effective force at the extreme images should already be zero, but
+        # we will manually remove any value
+        if ensure_zero_extrema:
+            self.G[:self.n_dofs_image] = 0
+            self.G[-self.n_dofs_image:] = 0
+
         # Should be the same if we project the gradient before, instead
         # of the total force
         # nebm_clib.project_images(self.G, y,
@@ -388,7 +394,7 @@ def compute_distances(self):
     #     else:
     #         return 0
 
-    def step_RHS(self, t, y):
+    def step_RHS(self, band):
         """
 
         This function is called on every iteration of the integrators in
@@ -400,7 +406,7 @@ def step_RHS(self, t, y):
 
         # Update the effective field, energies, spring forces and tangents
         # using the *y* array
-        self.nebm_step(y)
+        self.nebm_step(band)
 
         # Now set the RHS of the equation as the effective force on the energy
         # band, which is stored on the self.G array