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 <string_view>
#include <vcg/complex/algorithms/update/position.h>
#include <boost/algorithm/string.hpp>

#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 <= 7) {
        std::cerr << "Wrong number of input parameters. Expects at least 6 input parameters."
                     "Usage:\n"
                     "1)full pattern file path\n"
                     "2)reduced pattern file path\n"
                     "3)Number of optimizaion function calls\n"
                     "4)Translational error weight\n"
                     "5)Optimization results directory path\n"
                     "6)[optional]Optimization parameters\n"
                     "7)[optional]Intermediate results directory path\n"
                     "Exiting.."
                  << std::endl;
        return 1;
    }

    // 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 beamI2{"I2", 0.001, 1000};
    ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
    ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
    ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95};
    ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
    ReducedPatternOptimization::Settings settings_optimization;
    settings_optimization.parameterRanges
        = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
    settings_optimization.numberOfFunctionCalls = std::atoi(argv[3]);

    const bool input_optimizationParametersDefined = argc >= 8;
    if (input_optimizationParametersDefined) {
        const std::string optimizationParametersTag = argv[7];
        std::vector <std::string> split=Utilities::split(optimizationParametersTag,",");
        //parse parameter tag
        std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>
            optimizationParameters;
        std::vector<ReducedPatternOptimization::OptimizationParameterIndex> parameterGroup;
        for (std::string &s : split) {
//        std::cout<<s<<std::endl;
//        s=Utilities::trimLeftAndRightSpaces(s);
            if(boost::algorithm::contains(s,"{")){
//            std::cout<<"{ detected"<<std::endl;
                parameterGroup.clear();
                }
               if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::E)))){
                parameterGroup.push_back(ReducedPatternOptimization::E);
                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::A)))){
                parameterGroup.push_back(ReducedPatternOptimization::A);
                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::I2)))){
                parameterGroup.push_back(ReducedPatternOptimization::I2);
                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::I3)))){
                parameterGroup.push_back(ReducedPatternOptimization::I3);
                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::J)))){
                parameterGroup.push_back(ReducedPatternOptimization::J);
                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::R)))){
                parameterGroup.push_back(ReducedPatternOptimization::R);
                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::Theta)))){
                parameterGroup.push_back(ReducedPatternOptimization::Theta);
                }
                else{
                std::cerr << "Wrong optimization parameter input: " << optimizationParametersTag
                          << std::endl;

                }
                if(boost::algorithm::contains(s,"}")){
                optimizationParameters.push_back(parameterGroup);
                }

        }
        settings_optimization.optimizationVariables=optimizationParameters;
        #ifdef POLYSCOPE_DEFINED
        for(const std::vector<ReducedPatternOptimization::OptimizationParameterIndex>& v:settings_optimization.optimizationVariables){
        for(const ReducedPatternOptimization::OptimizationParameterIndex& i:v ){
        std::cout<<static_cast<int>(i)<<" ";
        }
        std::cout<<std::endl;
        }
        #endif

    } else {
        enum OptimizationParameterComparisonScenarioIndex {
            AllVar,
            GeoYM,
            noYM,
            YMMat_Geo,
            YM_MatGeo,
            MatGeo_YM,
            Geo_YM_Mat,
            YM_Geo_Mat,
            Geo_Mat,
            YMGeo_Mat,
            NumberOfScenarios
        };
        const std::vector<
            std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
            optimizationParameters = [&]() {
                std::vector<
                    std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
                    optimizationParameters(NumberOfScenarios);
                using namespace ReducedPatternOptimization;
                optimizationParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}};
                optimizationParameters[GeoYM] = {{R, Theta, E}};
                optimizationParameters[noYM] = {{A, I2, I3, J, R, Theta}};
                optimizationParameters[YMMat_Geo] = {{E, A, I2, I3, J}, {R, Theta}};
                optimizationParameters[YM_MatGeo] = {{E}, {A, I2, I3, J, R, Theta}};
                optimizationParameters[MatGeo_YM] = {{A, I2, I3, J, R, Theta}, {E}};
                optimizationParameters[Geo_YM_Mat] = {{R, Theta}, {E}, {A, I2, I3, J}};
                optimizationParameters[YM_Geo_Mat] = {{E}, {R, Theta}, {A, I2, I3, J}};
                optimizationParameters[Geo_Mat] = {{R, Theta}, {A, I2, I3, J}};
                optimizationParameters[YMGeo_Mat] = {{E, R, Theta}, {A, I2, I3, J}};
                return optimizationParameters;
            }();

        constexpr OptimizationParameterComparisonScenarioIndex scenario = YMGeo_Mat;
        settings_optimization.optimizationVariables = optimizationParameters[scenario];
    }
    if(settings_optimization.optimizationVariables.empty()){
    std::cerr<<"No optimization variables. Exiting.."<<std::endl;
    std::terminate();
    }

    const bool input_intermediateResultsDirectoryDefined = argc >= 7;
    if (input_intermediateResultsDirectoryDefined) {
        settings_optimization.intermediateResultsDirectoryPath = std::filesystem::path(argv[6]);
    }

    settings_optimization.normalizationStrategy
        = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
#ifdef POLYSCOPE_DEFINED
  settings_optimization.translationNormalizationParameter = 0;
  settings_optimization.rotationNormalizationParameter = 0;
#else
  settings_optimization.translationNormalizationParameter = 3e-4;
  settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
#endif
  settings_optimization.solverAccuracy = 1e-1;
  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.parameterRanges) {
      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 std::string optimizationResultsDirectory = argv[5];
  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.parameterRanges);
      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 = crashedJobsDirPath.string();
          return 1;
      }
      resultsOutputDir = convergedJobsDirPath.string();
      optimizationResults.save(resultsOutputDir, true);
      csvFile csv_resultsLocalFile(std::filesystem::path(resultsOutputDir).append("results.csv"),
                                   true);
      csvFile csv_results({}, false);
      std::vector<csvFile *> csvVector{&csv_resultsLocalFile, &csv_results};
      csv_results << "Name";
      csv_resultsLocalFile << "Name";
      optimizationResults.writeHeaderTo(csvVector);
      settings_optimization.writeHeaderTo(csv_results);
      csv_results << endrow;
      csv_resultsLocalFile << endrow;
      csv_results << std::to_string(fullPattern.EN()) + "#" + pairName;
      optimizationResults.writeResultsTo(csvVector);
      settings_optimization.writeSettingsTo(csv_results);
      csv_results << endrow;
      csv_resultsLocalFile << endrow;
  }

#ifdef POLYSCOPE_DEFINED
  std::vector<std::string> scenarioLabels(
      optimizationResults.objectiveValue.perSimulationScenario_total.size());
  const double colorAxial = 1;
  const double colorShear = 3;
  const double colorBending = 5;
  const double colorDome = 7;
  const double colorSaddle = 9;
  std::vector<double> colors(optimizationResults.objectiveValue.perSimulationScenario_total.size());
  for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
      scenarioLabels[scenarioIndex]
          = optimizationResults.reducedPatternSimulationJobs[scenarioIndex]->getLabel();
      if (scenarioLabels[scenarioIndex].rfind("Axial", 0) == 0) {
          colors[scenarioIndex] = colorAxial;

      } else if (scenarioLabels[scenarioIndex].rfind("Shear", 0) == 0) {
          colors[scenarioIndex] = colorShear;

      } else if (scenarioLabels[scenarioIndex].rfind("Bending", 0) == 0) {
          colors[scenarioIndex] = colorBending;

      } else if (scenarioLabels[scenarioIndex].rfind("Dome", 0) == 0) {
          colors[scenarioIndex] = colorDome;
      } else if (scenarioLabels[scenarioIndex].rfind("Saddle", 0) == 0) {
          colors[scenarioIndex] = colorSaddle;
      } else {
          std::cerr << "Label could not be identified" << std::endl;
      }
  }
  std::vector<double> y(optimizationResults.objectiveValue.perSimulationScenario_total.size());
  for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
      y[scenarioIndex]
          = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex]
            + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex];
      //      y[scenarioIndex] = optimizationResults.objectiveValue
      //                             .perSimulationScenario_total[scenarioIndex];
  }
  std::vector<double> x = matplot::linspace(0, y.size() - 1, y.size());
  std::vector<double> markerSizes(y.size(), 5);
  SimulationResultsReporter::createPlot("scenario index", "error", x, y, markerSizes, colors);
  //    optimizationResults.saveMeshFiles();
  optimizationResults.draw();
#endif

  return 0;
}