#include <getopt.h>
#include <limits>
#include <stereographic.hpp>
#include <unordered_map>
#include <list>

#include "Eigen/src/Eigenvalues/ComplexEigenSolver.h"
#include "Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h"
#include "complex_normal.hpp"
#include "p-spin.hpp"
#include "dynamics.hpp"
#include "stokes.hpp"

#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
#include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h"

#define PSPIN_P 3
const int p = PSPIN_P; // polynomial degree of Hamiltonian
using Complex = std::complex<Real>;
using ComplexVector = Vector<Complex>;
using ComplexMatrix = Matrix<Complex>;
using ComplexTensor = Tensor<Complex, p>;

using Rng = randutils::random_generator<pcg32>;

std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, Real δ) {
  std::list<std::array<ComplexVector, 2>> pairs;

  for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) {
    for (auto it2 = std::next(it1); it2 != saddles.end(); it2++) {
      if ((it1->second - it2->second).norm() < δ) {
        pairs.push_back({it1->second, it2->second});
      }
    }
  }

  return pairs;
}

template <class Generator>
std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, Real ε, Real Δ, Generator& r) {
  Real σ = sqrt(factorial(p) / (Real(2) * pow(z10.size(), p - 1)));
  complex_normal_distribution<Real> dJ(0, σ, 0);

  ComplexTensor J = J0;
  ComplexVector z1, z2;
  Real ΔH = abs(imag(getHamiltonian(J, z10) - getHamiltonian(J, z20))) / z10.size();
  Real γ = ΔH;

  std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
    [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
      Complex Ji = getJ<Complex, p>(JJ, ind);
      setJ<Complex, p>(JJ, ind, Ji + γ * dJ(r.engine()));
    };

  while (ΔH > ε) {
    ComplexTensor Jp = J;
    iterateOver<Complex, p>(Jp, perturbJ);

    try {
      z1 = findSaddle(Jp, z10, Δ);
      z2 = findSaddle(Jp, z20, Δ);

      Real Δz = (z1 - z2).norm();

      if (Δz > 1e-2) {
        Real ΔHNew = abs(imag(getHamiltonian(Jp, z1) - getHamiltonian(Jp, z2))) / z1.size();

        if (ΔHNew < ΔH) {
          J = Jp;
          ΔH = ΔHNew;
          γ = ΔH;

          std::cerr << "Match imaginary energies: Found couplings with ΔH = " << ΔH << std::endl;
        }
      }
    } catch (std::exception& e) {}
  }

  return {J, z1, z2};
}

int main(int argc, char* argv[]) {
  // model parameters
  unsigned N = 10; // number of spins
  Real T = 1;    // temperature
  Real Rκ = 0;   // real part of distribution parameter
  Real Iκ = 0;   // imaginary part of distribution parameter
  // simulation parameters
  Real ε = 1e-15;
  Real εJ = 1e-5;
  Real δ = 1;   // threshold for determining saddle
  Real Δ = 1e-3;
  Real γ = 1e-1;   // step size
  unsigned t = 1000; // number of Langevin steps
  unsigned M = 100;
  unsigned n = 100;

  int opt;

  while ((opt = getopt(argc, argv, "N:M:n:T:e:r:i:g:t:E:")) != -1) {
    switch (opt) {
    case 'N':
      N = (unsigned)atof(optarg);
      break;
    case 't':
      t = (unsigned)atof(optarg);
      break;
    case 'T':
      T = atof(optarg);
      break;
    case 'e':
      δ = atof(optarg);
      break;
    case 'E':
      ε = atof(optarg);
      break;
    case 'g':
      γ = atof(optarg);
    case 'r':
      Rκ = atof(optarg);
      break;
    case 'i':
      Iκ = atof(optarg);
      break;
    case 'n':
      n = (unsigned)atof(optarg);
      break;
    case 'M':
      M = (unsigned)atof(optarg);
      break;
    default:
      exit(1);
    }
  }

  Complex κ(Rκ, Iκ);
  Real σ = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));

  Rng r;

  Tensor<Real, p> ReJ = generateCouplings<Real, p>(N, std::normal_distribution<Real>(0, σ), r.engine());

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

  Vector<Real> zMin = randomMinimum(ReJ, Red, r, ε);
  auto [Hr, dHr, ddHr] = hamGradHess(ReJ, zMin);
  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - ((ddHr * zMin) * zMin.transpose()) / (Real)zMin.size());
  std::cout << eigenS.eigenvalues().transpose() << std::endl;
  getchar();

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

  ComplexTensor J = generateCouplings<Complex, p>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine());

  std::function<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
    (const ComplexTensor& J, const ComplexVector& z) {
      Real W;
      std::tie(W, std::ignore) = WdW(J, z);
      return W;
    };

  ComplexVector zSaddle = randomSaddle(J, d, r, ε);
  std::cerr << "Found saddle." << std::endl;

  ComplexVector zSaddleNext;
  bool foundSaddle = false;
  while (!foundSaddle) {
    ComplexVector z0 = normalize(zSaddle + δ * randomVector<Complex>(N, d, r.engine()));
    try {
      zSaddleNext = findSaddle(J, z0, ε);
      Real saddleDistance = (zSaddleNext - zSaddle).norm();
      if (saddleDistance / N > 1e-2) {
        std::cout << saddleDistance << std::endl;
        foundSaddle = true;
      }
    } catch (std::exception& e) {}
  }

  auto [H1, dH1, ddH1] = hamGradHess(J, zSaddle);
  auto [H2, dH2, ddH2] = hamGradHess(J, zSaddleNext);

  Eigen::ComplexEigenSolver<ComplexMatrix> ces;
  ces.compute(ddH1);

  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;
  Real smallestNorm = std::numeric_limits<Real>::infinity();
  for (Complex z : ces.eigenvalues()) {
    if (norm(z) < smallestNorm) {
      smallestNorm = norm(z);
    }
  }
  std::cerr << smallestNorm << std::endl;

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

  /*
  if (stokesLineTest<Real>(J, zSaddle, zSaddleNext, 10, 4)) {
    std::cerr << "Found a Stokes line" << std::endl;
//    stokesLineTestNew<Real>(J, zSaddle, zSaddleNext, 10, 3);
  } else {
    std::cerr << "Didn't find a Stokes line" << std::endl;
  }
  */

  /*
  J(0,0,0) = Complex(2, 3);
  J(1,1,1) = Complex(-2, 0.3);
  J(0,1,1) = Complex(4, 0.2);
  J(1,0,1) = Complex(4, 0.2);
  J(1,1,0) = Complex(4, 0.2);
  J(1,0,0) = Complex(0.7, 0.4);
  J(0,1,0) = Complex(0.7, 0.4);
  J(0,0,1) = Complex(0.7, 0.4);

  ComplexVector z0(2);;
  z0 << Complex(0.8, 0.3), Complex(0.7, 0.2);
  ComplexVector z1(2);
  z1 << Complex(-0.5, 0.2), Complex(1.0, 1.0);

  Cord test(J, z0, z1, 2);

  test.gs[0](0) = Complex(0.2, 0.2);
  test.gs[0](1) = Complex(0.4, 0.4);
  test.gs[1](0) = Complex(0.1, 0.2);
  test.gs[1](1) = Complex(0.3, 0.4);

  auto [dgs, ddgs] = test.gsGradHess(J, 0.7);

  std::cout << dgs << std::endl;
  std::cout << ddgs << std::endl;
  */

  Cord test(J, zSaddle, zSaddleNext, 5);
  test.relaxNewton(J, 25, 1, 1e4);

  std::cout << test.z0.transpose() << std::endl;
  std::cout << test.z1.transpose() << std::endl;

  for (Vector<Complex>& g : test.gs) {
    std::cout << g.transpose() << std::endl;
  }

  return 0;
}