From 76433eaa9dd5b5a075fe82690c7a7a2f94475b2c Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Thu, 20 Oct 2022 13:34:53 +0300 Subject: [PATCH] Refactoring --- reducedmodeloptimizer.cpp | 1014 ++++++++++++++++++++----------------- reducedmodeloptimizer.hpp | 23 +- 2 files changed, 582 insertions(+), 455 deletions(-) diff --git a/reducedmodeloptimizer.cpp b/reducedmodeloptimizer.cpp index a90df3c..d993018 100644 --- a/reducedmodeloptimizer.cpp +++ b/reducedmodeloptimizer.cpp @@ -8,7 +8,10 @@ #include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" +#include "csvfile.hpp" + #ifdef USE_ENSMALLEN +#include #include "ensmallen.hpp" #endif #include @@ -17,7 +20,7 @@ #include #include //#define USE_SCENARIO_WEIGHTS -#include +constexpr bool beVerbose = false; using namespace ReducedModelOptimization; @@ -190,12 +193,11 @@ double ReducedModelOptimizer::objective(const dlib::matrix& x) { return objective(std::vector(x.begin(), x.end()), g_optimizationState); } -#endif -// double ReducedModelOptimizer::objective(const double &xValue) -//{ -// return objective(std::vector({xValue})); -//} +double ReducedModelOptimizer::objective(const double& xValue) { + return objective(std::vector({xValue}), g_optimizationState); +} +#endif double ReducedModelOptimizer::objective( const std::vector& x, @@ -214,12 +216,12 @@ double ReducedModelOptimizer::objective( // e.secondBendingConstFactor // << std::endl; // const int n = x.size(); - std::shared_ptr& pReducedPatternSimulationEdgeMesh = + std::shared_ptr& pReducedModelSimulationEdgeMesh = optimizationState .simulationJobs_reducedModel[optimizationState .simulationScenarioIndices[0]] ->pMesh; - function_updateReducedPattern(x, pReducedPatternSimulationEdgeMesh); + function_updateReducedPattern(x, pReducedModelSimulationEdgeMesh); // optimizationState.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); // optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); // polyscope::show(); @@ -237,7 +239,7 @@ double ReducedModelOptimizer::objective( LinearSimulationModel::label) { pReducedModelSimulator = SimulationModelFactory::create( optimizationState.simulationModelLabel_reducedModel); - pReducedModelSimulator->setStructure(pReducedPatternSimulationEdgeMesh); + pReducedModelSimulator->setStructure(pReducedModelSimulationEdgeMesh); } std::for_each( @@ -339,21 +341,26 @@ double ReducedModelOptimizer::objective( simulationErrorsPerScenario.begin(), simulationErrorsPerScenario.end()); // std::cout << totalError << std::endl; + // if (optimizationState.numberOfFunctionCalls % 100 == 0) { // optimizationState.objectiveValueHistory.push_back(totalError); + // } // optimizationState.plotColors.push_back(10); ++optimizationState.numberOfFunctionCalls; if (totalError < optimizationState.minY) { optimizationState.minY = totalError; - // optimizationState.objectiveValueHistory.push_back(totalError); - // optimizationState.objectiveValueHistory_iteration.push_back(optimizationState.numberOfFunctionCalls); + optimizationState.objectiveValueHistory.push_back(totalError); + optimizationState.objectiveValueHistory_iteration.push_back( + optimizationState.numberOfFunctionCalls); optimizationState.minX.assign(x.begin(), x.end()); #ifdef POLYSCOPE_DEFINED - std::cout << "New best:" << totalError << std::endl; - std::cout.precision(17); - for (int i = 0; i < x.size(); i++) { - std::cout << x[i] << " "; + if (beVerbose) { + std::cout << "New best:" << totalError << std::endl; + std::cout.precision(17); + for (int i = 0; i < x.size(); i++) { + std::cout << x[i] << " "; + } + std::cout << std::endl; } - std::cout << std::endl; #endif // optimizationState.objectiveValueHistoryY.push_back(std::log(totalError)); // optimizationState.objectiveValueHistoryX.push_back(optimizationState.numberOfFunctionCalls); @@ -373,24 +380,10 @@ double ReducedModelOptimizer::objective( } #ifdef POLYSCOPE_DEFINED -#ifdef USE_ENSMALLEN - if (optimizationState.numberOfFunctionCalls % 1000 == 0) { + if (beVerbose && optimizationState.numberOfFunctionCalls % 1000 == 0) { std::cout << "Number of function calls:" << optimizationState.numberOfFunctionCalls << std::endl; } -#else - if (optimizationState.optimizationSettings.dlib.numberOfFunctionCalls >= - 100 && - optimizationState.numberOfFunctionCalls % - (optimizationState.optimizationSettings.dlib - .numberOfFunctionCalls / - 100) == - 0) { - std::cout << "Number of function calls:" - << optimizationState.numberOfFunctionCalls << std::endl; - } - -#endif #endif // compute error and return it return totalError; @@ -829,13 +822,15 @@ void ReducedModelOptimizer::initializeOptimizationParameters( optimizationState.parametersInitialValue[optimizationParameterIndex] = 0; optimizationState - .optimizationInitialValue[optimizationParameterIndex] = 0.5; + .optimizationInitialValue[optimizationParameterIndex] = + initialValue_R; break; case Theta: optimizationState.parametersInitialValue[optimizationParameterIndex] = 0; optimizationState - .optimizationInitialValue[optimizationParameterIndex] = 0; + .optimizationInitialValue[optimizationParameterIndex] = + initialValue_theta; break; } } @@ -1129,10 +1124,14 @@ void ReducedModelOptimizer::getResults( << results.objectiveValue.perSimulationScenario_total[i] << std::endl; std::cout << std::endl; - - reducedModelResults.registerForDrawing(Colors::reducedDeformed, true, 0.01); + const double drawingRadius_pattern = + optimizationState.fullPatternResults[simulationScenarioIndex] + .pJob->pMesh->elements[0] + .dimensions.getDrawingRadius(); + reducedModelResults.registerForDrawing(Colors::reducedDeformed, true, + drawingRadius_pattern); optimizationState.fullPatternResults[simulationScenarioIndex] - .registerForDrawing(Colors::patternDeformed, true, 0.01); + .registerForDrawing(Colors::patternDeformed, true /*, 0.01*/); ChronosEulerNonLinearSimulationModel debug_chronosModel; auto debug_chronosResults = debug_chronosModel.executeSimulation( optimizationState.pSimulationJobs_pattern[simulationScenarioIndex]); @@ -1165,71 +1164,193 @@ void ReducedModelOptimizer::getResults( // results.draw(); } +//#ifdef DLIB_DEFINED +// std::array +// ReducedModelOptimizer::computeFullPatternMaxSimulationForces( +// const std::vector& +// desiredBaseSimulationScenario) const { +// std::array +// fullPatternMaxSimulationForces{0}; +// for (const BaseSimulationScenario& scenario : +// desiredBaseSimulationScenario) { +// const double maxForce = +// computeFullPatternMaxSimulationForce(scenario); +// fullPatternMaxSimulationForces[scenario] = maxForce; +// } + +// return fullPatternMaxSimulationForces; +//} + +// std::array +// ReducedModelOptimizer::getFullPatternMaxSimulationForces( +// const std::vector& +// desiredBaseSimulationScenarioIndices, +// const std::filesystem::path& intermediateResultsDirectoryPath, +// const bool& recomputeForceMagnitudes) { +// std::array +// fullPatternSimulationScenarioMaxMagnitudes; +// //#ifdef POLYSCOPE_DEFINED +// const std::filesystem::path forceMagnitudesDirectoryPath( +// std::filesystem::path(intermediateResultsDirectoryPath) +// .append("ForceMagnitudes")); +// std::filesystem::path patternMaxForceMagnitudesFilePath( +// std::filesystem::path(forceMagnitudesDirectoryPath) +// .append(m_pSimulationEdgeMesh_pattern->getLabel() + ".json")); +// const bool fullPatternScenarioMagnitudesExist = +// std::filesystem::exists(patternMaxForceMagnitudesFilePath); +// if (fullPatternScenarioMagnitudesExist && !recomputeForceMagnitudes) { +// nlohmann::json json; +// std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); +// ifs >> json; +// fullPatternSimulationScenarioMaxMagnitudes = +// static_cast>( +// json.at("maxMagn")); +// return fullPatternSimulationScenarioMaxMagnitudes; +// } + +// //#endif +// fullPatternSimulationScenarioMaxMagnitudes = +// computeFullPatternMaxSimulationForces( +// desiredBaseSimulationScenarioIndices); + +// //#ifdef POLYSCOPE_DEFINED +// nlohmann::json json; +// json["maxMagn"] = fullPatternSimulationScenarioMaxMagnitudes; + +// std::filesystem::create_directories(forceMagnitudesDirectoryPath); +// std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string()); +// jsonFile << json; + +//#ifdef POLYSCOPE_DEFINED +// std::cout << "Saved base scenario max magnitudes to:" +// << patternMaxForceMagnitudesFilePath << std::endl; +//#endif + +// //#endif +// assert(fullPatternSimulationScenarioMaxMagnitudes.size() == +// desiredBaseSimulationScenarioIndices.size()); + +// return fullPatternSimulationScenarioMaxMagnitudes; +//} +//#endif + #ifdef DLIB_DEFINED -std::array -ReducedModelOptimizer::computeFullPatternMaxSimulationForces( - const std::vector& desiredBaseSimulationScenario) - const { - std::array - fullPatternMaxSimulationForces{0}; - for (const BaseSimulationScenario& scenario : desiredBaseSimulationScenario) { - const double maxForce = computeFullPatternMaxSimulationForce(scenario); - fullPatternMaxSimulationForces[scenario] = maxForce; + +ReducedModelOptimizer::FunctionEvaluation +ReducedModelOptimizer::optimizeWith_bobyqa( + const std::vector& parameterGroup, + const Settings& settings) { + g_optimizationState = optimizationState; + const size_t npt = 2 * parameterGroup.size() + 1; + // Set min max values + dlib::matrix xMin_dlib(parameterGroup.size()); + dlib::matrix xMax_dlib(parameterGroup.size()); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + + xMin_dlib(xIndex) = settings.variablesRanges[parameterIndex].min; + xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max; + } + // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; + const double rhobeg = 0.1; + // const double rhobeg = 10; + const double rhoend = settings.solver_accuracy * 1e-10; + // const size_t maxFun = 10 * (x.size() ^ 2); + const size_t bobyqa_maxFunctionCalls = settings.solver_maxIterations; + dlib::matrix x; + dlib::function_evaluation optimalResult_dlib; + optimalResult_dlib.x.set_size(parameterGroup.size()); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + optimalResult_dlib.x(xIndex) = + optimizationState.optimizationInitialValue[parameterIndex]; } - return fullPatternMaxSimulationForces; + double (*objF)(const dlib::matrix&) = &objective; + optimalResult_dlib.y = + dlib::find_min_bobyqa(objF, optimalResult_dlib.x, npt, xMin_dlib, + xMax_dlib, rhobeg, rhoend, bobyqa_maxFunctionCalls); + + FunctionEvaluation parameterGroup_optimalResult; + parameterGroup_optimalResult.x.resize(parameterGroup.size()); + std::copy(g_optimizationState.minX.begin(), g_optimizationState.minX.end(), + parameterGroup_optimalResult.x.begin()); + parameterGroup_optimalResult.y = optimalResult_dlib.y; + optimizationState = g_optimizationState; + + return parameterGroup_optimalResult; } -std::array -ReducedModelOptimizer::getFullPatternMaxSimulationForces( - const std::vector& - desiredBaseSimulationScenarioIndices, - const std::filesystem::path& intermediateResultsDirectoryPath, - const bool& recomputeForceMagnitudes) { - std::array - fullPatternSimulationScenarioMaxMagnitudes; - //#ifdef POLYSCOPE_DEFINED - const std::filesystem::path forceMagnitudesDirectoryPath( - std::filesystem::path(intermediateResultsDirectoryPath) - .append("ForceMagnitudes")); - std::filesystem::path patternMaxForceMagnitudesFilePath( - std::filesystem::path(forceMagnitudesDirectoryPath) - .append(m_pFullPatternSimulationEdgeMesh->getLabel() + ".json")); - const bool fullPatternScenarioMagnitudesExist = - std::filesystem::exists(patternMaxForceMagnitudesFilePath); - if (fullPatternScenarioMagnitudesExist && !recomputeForceMagnitudes) { - nlohmann::json json; - std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); - ifs >> json; - fullPatternSimulationScenarioMaxMagnitudes = - static_cast>( - json.at("maxMagn")); - return fullPatternSimulationScenarioMaxMagnitudes; - } - - //#endif - fullPatternSimulationScenarioMaxMagnitudes = - computeFullPatternMaxSimulationForces( - desiredBaseSimulationScenarioIndices); - - //#ifdef POLYSCOPE_DEFINED - nlohmann::json json; - json["maxMagn"] = fullPatternSimulationScenarioMaxMagnitudes; - - std::filesystem::create_directories(forceMagnitudesDirectoryPath); - std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string()); - jsonFile << json; - +ReducedModelOptimizer::FunctionEvaluation +ReducedModelOptimizer::optimizeWith_dlib( + const std::vector& parameterGroup, + const Settings& settings, + const int& totalNumberOfOptimizationParameters) { + g_optimizationState = optimizationState; #ifdef POLYSCOPE_DEFINED - std::cout << "Saved base scenario max magnitudes to:" - << patternMaxForceMagnitudesFilePath << std::endl; + std::cout << "Optimizing using dlib" << std::endl; #endif + // Set min max values + dlib::matrix xMin_dlib(parameterGroup.size()); + dlib::matrix xMax_dlib(parameterGroup.size()); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; - //#endif - assert(fullPatternSimulationScenarioMaxMagnitudes.size() == - desiredBaseSimulationScenarioIndices.size()); + xMin_dlib(xIndex) = settings.variablesRanges[parameterIndex].min; + xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max; + } + const double portionOfTotalFunctionCallsBudget = [&]() { + // if (settings.optimizationVariablesGroupsWeights.has_value()) { + // assert(settings.optimizationVariablesGroupsWeights->size() == + // settings.optimizationStrategy.size()); + // return settings.optimizationVariablesGroupsWeights + // .value()[parameterGroupIndex]; - return fullPatternSimulationScenarioMaxMagnitudes; + // } else { + return static_cast(parameterGroup.size()) / + totalNumberOfOptimizationParameters; + // } + }(); + const int dlib_maxNumberOfFunctionCalls = 100e3; + const int optimizationNumberOfFunctionCalls = + portionOfTotalFunctionCallsBudget * dlib_maxNumberOfFunctionCalls; + const dlib::function_evaluation optimalResult_dlib = [&]() { + if (parameterGroup.size() == 1) { + dlib::function_evaluation result; + double optimalX = + optimizationState.optimizationInitialValue[parameterGroup[0]]; + double (*singleVariableObjective)(const double&) = &objective; + try { + result.y = dlib::find_min_single_variable( + singleVariableObjective, optimalX, xMin_dlib(0), xMax_dlib(0), + settings.solver_accuracy, optimizationNumberOfFunctionCalls); + } catch (const std::exception& + e) { // reference to the base of a polymorphic object +#ifdef POLYSCOPE_DEFINED + std::cout << e.what() << std::endl; +#endif + } + + result.x.set_size(1); + result.x(0) = optimalX; + return result; + } + double (*objF)(const dlib::matrix&) = &objective; + return dlib::find_min_global( + objF, xMin_dlib, xMax_dlib, + dlib::max_function_calls(optimizationNumberOfFunctionCalls), + std::chrono::hours(24 * 365 * 290), settings.solver_accuracy); + }(); + + FunctionEvaluation parameterGroup_optimalResult; + parameterGroup_optimalResult.x.resize(parameterGroup.size()); + std::copy(optimalResult_dlib.x.begin(), optimalResult_dlib.x.end(), + parameterGroup_optimalResult.x.begin()); + parameterGroup_optimalResult.y = optimalResult_dlib.y; + + optimizationState = g_optimizationState; + + return parameterGroup_optimalResult; } #endif @@ -1240,10 +1361,10 @@ void ReducedModelOptimizer::runOptimization( settings.simulationModelLabel_reducedModel; optimizationState.objectiveValueHistory.clear(); optimizationState.objectiveValueHistory_iteration.clear(); - optimizationState.objectiveValueHistory.reserve( - settings.dlib.numberOfFunctionCalls / 100); - optimizationState.objectiveValueHistory_iteration.reserve( - settings.dlib.numberOfFunctionCalls / 100); + // optimizationState.objectiveValueHistory.reserve( + // settings.dlibSettings.numberOfFunctionCalls / 100); + // optimizationState.objectiveValueHistory_iteration.reserve( + // settings.dlibSettings.numberOfFunctionCalls / 100); optimizationState.minY = std::numeric_limits::max(); optimizationState.numOfSimulationCrashes = 0; optimizationState.numberOfFunctionCalls = 0; @@ -1323,43 +1444,118 @@ void ReducedModelOptimizer::runOptimization( xMin[xIndex] = settings.variablesRanges[parameterIndex].min; xMax[xIndex] = settings.variablesRanges[parameterIndex].max; } -#ifdef USE_ENSMALLEN -#ifdef POLYSCOPE_DEFINED - std::cout << "Optimizing using ensmallen" << std::endl; -#endif + //#ifdef POLYSCOPE_DEFINED + // std::cout << "Optimizing using ensmallen" << std::endl; + //#endif + + optimizationState.xMin = xMin; + optimizationState.xMax = xMax; +#ifndef DLIB_DEFINED arma::mat x(parameterGroup.size(), 1); for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; x(xIndex, 0) = optimizationState.optimizationInitialValue[parameterIndex]; } - - optimizationState.xMin = xMin; - optimizationState.xMax = xMax; - EnsmallenReducedModelOptimizationObjective optimizationFunction( optimizationState); -// optimizationFunction.optimizationState = optimizationState; -// Set min max values -// ens::CNE optimizer; -// ens::DE optimizer; -// ens::SPSA optimizer; -#ifdef USE_PSO - arma::mat xMin_arma(optimizationState.xMin); - arma::mat xMax_arma(optimizationState.xMax); - // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); - ens::LBestPSO optimizer(settings.pso.numberOfParticles, xMin_arma, - xMax_arma, 1e5, 500, settings.solverAccuracy, 2.05, - 2.35); -#else - // ens::SA optimizer; - ens::SA<> optimizer(ens::ExponentialSchedule(), 100000, 10000., 1000, 100, - settings.solverAccuracy, 3, 20, 0.3, 0.3); #endif - // ens::LBestPSO optimizer; - parameterGroup_optimalResult.y = - optimizer.Optimize(optimizationFunction, x); - // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + switch (settings.optimizationTool) { +#ifdef DLIB_DEFINED + // struct SettingsDlibGlobal { + // int numberOfFunctionCalls{100000}; + // } dlibSettings; + case ReducedModelOptimization::Settings::OptimizationTool:: + dlib_globalMin: { + parameterGroup_optimalResult = optimizeWith_dlib( + parameterGroup, settings, totalNumberOfOptimizationParameters); + break; + } + case ReducedModelOptimization::Settings::OptimizationTool::dlib_bobyqa: { + parameterGroup_optimalResult = + optimizeWith_bobyqa(parameterGroup, settings); + break; + } +#else + // case + // ReducedModelOptimization::Settings::OptimizationTool::ens_CNE: { + // ens::CNE optimizer(20000, 20000, 0.5, 0.2, 0.2, + // settings.solver_accuracy); + // parameterGroup_optimalResult.y = + // optimizer.Optimize(optimizationFunction, x); + // break; + // } + case ReducedModelOptimization::Settings::OptimizationTool::ens_DE: { + ens::DE optimizer(1000, settings.solver_maxIterations, 0.6, 0.6, + settings.solver_accuracy); + // 1e-20); + // ens::DE optimizer; + parameterGroup_optimalResult.y = + optimizer.Optimize(optimizationFunction, x); + break; + } + case ReducedModelOptimization::Settings::OptimizationTool::ens_PSO: { + constexpr int numberOfParticles = 100; + constexpr int horizonSize = 350; + ens::LBestPSO optimizer(numberOfParticles, xMin, xMax, + settings.solver_maxIterations, horizonSize, + settings.solver_accuracy); + parameterGroup_optimalResult.y = + optimizer.Optimize(optimizationFunction, x); + break; + } + case ReducedModelOptimization::Settings::OptimizationTool::ens_SA: { + ens::SA optimizer(ens::ExponentialSchedule(), + settings.solver_maxIterations, 10000., 1000, 100, + settings.solver_accuracy, 3, 20, 10, 0.3); + // ens::SA optimizer; + parameterGroup_optimalResult.y = + optimizer.Optimize(optimizationFunction, x); + break; + } + // case + // ReducedModelOptimization::Settings::OptimizationTool::ens_SPSA: + // { + // // ens::SPSA optimizer; + // // ens::SA sa_optimizer; + // // ens::SA sa_optimizer(ens::ExponentialSchedule(), 10, + // 10000., + // // 10000, 100, + // // settings.solverAccuracy, 3, 20, + // 0.3, + // // 0.3); + // // ens::SA optimizer; + // // parameterGroup_optimalResult.y = + // // sa_optimizer.Optimize(optimizationFunction, x); + // ens::SPSA optimizer(0.602, 0.101, 0.16, 0.3, + // settings.solver_maxIterations, + // settings.solver_accuracy); + // // arma::mat + // sa_x(parameterGroup_optimalResult.x); + // // std::cout << "Optimizing using spsa:" << x << + // std::endl; parameterGroup_optimalResult.y = + // optimizer.Optimize(optimizationFunction, x); + // // std::cout << "Results of Optimizing using spsa:" << + // x << + // // std::endl; + // break; + // } + std::terminate(); +#ifdef POLYSCOPE_DEFINED + std::cout << "optimal x:" + << "\n" + << x << std::endl; + std::cout << "optimal objective:" << optimizationFunction.Evaluate(x) + << std::endl; + std::cout << "optimal objective using state var:" + << optimizationFunction.Evaluate(optimizationState.minX) + << std::endl; +#endif +#endif + } + + // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) + // { // if (x[xIndex] > optimizationState.xMax[xIndex]) { // x[xIndex] = optimizationState.xMax[xIndex]; // } else if (x[xIndex] < optimizationState.xMin[xIndex]) { @@ -1367,118 +1563,12 @@ void ReducedModelOptimizer::runOptimization( // } // } -#ifdef POLYSCOPE_DEFINED - std::cout << "optimal x:" - << "\n" - << x << std::endl; - std::cout << "optimal objective:" << optimizationFunction.Evaluate(x) - << std::endl; - std::cout << "optimal objective using state var:" - << optimizationFunction.Evaluate(optimizationState.minX) - << std::endl; -#endif parameterGroup_optimalResult.x.clear(); parameterGroup_optimalResult.x.resize(parameterGroup.size()); // std::copy(x.begin(), x.end(), // parameterGroup_optimalResult.x.begin()); std::copy(optimizationState.minX.begin(), optimizationState.minX.end(), parameterGroup_optimalResult.x.begin()); -#else - g_optimizationState = optimizationState; -#ifdef POLYSCOPE_DEFINED - std::cout << "Optimizing using dlib" << std::endl; -#endif - // Set min max values - dlib::matrix xMin_dlib(parameterGroup.size()); - dlib::matrix xMax_dlib(parameterGroup.size()); - for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { - const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; - - xMin_dlib(xIndex) = settings.variablesRanges[parameterIndex].min; - xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max; - } - const double portionOfTotalFunctionCallsBudget = [&]() { - if (settings.optimizationVariablesGroupsWeights.has_value()) { - assert(settings.optimizationVariablesGroupsWeights->size() == - settings.optimizationStrategy.size()); - return settings.optimizationVariablesGroupsWeights - .value()[parameterGroupIndex]; - - } else { - return static_cast(parameterGroup.size()) / - totalNumberOfOptimizationParameters; - } - }(); - const int optimizationNumberOfFunctionCalls = - portionOfTotalFunctionCallsBudget * settings.dlib.numberOfFunctionCalls; - const dlib::function_evaluation optimalResult_dlib = [&]() { - if (parameterGroup.size() == 1) { - dlib::function_evaluation result; - double optimalX = - optimizationState.optimizationInitialValue[parameterGroup[0]]; - double (*singleVariableObjective)(const double&) = &objective; - try { - result.y = dlib::find_min_single_variable( - singleVariableObjective, optimalX, xMin_dlib(0), xMax_dlib(0), - 1e-5, optimizationNumberOfFunctionCalls); - } catch (const std::exception& - e) { // reference to the base of a polymorphic object -#ifdef POLYSCOPE_DEFINED - std::cout << e.what() << std::endl; -#endif - } - - result.x.set_size(1); - result.x(0) = optimalX; - return result; - } - double (*objF)(const dlib::matrix&) = &objective; - return dlib::find_min_global( - objF, xMin_dlib, xMax_dlib, - dlib::max_function_calls(optimizationNumberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), settings.solverAccuracy); - }(); - // constexpr bool useBOBYQA = false; - // if (useBOBYQA) { - // const size_t npt = 2 * - // optimizationState.numberOfOptimizationParameters; - // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; - // const double rhobeg = 0.1; - // // const double rhobeg = 10; - // const double rhoend = rhobeg * 1e-6; - // // const size_t maxFun = 10 * (x.size() ^ 2); - // const size_t bobyqa_maxFunctionCalls = 10000; - // dlib::matrix x; - // dlib::function_evaluation optimalResult_dlib; - // optimalResult_dlib.x.set_size(parameterGroup.size()); - // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { - // const OptimizationParameterIndex parameterIndex = - // parameterGroup[xIndex]; optimalResult_dlib.x(xIndex) = - // optimizationState.optimizationInitialValue[parameterIndex]; - - // std::cout << "xIndex:" << xIndex << std::endl; - // std::cout << "xInit:" << optimalResult_dlib.x(xIndex) << - // std::endl; - // } - - // optimalResult_dlib.y = dlib::find_min_bobyqa(objF, - // optimalResult_dlib.x, - // npt, - // xMin, - // xMax, - // rhobeg, - // rhoend, - // bobyqa_maxFunctionCalls); - // } - - parameterGroup_optimalResult.x.clear(); - parameterGroup_optimalResult.x.resize(parameterGroup.size()); - std::copy(optimalResult_dlib.x.begin(), optimalResult_dlib.x.end(), - parameterGroup_optimalResult.x.begin()); - parameterGroup_optimalResult.y = optimalResult_dlib.y; - - optimizationState = g_optimizationState; -#endif const auto firstNonNullReducedSimulationJobIt = std::find_if(optimizationState.simulationJobs_reducedModel.begin(), @@ -1493,18 +1583,16 @@ void ReducedModelOptimizer::runOptimization( ->pMesh); // TODO: Check if its ok to update only the changed // parameters #ifdef POLYSCOPE_DEFINED - std::cout << "Optimal x:" - << "\n"; - for (const auto& x : parameterGroup_optimalResult.x) { - std::cout << x << std::endl; + if (beVerbose) { + std::cout << "Optimal x:" + << "\n"; + for (const auto& x : parameterGroup_optimalResult.x) { + std::cout << x << std::endl; + } } std::cout << "optimal objective:" << objective(parameterGroup_optimalResult.x, optimizationState) << std::endl; - // std::cout << "Minima:" << minima << std::endl; - std::cout << "optimizationState MIN:" << optimizationState.minY - << std::endl; - std::cout << "opt res y:" << parameterGroup_optimalResult.y << std::endl; #endif optimization_optimalResult.y = parameterGroup_optimalResult.y; @@ -1629,21 +1717,22 @@ void ReducedModelOptimizer::constructDomeSimulationScenario( // Eigen::Vector3d(interfaceVector[0], // // interfaceVector[1], // // 0) - // // * - // 0.01 /** std::abs(forceMagnitude)*/; + // // * 0.01 /** std::abs(forceMagnitude)*/; // // * std::abs(forceMagnitude / // maxForceMagnitude_dome); - // // CoordType v = (pMesh->vert[viPair.first].cP() - - // pMesh->vert[viPair.second].cP()) - // // ^ CoordType(0, 0, -1).Normalize(); + // // CoordType v = (pMesh->vert[viPair.first].cP() + // - pMesh->vert[viPair.second].cP()) + // // ^ CoordType(0, 0, + // -1).Normalize(); // // nodalForces[viPair.first] = Vector6d({0, 0, 0, // v[0], v[1], 0}) * forceMagnitude // // * 0.0001; - // // nodalForces[viPair.second] = Vector6d({0, 0, 0, - // -v[0], -v[1], 0}) * forceMagnitude + // // nodalForces[viPair.second] = Vector6d({0, 0, + // 0, -v[0], -v[1], 0}) * forceMagnitude // // * 0.0001; - // job.constrainedVertices[viPair.first] = {0, 1, 2}; - // job.constrainedVertices[viPair.second] = {2}; + // job.constrainedVertices[viPair.first] = {0, 1, + // 2}; job.constrainedVertices[viPair.second] = + // {2}; // job.constrainedVertices[viPair.first] = // std::unordered_set{0, 1, 2}; // job.constrainedVertices[viPair.second] = @@ -1726,21 +1815,23 @@ void ReducedModelOptimizer::constructSaddleSimulationScenario( // simulator.executeSimulation(std::make_shared(job)); // // DRMSimulationModel simulator; // // DRMSimulationModel::Settings settings; -// // // settings.threshold_averageResidualToExternalForcesNorm = 1e-2; +// // // settings.threshold_averageResidualToExternalForcesNorm = +// 1e-2; // // settings.totalTranslationalKineticEnergyThreshold = 1e-10; // // // settings.viscousDampingFactor = 1e-2; // // // settings.useKineticDamping = true; // // settings.shouldDraw = false; // // settings.debugModeStep = 200; -// // // settings.threshold_averageResidualToExternalForcesNorm = 1e-5; +// // // settings.threshold_averageResidualToExternalForcesNorm = +// 1e-5; // // settings.maxDRMIterations = 200000; // // SimulationResults results = // simulator.executeSimulation(std::make_shared(job), // // settings); // const double &desiredRotationAngle // = g_baseScenarioMaxDisplacementHelperData -// .desiredMaxRotationAngle; //TODO: move from OptimizationState to -// a new struct +// .desiredMaxRotationAngle; //TODO: move from OptimizationState +// to a new struct // double error; // // job.pMesh->setLabel("initial"); // // job.pMesh->registerForDrawing(); @@ -1766,8 +1857,8 @@ void ReducedModelOptimizer::constructSaddleSimulationScenario( // debugSimulationSettings.threshold_averageResidualToExternalForcesNorm // = 1e-3; // // -// debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = -// true; +// debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold +// = true; // // auto debugResults = // simulator.executeSimulation(std::make_shared(job), // // debugSimulationSettings); @@ -1784,7 +1875,7 @@ void ReducedModelOptimizer::constructSaddleSimulationScenario( // } // // std::filesystem::path // outputPath(std::filesystem::path(saveJobToPath) -// // .append(job.pMesh->getLabel()) +// // .append(job.pMesh->getLabel()) // // .append("mag_" + // optimizationState.currentScenarioName)); // // std::filesystem::create_directories(outputPath); @@ -1850,10 +1941,13 @@ void search(const std::function objectiveFunction, // // settings.totalTranslationalKineticEnergyThreshold = 1e-8; // // settings.viscousDampingFactor = 5e-3; // // settings.useTranslationalKineticEnergyForKineticDamping = true; -// // // settings.threshold_averageResidualToExternalForcesNorm = 1e-5; +// // // settings.threshold_averageResidualToExternalForcesNorm = +// 1e-5; // // // settings.useAverage = true; -// // // settings.totalTranslationalKineticEnergyThreshold = 1e-10; -// // // settings.shouldUseTranslationalKineticEnergyThreshold = true; +// // // settings.totalTranslationalKineticEnergyThreshold = +// 1e-10; +// // // settings.shouldUseTranslationalKineticEnergyThreshold = +// true; // // // settings.shouldDraw = true; // // // settings.isDebugMode = true; // // // settings.drawingStep = 200000; @@ -1904,90 +1998,94 @@ struct ForceMagnitudeOptimization { }; #endif -#ifdef DLIB_DEFINED -double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( - const BaseSimulationScenario& scenario) const { - double forceMagnitude = 1; - double minimumError; - bool wasSuccessful = false; - g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction = - constructBaseScenarioFunctions[scenario]; - const double baseTriangleHeight = vcg::Distance( - baseTriangle.cP(0), (baseTriangle.cP(1) + baseTriangle.cP(2)) / 2); - std::function objectiveFunction; - double translationalOptimizationEpsilon{baseTriangleHeight * 0.001}; - double objectiveEpsilon = translationalOptimizationEpsilon; - objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective; - g_baseScenarioMaxDisplacementHelperData.interfaceViForComputingScenarioError = - m_fullPatternOppositeInterfaceViPairs[1].first; - g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = - 0.1 * baseTriangleHeight; - g_baseScenarioMaxDisplacementHelperData.currentScenarioName = - baseSimulationScenarioNames[scenario]; - g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationEdgeMesh = - m_pFullPatternSimulationEdgeMesh; - double forceMagnitudeEpsilon = 1e-4; +//#ifdef DLIB_DEFINED +// double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( +// const BaseSimulationScenario& scenario) const { +// double forceMagnitude = 1; +// double minimumError; +// bool wasSuccessful = false; +// g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction = +// constructBaseScenarioFunctions[scenario]; +// const double baseTriangleHeight = vcg::Distance( +// baseTriangle.cP(0), (baseTriangle.cP(1) + baseTriangle.cP(2)) / 2); +// std::function objectiveFunction; +// double translationalOptimizationEpsilon{baseTriangleHeight * 0.001}; +// double objectiveEpsilon = translationalOptimizationEpsilon; +// objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective; +// g_baseScenarioMaxDisplacementHelperData.interfaceViForComputingScenarioError +// = +// m_fullPatternOppositeInterfaceViPairs[1].first; +// g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = +// 0.1 * baseTriangleHeight; +// g_baseScenarioMaxDisplacementHelperData.currentScenarioName = +// baseSimulationScenarioNames[scenario]; +// g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationEdgeMesh = +// m_pFullPatternSimulationEdgeMesh; +// double forceMagnitudeEpsilon = 1e-4; - switch (scenario) { - case Axial: - g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = - 0.03 * baseTriangleHeight; - break; - case Shear: - g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = - 0.04 * baseTriangleHeight; - break; - case Bending: - g_baseScenarioMaxDisplacementHelperData - .interfaceViForComputingScenarioError = - m_fullPatternOppositeInterfaceViPairs[0].first; - g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = - 0.05 * baseTriangleHeight; - break; - case Dome: - g_baseScenarioMaxDisplacementHelperData.desiredMaxRotationAngle = - vcg::math::ToRad(20.0); - objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective; - forceMagnitudeEpsilon *= 1e-2; - objectiveEpsilon = vcg::math::ToRad(1.0); - forceMagnitude = 0.6; - break; - case Saddle: - g_baseScenarioMaxDisplacementHelperData - .interfaceViForComputingScenarioError = - m_fullPatternOppositeInterfaceViPairs[0].first; - g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = - 0.05 * baseTriangleHeight; - break; - case S: - g_baseScenarioMaxDisplacementHelperData - .interfaceViForComputingScenarioError = - m_fullPatternOppositeInterfaceViPairs[1].first; - g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = - 0.05 * baseTriangleHeight; - break; - default: - std::cerr << "Unrecognized base scenario" << std::endl; - break; - } +// switch (scenario) { +// case Axial: +// g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = +// 0.03 * baseTriangleHeight; +// break; +// case Shear: +// g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = +// 0.04 * baseTriangleHeight; +// break; +// case Bending: +// g_baseScenarioMaxDisplacementHelperData +// .interfaceViForComputingScenarioError = +// m_fullPatternOppositeInterfaceViPairs[0].first; +// g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = +// 0.05 * baseTriangleHeight; +// break; +// case Dome: +// g_baseScenarioMaxDisplacementHelperData.desiredMaxRotationAngle = +// vcg::math::ToRad(20.0); +// objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective; +// forceMagnitudeEpsilon *= 1e-2; +// objectiveEpsilon = vcg::math::ToRad(1.0); +// forceMagnitude = 0.6; +// break; +// case Saddle: +// g_baseScenarioMaxDisplacementHelperData +// .interfaceViForComputingScenarioError = +// m_fullPatternOppositeInterfaceViPairs[0].first; +// g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = +// 0.05 * baseTriangleHeight; +// break; +// case S: +// g_baseScenarioMaxDisplacementHelperData +// .interfaceViForComputingScenarioError = +// m_fullPatternOppositeInterfaceViPairs[1].first; +// g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = +// 0.05 * baseTriangleHeight; +// break; +// default: +// std::cerr << "Unrecognized base scenario" << std::endl; +// break; +// } - constexpr int maxIterations = 1000; - minimumError = - dlib::find_min_single_variable(objectiveFunction, forceMagnitude, 1e-4, - 1e4, forceMagnitudeEpsilon, maxIterations); +// constexpr int maxIterations = 1000; +// minimumError = +// dlib::find_min_single_variable(objectiveFunction, forceMagnitude, 1e-4, +// 1e4, forceMagnitudeEpsilon, +// maxIterations); #ifdef DLIB_DEFINED #else - // ens::SA<> optimizer; - // arma::vec lowerBound("0.00001"); - // arma::vec upperBound("10000"); - // ens::LBestPSO optimizer(64, lowerBound, upperBound, maxIterations, 350, - // forceMagnitudeEpsilon); ForceMagnitudeOptimization f(objectiveFunction); - // // Create function to be optimized. arma::mat - // forceMagnitude_mat({forceMagnitude}); minimumError = - // optimizer.Optimize(f, forceMagnitude_mat); std::cout << - // ReducedModelOptimization::baseSimulationScenarioNames[scenario] << ": " - // << optimalObjective << std::endl; - // forceMagnitude = forceMagnitude_mat(0, 0); +// ens::SA<> optimizer; +// arma::vec lowerBound("0.00001"); +// arma::vec upperBound("10000"); +// ens::LBestPSO optimizer(64, lowerBound, upperBound, maxIterations, +// 350, forceMagnitudeEpsilon); ForceMagnitudeOptimization +// f(objectiveFunction); +// // Create function to be optimized. arma::mat +// forceMagnitude_mat({forceMagnitude}); minimumError = +// optimizer.Optimize(f, forceMagnitude_mat); std::cout << +// ReducedModelOptimization::baseSimulationScenarioNames[scenario] << ": +// " +// << optimalObjective << std::endl; +// forceMagnitude = forceMagnitude_mat(0, 0); // search(objectiveFunction, // optimizationState.desiredMaxDisplacementValue, @@ -1997,36 +2095,37 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( // maxIterations); #endif - wasSuccessful = minimumError < objectiveEpsilon; +// wasSuccessful = minimumError < objectiveEpsilon; -#ifdef POLYSCOPE_DEFINED - std::cout << "Max " - << ReducedModelOptimization::baseSimulationScenarioNames[scenario] - << " magnitude:" << forceMagnitude << std::endl; - if (!wasSuccessful) { - SimulationJob job; - job.pMesh = optimizationState.pFullPatternSimulationEdgeMesh; - g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( - forceMagnitude, optimizationState.fullPatternOppositeInterfaceViPairs, - job); - std::cerr - << ReducedModelOptimization::baseSimulationScenarioNames[scenario] + - " max scenario magnitude was not succefully determined." - << std::endl; - std::filesystem::path outputPath( - std::filesystem::path("../nonConvergingJobs") - .append(m_pFullPatternSimulationEdgeMesh->getLabel()) - .append("magFinal_" + ReducedModelOptimization:: - baseSimulationScenarioNames[scenario])); - std::filesystem::create_directories(outputPath); - job.save(outputPath); - std::terminate(); - } -#endif +//#ifdef POLYSCOPE_DEFINED +// std::cout << "Max " +// << ReducedModelOptimization::baseSimulationScenarioNames[scenario] +// << " magnitude:" << forceMagnitude << std::endl; +// if (!wasSuccessful) { +// SimulationJob job; +// job.pMesh = optimizationState.pFullPatternSimulationEdgeMesh; +// g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( +// forceMagnitude, optimizationState.fullPatternOppositeInterfaceViPairs, +// job); +// std::cerr << ReducedModelOptimization::baseSimulationScenarioNames[scenario] +// + +// " max scenario magnitude was not succefully determined." +// << std::endl; +// std::filesystem::path outputPath( +// std::filesystem::path("../nonConvergingJobs") +// .append(m_pFullPatternSimulationEdgeMesh->getLabel()) +// .append( +// "magFinal_" + +// ReducedModelOptimization::baseSimulationScenarioNames[scenario])); +// std::filesystem::create_directories(outputPath); +// job.save(outputPath); +// std::terminate(); +//} +//#endif - return forceMagnitude; -} -#endif +// return forceMagnitude; +//} +//#endif void ReducedModelOptimizer::computeScenarioWeights( const std::vector& baseSimulationScenarios, @@ -2073,8 +2172,8 @@ void ReducedModelOptimizer::computeScenarioWeights( optimizationState.scenarioWeights.resize( totalNumberOfSimulationScenarios, 0); // This essentially holds the base scenario weights but I use - // totalNumberOfSimulationScenarios elements instead in order to have - // O(1) access time during the optimization + // totalNumberOfSimulationScenarios elements instead in order to + // have O(1) access time during the optimization for (const BaseSimulationScenario& baseScenario : baseSimulationScenarios) { const int baseSimulationScenarioIndexOffset = std::accumulate( simulationScenariosResolution.begin(), @@ -2153,7 +2252,9 @@ ReducedModelOptimizer::createFullPatternSimulationJobs( } #ifdef POLYSCOPE_DEFINED - std::cout << "Created pattern simulation jobs" << std::endl; + if (beVerbose) { + std::cout << "Created pattern simulation jobs" << std::endl; + } #endif return scenarios; } @@ -2172,8 +2273,8 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors( double translationalDisplacementNormSum = 0; for (auto interfaceViPair : optimizationState.reducedToFullInterfaceViMap) { const int fullPatternVi = interfaceViPair.second; - // If the full pattern vertex is translationally constrained dont take it - // into account + // If the full pattern vertex is translationally constrained dont take + // it into account if (optimizationState.pSimulationJobs_pattern[simulationScenarioIndex] ->constrainedVertices.contains(fullPatternVi)) { const std::unordered_set constrainedDof = @@ -2363,14 +2464,14 @@ void ReducedModelOptimizer::optimize( // .append("DRMSettings") // .append("defaultDRMSimulationSettings.json")); - // simulationSettings.threshold_averageResidualToExternalForcesNorm = 1e-5; - // simulationSettings.viscousDampingFactor = 1e-3; + // simulationSettings.threshold_averageResidualToExternalForcesNorm = + // 1e-5; simulationSettings.viscousDampingFactor = 1e-3; // simulationSettings.beVerbose = true; // simulationSettings.shouldDraw = true; // simulationSettings.isDebugMode = true; // simulationSettings.debugModeStep = 100000; //#ifdef POLYSCOPE_DEFINED - constexpr bool drawPatternSimulationResults = true; + constexpr bool drawPatternSimulationResults = false; auto t1 = std::chrono::high_resolution_clock::now(); //#endif results.wasSuccessful = true; @@ -2380,10 +2481,10 @@ void ReducedModelOptimizer::optimize( std::atomic numberOfComputedJobs = 0; std::for_each( -// std::execution::unseq, -#ifndef POLYSCOPE_DEFINED + // std::execution::unseq, + // #ifndef POLYSCOPE_DEFINED std::execution::par_unseq, -#endif + // #endif optimizationState.simulationScenarioIndices.begin(), optimizationState.simulationScenarioIndices.end(), [&](const int& simulationScenarioIndex) { @@ -2392,78 +2493,82 @@ void ReducedModelOptimizer::optimize( // std::cout << "Executing " << // pFullPatternSimulationJob->getLabel() << std::endl; - // std::filesystem::path jobResultsDirectoryPath( - // std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath) - // .append("FullPatternResults") - // .append(m_pFullPatternSimulationEdgeMesh->getLabel()) - // .append(pPatternSimulationJob->getLabel())); std::unique_ptr pPatternSimulator = SimulationModelFactory::create( optimizationSettings.simulationModelLabel_groundTruth); pPatternSimulator->setStructure(m_pSimulationEdgeMesh_pattern); SimulationResults simulationResults_pattern; - SimulationResults debug_chronosResults; - if (optimizationSettings.simulationModelLabel_groundTruth == - DRMSimulationModel::label) { + // SimulationResults debug_chronosResults; + std::filesystem::path jobResultsDirectoryPath( + std::filesystem::path( + optimizationSettings.intermediateResultsDirectoryPath) + .append("FullPatternResults") + .append(m_pSimulationEdgeMesh_pattern->getLabel()) + .append(pSimulationJob_pattern->getLabel())); + if (!optimizationSettings.intermediateResultsDirectoryPath.empty() && + std::filesystem::exists(jobResultsDirectoryPath)) { + simulationResults_pattern.load(jobResultsDirectoryPath, + pSimulationJob_pattern); + } else { + if (optimizationSettings.simulationModelLabel_groundTruth == + DRMSimulationModel::label) { #ifdef POLYSCOPE_DEFINED - if (drawPatternSimulationResults) { - ChronosEulerNonLinearSimulationModel debug_chronosModel; - debug_chronosResults = - debug_chronosModel.executeSimulation(pSimulationJob_pattern); - debug_chronosResults.registerForDrawing(RGBColor({0, 0, 255}), true, - 0.0001, true); - } + if (drawPatternSimulationResults) { + ChronosEulerNonLinearSimulationModel debug_chronosModel; + const auto debug_chronosResults = + debug_chronosModel.executeSimulation(pSimulationJob_pattern); + debug_chronosResults.registerForDrawing(RGBColor({0, 0, 255}), + true, 0.0001, true); + } #endif - DRMSimulationModel::Settings drmSettings; - drmSettings.threshold_averageResidualToExternalForcesNorm = 1e-5; - drmSettings.threshold_totalTranslationalKineticEnergy = 1e-14; - // drmSettings.threshold_residualForcesMovingAverageDerivativeNorm - // = - // 1e-11; - drmSettings.maxDRMIterations = 10000; - drmSettings.Dtini = 1; - // drmSettings.beVerbose = true; - // drmSettings.debugModeStep = 1e4; - // drmSettings.shouldCreatePlots = true; - // drmSettings.shouldDraw = true; - // drmSettings.threshold_totalTranslationalKineticEnergy = - // 1e-8; - pSimulationJob_pattern->pMesh->getLabel(); - simulationResults_pattern = - dynamic_cast(pPatternSimulator.get()) - ->executeSimulation(pSimulationJob_pattern, drmSettings); - if (!simulationResults_pattern.converged) { -#ifdef POLYSCOPE_DEFINED - std::filesystem::path outputPath( - std::filesystem::path("../nonConvergingJobs") - .append(m_pSimulationEdgeMesh_pattern->getLabel()) - .append("final_" + pSimulationJob_pattern->getLabel())); - std::filesystem::create_directories(outputPath); - pSimulationJob_pattern->save(outputPath); - std::cerr << m_pSimulationEdgeMesh_pattern->getLabel() - << " sim scenario " << pSimulationJob_pattern->getLabel() - << " did not converge on first try." << std::endl; -#endif - - // return; - DRMSimulationModel::Settings drmSettings_secondTry = drmSettings; - drmSettings_secondTry.linearGuessForceScaleFactor = 1; - // drmSettings_secondTry.viscousDampingFactor = 4e-8; - // *drmSettings_secondTry.maxDRMIterations *= 5; - drmSettings_secondTry.maxDRMIterations = 2e6; - drmSettings_secondTry.threshold_totalTranslationalKineticEnergy = - 1e-13; - + DRMSimulationModel::Settings drmSettings; + drmSettings.threshold_averageResidualToExternalForcesNorm = 1e-5; + drmSettings.threshold_totalTranslationalKineticEnergy = 1e-14; + drmSettings.maxDRMIterations = 10000; + // drmSettings.beVerbose = true; + // drmSettings.debugModeStep = 1e4; + // drmSettings.shouldCreatePlots = true; + // drmSettings.shouldDraw = true; + pSimulationJob_pattern->pMesh->getLabel(); simulationResults_pattern = dynamic_cast(pPatternSimulator.get()) - ->executeSimulation(pSimulationJob_pattern, - drmSettings_secondTry); - } + ->executeSimulation(pSimulationJob_pattern, drmSettings); + if (!simulationResults_pattern.converged) { +#ifdef POLYSCOPE_DEFINED + std::filesystem::path outputPath( + std::filesystem::path("../nonConvergingJobs") + .append(m_pSimulationEdgeMesh_pattern->getLabel()) + .append("final_" + pSimulationJob_pattern->getLabel())); + std::filesystem::create_directories(outputPath); + pSimulationJob_pattern->save(outputPath); + if (optimizationSettings.beVerbose) { + std::cerr << m_pSimulationEdgeMesh_pattern->getLabel() + << " sim scenario " + << pSimulationJob_pattern->getLabel() + << " did not converge on first try." << std::endl; + } +#endif - } else { - simulationResults_pattern = - pPatternSimulator->executeSimulation(pSimulationJob_pattern); + // return; + DRMSimulationModel::Settings drmSettings_secondTry = drmSettings; + // drmSettings_secondTry.linearGuessForceScaleFactor = + // 1; drmSettings_secondTry.viscousDampingFactor = + // 4e-8; *drmSettings_secondTry.maxDRMIterations *= 5; + drmSettings_secondTry.maxDRMIterations = 2e6; + drmSettings_secondTry.threshold_totalTranslationalKineticEnergy = + 1e-13; + + simulationResults_pattern = + dynamic_cast(pPatternSimulator.get()) + ->executeSimulation(pSimulationJob_pattern, + drmSettings_secondTry); + } + } else { + simulationResults_pattern = + pPatternSimulator->executeSimulation(pSimulationJob_pattern); + } + simulationResults_pattern.save(jobResultsDirectoryPath); } // ChronosEulerSimulationModel simulationModel_pattern; @@ -2478,8 +2583,8 @@ void ReducedModelOptimizer::optimize( // if // (!patternSimulationResults.isEqual(debug_eulerPatternSimulationResults, // error)) { - // std::cerr << "The computed ground truth does not match - // the IGA model results. " + // std::cerr << "The computed ground truth does not + // match the IGA model results. " // "Error is:" // << error << std::endl; // } @@ -2519,11 +2624,13 @@ void ReducedModelOptimizer::optimize( optimizationState.simulationJobs_reducedModel[simulationScenarioIndex] = std::make_shared(simulationJob_reducedModel); #ifdef POLYSCOPE_DEFINED - std::cout << "Ran sim scenario:" << pSimulationJob_pattern->getLabel() - << ". "; - std::cout << ++numberOfComputedJobs << "/" - << optimizationState.simulationScenarioIndices.size() - << std::endl; + if (beVerbose) { + std::cout << "Ran sim scenario:" << pSimulationJob_pattern->getLabel() + << ". "; + std::cout << ++numberOfComputedJobs << "/" + << optimizationState.simulationScenarioIndices.size() + << std::endl; + } if (drawPatternSimulationResults) { optimizationState.pSimulationJobs_pattern[0] @@ -2536,7 +2643,8 @@ void ReducedModelOptimizer::optimize( simulationResults_pattern.registerForDrawing( ReducedModelOptimization::Colors::patternDeformed, true, 0.0001, true); - // SimulationResults fullPatternResults_linear + // SimulationResults + // fullPatternResults_linear // = // linearSimulator.executeSimulation(pFullPatternSimulationJob); // patternDrmResults.registerForDrawing(std::arraypMesh->unregister(); - // debug_eulerPatternSimulationResults.unregister(); + // debug_chronosResults.unregister(); } #endif }); @@ -2564,8 +2672,10 @@ void ReducedModelOptimizer::optimize( auto t2 = std::chrono::high_resolution_clock::now(); auto s_int = std::chrono::duration_cast(t2 - t1); #ifdef POLYSCOPE_DEFINED - std::cout << "Computed ground of truth pattern results in:" << s_int.count() - << " seconds." << std::endl; + if (beVerbose) { + std::cout << "Computed ground of truth pattern results in:" << s_int.count() + << " seconds." << std::endl; + } if (drawPatternSimulationResults) { optimizationState.pSimulationJobs_pattern[0]->pMesh->unregister(); } @@ -2577,7 +2687,9 @@ void ReducedModelOptimizer::optimize( computeObjectiveValueNormalizationFactors(optimizationSettings); // } #ifdef POLYSCOPE_DEFINED - std::cout << "Running reduced model optimization" << std::endl; + if (beVerbose) { + std::cout << "Running reduced model optimization" << std::endl; + } #endif runOptimization(optimizationSettings, results); } diff --git a/reducedmodeloptimizer.hpp b/reducedmodeloptimizer.hpp index a736970..3f60fbe 100644 --- a/reducedmodeloptimizer.hpp +++ b/reducedmodeloptimizer.hpp @@ -124,7 +124,11 @@ class ReducedModelOptimizer { constexpr static std:: array simulationScenariosResolution = {12, 12, 22, 22, 22, 22}; + // simulationScenariosResolution = {8, 8, 16, 16, 16, 16}; + // simulationScenariosResolution = {4, 4, 8, 8, 8, 8}; + // simulationScenariosResolution = {2, 2, 4, 4, 4, 4}; // simulationScenariosResolution = {2, 2, 2, 22, 22, 22}; + // simulationScenariosResolution = {1, 1, 1, 1, 1, 1}; constexpr static std:: array baseScenarioWeights = {1, 1, 2, 2, 2}; @@ -133,7 +137,7 @@ class ReducedModelOptimizer { simulationScenariosResolution.end(), 0); inline static int fanCardinality{6}; - inline static double initialValue_R{0.3}; + inline static double initialValue_R{0.5}; inline static double initialValue_theta{0}; int interfaceNodeIndex; double operator()(const Eigen::VectorXd& x, Eigen::VectorXd&) const; @@ -284,7 +288,7 @@ class ReducedModelOptimizer { const double& newJ, std::shared_ptr& pReducedPatternSimulationEdgeMesh)> function_updateReducedPattern_material_J; - static double objective(const std::vector& x); + // static double objective(const std::vector& x); void initializeUpdateReducedPatternFunctions(); // static double objective(const double &xValue); @@ -340,11 +344,22 @@ class ReducedModelOptimizer { void getResults(const FunctionEvaluation& optimalObjective, const ReducedModelOptimization::Settings& settings, ReducedModelOptimization::Results& results); - double computeFullPatternMaxSimulationForce( - const ReducedModelOptimization::BaseSimulationScenario& scenario) const; + // double computeFullPatternMaxSimulationForce( + // const ReducedModelOptimization::BaseSimulationScenario& scenario) + // const; #ifdef DLIB_DEFINED static double objective(const dlib::matrix& x); + static double objective(const double& xValue); + ReducedModelOptimizer::FunctionEvaluation optimizeWith_dlib( + const std::vector& + parameterGroup, + const ReducedModelOptimization::Settings& settings, + const int& totalNumberOfOptimizationParameters); + ReducedModelOptimizer::FunctionEvaluation optimizeWith_bobyqa( + const std::vector& + parameterGroup, + const ReducedModelOptimization::Settings& settings); #endif std::array computeFullPatternMaxSimulationForces(