Merged branched containing optional polyscope visualization

CMakeLists.txt

@@ -37,6 +37,17 @@ download_project(PROJ               MYSOURCES
 
 endif()
 
+
+#Polyscope
+download_project(PROJ               POLYSCOPE
+                 GIT_REPOSITORY     https://github.com/nmwsharp/polyscope.git
+                 GIT_TAG	    master
+                 PREFIX             ${CMAKE_CURRENT_SOURCE_DIR}/build/external/
+                 ${UPDATE_DISCONNECTED_IF_AVAILABLE}
+)
+add_subdirectory(${POLYSCOPE_SOURCE_DIR})
+add_compile_definitions(POLYSCOPE_DEFINED)
+
 #dlib
 download_project(PROJ               DLIB
                  GIT_REPOSITORY     https://github.com/davisking/dlib.git

src/main.cpp

@@ -1,8 +1,10 @@
 #include "beamformfinder.hpp"
 #include "csvfile.hpp"
 #include "edgemesh.hpp"
-#include "externvariables.hpp"
 #include "flatpattern.hpp"
+#include "polyscope/curve_network.h"
+#include "polyscope/point_cloud.h"
+#include "polyscope/polyscope.h"
 #include "reducedmodeloptimizer.hpp"
 #include "simulationhistoryplotter.hpp"
 #include "trianglepattterntopology.hpp"
@@ -14,30 +16,23 @@
 #include <string>
 #include <vcg/complex/algorithms/update/position.h>
 
-bool printDebugInfo;
-
 int main(int argc, char *argv[]) {
-  if (argc < 5) {
+  if (argc < 3) {
-    std::cerr << "Specify at least if D(ebug) or R(elease) and the pattern "
+    std::cerr << "Specify at least the two pattern filepaths to be "
-                 "pair to be optimized."
+                 "optimized.Exiting.."
               << std::endl;
     std::terminate();
   }
 
-  if (argv[1] == "D") {
-    printDebugInfo = true;
-  } else {
-    printDebugInfo = false;
-  }
   // Populate the pattern pair to be optimized
   ////Full pattern
-  const std::string filepath_fullPattern = argv[2];
+  const std::string filepath_fullPattern = argv[1];
   FlatPattern fullPattern(filepath_fullPattern);
   fullPattern.setLabel(
       std::filesystem::path(filepath_fullPattern).stem().string());
   fullPattern.scale(0.03);
   ////Reduced pattern
-  const std::string filepath_reducedPattern = argv[3];
+  const std::string filepath_reducedPattern = argv[2];
   FlatPattern reducedPattern(filepath_reducedPattern);
   reducedPattern.setLabel(
       std::filesystem::path(filepath_reducedPattern).stem().string());
@@ -59,9 +54,7 @@ int main(int argc, char *argv[]) {
   // Optimize pair
   const std::string pairName =
       fullPattern.getLabel() + "@" + reducedPattern.getLabel();
-  if (printDebugInfo) {
+
-    std::cout << "Optimizing " << pairName << std::endl;
-  }
   const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
   ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
   optimizer.initializePatterns(fullPattern, reducedPattern, {});
@@ -72,28 +65,20 @@ int main(int argc, char *argv[]) {
   const bool input_resultDirectoryDefined = argc >= 6;
   std::string optimizationResultsDirectory =
       input_resultDirectoryDefined ? argv[5] : "OptimizationResults";
-  //// Get current date for creating the results folder
-  std::time_t now = time(0);
-  std::tm *ltm = std::localtime(&now);
-  std::string currentDate = std::to_string(ltm->tm_mday) + "_" +
-                            std::to_string(1 + ltm->tm_mon) + "_" +
-                            std::to_string(1900 + ltm->tm_year);
-  std::filesystem::path optimizationResultsDirectoryDatePath(
-      std::filesystem::path(optimizationResultsDirectory).append(currentDate));
   if (optimizationResults.numberOfSimulationCrashes != 0) {
     const auto crashedJobsDirPath =
-        std::filesystem::path(optimizationResultsDirectoryDatePath)
+        std::filesystem::path(optimizationResultsDirectory)
             .append("CrashedJobs")
             .append(pairName);
     std::filesystem::create_directories(crashedJobsDirPath);
     optimizationResults.save(crashedJobsDirPath.string());
   } else {
-    std::filesystem::path dirPath_thisOptimization(
+    std::filesystem::path convergedJobsDirPath(
-        std::filesystem::path(optimizationResultsDirectoryDatePath)
+        std::filesystem::path(optimizationResultsDirectory)
             .append("ConvergedJobs")
             .append(pairName));
-    std::filesystem::create_directories(dirPath_thisOptimization);
+    std::filesystem::create_directories(convergedJobsDirPath);
-    optimizationResults.save(dirPath_thisOptimization.string());
+    optimizationResults.save(convergedJobsDirPath.string());
   }
   csvFile csv_results({}, false);
   //    csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)

src/reducedmodeloptimizer.cpp

@@ -10,7 +10,7 @@ struct GlobalOptimizationVariables {
   std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
   std::vector<std::vector<Vector6d>> fullPatternDisplacements;
   std::vector<double> fullPatternDisplacementNormSum;
-  std::vector<SimulationJob> g_fullPatternSimulationJob;
+  std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
   std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
   std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
       reducedToFullInterfaceViMap;
@@ -34,84 +34,6 @@ struct GlobalOptimizationVariables {
   int numberOfOptimizationParameters{3};
   ReducedModelOptimizer::Settings optimizationSettings;
 } global;
-//#pragma omp threadprivate(global)
-
-// struct OptimizationCallback {
-//  double operator()(const size_t &iterations, const Eigen::VectorXd &x,
-//                    const double &fval, Eigen::VectorXd &gradient) const {
-//    // run simulation
-//    //    SimulationResults reducedModelResults =
-//    //        simulator.executeSimulation(reducedModelSimulationJob);
-//    //    reducedModelResults.draw(reducedModelSimulationJob);
-//    gObjectiveValueHistory.push_back(fval);
-//    auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
-//                                   gObjectiveValueHistory.size());
-//    gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
-//    //    const std::string plotImageFilename = "objectivePlot.png";
-//    //    matplot::save(plotImageFilename);
-//    //    if (numberOfOptimizationRounds % 30 == 0) {
-//    //      std::filesystem::copy_file(
-//    //          std::filesystem::path(plotImageFilename),
-//    //          std::filesystem::path("objectivePlot_copy.png"));
-//    //    }
-//    //    std::stringstream ss;
-//    //    ss << x;
-//    //    reducedModelResults.simulationLabel = ss.str();
-//    //    SimulationResultsReporter resultsReporter;
-//    //    resultsReporter.reportResults(
-//    //        {reducedModelResults},
-//    //        std::filesystem::current_path().append("Results"));
-
-//    return true;
-//  }
-//};
-
-// struct Objective {
-//  double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const {
-//    assert(x.rows() == 4);
-
-//    //      drawSimulationJob(simulationJob);
-//    // Set mesh from x
-//    std::shared_ptr<SimulationMesh> reducedModel =
-//        g_reducedPatternSimulationJob.mesh;
-//    for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) {
-//      if (g_reducedPatternExludedEdges.contains(ei)) {
-//        continue;
-//      }
-//      Element &e = reducedModel->elements[ei];
-//      e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * x(0);
-//      e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * x(1);
-//      e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * x(2);
-//      e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * x(3);
-//    }
-//    // run simulation
-//    SimulationResults reducedModelResults =
-//        simulator.executeSimulation(g_reducedPatternSimulationJob);
-//    //    std::stringstream ss;
-//    //    ss << x;
-//    //    reducedModelResults.simulationLabel = ss.str();
-//    //    SimulationResultsReporter resultsReporter;
-//    //    resultsReporter.reportResults(
-//    //        {reducedModelResults},
-//    //        std::filesystem::current_path().append("Results"));
-//    // compute error and return it
-//    double error = 0;
-//    for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
-//      VertexIndex reducedModelVi = reducedFullViPair.first;
-//      Eigen::Vector3d vertexDisplacement(
-//          reducedModelResults.displacements[reducedModelVi][0],
-//          reducedModelResults.displacements[reducedModelVi][1],
-//          reducedModelResults.displacements[reducedModelVi][2]);
-//      Eigen::Vector3d errorVector =
-//          Eigen::Vector3d(
-//              g_optimalReducedModelDisplacements.row(reducedModelVi)) -
-//          vertexDisplacement;
-//      error += errorVector.norm();
-//    }
-
-//    return error;
-//  }
-//};
 
 double ReducedModelOptimizer::computeError(
     const std::vector<Vector6d> &reducedPatternDisplacements,
@@ -204,12 +126,10 @@ void updateMesh(long n, const double *x) {
       pReducedPatternSimulationMesh->vert[vi].P() =
           global.g_innerHexagonVectors[vi / 2] * x[n - 1];
     }
-
-    pReducedPatternSimulationMesh->reset();
-    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
-    //    pReducedPatternSimulationMesh->registerForDrawing("Optimized
-    //    hexagon"); polyscope::show();
   }
+
+  pReducedPatternSimulationMesh->reset();
+  pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
 }
 
 double ReducedModelOptimizer::objective(double b, double h, double E) {
@@ -252,9 +172,12 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
         simulationSettings);
     std::string filename;
     if (!reducedModelResults.converged) {
+
       error += std::numeric_limits<double>::max();
       global.numOfSimulationCrashes++;
-    } else {
+if(beVerbose){
+      std::cout << "Failed simulation" << std::endl;
+    } }else {
       error += computeError(
           reducedModelResults.displacements,
           global.fullPatternDisplacements[simulationScenarioIndex],
@@ -282,6 +205,13 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
         << " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] +
                "\n\n";
     out.close();
+#ifdef POLYSCOPE_DEFINED
+    ReducedModelOptimizer::visualizeResults(
+        global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs,
+        std::vector<SimulationScenario>{
+            static_cast<SimulationScenario>(simulationScenarioIndex)},
+        global.reducedToFullInterfaceViMap);
+#endif // POLYSCOPE_DEFINED
   }
   //  std::cout << error << std::endl;
   if (error < global.minY) {
@@ -414,6 +344,61 @@ void ReducedModelOptimizer::computeMaps(
     assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
     fullPatternOppositeInterfaceViMap[vi0] = vi1;
   }
+
+  const bool debugMapping = false;
+  if (debugMapping) {
+#if POLYSCOPE_DEFINED
+    reducedPattern.registerForDrawing();
+    std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
+        reducedPattern.EN(), glm::vec3(0, 0, 0));
+    for (const size_t ei : global.reducedPatternExludedEdges) {
+      colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
+    }
+    const std::string label = reducedPattern.getLabel();
+    polyscope::getCurveNetwork(label)
+        ->addEdgeColorQuantity("Excluded edges",
+                               colors_reducedPatternExcludedEdges)
+        ->setEnabled(true);
+    polyscope::show();
+
+    std::vector<glm::vec3> nodeColorsOpposite(fullPattern.VN(),
+                                              glm::vec3(0, 0, 0));
+    for (const std::pair<size_t, size_t> oppositeVerts :
+         fullPatternOppositeInterfaceViMap) {
+      auto color = polyscope::getNextUniqueColor();
+      nodeColorsOpposite[oppositeVerts.first] = color;
+      nodeColorsOpposite[oppositeVerts.second] = color;
+    }
+    fullPattern.registerForDrawing();
+    polyscope::getCurveNetwork(fullPattern.getLabel())
+        ->addNodeColorQuantity("oppositeMap", nodeColorsOpposite)
+        ->setEnabled(true);
+    polyscope::show();
+
+    std::vector<glm::vec3> nodeColorsReducedToFull_reduced(reducedPattern.VN(),
+                                                           glm::vec3(0, 0, 0));
+    std::vector<glm::vec3> nodeColorsReducedToFull_full(fullPattern.VN(),
+                                                        glm::vec3(0, 0, 0));
+    for (size_t vi = 0; vi < reducedPattern.VN(); vi++) {
+      if (global.reducedToFullInterfaceViMap.contains(vi)) {
+
+        auto color = polyscope::getNextUniqueColor();
+        nodeColorsReducedToFull_reduced[vi] = color;
+        nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] =
+            color;
+      }
+    }
+    polyscope::getCurveNetwork(reducedPattern.getLabel())
+        ->addNodeColorQuantity("reducedToFull_reduced",
+                               nodeColorsReducedToFull_reduced)
+        ->setEnabled(true);
+    polyscope::getCurveNetwork(fullPattern.getLabel())
+        ->addNodeColorQuantity("reducedToFull_full",
+                               nodeColorsReducedToFull_full)
+        ->setEnabled(true);
+    polyscope::show();
+#endif
+  }
 }
 
 void ReducedModelOptimizer::computeMaps(
@@ -438,6 +423,9 @@ void ReducedModelOptimizer::initializePatterns(
   // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
   assert(fullPattern.VN() == reducedPattern.VN() &&
          fullPattern.EN() >= reducedPattern.EN());
+#if POLYSCOPE_DEFINED
+  polyscope::removeAllStructures();
+#endif
   // Create copies of the input models
   FlatPattern copyFullPattern;
   FlatPattern copyReducedPattern;
@@ -543,6 +531,74 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
   //      reducedModelNodalForcedNormals;
 }
 
+#if POLYSCOPE_DEFINED
+void ReducedModelOptimizer::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) {
+  FormFinder simulator;
+  std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
+      fullPatternSimulationJobs[0]->pMesh;
+  pFullPatternSimulationMesh->registerForDrawing();
+  pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() +
+                                      "_undeformed.ply");
+  reducedPatternSimulationJobs[0]->pMesh->savePly(
+      reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply");
+  double totalError = 0;
+  for (const int simulationScenarioIndex : simulationScenarios) {
+    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
+        fullPatternSimulationJobs[simulationScenarioIndex];
+    pFullPatternSimulationJob->registerForDrawing(
+        pFullPatternSimulationMesh->getLabel());
+    SimulationResults fullModelResults =
+        simulator.executeSimulation(pFullPatternSimulationJob);
+    fullModelResults.registerForDrawing();
+    // fullModelResults.saveDeformedModel();
+    const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+        reducedPatternSimulationJobs[simulationScenarioIndex];
+    SimulationResults reducedModelResults =
+        simulator.executeSimulation(pReducedPatternSimulationJob);
+    double interfaceDisplacementNormSum = 0;
+    for (const auto &interfaceViPair : reducedToFullInterfaceViMap) {
+      const int fullPatternInterfaceIndex = interfaceViPair.second;
+      Eigen::Vector3d fullPatternDisplacementVector(
+          fullModelResults.displacements[fullPatternInterfaceIndex][0],
+          fullModelResults.displacements[fullPatternInterfaceIndex][1],
+          fullModelResults.displacements[fullPatternInterfaceIndex][2]);
+      interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
+    }
+    reducedModelResults.saveDeformedModel();
+    fullModelResults.saveDeformedModel();
+    double error = computeError(
+        reducedModelResults.displacements, fullModelResults.displacements,
+        interfaceDisplacementNormSum, reducedToFullInterfaceViMap);
+    std::cout << "Error of simulation scenario "
+              << simulationScenarioStrings[simulationScenarioIndex] << " is "
+              << error << std::endl;
+    totalError += error;
+    reducedModelResults.registerForDrawing();
+    //    firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
+    //    registerWorldAxes();
+    const std::string screenshotFilename =
+        "/home/iason/Coding/Projects/Approximating shapes with flat "
+        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
+        "Images/" +
+        pFullPatternSimulationMesh->getLabel() + "_" +
+        simulationScenarioStrings[simulationScenarioIndex];
+    polyscope::show();
+    polyscope::screenshot(screenshotFilename, false);
+    fullModelResults.unregister();
+    reducedModelResults.unregister();
+    //    firstOptimizationRoundResults[simulationScenarioIndex].unregister();
+  }
+  std::cout << "Total error:" << totalError << std::endl;
+}
+#endif // POLYSCOPE_DEFINED
+
 void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
     const SimulationResults &fullModelResults,
     const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
@@ -947,8 +1003,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
         SimulationScenario::Saddle};
   }
 
-  std::vector<std::shared_ptr<SimulationJob>> simulationJobs =
-      createScenarios(m_pFullPatternSimulationMesh);
   global.g_optimalReducedModelDisplacements.resize(6);
   global.reducedPatternSimulationJobs.resize(6);
   global.fullPatternDisplacements.resize(6);
@@ -958,13 +1012,15 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   global.numOfSimulationCrashes = 0;
   global.numberOfFunctionCalls = 0;
   global.optimizationSettings = optimizationSettings;
+  global.fullPatternSimulationJobs =
+      createScenarios(m_pFullPatternSimulationMesh);
   //  polyscope::removeAllStructures();
 
   FormFinder::Settings settings;
   //  settings.shouldDraw = true;
   for (int simulationScenarioIndex : global.simulationScenarioIndices) {
     const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
-        simulationJobs[simulationScenarioIndex];
+        global.fullPatternSimulationJobs[simulationScenarioIndex];
     SimulationResults fullModelResults =
         simulator.executeSimulation(pFullPatternSimulationJob, settings);
     global.fullPatternDisplacements[simulationScenarioIndex] =
@@ -997,7 +1053,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   Results optResults = runOptimization(optimizationSettings);
   for (int simulationScenarioIndex : global.simulationScenarioIndices) {
     optResults.fullPatternSimulationJobs.push_back(
-        simulationJobs[simulationScenarioIndex]);
+        global.fullPatternSimulationJobs[simulationScenarioIndex]);
     optResults.reducedPatternSimulationJobs.push_back(
         global.reducedPatternSimulationJobs[simulationScenarioIndex]);
   }

src/reducedmodeloptimizer.hpp

@@ -248,7 +248,52 @@ struct ReducedModelOptimizer::Results {
       }
     }
   }
+#if POLYSCOPE_DEFINED
+  void draw() const {
+    initPolyscope();
+    FormFinder simulator;
+    assert(fullPatternSimulationJobs.size() ==
+           reducedPatternSimulationJobs.size());
+    fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
+    reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
 
+    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();
+      // Drawing of reduced pattern results
+      const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+          reducedPatternSimulationJobs[simulationJobIndex];
+      SimulationResults reducedModelResults =
+          simulator.executeSimulation(pReducedPatternSimulationJob);
+      reducedModelResults.registerForDrawing();
+      polyscope::show();
+      // Save a screensh
+      //    const std::string screenshotFilename =
+      //        "/home/iason/Coding/Projects/Approximating shapes with flat "
+      //        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
+      //        + m_pFullPatternSimulationMesh->getLabel() + "_" +
+      //        simulationScenarioStrings[simulationScenarioIndex];
+      //    polyscope::screenshot(screenshotFilename, false);
+      fullModelResults.unregister();
+      reducedModelResults.unregister();
+      //    double error = computeError(
+      //        reducedModelResults,
+      //        global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
+      //    std::cout << "Error of simulation scenario "
+      //              << simulationScenarioStrings[simulationScenarioIndex] << "
+      //              is "
+      //              << error << std::endl;
+    }
+  }
+#endif // POLYSCOPE_DEFINED
   void saveMeshFiles() const {
     const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
     assert(numberOfSimulationJobs != 0 &&