#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 writeHeaderTo(csvFile &os) const {
      if (!xRanges.empty()) {
        for (const xRange &range : xRanges) {
          os << range.label + " max";
          os << range.label + " min";
        }
      }
      os << "Function Calls";
      os << "Solution Accuracy";
      //        os << std::endl;
    }

    void writeSettingsTo(csvFile &os) const {
      if (!xRanges.empty()) {
        for (const xRange &range : xRanges) {
          os << range.max;
          os << range.min;
        }
      }
      os << numberOfFunctionCalls;
      os << solutionAccuracy;
    }
  };

  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 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();
      if (!std::filesystem::is_directory(simulationScenarioPath)) {
        continue;
      }
      // 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));
        }
      }
    }
  }

  void saveMeshFiles() const {
    const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
    assert(numberOfSimulationJobs != 0 &&
           fullPatternSimulationJobs.size() ==
               reducedPatternSimulationJobs.size());

    fullPatternSimulationJobs[0]->pMesh->savePly("FullPattern_undeformed.ply");
    reducedPatternSimulationJobs[0]->pMesh->savePly(
        "ReducedPattern_undeformed.ply");
    FormFinder simulator;
    for (int simulationJobIndex = 0;
         simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
      // Drawing of full pattern results
      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
          fullPatternSimulationJobs[simulationJobIndex];
      SimulationResults fullModelResults =
          simulator.executeSimulation(pFullPatternSimulationJob);
      fullModelResults.saveDeformedModel();

      // Drawing of reduced pattern results
      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
          reducedPatternSimulationJobs[simulationJobIndex];
      SimulationResults reducedModelResults =
          simulator.executeSimulation(pReducedPatternSimulationJob);
      reducedModelResults.saveDeformedModel();
    }
  }

  void
  writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization,
                csvFile &os) {
    os << "Obj value";
    for (const ReducedModelOptimizer::xRange &range :
         settings_optimization.xRanges) {
      os << range.label;
    }
    os << "Time(s)";
    os << "#Crashes";
  }

  void
  writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization,
                 csvFile &os) const {
    os << objectiveValue;
    for (const double &optimalX : x) {
      os << optimalX;
    }

    for (int unusedXVarCounter = 0;
         unusedXVarCounter < settings_optimization.xRanges.size() - x.size();
         unusedXVarCounter++) {
      os << "-";
    }

    os << time;
    if (numberOfSimulationCrashes == 0) {
      os << "-";
    } else {
      os << numberOfSimulationCrashes;
    }
  }
};

#endif // REDUCEDMODELOPTIMIZER_HPP