summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2021-01-06 11:16:20 +0100
committerJaron Kent-Dobias <jaron@kent-dobias.com>2021-01-06 11:16:20 +0100
commit1094eaa434d791d3a511f853df9d6d66d169b2ab (patch)
tree18bd0248f207a45bb3fb977a86af350bb02c1a55
parent252ad65ae0af4249bc089cd18aad4ab739b80d70 (diff)
downloadcode-1094eaa434d791d3a511f853df9d6d66d169b2ab.tar.gz
code-1094eaa434d791d3a511f853df9d6d66d169b2ab.tar.bz2
code-1094eaa434d791d3a511f853df9d6d66d169b2ab.zip
Lots of work, and fixed bug stemming from choice of matrix algorithm.
-rw-r--r--langevin.cpp5
-rw-r--r--stereographic.hpp20
2 files changed, 14 insertions, 11 deletions
diff --git a/langevin.cpp b/langevin.cpp
index 3bc12ed..05b436e 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -7,6 +7,8 @@
#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
+#include <eigen3/Eigen/LU>
+
#include "complex_normal.hpp"
#include "p-spin.hpp"
#include "stereographic.hpp"
@@ -57,7 +59,8 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, doub
while (W > ε) {
double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
- Vector ζNew = ζ - ddH.ldlt().solve(dH);
+ Vector dζ = ddH.partialPivLu().solve(dH);
+ Vector ζNew = ζ - dζ;
double WNew;
Vector dWNew;
diff --git a/stereographic.hpp b/stereographic.hpp
index 59432b7..4c3c831 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -58,25 +58,25 @@ Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix
unsigned N = ζ.size() + 1;
Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N);
- Scalar b = ζ.transpose() * ζ;
+ Scalar b = 1.0 + (Scalar)(ζ.transpose() * ζ);
for (unsigned i = 0; i < N - 1; i++) {
for (unsigned j = 0; j < N - 1; j++) {
for (unsigned k = 0; k < N - 1; k++) {
- dJ(i, j, k) = 16 * sqrt(N) * ζ(i) * ζ(j) * ζ(k) / pow(1.0 + b, 3);
+ dJ(i, j, k) = 16 * sqrt(N) * ζ(i) * ζ(j) * ζ(k) / pow(b, 3);
if (i == j) {
- dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(k) / pow(1.0 + b, 3);
+ dJ(i, j, k) -= 4 * sqrt(N) * ζ(k) / pow(b, 2);
}
if (i == k) {
- dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(j) / pow(1.0 + b, 3);
+ dJ(i, j, k) -= 4 * sqrt(N) * ζ(j) / pow(b, 2);
}
if (j == k) {
- dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
+ dJ(i, j, k) -= 4 * sqrt(N) * ζ(i) / pow(b, 2);
}
}
- dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 3);
+ dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(b, 3);
if (i == j) {
- dJ(i, j, N - 1) += 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
+ dJ(i, j, N - 1) += 4 * sqrt(N) * ζ(i) / pow(b, 2);
}
}
}
@@ -103,13 +103,13 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con
grad = gInvJac * gradZ;
- Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac.conjugate());
+ Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac);
Vector df = jacobian * gradZ;
Eigen::Tensor<Scalar, 1> dfT = Eigen::TensorMap<Eigen::Tensor<Scalar, 1>>(df.data(), {df.size()});
- std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(0, 0)};
+ std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(1, 0)};
Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip);
- hess = jacobian * hessZ * jacobian.transpose() + (Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size())).transpose();
+ hess = jacobian * hessZ * jacobian.transpose() + Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size());
return {hamiltonian, grad, hess};
}