summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--langevin.cpp23
-rw-r--r--p-spin.hpp109
-rw-r--r--tensor.hpp24
3 files changed, 77 insertions, 79 deletions
diff --git a/langevin.cpp b/langevin.cpp
index 8ecd69c..0e5a616 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -81,17 +81,10 @@ int main(int argc, char* argv[]) {
}
}
- Complex κ(Rκ, Iκ);
- Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));
- Real σ2 = sqrt(factorial(2) / ((Real)2 * pow(N, 2 - 1)));
- Real σ4 = sqrt(factorial(4) / ((Real)2 * pow(N, 4 - 1)));
-
Rng r;
- pSpinModel<Real, 2, 4> pSpin;
-
- std::get<0>(pSpin.Js) = generateCouplings<Real, 2>(N, std::normal_distribution<Real>(0, σ2), r.engine());
- std::get<1>(pSpin.Js) = generateCouplings<Real, 4>(N, std::normal_distribution<Real>(0, σ4 / 10), r.engine());
+// pSpinModel<Real, 2, 4> pSpin({J2, J4});
+ pSpinModel<Real, 2, 4>pSpin(N, r.engine(), 1, 0.1);
std::normal_distribution<Real> Red(0, 1);
@@ -100,18 +93,17 @@ int main(int argc, char* argv[]) {
RealVector dHr;
RealMatrix ddHr;
std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin);
- Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + 2.0 * zMin * zMin.transpose());
- std::cout << eigenS.eigenvalues().transpose() << std::endl;
+ Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + zMin * zMin.transpose() / (Real)N);
for (unsigned i = 0; i < N; i++) {
- RealVector zNew = normalize(zMin + 0.01 * eigenS.eigenvectors().col(i).real());
- std::cout << pSpin.getHamiltonian(zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl;
+ RealVector zNew = normalize(zMin + 1e-3 * eigenS.eigenvectors().col(i).real());
+ std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-6 << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << " " << eigenS.eigenvectors().col(0).dot(eigenS.eigenvectors().col(i)) << std::endl;
}
std::cout << std::endl;
getchar();
complex_normal_distribution<Real> d(0, 1, 0);
- pSpinModel<Complex, 2, 4> complexPSpin = pSpin.cast<Complex>();;
+ pSpinModel complexPSpin = pSpin.cast<Complex>();;
ComplexVector zSaddle = zMin.cast<Complex>();
ComplexVector zSaddleNext;
@@ -134,8 +126,7 @@ int main(int argc, char* argv[]) {
Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl;
- std::get<0>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<0>(complexPSpin.Js);
- std::get<1>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<1>(complexPSpin.Js);
+ complexPSpin *= exp(Complex(0, φ));
Cord test(complexPSpin, zSaddle, zSaddleNext, 3);
test.relaxNewton(10, 1, 1e4);
diff --git a/p-spin.hpp b/p-spin.hpp
index 7a49bca..544abf1 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -12,73 +12,92 @@ Vector<typename Derived::Scalar> normalize(const Eigen::MatrixBase<Derived>& z)
return z * sqrt((Real)z.size() / (typename Derived::Scalar)(z.transpose() * z));
}
-template <class Scalar>
-std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, 2>& J) {
- Tensor<Scalar, 3> J3(z.size(), z.size(), z.size());;
- J3.setZero();
- Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());
-
- return {Jz, J3};
-}
+template <typename Scalar, int... ps>
+class pSpinModel {
+private:
+ std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, 2>& J) const {
+ Tensor<Scalar, 3> J3(z.size(), z.size(), z.size());;
+ J3.setZero();
+ Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());
-template <class Scalar, int p>
-std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J) {
- Tensor<Scalar, 3> J3 = contractDown(J, z);
- Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
- Tensor<Scalar, 2> J3zT = J3.contract(zT, ip00);
- Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J3zT.data(), z.size(), z.size());
+ return {Jz, J3};
+ }
- return {Jz, J3};
-}
+ template <int p>
+ std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J) const {
+ Tensor<Scalar, 3> J3 = contractDown(J, z);
+ Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
+ Tensor<Scalar, 2> J3zT = J3.contract(zT, ip00);
+ Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J3zT.data(), z.size(), z.size());
-template <class Scalar, int p, int... ps>
-std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> hamGradHessHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J, const Tensor<Scalar, ps>& ...Js) {
- auto [Jz, J3] = hamGradTensorHelper(z, J);
-
- Vector<Scalar> Jzz = Jz * z;
- Scalar Jzzz = Jzz.transpose() * z;
-
- Real pBang = factorial(p);
-
- Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / pBang) * J3;
- Matrix<Scalar> ddH = ((p - 1) * p / pBang) * Jz;
- Vector<Scalar> dH = (p / pBang) * Jzz;
- Scalar H = Jzzz / pBang;
-
- if constexpr (sizeof...(Js) > 0) {
- auto [Hs, dHs, ddHs, dddHs] = hamGradHessHelper(z, Js...);
- H += Hs;
- dH += dHs;
- ddH += ddHs;
- dddH += dddHs;
+ return {Jz, J3};
}
- return {H, dH, ddH, dddH};
-}
+ template <int p, int... qs>
+ std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> hamGradHessHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J, const Tensor<Scalar, qs>& ...Js) const {
+ auto [Jz, J3] = hamGradTensorHelper(z, J);
+
+ Vector<Scalar> Jzz = Jz * z;
+ Scalar Jzzz = Jzz.transpose() * z;
+
+ Real pBang = factorial(p);
+
+ Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / pBang) * J3;
+ Matrix<Scalar> ddH = ((p - 1) * p / pBang) * Jz;
+ Vector<Scalar> dH = (p / pBang) * Jzz;
+ Scalar H = Jzzz / pBang;
+
+ if constexpr (sizeof...(Js) > 0) {
+ auto [Hs, dHs, ddHs, dddHs] = hamGradHessHelper(z, Js...);
+ H += Hs;
+ dH += dHs;
+ ddH += ddHs;
+ dddH += dddHs;
+ }
+
+ return {H, dH, ddH, dddH};
+ }
-template <typename Scalar, int... ps>
-class pSpinModel {
public:
std::tuple<Tensor<Scalar, ps>...> Js;
+ pSpinModel() {}
+
+ pSpinModel(const std::tuple<Tensor<Scalar, ps>...>& Js) : Js(Js) {}
+
+ template <class Generator, typename... T>
+ pSpinModel<Real>(unsigned N, Generator& r, T... μs) {
+ Js = std::make_tuple(μs * generateRealPSpinCouplings<Real, ps>(N, r)...);
+ }
+
unsigned dimension() const {
return std::get<0>(Js).dimension(0);
}
template <typename NewScalar>
pSpinModel<NewScalar, ps...> cast() const {
- pSpinModel<NewScalar, ps...> M;
- M.Js = std::apply(
+ return std::apply(
[] (const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Tensor<NewScalar, ps>...> {
return std::make_tuple(Ks.template cast<NewScalar>()...);
}, Js
- );
+ );
+ }
+
+ template <typename T>
+ pSpinModel<Scalar, ps...>& operator*=(T x) {
+ std::tuple<Tensor<Scalar, ps>...> newJs = std::apply(
+ [x] (const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Tensor<Scalar, ps>...> {
+ return std::make_tuple((x * Ks)...);
+ }, Js
+ );
+
+ Js = newJs;
- return M;
+ return *this;
}
std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> hamGradHess(const Vector<Scalar>& z) const {
- return std::apply([&z](const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> { return hamGradHessHelper(z, Ks...); }, Js);
+ return std::apply([&z, this](const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> { return hamGradHessHelper(z, Ks...); }, Js);
}
Scalar getHamiltonian(const Vector<Scalar>& z) const {
diff --git a/tensor.hpp b/tensor.hpp
index fc99042..1212bda 100644
--- a/tensor.hpp
+++ b/tensor.hpp
@@ -5,6 +5,8 @@
#include <eigen3/unsupported/Eigen/CXX11/Tensor>
+#include "factorial.hpp"
+
template <class Scalar>
using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
@@ -89,25 +91,11 @@ Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
return J;
}
-template <class Scalar, int p, class Distribution, class Generator>
-Tensor<Scalar, p> plantState(const Tensor<Scalar, p>& J, const Vector<Scalar>& z, double β) {
- Tensor<Scalar, p> JPlanted = J;
-
- std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> plant =
- [&z, β] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
- Scalar Ji = getJ<Scalar, p>(JJ, ind);
- Scalar zzz = 1;
-
- for (unsigned i : ind) {
- zzz *= z(i);
- }
-
- setJ<Scalar, p>(JJ, ind, Ji - β * zzz / pow(zzz.size(), p));
- };
-
- iterateOver<Scalar, p>(JPlanted, plant);
+template <class Real, int p, class Generator>
+Tensor<Real, p> generateRealPSpinCouplings(unsigned N, Generator& r) {
+ Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));
- return JPlanted;
+ return generateCouplings<Real, p>(N, std::normal_distribution<Real>(0, σp), r);
}
template <class Scalar>