summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--dynamics.hpp18
-rw-r--r--langevin.cpp34
-rw-r--r--p-spin.hpp9
-rw-r--r--tensor.hpp49
4 files changed, 69 insertions, 41 deletions
diff --git a/dynamics.hpp b/dynamics.hpp
index 85eef71..b373659 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -92,26 +92,26 @@ Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
}
template <class Scalar, int p, class Distribution, class Generator>
-std::tuple<double, Vector<Scalar>> langevin(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double T, double γ, unsigned N, Distribution d, Generator& r) {
+std::tuple<double, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
+ std::function<double(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy,
+ double T, double γ, unsigned N, Distribution d, Generator& r) {
Vector<Scalar> z = z0;
- double W;
- std::tie(W, std::ignore) = WdW(J, z);
+ double E = energy(J, z);
std::uniform_real_distribution<double> D(0, 1);
for (unsigned i = 0; i < N; i++) {
- Vector<Scalar> zNewTmp = z + randomVector<Scalar>(z.size(), d, r);
+ Vector<Scalar> zNewTmp = z + γ * randomVector<Scalar>(z.size(), d, r);
Vector<Scalar> zNew = normalize(zNewTmp);
- double WNew;
- std::tie(WNew, std::ignore) = WdW(J, zNew);
+ double ENew = energy(J, zNew);
- if (exp((W - WNew) / T) > D(r)) {
+ if (E - ENew > T * log(D(r))) {
z = zNew;
- W = WNew;
+ E = ENew;
}
}
- return {W, z};
+ return {E, z};
}
diff --git a/langevin.cpp b/langevin.cpp
index 870879b..5ca2d72 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -83,11 +83,41 @@ int main(int argc, char* argv[]) {
ComplexVector z0 = normalize(randomVector<Complex>(N, d, r.engine()));
ComplexVector zSaddle = findSaddle(J, z0, ε);
- ComplexVector zSaddlePrev = ComplexVector::Zero(N);
ComplexVector z = zSaddle;
+ std::function<double(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
+ (const ComplexTensor& J, const ComplexVector& z) {
+ double W;
+ std::tie(W, std::ignore) = WdW(J, z);
+ 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...
- std::tie(std::ignore, z) = langevin(J, z, T, γ, M, d, r.engine());
+ std::tie(std::ignore, z) = metropolis(J, z, energyNormGrad, T, γ, M, d, r.engine());
try {
ComplexVector zSaddleNext = findSaddle(J, z, ε);
if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles.
diff --git a/p-spin.hpp b/p-spin.hpp
index 4621db6..e5d4f94 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -5,15 +5,6 @@
#include "tensor.hpp"
#include "factorial.hpp"
-template <class Scalar>
-using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
-
-template <class Scalar>
-using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
-
-template <class Scalar, int p>
-using Tensor = Eigen::Tensor<Scalar, p>;
-
template <class Scalar, int p>
std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
Matrix<Scalar> Jz = contractDown(J, z); // Contracts J into p - 2 copies of z.
diff --git a/tensor.hpp b/tensor.hpp
index 21f2a89..6a5b2c2 100644
--- a/tensor.hpp
+++ b/tensor.hpp
@@ -5,43 +5,52 @@
#include <eigen3/unsupported/Eigen/CXX11/Tensor>
+template <class Scalar>
+using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+template <class Scalar>
+using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+
+template <class Scalar, int p>
+using Tensor = Eigen::Tensor<Scalar, p>;
+
template <class Scalar, int p, std::size_t... Indices>
-Eigen::Tensor<Scalar, p> initializeJHelper(unsigned N, std::index_sequence<Indices...>) {
+Tensor<Scalar, p> initializeJHelper(unsigned N, std::index_sequence<Indices...>) {
std::array<unsigned, p> Ns;
std::fill_n(Ns.begin(), p, N);
- return Eigen::Tensor<Scalar, p>(std::get<Indices>(Ns)...);
+ return Tensor<Scalar, p>(std::get<Indices>(Ns)...);
}
template <class Scalar, int p>
-Eigen::Tensor<Scalar, p> initializeJ(unsigned N) {
+Tensor<Scalar, p> initializeJ(unsigned N) {
return initializeJHelper<Scalar, p>(N, std::make_index_sequence<p>());
}
template <class Scalar, int p, std::size_t... Indices>
-void setJHelper(Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, Scalar val, std::index_sequence<Indices...>) {
+void setJHelper(Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, Scalar val, std::index_sequence<Indices...>) {
J(std::get<Indices>(ind)...) = val;
}
template <class Scalar, int p>
-void setJ(Eigen::Tensor<Scalar, p>& J, std::array<unsigned, p> ind, Scalar val) {
+void setJ(Tensor<Scalar, p>& J, std::array<unsigned, p> ind, Scalar val) {
do {
setJHelper<Scalar, p>(J, ind, val, std::make_index_sequence<p>());
} while (std::next_permutation(ind.begin(), ind.end()));
}
template <class Scalar, int p, std::size_t... Indices>
-Scalar getJHelper(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) {
+Scalar getJHelper(const Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) {
return J(std::get<Indices>(ind)...);
}
template <class Scalar, int p>
-Scalar getJ(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind) {
+Scalar getJ(const Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind) {
return getJHelper<Scalar, p>(J, ind, std::make_index_sequence<p>());
}
template <class Scalar, int p>
-void iterateOverHelper(Eigen::Tensor<Scalar, p>& J,
- std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f,
+void iterateOverHelper(Tensor<Scalar, p>& J,
+ std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)>& f,
unsigned l, std::array<unsigned, p> is) {
if (l == 0) {
f(J, is);
@@ -61,17 +70,17 @@ void iterateOverHelper(Eigen::Tensor<Scalar, p>& J,
}
template <class Scalar, int p>
-void iterateOver(Eigen::Tensor<Scalar, p>& J, std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f) {
+void iterateOver(Tensor<Scalar, p>& J, std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)>& f) {
std::array<unsigned, p> is;
iterateOverHelper<Scalar, p>(J, f, p, is);
}
template <class Scalar, int p, class Distribution, class Generator>
-Eigen::Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
- Eigen::Tensor<Scalar, p> J = initializeJ<Scalar, p>(N);
+Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
+ Tensor<Scalar, p> J = initializeJ<Scalar, p>(N);
- std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)> setRandom =
- [&d, &r] (Eigen::Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
+ std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> setRandom =
+ [&d, &r] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
setJ<Scalar, p>(JJ, ind, d(r));
};
@@ -81,17 +90,15 @@ Eigen::Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator
}
template <class Scalar>
-Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>
-contractDown(const Eigen::Tensor<Scalar, 2>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) {
- return Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>(J.data(), z.size(), z.size());
+Matrix<Scalar> contractDown(const Tensor<Scalar, 2>& J, const Vector<Scalar>& z) {
+ return Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());
}
const std::array<Eigen::IndexPair<int>, 1> ip00 = {Eigen::IndexPair<int>(0, 0)};
template <class Scalar, int r>
-Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>
-contractDown(const Eigen::Tensor<Scalar, r>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) {
- Eigen::Tensor<Scalar, 1> zT = Eigen::TensorMap<Eigen::Tensor<const Scalar, 1>>(z.data(), {z.size()});
- Eigen::Tensor<Scalar, r - 1> Jz = J.contract(zT, ip00);
+Matrix<Scalar> contractDown(const Tensor<Scalar, r>& J, const Vector<Scalar>& z) {
+ Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
+ Tensor<Scalar, r - 1> Jz = J.contract(zT, ip00);
return contractDown(Jz, z);
}