From 61ee450ddf82c6206ac3e55dd5e064d5002e57fe Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 23 Nov 2021 16:53:17 +0100 Subject: Fixed some problems with changing the Real type, and went back to geometric increases in fit. --- collectStokesData.hpp | 8 ++++---- p-spin.hpp | 6 +++--- stokes.hpp | 23 +++++++++++------------ 3 files changed, 18 insertions(+), 19 deletions(-) diff --git a/collectStokesData.hpp b/collectStokesData.hpp index 42500c0..b48c6bc 100644 --- a/collectStokesData.hpp +++ b/collectStokesData.hpp @@ -5,8 +5,8 @@ using Complex = std::complex; template -void collectStokesData(std::string tag, unsigned N, Generator& r, double ε, double δz, bool minimum, T... μs) { - unsigned nIts = 7; +void collectStokesData(std::string tag, unsigned N, Generator& r, Real ε, Real δz, bool minimum, T... μs) { + unsigned nIts = 5; unsigned nGs = 4; unsigned coDim = 16; Real newSaddleThres = 1e-4; @@ -83,11 +83,11 @@ void collectStokesData(std::string tag, unsigned N, Generator& r, double ε, dou std::vector> oldGs; for (int j = 0; j < nIts; j++) { - Cord c(M, zMin, zSaddle, nGs + j * 2); + Cord c(M, zMin, zSaddle, nGs + pow(2, j)); for (unsigned i = 0; i < oldGs.size(); i++) { c.gs[i] = oldGs[i]; } - c.relaxNewton(c.gs.size() * 4, 1, 1e-10, 1e4); + c.relaxNewton(c.gs.size() * 2, 1, 1e-10, 1e4); oldGs = c.gs; Real reConstraintError = 0; diff --git a/p-spin.hpp b/p-spin.hpp index 544abf1..dee6273 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -150,7 +150,7 @@ Tensor ddzDot(const Vector& z, const Vector& dH) { return - zT.conjugate().contract(dHT.conjugate(), ei).contract(zT.conjugate(), ei) / pow(z², 2) - dHT.conjugate().contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) / pow(z², 2) - + zT.conjugate().contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * (2.0 * dH.dot(z) / pow(z², 3)); + + zT.conjugate().contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * ((Real)2 * dH.dot(z) / pow(z², 3)); } template @@ -172,7 +172,7 @@ Tensor dcdzDot(const Vector& z, const Vector& dH, con -ddHT.conjugate().contract(zT, ip00).contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) / pow(z², 2) -IT.contract(zT.conjugate(), ei) * (dH.dot(z) / pow(z², 2)) -IT.contract(zT.conjugate(), ei).shuffle((std::array){0,2,1}) * (dH.dot(z) / pow(z², 2)) - +zT.contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * (2.0 * dH.dot(z) / pow(z², 3)) + +zT.contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * ((Real)2 * dH.dot(z) / pow(z², 3)) ; } @@ -197,5 +197,5 @@ Tensor ddzDotConjugate(const Vector& z, const Vector& - zT.conjugate().contract(IT, ei) * (z.dot(dH) / pow(z², 2)) - IT.contract(zT.conjugate(), ei).shuffle((std::array){0, 2, 1}) * (z.dot(dH) / pow(z², 2)) - ddHT.contract(zT.conjugate(), ip00).contract(zT.conjugate(), ei).contract(zT, ei) / pow(z², 2) - + zT.conjugate().contract(zT.conjugate(), ei).contract(zT, ei) * (2.0 * z.dot(dH) / pow(z², 3)); + + zT.conjugate().contract(zT.conjugate(), ei).contract(zT, ei) * ((Real)2 * z.dot(dH) / pow(z², 3)); } diff --git a/stokes.hpp b/stokes.hpp index 76b950f..bcb6608 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -34,7 +34,7 @@ Real dJacobi(unsigned α, unsigned β, unsigned n, Real x) { f += pow(y, n - s) * pow(z, s - 1) / (tgamma(s - 1 + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s + 1)); } - return 0.5 * tgamma(n + α + 1) * tgamma(n + β + 1) * f; + return tgamma(n + α + 1) * tgamma(n + β + 1) * f / 2; } template @@ -44,7 +44,6 @@ public: std::vector> gs; Vector z0; Vector z1; - double φ; Cord(const pSpinModel& M, const Vector& z2, const Vector& z3, unsigned ng) : M(M), gs(ng, Vector::Zero(z2.size())) { Scalar H2 = M.getHamiltonian(z2); @@ -148,7 +147,7 @@ public: for (unsigned i = 0; i < gs.size(); i++) { Real fdgn = gCoeff(i, t); Real dfdgn = dgCoeff(i, t); - Vector dC = - 0.5 / ż.norm() / dz.norm() * ( + Vector dC = - 1 / ż.norm() / dz.norm() / 2 * ( dfdgn * ż.conjugate() + fdgn * dżc * dz + fdgn * dż * dz.conjugate() - std::real(dz.dot(ż)) * ( dfdgn * dz.conjugate() / dz.squaredNorm() + @@ -163,7 +162,7 @@ public: Real fdgm = gCoeff(j, t); Real dfdgm = dgCoeff(j, t); Scalar CC = - std::real(dz.dot(ż)) / ż.norm() / dz.norm(); - Vector dCm = - 0.5 / ż.norm() / dz.norm() * ( + Vector dCm = - 1 / ż.norm() / dz.norm() / 2 * ( dfdgm * ż.conjugate() + fdgm * dżc * dz + fdgm * dż * dz.conjugate() - std::real(dz.dot(ż)) * ( dfdgm * dz.conjugate() / dz.squaredNorm() + @@ -171,7 +170,7 @@ public: ) ); - Matrix ddC = 0.5 / ż.norm() / dz.norm() * ( + Matrix ddC = 1 / ż.norm() / dz.norm() / 2 * ( - ( dfdgn * fdgm * dżc.transpose() + dfdgm * fdgn * dżc + fdgn * fdgm * ddżcdz + fdgn * fdgm * ddżdz @@ -187,18 +186,18 @@ public: - CC * dz.norm() / ż.norm() * fdgn * fdgm * ( ddżżc + ddżcż + dżc * dż.transpose() + dż * dżc.transpose() ) - - 0.5 * CC / dz.norm() / ż.norm() * ( + - CC / dz.norm() / ż.norm() / (Real)2 * ( dfdgn * fdgm * dz.conjugate() * (dżc * ż + dż * ż.conjugate()).transpose() + dfdgm * fdgn * (dżc * ż + dż * ż.conjugate()) * dz.adjoint() ) - + 0.5 * CC * ż.norm() / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.adjoint() - + 0.5 * CC * dz.norm() / pow(ż.norm(), 3) + + CC * ż.norm() / (Real)2 / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.adjoint() + + CC * dz.norm() / (Real)2 / pow(ż.norm(), 3) * fdgn * (dżc * ż + dż * ż.conjugate()) * fdgm * (dżc * ż + dż * ż.conjugate()).transpose() ) ; - Matrix dcdC = 0.5 / ż.norm() / dz.norm() * ( + Matrix dcdC = 1 / ż.norm() / dz.norm() / 2 * ( - ( dfdgn * fdgm * dż.adjoint() + dfdgm * fdgn * dż + fdgn * fdgm * (dcdżcdz + dcdżcdz.adjoint()) @@ -215,12 +214,12 @@ public: - CC * dz.norm() / ż.norm() * fdgn * fdgm * ( dcdżcż + dcdżcż.adjoint() + dżc * dżc.adjoint() + dż * dż.adjoint() ) - - 0.5 * CC / dz.norm() / ż.norm() * ( + - CC / dz.norm() / ż.norm() / (Real)2 * ( dfdgn * fdgm * dz.conjugate() * (dżc * ż + dż * ż.conjugate()).adjoint() + dfdgm * fdgn * (dżc * ż + dż * ż.conjugate()) * dz.transpose() ) - + 0.5 * CC * ż.norm() / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.transpose() - + 0.5 * CC * dz.norm() / pow(ż.norm(), 3) * fdgn * fdgm + + CC * ż.norm() / (Real)2 / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.transpose() + + CC * dz.norm() / (Real)2 / pow(ż.norm(), 3) * fdgn * fdgm * (dżc * ż + dż * ż.conjugate()) * (dżc * ż + dż * ż.conjugate()).adjoint() ) -- cgit v1.2.3-70-g09d2