summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--langevin.cpp137
-rw-r--r--p-spin.hpp39
2 files changed, 134 insertions, 42 deletions
diff --git a/langevin.cpp b/langevin.cpp
index e5b370c..dc7c5cf 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -41,12 +41,9 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
complex_normal_distribution<Real> dJ(0, σ, 0);
ComplexTensor J = J0;
- Complex H1, H2;
ComplexVector z1, z2;
- std::tie(H1, std::ignore, std::ignore) = hamGradHess(J, z10);
- std::tie(H2, std::ignore, std::ignore) = hamGradHess(J, z20);
- Real prevdist = abs(imag(H1-H2));
- Real γ = prevdist / 10;
+ Real ΔH = abs(imag(getHamiltonian(J, z10) - getHamiltonian(J, z20))) / z10.size();
+ Real γ = ΔH;
std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
[&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
@@ -54,7 +51,7 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
setJ<Complex, p>(JJ, ind, Ji + γ * dJ(r.engine()));
};
- while (true) {
+ while (ΔH > ε) {
ComplexTensor Jp = J;
iterateOver<Complex, p>(Jp, perturbJ);
@@ -62,18 +59,17 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
z1 = findSaddle(Jp, z10, Δ);
z2 = findSaddle(Jp, z20, Δ);
- Real dist = (z1 - z2).norm();
+ Real Δz = (z1 - z2).norm();
- std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, z1);
- std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, z2);
+ if (Δz > 1e-2) {
+ Real ΔHNew = abs(imag(getHamiltonian(Jp, z1) - getHamiltonian(Jp, z2))) / z1.size();
- if (abs(imag(H1 - H2)) < prevdist && dist > Real(1) / 100) {
- J = Jp;
- prevdist = abs(imag(H1 - H2));
- γ = prevdist / 10;
+ if (ΔHNew < ΔH) {
+ J = Jp;
+ ΔH = ΔHNew;
+ γ = ΔH;
- if (abs(imag(H1 - H2)) < ε && dist > Real(1) / 100) {
- break;
+ std::cerr << "Match imaginary energies: Found couplings with ΔH = " << ΔH << std::endl;
}
}
} catch (std::exception& e) {}
@@ -82,20 +78,19 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
return {J, z1, z2};
}
-
int main(int argc, char* argv[]) {
// model parameters
unsigned N = 10; // number of spins
Real T = 1; // temperature
- Real Rκ = 1; // real part of distribution parameter
+ Real Rκ = 0; // real part of distribution parameter
Real Iκ = 0; // imaginary part of distribution parameter
// simulation parameters
- Real ε = 1e-12;
+ Real ε = 1e-15;
Real εJ = 1e-5;
Real δ = 1e-2; // threshold for determining saddle
Real Δ = 1e-3;
- Real γ = 1e-2; // step size
+ Real γ = 1e-1; // step size
unsigned t = 1000; // number of Langevin steps
unsigned M = 100;
unsigned n = 100;
@@ -146,7 +141,6 @@ int main(int argc, char* argv[]) {
complex_normal_distribution<Real> d(0, 1, 0);
ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine());
- ComplexVector z = normalize(randomVector<Complex>(N, d, r.engine()));
std::function<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
(const ComplexTensor& J, const ComplexVector& z) {
@@ -155,45 +149,104 @@ int main(int argc, char* argv[]) {
return W;
};
- std::unordered_map<uint64_t, ComplexVector> saddles;
- std::list<std::array<ComplexVector, 2>> nearbySaddles;
- while (true) { // Until we find two saddles sufficiently close...
- std::tie(std::ignore, z) = metropolis(J, z, energyNormGrad, T, γ, M, d, r.engine());
- try {
- ComplexVector zSaddleNext = findSaddle(J, z, ε);
- uint64_t saddleKey = round(1e2 * real(zSaddleNext(0)));
- auto got = saddles.find(saddleKey);
-
- if (got == saddles.end()) {
- saddles[saddleKey] = zSaddleNext;
- nearbySaddles = saddlesCloserThan(saddles, δ);
- if (nearbySaddles.size() > 0) {
- break;
- }
- std::cerr << "Found " << saddles.size() << " distinct saddles." << std::endl;
+ std::function<Real(const ComplexTensor&, const ComplexVector&)> energyThres = [N, κ]
+ (const ComplexTensor& J, const ComplexVector& z) {
+ Real a = z.squaredNorm() / (Real)N;
+ /*
+ Complex ε = getHamiltonian(J, z) / (Real)N;
+ Real t = norm(ε) / getThresholdEnergyDensity(p, κ, ε, a);
+ */
+
+ Real amin = 1.1547;
+ Real amax = 1.46806;
+
+ /*
+ Real tmin = 1;
+ Real tmax = 1.17514;
+ */
+
+ Real E = 0;
+
+ if (a > amax) {
+ E += pow(a - amax, 2);
+ } else if (a < amin) {
+ E += pow(a - amin, 2);
}
- } catch (std::exception& e) {
-// std::cerr << "Could not find a saddle: " << e.what() << std::endl;
+
+ /*
+ if (t > tmax) {
+ E += pow(t - tmax, 2);
+ } else if (t < tmin) {
+ E += pow(t - tmin, 2);
+ }
+ */
+
+ return E;
+ };
+
+ Real largestThreshold = 0;
+ ComplexVector saddlePastThreshold;
+ bool foundSaddle = false;
+
+#pragma omp parallel default(none) shared(largestThreshold, saddlePastThreshold, foundSaddle, M, J, energyThres, d, N, γ, ε, κ)
+ {
+ Rng rPrivate;
+ ComplexVector z = normalize(randomVector<Complex>(N, d, rPrivate.engine()));
+
+ while (largestThreshold < 1) { // Until we find a saddle past the threshold...
+ std::tie(std::ignore, z) = metropolis(J, z, energyThres, (Real)0.1, γ, M, d, rPrivate.engine());
+ try {
+ ComplexVector latestSaddle = findSaddle(J, z, ε);
+ Real threshold = getProportionOfThreshold(κ, J, latestSaddle);
+
+#pragma omp critical
+ {
+ if (threshold > largestThreshold + 1e-6) {
+ largestThreshold = threshold;
+ saddlePastThreshold = latestSaddle;
+ std::cerr << "Found saddle with threshold porportion " << largestThreshold << std::endl;;
+ }
+ }
+ } catch (std::exception& e) {}
}
+ }
+
+ std::cerr << "Found saddle with energy beyond threshold." << std::endl;
+
+ ComplexVector zSaddleNext;
+ Real closestSaddle = std::numeric_limits<Real>::infinity();
+ ComplexVector z = saddlePastThreshold;
+ while (closestSaddle > δ) { // Until we find two saddles sufficiently close...
+ std::tie(std::ignore, z) = metropolis(J, saddlePastThreshold, energyNormGrad, T, γ, 1e4, d, r.engine());
+ try {
+ zSaddleNext = findSaddle(J, z, ε);
+ Real saddleDistance = (zSaddleNext - saddlePastThreshold).norm();
+ if (saddleDistance < closestSaddle && saddleDistance > 1e-2) {
+ closestSaddle = saddleDistance;
+ std::cerr << "Nearby saddle: found saddle at distance " << saddleDistance << std::endl;
+ }
+ } catch (std::exception& e) {}
}
std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl;
- ComplexVector z1 = nearbySaddles.front()[0];
- ComplexVector z2 = nearbySaddles.front()[1];
+ ComplexVector z1 = saddlePastThreshold;
+ ComplexVector z2 = zSaddleNext;
- std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-15, ε, r);
+ std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r);
std::cerr << "Im H is now sufficently close, starting to relax rope." << std::endl;
+ std::cerr << "Threshold proportions of saddles are " << getProportionOfThreshold(κ, J, z1) << " and " << getProportionOfThreshold(κ, J, z2) << std::endl;
+
Rope stokes(10, z1, z2, J);
for (unsigned i = 0; i < 9; i++) {
stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0);
- std::cout << stokes.n() << " " << stokes.cost(J) << std::endl;
+ std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << std::endl;
stokes = stokes.interpolate();
}
diff --git a/p-spin.hpp b/p-spin.hpp
index f1dc07f..b2bdf07 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -26,6 +26,45 @@ std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scal
return {hamiltonian, gradient, hessian};
}
+template <class Scalar, int p>
+Scalar getHamiltonian(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+ Scalar H;
+ std::tie(H, std::ignore, std::ignore) = hamGradHess(J, z);
+ return H;
+}
+
+template <class Scalar, int p>
+Vector<Scalar> getGradient(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+ Vector<Scalar> dH;
+ std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z);
+ return dH;
+}
+
+template <class Scalar, int p>
+Matrix<Scalar> getHessian(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+ Matrix<Scalar> ddH;
+ std::tie(std::ignore, std::ignore, ddH) = hamGradHess(J, z);
+ return ddH;
+}
+
+template <class Scalar>
+Real getThresholdEnergyDensity(unsigned p, Scalar κ, Scalar ε, Real a) {
+ Real apm2 = pow(a, p - 2);
+ Scalar δ = κ / apm2;
+ Real θ = arg(κ) + 2 * arg(ε);
+
+ return (p - 1) * apm2 / (2 * p) * pow(1 - norm(δ), 2) / (1 + norm(δ) - 2 * abs(δ) * cos(θ));
+}
+
+template <class Scalar, int p>
+Real getProportionOfThreshold(Scalar κ, const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+ Real N = z.size();
+ Scalar ε = getHamiltonian(J, z) / N;
+ Real a = z.squaredNorm() / N;
+
+ return norm(ε) / getThresholdEnergyDensity(p, κ, ε, a);
+}
+
template <class Scalar>
Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate();