ReducedModelOptimization/src/reducedmodeloptimizer.hpp

#ifndef REDUCEDMODELOPTIMIZER_HPP
#define REDUCEDMODELOPTIMIZER_HPP

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

#ifdef POLYSCOPE_DEFINED
#include "polyscope/color_management.h"
#endif // POLYSCOPE_DEFINED
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;
#ifdef POLYSCOPE_DEFINED
  struct StaticColors {
    glm::vec3 fullInitial;
    glm::vec3 fullDeformed;
    glm::vec3 reducedInitial;
    glm::vec3 reducedDeformed;
    StaticColors() {
      fullInitial = {0.416666, 0.109804, 0.890196};
      fullDeformed = {0.583333, 0.890196, 0.109804};
      reducedInitial = {0.890196, 0.109804, 0.193138};
      reducedDeformed = {0.109804, 0.890196, 0.806863};
    }
  };
  inline static StaticColors colors;
#endif // POLYSCOPE_DEFINED

public:
  inline static int fanSize{6};
  inline static VectorType patternPlaneNormal{0, 0, 1};
  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 {
    enum NormalizationStrategy {
      NonNormalized,
      Epsilon,
      MaxDisplacement,
      EqualDisplacements
    };
    inline static vector<std::string> normalizationStrategyStrings{
        "NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"};
    std::vector<xRange> xRanges;
    int numberOfFunctionCalls{100};
    double solutionAccuracy{1e-2};
    NormalizationStrategy normalizationStrategy{Epsilon};
    double normalizationParameter{0.003};

    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=" + normalizationStrategyStrings[normalizationStrategy];

      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 << "Normalization strategy";
      //        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;
      os << normalizationStrategyStrings[normalizationStrategy] + "_" +
                std::to_string(normalizationParameter);
    }
  };
  struct Results {
    double time{-1};
    int numberOfSimulationCrashes{0};
    std::vector<double> x;
    double objectiveValue;
    double rawObjectiveValue;
    std::vector<double> objectiveValuePerSimulationScenario;
    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));
          }
        }
      }
    }
#if POLYSCOPE_DEFINED
    void draw() const {
      initPolyscope();
      FormFinder simulator;
      LinearSimulationModel linearSimulator;
      assert(fullPatternSimulationJobs.size() ==
             reducedPatternSimulationJobs.size());
      fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
          colors.fullInitial);
      reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(
          colors.reducedInitial, false);

      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(colors.fullDeformed);
        SimulationResults fullModelLinearResults =
            linearSimulator.executeSimulation(pFullPatternSimulationJob);
        fullModelLinearResults.setLabelPrefix("linear");
        fullModelLinearResults.registerForDrawing(colors.fullDeformed, false);
        // Drawing of reduced pattern results
        const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
            reducedPatternSimulationJobs[simulationJobIndex];
        SimulationResults reducedModelResults =
            simulator.executeSimulation(pReducedPatternSimulationJob);
        reducedModelResults.registerForDrawing(colors.reducedDeformed);
        SimulationResults reducedModelLinearResults =
            linearSimulator.executeSimulation(pReducedPatternSimulationJob);
        reducedModelLinearResults.setLabelPrefix("linear");
        reducedModelLinearResults.registerForDrawing(colors.reducedDeformed,
                                                     false);
        polyscope::options::programName =
            fullPatternSimulationJobs[0]->pMesh->getLabel();
        polyscope::view::resetCameraToHomeView();
        polyscope::show();
        // Save a screensh
        const std::string screenshotFilename =
            "/home/iason/Coding/Projects/Approximating shapes with flat "
            "patterns/RodModelOptimizationForPatterns/Results/Images/" +
            fullPatternSimulationJobs[0]->pMesh->getLabel() + "_" +
            pFullPatternSimulationJob->getLabel();
        polyscope::screenshot(screenshotFilename, false);
        fullModelResults.unregister();
        reducedModelResults.unregister();
        reducedModelLinearResults.unregister();
        fullModelLinearResults.unregister();
        //                    double error = computeError(
        //                        reducedModelResults.displacements,fullModelResults.displacements,
        //                        );
        //                    std::cout << "Error of simulation scenario "
        //                              <<
        //                              simula         simulationScenarioStrings[simulationScenarioIndex]
        //                              << " is "
        //                              << error << std::endl;
      }
    }
#endif // POLYSCOPE_DEFINED
    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 << "Total raw Obj value";
      os << "Total Obj value";
      for (int simulationScenarioIndex = 0;
           simulationScenarioIndex <
           SimulationScenario::NumberOfSimulationScenarios;
           simulationScenarioIndex++) {
        os << "Obj value " + simulationScenarioStrings[simulationScenarioIndex];
      }
      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 << rawObjectiveValue;
      os << objectiveValue;
      for (double scenarioObjValue : objectiveValuePerSimulationScenario) {
        os << scenarioObjValue;
      }
      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;
      }
    }
  };

  inline static const std::string simulationScenarioStrings[] = {
      "Axial", "Shear", "Bending", "Dome", "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(
      PatternGeometry &fullPattern, PatternGeometry &reducedPatterm,
      const std::unordered_set<size_t> &reducedModelExcludedEges);

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

  static double objective(double x0, double x1, double x2, double x3,
                          double innerHexagonRotationAngle);
  static double objective(double b, double r, 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(
      PatternGeometry &fullModel, PatternGeometry &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,
      PatternGeometry &fullPattern, PatternGeometry &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 void registerResultsForDrawing(
      const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
      const std::unordered_map<ReducedPatternVertexIndex,
                               FullPatternVertexIndex>
          &reducedToFullInterfaceViMap);

  static double
  computeRawError(const std::vector<Vector6d> &reducedPatternDisplacements,
                  const std::vector<Vector6d> &fullPatternDisplacements,
                  const std::unordered_map<ReducedPatternVertexIndex,
                                           FullPatternVertexIndex>
                      &reducedToFullInterfaceViMap);

  static double
  computeError(const std::vector<Vector6d> &reducedPatternDisplacements,
               const std::vector<Vector6d> &fullPatternDisplacements,
               const std::unordered_map<ReducedPatternVertexIndex,
                                        FullPatternVertexIndex>
                   &reducedToFullInterfaceViMap,
               const double &normalizationFactor);

private:
    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(PatternGeometry &fullModel, PatternGeometry &reducedPattern,
                   const std::unordered_set<size_t> &reducedModelExcludedEges);
  void createSimulationMeshes(PatternGeometry &fullModel,
                              PatternGeometry &reducedModel);
  static void
  initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh);

  static double objective(long n, const double *x);
  FormFinder simulator;
  void computeObjectiveValueNormalizationFactors();
};
#endif // REDUCEDMODELOPTIMIZER_HPP