Normalizing using an estimate of the maximum displacement. Removed the R and D command line argument.Refactoring

src/reducedmodeloptimizer.cpp

@@ -9,7 +9,7 @@
 struct GlobalOptimizationVariables {
   std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
   std::vector<std::vector<Vector6d>> fullPatternDisplacements;
-  std::vector<double> fullPatternDisplacementNormSum;
+  std::vector<double> reducedPatternMaximumDisplacementNormSum;
   std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
   std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
   std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@@ -35,26 +35,49 @@ struct GlobalOptimizationVariables {
   ReducedModelOptimizer::Settings optimizationSettings;
 } global;
 
+std::vector<SimulationJob> reducedPatternMaximumDisplacementSimulationJobs;
+
 double ReducedModelOptimizer::computeError(
     const std::vector<Vector6d> &reducedPatternDisplacements,
     const std::vector<Vector6d> &fullPatternDisplacements,
-    const double &interfaceDisplacementsNormSum,
+    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
+        &reducedToFullInterfaceViMap,
+    const double &normalizationFactor) {
+  const double rawError =
+      computeRawError(reducedPatternDisplacements, fullPatternDisplacements,
+                      reducedToFullInterfaceViMap);
+  if (global.optimizationSettings.shouldNormalizeObjectiveValue) {
+    assert(rawError <= normalizationFactor);
+    if (rawError > normalizationFactor) {
+      std::cout << "Normalized error is bigger than one." << std::endl;
+    }
+    return rawError /
+           normalizationFactor; // std::max(interfaceDisplacementsNormSum,
+                                // 0.00003);
+  }
+
+  return rawError;
+}
+
+double ReducedModelOptimizer::computeRawError(
+    const std::vector<Vector6d> &reducedPatternDisplacements,
+    const std::vector<Vector6d> &fullPatternDisplacements,
     const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
         &reducedToFullInterfaceViMap) {
   double error = 0;
   for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
-    VertexIndex reducedModelVi = reducedFullViPair.first;
     //    const auto pos =
     //        g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
     //    std::cout << "Interface vi " << reducedModelVi << " is at position "
     //              << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
     Eigen::Vector3d reducedVertexDisplacement(
-        reducedPatternDisplacements[reducedModelVi][0],
+        reducedPatternDisplacements[reducedFullViPair.first][0],
-        reducedPatternDisplacements[reducedModelVi][1],
+        reducedPatternDisplacements[reducedFullViPair.first][1],
-        reducedPatternDisplacements[reducedModelVi][2]);
+        reducedPatternDisplacements[reducedFullViPair.first][2]);
     if (!std::isfinite(reducedVertexDisplacement[0]) ||
         !std::isfinite(reducedVertexDisplacement[1]) ||
         !std::isfinite(reducedVertexDisplacement[2])) {
+      std::cout << "Displacements are not finite" << std::endl;
       std::terminate();
     }
     Eigen::Vector3d fullVertexDisplacement(
@@ -67,9 +90,6 @@ double ReducedModelOptimizer::computeError(
     error += errorVector.norm();
   }
 
-  if (global.optimizationSettings.normalizeObjectiveValue) {
-    return error / std::max(interfaceDisplacementsNormSum, 0.00003);
-  }
   return error;
 }
 
@@ -132,8 +152,8 @@ void updateMesh(long n, const double *x) {
   pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
 }
 
-double ReducedModelOptimizer::objective(double b, double h, double E) {
+double ReducedModelOptimizer::objective(double b, double r, double E) {
-  std::vector<double> x{b, h, E};
+  std::vector<double> x{b, r, E};
   return ReducedModelOptimizer::objective(x.size(), x.data());
 }
 
@@ -146,12 +166,9 @@ double ReducedModelOptimizer::objective(double b, double h, double E,
 double ReducedModelOptimizer::objective(long n, const double *x) {
   //  std::cout.precision(17);
 
-  //  for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
+  //  const Element &e =
-  //    std::cout << "x[" + std::to_string(parameterIndex) + "]="
+  //  global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout <<
-  //              << x[parameterIndex] << std::endl;
+  //  e.axialConstFactor << " " << e.torsionConstFactor << " "
-  //  }
-  const Element &e = global.reducedPatternSimulationJobs[0]->pMesh->elements[0];
-  //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
   //            << e.firstBendingConstFactor << " " <<
   //            e.secondBendingConstFactor
   //            << std::endl;
@@ -167,51 +184,79 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
   FormFinder::Settings simulationSettings;
   //    simulationSettings.shouldDraw = true;
   for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
+    //    std::cout << "Executing "
+    //              << simulationScenarioStrings[simulationScenarioIndex]
+    //              << std::endl;
     SimulationResults reducedModelResults = simulator.executeSimulation(
         global.reducedPatternSimulationJobs[simulationScenarioIndex],
         simulationSettings);
-    std::string filename;
+    //    std::string filename;
     if (!reducedModelResults.converged) {
-
       error += std::numeric_limits<double>::max();
       global.numOfSimulationCrashes++;
-if(beVerbose){
+      bool beVerbose = false; // TODO: use an extern or data member for that
-      std::cout << "Failed simulation" << std::endl;
+      if (beVerbose) {
-    } }else {
+        std::cout << "Failed simulation" << std::endl;
-      error += computeError(
+      }
-          reducedModelResults.displacements,
+    } else {
-          global.fullPatternDisplacements[simulationScenarioIndex],
+      const double thisSimulationScenarioError =
-          global.fullPatternDisplacementNormSum[simulationScenarioIndex],
+          computeError(reducedModelResults.displacements,
-          global.reducedToFullInterfaceViMap);
+                       global.fullPatternDisplacements[simulationScenarioIndex],
-    }
+                       global.reducedToFullInterfaceViMap,
-    std::ofstream out(filename, std::ios_base::app);
+                       global.reducedPatternMaximumDisplacementNormSum
-    auto pMesh =
+                           [simulationScenarioIndex]);
-        global
+      if (thisSimulationScenarioError > 1) {
-            .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
+        std::cout << "Simulation scenario "
-            ->pMesh;
+                  << simulationScenarioStrings[simulationScenarioIndex]
+                  << " results in an error bigger than one." << std::endl;
+        for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
+          std::cout << "x[" + std::to_string(parameterIndex) + "]="
+                    << x[parameterIndex] << std::endl;
+        }
 
-    for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
-      out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex]
-          << std::endl;
-    }
-    out << pMesh->elements[0].dimensions.toString() + "\n" +
-               pMesh->elements[0].material.toString() + " \nA="
-        << pMesh->elements[0].A << " \nratio="
-        << pMesh->elements[0].dimensions.b / pMesh->elements[0].dimensions.h
-        << " \naxialRig:" << pMesh->elements[0].rigidity.axial
-        << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional
-        << " \nfirstBendingRig:" << pMesh->elements[0].rigidity.firstBending
-        << " \nsecondBendingRig:" << pMesh->elements[0].rigidity.secondBending
-        << " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] +
-               "\n\n";
-    out.close();
 #ifdef POLYSCOPE_DEFINED
-    ReducedModelOptimizer::visualizeResults(
+        ReducedModelOptimizer::visualizeResults(
-        global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs,
+            global.fullPatternSimulationJobs[simulationScenarioIndex],
-        std::vector<SimulationScenario>{
+            global.reducedPatternSimulationJobs[simulationScenarioIndex],
-            static_cast<SimulationScenario>(simulationScenarioIndex)},
+            global.reducedToFullInterfaceViMap, false);
-        global.reducedToFullInterfaceViMap);
+
+        ReducedModelOptimizer::visualizeResults(
+            global.fullPatternSimulationJobs[simulationScenarioIndex],
+            std::make_shared<SimulationJob>(
+                reducedPatternMaximumDisplacementSimulationJobs
+                    [simulationScenarioIndex]),
+            global.reducedToFullInterfaceViMap, true);
+        polyscope::removeAllStructures();
 #endif // POLYSCOPE_DEFINED
+      }
+      error += thisSimulationScenarioError;
+    }
+    //    std::ofstream out(filename, std::ios_base::app);
+    //    auto pMesh =
+    //        global
+    //            .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
+    //            ->pMesh;
+    //    for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++)
+    //    {
+    //      out << "x[" + std::to_string(parameterIndex) + "]=" <<
+    //      x[parameterIndex]
+    //          << std::endl;
+    //    }
+    //    out << pMesh->elements[0].dimensions.toString() + "\n" +
+    //               pMesh->elements[0].material.toString() + " \nA="
+    //        << pMesh->elements[0].A << " \nratio="
+    //        << pMesh->elements[0].dimensions.b /
+    //        pMesh->elements[0].dimensions.h
+    //        << " \naxialRig:" << pMesh->elements[0].rigidity.axial
+    //        << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional
+    //        << " \nfirstBendingRig:" <<
+    //        pMesh->elements[0].rigidity.firstBending
+    //        << " \nsecondBendingRig:" <<
+    //        pMesh->elements[0].rigidity.secondBending
+    //        << " \nscenario:" +
+    //        simulationScenarioStrings[simulationScenarioIndex] +
+    //               "\n\n";
+    //    out.close();
   }
   //  std::cout << error << std::endl;
   if (error < global.minY) {
@@ -533,69 +578,50 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
 
 #if POLYSCOPE_DEFINED
 void ReducedModelOptimizer::visualizeResults(
-    const std::vector<std::shared_ptr<SimulationJob>>
+    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
-        &fullPatternSimulationJobs,
+    const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
-    const std::vector<std::shared_ptr<SimulationJob>>
-        &reducedPatternSimulationJobs,
-    const std::vector<SimulationScenario> &simulationScenarios,
     const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
-        &reducedToFullInterfaceViMap) {
+        &reducedToFullInterfaceViMap,
+    const bool draw = true) {
   FormFinder simulator;
   std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
-      fullPatternSimulationJobs[0]->pMesh;
+      pFullPatternSimulationJob->pMesh;
   pFullPatternSimulationMesh->registerForDrawing();
-  pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() +
+  //  pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel()
-                                      "_undeformed.ply");
+  //  +
-  reducedPatternSimulationJobs[0]->pMesh->savePly(
+  //                                      "_undeformed.ply");
-      reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply");
+  //  reducedPatternSimulationJobs[0]->pMesh->savePly(
-  double totalError = 0;
+  //      reducedPatternSimulationJobs[0]->pMesh->getLabel() +
-  for (const int simulationScenarioIndex : simulationScenarios) {
+  //      "_undeformed.ply");
-    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
+  pFullPatternSimulationJob->registerForDrawing(
-        fullPatternSimulationJobs[simulationScenarioIndex];
+      pFullPatternSimulationMesh->getLabel());
-    pFullPatternSimulationJob->registerForDrawing(
+  SimulationResults fullModelResults =
-        pFullPatternSimulationMesh->getLabel());
+      simulator.executeSimulation(pFullPatternSimulationJob);
-    SimulationResults fullModelResults =
+  fullModelResults.registerForDrawing();
-        simulator.executeSimulation(pFullPatternSimulationJob);
+  // fullModelResults.saveDeformedModel();
-    fullModelResults.registerForDrawing();
+  SimulationResults reducedModelResults =
-    // fullModelResults.saveDeformedModel();
+      simulator.executeSimulation(pReducedPatternSimulationJob);
-    const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
+  //    reducedModelResults.saveDeformedModel();
-        reducedPatternSimulationJobs[simulationScenarioIndex];
+  //    fullModelResults.saveDeformedModel();
-    SimulationResults reducedModelResults =
+  double error = computeRawError(
-        simulator.executeSimulation(pReducedPatternSimulationJob);
+      reducedModelResults.displacements, fullModelResults.displacements,
-    double interfaceDisplacementNormSum = 0;
+      reducedToFullInterfaceViMap/*,
-    for (const auto &interfaceViPair : reducedToFullInterfaceViMap) {
+      global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/);
-      const int fullPatternInterfaceIndex = interfaceViPair.second;
+  std::cout << "Error is " << error << std::endl;
-      Eigen::Vector3d fullPatternDisplacementVector(
+  reducedModelResults.registerForDrawing();
-          fullModelResults.displacements[fullPatternInterfaceIndex][0],
+  if (draw) {
-          fullModelResults.displacements[fullPatternInterfaceIndex][1],
+    polyscope::show();
-          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() + "_" +
+        pFullPatternSimulationMesh->getLabel() + "_" + "noScenarioName";
-        simulationScenarioStrings[simulationScenarioIndex];
+    //        simulationScenarioStrings[simulationScenarioIndex];
-    polyscope::show();
     polyscope::screenshot(screenshotFilename, false);
     fullModelResults.unregister();
     reducedModelResults.unregister();
-    //    firstOptimizationRoundResults[simulationScenarioIndex].unregister();
+    pFullPatternSimulationMesh->unregister();
   }
-  std::cout << "Total error:" << totalError << std::endl;
 }
 #endif // POLYSCOPE_DEFINED
 
@@ -694,7 +720,19 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
   results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
   results.x = global.minX;
   results.objectiveValue = global.minY;
+  if (settings.shouldNormalizeObjectiveValue) {
+    global.optimizationSettings.shouldNormalizeObjectiveValue = false;
+    if (global.optimizeInnerHexagonSize) {
+      results.objectiveValue = objective(4, results.x.data());
+
+    } else {
+      results.objectiveValue =
+          objective(results.x[0], results.x[1], results.x[2]);
+    }
+    global.optimizationSettings.shouldNormalizeObjectiveValue = true;
+  }
   results.time = elapsed.count() / 1000.0;
+  const bool printDebugInfo = false;
   if (printDebugInfo) {
     std::cout << "Finished optimizing." << endl;
     //  std::cout << "Solution x:" << endl;
@@ -1006,7 +1044,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   global.g_optimalReducedModelDisplacements.resize(6);
   global.reducedPatternSimulationJobs.resize(6);
   global.fullPatternDisplacements.resize(6);
-  global.fullPatternDisplacementNormSum.resize(6);
+  global.reducedPatternMaximumDisplacementNormSum.resize(6);
   global.g_firstRoundIterationIndex = 0;
   global.minY = std::numeric_limits<double>::max();
   global.numOfSimulationCrashes = 0;
@@ -1014,34 +1052,18 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   global.optimizationSettings = optimizationSettings;
   global.fullPatternSimulationJobs =
       createScenarios(m_pFullPatternSimulationMesh);
+  reducedPatternMaximumDisplacementSimulationJobs.resize(6);
   //  polyscope::removeAllStructures();
 
-  FormFinder::Settings settings;
+  FormFinder::Settings simulationSettings;
   //  settings.shouldDraw = true;
   for (int simulationScenarioIndex : global.simulationScenarioIndices) {
     const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
         global.fullPatternSimulationJobs[simulationScenarioIndex];
-    SimulationResults fullModelResults =
+    SimulationResults fullModelResults = simulator.executeSimulation(
-        simulator.executeSimulation(pFullPatternSimulationJob, settings);
+        pFullPatternSimulationJob, simulationSettings);
     global.fullPatternDisplacements[simulationScenarioIndex] =
         fullModelResults.displacements;
-    double interfaceDisplacementNormSum = 0;
-    for (const auto &interfaceViPair : global.reducedToFullInterfaceViMap) {
-      const int fullPatternInterfaceIndex = interfaceViPair.second;
-      Eigen::Vector3d fullPatternDisplacementVector(
-          fullModelResults.displacements[fullPatternInterfaceIndex][0],
-          fullModelResults.displacements[fullPatternInterfaceIndex][1],
-          fullModelResults.displacements[fullPatternInterfaceIndex][2]);
-      interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
-    }
-    global.fullPatternDisplacementNormSum[simulationScenarioIndex] =
-        interfaceDisplacementNormSum;
-    // global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize(
-    //    m_pReducedPatternSimulationMesh->VN(), 3);
-    // computeDesiredReducedModelDisplacements(
-    //    fullModelResults, global.reducedToFullInterfaceViMap,
-    //    global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);
-
     SimulationJob reducedPatternSimulationJob;
     reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
     computeReducedModelSimulationJob(*pFullPatternSimulationJob,
@@ -1049,6 +1071,54 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
                                      reducedPatternSimulationJob);
     global.reducedPatternSimulationJobs[simulationScenarioIndex] =
         std::make_shared<SimulationJob>(reducedPatternSimulationJob);
+
+    if (optimizationSettings.shouldNormalizeObjectiveValue) {
+      // Compute the quantities for normalizing the obj value
+      const double minB = optimizationSettings.xRanges[0].min;
+      const double maxRatio = optimizationSettings.xRanges[1].max;
+      const double minE = optimizationSettings.xRanges[2].min;
+      const double minHS = 0.3;
+      std::vector<double> mostFlexibleOptimizationParameters{minB, maxRatio,
+                                                             minE, minHS};
+      if (global.optimizeInnerHexagonSize) {
+        updateMesh(4, mostFlexibleOptimizationParameters.data());
+      } else {
+        updateMesh(3, mostFlexibleOptimizationParameters.data());
+      }
+      reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] =
+          global.reducedPatternSimulationJobs[simulationScenarioIndex]
+              ->getCopy();
+      reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex]
+          .pMesh->setLabel("reduced_maxDisplacement");
+      SimulationResults reducedModelResults = simulator.executeSimulation(
+          global.reducedPatternSimulationJobs[simulationScenarioIndex],
+          simulationSettings);
+      const double errorOfMaxDisplacedReduced = computeRawError(
+          reducedModelResults.displacements, fullModelResults.displacements,
+          global.reducedToFullInterfaceViMap);
+      const double errorOfNonDisplacedReduced = computeRawError(
+          std::vector<Vector6d>(reducedModelResults.displacements.size(),
+                                Vector6d(0)),
+          fullModelResults.displacements, global.reducedToFullInterfaceViMap);
+      const double displacementMultiplier = 2;
+      //      if (global.optimizeInnerHexagonSize) {
+      //        displacementMultiplier = 2;
+      //      } else {
+      //        displacementMultiplier = 2;
+      //      }
+      global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex] =
+          displacementMultiplier *
+          std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced);
+      //      if (errorOfMaxDisplacedReduced > errorOfNonDisplacedReduced) {
+      //        std::cout << "Max disp results in a bigger error for scenario "
+      //                  << simulationScenarioStrings[simulationScenarioIndex]
+      //                  << std::endl;
+      //      } else {
+      //        std::cout << "Zero disp results in a bigger error for scenario "
+      //                  << simulationScenarioStrings[simulationScenarioIndex]
+      //                  << std::endl;
+      //      }
+    }
   }
   Results optResults = runOptimization(optimizationSettings);
   for (int simulationScenarioIndex : global.simulationScenarioIndices) {
@@ -1061,7 +1131,10 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
   //  optResults.draw();
 
   //  visualizeResults(simulationJobs, global.simulationScenarioIndices);
-  // visualizeResults(simulationJobs, global.reducedPatternSimulationJobs,
+#ifdef POLYSCOPE_DEFINED
-  //                 global.simulationScenarioIndices,global.reducedToFullInterfaceViMap);
+//  visualizeResults(
+//      global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs,
+//      global.simulationScenarioIndices, global.reducedToFullInterfaceViMap);
+#endif // POLYSCOPE_DEFINED
   return optResults;
 }