#ifndef REDUCEDMODELOPTIMIZER_HPP
#define REDUCEDMODELOPTIMIZER_HPP

#include "beamformfinder.hpp"
#include "csvfile.hpp"
#include "edgemesh.hpp"
#include "elementalmesh.hpp"
#include "matplot/matplot.h"
#include <Eigen/Dense>

using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;

class ReducedModelOptimizer {
  std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
  std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
  std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
      m_fullToReducedInterfaceViMap;
  std::unordered_map<FullPatternVertexIndex, FullPatternVertexIndex>
      m_fullPatternOppositeInterfaceViMap;
  std::unordered_map<size_t, size_t> nodeToSlot;
  std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
  std::vector<double> initialGuess;

public:
  enum SimulationScenario {
    Axial,
    Shear,
    Bending,
    Dome,
    Saddle,
    NumberOfSimulationScenarios
  };

  struct xRange {
    std::string label;
    double min;
    double max;
    std::string toString() const {
      return label + "=[" + std::to_string(min) + "," + std::to_string(max) +
             "]";
    }
  };
  struct Results;

  struct Settings {
    std::vector<xRange> xRanges;
    int numberOfFunctionCalls{100};
    double solutionAccuracy{1e-2};
    bool normalizeObjectiveValue{ true };

    std::string toString() const {
      std::string settingsString;
      if (!xRanges.empty()) {
        std::string xRangesString;
        for (const xRange &range : xRanges) {
          xRangesString += range.toString() + " ";
        }
        settingsString += xRangesString;
      }
      settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) +
                        " Accuracy=" + std::to_string(solutionAccuracy)+" Norm="+(normalizeObjectiveValue ? "yes" : "no");

      return settingsString;
    }

    void writeTo(csvFile &csv) const {
      // Create settings csv header
      if (!xRanges.empty()) {
        for (const xRange &range : xRanges) {
          csv << range.label + " max";
          csv << range.label + " min";
        }
      }
      csv << "Function Calls";
      csv << "Solution Accuracy";
      csv << "Normalize obj value";
      csv << endrow;

      if (!xRanges.empty()) {
        for (const xRange &range : xRanges) {
          csv << range.max;
          csv << range.min;
        }
      }
      csv << numberOfFunctionCalls;
      csv << solutionAccuracy;
      csv <<(normalizeObjectiveValue ? "yes" : "no");
      csv << endrow;
    }
  };

  inline static const std::string simulationScenarioStrings[] = {
      "Axial", "Shear", "Bending", "Double", "Saddle"};
  Results optimize(const Settings &xRanges,
                   const std::vector<SimulationScenario> &simulationScenarios =
                       std::vector<SimulationScenario>());
  double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;

  ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
  static void computeReducedModelSimulationJob(
      const SimulationJob &simulationJobOfFullModel,
      const std::unordered_map<size_t, size_t> &simulationJobFullToReducedMap,
      SimulationJob &simulationJobOfReducedModel);

  SimulationJob
  getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);

  void initializePatterns(
      FlatPattern &fullPattern, FlatPattern &reducedPatterm,
      const std::unordered_set<size_t> &reducedModelExcludedEges);

  void setInitialGuess(std::vector<double> v);

  static void runBeamOptimization();

  static void runSimulation(const std::string &filename,
                            std::vector<double> &x);

  static double objective(double x0, double x1, double x2, double x3);
  static double objective(double b, double h, double E);

  static std::vector<std::shared_ptr<SimulationJob>>
  createScenarios(const std::shared_ptr<SimulationMesh> &pMesh,
                  const std::unordered_map<size_t, size_t>
                      &fullPatternOppositeInterfaceViMap);
  static void createSimulationMeshes(
      FlatPattern &fullModel, FlatPattern &reducedModel,
      std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
      std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
  static void computeMaps(
      const std::unordered_set<size_t> &reducedModelExcludedEdges,
      const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
      FlatPattern &fullPattern, FlatPattern &reducedPattern,
      std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
          &reducedToFullInterfaceViMap,
      std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
          &fullToReducedInterfaceViMap,
      std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
          &fullPatternOppositeInterfaceViMap);
  static void
  visualizeResults(const std::vector<std::shared_ptr<SimulationJob>>
                       &fullPatternSimulationJobs,
                   const std::vector<std::shared_ptr<SimulationJob>>
                       &reducedPatternSimulationJobs,
                   const std::vector<SimulationScenario> &simulationScenarios,
                   const std::unordered_map<ReducedPatternVertexIndex,
                                            FullPatternVertexIndex>
                       &reducedToFullInterfaceViMap);

  static double computeError(const std::vector<Vector6d> &reducedPatternResults,
                             const std::vector<Vector6d> &fullPatternResults,const double& interfaceDisplacementNormSum,
                             const std::unordered_map<ReducedPatternVertexIndex,
                                                      FullPatternVertexIndex>
                                 &reducedToFullInterfaceViMap);

private:
  void
  visualizeResults(const std::vector<std::shared_ptr<SimulationJob>>
                       &fullPatternSimulationJobs,
                   const std::vector<SimulationScenario> &simulationScenarios);
  static void computeDesiredReducedModelDisplacements(
      const SimulationResults &fullModelResults,
      const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
      Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
  static Results runOptimization(const Settings &settings);
  std::vector<std::shared_ptr<SimulationJob>>
  createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
  void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern,
                   const std::unordered_set<size_t> &reducedModelExcludedEges);
  void createSimulationMeshes(FlatPattern &fullModel,
                              FlatPattern &reducedModel);
  static void
  initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh);

  static double
  computeError(const SimulationResults &reducedPatternResults,
               const Eigen::MatrixX3d &optimalReducedPatternDisplacements);
  static double objective(long n, const double *x);
  FormFinder simulator;
};
struct ReducedModelOptimizer::Results {
  double time{-1};
  int numberOfSimulationCrashes{0};
  std::vector<double> x;
  double objectiveValue;
  std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
  std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;

  void draw() const {
    initPolyscope();
    FormFinder simulator;
    assert(fullPatternSimulationJobs.size() ==
           reducedPatternSimulationJobs.size());
    fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
    reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();

    const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
    for (int simulationJobIndex = 0;
         simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
      // Drawing of full pattern results
      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
          fullPatternSimulationJobs[simulationJobIndex];
      pFullPatternSimulationJob->registerForDrawing(
          fullPatternSimulationJobs[0]->pMesh->getLabel());
      SimulationResults fullModelResults =
          simulator.executeSimulation(pFullPatternSimulationJob);
      fullModelResults.registerForDrawing();
      // Drawing of reduced pattern results
      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
          reducedPatternSimulationJobs[simulationJobIndex];
      SimulationResults reducedModelResults =
          simulator.executeSimulation(pReducedPatternSimulationJob);
      reducedModelResults.registerForDrawing();
      polyscope::show();
      // Save a screensh
      //    const std::string screenshotFilename =
      //        "/home/iason/Coding/Projects/Approximating shapes with flat "
      //        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
      //        + m_pFullPatternSimulationMesh->getLabel() + "_" +
      //        simulationScenarioStrings[simulationScenarioIndex];
      //    polyscope::screenshot(screenshotFilename, false);
      fullModelResults.unregister();
      reducedModelResults.unregister();
      //    double error = computeError(
      //        reducedModelResults,
      //        global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
      //    std::cout << "Error of simulation scenario "
      //              << simulationScenarioStrings[simulationScenarioIndex] << "
      //              is "
      //              << error << std::endl;
    }
  }
  void save(const string &saveToPath) const {
    assert(std::filesystem::is_directory(saveToPath));

    const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
    for (int simulationJobIndex = 0;
         simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
          fullPatternSimulationJobs[simulationJobIndex];
      std::filesystem::path simulationJobFolderPath(
          std::filesystem::path(saveToPath)
              .append(pFullPatternSimulationJob->label));
      std::filesystem::create_directory(simulationJobFolderPath);
      const auto fullPatternDirectoryPath =
          std::filesystem::path(simulationJobFolderPath).append("Full");
      std::filesystem::create_directory(fullPatternDirectoryPath);
      pFullPatternSimulationJob->save(fullPatternDirectoryPath.string());
      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
          reducedPatternSimulationJobs[simulationJobIndex];
      const auto reducedPatternDirectoryPath =
          std::filesystem::path(simulationJobFolderPath)
              .append("Reduced");
      if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
          std::filesystem::create_directory(reducedPatternDirectoryPath);
      }
      pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string());
    }
  }
  void load(const string &loadFromPath) {
    assert(std::filesystem::is_directory(loadFromPath));

    for (const auto &directoryEntry :
         filesystem::directory_iterator(loadFromPath)) {
      const auto simulationScenarioPath = directoryEntry.path();
      // Load reduced pattern files
      for (const auto &fileEntry : filesystem::directory_iterator(
               std::filesystem::path(simulationScenarioPath)
                   .append("Full"))) {
        const auto filepath = fileEntry.path();
        if (filepath.extension() == ".json") {
          SimulationJob job;
          job.load(filepath.string());
          fullPatternSimulationJobs.push_back(
              std::make_shared<SimulationJob>(job));
        }
      }

      // Load full pattern files
      for (const auto &fileEntry : filesystem::directory_iterator(
               std::filesystem::path(simulationScenarioPath)
                   .append("Reduced"))) {
        const auto filepath = fileEntry.path();
        if (filepath.extension() == ".json") {
          SimulationJob job;
          job.load(filepath.string());
          reducedPatternSimulationJobs.push_back(
              std::make_shared<SimulationJob>(job));
        }
      }
    }
  }
};

#endif // REDUCEDMODELOPTIMIZER_HPP