diff --git a/src/main.cpp b/src/main.cpp index 6cf66d0..037610c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -52,7 +52,7 @@ int main(int argc, char *argv[]) { settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; settings_optimization.shouldNormalizeObjectiveValue = true; - settings_optimization.solutionAccuracy = 0.1; + settings_optimization.solutionAccuracy = 0.01; // Optimize pair const std::string pairName = diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 34d904c..402e1e1 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -9,7 +9,7 @@ struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; std::vector> fullPatternDisplacements; - std::vector reducedPatternMaximumDisplacementNormSum; + std::vector fullPatternDisplacementNormSum; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; std::unordered_map @@ -47,13 +47,7 @@ double ReducedModelOptimizer::computeError( 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 / std::max(normalizationFactor, 0.0003); } return rawError; @@ -199,36 +193,11 @@ double ReducedModelOptimizer::objective(long n, const double *x) { std::cout << "Failed simulation" << std::endl; } } else { - const double thisSimulationScenarioError = - computeError(reducedModelResults.displacements, - global.fullPatternDisplacements[simulationScenarioIndex], - global.reducedToFullInterfaceViMap, - global.reducedPatternMaximumDisplacementNormSum - [simulationScenarioIndex]); - if (thisSimulationScenarioError > 1) { - std::cout << "Simulation scenario " - << 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; - } - -#ifdef POLYSCOPE_DEFINED - ReducedModelOptimizer::visualizeResults( - global.fullPatternSimulationJobs[simulationScenarioIndex], - global.reducedPatternSimulationJobs[simulationScenarioIndex], - global.reducedToFullInterfaceViMap, false); - - ReducedModelOptimizer::visualizeResults( - global.fullPatternSimulationJobs[simulationScenarioIndex], - std::make_shared( - reducedPatternMaximumDisplacementSimulationJobs - [simulationScenarioIndex]), - global.reducedToFullInterfaceViMap, true); - polyscope::removeAllStructures(); -#endif // POLYSCOPE_DEFINED - } + const double thisSimulationScenarioError = computeError( + reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, + global.fullPatternDisplacementNormSum[simulationScenarioIndex]); error += thisSimulationScenarioError; } // std::ofstream out(filename, std::ios_base::app); @@ -269,7 +238,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { //} // compute error and return it - global.gObjectiveValueHistory.push_back(error); + // global.gObjectiveValueHistory.push_back(error); // auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), // gObjectiveValueHistory.size()); // std::vector colors(gObjectiveValueHistory.size(), 2); @@ -577,7 +546,7 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( } #if POLYSCOPE_DEFINED -void ReducedModelOptimizer::visualizeResults( +void ReducedModelOptimizer::registerResultsForDrawing( const std::shared_ptr &pFullPatternSimulationJob, const std::shared_ptr &pReducedPatternSimulationJob, const std::unordered_map @@ -609,19 +578,6 @@ void ReducedModelOptimizer::visualizeResults( global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); std::cout << "Error is " << error << std::endl; reducedModelResults.registerForDrawing(); - if (draw) { - polyscope::show(); - const std::string screenshotFilename = - "/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" - "Images/" + - pFullPatternSimulationMesh->getLabel() + "_" + "noScenarioName"; - // simulationScenarioStrings[simulationScenarioIndex]; - polyscope::screenshot(screenshotFilename, false); - fullModelResults.unregister(); - reducedModelResults.unregister(); - pFullPatternSimulationMesh->unregister(); - } } #endif // POLYSCOPE_DEFINED @@ -720,16 +676,37 @@ 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()); + // 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; + // } else { + // results.objectiveValue = + // objective(results.x[0], results.x[1], results.x[2]); + // } + // global.optimizationSettings.shouldNormalizeObjectiveValue = true; + // } + // Compute obj value per simulation scenario + updateMesh(results.x.size(), results.x.data()); + results.objectiveValuePerSimulationScenario.resize( + global.simulationScenarioIndices.size()); + FormFinder::Settings simulationSettings; + FormFinder simulator; + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < + SimulationScenario::NumberOfSimulationScenarios; + simulationScenarioIndex++) { + SimulationResults reducedModelResults = simulator.executeSimulation( + global.reducedPatternSimulationJobs[simulationScenarioIndex], + simulationSettings); + const double error = computeError( + reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, + global.fullPatternDisplacementNormSum[simulationScenarioIndex]); + results.objectiveValuePerSimulationScenario[simulationScenarioIndex] = + error; } results.time = elapsed.count() / 1000.0; const bool printDebugInfo = false; @@ -1044,7 +1021,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.g_optimalReducedModelDisplacements.resize(6); global.reducedPatternSimulationJobs.resize(6); global.fullPatternDisplacements.resize(6); - global.reducedPatternMaximumDisplacementNormSum.resize(6); global.g_firstRoundIterationIndex = 0; global.minY = std::numeric_limits::max(); global.numOfSimulationCrashes = 0; @@ -1052,7 +1028,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.optimizationSettings = optimizationSettings; global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh); - reducedPatternMaximumDisplacementSimulationJobs.resize(6); + global.fullPatternDisplacementNormSum.resize(6); // polyscope::removeAllStructures(); FormFinder::Settings simulationSettings; @@ -1073,51 +1049,17 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( std::make_shared(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 mostFlexibleOptimizationParameters{minB, maxRatio, - minE, minHS}; - if (global.optimizeInnerHexagonSize) { - updateMesh(4, mostFlexibleOptimizationParameters.data()); - } else { - updateMesh(3, mostFlexibleOptimizationParameters.data()); + double displacementNormSum = 0; + for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + const int interfaceFullPatternVi = interfaceViPair.second; + const Eigen::Vector3d displacement( + fullModelResults.displacements[interfaceFullPatternVi][0], + fullModelResults.displacements[interfaceFullPatternVi][1], + fullModelResults.displacements[interfaceFullPatternVi][2]); + displacementNormSum += displacement.norm(); } - 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(reducedModelResults.displacements.size(), - Vector6d(0)), - fullModelResults.displacements, global.reducedToFullInterfaceViMap); - 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; - // } + global.fullPatternDisplacementNormSum[simulationScenarioIndex] = + displacementNormSum; } } Results optResults = runOptimization(optimizationSettings); @@ -1132,9 +1074,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( // visualizeResults(simulationJobs, global.simulationScenarioIndices); #ifdef POLYSCOPE_DEFINED -// visualizeResults( -// global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs, -// global.simulationScenarioIndices, global.reducedToFullInterfaceViMap); + optResults.draw(); #endif // POLYSCOPE_DEFINED return optResults; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index f3f1f81..4709a29 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -95,6 +95,7 @@ public: int numberOfSimulationCrashes{0}; std::vector x; double objectiveValue; + std::vector objectiveValuePerSimulationScenario; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; @@ -189,21 +190,22 @@ public: 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); + const std::string screenshotFilename = + "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/Results/Images/" + + pFullPatternSimulationJob->pMesh->getLabel() + "_" + + simulationScenarioStrings[simulationJobIndex]; + 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; + // double error = computeError( + // reducedModelResults.displacements,fullModelResults.displacements, + // ); + // std::cout << "Error of simulation scenario " + // << + // simulationScenarioStrings[simulationScenarioIndex] + // << " is " + // << error << std::endl; } } #endif // POLYSCOPE_DEFINED @@ -239,7 +241,13 @@ public: void writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization, csvFile &os) { - os << "Obj value"; + os << "Total Obj value"; + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < + SimulationScenario::NumberOfSimulationScenarios; + simulationScenarioIndex++) { + os << "Obj value " + simulationScenarioStrings[simulationScenarioIndex]; + } for (const ReducedModelOptimizer::xRange &range : settings_optimization.xRanges) { os << range.label; @@ -252,6 +260,9 @@ public: writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization, csvFile &os) const { os << objectiveValue; + for (double scenarioObjValue : objectiveValuePerSimulationScenario) { + os << scenarioObjValue; + } for (const double &optimalX : x) { os << optimalX; } @@ -319,7 +330,7 @@ public: &fullToReducedInterfaceViMap, std::unordered_map &fullPatternOppositeInterfaceViMap); - static void visualizeResults( + static void registerResultsForDrawing( const std::shared_ptr &pFullPatternSimulationJob, const std::shared_ptr &pReducedPatternSimulationJob, const std::unordered_map> - &fullPatternSimulationJobs, - const std::vector &simulationScenarios); + void registerResultsForDrawing( + const std::vector> + &fullPatternSimulationJobs, + const std::vector &simulationScenarios); static void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap,