From 233694563ae9ea7a699d818be47f36b8b923ecb5 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Mon, 3 Jan 2022 15:30:13 +0200 Subject: [PATCH] Added per base scenario weights --- src/main.cpp | 347 ++++++++++++++++++++-------------- src/reducedmodeloptimizer.cpp | 287 +++++++++++++++++----------- src/reducedmodeloptimizer.hpp | 35 ++-- 3 files changed, 407 insertions(+), 262 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 3678ddf..37bf2e2 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -4,6 +4,7 @@ #include "reducedmodeloptimizer.hpp" #include "simulationhistoryplotter.hpp" #include "trianglepattterntopology.hpp" +#include #include #include #include @@ -12,7 +13,6 @@ #include #include #include -#include #ifdef POLYSCOPE_DEFINED #include "polyscope/curve_network.h" @@ -21,7 +21,7 @@ #endif int main(int argc, char *argv[]) { - if (argc <= 7) { + if (argc <= 6) { std::cerr << "Wrong number of input parameters. Expects at least 6 input parameters." "Usage:\n" "1)full pattern file path\n" @@ -29,8 +29,8 @@ int main(int argc, char *argv[]) "3)Number of optimizaion function calls\n" "4)Translational error weight\n" "5)Optimization results directory path\n" - "6)[optional]Optimization parameters\n" - "7)[optional]Intermediate results directory path\n" + "6)[optional]Intermediate results directory path\n" + "7)[optional]Optimization parameters\n" "Exiting.." << std::endl; return 1; @@ -49,6 +49,12 @@ int main(int argc, char *argv[]) reducedPattern.scale(0.03, interfaceNodeIndex); // Set the optization settings + const std::filesystem::path optimizationSettingsFilePath = argv[3]; + if (!std::filesystem::exists(optimizationSettingsFilePath)) { + std::cerr << "Input optimization settings file does not exist:" + << optimizationSettingsFilePath << std::endl; + } + ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000}; ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000}; ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000}; @@ -57,65 +63,77 @@ int main(int argc, char *argv[]) ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95}; ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0}; ReducedPatternOptimization::Settings settings_optimization; - settings_optimization.parameterRanges + settings_optimization.variablesRanges = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle}; settings_optimization.numberOfFunctionCalls = std::atoi(argv[3]); const bool input_optimizationParametersDefined = argc >= 8; if (input_optimizationParametersDefined) { const std::string optimizationParametersTag = argv[7]; - std::vector split=Utilities::split(optimizationParametersTag,","); + std::vector split = Utilities::split(optimizationParametersTag, ","); //parse parameter tag std::vector> optimizationParameters; std::vector parameterGroup; for (std::string &s : split) { -// std::cout<(ReducedPatternOptimization::E)))){ + } + if (boost::algorithm::contains(s, + std::to_string( + static_cast(ReducedPatternOptimization::E)))) { parameterGroup.push_back(ReducedPatternOptimization::E); - }else if(boost::algorithm::contains(s,std::to_string(static_cast(ReducedPatternOptimization::A)))){ + } else if (boost::algorithm::contains(s, + std::to_string(static_cast( + ReducedPatternOptimization::A)))) { parameterGroup.push_back(ReducedPatternOptimization::A); - }else if(boost::algorithm::contains(s,std::to_string(static_cast(ReducedPatternOptimization::I2)))){ + } else if (boost::algorithm::contains(s, + std::to_string(static_cast( + ReducedPatternOptimization::I2)))) { parameterGroup.push_back(ReducedPatternOptimization::I2); - }else if(boost::algorithm::contains(s,std::to_string(static_cast(ReducedPatternOptimization::I3)))){ + } else if (boost::algorithm::contains(s, + std::to_string(static_cast( + ReducedPatternOptimization::I3)))) { parameterGroup.push_back(ReducedPatternOptimization::I3); - }else if(boost::algorithm::contains(s,std::to_string(static_cast(ReducedPatternOptimization::J)))){ + } else if (boost::algorithm::contains(s, + std::to_string(static_cast( + ReducedPatternOptimization::J)))) { parameterGroup.push_back(ReducedPatternOptimization::J); - }else if(boost::algorithm::contains(s,std::to_string(static_cast(ReducedPatternOptimization::R)))){ + } else if (boost::algorithm::contains(s, + std::to_string(static_cast( + ReducedPatternOptimization::R)))) { parameterGroup.push_back(ReducedPatternOptimization::R); - }else if(boost::algorithm::contains(s,std::to_string(static_cast(ReducedPatternOptimization::Theta)))){ + } else if (boost::algorithm::contains(s, + std::to_string(static_cast( + ReducedPatternOptimization::Theta)))) { parameterGroup.push_back(ReducedPatternOptimization::Theta); - } - else{ + } else { std::cerr << "Wrong optimization parameter input: " << optimizationParametersTag << std::endl; - - } - if(boost::algorithm::contains(s,"}")){ + } + if (boost::algorithm::contains(s, "}")) { optimizationParameters.push_back(parameterGroup); - } - + } } - settings_optimization.optimizationVariables=optimizationParameters; - #ifdef POLYSCOPE_DEFINED - for(const std::vector& v:settings_optimization.optimizationVariables){ - for(const ReducedPatternOptimization::OptimizationParameterIndex& i:v ){ - std::cout<(i)<<" "; + settings_optimization.optimizationStrategy = optimizationParameters; +#ifdef POLYSCOPE_DEFINED + for (const std::vector &v : + settings_optimization.optimizationStrategy) { + for (const ReducedPatternOptimization::OptimizationParameterIndex &i : v) { + std::cout << static_cast(i) << " "; + } + std::cout << std::endl; } - std::cout<= 7; - if (input_intermediateResultsDirectoryDefined) { - settings_optimization.intermediateResultsDirectoryPath = std::filesystem::path(argv[6]); + if (settings_optimization.optimizationStrategy.empty()) { + std::cerr << "No optimization variables. Exiting.." << std::endl; + std::terminate(); } settings_optimization.normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon; -#ifdef POLYSCOPE_DEFINED - settings_optimization.translationNormalizationParameter = 0; - settings_optimization.rotationNormalizationParameter = 0; -#else - settings_optimization.translationNormalizationParameter = 3e-4; - settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0); -#endif - settings_optimization.solverAccuracy = 1e-1; - settings_optimization.objectiveWeights.translational = std::atof(argv[4]); - settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]); + //#ifdef POLYSCOPE_DEFINED + //#else + settings_optimization.translationNormalizationParameter = 4e-3; + settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0); - std::string xConcatNames; - for (const auto &x : settings_optimization.parameterRanges) { - xConcatNames.append(x.label + "_"); - } - xConcatNames.pop_back(); + // settings_optimization.translationNormalizationParameter = 0; + // settings_optimization.rotationNormalizationParameter = 0; - // Optimize pairthere - const std::string pairName = fullPattern.getLabel(); // + "@" + reducedPattern.getLabel(); - const std::string optimizationName = pairName + "(" - + std::to_string(settings_optimization.numberOfFunctionCalls) - + "_" - + to_string_with_precision( - settings_optimization.objectiveWeights.translational) - + ")" + "_" + xConcatNames; - const std::string optimizationResultsDirectory = argv[5]; - std::string resultsOutputDir; - bool optimizationResultFolderExists = false; - const std::filesystem::path crashedJobsDirPath(std::filesystem::path(optimizationResultsDirectory) - .append("CrashedJobs") - .append(optimizationName)); - if (std::filesystem::exists(crashedJobsDirPath)) { - resultsOutputDir = crashedJobsDirPath.string(); - optimizationResultFolderExists = true; - } - const std::filesystem::path convergedJobsDirPath( - std::filesystem::path(optimizationResultsDirectory) - .append("ConvergedJobs") - .append(optimizationName)); - if (std::filesystem::exists(convergedJobsDirPath)) { - resultsOutputDir = convergedJobsDirPath.string(); - optimizationResultFolderExists = true; - } + // settings_optimization.translationNormalizationParameter = 1e-3; + // settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0); - ReducedPatternOptimization::Results optimizationResults; - constexpr bool shouldReoptimize = true; - bool optimizationAlreadyComputed = false; - if (!shouldReoptimize && optimizationResultFolderExists) { - const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir); - if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) { - optimizationAlreadyComputed = true; - } - } + // settings_optimization.translationNormalizationParameter = 5e-2; + // settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0); - if (!optimizationAlreadyComputed) { - auto start = std::chrono::system_clock::now(); - const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; - assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]); - ReducedModelOptimizer optimizer(numberOfNodesPerSlot); - optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.parameterRanges); - optimizer.optimize(settings_optimization, optimizationResults); - optimizationResults.label = optimizationName; - optimizationResults.baseTriangleFullPattern.copy(fullPattern); - optimizationResults.settings = settings_optimization; - auto end = std::chrono::system_clock::now(); - auto elapsed = std::chrono::duration_cast(end - start); - optimizationResults.time = elapsed.count() / 1000.0; + //#endif + settings_optimization.solverAccuracy = 1e-2; + //TODO refactor that + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial] + .translational + = 1.5; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].rotational + = 2 + - settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial] + .translational; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear] + .translational + = 1.5; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].rotational + = 2 + - settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear] + .translational; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending] + .translational + = 1.2; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending] + .rotational + = 2 + - settings_optimization + .perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending] + .translational; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome] + .translational + = 1.95; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].rotational + = 2 + - settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome] + .translational; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle] + .translational + = 1.2; + settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle] + .rotational + = 2 + - settings_optimization + .perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle] + .translational; - if (!optimizationResults.wasSuccessful) { - // resultsOutputDir = crashedJobsDirPath.string(); - return 1; - } - resultsOutputDir = convergedJobsDirPath.string(); - optimizationResults.save(resultsOutputDir, true); - csvFile csv_resultsLocalFile(std::filesystem::path(resultsOutputDir).append("results.csv"), - true); - csvFile csv_results({}, false); - std::vector csvVector{&csv_resultsLocalFile, &csv_results}; - csv_results << "Name"; - csv_resultsLocalFile << "Name"; - optimizationResults.writeHeaderTo(csvVector); - settings_optimization.writeHeaderTo(csv_results); - csv_results << endrow; - csv_resultsLocalFile << endrow; - csv_results << std::to_string(fullPattern.EN()) + "#" + pairName; - optimizationResults.writeResultsTo(csvVector); - settings_optimization.writeSettingsTo(csv_results); - csv_results << endrow; - csv_resultsLocalFile << endrow; - } + // Optimize pairthere + const std::string pairName = fullPattern.getLabel(); // + "@" + reducedPattern.getLabel(); + const std::string optimizationName + = pairName + "(" + std::to_string(settings_optimization.numberOfFunctionCalls) + "_" + + to_string_with_precision( + settings_optimization + .perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial] + .translational) + + ")"; + const std::string optimizationResultsDirectory = argv[5]; + std::string resultsOutputDir; + bool optimizationResultFolderExists = false; + const std::filesystem::path crashedJobsDirPath( + std::filesystem::path(optimizationResultsDirectory) + .append("CrashedJobs") + .append(optimizationName)); + if (std::filesystem::exists(crashedJobsDirPath)) { + resultsOutputDir = crashedJobsDirPath.string(); + optimizationResultFolderExists = true; + } + const std::filesystem::path convergedJobsDirPath( + std::filesystem::path(optimizationResultsDirectory) + .append("ConvergedJobs") + .append(optimizationName)); + if (std::filesystem::exists(convergedJobsDirPath)) { + resultsOutputDir = convergedJobsDirPath.string(); + optimizationResultFolderExists = true; + } + + ReducedPatternOptimization::Results optimizationResults; + constexpr bool shouldReoptimize = true; + bool optimizationAlreadyComputed = false; + if (!shouldReoptimize && optimizationResultFolderExists) { + const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir); + if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) { + optimizationAlreadyComputed = true; + } + } + + if (!optimizationAlreadyComputed) { + auto start = std::chrono::system_clock::now(); + const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; + assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]); + ReducedModelOptimizer optimizer(numberOfNodesPerSlot); + const bool input_intermediateResultsDirectoryDefined = argc >= 7; + if (input_intermediateResultsDirectoryDefined) { + optimizer.setIntermediateResultsDirectoryPath(std::filesystem::path(argv[6])); + } + optimizer.initializePatterns(fullPattern, + reducedPattern, + settings_optimization.variablesRanges); + optimizer.optimize(settings_optimization, optimizationResults); + optimizationResults.label = optimizationName; + optimizationResults.baseTriangleFullPattern.copy(fullPattern); + optimizationResults.settings = settings_optimization; + auto end = std::chrono::system_clock::now(); + auto elapsed = std::chrono::duration_cast(end - start); + optimizationResults.time = elapsed.count() / 1000.0; + + if (!optimizationResults.wasSuccessful) { + // resultsOutputDir = crashedJobsDirPath.string(); + return 1; + } + resultsOutputDir = convergedJobsDirPath.string(); + optimizationResults.save(resultsOutputDir, true); + csvFile csv_resultsLocalFile(std::filesystem::path(resultsOutputDir).append("results.csv"), + true); + csvFile csv_results({}, false); + std::vector csvVector{&csv_resultsLocalFile, &csv_results}; + csv_results << "Name"; + csv_resultsLocalFile << "Name"; + optimizationResults.writeHeaderTo(csvVector); + settings_optimization.writeHeaderTo(csv_results); + csv_results << endrow; + csv_resultsLocalFile << endrow; + csv_results << std::to_string(fullPattern.EN()) + "#" + pairName; + csv_resultsLocalFile << std::to_string(fullPattern.EN()) + "#" + pairName; + optimizationResults.writeResultsTo(csvVector); + settings_optimization.writeSettingsTo(csv_results); + csv_results << endrow; + csv_resultsLocalFile << endrow; + } #ifdef POLYSCOPE_DEFINED std::vector scenarioLabels( @@ -257,8 +322,8 @@ int main(int argc, char *argv[]) const double colorAxial = 1; const double colorShear = 3; const double colorBending = 5; - const double colorDome = 7; - const double colorSaddle = 9; + const double colorDome = 0.1; + const double colorSaddle = 0; std::vector colors(optimizationResults.objectiveValue.perSimulationScenario_total.size()); for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) { scenarioLabels[scenarioIndex] @@ -283,16 +348,24 @@ int main(int argc, char *argv[]) std::vector y(optimizationResults.objectiveValue.perSimulationScenario_total.size()); for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) { y[scenarioIndex] - = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex] - + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex]; - // y[scenarioIndex] = optimizationResults.objectiveValue - // .perSimulationScenario_total[scenarioIndex]; + // = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex] + // + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex]; + = optimizationResults.objectiveValue.perSimulationScenario_total_unweighted[scenarioIndex]; } std::vector x = matplot::linspace(0, y.size() - 1, y.size()); std::vector markerSizes(y.size(), 5); - SimulationResultsReporter::createPlot("scenario index", "error", x, y, markerSizes, colors); + SimulationResultsReporter::createPlot("scenario index", + "error", + x, + y, + markerSizes, + colors, + std::filesystem::path(resultsOutputDir) + .append("perScenarioObjectiveValues.png")); // optimizationResults.saveMeshFiles(); + std::cout << "Saved results to:" << resultsOutputDir << std::endl; optimizationResults.draw(); + #endif return 0; diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 89c0962..e7d3d1d 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -17,6 +17,7 @@ using namespace ReducedPatternOptimization; struct GlobalOptimizationVariables { + // int totalNumberOfSimulationScenarios; std::vector fullPatternResults; std::vector translationalDisplacementNormalizationValues; std::vector rotationalDisplacementNormalizationValues; @@ -30,10 +31,10 @@ struct GlobalOptimizationVariables std::vector objectiveValueHistory; std::vector plotColors; std::array + ReducedPatternOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables> parametersInitialValue; std::array + ReducedPatternOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables> optimizationInitialValue; std::vector simulationScenarioIndices; double minY{DBL_MAX}; @@ -60,6 +61,7 @@ struct GlobalOptimizationVariables std::vector xMin; std::vector xMax; std::vector scenarioWeights; + std::vector objectiveWeights; } global; double ReducedModelOptimizer::computeDisplacementError( @@ -136,7 +138,8 @@ double ReducedModelOptimizer::computeError( &reducedToFullInterfaceViMap, const double &normalizationFactor_translationalDisplacement, const double &normalizationFactor_rotationalDisplacement, - const double &scenarioWeight) + const double &scenarioWeight, + const Settings::ObjectiveWeights &objectiveWeights) { const double translationalError = computeDisplacementError(simulationResults_fullPattern.displacements, @@ -149,21 +152,21 @@ double ReducedModelOptimizer::computeError( // reducedToFullInterfaceViMap); // std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl; -// std::cout << "trans error:" << translationalError << std::endl; + // std::cout << "trans error:" << translationalError << std::endl; const double rotationalError = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, simulationResults_reducedPattern.rotationalDisplacementQuaternion, reducedToFullInterfaceViMap, normalizationFactor_rotationalDisplacement); - // computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, - // simulationResults_reducedPattern.rotationalDisplacementQuaternion, - // reducedToFullInterfaceViMap); + // const double rotationalError + // = computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, + // simulationResults_reducedPattern.rotationalDisplacementQuaternion, + // reducedToFullInterfaceViMap); // std::cout << "rot error:" << rotationalError<< std::endl; - const double simulationError = global.optimizationSettings.objectiveWeights.translational - * translationalError - + global.optimizationSettings.objectiveWeights.rotational - * rotationalError; - return simulationError * simulationError * scenarioWeight; + const double simulationError = objectiveWeights.translational * translationalError + + objectiveWeights.rotational * rotationalError; + assert(!std::isnan(simulationError)); + return simulationError * simulationError * scenarioWeight * scenarioWeight; } //double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, @@ -236,7 +239,7 @@ double ReducedModelOptimizer::objective(const std::vector &x) // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); // run simulations - std::vector simulationErrorsPerScenario(global.simulationScenarioIndices.size()); + std::vector simulationErrorsPerScenario(totalNumberOfSimulationScenarios, 0); LinearSimulationModel simulator; simulator.setStructure(pReducedPatternSimulationMesh); // simulator.initialize(); @@ -287,21 +290,22 @@ double ReducedModelOptimizer::objective(const std::vector &x) global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex], global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], - scenarioWeight); + scenarioWeight, + global.objectiveWeights[simulationScenarioIndex]); simulationErrorsPerScenario[simulationScenarioIndex] = simulationScenarioError; // } -// #ifdef POLYSCOPE_DEFINED -// reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); -// reducedModelResults.registerForDrawing(Colors::reducedDeformed); -// global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); -// global.fullPatternResults[simulationScenarioIndex].registerForDrawing( -// Colors::fullDeformed); -// polyscope::show(); -// reducedModelResults.unregister(); -// global.pFullPatternSimulationMesh->unregister(); -// global.fullPatternResults[simulationScenarioIndex].unregister(); -// #endif + // #ifdef POLYSCOPE_DEFINED + // reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); + // reducedModelResults.registerForDrawing(Colors::reducedDeformed); + // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); + // global.fullPatternResults[simulationScenarioIndex].registerForDrawing( + // Colors::fullDeformed); + // polyscope::show(); + // reducedModelResults.unregister(); + // global.pFullPatternSimulationMesh->unregister(); + // global.fullPatternResults[simulationScenarioIndex].unregister(); + // #endif } }); // double totalError = std::accumulate(simulationErrorsPerScenario.begin(), @@ -361,7 +365,7 @@ void ReducedModelOptimizer::createSimulationMeshes( { if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { std::cerr << "Error: A rectangular cross section is expected." << std::endl; - terminate(); + std::terminate(); } // Full pattern pFullPatternSimulationMesh = std::make_shared(fullModel); @@ -542,9 +546,10 @@ ReducedModelOptimizer::ReducedModelOptimizer(const std::vector &numberOf scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = true; } -void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern, - PatternGeometry &reducedPattern, - const std::vector &optimizationParameters) +void ReducedModelOptimizer::initializePatterns( + PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, + const std::array &optimizationParameters) { assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN()); fullPatternNumberOfEdges = fullPattern.EN(); @@ -606,6 +611,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() global.functions_updateReducedPatternParameter[E] = [](const double &newE, std::shared_ptr &pReducedPatternSimulationMesh) { + // std::cout << "Updating E with new value:" << newE << std::endl; for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE)); @@ -613,6 +619,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() }; global.functions_updateReducedPatternParameter[A] = [](const double &newA, std::shared_ptr &pReducedPatternSimulationMesh) { + assert(pReducedPatternSimulationMesh != nullptr); const double beamWidth = std::sqrt(newA); const double beamHeight = beamWidth; for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { @@ -648,15 +655,18 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() } void ReducedModelOptimizer::initializeOptimizationParameters( - const std::shared_ptr &mesh, const std::vector &optimizationParameters) + const std::shared_ptr &mesh, + const std::array + &optimizationParameters) { - global.numberOfOptimizationParameters = NumberOfOptimizationParameters; + global.numberOfOptimizationParameters = NumberOfOptimizationVariables; for (int optimizationParameterIndex = 0; optimizationParameterIndex < optimizationParameters.size(); optimizationParameterIndex++) { for (int optimizationParameterIndex = E; - optimizationParameterIndex != NumberOfOptimizationParameters; + optimizationParameterIndex != NumberOfOptimizationVariables; optimizationParameterIndex++) { // const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex]; switch (optimizationParameterIndex) { @@ -764,24 +774,19 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv << std::endl; } //Optimal X values - results.optimalXNameValuePairs.resize(NumberOfOptimizationParameters); + results.optimalXNameValuePairs.resize(NumberOfOptimizationVariables); for (int optimizationParameterIndex = E; - optimizationParameterIndex != NumberOfOptimizationParameters; + optimizationParameterIndex != NumberOfOptimizationVariables; optimizationParameterIndex++) { - if (settings.parameterRanges[optimizationParameterIndex].max - - settings.parameterRanges[optimizationParameterIndex].min - < 1e-10) { - continue; - } if (optimizationParameterIndex != OptimizationParameterIndex::R && optimizationParameterIndex != OptimizationParameterIndex::Theta ) { results.optimalXNameValuePairs[optimizationParameterIndex] - = std::make_pair(settings.parameterRanges[optimizationParameterIndex].label, + = std::make_pair(settings.variablesRanges[optimizationParameterIndex].label, optimalObjective.x[optimizationParameterIndex] * global.parametersInitialValue[optimizationParameterIndex]); } else { //Hex size and angle are pure values (not multipliers of the initial values) results.optimalXNameValuePairs[optimizationParameterIndex] - = std::make_pair(settings.parameterRanges[optimizationParameterIndex].label, + = std::make_pair(settings.variablesRanges[optimizationParameterIndex].label, optimalObjective.x[optimizationParameterIndex]); } @@ -815,6 +820,8 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv global.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_total.resize( global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_total_unweighted.resize( + global.simulationScenarioIndices.size()); //#ifdef POLYSCOPE_DEFINED // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); //#endif @@ -840,7 +847,17 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex], global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], - scenarioWeight); + scenarioWeight, + global.objectiveWeights[simulationScenarioIndex]); + + results.objectiveValue.perSimulationScenario_total_unweighted[i] = computeError( + global.fullPatternResults[simulationScenarioIndex], + reducedModelResults, + global.reducedToFullInterfaceViMap, + global.translationalDisplacementNormalizationValues[simulationScenarioIndex], + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + 1, + global.objectiveWeights[simulationScenarioIndex]); // results.objectiveValue.total += results.objectiveValue.perSimulationScenario_total[i]; //Raw translational @@ -995,9 +1012,9 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, #if POLYSCOPE_DEFINED // global.plotColors.reserve(settings.numberOfFunctionCalls); #endif - assert(!settings.optimizationVariables.empty()); + assert(!settings.optimizationStrategy.empty()); const std::vector> &optimizationParametersGroups - = settings.optimizationVariables; + = settings.optimizationStrategy; #ifndef USE_ENSMALLEN const int totalNumberOfOptimizationParameters @@ -1011,15 +1028,16 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, #endif FunctionEvaluation optimization_optimalResult; - optimization_optimalResult.x.resize(NumberOfOptimizationParameters,0); - for (int optimizationParameterIndex = E; - optimizationParameterIndex != NumberOfOptimizationParameters; - optimizationParameterIndex++) { - optimization_optimalResult.x[optimizationParameterIndex] - = global.parametersInitialValue[optimizationParameterIndex]; - } - for (const std::vector ¶meterGroup : - optimizationParametersGroups) { + optimization_optimalResult.x.resize(NumberOfOptimizationVariables, 0); + for (int optimizationParameterIndex = E; + optimizationParameterIndex != NumberOfOptimizationVariables; + optimizationParameterIndex++) { + optimization_optimalResult.x[optimizationParameterIndex] + = global.optimizationInitialValue[optimizationParameterIndex]; + } + for (int parameterGroupIndex=0;parameterGroupIndex ¶meterGroup = + optimizationParametersGroups[parameterGroupIndex]; FunctionEvaluation parameterGroup_optimalResult; //Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup function_updateReducedPattern = [&](const std::vector &x, @@ -1036,8 +1054,8 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, = global.parametersInitialValue[parameterIndex]; return x[xIndex] * parameterInitialValue; }(); - // std::cout << "Optimization parameter:" << parameterIndex << std::endl; - // std::cout << "New value:" << parameterNewValue << std::endl; + // std::cout << "Optimization parameter:" << parameterIndex << std::endl; + // std::cout << "New value:" << parameterNewValue << std::endl; global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, pMesh); } @@ -1051,8 +1069,8 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; - xMin[xIndex] = settings.parameterRanges[parameterIndex].min; - xMax[xIndex] = settings.parameterRanges[parameterIndex].max; + xMin[xIndex] = settings.variablesRanges[parameterIndex].min; + xMax[xIndex] = settings.variablesRanges[parameterIndex].max; } #ifdef USE_ENSMALLEN arma::mat x(parameterGroup.size(), 1); @@ -1070,14 +1088,21 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, // that is able to handle arbitrary functions. EnsmallenOptimizationObjective optimizationFunction; //Set min max values - // ens::SA optimizer; - // ens::CNE optimizer; - // ens::DE optimizer; - // ens::SPSA optimizer; + // ens::SA optimizer; + // ens::CNE optimizer; + // ens::DE optimizer; + // ens::SPSA optimizer; // arma::mat xMin_arma(global.xMin); // arma::mat xMax_arma(global.xMax); - // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); - ens::LBestPSO optimizer(64, 1, 1, 3000, 350, 1e-1, 2.05, 2.15); + // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); + ens::LBestPSO optimizer(64, + 1, + 1, + settings.numberOfFunctionCalls, + 500, + settings.solverAccuracy, + 2.05, + 2.35); // ens::LBestPSO optimizer; const double minima = optimizer.Optimize(optimizationFunction, x); // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { @@ -1106,12 +1131,22 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; - xMin_dlib(xIndex) = settings.parameterRanges[parameterIndex].min; - xMax_dlib(xIndex) = settings.parameterRanges[parameterIndex].max; + xMin_dlib(xIndex) = settings.variablesRanges[parameterIndex].min; + xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max; } - const int optimizationNumberOfFunctionCalls = static_cast( - (static_cast(parameterGroup.size()) / totalNumberOfOptimizationParameters) - * settings.numberOfFunctionCalls); + 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.numberOfFunctionCalls; const dlib::function_evaluation optimalResult_dlib = [&]() { if (parameterGroup.size() == 1) { dlib::function_evaluation result; @@ -1181,9 +1216,16 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, #endif + const auto firstNonNullReducedSimulationJobIt + = std::find_if(global.reducedPatternSimulationJobs.begin(), + global.reducedPatternSimulationJobs.end(), + [](const std::shared_ptr &pJob) { + return pJob != nullptr; + }); function_updateReducedPattern( - std::vector(parameterGroup_optimalResult.x.begin(), parameterGroup_optimalResult.x.end()), - global.reducedPatternSimulationJobs[0] + std::vector(parameterGroup_optimalResult.x.begin(), + parameterGroup_optimalResult.x.end()), + (*firstNonNullReducedSimulationJobIt) ->pMesh); //TODO: Check if its ok to update only the changed parameters #ifdef POLYSCOPE_DEFINED std::cout << "Optimal x:" @@ -1203,7 +1245,7 @@ optimization_optimalResult.y=parameterGroup_optimalResult.y; optimization_optimalResult.x[parameterIndex] = parameterGroup_optimalResult.x[xIndex]; } } - getResults(optimization_optimalResult, settings, results); + getResults(optimization_optimalResult, settings, results); } void ReducedModelOptimizer::constructAxialSimulationScenario( @@ -1359,7 +1401,7 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni const double &desiredRotationAngle = global.desiredMaxRotationAngle; double error; #ifdef POLYSCOPE_DEFINED - job.pMesh->setLabel("initial"); + // job.pMesh->setLabel("initial"); // job.pMesh->registerForDrawing(); // results.registerForDrawing(); // polyscope::show(); @@ -1587,55 +1629,66 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( } void ReducedModelOptimizer::computeScenarioWeights( - const std::vector &baseSimulationScenarios, - const std::array &maxForceMagnitudePerBaseScenario) + const std::vector &baseSimulationScenarios) { - assert(maxForceMagnitudePerBaseScenario.size() == NumberOfBaseSimulationScenarios); - const double userWeightsSum = std::accumulate(baseScenarioWeights.begin(), - baseScenarioWeights.end(), - 0); - const int numberOfBaseScenarios = baseSimulationScenarios.size(); - - Eigen::VectorXd baseScenario_equalizationWeights(numberOfBaseScenarios); - Eigen::VectorXd baseScenario_userWeights(numberOfBaseScenarios); - Eigen::Matrix baseScenario_maxMagnitudeWeights( - maxForceMagnitudePerBaseScenario.data()); + std::array baseScenario_equalizationWeights; + baseScenario_equalizationWeights.fill(0); + std::array baseScenario_userWeights; + baseScenario_userWeights.fill(0); + //Fill only used base scenario weights for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { - const double baseScenarioEqualizationWeight - = totalNumberOfSimulationScenarios - / static_cast(simulationScenariosResolution[baseScenario]); - baseScenario_equalizationWeights(baseScenario) = baseScenarioEqualizationWeight; - const double baseScenarioUserWeight = static_cast(baseScenarioWeights[baseScenario]) - / userWeightsSum; - baseScenario_userWeights(baseScenario) = baseScenarioUserWeight; + baseScenario_equalizationWeights[baseScenario] = static_cast( + simulationScenariosResolution[baseScenario]); + baseScenario_userWeights[baseScenario] = baseScenarioWeights[baseScenario]; } - baseScenario_maxMagnitudeWeights.normalize(); - baseScenario_equalizationWeights.normalize(); - baseScenario_userWeights.normalize(); + Utilities::normalize(baseScenario_userWeights.begin(), baseScenario_userWeights.end()); + Utilities::normalize(baseScenario_equalizationWeights.begin(), + baseScenario_equalizationWeights.end()); + std::transform(baseScenario_equalizationWeights.begin(), + baseScenario_equalizationWeights.end(), + baseScenario_equalizationWeights.begin(), + [](const double &d) { + if (d == 0) { + return d; + } + return 1 / d; + }); + Utilities::normalize(baseScenario_equalizationWeights.begin(), + baseScenario_equalizationWeights.end()); + + std::array baseScenarios_weights; + baseScenarios_weights.fill(0); + for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { + baseScenarios_weights[baseScenario] = baseScenario_equalizationWeights[baseScenario] + + 2 * baseScenario_userWeights[baseScenario]; + } + Utilities::normalize(baseScenarios_weights.begin(), baseScenarios_weights.end()); + + global.objectiveWeights.resize(totalNumberOfSimulationScenarios); global.scenarioWeights.resize( - totalNumberOfSimulationScenarios); //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, + 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 for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { const int baseSimulationScenarioIndexOffset = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.begin() + baseScenario, 0); - const double baseScenarioEqualizationWeight = baseScenario_equalizationWeights(baseScenario); - // const double baseScenarioUserWeight = baseScenario_userWeights(baseScenario); - const double baseScenarioUserWeight = baseScenario_userWeights(baseScenario); - - const double baseScenarioMaxMagnitudeWeight = baseScenario_maxMagnitudeWeights(baseScenario); for (int simulationScenarioIndex = 0; simulationScenarioIndex < simulationScenariosResolution[baseScenario]; simulationScenarioIndex++) { const int globalScenarioIndex = baseSimulationScenarioIndexOffset + simulationScenarioIndex; - const double scenarioWeight = /*baseScenarioUserWeight +*/ baseScenarioEqualizationWeight - + baseScenarioMaxMagnitudeWeight; - global.scenarioWeights[globalScenarioIndex] = scenarioWeight; + global.scenarioWeights[globalScenarioIndex] = baseScenarios_weights[baseScenario]; + + global.objectiveWeights[globalScenarioIndex] + = global.optimizationSettings.perBaseScenarioObjectiveWeights[baseScenario]; } } + + //Dont normalize since we want the base scenarios to be normalized not the "global scenarios" + // Utilities::normalize(global.scenarioWeights.begin(), global.scenarioWeights.end()); } std::vector> ReducedModelOptimizer::createFullPatternSimulationJobs( @@ -1771,12 +1824,25 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() } } +void ReducedModelOptimizer::setIntermediateResultsDirectoryPath( + const std::filesystem::path &newIntermediateResultsDirectoryPath) +{ + intermediateResultsDirectoryPath = newIntermediateResultsDirectoryPath; +} + void ReducedModelOptimizer::optimize( const Settings &optimizationSettings, ReducedPatternOptimization::Results &results, const std::vector &desiredBaseSimulationScenarioIndices) { - assert(!optimizationSettings.optimizationVariables.empty()); + assert(!optimizationSettings.optimizationStrategy.empty()); + assert(!optimizationSettings.optimizationVariablesGroupsWeights.has_value() + || (optimizationSettings.optimizationStrategy.size() + == optimizationSettings.optimizationVariablesGroupsWeights->size() + && std::accumulate(optimizationSettings.optimizationVariablesGroupsWeights->begin(), + optimizationSettings.optimizationVariablesGroupsWeights->end(), + 0))); + for (int baseSimulationScenarioIndex : desiredBaseSimulationScenarioIndices) { //Increase the size of the vector holding the simulation scenario indices global.simulationScenarioIndices.resize( @@ -1792,6 +1858,7 @@ void ReducedModelOptimizer::optimize( 0)); } + // global.totalNumberOfSimulationScenarios = totalNumberOfSimulationScenarios; global.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios); global.fullPatternResults.resize(totalNumberOfSimulationScenarios); global.translationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios); @@ -1802,9 +1869,6 @@ void ReducedModelOptimizer::optimize( global.optimizationSettings = optimizationSettings; global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh; - const std::filesystem::path intermediateResultsDirectoryPath( - optimizationSettings.intermediateResultsDirectoryPath); - std::array fullPatternSimulationScenarioMaxMagnitudes = getFullPatternMaxSimulationForces(desiredBaseSimulationScenarioIndices, intermediateResultsDirectoryPath); @@ -1813,14 +1877,14 @@ void ReducedModelOptimizer::optimize( fullPatternSimulationScenarioMaxMagnitudes); // polyscope::removeAllStructures(); - computeScenarioWeights(desiredBaseSimulationScenarioIndices, - fullPatternSimulationScenarioMaxMagnitudes); + // ComputeScenarioWeights({Axial, Shear, Bending, Dome, Saddle}, + computeScenarioWeights(desiredBaseSimulationScenarioIndices); results.baseTriangle = global.baseTriangle; DRMSimulationModel::Settings simulationSettings; simulationSettings.maxDRMIterations = 200000; - simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; + simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-10; simulationSettings.viscousDampingFactor = 5e-3; simulationSettings.useKineticDamping = true; @@ -1848,8 +1912,9 @@ void ReducedModelOptimizer::optimize( .append(m_pFullPatternSimulationMesh->getLabel()) .append(pFullPatternSimulationJob->getLabel())); // .append(pFullPatternSimulationJob->getLabel() + ".json") + constexpr bool recomputeFullPatternResults = false; SimulationResults fullPatternResults; - if (std::filesystem::exists(jobResultsDirectoryPath)) { + if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) { fullPatternResults.load(std::filesystem::path(jobResultsDirectoryPath).append("Results"), pFullPatternSimulationJob); } else { diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 9376e5e..cdc1e22 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -41,6 +41,7 @@ class ReducedModelOptimizer std::vector scenarioIsSymmetrical; int fullPatternNumberOfEdges; constexpr static double youngsModulus{1 * 1e9}; + std::filesystem::path intermediateResultsDirectoryPath; public: struct FunctionEvaluation @@ -64,13 +65,12 @@ public: // }; // inline constexpr static ParameterLabels parameterLabels(); - inline static std::array + inline static std::array parameterLabels = {"R", "A", "I2", "I3", "J", "Theta", "R"}; - constexpr static std::array + constexpr static std::array simulationScenariosResolution = {11, 11, 20, 20, 20}; - constexpr static std::array - baseScenarioWeights = {1, 1, 100, 50, 50}; - // constexpr static std::array simulationScenariosResolution = {3, 3, 3, 3, 3}; + constexpr static std::array + baseScenarioWeights = {1, 1, 2, 3, 2}; inline static int totalNumberOfSimulationScenarios = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.end(), @@ -98,11 +98,6 @@ public: SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); - void initializePatterns( - PatternGeometry &fullPattern, - PatternGeometry &reducedPatterm, - const std::vector &optimizationParameters); - static void runSimulation(const std::string &filename, std::vector &x); static std::vector> createFullPatternSimulationJobs( @@ -170,7 +165,8 @@ public: &reducedToFullInterfaceViMap, const double &normalizationFactor_translationalDisplacement, const double &normalizationFactor_rotationalDisplacement, - const double &scenarioWeight); + const double &scenarioWeight, + const ReducedPatternOptimization::Settings::ObjectiveWeights &objectiveWeights); static void constructAxialSimulationScenario( const double &forceMagnitude, const std::vector> @@ -225,6 +221,16 @@ public: static void initializeUpdateReducedPatternFunctions(); static double objective(const double &xValue); + void initializePatterns( + PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, + const std::array + &optimizationParameters); + + void setIntermediateResultsDirectoryPath( + const std::filesystem::path &newIntermediateResultsDirectoryPath); + private: static void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, @@ -236,7 +242,9 @@ private: void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); static void initializeOptimizationParameters( const std::shared_ptr &mesh, - const std::vector &optimizationParamters); + const std::array + &optimizationParamters); DRMSimulationModel simulator; void computeObjectiveValueNormalizationFactors(); @@ -256,8 +264,7 @@ private: &desiredBaseSimulationScenarioIndices, const std::filesystem::path &intermediateResultsDirectoryPath); void computeScenarioWeights(const std::vector - &baseSimulationScenarios, - const std::array &maxForceMagnitudePerBaseScenario); + &baseSimulationScenarios); std::array computeFullPatternMaxSimulationForces( const std::vector