diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-11-09 10:13:59 +0100 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-11-09 10:13:59 +0100 |
commit | 7bf6e952b53699977f5091a78f0f9f48f7b359c5 (patch) | |
tree | 4d5637c53daf33e10fc58334cb3b15b8aa07c0d8 /stokes.hpp | |
parent | d448af5010b664025c816dc2c6e383ac7bea3491 (diff) | |
download | code-7bf6e952b53699977f5091a78f0f9f48f7b359c5.tar.gz code-7bf6e952b53699977f5091a78f0f9f48f7b359c5.tar.bz2 code-7bf6e952b53699977f5091a78f0f9f48f7b359c5.zip |
Generalized code to easily allow mixed p-spins.
Diffstat (limited to 'stokes.hpp')
-rw-r--r-- | stokes.hpp | 70 |
1 files changed, 18 insertions, 52 deletions
@@ -4,17 +4,18 @@ #include "complex_normal.hpp" #include "dynamics.hpp" -template <class Scalar> +template <class Scalar, int ...ps> class Cord { public: + const pSpinModel<Scalar, ps...>& M; std::vector<Vector<Scalar>> gs; Vector<Scalar> z0; Vector<Scalar> z1; + double φ; - 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); + Cord(const pSpinModel<Scalar, ps...>& M, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : M(M), gs(ng, Vector<Scalar>::Zero(z2.size())) { + Scalar H2 = M.getHamiltonian(z2); + Scalar H3 = M.getHamiltonian(z3); if (real(H2) > real(H3)) { z0 = z2; @@ -53,64 +54,31 @@ public: return z; } - template <int p> - Real cost(const Tensor<Scalar, p>& J, Real t) const { + Real cost(Real t) const { Vector<Scalar> z = f(t); Scalar H; Vector<Scalar> dH; - std::tie(H, dH, std::ignore) = hamGradHess(J, z); + std::tie(H, dH, std::ignore, std::ignore) = M.hamGradHess(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 totalCost(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); + tc += cost(t); } return tc; } - template <int p> - std::vector<Vector<Scalar>> dgs(const Tensor<Scalar, p>& J, Real t) const { + std::pair<Vector<Scalar>, Matrix<Scalar>> gsGradHess(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<Vector<Scalar>, Matrix<Scalar>> gsGradHess(const Tensor<Scalar, p>& J, Real t) const { - Vector<Scalar> z = f(t); - auto [H, dH, ddH] = hamGradHess(J, z); - Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / (Real)factorial(p)) * J; // BUG IF P != 3 + auto [H, dH, ddH, dddH] = M.hamGradHess(z); Vector<Scalar> ż = zDot(z, dH); Tensor<Scalar, 1> żT = Eigen::TensorMap<Tensor<const Scalar, 1>>(ż.data(), {z.size()}); Vector<Scalar> dz = df(t); @@ -241,14 +209,13 @@ public: return {x, M}; } - template <int p> - std::pair<Real, Real> relaxStepNewton(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀) { + std::pair<Real, Real> relaxStepNewton(unsigned nt, Real δ₀) { Vector<Scalar> dgsTot = Vector<Scalar>::Zero(2 * z0.size() * gs.size()); Matrix<Scalar> HessTot = Matrix<Scalar>::Zero(2 * z0.size() * gs.size(), 2 * z0.size() * gs.size()); for (unsigned i = 0; i < nt; i++) { Real t = (i + 1.0) / (nt + 1.0); - auto [dgsi, Hessi] = gsGradHess(J, t); + auto [dgsi, Hessi] = gsGradHess(t); dgsTot += dgsi / nt; HessTot += Hessi / nt; @@ -259,7 +226,7 @@ public: Cord cNew(*this); Real δ = δ₀; - Real oldCost = totalCost(J, nt); + Real oldCost = totalCost(nt); Real newCost = std::numeric_limits<Real>::infinity(); Matrix<Scalar> I = abs(HessTot.diagonal().array()).matrix().asDiagonal(); @@ -273,7 +240,7 @@ public: } } - newCost = cNew.totalCost(J, nt); + newCost = cNew.totalCost(nt); δ *= 2; } @@ -283,13 +250,12 @@ public: return {δ / 2, stepSize}; } - template <int p> - void relaxNewton(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀, unsigned maxSteps) { + void relaxNewton(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) = relaxStepNewton(J, nt, δ); + std::tie(δ, stepSize) = relaxStepNewton(nt, δ); δ /= 3; steps++; } |