diff options
Diffstat (limited to 'stokes.hpp')
| -rw-r--r-- | stokes.hpp | 372 |
1 files changed, 211 insertions, 161 deletions
@@ -4,6 +4,8 @@ #include "complex_normal.hpp" #include "dynamics.hpp" +#include <iostream> + class ropeRelaxationStallException: public std::exception { virtual const char* what() const throw() { return "Gradient descent stalled."; @@ -11,49 +13,6 @@ class ropeRelaxationStallException: public std::exception { }; 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) { - Real z² = z.squaredNorm(); - Real z´² = z´.squaredNorm(); - - Vector<Scalar> ż = zDot(z, dH); - Real ż² = ż.squaredNorm(); - - Real Reż·z´ = real(ż.dot(z´)); - - 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 = - ( - dżc * z´ + dż * z´.conjugate() - (dż * ż.conjugate() + dżc * ż) * Reż·z´ / ż² - ) / sqrt(ż² * z´²) / 2; - - 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²; - - Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´)); - - Vector<Scalar> ddtdLdz´ = - ( - ( - ż´.conjugate() - ( - Reż·z´ * z´´.conjugate() + dReż·z´ * z´.conjugate() - - (Reż·z´ / z´²) * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() - ) / z´² - ) - - 0.5 * ( - (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´² - ) * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate()) - ) / sqrt(ż² * z´²) / 2; - - return dLdz - ddtdLdz´; -} - -template <class Scalar> class Rope { public: std::vector<Vector<Scalar>> z; @@ -93,65 +52,10 @@ class Rope { return l; } - template <int p> - Real error(const Tensor<Scalar, p>& J) const { - Scalar H0, HN; - std::tie(H0, std::ignore, std::ignore) = hamGradHess(J, z.front()); - std::tie(HN, std::ignore, std::ignore) = hamGradHess(J, z.back()); - - Real ImH = imag((H0 + HN) / 2.0); - - Real err = 0; - - for (unsigned i = 1; i < z.size() - 1; i++) { - Scalar Hi; - std::tie(Hi, std::ignore, std::ignore) = hamGradHess(J, z[i]); - - err += pow(imag(Hi) - ImH, 2); - } - - return sqrt(err); - } - Vector<Scalar> dz(unsigned i) const { return z[i + 1] - z[i - 1]; } - Vector<Scalar> ddz(unsigned i) const { - return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]); - } - - template <int p> - std::vector<Vector<Scalar>> generateGradientδz(const Tensor<Scalar, p>& J) const { - std::vector<Vector<Scalar>> δz(z.size()); - -#pragma omp parallel for - 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]); - - δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH); - } - - for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate(); -// δz[i] = δz[i] - ((δz[i].conjugate().dot(dz(i))) / dz(i).squaredNorm()) * dz(i).conjugate(); - } - - // We return a δz with average norm of one. - Real mag = 0; - for (unsigned i = 1; i < z.size() - 1; i++) { - mag += δz[i].norm(); - } - - for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] /= mag / n(); - } - - return δz; - } - template <int p> std::vector<Vector<Scalar>> generateDiscreteGradientδz(const Tensor<Scalar, p>& J, Real γ) const { std::vector<Vector<Scalar>> δz(z.size()); @@ -185,29 +89,38 @@ class Rope { dC += 0.5 * (ż[i + 1].conjugate() - dz(i + 1).conjugate() * real(ż[i + 1].dot(dz(i + 1))) / dz(i + 1).squaredNorm()) / (ż[i + 1].norm() * dz(i + 1).norm()); } - dC += - γ * (z[i - 1] + z[i + 1]).conjugate(); + dC += γ * (2 * z[i] - z[i - 1] - z[i + 1]).conjugate(); - δz[i] = dC; - } + δz[i] = dC.conjugate(); - Real size = 0; - for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate(); + δz[i] -= z[i].conjugate() * z[i].conjugate().dot(δz[i]) / z²; } return δz; } - template<class Gen> - std::vector<Vector<Scalar>> generateRandomδz(Gen& r) const { - std::vector<Vector<Scalar>> δz(z.size()); + void spread() { + Real l = length(); + + Real a = 0; + unsigned pos = 0; + + std::vector<Vector<Scalar>> zNew = z; - complex_normal_distribution<> d(0, 1, 0); for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] = randomVector<Scalar>(z[0].size(), d, r); + Real b = i * l / (z.size() - 1); + + while (b > a) { + pos++; + a += (z[pos] - z[pos - 1]).norm(); + } + + Vector<Scalar> δz = z[pos] - z[pos - 1]; + + zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz); } - return δz; + z = zNew; } template<int p> @@ -222,8 +135,6 @@ class Rope { rNew.z[i] = normalize(z[i] - (δ * Δl) * δz[i]); } - rNew.spread(); - if (rNew.cost(J, γ) < cost(J, γ)) { break; } else { @@ -235,71 +146,33 @@ class Rope { } } +// rNew.spread(); + z = rNew.z; return δ; } - void spread() { - Real l = length(); - - Real a = 0; - unsigned pos = 0; - - std::vector<Vector<Scalar>> zNew = z; - - for (unsigned i = 1; i < z.size() - 1; i++) { - Real b = i * l / (z.size() - 1); - - while (b > a) { - pos++; - a += (z[pos] - z[pos - 1]).norm(); - } - - Vector<Scalar> δz = z[pos] - z[pos - 1]; - - zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz); - } - - z = zNew; - } - - template <int p> - void relaxGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0) { - Real δ = δ0; - try { - for (unsigned i = 0; i < N; i++) { - std::vector<Vector<Scalar>> δz = generateGradientδz(J); - δ = 1.1 * perturb(J, δ, δz); - } - } catch (std::exception& e) { - } - } - template <int p> void relaxDiscreteGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Real γ) { Real δ = δ0; try { for (unsigned i = 0; i < N; i++) { std::vector<Vector<Scalar>> δz = generateDiscreteGradientδz(J, γ); - δ = 1.1 * perturb(J, δ, δz, γ); + double stepSize = 0; + for (const Vector<Scalar>& v : δz) { + stepSize += v.norm(); + } + if (stepSize / δz.size() < 1e-6) { + break; + } + std::cout << cost(J) << " " << stepSize / δz.size() << std::endl; + δ = 2 * perturb(J, δ, δz, γ); } } catch (std::exception& e) { } } - template <int p, class Gen> - void relaxRandom(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Gen& r) { - Real δ = δ0; - for (unsigned i = 0; i < N; i++) { - try { - std::vector<Vector<Scalar>> δz = generateRandomδz(r); - δ = 1.1 * perturb(J, δ, δz); - } catch (std::exception& e) { - } - } - } - template <int p> Real cost(const Tensor<Scalar, p>& J, Real γ = 0) const { Real c = 0; @@ -334,3 +207,180 @@ class Rope { return r; } }; + +template <class Real, class Scalar, int p> +bool stokesLineTest(const Tensor<Scalar, p>& J, const Vector<Scalar>& z1, const Vector<Scalar>& z2, unsigned n0, unsigned steps) { + Rope stokes(n0, z1, z2, J); + + Real oldCost = stokes.cost(J); + + for (unsigned i = 0; i < steps; i++) { + stokes.relaxDiscreteGradient(J, 1e6, 1, pow(2, steps)); + + Real newCost = stokes.cost(J); + + if (newCost > oldCost) { + return false; + } + + oldCost = newCost; + + stokes = stokes.interpolate(); + } + return true; +} + +template <class Scalar> +class Cord { +public: + std::vector<Vector<Scalar>> gs; + Vector<Scalar> z0; + Vector<Scalar> z1; + + template <int p> + Cord(const Tensor<Scalar, p>& J, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : gs(ng, Vector<Scalar>::Zero(z2.size())) { + Scalar H2 = getHamiltonian(J, z2); + Scalar H3 = getHamiltonian(J, z3); + + if (real(H2) > real(H3)) { + z0 = z2; + z1 = z3; + } else { + z0 = z3; + z1 = z2; + } + } + + Real gCoeff(unsigned i, Real t) const { + return (1 - t) * t * pow(t, i); + } + + Real dgCoeff(unsigned i, Real t) const { + return (i + 1) * (1 - t) * pow(t, i) - pow(t, i + 1); + } + + Vector<Scalar> f(Real t) const { + Vector<Scalar> z = (1 - t) * z0 + t * z1; + + for (unsigned i = 0; i < gs.size(); i++) { + z += gCoeff(i, t) * gs[i]; + } + + return z; + } + + Vector<Scalar> df(Real t) const { + Vector<Scalar> z = z1 - z0; + + for (unsigned i = 0; i < gs.size(); i++) { + z += dgCoeff(i, t) * gs[i]; + } + + return z; + } + + template <int p> + Real cost(const Tensor<Scalar, p>& J, Real t) const { + Vector<Scalar> z = f(t); + Scalar H; + Vector<Scalar> dH; + std::tie(H, dH, std::ignore) = hamGradHess(J, z); + Vector<Scalar> ż = zDot(z, dH); + Vector<Scalar> dz = df(t); + + return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm(); + } + + template <int p> + Real totalCost(const Tensor<Scalar, p>& J, unsigned nt) const { + Real tc = 0; + + for (unsigned i = 0; i < nt; i++) { + Real t = (i + 1.0) / (nt + 1.0); + tc += cost(J, t); + } + + return tc; + } + + template <int p> + std::vector<Vector<Scalar>> dgs(const Tensor<Scalar, p>& J, Real t) const { + Vector<Scalar> z = f(t); + auto [H, dH, ddH] = hamGradHess(J, z); + Vector<Scalar> ż = zDot(z, dH); + Vector<Scalar> dz = df(t); + Matrix<Scalar> dż = dzDot(z, dH); + Matrix<Scalar> dżc = dzDotConjugate(z, dH, ddH); + + std::vector<Vector<Scalar>> x; + x.reserve(gs.size()); + + for (unsigned i = 0; i < gs.size(); i++) { + Real fdg = gCoeff(i, t); + Real dfdg = dgCoeff(i, t); + Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * ( + dfdg * ż.conjugate() + fdg * dżc * dz + fdg * dż * dz.conjugate() + - real(dz.dot(ż)) * ( + dfdg * dz.conjugate() / dz.squaredNorm() + + fdg * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() + ) + ); + + x.push_back(dC.conjugate()); + } + + return x; + } + + template <int p> + std::pair<Real, Real> relaxStep(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀) { + std::vector<Vector<Scalar>> dgsTot(gs.size(), Vector<Scalar>::Zero(z0.size())); + + for (unsigned i = 0; i < nt; i++) { + Real t = (i + 1.0) / (nt + 1.0); + std::vector<Vector<Scalar>> dgsI = dgs(J, t); + + for (unsigned j = 0; j < gs.size(); j++) { + dgsTot[j] += dgsI[j] / nt; + } + } + + Real stepSize = 0; + for (const Vector<Scalar>& dgi : dgsTot) { + stepSize += dgi.squaredNorm(); + } + stepSize = sqrt(stepSize); + + Cord cNew(*this); + + Real δ = δ₀; + Real oldCost = totalCost(J, nt); + Real newCost = std::numeric_limits<Real>::infinity(); + + while (newCost > oldCost) { + for (unsigned i = 0; i < gs.size(); i++) { + cNew.gs[i] = gs[i] - δ * dgsTot[i]; + } + + newCost = cNew.totalCost(J, nt); + + δ /= 2; + } + + gs = cNew.gs; + + return {2 * δ, stepSize}; + } + + template <int p> + void relax(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀, unsigned maxSteps) { + Real δ = δ₀; + Real stepSize = std::numeric_limits<Real>::infinity(); + unsigned steps = 0; + while (stepSize > 1e-7 && steps < maxSteps) { + std::tie(δ, stepSize) = relaxStep(J, nt, δ); + std::cout << totalCost(J, nt) << " " << δ << " " << stepSize << std::endl; + steps++; + } + } +}; |