From 61ee450ddf82c6206ac3e55dd5e064d5002e57fe Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
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<Real>;
template<int ...ps, class Generator, typename... T>
-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<Vector<Complex>> 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<Scalar, 3> ddzDot(const Vector<Scalar>& z, const Vector<Scalar>& 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 <class Scalar>
@@ -172,7 +172,7 @@ Tensor<Scalar, 3> dcdzDot(const Vector<Scalar>& z, const Vector<Scalar>& 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<int, 3>){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<Scalar, 3> ddzDotConjugate(const Vector<Scalar>& z, const Vector<Scalar>&
- zT.conjugate().contract(IT, ei) * (z.dot(dH) / pow(z², 2))
- IT.contract(zT.conjugate(), ei).shuffle((std::array<int, 3>){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 <class Scalar, int ...ps>
@@ -44,7 +44,6 @@ public:
std::vector<Vector<Scalar>> gs;
Vector<Scalar> z0;
Vector<Scalar> z1;
- double φ;
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);
@@ -148,7 +147,7 @@ public:
for (unsigned i = 0; i < gs.size(); i++) {
Real fdgn = gCoeff(i, t);
Real dfdgn = dgCoeff(i, t);
- Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * (
+ Vector<Scalar> 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<Scalar> dCm = - 0.5 / ż.norm() / dz.norm() * (
+ Vector<Scalar> 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<Scalar> ddC = 0.5 / ż.norm() / dz.norm() * (
+ Matrix<Scalar> 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<Scalar> dcdC = 0.5 / ż.norm() / dz.norm() * (
+ Matrix<Scalar> 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