#include <iostream>
#include <Kokkos_Core.hpp>
#include <RevisionInfo.hpp>
#include <rang.hpp>
#include <FPEManager.hpp>
#include <SignalManager.hpp>
#include <ConsoleManager.hpp>

// #include <RawKokkosAcousticSolver.hpp>
// #include <MeshLessAcousticSolver.hpp>
// #include <AcousticSolverClass.hpp>
// #include <AcousticSolverTest.hpp>

#include <Connectivity1D.hpp>
#include <Mesh.hpp>
#include <BoundaryCondition.hpp>
#include <AcousticSolver.hpp>

#include <TinyVector.hpp>
#include <TinyMatrix.hpp>

#include <GmshReader.hpp>

#include <CLI/CLI.hpp>
#include <cassert>
#include <limits>
#include <map>

int main(int argc, char *argv[])
{
  long unsigned number = 10;
  std::string filename;
  {
    CLI::App app{"Pastis help"};

    app.add_option("number,-n,--number", number, "Number of cells");//->required();

    app.add_option("filename,-f,--filename", filename, "gmsh file");//->required();

    int threads=-1;
    app.add_option("--threads", threads, "Number of Kokkos threads")->check(CLI::Range(1,std::numeric_limits<decltype(threads)>::max()));

    std::string colorize="auto";
    app.add_set("--colorize", colorize, {"auto", "yes", "no"}, "Colorize console output", true);

    bool disable_fpe = false;
    app.add_flag("--no-fpe", disable_fpe, "Do not trap floating point exceptions");
    bool disable_signals = false;
    app.add_flag("--no-signal", disable_signals, "Do not catches signals");

    std::string pause_on_error="auto";
    app.add_set("--pause-on-error", pause_on_error, {"auto", "yes", "no"}, "Pause for debugging on unexpected error", true);

    std::atexit([](){std::cout << rang::style::reset;});
    try {
      app.parse(argc, argv);
    } catch (const CLI::ParseError &e) {
      return app.exit(e);
    }

    ConsoleManager::init(colorize);
    FPEManager::init(not disable_fpe);
    SignalManager::setPauseForDebug(pause_on_error);
    SignalManager::init(not disable_signals);
  }
  
  std::cout << "Code version: "
	    << rang::style::bold << RevisionInfo::version() << rang::style::reset << '\n';

  std::cout << "-------------------- "
	    << rang::fg::green
	    << "git info"
	    << rang::fg::reset
	    <<" -------------------------"
	    << '\n';
  std::cout << "tag:  " << rang::fg::reset
	    << rang::style::bold << RevisionInfo::gitTag() << rang::style::reset << '\n';
  std::cout << "HEAD: " << rang::style::bold << RevisionInfo::gitHead() << rang::style::reset << '\n';
  std::cout << "hash: " << rang::style::bold << RevisionInfo::gitHash() << rang::style::reset << "  (";

  if (RevisionInfo::gitIsClean()) {
    std::cout << rang::fgB::green << "clean" << rang::fg::reset;
  } else {
    std::cout << rang::fgB::red << "dirty" << rang::fg::reset; 
  }
  std::cout << ")\n";
  std::cout << "-------------------------------------------------------\n";

  Kokkos::initialize(argc, argv);
  Kokkos::DefaultExecutionSpace::print_configuration(std::cout);

  std::map<std::string, double> method_cost_map;

  // { // Basic function based acoustic solver
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   RawKokkos::AcousticSolver(number);
  //   method_cost_map["RawKokkos"] = timer.seconds();
  // }

  // { // class for acoustic solver (mesh less)
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   MeshLessAcousticSolver acoustic_solver(number);
  //   method_cost_map["MeshLessAcousticSolver"] = timer.seconds();
  // }

  // { // class for acoustic solver
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   AcousticSolverClass acoustic_solver(number);
  //   method_cost_map["AcousticSolverClass"] = timer.seconds();
  // }

  // { // class for acoustic solver test
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   AcousticSolverTest acoustic_solver(number);
  //   method_cost_map["AcousticSolverTest"] = timer.seconds();
  // }

  if (filename != "") {
    std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
    GmshReader gmsh_reader(filename);

    typedef Mesh<Connectivity2D> MeshType;
    typedef MeshData<MeshType> MeshDataType;
    typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;

    MeshType& mesh = gmsh_reader.mesh();

    Kokkos::Timer timer;
    timer.reset();
    MeshDataType mesh_data(mesh);

    std::vector<BoundaryConditionHandler> bc_list;
    { // quite dirty!
      for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) {
	Connectivity2D::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
	std::cout << i_boundary
		  << " -> "
		  << mesh.connectivity().nodesBoundary(i_boundary).second.extent(0)
		  << '\n';
	unsigned int ref = nodes_boundary.first;
	TinyVector<2> normal(0,0);
	if ((ref == 3) or (ref == 4)) {
	  normal = TinyVector<2>(0,1);
	} else {
	  normal = TinyVector<2>(1,0);
	}
	std::cout << "boundary=" << i_boundary << " ref=" << ref << " n=" << normal << '\n';
	const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second;
	std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0));
	for (size_t r=0; r<nodes_ids.extent(0); ++r) {
	  node_boundary_vector[r] = nodes_ids[r];
	}
	SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
	  = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal);
	std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
	bc_list.push_back(BoundaryConditionHandler(bc));
      }

      // PressureBoundaryCondition* pres_bc1
      // 	= new PressureBoundaryCondition(0.1,
      // 					std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
      // std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
      // bc_list.push_back(BoundaryConditionHandler(bc1));
    }

    UnknownsType unknowns(mesh_data);

    unknowns.initializeSod();

    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);

    typedef TinyVector<MeshType::dimension> Rd;

    const Kokkos::View<const double*> Vj = mesh_data.Vj();
    const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();

    const double tmax=0.2;
    double t=0;

    int itermax=std::numeric_limits<int>::max();
    int iteration=0;

    Kokkos::View<double*> rhoj = unknowns.rhoj();
    Kokkos::View<double*> ej = unknowns.ej();
    Kokkos::View<double*> pj = unknowns.pj();
    Kokkos::View<double*> gammaj = unknowns.gammaj();
    Kokkos::View<double*> cj = unknowns.cj();

    BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);

    while((t<tmax) and (iteration<itermax)) {
      double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
      if (t+dt>tmax) {
	dt=tmax-t;
      }
      acoustic_solver.computeNextStep(t,dt, unknowns);

      block_eos.updatePandCFromRhoE();    
    
      t += dt;
      ++iteration;
    }

    std::ofstream gnuplot("sol.gnu");
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      for (int r=0; r<3; ++r) {
	const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,r)];
	gnuplot << x[0] << ' ' << x[1] << '\n';
      }
      const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,0)];
      gnuplot << x[0] << ' ' << x[1] << "\n\n";
    }
    gnuplot.close();

    std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
	      << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";

    method_cost_map["AcousticSolverWithMesh"] = timer.seconds();

  } else {
    // class for acoustic solver test
    Kokkos::Timer timer;
    timer.reset();
    Connectivity1D connectivity(number);
    typedef Mesh<Connectivity1D> MeshType;
    typedef MeshData<MeshType> MeshDataType;
    typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;

    MeshType mesh(connectivity);
    MeshDataType mesh_data(mesh);

    std::vector<BoundaryConditionHandler> bc_list;
    { // quite dirty!
      SymmetryBoundaryCondition<MeshType::dimension>* sym_bc0
	= new SymmetryBoundaryCondition<MeshType::dimension>(std::vector<unsigned int>({0u}),
							     TinyVector<1>(-1));
      std::shared_ptr<SymmetryBoundaryCondition<1>> bc0(sym_bc0);
      bc_list.push_back(BoundaryConditionHandler(bc0));

      PressureBoundaryCondition* pres_bc1
	= new PressureBoundaryCondition(0.1,
					std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
      std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
      bc_list.push_back(BoundaryConditionHandler(bc1));
    }

    UnknownsType unknowns(mesh_data);

    unknowns.initializeSod();

    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);

    typedef TinyVector<MeshType::dimension> Rd;

    const Kokkos::View<const double*> Vj = mesh_data.Vj();
    const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();

    const double tmax=0.2;
    double t=0;

    int itermax=std::numeric_limits<int>::max();
    int iteration=0;

    Kokkos::View<double*> rhoj = unknowns.rhoj();
    Kokkos::View<double*> ej = unknowns.ej();
    Kokkos::View<double*> pj = unknowns.pj();
    Kokkos::View<double*> gammaj = unknowns.gammaj();
    Kokkos::View<double*> cj = unknowns.cj();

    BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);

    while((t<tmax) and (iteration<itermax)) {
      double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
      if (t+dt>tmax) {
	dt=tmax-t;
      }

      acoustic_solver.computeNextStep(t,dt, unknowns);

      block_eos.updatePandCFromRhoE();    
    
      t += dt;
      ++iteration;
    }

    std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
	      << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";

    method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
    
    { // gnuplot output for density
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> rhoj = unknowns.rhoj();
      std::ofstream fout("rho");
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
        fout << xj[j][0] << ' ' << rhoj[j] << '\n';
      }
    }

  }

  Kokkos::finalize();

  std::string::size_type size=0;
  for (const auto& method_cost : method_cost_map) {
    size = std::max(size, method_cost.first.size());
  }

  for (const auto& method_cost : method_cost_map) {
    std::cout << "* ["
      	      << rang::fgB::cyan
	      << std::setw(size) << std::left
	      << method_cost.first
	      << rang::fg::reset
	      << "] Execution time: "
	      << rang::style::bold
	      << method_cost.second
	      << rang::style::reset << '\n';
  }
  
  return 0;
}