ReducedModelOptimization/src/main.cpp

#include "csvfile.hpp"
#include "drmsimulationmodel.hpp"
#include "edgemesh.hpp"
#include "reducedmodeloptimizer.hpp"
#include "simulationhistoryplotter.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>
#include <filesystem>
#include <iostream>
#include <iterator>
#include <stdexcept>
#include <string>
#include <vcg/complex/algorithms/update/position.h>

#ifdef POLYSCOPE_DEFINED
#include "polyscope/curve_network.h"
#include "polyscope/point_cloud.h"
#include "polyscope/polyscope.h"
#endif
int main(int argc, char *argv[]) {
  if (argc < 3) {
    std::cerr << "Specify at least the two pattern filepaths to be "
                 "optimized.Exiting.."
              << std::endl;
    std::terminate();
  }

  // Populate the pattern pair to be optimized
  const int interfaceNodeIndex=3;
  ////Full pattern
  const std::string filepath_fullPattern = argv[1];
  PatternGeometry fullPattern(filepath_fullPattern);
  //  fullPattern.prependToLabel(std::to_string(fullPattern.EN()) + "#");
  fullPattern.scale(0.03, interfaceNodeIndex);
  ////Reduced pattern
  const std::string filepath_reducedPattern = argv[2];
  PatternGeometry reducedPattern(filepath_reducedPattern);
  reducedPattern.scale(0.03, interfaceNodeIndex);

  // Set the optization settings
  ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
  ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
  ReducedPatternOptimization::xRange beamI{"I", 0.001, 1000};
  ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
  ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
  ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
  ReducedPatternOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95};
  ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0};
  ReducedPatternOptimization::Settings settings_optimization;
  settings_optimization.xRanges
      = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
  const bool input_numberOfFunctionCallsDefined = argc >= 4;
  settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined
                                                    ? std::atoi(argv[3])
                                                    : 100;
  settings_optimization.normalizationStrategy
      = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;

  settings_optimization.translationNormalizationParameter = 1e-3;
  settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
  //  settings_optimization.translationNormalizationParameter = 1e-15;
  //  settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(1e-15);
  //  settings_optimization.solverAccuracy = 1e-3;
  settings_optimization.solverAccuracy = 1e-3;
  settings_optimization.objectiveWeights.translational = std::atof(argv[4]);
  settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]);

  std::string xConcatNames;
  for (const auto &x : settings_optimization.xRanges) {
      xConcatNames.append(x.label + "_");
  }
  xConcatNames.pop_back();

  // Optimize pairthere
  const std::string pairName = fullPattern.getLabel() + "@" + reducedPattern.getLabel();
  const std::string optimizationName = pairName + "("
                                       + std::to_string(settings_optimization.numberOfFunctionCalls)
                                       + "_"
                                       + to_string_with_precision(
                                           settings_optimization.objectiveWeights.translational)
                                       + ")" + "_" + xConcatNames;
  const bool input_resultDirectoryDefined = argc >= 6;
  const std::string optimizationResultsDirectory = input_resultDirectoryDefined
                                                       ? argv[5]
                                                       : std::filesystem::current_path().append(
                                                           "OptimizationResults").string();
  std::string resultsOutputDir;
  bool optimizationResultFolderExists = false;
  const std::filesystem::path crashedJobsDirPath(std::filesystem::path(optimizationResultsDirectory)
                                                     .append("CrashedJobs")
                                                     .append(optimizationName));
  if (std::filesystem::exists(crashedJobsDirPath)) {
      resultsOutputDir = crashedJobsDirPath.string();
      optimizationResultFolderExists = true;
  }
  const std::filesystem::path convergedJobsDirPath(
      std::filesystem::path(optimizationResultsDirectory)
          .append("ConvergedJobs")
          .append(optimizationName));
  if (std::filesystem::exists(convergedJobsDirPath)) {
      resultsOutputDir = convergedJobsDirPath.string();
      optimizationResultFolderExists = true;
  }

  ReducedPatternOptimization::Results optimizationResults;
  constexpr bool shouldReoptimize = true;
  bool optimizationAlreadyComputed = false;
  if (!shouldReoptimize && optimizationResultFolderExists) {
      const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir);
      if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) {
          optimizationAlreadyComputed = true;
      }
  }

  if (!optimizationAlreadyComputed) {
      auto start = std::chrono::system_clock::now();
      const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
      assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]);
      ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
      optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges);
      optimizer.optimize(settings_optimization, optimizationResults);
      optimizationResults.label = optimizationName;
      optimizationResults.baseTriangleFullPattern.copy(fullPattern);
      optimizationResults.settings = settings_optimization;
      auto end = std::chrono::system_clock::now();
      auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
      optimizationResults.time = elapsed.count() / 1000.0;

      if (optimizationResults.wasSuccessful) {
          resultsOutputDir = convergedJobsDirPath.string();
          csvFile csv_results({}, false);
          //          csvFile csv_results(std::filesystem::path(resultsOutputDir).append("optimizationDistances.csv"), false);
          csv_results << "Name";
          optimizationResults.writeHeaderTo(csv_results);
          settings_optimization.writeHeaderTo(csv_results);
          csv_results << endrow;
          csv_results << std::to_string(fullPattern.EN()) + "#" + pairName;
          optimizationResults.writeResultsTo(csv_results);
          settings_optimization.writeSettingsTo(csv_results);
          csv_results << endrow;
      } else {
          resultsOutputDir = crashedJobsDirPath.string();
      }
      optimizationResults.save(resultsOutputDir, true);
  }

#ifdef POLYSCOPE_DEFINED
  //    optimizationResults.saveMeshFiles();
  std::cout << "E:"
            << optimizationResults.reducedPatternSimulationJobs[0]
                   ->pMesh->elements[0]
                   .material.youngsModulus
            << std::endl;
  optimizationResults.draw();
#endif

  return 0;
}