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

#include "Eigen/Dense"
#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 RealVector = Vector<Real>;
using RealMatrix = Matrix<Real>;
using RealTensor = Tensor<Real, p>;
using ComplexVector = Vector<Complex>;
using ComplexMatrix = Matrix<Complex>;
using ComplexTensor = Tensor<Complex, p>;

using Rng = randutils::random_generator<pcg32>;

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);
    }
  }

  Rng r;

//  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);

  RealVector zMin = randomMinimum(pSpin, Red, r, ε);
  Real Hr;
  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() + zMin * zMin.transpose() / (Real)N);
  for (unsigned i = 0; i < N; i++) {
    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 complexPSpin = pSpin.cast<Complex>();;
  ComplexVector zSaddle = zMin.cast<Complex>();

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

  Complex H1 = complexPSpin.getHamiltonian(zSaddle);
  Complex H2 = complexPSpin.getHamiltonian(zSaddleNext);

  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;

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

  Cord test(complexPSpin, zSaddle, zSaddleNext, 3);
  test.relaxNewton(10, 1, 1e4);

  std::cout << test.totalCost(10) << " " << test.totalCost(20) << std::endl;

  return 0;
}