summaryrefslogtreecommitdiff
path: root/topology.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'topology.cpp')
-rw-r--r--topology.cpp217
1 files 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 <getopt.h>
#include <iomanip>
#include <random>
+#include <list>
#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
@@ -19,6 +20,10 @@ class Tensor : public Eigen::Tensor<Real, 3> {
using Eigen::Tensor<Real, 3>::Tensor;
public:
+ Tensor(const Matrix& A) {
+ *this = Eigen::TensorMap<const Eigen::Tensor<Real, 3>>(A.data(), 1, A.rows(), A.cols());
+ }
+
Matrix operator*(const Vector& x) const {
std::array<Eigen::IndexPair<int>, 1> ip20 = {Eigen::IndexPair<int>(2, 0)};
const Eigen::Tensor<Real, 1> xT = Eigen::TensorMap<const Eigen::Tensor<Real, 1>>(x.data(), x.size());
@@ -27,52 +32,43 @@ 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());
}
-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<0,1>, 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<Real, std::normal_distribution>(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<Real, std::normal_distribution>();
}
- }
- std::tuple<Real, Vector, Matrix> 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<Real, std::normal_distribution>();
- }
-
- x0 = normalize(x0);
+ virtual std::tuple<Vector, Matrix, Tensor> 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<Vector, Matrix> ∂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<0,1>, 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<Real, std::normal_distribution>(0, sqrt(12) / N);
+ }
+ }
+ }
+ }
+
+ std::tuple<Vector, Matrix, Tensor> 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<Eigen::Symmetry<1,2>> 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<Real, std::normal_distribution>(0, 1.0 / N);
+ }
+ }
+ }
+ }
+
+ std::tuple<Vector, Matrix, Tensor> 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<Real, std::normal_distribution>(),
- r.variate<Real, std::normal_distribution>();
- 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<Real, std::normal_distribution>();
+ }
+ 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<Real, std::normal_distribution>(),
+ r.variate<Real, std::normal_distribution>();
+ 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<Vector> points;
+
+ for (unsigned sample = 0; sample < samples; sample++) {
+ Vector v0(N + 2);
+ for (Real& v0i : v0) {
+ v0i = r.variate<Real, std::normal_distribution>();
+ }
+ 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;