diff options
| -rw-r--r-- | topology.cpp | 64 |
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; } |