From 11cc46f5877fa814c891f715aadee80057522095 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 18 Aug 2024 13:15:43 +0200 Subject: Generalized topology code to extensively many constraints. --- topology.cpp | 217 +++++++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 159 insertions(+), 58 deletions(-) diff --git a/topology.cpp b/topology.cpp index 9483c48..9e755c0 100644 --- a/topology.cpp +++ b/topology.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" @@ -19,6 +20,10 @@ class Tensor : public Eigen::Tensor { using Eigen::Tensor::Tensor; public: + Tensor(const Matrix& A) { + *this = Eigen::TensorMap>(A.data(), 1, A.rows(), A.cols()); + } + Matrix operator*(const Vector& x) const { std::array, 1> ip20 = {Eigen::IndexPair(2, 0)}; const Eigen::Tensor xT = Eigen::TensorMap>(x.data(), x.size()); @@ -27,52 +32,43 @@ public: } }; +Matrix operator*(const Eigen::Matrix& x, const Tensor& J) { + std::array, 1> ip00 = {Eigen::IndexPair(0, 0)}; + const Eigen::Tensor xT = Eigen::TensorMap>(x.data(), x.size()); + const Eigen::Tensor JxT = J.contract(xT, ip00); + return Eigen::Map(JxT.data(), J.dimension(1), J.dimension(2)); +} + Vector normalize(const Vector& x) { return x * sqrt(x.size() / x.squaredNorm()); } -class Spherical3Spin { -private: - Tensor J; +class ConstraintModel { +protected: + Vector minusE(const Vector& H) const { + return H - Vector::Constant(M, E); + } public: + Vector x0; + Real E; unsigned N; + unsigned M; - Spherical3Spin(unsigned N, Rng& r) : J(N, N, N), N(N) { - Eigen::StaticSGroup, Eigen::Symmetry<1,2>> sym123; - - 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(0, sqrt(12) / N); - } - } + ConstraintModel(unsigned N, unsigned M, Real E, Rng& r) : x0(N), E(E), N(N), M(M) { + for (Real& x0i : x0) { + x0i = r.variate(); } - } - std::tuple H_∂H_∂∂H(const Vector& x) const { - Matrix ∂∂H = J * x; - Vector ∂H = ∂∂H * x / 2; - Real H = ∂H.dot(x) / 6; - return {H, ∂H, ∂∂H}; + x0 = normalize(x0); } -}; -class ConstrainedHeight { - private: - Vector x0; - Real E; - - public: - Spherical3Spin S; - unsigned N; - - ConstrainedHeight(unsigned N, Real E, Rng& r) : x0(N), E(E), S(N, r), N(N) { - for (Real& x0ᵢ : x0) { - x0ᵢ = r.variate(); - } - - x0 = normalize(x0); + virtual std::tuple H_∂H_∂∂H(const Vector& x) const { + Tensor t(M, N, N); + Tensor ∂∂H = t.setConstant(0); + Matrix ∂H = Matrix::Zero(M, N); + Vector H = minusE(Vector::Zero(M)); + return {H, ∂H, ∂∂H}; } Real overlap(const Vector& v) const { @@ -82,31 +78,31 @@ class ConstrainedHeight { std::tuple ∂L_∂∂L(const Vector& v) const { Vector x = v.head(N); Real ω₀ = v(N); - Real ω₁ = v(N + 1); + Vector ω₁ = v.tail(M); - auto [H, ∂H∂x, ∂²H∂²x] = S.H_∂H_∂∂H(x); + auto [H, ∂H∂x, ∂∂H∂∂x] = H_∂H_∂∂H(x); - Vector ∂L∂x = x0 + ω₀ * x + ω₁ * ∂H∂x; + Vector ∂L∂x = x0 + ω₀ * x + ∂H∂x.transpose() * ω₁; Real ∂L∂ω₀ = 0.5 * (x.squaredNorm() - N); - Real ∂L∂ω₁ = H - N * E; + Vector ∂L∂ω₁ = H; - Vector ∂L(N + 2); + Vector ∂L(N + 1 + M); ∂L << ∂L∂x, ∂L∂ω₀, ∂L∂ω₁; - Matrix ∂L²∂x² = ω₀ * Matrix::Identity(N, N) + ω₁ * ∂²H∂²x; + Matrix ∂L²∂x² = ω₁ * ∂∂H∂∂x + ω₀ * Matrix::Identity(N, N); Vector ∂²L∂x∂ω₀ = x; - Vector ∂²L∂x∂ω₁ = ∂H∂x; + Matrix ∂²L∂x∂ω₁ = ∂H∂x; - Matrix ∂∂L(N + 2, N + 2); + Matrix ∂∂L(N + 1 + M, N + 1 + M); ∂∂L << - ∂L²∂x², ∂²L∂x∂ω₀, ∂²L∂x∂ω₁, - ∂²L∂x∂ω₀.transpose(), 0, 0, - ∂²L∂x∂ω₁.transpose(), 0, 0; + ∂L²∂x², ∂²L∂x∂ω₀, ∂²L∂x∂ω₁.transpose(), + ∂²L∂x∂ω₀.transpose(), Vector::Zero(M + 1).transpose(), + ∂²L∂x∂ω₁, Matrix::Zero(M, M + 1); return {∂L, ∂∂L}; } - Vector newtonMethod(const Vector& v0, Real γ = 1) { + Vector newtonMethod(const Vector& v0, Real γ = 1) const { Vector v = v0; Vector ∂L; @@ -121,20 +117,76 @@ class ConstrainedHeight { } }; +class Spherical3Spin : public ConstraintModel { +private: + Tensor J; + +public: + Spherical3Spin(unsigned N, Real E, Rng& r) : ConstraintModel(N, 1, N * E, r), J(N, N, N) { + Eigen::StaticSGroup, Eigen::Symmetry<1,2>> sym123; + + 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(0, sqrt(12) / N); + } + } + } + } + + std::tuple H_∂H_∂∂H(const Vector& x) const override { + Tensor ∂∂H = J * x; + Matrix ∂H = ∂∂H * x / 2; + Vector H = minusE(∂H * x / 6); + return {H, ∂H, ∂∂H}; + } +}; + +class ExtensiveQuadratic : public ConstraintModel { +private: + Tensor J; + +public: + ExtensiveQuadratic(unsigned N, unsigned M, Real E, Rng& r) : ConstraintModel(N, M, E, r), J(M, N, N) { + Eigen::StaticSGroup> sym23; + + for (unsigned i = 0; i < M; i++) { + for (unsigned j = 0; j < N; j++) { + for (unsigned k = j; k < N; k++) { + sym23(J, i, j, k) = r.variate(0, 1.0 / N); + } + } + } + } + + std::tuple H_∂H_∂∂H(const Vector& x) const override { + Tensor ∂∂H = J; + Matrix ∂H = ∂∂H * x; + Vector H = minusE(∂H * x / 2); + return {H, ∂H, ∂∂H}; + } +}; int main(int argc, char* argv[]) { unsigned N = 10; + double α = 0.25; double E = 0; double γ = 1; unsigned samples = 10; + bool count = false; + bool quadratic = false; + int opt; - while ((opt = getopt(argc, argv, "N:E:g:n:")) != -1) { + while ((opt = getopt(argc, argv, "N:a:E:g:n:cq")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); break; + case 'a': + α = atof(optarg); + break; case 'E': E = atof(optarg); break; @@ -144,26 +196,75 @@ int main(int argc, char* argv[]) { case 'n': samples = atoi(optarg); break; + case 'c': + count = true; + break; + case 'q': + quadratic = true; + break; default: exit(1); } } - Rng r; + unsigned M = std::round(α * N); - Vector x = Vector::Zero(N); - x(0) = sqrt(N); + Rng r; std::cout << std::setprecision(15); - for (unsigned sample = 0; sample < samples; sample++) { - Vector v0(N + 2); - v0 << x, - r.variate(), - r.variate(); - ConstrainedHeight M(N, E, r); - Vector v = M.newtonMethod(v0, γ); - std::cout << M.overlap(v) << std::endl; + if (quadratic) { + Vector x = Vector::Zero(N); + x(0) = sqrt(N); + + for (unsigned sample = 0; sample < samples; sample++) { + Vector v0(N + M + 1); + v0.head(N) = x; + for (Real& v0ᵢ : v0.tail(M + 1)) { + v0ᵢ = r.variate(); + } + ExtensiveQuadratic S(N, M, E, r); + Vector v = S.newtonMethod(v0, γ); + std::cout << S.overlap(v) << std::endl; + } + } else if (!count) { + Vector x = Vector::Zero(N); + x(0) = sqrt(N); + + for (unsigned sample = 0; sample < samples; sample++) { + Vector v0(N + 2); + v0 << x, + r.variate(), + r.variate(); + Spherical3Spin S(N, E, r); + Vector v = S.newtonMethod(v0, γ); + std::cout << S.overlap(v) << std::endl; + } + } else { + Spherical3Spin S(N, E, r); + std::list points; + + for (unsigned sample = 0; sample < samples; sample++) { + Vector v0(N + 2); + for (Real& v0i : v0) { + v0i = r.variate(); + } + v0.head(N) = normalize(v0.head(N)); + + Vector v = S.newtonMethod(v0, γ); + if (v(N) == v(N)) { + bool inList = false; + for (const Vector& u : points) { + if ((u - v).head(N).squaredNorm() < 1e-8) { + inList = true; + } + } + if (!inList) { + points.push_back(v); + } + } + } + std::cout << points.size() << std::endl; } return 0; -- cgit v1.2.3-70-g09d2