diff options
Diffstat (limited to 'p-spin.hpp')
-rw-r--r-- | p-spin.hpp | 39 |
1 files changed, 39 insertions, 0 deletions
@@ -1,6 +1,7 @@ #pragma once #include <eigen3/Eigen/Dense> +#include <iterator> #include "types.hpp" #include "tensor.hpp" @@ -104,3 +105,41 @@ Matrix<Scalar> dzDotConjugate(const Vector<Scalar>& z, const Vector<Scalar>& dH, Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z² ); } + +template <class Scalar> +Tensor<Scalar, 3> ddzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) { + Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()}); + Tensor<Scalar, 1> dHT = Eigen::TensorMap<Tensor<const Scalar, 1>>(dH.data(), {dH.size()}); + + Eigen::array<Eigen::IndexPair<int>, 0> ei = {}; + + Scalar z² = z.squaredNorm(); + + 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)); +} + +template <class Scalar> +Tensor<Scalar, 3> ddzDotConjugate(const Vector<Scalar>& z, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH, const Tensor<Scalar, 3>& dddH) { + Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()}); + Tensor<Scalar, 1> dHT = Eigen::TensorMap<Tensor<const Scalar, 1>>(dH.data(), {dH.size()}); + Tensor<Scalar, 2> ddHT = Eigen::TensorMap<Tensor<const Scalar, 2>>(ddH.data(), {z.size(), z.size()}); + + Matrix<Scalar> I = Matrix<Real>::Identity(z.size(), z.size()); + Tensor<Scalar, 2> IT = Eigen::TensorMap<Tensor<const Scalar, 2>>(I.data(), {z.size(), z.size()}); + + Eigen::array<Eigen::IndexPair<int>, 0> ei = {}; + Eigen::array<Eigen::IndexPair<int>, 1> ip00 = {Eigen::IndexPair<int>(0, 0)}; + + Scalar z² = z.squaredNorm(); + + return - dddH + dddH.contract(zT.conjugate(), ip00).contract(zT, ei) / z² + + IT.contract(ddHT.contract(zT.conjugate(), ip00), ei).shuffle((std::array<int, 3>){0, 2, 1}) / z² + + ddHT.contract(zT.conjugate(), ip00).contract(IT, ei) / z² + - zT.conjugate().contract(ddHT.contract(zT.conjugate(), ip00), ei).contract(zT, ei) / pow(z², 2) + - 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)); +} |