From 02aeb1dc3a01818b6d737514251369dab7914075 Mon Sep 17 00:00:00 2001 From: LemurPwned Date: Sat, 14 Jan 2023 18:53:56 +0100 Subject: [PATCH] nedler-mead and skip grad in sb model --- cmtj/models/sb.py | 209 +++++++++++++++++++++++++++---- docs/api/models/dw-reference.md | 1 + docs/api/models/sb-reference.md | 1 + examples/overview/examples.ipynb | 2 +- 4 files changed, 187 insertions(+), 26 deletions(-) create mode 100644 docs/api/models/dw-reference.md create mode 100644 docs/api/models/sb-reference.md diff --git a/cmtj/models/sb.py b/cmtj/models/sb.py index f48c5c1..e356445 100644 --- a/cmtj/models/sb.py +++ b/cmtj/models/sb.py @@ -102,7 +102,7 @@ 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 + @@ -110,7 +110,8 @@ def grad_ext_field(self, Hinplane: VectorObj): 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) diff --git a/docs/api/models/dw-reference.md b/docs/api/models/dw-reference.md new file mode 100644 index 0000000..9e35199 --- /dev/null +++ b/docs/api/models/dw-reference.md @@ -0,0 +1 @@ +::: cmtj.models.domain_dynamics diff --git a/docs/api/models/sb-reference.md b/docs/api/models/sb-reference.md new file mode 100644 index 0000000..71fd541 --- /dev/null +++ b/docs/api/models/sb-reference.md @@ -0,0 +1 @@ +::: cmtj.models.sb diff --git a/examples/overview/examples.ipynb b/examples/overview/examples.ipynb index 680074b..1cf94e1 100644 --- a/examples/overview/examples.ipynb +++ b/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." ] }, {