Refactoring

src/main.cpp

@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
   ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9};
   ReducedModelOptimizer::xRange innerHexagonSize{"HexagonSize", 0.1, 0.9};
   // Test set of full patterns
-  std::string fullPatternsTestSetDirectory = "../TestSet";
+  std::string fullPatternsTestSetDirectory = "../../TestSet";
   if (!std::filesystem::exists(
           std::filesystem::path(fullPatternsTestSetDirectory))) {
     std::cerr << "Full pattern directory does not exist: "
@@ -70,130 +70,108 @@ int main(int argc, char *argv[]) {
       continue;
     }
 
-    FlatPattern *pFullPattern = new FlatPattern(filepathString);
+    for (int reducedPatternIndex = 0;
-    pFullPattern->setLabel(filepath.stem().string());
+         reducedPatternIndex < reducedModels.size(); reducedPatternIndex++) {
-    pFullPattern->scale(0.03);
+      FlatPattern *pFullPattern = new FlatPattern(filepathString);
-    FlatPattern *pReducedPattern = new FlatPattern();
+      pFullPattern->setLabel(filepath.stem().string());
-    pReducedPattern->copy(*reducedModels[1]);
+      pFullPattern->scale(0.03);
-    patternPairs.push_back(std::make_pair(pFullPattern, pReducedPattern));
+      FlatPattern *pReducedPattern = new FlatPattern();
+      pReducedPattern->copy(*reducedModels[reducedPatternIndex]);
+      patternPairs.push_back(std::make_pair(pFullPattern, pReducedPattern));
+    }
   }
 
   //  for (double rangeOffset = 0.15; rangeOffset <= 0.95; rangeOffset += 0.05)
   //  {
-  ReducedModelOptimizer::Settings settings;
+  ReducedModelOptimizer::Settings settings_optimization;
-  for (settings.maxSimulations = 100; settings.maxSimulations < 5000;
+  settings_optimization.xRanges = {beamWidth, beamDimensionsRatio, beamE,
-       settings.maxSimulations += 100) {
+                                   innerHexagonSize};
-    std::string xRangesString = beamWidth.toString() + " " +
-                                beamDimensionsRatio.toString() + " " +
-                                beamE.toString();
-    std::cout << xRangesString << std::endl;
-    settings.xRanges = {beamWidth, beamDimensionsRatio, beamE,
-                        innerHexagonSize};
-    //    std::filesystem::path thisOptimizationDirectory(
-    //        std::filesystem::path("../OptimizationResults").append(xRangesString));
-    //    std::filesystem::create_directories(thisOptimizationDirectory);
-    //    csvfile thisOptimizationStatistics(
-    //        std::filesystem::path(thisOptimizationDirectory)
-    //            .append("statistics.csv")
-    //            .string(),
-    //        true);
 
-    double totalError = 0;
+  //  for (settings_optimization.numberOfFunctionCalls = 100;
-    int totalNumberOfSimulationCrashes = 0;
+  //       settings_optimization.numberOfFunctionCalls < 5000;
-    std::vector<std::pair<std::string, ReducedModelOptimizer::Results>>
+  //       settings_optimization.numberOfFunctionCalls += 100) {
-        resultsPerPattern(patternPairs.size());
+  settings_optimization.numberOfFunctionCalls = 1000;
-    auto start = std::chrono::high_resolution_clock::now();
+  const std::string optimizationSettingsString =
-
+      settings_optimization.toString();
-    //#pragma omp parallel for
+  std::filesystem::path thisOptimizationDirectory(
-    for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
+      std::filesystem::path("../OptimizationResults")
-         patternPairIndex++) {
+          .append(optimizationSettingsString));
-      // const auto filepathString = filepath.string();
+  std::filesystem::create_directories(thisOptimizationDirectory);
-      // Use only the base triangle version
+  std::cout << optimizationSettingsString << std::endl;
-      //    std::cout << "Full pattern:" << filepathString << std::endl;
+  csvFile csv_settings(std::filesystem::path(thisOptimizationDirectory)
-      //    for (int reducedPatternIndex = 0;
+                           .append("settings.csv")
-      //         reducedPatternIndex < reducedModels.size();
-      //         reducedPatternIndex++) {
-      //          FlatPattern *pReducedModel =
-      // reducedModels[reducedPatternIndex];
-      std::unordered_set<size_t> optimizationExcludedEi;
-      //    if (pReducedModel !=
-      //        reducedModels[0]) { // assumes that the singleBar reduced model
-      //                            //              is
-      //      // the
-      //      // first in the reducedModels vector
-      //      optimizationExcludedEi.insert(0);
-      //    }
-      //    FlatPattern cp;
-      //    cp.copy(*reducedModels[0]);
-      const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
-      ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
-      optimizer.initializePatterns(*patternPairs[patternPairIndex].first,
-                                   *patternPairs[patternPairIndex].second,
-                                   optimizationExcludedEi);
-      //    optimizer.optimize({ReducedModelOptimizer::Axial});
-      ReducedModelOptimizer::Results optimizationResults =
-          optimizer.optimize(settings);
-
-      totalError += optimizationResults.objectiveValue;
-      resultsPerPattern[patternPairIndex] =
-          std::make_pair(patternPairs[patternPairIndex].first->getLabel(),
-                         optimizationResults);
-      totalNumberOfSimulationCrashes +=
-          optimizationResults.numberOfSimulationCrashes;
-      //      std::cout << "Have optimized " << ++patternsOptimized << "/"
-      //                << static_cast<int>(
-      //                       std::distance(std::filesystem::directory_iterator(
-      //                                         fullPatternsTestSetDirectory),
-      //                                     std::filesystem::directory_iterator()))
-      //                << " patterns." << std::endl;
-      //    }
-    }
-    auto end = std::chrono::high_resolution_clock::now();
-    auto runtime_ms =
-        std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
-
-    //    for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
-    //         patternPairIndex++) {
-    //        std::filesystem::path
-    //        saveToPath(std::filesystem::path("../OptimizationResults")
-    //                           .append(resultsPerPattern[patternPairIndex].first));
-    //        std::filesystem::create_directory(std::filesystem::path(saveToPath));
-
-    //      resultsPerPattern[patternPairIndex].second.save(saveToPath);
-    //    }
-
-    if (!std::filesystem::exists(
-            std::filesystem::path("OptimizationResults/"))) {
-      std::filesystem::create_directory(
-          std::filesystem::path("OptimizationResults"));
-    }
-
-    csvfile statistics(std::filesystem::path("OptimizationResults")
-                           .append("statistics.csv")
                            .string(),
-                       false);
+                       true);
-    for (const auto &range : settings.xRanges) {
+  settings_optimization.writeTo(csv_settings);
-      statistics << range.min << range.max;
-    }
-    statistics << settings.maxSimulations;
-    if (totalNumberOfSimulationCrashes == 0) {
-      statistics << "No crashes";
-    } else {
-      statistics << totalNumberOfSimulationCrashes;
-    }
 
-    statistics << totalError;
+  double totalError = 0;
-    for (const auto &patternObjectiveValue : resultsPerPattern) {
+  int totalNumberOfSimulationCrashes = 0;
-      statistics << patternObjectiveValue.first;
+  std::vector<ReducedModelOptimizer::Results> optimizationResults_testSet(
-      statistics << patternObjectiveValue.second.objectiveValue;
+      patternPairs.size());
-      for (const double &optimalX : patternObjectiveValue.second.x) {
+  auto start = std::chrono::high_resolution_clock::now();
-        statistics << optimalX;
-      }
-    }
-    statistics << runtime_ms.count() / 1000.0;
 
-    statistics << endrow;
+  //#pragma omp parallel for
+  for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
+       patternPairIndex++) {
+    const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
+    ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
+    optimizer.initializePatterns(*patternPairs[patternPairIndex].first,
+                                 *patternPairs[patternPairIndex].second, {});
+    ReducedModelOptimizer::Results optimizationResults =
+        optimizer.optimize(settings_optimization);
+
+    totalError += optimizationResults.objectiveValue;
+    optimizationResults_testSet[patternPairIndex] = optimizationResults;
+    totalNumberOfSimulationCrashes +=
+        optimizationResults.numberOfSimulationCrashes;
   }
+  auto end = std::chrono::high_resolution_clock::now();
+  auto runtime_ms =
+      std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
+
+  //  for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
+  //       patternPairIndex++) {
+  //    std::filesystem::path fullPatternPath(
+  //        std::filesystem::path(thisOptimizationDirectory)
+  //            .append(patternPairs[patternPairIndex].first->getLabel()));
+  //    std::filesystem::create_directories(std::filesystem::path(fullPatternPath));
+
+  //    std::filesystem::path reducedPatternPath(fullPatternPath.append(
+  //        patternPairs[patternPairIndex].second->getLabel()));
+  //    std::filesystem::create_directories(
+  //        std::filesystem::path(reducedPatternPath));
+
+  //    optimizationResults_testSet[patternPairIndex].save(reducedPatternPath);
+  //  }
+  csvFile statistics(std::filesystem::path(thisOptimizationDirectory)
+                         .append("statistics.csv")
+                         .string(),
+                     false);
+
+  for (const auto &range : settings_optimization.xRanges) {
+    statistics << range.min << range.max;
+  }
+  statistics << settings_optimization.numberOfFunctionCalls;
+  if (totalNumberOfSimulationCrashes == 0) {
+    statistics << "No crashes";
+  } else {
+    statistics << totalNumberOfSimulationCrashes;
+  }
+
+  statistics << totalError;
+  for (int patternPairIndex = 0; patternPairIndex < patternPairs.size();
+       patternPairIndex++) {
+    statistics << patternPairs[patternPairIndex].first->getLabel();
+    statistics << optimizationResults_testSet[patternPairIndex].objectiveValue;
+    for (const double &optimalX :
+         optimizationResults_testSet[patternPairIndex].x) {
+      statistics << optimalX;
+    }
+  }
+  statistics << runtime_ms.count() / 1000.0;
+
+  statistics << endrow;
+  //  }
 
   for (auto patternPair : patternPairs) {
     delete patternPair.first;

src/reducedmodeloptimizer.cpp

@@ -690,12 +690,14 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(
   if (global.optimizeInnerHexagonSize) {
     double (*objF)(double, double, double, double) = &objective;
     result = dlib::find_min_global(
-        objF, xMin, xMax, dlib::max_function_calls(settings.maxSimulations),
+        objF, xMin, xMax,
+        dlib::max_function_calls(settings.numberOfFunctionCalls),
         std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
   } else {
     double (*objF)(double, double, double) = &objective;
     result = dlib::find_min_global(
-        objF, xMin, xMax, dlib::max_function_calls(settings.maxSimulations),
+        objF, xMin, xMax,
+        dlib::max_function_calls(settings.numberOfFunctionCalls),
         std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
   }
   auto end = std::chrono::system_clock::now();
@@ -1026,7 +1028,8 @@ void ReducedModelOptimizer::visualizeResults(
     //    registerWorldAxes();
     const std::string screenshotFilename =
         "/home/iason/Coding/Projects/Approximating shapes with flat "
-        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" +
+        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
+        "Images" +
         m_pFullPatternSimulationMesh->getLabel() + "_" +
         simulationScenarioStrings[simulationScenarioIndex];
     polyscope::show();

src/reducedmodeloptimizer.hpp

@@ -2,6 +2,7 @@
 #define REDUCEDMODELOPTIMIZER_HPP
 
 #include "beamformfinder.hpp"
+#include "csvfile.hpp"
 #include "edgemesh.hpp"
 #include "elementalmesh.hpp"
 #include "matplot/matplot.h"
@@ -47,8 +48,46 @@ public:
 
   struct Settings {
     std::vector<xRange> xRanges;
-    int maxSimulations{100};
+    int numberOfFunctionCalls{100};
     double solutionAccuracy{1e-5};
+
+    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);
+
+      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 << endrow;
+
+      if (!xRanges.empty()) {
+        for (const xRange &range : xRanges) {
+          csv << range.max;
+          csv << range.min;
+        }
+      }
+      csv << numberOfFunctionCalls;
+      csv << solutionAccuracy;
+      csv << endrow;
+    }
   };
 
   inline static const std::string simulationScenarioStrings[] = {

src/reducedmodeloptimizer.hpp.autosave

@@ -0,0 +1,151 @@
+#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 Results {
+    int numberOfSimulationCrashes{0};
+    std::vector<double> x;
+    double objectiveValue;
+  };
+  struct xRange {
+    std::string label;
+    double min;
+    double max;
+    std::string toString() const {
+      return label + "=[" + std::to_string(min) + "," + std::to_string(max) +
+             "]";
+    }
+  };
+
+  struct Settings {
+    std::vector<xRange> xRanges;
+    int numberOfFunctionCalls{100};
+    double solutionAccuracy{1e-5};
+
+    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);
+
+      return settingsString;
+    }
+
+    void writeTo(csvFile &csv) const {
+      // Create settings csv header
+      if (!xRanges.empty()) {
+        for (const xRange &range : xRanges) {
+          csv << range.label + " min";
+          csv << range.label + " max";
+        }
+      }
+      csv << "Function Calls";
+      csv << "Solution Accuracy";
+      csv << endrow;
+
+      if (!xRanges.empty()) {
+        for (const xRange &range : xRanges) {
+          csv << range.min;
+          csv << range.max;
+        }
+      }
+      csv << numberOfFunctionCalls;
+      csv << solutionAccuracy;
+      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);
+
+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,
+                                 double (*pObjectiveFunction)(long,
+                                                              const double *));
+  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;
+};
+
+#endif // REDUCEDMODELOPTIMIZER_HPP