Work to find beyond-threshold points.

1 files changed, 39 insertions, 0 deletions

diff --git a/p-spin.hpp b/p-spin.hpp
index f1dc07f..b2bdf07 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -26,6 +26,45 @@ std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scal
   return {hamiltonian, gradient, hessian};
 }
 
+template <class Scalar, int p>
+Scalar getHamiltonian(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+  Scalar H;
+  std::tie(H, std::ignore, std::ignore) = hamGradHess(J, z);
+  return H;
+}
+
+template <class Scalar, int p>
+Vector<Scalar> getGradient(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+  Vector<Scalar> dH;
+  std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z);
+  return dH;
+}
+
+template <class Scalar, int p>
+Matrix<Scalar> getHessian(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+  Matrix<Scalar> ddH;
+  std::tie(std::ignore, std::ignore, ddH) = hamGradHess(J, z);
+  return ddH;
+}
+
+template <class Scalar>
+Real getThresholdEnergyDensity(unsigned p, Scalar κ, Scalar ε, Real a) {
+  Real apm2 = pow(a, p - 2);
+  Scalar δ = κ / apm2;
+  Real θ = arg(κ) + 2 * arg(ε);
+
+  return (p - 1) * apm2 / (2 * p) * pow(1 - norm(δ), 2) / (1 + norm(δ) - 2 * abs(δ) * cos(θ));
+}
+
+template <class Scalar, int p>
+Real getProportionOfThreshold(Scalar κ, const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+  Real N = z.size();
+  Scalar ε = getHamiltonian(J, z) / N;
+  Real a = z.squaredNorm() / N;
+
+  return norm(ε) / getThresholdEnergyDensity(p, κ, ε, a);
+}
+
 template <class Scalar>
 Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
   return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate();