Merge branch 'nebm-kf' of github.com:computationalmodelling/fidimag into nebm-kf

Author: davidcortesortuno

.readthedocs.yml

@@ -0,0 +1,13 @@
+version: 2
+
+build:
+  os: "ubuntu-20.04"
+  tools:
+    python: "mambaforge-22.9"
+
+conda:
+  environment: doc/environment.yml
+
+sphinx:
+  # Path to your Sphinx configuration file.
+  configuration: doc/conf.py

doc/conf.py

@@ -78,7 +78,7 @@
 
 html_theme_options = {
     "github_url": "https://github.com/computationalmodelling/fidimag",
-    "use_edit_page_button": True,
+    # "use_edit_page_button": True,
     "show_toc_level": 1,
     # "navbar_align": "right",  # For testing that the navbar items align properly
 }

fidimag/common/chain_method_base.py

@@ -199,10 +199,7 @@ def __init__(self, sim,
         self.n_dofs_image_material = np.sum(self._material)
 
         # VTK saver for the magnetisation/spin field --------------------------
-        self.VTK = VTK(self.mesh,
-                       directory='vtks'.format(self.name),
-                       filename='image'
-                       )
+        self.VTK = VTK(self.mesh, directory='vtks', filename='image')
 
         # Functions to convert the energy band coordinates to Cartesian
         # coordinates when saving VTK and NPY files We assume Cartesian
@@ -303,6 +300,9 @@ def __init__(self, sim,
 
         self.G_log = []
 
+        # Make sure the sim object does not have a driver: (see sim_base.py -> set_m)
+        sim.driver = None
+
     # TODO: Move this property to the NEBM classes because they are only
     # relevant to the NEBM and not the string method
     @property
@@ -403,21 +403,17 @@ def save_npys(self, coordinates_function=None):
                 np.save(name, self.band[i])
         self.band.shape = (-1)
 
-    def initialise_integrator(self,
-                              integrator='sundials',
-                              rtol=1e-6, atol=1e-6):
+    def initialise_integrator(self, integrator='sundials', rtol=1e-6, atol=1e-6):
         self.t = 0
         self.iterations = 0
         self.ode_count = 1
 
         if integrator == 'sundials':
             if not self.openmp:
-                self.integrator = cvode.CvodeSolver(self.band,
-                                                    self.Sundials_RHS)
+                self.integrator = cvode.CvodeSolver(self.band, self.Sundials_RHS)
                 self.integrator.set_options(rtol, atol)
             else:
-                self.integrator = cvode.CvodeSolver_OpenMP(self.band,
-                                                           self.Sundials_RHS)
+                self.integrator = cvode.CvodeSolver_OpenMP(self.band, self.Sundials_RHS)
                 self.integrator.set_options(rtol, atol)
         # elif integrator == 'scipy':
         #     self.integrator = ScipyIntegrator(self.band, self.step_RHS)
@@ -440,8 +436,7 @@ def initialise_integrator(self,
             # In Verlet algorithm we only use the total force G and not YxYxG:
             self._llg_evolve = False
         else:
-            raise Exception('No valid integrator specified. Available: '
-                            '"sundials", "scipy"')
+            raise Exception('No valid integrator specified. Available: "sundials", "scipy"')
 
     def create_tablewriter(self):
         entities_energy = {
@@ -763,11 +758,15 @@ def compute_polynomial_factors(self, compute_fields=True):
             self.compute_tangents(self.band)
             self.compute_distances()
 
+        self.gradientE.shape = (self.n_images, -1)
+        self.tangents.shape = (self.n_images, -1)
+
         deltas = np.zeros(self.n_images)
         for i in range(self.n_images):
-            deltas[i] = np.dot(self.scale * (self.gradientE).reshape(self.n_images, -1)[i],
-                               self.tangents.reshape(self.n_images, -1)[i]
-                               )
+            deltas[i] = np.dot(self.scale * self.gradientE[i], self.tangents[i])
+
+        self.gradientE.shape = (-1)
+        self.tangents.shape = (-1)
 
         ds = self.path_distances
         E = self.energies

fidimag/common/chain_method_integrators.py

@@ -160,8 +160,8 @@ def __init__(self, ChainObj,
                  # band, forces, distances, rhs_fun, action_fun,
                  # n_images, n_dofs_image,
                  maxSteps=1000,
-                 maxCreep=6, actionTol=1e-2, forcesTol=1e-6,
-                 etaScale=1e-6, dEta=4., minEta=1e-6,
+                 maxCreep=4, actionTol=1e-2, forcesTol=1e-8,
+                 etaScale=1e4, dEta=4., minEta=1e-6,
                  # perturbSeed=42, perturbFactor=0.1,
                  nTrail=13, resetMax=20
                  ):
@@ -188,13 +188,17 @@ def __init__(self, ChainObj,
         self.n_dofs_image = self.ChainObj.n_dofs_image
         # self.forces_prev = np.zeros_like(band).reshape(n_images, -1)
         # self.G :
-        self.forces = self.ChainObj.G
+        self.forces = -self.ChainObj.G
         self.distances = self.ChainObj.distances
         self.forces_old = np.zeros_like(self.ChainObj.G)
 
         # self.band should be just a reference to the band in the ChainObj
-        self.band = self.ChainObj.band
-        self.band_old = np.copy(self.band)
+        # self.band = self.ChainObj.band
+        self.band_old = np.copy(self.ChainObj.band)
+
+        # CHECK
+        self.ChainObj.k[:] = 0.0
+        self.ChainObj.spring_force[:] = 0.0
 
     def run_for(self, n_steps):
 
@@ -218,43 +222,67 @@ def run_for(self, n_steps):
         self.ChainObj.tablewriter_dm.save()
 
         np.save(self.ChainObj.name + '_init.npy', self.ChainObj.band)
+        self.action_old = 0.0
 
         while not exitFlag:
 
             if totalRestart:
                 if self.i_step > 0:
                     print('Restarting')
-                self.band[:] = self.band_old
+                self.ChainObj.band[:] = self.band_old
 
                 # Compute from self.band. Do not update the step at this stage:
                 # This step updates the forces and distances in the G array of the nebm module,
                 # using the current band state self.y
                 # TODO: make a specific function to update G??
                 print('Computing forces')
-                self.ChainObj.nebm_step(self.band, ensure_zero_extrema=True)
+                self.ChainObj.nebm_step(self.ChainObj.band, ensure_zero_extrema=True)
                 self.action = self.ChainObj.compute_action()
 
                 # self.step += 1
-                self.forces_old[:] = self.forces  # Scale field??
+                self.forces_old[:] = self.ChainObj.G  # Scale field??
                 # self.gradE_last[~_material] = 0.0
                 # self.gradE[:] = self.gradE_last
                 self.action_old = self.action
                 self.trailAction[nStart] = self.action
                 nStart = next(trailPool)
                 eta = 1.0
                 totalRestart = False
-
             creepCount = 0
 
+
             # Creep stage: minimise with a fixed eta
             while creepCount < self.maxCreep:
                 # Update spin. Avoid pinned or zero-Ms sites
-                self.band[:] = self.band_old - eta * self.etaScale * self.forces_old
-                normalise_spins(self.band)
+                self.ChainObj.band[:] = self.band_old - eta * self.etaScale * self.forces_old
+                normalise_spins(self.ChainObj.band)
+
+                # self.refine_path(self.ChainObj.path_distances, self.ChainObj.band)  # Resets the path to equidistant structures (smoothing  kinks?)
+                # normalise_spins(self.ChainObj.band)
+
+                self.trailAction[nStart] = self.action
+                nStart = next(trailPool)
 
-                self.trailAction
+                self.ChainObj.nebm_step(self.ChainObj.band, ensure_zero_extrema=True)
+
+                # gradE = self.ChainObj.gradientE.reshape(self.n_images, -1)
+                # band = self.ChainObj.band.reshape(self.n_images, -1)
+                # tgts = self.ChainObj.tangents.reshape(self.n_images, -1)
+                # forces = self.ChainObj.G.reshape(self.n_images, -1)
+                # for im_idx in range(gradE.shape[0]):
+                #         gradE_dot_tgt = np.sum(forces[im_idx] * tgts[im_idx])
+                #         print(f'Im {im_idx}', forces[im_idx], tgts[im_idx], gradE_dot_tgt)
+                #         print('-----', np.linalg.norm(tgts[im_idx]))
+                # gradE = self.ChainObj.gradientE.reshape(self.n_images, -1)
+                # band = self.ChainObj.band.reshape(self.n_images, -1)
+                # tgts = self.ChainObj.tangents.reshape(self.n_images, -1)
+                # forces = self.ChainObj.G.reshape(self.n_images, -1)
+
+                # for im_idx in range(gradE.shape[0]):
+                #     gradE_dot_tgt = np.sum(forces[im_idx] * tgts[im_idx])
+                #     print(np.linalg.norm(tgts[im_idx]), gradE_dot_tgt)
+                # print('--')
 
-                self.ChainObj.nebm_step(self.band, ensure_zero_extrema=True)
                 self.action = self.ChainObj.compute_action()
 
                 self.trailAction[nStart] = self.action
@@ -269,9 +297,9 @@ def run_for(self, n_steps):
                 # 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, although G is zero at the extrema
-                Gnorms2 = np.sum(self.forces[INNER_DOFS].reshape(-1, 3)**2, axis=1)
+                Gnorms2 = np.sum(self.ChainObj.G.reshape(-1, 3)**2, axis=1)
                 # Compute the root mean square per image
-                rms_G_norms_per_image = np.sum(Gnorms2.reshape(self.n_images - 2, -1), axis=1) / self.ChainObj.n_spins
+                rms_G_norms_per_image = np.sum(Gnorms2.reshape(self.n_images, -1), axis=1) / self.ChainObj.n_images
                 rms_G_norms_per_image = np.sqrt(rms_G_norms_per_image)
                 mean_rms_G_norms_per_image = np.mean(rms_G_norms_per_image)
 
@@ -314,7 +342,7 @@ def run_for(self, n_steps):
                         resetCount += 1
                         # bestAction = self.action_old
                         print('Refining path')
-                        self.refine_path(self.ChainObj.path_distances, self.band)  # Resets the path to equidistant structures (smoothing  kinks?)
+                        self.refine_path(self.ChainObj.path_distances, self.ChainObj.band)  # Resets the path to equidistant structures (smoothing  kinks?)
                         # PathChanged[:] = True
 
                         if resetCount > self.resetMax:
@@ -328,14 +356,14 @@ def run_for(self, n_steps):
                     # Otherwise, just start again with smaller alpha 
                     else:
                         print('Decreasing alpha')
-                        self.band[:] = self.band_old
-                        self.forces[:] = self.forces_old
+                        self.ChainObj.band[:] = self.band_old
+                        self.ChainObj.G[:] = self.forces_old
                 # If action decreases, move to next creep step
                 else:
                     creepCount += 1
                     self.action_old = self.action
-                    self.band_old[:] = self.band
-                    self.forces_old[:] = self.forces
+                    self.band_old[:] = self.ChainObj.band
+                    self.forces_old[:] = self.ChainObj.G
 
                     if (mean_rms_G_norms_per_image < self.forcesTol):
                         print('Change of mean of the force RMSquares negligible')

fidimag/common/nebm_FS.py

@@ -126,6 +126,8 @@ def __init__(self, sim,
                                       openmp=openmp
                                       )
 
+        # The gradient = - H_eff
+        self.gradientE = np.zeros_like(self.band)
         # We need the gradient norm to calculate the action
         self.gradientENorm = np.zeros(self.n_images)
 
@@ -243,7 +245,7 @@ def compute_effective_field_and_energy(self, y):
         # Do not update the extreme images
         for i in range(1, len(y) - 1):
 
-            self.sim.set_m(y[i])
+            self.sim.spin[:] = y[i]
             # elif self.coordinates == 'Cartesian':
             #     self.sim.set_m(self.band[i])
 
@@ -263,9 +265,9 @@ def compute_tangents(self, y):
         nebm_clib.project_images(self.tangents, y,
                                  self.n_images, self.n_dofs_image
                                  )
-        # nebm_clib.normalise_images(self.tangents,
-        #                            self.n_images, self.n_dofs_image
-        #                            )
+        nebm_clib.normalise_images(self.tangents,
+                                   self.n_images, self.n_dofs_image
+                                   )
 
     def compute_spring_force(self, y):
         """
@@ -321,13 +323,18 @@ def compute_action(self):
         #                                       )
 
         # NOTE: Gradient here is projected in the S2^N tangent space
+        # nebm_clib.project_images(self.gradientE, self.band, self.n_images, self.n_dofs_image)
         self.gradientE.shape = (self.n_images, -1)
-        # NOTE: HEre we have to divide by the number of spins per image,not n_images:
+        # # NOTE: HEre we have to divide by the number of spins per image,not n_images:
         Gnorms2 = np.sum(self.gradientE**2, axis=1) / self.n_spins
+        # self.G.shape = (self.n_images, -1)
+        # Gnorms2 = np.sum(self.G**2, axis=1) / self.n_spins
 
         # Compute the root mean square per image
-        self.gradientENorm[:] = np.sqrt(Gnorms2)
+        # self.gradientENorm[:] = np.sqrt(Gnorms2)
+        self.gradientENorm[:] = np.linalg.norm(self.gradientE, axis=1) / self.n_spins
         self.gradientE.shape = (-1)
+        # self.G.shape = (-1)
 
         # DEBUG:
         # print('gradEnorm', self.gradientENorm)
@@ -374,8 +381,11 @@ def nebm_step(self, y, ensure_zero_extrema=False):
 
         self.compute_effective_field_and_energy(y)
         nebm_clib.project_images(self.gradientE, y, self.n_images, self.n_dofs_image)
-        self.compute_tangents(y)
-        nebm_clib.normalise_spins(self.tangents, self.n_images, self.n_dofs_image)
+        # print('GRADE', self.gradientE)
+        self.compute_tangents(y)  # In the function: tangents -> project -> normalise
+        # self.tangents.shape = (self.n_images, -1)
+        # print(f'tangent norms', np.linalg.norm(self.tangents, axis=1))
+        # self.tangents.shape = (-1)
         self.compute_spring_force(y)
 
         nebm_clib.compute_effective_force(self.G,
@@ -386,6 +396,27 @@ def nebm_step(self, y, ensure_zero_extrema=False):
                                           self.n_images,
                                           self.n_dofs_image
                                           )
+        # self.G.shape = (self.n_images, -1, 3)
+        # self.tangents.shape = (self.n_images, -1, 3)
+        # self.gradientE.shape = (self.n_images, -1, 3)
+        # for im_idx in range(self.G.shape[0]):
+        #     for s_idx in range(self.G.shape[1]):
+        #         gradE_dot_t = np.dot(self.gradientE[im_idx, s_idx, :], self.tangents[im_idx, s_idx, :])
+        #         self.G[im_idx, s_idx, :] = -self.gradientE[im_idx, s_idx, :] + gradE_dot_t * self.tangents[im_idx, s_idx, :]
+        # self.G.shape = (-1)
+        # self.tangents.shape = (-1)
+        # self.gradientE.shape = (-1)
+
+
+        # tgts = self.tangents.reshape(self.n_images, -1)
+        # forces = self.G.reshape(self.n_images, -1)
+        # for im_idx in range(forces.shape[0]):
+        #     gradE_dot_tgt = np.sum(forces[im_idx] * tgts[im_idx])
+        #     print(f'Im {im_idx}', forces[im_idx], tgts[im_idx], gradE_dot_tgt)
+        #     print('-----', np.linalg.norm(tgts[im_idx]))
+        #     print(self.spring_force.reshape(self.n_images, -1)[im_idx])
+        # tgts = self.tangents.reshape(self.n_images, -1)
+        # forces = self.G.reshape(self.n_images, -1)
 
         # The effective force at the extreme images should already be zero, but
         # we will manually remove any value