summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--langevin.cpp92
-rw-r--r--stokes.hpp129
2 files changed, 154 insertions, 67 deletions
diff --git a/langevin.cpp b/langevin.cpp
index 5ca2d72..5d9acc9 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -1,8 +1,10 @@
#include <getopt.h>
+#include <stereographic.hpp>
#include "complex_normal.hpp"
#include "p-spin.hpp"
#include "dynamics.hpp"
+#include "stokes.hpp"
#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
@@ -26,7 +28,7 @@ int main(int argc, char* argv[]) {
// simulation parameters
double ε = 1e-4;
- double εJ = 1e-2;
+ double εJ = 1e-5;
double δ = 1e-2; // threshold for determining saddle
double Δ = 1e-3;
double γ = 1e-2; // step size
@@ -92,28 +94,6 @@ int main(int argc, char* argv[]) {
return W;
};
- double aGoal = 1e3;
-
- std::function<double(const ComplexTensor&, const ComplexVector&)> energyInvA = [aGoal]
- (const ComplexTensor& J, const ComplexVector& z) {
- double a = z.squaredNorm();
- if (a > aGoal) {
- return -aGoal;
- } else {
- return -z.squaredNorm();
- }
- };
-
- while (zSaddle.squaredNorm() < aGoal) {
- std::tie(std::ignore, z) = metropolis(J, z, energyInvA, T, γ, 100, d, r.engine());
- try {
- std::cerr << "Starting descent from " << z.squaredNorm() << "." << std::endl;
- zSaddle = findSaddle(J, z, ε);
- } catch (std::exception& e) {
- }
- std::cerr << "Current saddle is of size " << zSaddle.squaredNorm() << "." << std::endl;
- }
-
ComplexVector zSaddlePrev = ComplexVector::Zero(N);
while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close...
@@ -131,7 +111,11 @@ int main(int argc, char* argv[]) {
std::cerr << "Current saddles are " << (zSaddle - zSaddlePrev).norm() << " apart." << std::endl;
}
- std::cerr << "Found sufficiently nearby saddles, perturbing J." << std::endl;
+ Rope<Complex> stokest(20, zSaddle, zSaddlePrev);
+
+ stokest.relax(J, 10000, 0.0001);
+
+ std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl;
complex_normal_distribution<> dJ(0, εJ * σ, 0);
@@ -141,20 +125,68 @@ int main(int argc, char* argv[]) {
setJ<Complex, p>(JJ, ind, Ji + dJ(r.engine()));
};
- for (unsigned i = 0; i < n; i++) {
- ComplexTensor Jp = J;
+ ComplexTensor Jp = J;
+ ComplexVector z1, z2;
+ Complex H1, H2;
+ double prevdist = 100;
- iterateOver<Complex, p>(Jp, perturbJ);
+ while (true) {
+ ComplexTensor Jpp = Jp;
+ iterateOver<Complex, p>(Jpp, perturbJ);
try {
- ComplexVector zSaddleNew = findSaddle(Jp, zSaddle, ε);
- ComplexVector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε);
+ z1 = findSaddle(Jpp, zSaddle, ε);
+ z2 = findSaddle(Jpp, zSaddlePrev, ε);
+
+ double dist = (z1 - z2).norm();
+
+ std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jpp, z1);
+ std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jpp, z2);
+
+ if (abs(imag(H1 - H2)) < prevdist && dist > 1e-2) {
+ Jp = Jpp;
+ prevdist = abs(imag(H1 - H2));
+
+ std::cout << abs(imag(H1 - H2)) << " " << dist << std::endl;
+ if (abs(imag(H1 - H2)) < 1e-8 && dist > 1e-2) {
+ std::cout << "Found distinct critical points with sufficiently similar Im H." << std::endl;
+ break;
+ }
+ }
+
+
- std::cout << zSaddleNew.transpose() << " " << zSaddlePrevNew.transpose() << std::endl;
} catch (std::exception& e) {
std::cerr << "Couldn't find a saddle with new couplings, skipping." << std::endl;
}
}
+ Rope<Complex> stokes(20, z1, z2);
+
+ std::cout << stokes.cost(Jp) << std::endl;
+
+ stokes.relax(Jp, 10000, 0.0001);
+
+ std::cout << stokes.cost(Jp) << std::endl;
+
+ Rope<Complex> stokes2 = stokes.interpolate();
+
+ stokes2.relax(Jp, 10000, 0.001);
+
+ std::cout << stokes2.cost(Jp) << std::endl;
+
+ stokes2 = stokes2.interpolate();
+
+ stokes2.relax(Jp, 10000, 0.001);
+
+ std::cout << stokes2.cost(Jp) << std::endl;
+
+ stokes2 = stokes2.interpolate();
+
+ stokes2.relax(Jp, 10000, 0.001);
+
+ std::cout << stokes2.cost(Jp) << std::endl;
+
+
return 0;
}
diff --git a/stokes.hpp b/stokes.hpp
index b60ae0e..43a69e3 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,47 +1,96 @@
-#include "stereographic.hpp"
+#include "p-spin.hpp"
+#include <iostream>
+#include "dynamics.hpp"
template <class Scalar>
-class Rope {
- private:
- std::vector<Vector<Scalar>> z;
+Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
+ return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate();
+}
- template <int p>
- Vector<Scalar> δz(const Tensor<Scalar, p>& J) const {
- Vector<Scalar> dz(z.size());
-
- for (unsigned i = 1; i < z.size() - 1; i++) {
- Vector<Scalar> z12 = (z[i] + z[i - 1]) / 2.0;
- Vector<Scalar> g;
- std::tie(std::ignore, g, std::ignore) = stereographicHamGradHessHess(J, z12);
+template <class Scalar>
+double segmentCost(const Vector<Scalar>& z, const Vector<Scalar>& dz, const Vector<Scalar>& dH) {
+ Vector<Scalar> zD = zDot(z, dH);
+ return 1.0 - pow(real(zD.dot(dz)), 2) / zD.squaredNorm() / dz.squaredNorm();
+}
- Vector<Scalar> Δ = z[i] - z[i - 1];
+template <class Scalar>
+Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, const Vector<Scalar>& z´´, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) {
+ double z² = z.squaredNorm();
+ double z´² = z´.squaredNorm();
+
+ Vector<Scalar> ż = zDot(z, dH);
+ double ż² = ż.squaredNorm();
+
+ double Reż·z´² = real(pow(ż.dot(z´), 2));
+
+ Matrix<Scalar> dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²;
+ Matrix<Scalar> dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z²
+ + (z.dot(dH) / z²) * (Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z²);
+
+ Vector<Scalar> dLdz = - 0.5 * (
+ ż.dot(z´) * (dżc * z´) + z´.dot(ż) * (dż * z´.conjugate())
+ - (dż * ż.conjugate() + dżc * ż) * Reż·z´² / ż²
+ ) / (ż² * z´²);
+
+ Vector<Scalar> ż´ = -(ddH * z´).conjugate() + ((ddH * z´).dot(z) / z²) * z.conjugate() + (
+ dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - (dH.dot(z) * (z´.dot(z) + z.dot(z´)) / z²) * z.conjugate()
+ ) / z²;
+
+ Vector<Scalar> ddtdLdz´ = -0.5 * (
+ ż´.dot(z´) * ż.conjugate() + ż.dot(z´´) * ż.conjugate() + ż.dot(z´) * ż´.conjugate()
+ - (
+ Reż·z´² * z´´.conjugate() + real(2.0 * ż.dot(z´) * (ż.dot(z´´) + ż´.dot(z´))) * z´.conjugate()
+ - Reż·z´² * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() / z´²
+ ) / z´²
+ - (
+ (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´²
+ ) * (
+ ż.dot(z´) * ż.conjugate() - Reż·z´² / z´² * z´.conjugate()
+ )
+ ) / (ż² * z´²);
+
+ return dLdz - ddtdLdz´;
+}
- g.normalize();
- Δ.normalize();
+template <class Scalar>
+class Rope {
+ public:
+ std::vector<Vector<Scalar>> z;
- if (abs(arg(g.dot(Δ))) < M_PI / 2) {
- dz[i] = g - Δ;
- } else {
- dz[i] = -g - Δ;
- }
+ Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2) : z(N + 2) {
+ for (unsigned i = 0; i < N + 2; i++) {
+ z[i] = z1 + (z2 - z1) * ((double)i / (N + 1.0));
+ z[i] = normalize(z[i]);
}
+ }
+
+ Vector<Scalar> dz(unsigned i) const {
+ return z[i + 1] - z[i - 1];
+ }
- return dz;
+ Vector<Scalar> ddz(unsigned i) const {
+ return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]);
}
template <int p>
void perturb(const Tensor<Scalar, p>& J, double δ) {
- z += δ * this->δz(J);
+ for (unsigned i = 1; i < z.size() - 1; i++) {
+ Vector<Scalar> dH;
+ Matrix<Scalar> ddH;
+ std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]);
+
+ Vector<Scalar> δz = variation(z[i], this->dz(i), this->ddz(i), dH, ddH);
+
+ z[i] -= δ * δz.conjugate();
+ z[i] = normalize(z[i]);
+ }
}
void spread() {
- std::vector<double> d(z.size() - 1);
double l = 0;
for (unsigned i= 0; i < z.size() - 1; i++) {
- double Δ = (z[i + 1] - z[i]).norm();
- d[i] = Δ;
- l += Δ;
+ l += (z[i + 1] - z[i]).norm();
}
double a = 0;
@@ -49,19 +98,16 @@ class Rope {
for (unsigned i = 1; i < z.size() - 1; i++) {
double b = i * l / (z.size() - 1);
- while (b >= a) {
- a += d[pos];
+
+ while (b > a) {
pos++;
+ a += (z[pos] - z[pos - 1]).norm();
}
- z[i] = z[pos] - (a - b) * (z[pos] - z[pos - 1]).normalized();
- }
- }
+ Vector<Scalar> δz = z[pos] - z[pos - 1];
- public:
- Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2) : z(N + 2) {
- for (unsigned i = 0; i < N + 2; i++) {
- z[i] = z1 + (z2 - z1) * ((double)i / (N + 1.0));
+ z[i] = z[pos] - (a - b) / δz.norm() * δz;
+ z[i] = normalize(z[i]);
}
}
@@ -73,7 +119,16 @@ class Rope {
template <int p>
double cost(const Tensor<Scalar, p>& J) const {
- return this->δz(J).norm();
+ double c = 0;
+
+ for (unsigned i = 1; i < z.size() - 1; i++) {
+ Vector<Scalar> dH;
+ std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z[i]);
+
+ c += segmentCost(z[i], this->dz(i), dH);
+ }
+
+ return c;
}
template <int p>
@@ -84,7 +139,7 @@ class Rope {
}
Rope<Scalar> interpolate() const {
- Rope<Scalar> r(2 * z.size() - 1);
+ Rope<Scalar> r(2 * z.size() - 1, z[0], z[z.size() - 1]);
for (unsigned i = 0; i < z.size(); i++) {
r.z[2 * i] = z[i];