#include "stokes.hpp"
#include <iostream>

using Complex = std::complex<Real>;

template<int ...ps, class Generator, typename... T>
void collectStokesData(unsigned N, Generator& r, double ε, T... μs) {
  pSpinModel<Real, ps...> ReM(N, r, μs...);

  std::normal_distribution<Real> Red(0, 1);

  Vector<Real> xMin = randomMinimum(ReM, Red, r, ε);

  Real Hx;
  Vector<Real> dHx;
  Matrix<Real> ddHx;
  std::tie(Hx, dHx, ddHx, std::ignore) = ReM.hamGradHess(xMin);
  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHx - xMin.dot(dHx) * Matrix<Real>::Identity(xMin.size(), xMin.size()) / Real(N));

  complex_normal_distribution<Real> d(0, 1, 0);

  pSpinModel M = ReM.template cast<Complex>();;
  Vector<Complex> zMin = xMin.cast<Complex>();

  Vector<Complex> zSaddle = zMin;

  while ((zSaddle - zMin).norm() < 1e-3 * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) {
    Vector<Complex> z0 = normalize(zSaddle + 0.5 * randomVector<Complex>(N, d, r));
    zSaddle = findSaddle(M, z0, ε);
  }

  Complex H1 = M.getHamiltonian(zMin);
  Complex H2 = M.getHamiltonian(zSaddle);

  Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());

  M *= exp(Complex(0, φ));

  Cord c(M, zMin, zSaddle, 3);
  c.relaxNewton(10, 1, 1e4);

  ((std::cout << ps), ...) << std::endl;;
  ((std::cout << μs), ...) << std::endl;;

  std::apply([](const Tensor<Real, ps>&... Js) -> void {
        ((std::cout << Js << std::endl), ...);
      } , ReM.Js);

  std::cout << xMin.transpose() << std::endl;
  std::cout << Hx << std::endl;
  std::cout << eigenS.eigenvalues().real().transpose() << std::endl;
}