1 files changed, 24 insertions, 40 deletions

diff --git a/topology.cpp b/topology.cpp
index 6048e1e..9483c48 100644
--- a/topology.cpp
+++ b/topology.cpp
@@ -27,25 +27,10 @@ public:
   }
 };
 
-Matrix operator*(const Eigen::Matrix<Real, 1, Eigen::Dynamic>& x, const Tensor& J) {
-  std::array<Eigen::IndexPair<int>, 1> ip00 = {Eigen::IndexPair<int>(0, 0)};
-  const Eigen::Tensor<Real, 1> xT = Eigen::TensorMap<const Eigen::Tensor<Real, 1>>(x.data(), x.size());
-  const Eigen::Tensor<Real, 2> JxT = J.contract(xT, ip00);
-  return Eigen::Map<const Matrix>(JxT.data(), J.dimension(1), J.dimension(2));
-}
-
 Vector normalize(const Vector& x) {
   return x * sqrt(x.size() / x.squaredNorm());
 }
 
-Vector ∂HFromV∂V(const Vector& V, const Matrix& ∂V) {
-  return V.transpose() * ∂V;
-}
-
-Vector VFromABJx(const Vector& b, const Matrix& A, const Matrix& Jx, const Vector& x) {
-  return b + (A + 0.5 * Jx) * x;
-}
-
 class Spherical3Spin {
 private:
   Tensor J;
@@ -59,16 +44,16 @@ public:
     for (unsigned i = 0; i < N; i++) {
       for (unsigned j = i; j < N; j++) {
         for (unsigned k = j; k < N; k++) {
-          sym123(J, i, j, k) = r.variate<Real, std::normal_distribution>(0, sqrt(2) / N);
+          sym123(J, i, j, k) = r.variate<Real, std::normal_distribution>(0, sqrt(12) / N);
         }
       }
     }
   }
 
   std::tuple<Real, Vector, Matrix> H_∂H_∂∂H(const Vector& x) const {
-    Matrix ∂∂H = J * x / (6 * N);
-    Vector ∂H = ∂∂H * x;
-    Real   H = ∂H.dot(x);
+    Matrix ∂∂H = J * x;
+    Vector ∂H = ∂∂H * x / 2;
+    Real   H = ∂H.dot(x) / 6;
     return {H, ∂H, ∂∂H};
   }
 };
@@ -77,9 +62,9 @@ class ConstrainedHeight {
   private:
     Vector x0;
     Real E;
-    Spherical3Spin S;
 
   public:
+    Spherical3Spin S;
     unsigned N;
 
   ConstrainedHeight(unsigned N, Real E, Rng& r) : x0(N), E(E), S(N, r), N(N) {
@@ -99,14 +84,11 @@ class ConstrainedHeight {
     Real ω₀ = v(N);
     Real ω₁ = v(N + 1);
 
-    Real H;
-    Vector ∂H∂x;
-    Matrix ∂²H∂²x;
-    std::tie(H, ∂H∂x, ∂²H∂²x) = S.H_∂H_∂∂H(x);
+    auto [H, ∂H∂x, ∂²H∂²x] = S.H_∂H_∂∂H(x);
 
     Vector ∂L∂x = x0 + ω₀ * x + ω₁ * ∂H∂x;
     Real ∂L∂ω₀ = 0.5 * (x.squaredNorm() - N);
-    Real ∂L∂ω₁ = H - E;
+    Real ∂L∂ω₁ = H - N * E;
 
     Vector ∂L(N + 2);
     ∂L << ∂L∂x, ∂L∂ω₀, ∂L∂ω₁;
@@ -116,7 +98,10 @@ class ConstrainedHeight {
     Vector ∂²L∂x∂ω₁ = ∂H∂x;
 
     Matrix ∂∂L(N + 2, N + 2);
-    ∂∂L << ∂L²∂x², ∂²L∂x∂ω₀, ∂²L∂x∂ω₁, ∂²L∂x∂ω₀.transpose(), 0, 0, ∂²L∂x∂ω₁.transpose(), 0, 0;
+    ∂∂L <<
+      ∂L²∂x², ∂²L∂x∂ω₀, ∂²L∂x∂ω₁,
+      ∂²L∂x∂ω₀.transpose(), 0, 0,
+      ∂²L∂x∂ω₁.transpose(), 0, 0;
 
     return {∂L, ∂∂L};
   }
@@ -129,6 +114,7 @@ class ConstrainedHeight {
 
     while (std::tie(∂L, ∂∂L) = ∂L_∂∂L(v), ∂L.squaredNorm() > 1e-10) {
       v -= γ * ∂∂L.partialPivLu().solve(∂L);
+      v.head(N) = normalize(v.head(N)); // might as well stay on the sphere
     }
 
     return v;
@@ -139,11 +125,12 @@ class ConstrainedHeight {
 int main(int argc, char* argv[]) {
   unsigned N = 10;
   double E = 0;
+  double γ = 1;
   unsigned samples = 10;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:E:n:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:E:g:n:")) != -1) {
     switch (opt) {
     case 'N':
       N = (unsigned)atof(optarg);
@@ -151,6 +138,9 @@ int main(int argc, char* argv[]) {
     case 'E':
       E = atof(optarg);
       break;
+    case 'g':
+      γ = atof(optarg);
+      break;
     case 'n':
       samples = atoi(optarg);
       break;
@@ -161,24 +151,18 @@ int main(int argc, char* argv[]) {
 
   Rng r;
 
-  ConstrainedHeight M(N, E, r);
-
-  Vector x(N);
-  for (Real& xᵢ : x) {
-    xᵢ = r.variate<Real, std::normal_distribution>();
-  }
-  x = normalize(x);
-
-  Vector v0(N + 2);
-  v0.head(N) = x;
-  v0(N) = r.variate<Real, std::normal_distribution>();
-  v0(N+1) = r.variate<Real, std::normal_distribution>();
+  Vector x = Vector::Zero(N);
+  x(0) = sqrt(N);
 
   std::cout << std::setprecision(15);
 
   for (unsigned sample = 0; sample < samples; sample++) {
+    Vector v0(N + 2);
+    v0 << x,
+      r.variate<Real, std::normal_distribution>(),
+      r.variate<Real, std::normal_distribution>();
     ConstrainedHeight M(N, E, r);
-    Vector v = M.newtonMethod(v0, 0.05);
+    Vector v = M.newtonMethod(v0, γ);
     std::cout << M.overlap(v) << std::endl;
   }