summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--stokes.hpp52
1 files changed, 26 insertions, 26 deletions
diff --git a/stokes.hpp b/stokes.hpp
index 117f4de..abf0c50 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,11 +1,5 @@
#include "p-spin.hpp"
-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();
-}
-
class ropeRelaxationStallException: public std::exception {
virtual const char* what() const throw() {
return "Gradient descent stalled.";
@@ -61,6 +55,10 @@ class Rope {
public:
std::vector<Vector<Scalar>> z;
+ unsigned n() const {
+ return z.size() - 2;
+ }
+
Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2) : z(N + 2) {
for (unsigned i = 0; i < N + 2; i++) {
z[i] = normalize(z1 + (z2 - z1) * ((double)i / (N + 1.0)));
@@ -85,7 +83,7 @@ class Rope {
Matrix<Scalar> ddH;
std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]);
- δz[i] = variation(z[i], this->dz(i), this->ddz(i), dH, ddH);
+ δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH);
}
return δz;
@@ -93,12 +91,12 @@ class Rope {
template <int p>
double perturb(const Tensor<Scalar, p>& J, double δ0) {
- std::vector<Vector<Scalar>> Δz = this->δz(J);
+ std::vector<Vector<Scalar>> Δz = δz(J);
Rope rNew = *this;
double δ = δ0;
- while (rNew.cost(J) >= this->cost(J)) {
+ while (rNew.cost(J) >= cost(J)) {
for (unsigned i = 1; i < z.size() - 1; i++) {
rNew.z[i] = normalize(z[i] - δ * Δz[i].conjugate());
}
@@ -144,11 +142,22 @@ class Rope {
}
template <int p>
- double step(const Tensor<Scalar, p>& J, double δ) {
- double δNew = this->perturb(J, δ);
- this->spread();
+ double step(const Tensor<Scalar, p>& J, double δ0) {
+ double δ = perturb(J, δ0);
+ spread();
- return δNew;
+ return δ;
+ }
+
+ template <int p>
+ void relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) {
+ double δ = δ0;
+ try {
+ for (unsigned i = 0; i < N; i++) {
+ δ = 2 * step(J, δ);
+ }
+ } catch (std::exception& e) {
+ }
}
template <int p>
@@ -159,25 +168,16 @@ class Rope {
Vector<Scalar> dH;
std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z[i]);
- c += segmentCost(z[i], this->dz(i), dH);
+ Vector<Scalar> zD = zDot(z[i], dH);
+
+ c += 1.0 - pow(real(zD.dot(dz(i))), 2) / zD.squaredNorm() / dz(i).squaredNorm();
}
return c;
}
- template <int p>
- void relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) {
- double δ = δ0;
- try {
- for (unsigned i = 0; i < N; i++) {
- δ = 2 * this->step(J, δ);
- }
- } catch (std::exception& e) {
- }
- }
-
Rope<Scalar> interpolate() const {
- Rope<Scalar> r(2 * (z.size() - 2) + 1, z[0], z[z.size() - 1]);
+ Rope<Scalar> r(2 * n() + 1, z[0], z[z.size() - 1]);
for (unsigned i = 0; i < z.size(); i++) {
r.z[2 * i] = z[i];