nedler-mead and skip grad in sb model

cmtj/models/sb.py

@@ -102,15 +102,16 @@ def ext_field(self, Hinplane: VectorObj):
         return -self.Ms * mu0 * (hx * self.stheta * self.cphi + hy *
                                  self.sphi * self.ctheta + hz * self.ctheta)
 
-    def grad_ext_field(self, Hinplane: VectorObj):
+    def grad_ext_field(self, Hinplane: VectorObj, full_grad: bool):
         hx, hy, hz = Hinplane.get_cartesian()
         pref = -self.Ms * mu0
         dEdtheta = pref * (hx * self.cphi * self.ctheta +
                            hy * self.sphi * self.ctheta - hz * self.stheta)
 
         dEdphi = pref * (-hx * self.sphi * self.stheta +
                          hy * self.stheta * self.cphi)
-
+        if not full_grad:
+            return [dEdtheta, dEdphi, 0, 0, 0]
         d2Edtheta2 = pref * (-hx * self.stheta * self.cphi -
                              hy * self.sphi * self.stheta - hz * self.ctheta)
         d2Edphi2 = pref * (-hx * self.stheta * self.cphi -
@@ -131,19 +132,23 @@ def iec_interaction(self, J, layer: "LayerSB"):
         return -constJ * (m1x * m2x + m1y * m2y + m1z * m2z)
 
     def grad_iec_interaction(self, J_top, J_bottom, top_layer: "LayerSB",
-                             bottom_layer: "LayerSB"):
+                             bottom_layer: "LayerSB", full_grad: bool):
         """
         IEC gradient is not symmetric and is computed per layer
         Compute joint interaction on the layers from the top and bottom
         """
-        grad1 = self.__any_grad_iec_interaction(J_top, top_layer)
-        grad2 = self.__any_grad_iec_interaction(J_bottom, bottom_layer)
+        grad1 = self.__any_grad_iec_interaction(J_top,
+                                                top_layer,
+                                                full_grad=full_grad)
+        grad2 = self.__any_grad_iec_interaction(J_bottom,
+                                                bottom_layer,
+                                                full_grad=full_grad)
         return [
             grad1[0] + grad2[0], grad1[1] + grad2[1], grad1[2] + grad2[2],
             grad1[3] + grad2[3], grad1[4] + grad2[4]
         ]
 
-    def __any_grad_iec_interaction(self, J, layer: "LayerSB"):
+    def __any_grad_iec_interaction(self, J, layer: "LayerSB", full_grad: bool):
         # from any of the layers
         if (layer is None) or J == 0:
             return [0, 0, 0, 0, 0]
@@ -158,7 +163,8 @@ def __any_grad_iec_interaction(self, J, layer: "LayerSB"):
         dEdphi = -constJ * (-self.sphi * self.stheta * math.sin(
             layer.m.theta) * math.cos(layer.m.phi) + math.sin(layer.m.phi) *
                             self.stheta * math.sin(layer.m.theta) * self.cphi)
-
+        if not full_grad:
+            return [dEdtheta, dEdphi, 0, 0, 0]
         d2Edtheta2 = -constJ * (
             -self.sphi * layer.sphi * self.stheta * layer.stheta -
             self.stheta * layer.stheta * self.cphi * layer.cphi -
@@ -185,19 +191,23 @@ def dmi_interaction(self, D, layer: "LayerSB"):
         return D * (m1x * m2y - m1y * m2x)
 
     def grad_dmi_interaction(self, D_top, D_bottom, top_layer: "LayerSB",
-                             bottom_layer: "LayerSB"):
+                             bottom_layer: "LayerSB", full_grad: bool):
         """
         DMI energy and its gradient are both not symmetric and is computed per layer
         Compute joint interaction on the layers from the top and bottom
         """
-        grad1 = self.__any_grad_dmi_interaction(D_top, top_layer)
-        grad2 = self.__any_grad_dmi_interaction(D_bottom, bottom_layer)
+        grad1 = self.__any_grad_dmi_interaction(D_top,
+                                                top_layer,
+                                                full_grad=full_grad)
+        grad2 = self.__any_grad_dmi_interaction(D_bottom,
+                                                bottom_layer,
+                                                full_grad=full_grad)
         return [
             grad1[0] + grad2[0], grad1[1] + grad2[1], grad1[2] + grad2[2],
             grad1[3] + grad2[3], grad1[4] + grad2[4]
         ]
 
-    def __any_grad_dmi_interaction(self, D, layer: "LayerSB"):
+    def __any_grad_dmi_interaction(self, D, layer: "LayerSB", full_grad: bool):
         if (layer is None) or D == 0:
             return [0, 0, 0, 0, 0]
         constD = D
@@ -207,7 +217,8 @@ def __any_grad_dmi_interaction(self, D, layer: "LayerSB"):
 
         dEdphi = -constD * self.stheta * math.sin(
             layer.m.theta) * math.cos(self.m.phi - layer.m.phi)
-
+        if not full_grad:
+            return [dEdtheta, dEdphi, 0, 0, 0]
         d2Edtheta2 = constD * self.stheta * math.sin(self.m.phi - layer.m.phi)
 
         d2Edphi2 = d2Edtheta2
@@ -220,7 +231,7 @@ def __any_grad_dmi_interaction(self, D, layer: "LayerSB"):
     def surface_anisotropy(self):
         return (-self.Ks + 0.5 * self.Ms * mu0) * self.ctheta**2
 
-    def grad_surface_anisotropy(self):
+    def grad_surface_anisotropy(self, full_grad: bool):
         msq = math.pow(self.Ms, 2) * mu0
         dEdphi = 0
         d2Edphi2 = 0
@@ -238,7 +249,7 @@ def volume_anisotropy(self):
         mx, my, _ = self.m.get_cartesian()
         return -self.Kv.mag * (mx * ax + my * ay)**2
 
-    def grad_volume_anisotropy(self):
+    def grad_volume_anisotropy(self, full_grad: bool):
         """
         E_k/dtheta = ~ (a_x * m_x) * m_x *dm_x/dtheta
         """
@@ -247,6 +258,8 @@ def grad_volume_anisotropy(self):
         dEdtheta = -2 * Kv * self.stheta * self.ctheta * (math.cos(angdiff)**2)
         dEdphi = -2 * Kv * (self.stheta**
                             2) * math.sin(angdiff) * math.cos(angdiff)
+        if not full_grad:
+            return [dEdtheta, dEdphi, 0, 0, 0]
         d2Edtheta2 = -2 * Kv * math.cos(2 * self.m.theta) * (math.cos(angdiff)
                                                              **2)
         d2Edphi2 = 2 * Kv * (self.stheta**2) * math.cos(2 * angdiff)
@@ -255,14 +268,22 @@ def grad_volume_anisotropy(self):
             math.cos(2 * self.Kv.phi - 2 * self.m.phi + 2 * self.m.theta))
         return [dEdtheta, dEdphi, d2Edtheta2, d2Edphi2, d2Edphidtheta]
 
-    def compute_grad_energy(self, Hinplane: VectorObj, Jtop: float,
-                            Jbottom: float, Dtop: float, Dbottom: float,
-                            top_layer: "LayerSB", bottom_layer: "LayerSB"):
-        g1 = self.grad_ext_field(Hinplane)
-        g2 = self.grad_surface_anisotropy()
-        g3 = self.grad_volume_anisotropy()
-        g4 = self.grad_iec_interaction(Jtop, Jbottom, top_layer, bottom_layer)
-        g5 = self.grad_dmi_interaction(Dtop, Dbottom, top_layer, bottom_layer)
+    def compute_grad_energy(self,
+                            Hinplane: VectorObj,
+                            Jtop: float,
+                            Jbottom: float,
+                            Dtop: float,
+                            Dbottom: float,
+                            top_layer: "LayerSB",
+                            bottom_layer: "LayerSB",
+                            full_grad: bool = False):
+        g1 = self.grad_ext_field(Hinplane, full_grad)
+        g2 = self.grad_surface_anisotropy(full_grad)
+        g3 = self.grad_volume_anisotropy(full_grad)
+        g4 = self.grad_iec_interaction(Jtop, Jbottom, top_layer, bottom_layer,
+                                       full_grad)
+        g5 = self.grad_dmi_interaction(Dtop, Dbottom, top_layer, bottom_layer,
+                                       full_grad)
         return [g1[i] + g2[i] + g3[i] + g4[i] + g5[i] for i in range(5)]
 
     def compute_energy(self, Hinplane: VectorObj, Jtop: float, Jbottom: float,
@@ -296,9 +317,14 @@ def compute_frequency_at_equilibrum(self, Hinplane: VectorObj, Jtop: float,
         :param top_layer: LayerSB definition of the layer above the current one.
         :param bottom layer: LayerSB definition of the layer below the current one."""
         (_, _, d2Edtheta2, d2Edphi2,
-         d2Edphidtheta) = self.compute_grad_energy(Hinplane, Jtop, Jbottom,
-                                                   Dtop, Dbottom, top_layer,
-                                                   bottom_layer)
+         d2Edphidtheta) = self.compute_grad_energy(Hinplane,
+                                                   Jtop,
+                                                   Jbottom,
+                                                   Dtop,
+                                                   Dbottom,
+                                                   top_layer,
+                                                   bottom_layer,
+                                                   full_grad=True)
 
         if self.stheta != 0.:
             fmr = (d2Edtheta2 * d2Edphi2 - d2Edphidtheta**2) / math.pow(
@@ -385,6 +411,28 @@ def compute_energy_step(self):
             self.history[f"theta_{i}"].append(layer.m.theta)
             self.history[f"phi_{i}"].append(layer.m.phi)
 
+    def transform_flat_position_to_layer_position(self, position_vector):
+        position_ = [
+            position_vector[n:n + 2] for n in range(0, len(position_vector), 2)
+        ]
+        return position_
+
+    def compute_energy(self, position_vector):
+        position_ = self.transform_flat_position_to_layer_position(
+            position_vector=position_vector)
+        self.update_layers(position_)
+        total_energy = 0
+        for i, layer in enumerate(self.layers):
+            Jtop, Jbottom, top_layer_handle, bottom_layer_handle = self.__get_interaction_constant(
+                self.layers, i, self.J)
+            Dtop, Dbottom, _, _ = self.__get_interaction_constant(
+                self.layers, i, self.D)
+            evec = layer.compute_energy(self.Hext, Jtop, Jbottom, Dtop,
+                                        Dbottom, top_layer_handle,
+                                        bottom_layer_handle)
+            total_energy += sum(evec)
+        return total_energy
+
     def update_layers(self, position_vector):
         # update layers
         for l, layer in enumerate(self.layers):
@@ -470,3 +518,114 @@ def adam_gradient_descent(self,
             self.update_layers(new_position)
         if not self.silent:
             self.print_summary(step)
+
+    def nelder_mead_method(self,
+                           max_steps: int,
+                           tol: float = 1e-8,
+                           alpha: float = 1.,
+                           gamma: float = 2.,
+                           rho: float = -0.5,
+                           sigma: float = 0.5):
+        """
+        A naive implementation of the Nedler-Mead method.
+        See: https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
+        Defaults taken from the wikipedia page.
+        :param max_steps: maximum number of steps.
+        :param tol: tolerance of the solution.
+        """
+
+        def simplex_centroid(points):
+            return np.mean(points, axis=0)
+
+        def reflection(point, centroid):
+            return centroid + alpha * (centroid - point)
+
+        def expansion(point, centroid):
+            return centroid + gamma * (point - centroid)
+
+        def contraction(point, centroid):
+            return centroid + rho * (point - centroid)
+
+        def shrink(points):
+            return np.asarray(
+                [points[0], *(points[0] + sigma * (points[1:] - points[0]))])
+
+        def point_energies(points):
+            return np.asarray([self.compute_energy(p) for p in points])
+
+        def sort_points(points, energies):
+            return points[np.argsort(energies)]
+
+        steps = 0
+        flat_coords = self.current_position.flatten()
+        points = [flat_coords]
+        for i in range(2):
+            points.append(flat_coords +
+                          np.random.normal(size=flat_coords.shape))
+
+        points = np.asarray(points)
+
+        while True:
+            # sort points by energy
+            energies = point_energies(points)
+            # sort points by energy, lowest first
+            points = sort_points(points, energies=energies)
+            # print("Best point", points[0])
+            position_ = self.transform_flat_position_to_layer_position(
+                position_vector=points[0])
+            self.update_layers(position_)
+            if self.position_difference(points[0]) < tol:
+                break
+            # do not allow more than max_steps
+            if steps > max_steps:
+                break
+            steps += 1
+            # compute excluding the worst point
+            centroid = simplex_centroid(points[:-1])
+            # reflection
+            reflected_point = reflection(points[-1], centroid)
+            reflected_point_energy = self.compute_energy(reflected_point)
+            # if reflected is better than the second worst point, replace worst
+            if (reflected_point_energy <
+                    energies[-2]) and (reflected_point_energy >= energies[0]):
+                points[-1] = reflected_point
+                continue
+
+            # if reflected point is better than the best point, expand
+            if reflected_point_energy < energies[0]:
+                expansion_point = expansion(reflected_point, centroid)
+                expansion_point_energy = self.compute_energy(expansion_point)
+                if expansion_point_energy < reflected_point_energy:
+                    points[-1] = expansion_point
+                else:
+                    points[-1] = reflected_point
+                continue
+
+            # if reflected point is better than the worst point, contract
+            if reflected_point_energy < energies[-1]:
+                contracted_point = contraction(reflected_point, centroid)
+                contracted_point_energy = self.compute_energy(contracted_point)
+                if contracted_point_energy < reflected_point_energy:
+                    points[-1] = contracted_point
+                else:
+                    points = shrink(points)
+            # if reflected point is worse than the worst point, contract
+            else:
+                contracted_point = contraction(points[-1], centroid)
+                contracted_point_energy = self.compute_energy(contracted_point)
+                if contracted_point_energy < energies[-1]:
+                    points[-1] = contracted_point
+                else:
+                    points = shrink(points)
+
+        # print("Best point", points[0])
+        position_ = self.transform_flat_position_to_layer_position(
+            position_vector=points[0])
+        self.update_layers(position_)
+        for i in range(len(self.layers)):
+            self.history[f"energy_{i}"].append(energies[0])
+            self.history[f"theta_{i}"].append(points[0][2 * i])
+            self.history[f"phi_{i}"].append(points[0][2 * i + 1])
+
+        if not self.silent:
+            self.print_summary(steps)

docs/api/models/dw-reference.md

@@ -0,0 +1 @@
+::: cmtj.models.domain_dynamics

docs/api/models/sb-reference.md

@@ -0,0 +1 @@
+::: cmtj.models.sb

examples/overview/examples.ipynb

@@ -576,7 +576,7 @@
    "metadata": {},
    "source": [
     "# CIMS\n",
-    "In this section we focus on generating switching currents. We are going to excite the system with a range of different pulse strengths. In this case, unlike in the PIMM, we excite the system directly with the current density pulse (which takes on the form of trapezoidal impulse) in order to evoke the torque response. After each simulation the system is allowed to relax to a steady state. We then pick up the steady state to see in which state (up/down) the magnetisation has settled. Based on the change of the hysteresis change, we can infer the critical current for a given value of the external field"
+    "In this section we focus on generating switching currents. We are going to excite the system with a range of different pulse strengths. In this case, unlike in the PIMM, we excite the system directly with the current density pulse (which takes on the form of trapezoidal impulse) in order to evoke the torque response. After each simulation the system is allowed to relax to a steady state. We then pick up the steady state to see in which state (up/down) the magnetisation has settled. Based on the change of the hysteresis change, we can infer the critical current for a given value of the external field."
    ]
   },
   {