From 535dd8a44d74bee894391787dd86a68b26bf2ae6 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Wed, 15 Dec 2021 15:21:21 +0200 Subject: [PATCH] Added ensmallen,armadillo. Parametrized optimization with choice of optimization parameters --- CMakeLists.txt | 73 ++- main.py | 130 ++++ src/main.cpp | 33 +- src/reducedmodeloptimizer.cpp | 1056 +++++++++++++-------------------- src/reducedmodeloptimizer.hpp | 43 +- 5 files changed, 612 insertions(+), 723 deletions(-) create mode 100644 main.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 253e82a..ed623b3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,15 +25,48 @@ endif() ##Create directory for the external libraries file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) + +#Add the project sources +file(GLOB SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/*.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp) + +if(MSVC) + add_compile_definitions(_HAS_STD_BYTE=0) +endif(MSVC) + +add_executable(${PROJECT_NAME} ${SOURCES}) + if(${CMAKE_BUILD_TYPE} STREQUAL "Release") - set(USE_POLYSCOPE FALSE) - set(MYSOURCES_STATIC_LINK TRUE) + set(USE_POLYSCOPE false) + set(MYSOURCES_STATIC_LINK true) + set(USE_ENSMALLEN false) else() - set(USE_POLYSCOPE TRUE) + set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE) + set(USE_POLYSCOPE true) add_compile_definitions(POLYSCOPE_DEFINED) - set(MYSOURCES_STATIC_LINK FALSE) + set(MYSOURCES_STATIC_LINK false) + set(USE_ENSMALLEN true) + if(${USE_ENSMALLEN}) + add_compile_definitions(USE_ENSMALLEN) + endif() endif() +#dlib +set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin) +file(MAKE_DIRECTORY ${DLIB_BIN_DIR}) +download_project(PROJ DLIB + GIT_REPOSITORY https://github.com/davisking/dlib.git + GIT_TAG master + BINARY_DIR ${DLIB_BIN_DIR} + PREFIX ${EXTERNAL_DEPS_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR}) +if(${MYSOURCES_STATIC_LINK}) + target_link_libraries(${PROJECT_NAME} "-static" dlib::dlib) +else() + target_link_libraries(${PROJECT_NAME} dlib::dlib) +endif() +add_compile_definitions(DLIB_DEFINED) set(MYSOURCES_SOURCE_DIR "/home/iason/Coding/Libraries/MySources") if (EXISTS ${MYSOURCES_SOURCE_DIR}) @@ -48,36 +81,16 @@ download_project(PROJ MYSOURCES endif() add_subdirectory(${MYSOURCES_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}/MySourcesBinDir) -#dlib -set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin) -file(MAKE_DIRECTORY ${DLIB_BIN_DIR}) -download_project(PROJ DLIB - GIT_REPOSITORY https://github.com/davisking/dlib.git - GIT_TAG master - BINARY_DIR ${DLIB_BIN_DIR} - PREFIX ${EXTERNAL_DEPS_DIR} - ${UPDATE_DISCONNECTED_IF_AVAILABLE} -) -add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR}) - -##Eigen 3 NOTE: Eigen is required on the system the code is ran -find_package(Eigen3 3.3 REQUIRED) - -#Add the project sources -file(GLOB SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/*.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp) - -if(MSVC) - add_compile_definitions(_HAS_STD_BYTE=0) -endif(MSVC) - -add_executable(${PROJECT_NAME} ${SOURCES}) - target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) target_include_directories(${PROJECT_NAME} PUBLIC ${MYSOURCES_SOURCE_DIR} ) if(${MYSOURCES_STATIC_LINK}) - target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen dlib::dlib MySources) + target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen MySources) else() - target_link_libraries(${PROJECT_NAME} Eigen3::Eigen dlib::dlib MySources) + target_link_libraries(${PROJECT_NAME} Eigen3::Eigen MySources) endif() + +##Eigen 3 NOTE: Eigen is required on the system the code is ran +find_package(Eigen3 3.3 REQUIRED) + diff --git a/main.py b/main.py new file mode 100644 index 0000000..b60599a --- /dev/null +++ b/main.py @@ -0,0 +1,130 @@ +# This is a sample Python script. + +# Press Shift+F10 to execute it or replace it with your code. +# Press Double Shift to search everywhere for classes, files, tool windows, actions, and settings. + + +import multiprocessing +import subprocess +import os +from os import listdir +from os.path import isfile, join +import itertools +import shutil +import sys +from datetime import datetime +from subprocess import check_output + +numberOfOptimizedPatterns=0 +numberOfSkippedPatterns=0 +numberOfFunctionCalls=0 +start_time = datetime.now() +def listener(q): + #print("Entered listener") + with open(os.path.join(resultsDir,'results.csv'), 'a') as f: + while 1: + m = q.get() + if m == 'kill': + # f.write('killed') + break + global numberOfOptimizedPatterns + numberOfOptimizedPatterns = numberOfOptimizedPatterns + 1 + print("Optimized patterns:" + str(numberOfOptimizedPatterns)) + if not m: + global numberOfSkippedPatterns + numberOfSkippedPatterns=numberOfSkippedPatterns+1 + continue + f.write(m) + f.flush() + global start_time + #print("Before") + if numberOfSkippedPatterns != numberOfOptimizedPatterns: + averageTimePerPattern_min=(datetime.now()-start_time).total_seconds()/(60*(numberOfOptimizedPatterns-numberOfSkippedPatterns)) + print("Average minutes/pattern:"+str(averageTimePerPattern_min)) + #print("After") + # global totalNumberOfPatterns + # print(totalNumberOfPatterns) + # global totalNumberOfPatterns + # print(totalNumberOfPatterns) + # completionPercentage=numberOfOptimizedPatterns/totalNumberOfPatterns + # print("Optimized patterns:" + str(completionPercentage)) + +def optimize(fullPatternFilepath, reducedPatternFilepath,translationalObjectiveWeight): + """Call run(), catch exceptions.""" + # dirname = os.path.abspath(os.path.dirname(__file__)) + dirname="/home/iason/Coding/build/ReducedModelOptimization/Release" + executableFilepath = os.path.join(dirname, 'ReducedModelOptimization') + #putResultsTo=os.path.join(resultsDir,str(translationalObjectiveWeight)) + #putResultsTo=os.path.join(putResultsTo,os.path.basename(os.path.dirname(os.path.dirname(fullPatternFilepath)))) + putResultsTo=os.path.join(resultsDir,os.path.basename(os.path.dirname(os.path.dirname(fullPatternFilepath)))) + try: + # global resultsDir + p = args = (executableFilepath, + fullPatternFilepath, + reducedPatternFilepath, + str(numberOfFunctionCalls),str(translationalObjectiveWeight), + #os.path.join(resultsDir,os.path.basename(os.path.dirname(fullPatternFilepath)))) + putResultsTo) + + patternStartTime=datetime.now() + #print("Optimizing " + fullPatternFilepath+" at "+str(datetime.now())) + popen = subprocess.Popen(args, stdout=subprocess.PIPE,bufsize=1005024) + #popen.wait() + #output = popen.stdout.read() + output = check_output(args).decode("utf-8") + #output,error=popen.communicate() + duration_min=(datetime.now() - patternStartTime).total_seconds()/60 + #print("Optimized " + fullPatternFilepath+" in "+str(duration_min)+" minutes") + q.put(output) + + except Exception as e: + print("error: %s run(*%r, **%r)" % (e, fullPatternFilepath,reducedPatternFilepath)) + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + dirOfThisFile = os.path.abspath(os.path.dirname(__file__)) + fullPatternDirectory= "/home/iason/Coding/Projects/Approximating shapes with flat patterns/ReducedModelOptimization/TestSet/FullPatterns/selectionOfPatterns" + #resultsDir=os.path.join(dirOfThisFile,os.path.join('Results/OptimizationResults/',os.path.basename(fullPatternDirectory))) + numberOfFunctionCalls=100000 + optimizationBatchName='variableComparison_allVars_joint_forConfirmation'+'_'+str(int(numberOfFunctionCalls/1000))+'k' + resultsDir=os.path.join(dirOfThisFile,os.path.join('Results/OptimizationResults/',optimizationBatchName)) + #print(resultsDir) + if not os.path.exists(resultsDir): + os.makedirs(resultsDir) + #shutil.rmtree(resultsDir) + manager = multiprocessing.Manager() + q = manager.Queue() + reducedPatternFilepath= "TestSet/ReducedPatterns/single_reduced.ply" + fullPatternFilepaths=[] + pool=multiprocessing.Pool(11) + watcher = pool.apply_async(listener, (q,)) + for subdir, dirs, files in os.walk(fullPatternDirectory): + #print(subdir) + #print(dirs) + #print(files) + #if !os.path.exists(resultsDir): + # Write results + #fullPatternFilepaths.extend([join(subdir, fullPatternFilepath) for fullPatternFilepath in listdir(subdir) if isfile(join(subdir, fullPatternFilepath))]) + fullPatternFilepaths.extend(os.path.join(subdir, fullPatternFile) for fullPatternFile in files if fullPatternFile.endswith(".ply")) + #fullPatternFilepaths.extend([os.path.join(fullPatternDirect]) + #print(fullPatternFilepaths) + + #print(optimizationPairs) + jobs=[] + translationalObjectiveWeights=[1.2] #[x/10 for x in range(2,18,2)] + jobs.extend(list(itertools.product(fullPatternFilepaths,[reducedPatternFilepath],translationalObjectiveWeights))) + #print(optimizationPairs) + totalNumberOfPatterns=len(jobs) + print("Runnning:",optimizationBatchName) + print("Number of function calls:",numberOfFunctionCalls) + print("Full pattern test set directory:"+fullPatternDirectory) + print("Total number of optimization jobs:"+ str(totalNumberOfPatterns)) + + print("Start time:", start_time) + pool.starmap(optimize,jobs) + print("Completed") + # f.close() + q.put('kill') + pool.close() + pool.join() + diff --git a/src/main.cpp b/src/main.cpp index f4a9f19..e68b54b 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -126,26 +126,27 @@ int main(int argc, char *argv[]) { auto elapsed = std::chrono::duration_cast(end - start); optimizationResults.time = elapsed.count() / 1000.0; - if (optimizationResults.wasSuccessful) { - resultsOutputDir = convergedJobsDirPath.string(); - 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_results << std::to_string(fullPattern.EN()) + "#" + pairName; - optimizationResults.writeResultsTo(csvVector); - settings_optimization.writeSettingsTo(csv_results); - csv_results << endrow; - } else { + 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; } #ifdef POLYSCOPE_DEFINED diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index c64d26b..5e2d933 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -3,7 +3,9 @@ #include "simulationhistoryplotter.hpp" #include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" +#include #include +#include #include using namespace ReducedPatternOptimization; @@ -24,6 +26,7 @@ struct GlobalOptimizationVariables std::vector objectiveValueHistory; std::vector plotColors; std::array parametersInitialValue; + std::array optimizationInitialValue; std::vector simulationScenarioIndices; double minY{DBL_MAX}; std::vector minX; @@ -46,7 +49,10 @@ struct GlobalOptimizationVariables std::shared_ptr &pReducedPatternSimulationMesh)>, 7> functions_updateReducedPatternParameter; - + std::vector xMin; + std::vector xMax; + std::vector scenarioEqualizationWeight; + std::vector scenarioUserWeights; } global; double ReducedModelOptimizer::computeDisplacementError( @@ -148,18 +154,47 @@ double ReducedModelOptimizer::computeError( // return ReducedModelOptimizer::objective(x.size(), x.data()); //} -double ReducedModelOptimizer::objective(const double &xValue) +#ifdef USE_ENSMALLEN +struct EnsmallenOptimizationObjective { - dlib::matrix x{xValue}; - return objective(x); -} + static double Evaluate(const arma::mat &x_arma) + { + std::vector x(x_arma.begin(), x_arma.end()); + for (int xi = 0; xi < x.size(); xi++) { + if (x[xi] > global.xMax[xi]) { + return std::numeric_limits::max(); + x[xi] = global.xMax[xi]; + // std::cout << "Out of range" << std::endl; + } else if (x[xi] < global.xMin[xi]) { + return std::numeric_limits::max(); + x[xi] = global.xMin[xi]; + // std::cout << "Out of range" << std::endl; + } + } + return ReducedModelOptimizer::objective(x); + } +}; +#endif + +#ifdef DLIB_DEFINED double ReducedModelOptimizer::objective(const dlib::matrix &x) { - std::cout.precision(17); - for (int i = 0; i < x.size(); i++) { - std::cout << x(i) << " "; - } - std::cout << std::endl; + return objective(std::vector(x.begin(), x.end())); +} +#endif + +double ReducedModelOptimizer::objective(const double &xValue) +{ + return objective({xValue}); +} + +double ReducedModelOptimizer::objective(const std::vector &x) +{ + // std::cout.precision(17); + // for (int i = 0; i < x.size(); i++) { + // std::cout << x[i] << " "; + // } + // std::cout << std::endl; // std::cout << x(x.size() - 2) << " " << x(x.size() - 1) << std::endl; // const Element &e = @@ -168,7 +203,7 @@ double ReducedModelOptimizer::objective(const dlib::matrix &x) // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor // << std::endl; - const int n = x.size(); + // const int n = x.size(); std::shared_ptr &pReducedPatternSimulationMesh = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; function_updateReducedPattern(x, pReducedPatternSimulationMesh); @@ -185,80 +220,92 @@ double ReducedModelOptimizer::objective(const dlib::matrix &x) // run simulations double totalError = 0; LinearSimulationModel simulator; + simulator.setStructure(pReducedPatternSimulationMesh); + // simulator.initialize(); // FormFinder simulator; - for (const int simulationScenarioIndex : global.simulationScenarioIndices) { - const std::shared_ptr &reducedJob - = global.reducedPatternSimulationJobs[simulationScenarioIndex]; - //#ifdef POLYSCOPE_DEFINED - // std::cout << reducedJob->getLabel() << ":" << std::endl; - //#endif - SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob); - // std::string filename; - if (!reducedModelResults.converged) { - totalError += std::numeric_limits::max(); - global.numOfSimulationCrashes++; + + std::for_each( + std::execution::par_unseq, + global.simulationScenarioIndices.begin(), + global.simulationScenarioIndices.end(), + [&](const int &simulationScenarioIndex) { + // for (const int simulationScenarioIndex : global.simulationScenarioIndices) { + const std::shared_ptr &reducedJob + = global.reducedPatternSimulationJobs[simulationScenarioIndex]; + //#ifdef POLYSCOPE_DEFINED + // std::cout << reducedJob->getLabel() << ":" << std::endl; + //#endif + SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob); + // std::string filename; + if (!reducedModelResults.converged) { + totalError += std::numeric_limits::max(); + global.numOfSimulationCrashes++; #ifdef POLYSCOPE_DEFINED - std::cout << "Failed simulation" << std::endl; + std::cout << "Failed simulation" << std::endl; #endif - } else { - constexpr bool usePotentialEnergy = false; - double simulationScenarioError; - if (usePotentialEnergy) { - simulationScenarioError = std::abs( - reducedModelResults.internalPotentialEnergy - - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); } else { - simulationScenarioError = computeError( + // double simulationScenarioError; + // constexpr bool usePotentialEnergy = false; + // if (usePotentialEnergy) { + // simulationScenarioError = std::abs( + // reducedModelResults.internalPotentialEnergy + // - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); + // } else { + const double simulationScenarioError = computeError( global.fullPatternResults[simulationScenarioIndex], reducedModelResults, global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex], global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + // } + //#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 + // if (global.optimizationSettings.normalizationStrategy != + // NormalizationStrategy::Epsilon && + // simulationScenarioError > 1) { + // std::cout << "Simulation scenario " + // << + // simulationScenarioStrings[simulationScenarioIndex] + // << " results in an error bigger than one." << + // std::endl; + // for (size_t parameterIndex = 0; parameterIndex < n; + // parameterIndex++) { + // std::cout << "x[" + std::to_string(parameterIndex) + "]=" + // << x[parameterIndex] << std::endl; + // } + // } + + // #ifdef POLYSCOPE_DEFINED + // ReducedModelOptimizer::visualizeResults( + // global.fullPatternSimulationJobs[simulationScenarioIndex], + // global.reducedPatternSimulationJobs[simulationScenarioIndex], + // global.reducedToFullInterfaceViMap, false); + + // ReducedModelOptimizer::visualizeResults( + // global.fullPatternSimulationJobs[simulationScenarioIndex], + // std::make_shared( + // reducedPatternMaximumDisplacementSimulationJobs + // [simulationScenarioIndex]), + // global.reducedToFullInterfaceViMap, true); + // polyscope::removeAllStructures(); + // #endif // POLYSCOPE_DEFINED + totalError += /*global.scenarioUserWeights[simulationScenarioIndex] + **/ + /*global.scenarioEqualizationWeight[simulationScenarioIndex] + **/ + simulationScenarioError * 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 - // if (global.optimizationSettings.normalizationStrategy != - // NormalizationStrategy::Epsilon && - // simulationScenarioError > 1) { - // std::cout << "Simulation scenario " - // << - // simulationScenarioStrings[simulationScenarioIndex] - // << " results in an error bigger than one." << - // std::endl; - // for (size_t parameterIndex = 0; parameterIndex < n; - // parameterIndex++) { - // std::cout << "x[" + std::to_string(parameterIndex) + "]=" - // << x[parameterIndex] << std::endl; - // } - // } - - // #ifdef POLYSCOPE_DEFINED - // ReducedModelOptimizer::visualizeResults( - // global.fullPatternSimulationJobs[simulationScenarioIndex], - // global.reducedPatternSimulationJobs[simulationScenarioIndex], - // global.reducedToFullInterfaceViMap, false); - - // ReducedModelOptimizer::visualizeResults( - // global.fullPatternSimulationJobs[simulationScenarioIndex], - // std::make_shared( - // reducedPatternMaximumDisplacementSimulationJobs - // [simulationScenarioIndex]), - // global.reducedToFullInterfaceViMap, true); - // polyscope::removeAllStructures(); - // #endif // POLYSCOPE_DEFINED - totalError += simulationScenarioError * simulationScenarioError; - } - } - // std::cout << error << std::endl; + }); + // std::cout << totalError << std::endl; // global.objectiveValueHistory.push_back(totalError); // global.plotColors.push_back(10); ++global.numberOfFunctionCalls; @@ -510,41 +557,6 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern, void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() { - // //TODO:Maybe I should split theta from R - // function_updateReducedPattern_geometry = [&](const dlib::matrix &x, - // std::shared_ptr - // &pReducedPatternSimulationMesh) { - // const int n = x.size(); - // assert(n >= 2); - - // // CoordType center_barycentric(1, 0, 0); - // // CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5); - // // CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric) - // // * x[n - 2]); - - // CoordType movableVertex_barycentric(1 - x(n - 2), x(n - 2) / 2, x(n - 2) / 2); - // CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0) - // * movableVertex_barycentric[0] - // + global.baseTriangle.cP(1) - // * movableVertex_barycentric[1] - // + global.baseTriangle.cP(2) - // * movableVertex_barycentric[2]; - // baseTriangleMovableVertexPosition - // = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, - // vcg::math::ToRad(x(n - 1))) - // * baseTriangleMovableVertexPosition; - - // for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; - // rotationCounter++) { - // pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() - // = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, - // vcg::math::ToRad(60.0 * rotationCounter)) - // * baseTriangleMovableVertexPosition; - // } - // }; - - int i = 0; - i++; global.functions_updateReducedPatternParameter[R] = [](const double &newR, std::shared_ptr &pReducedPatternSimulationMesh) { const CoordType barycentricCoordinates_hexagonBaseTriangleVertex(1 - newR, @@ -599,19 +611,13 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); } }; - // global.functions_updateReducedPatternParameter[I] = - // [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI) { - // for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - // Element &e = pReducedPatternSimulationMesh->elements[ei]; - // e.inertia.I2 = newI; - // e.inertia.I3 = newI; - // } - // }; + global.functions_updateReducedPatternParameter[I2] = [](const double &newI2, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.inertia.I2 = newI2; + e.updateRigidity(); } }; global.functions_updateReducedPatternParameter[I3] = @@ -619,6 +625,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.inertia.I3 = newI3; + e.updateRigidity(); } }; global.functions_updateReducedPatternParameter[J] = @@ -626,6 +633,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.inertia.J = newJ; + e.updateRigidity(); } }; } @@ -646,86 +654,38 @@ void ReducedModelOptimizer::initializeOptimizationParameters( case E: global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].material.youngsModulus; + global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case A: - global.parametersInitialValue[optimizationParameterIndex] - = mesh->elements[0].material.youngsModulus; + global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].A; + global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case I2: global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] .inertia.I2; + global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case I3: global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] .inertia.I3; + global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case J: global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] .inertia.J; + global.optimizationInitialValue[optimizationParameterIndex] = 1; break; case R: global.parametersInitialValue[optimizationParameterIndex] = 0; + global.optimizationInitialValue[optimizationParameterIndex] = 0.5; break; case Theta: global.parametersInitialValue[optimizationParameterIndex] = 0; + global.optimizationInitialValue[optimizationParameterIndex] = 0; break; } - - // if (xOptimizationParameter.label == "E") { - // global.updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_E); - // } else if (xOptimizationParameter.label == "A") { - // global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].A; - // global.updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_A); - // } else if (xOptimizationParameter.label == "I") { - // global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].inertia.I2; - // global.updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_I); - // } else if (xOptimizationParameter.label == "I2") { - // global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].inertia.I2; - // global.updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_I2); - // } else if (xOptimizationParameter.label == "I3") { - // global.parametersInitialValue(optimizationParameterIndex) = mesh->elements[0].inertia.I3; - // global.updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_I3); - // } else if (xOptimizationParameter.label == "J") { - // global.parametersInitialValue(optimizationParameterIndex) = mesh->elements[0].inertia.J; - // global.updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_J); - // } } - //NOTE:Assuming that I2,I3 and J are be passed after E and A because the update of E and A changes I2,I3,J. - // case "HexSize": - // updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_E); - // break; - // case "HexAngle": - // updateReducedPatternFunctions_material.push_back( - // function_updateReducedPattern_material_E); - // break; } - - // if (global.updateReducedPatternFunctions_material.size() != 0) { - // function_updateReducedPattern_material = - // [&](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // for (int optimizationParameterIndex = 0; - // optimizationParameterIndex - // < global.updateReducedPatternFunctions_material.size(); - // optimizationParameterIndex++) { - // global.updateReducedPatternFunctions_material[optimizationParameterIndex]( - // pReducedPatternSimulationMesh, - // global.parametersInitialValue(optimizationParameterIndex) - // * x(optimizationParameterIndex)); - // } - // }; - // } else { - // function_updateReducedPattern_material = - // [](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) {}; - // } } void ReducedModelOptimizer::computeReducedModelSimulationJob( @@ -758,99 +718,6 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements; } -//#if POLYSCOPE_DEFINED -//void ReducedModelOptimizer::visualizeResults( -// const std::vector> &fullPatternSimulationJobs, -// const std::vector> &reducedPatternSimulationJobs, -// const std::vector &simulationScenarios, -// const std::unordered_map -// &reducedToFullInterfaceViMap) -//{ -// DRMSimulationModel simulator; -// std::shared_ptr pFullPatternSimulationMesh = fullPatternSimulationJobs[0]->pMesh; -// pFullPatternSimulationMesh->registerForDrawing(); -// pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply"); -// reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel() -// + "_undeformed.ply"); -// double totalError = 0; -// for (const int simulationScenarioIndex : simulationScenarios) { -// const std::shared_ptr &pFullPatternSimulationJob -// = fullPatternSimulationJobs[simulationScenarioIndex]; -// pFullPatternSimulationJob->registerForDrawing(pFullPatternSimulationMesh->getLabel()); -// SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob); -// fullModelResults.registerForDrawing(); -// // fullModelResults.saveDeformedModel(); -// const std::shared_ptr &pReducedPatternSimulationJob -// = reducedPatternSimulationJobs[simulationScenarioIndex]; -// SimulationResults reducedModelResults = simulator.executeSimulation( -// pReducedPatternSimulationJob); -// double normalizationFactor = 1; -// if (global.optimizationSettings.normalizationStrategy -// != ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) { -// normalizationFactor -// = global.translationalDisplacementNormalizationValues[simulationScenarioIndex]; -// } -// reducedModelResults.saveDeformedModel(); -// fullModelResults.saveDeformedModel(); -// double error = computeDisplacementError(reducedModelResults.displacements, -// fullModelResults.displacements, -// reducedToFullInterfaceViMap, -// normalizationFactor); -// std::cout << "Error of simulation scenario " -// < baseSimulationScenarioNames[simulationScenarioIndex] << " is " << error << std::endl; -// totalError += error; -// reducedModelResults.registerForDrawing(); -// // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); -// // registerWorldAxes(); -// const std::string screenshotFilename -// = "/home/iason/Coding/Projects/Approximating shapes with flat " -// "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" -// "Images/" -// + pFullPatternSimulationMesh->getLabel() + "_" -// + baseSimulationScenarioNames[simulationScenarioIndex]; -// polyscope::show(); -// polyscope::screenshot(screenshotFilename, false); -// fullModelResults.unregister(); -// reducedModelResults.unregister(); -// // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); -// } -// std::cout << "Total error:" << totalError << std::endl; -//} - -//void ReducedModelOptimizer::registerResultsForDrawing( -// const std::shared_ptr &pFullPatternSimulationJob, -// const std::shared_ptr &pReducedPatternSimulationJob, -// const std::unordered_map -// &reducedToFullInterfaceViMap) { -// DRMSimulationModel simulator; -// std::shared_ptr pFullPatternSimulationMesh = -// pFullPatternSimulationJob->pMesh; -// pFullPatternSimulationMesh->registerForDrawing(); -// // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() -// // + -// // "_undeformed.ply"); -// // reducedPatternSimulationJobs[0]->pMesh->savePly( -// // reducedPatternSimulationJobs[0]->pMesh->getLabel() + -// // "_undeformed.ply"); -// pFullPatternSimulationJob->registerForDrawing( -// pFullPatternSimulationMesh->getLabel()); -// SimulationResults fullModelResults = -// simulator.executeSimulation(pFullPatternSimulationJob); -// fullModelResults.registerForDrawing(); -// // fullModelResults.saveDeformedModel(); -// SimulationResults reducedModelResults = -// simulator.executeSimulation(pReducedPatternSimulationJob); -// // reducedModelResults.saveDeformedModel(); -// // fullModelResults.saveDeformedModel(); -// double error = computeRawDisplacementError( -// reducedModelResults.displacements, fullModelResults.displacements, -// reducedToFullInterfaceViMap/*, -// global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); -// std::cout << "Error is " << error << std::endl; -// reducedModelResults.registerForDrawing(); -//} -//#endif // POLYSCOPE_DEFINED - void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, @@ -871,35 +738,46 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( } } -void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimizationResult_dlib, +void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjective, const Settings &settings, ReducedPatternOptimization::Results &results) { //Number of crashes // results.numberOfSimulationCrashes = global.numOfSimulationCrashes; //Value of optimal objective Y - results.objectiveValue.total = optimizationResult_dlib.y; + results.objectiveValue.total = optimalObjective.y; + // results.objectiveValue.total = 0; + if (optimalObjective.y != global.minY) { + std::cout << "Different optimal y:" << optimalObjective.y << " " << global.minY + << std::endl; + } //Optimal X values results.optimalXNameValuePairs.resize(settings.parameterRanges.size()); std::vector optimalX(settings.parameterRanges.size()); - for (int xVariableIndex = 0; xVariableIndex < settings.parameterRanges.size(); - xVariableIndex++) { - if (xVariableIndex < settings.parameterRanges.size() - 2) { - results.optimalXNameValuePairs[xVariableIndex] - = std::make_pair(settings.parameterRanges[xVariableIndex].label, - optimizationResult_dlib.x(xVariableIndex) - * global.parametersInitialValue[xVariableIndex]); + for (int optimizationParameterIndex = E; + optimizationParameterIndex != NumberOfOptimizationParameters; + optimizationParameterIndex++) { + if (settings.parameterRanges[optimizationParameterIndex].max + - settings.parameterRanges[optimizationParameterIndex].min + < 1e-10) { + continue; + } + if (optimizationParameterIndex < settings.parameterRanges.size() - 2) { + results.optimalXNameValuePairs[optimizationParameterIndex] + = std::make_pair(settings.parameterRanges[optimizationParameterIndex].label, + optimalObjective.x[optimizationParameterIndex] + * global.parametersInitialValue[optimizationParameterIndex]); } else { //Hex size and angle are pure values (not multipliers of the initial values) - results.optimalXNameValuePairs[xVariableIndex] - = std::make_pair(settings.parameterRanges[xVariableIndex].label, - optimizationResult_dlib.x(xVariableIndex)); + results.optimalXNameValuePairs[optimizationParameterIndex] + = std::make_pair(settings.parameterRanges[optimizationParameterIndex].label, + optimalObjective.x[optimizationParameterIndex]); } - optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex); + optimalX[optimizationParameterIndex] = optimalObjective.x[optimizationParameterIndex]; #ifdef POLYSCOPE_DEFINED - std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":" - << optimalX[xVariableIndex] << " "; + std::cout << results.optimalXNameValuePairs[optimizationParameterIndex].first << ":" + << optimalX[optimizationParameterIndex] << " "; #endif } #ifdef POLYSCOPE_DEFINED @@ -911,7 +789,9 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza // updateMeshFunction(optimalX); std::shared_ptr &pReducedPatternSimulationMesh = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; - function_updateReducedPattern(optimizationResult_dlib.x, pReducedPatternSimulationMesh); + function_updateReducedPattern(std::vector(optimalObjective.x.begin(), + optimalObjective.x.end()), + pReducedPatternSimulationMesh); // results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios); //TODO:use push_back it will make the code more readable LinearSimulationModel simulator; @@ -943,6 +823,7 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza global.translationalDisplacementNormalizationValues[simulationScenarioIndex], global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + // results.objectiveValue.total += results.objectiveValue.perSimulationScenario_total[i]; //Raw translational const double rawTranslationalError = computeRawTranslationalError( global.fullPatternResults[simulationScenarioIndex].displacements, @@ -1102,13 +983,18 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, #if POLYSCOPE_DEFINED // global.plotColors.reserve(settings.numberOfFunctionCalls); #endif +#ifdef USE_ENSMALLEN +#else double (*objF)(const dlib::matrix &) = &objective; +#endif enum OptimizationParameterComparisonScenarioIndex { AllVar, + YMAndGeo, noYM, SplittedMatGeo, SplittedYM_before, SplittedYM_after, + SplittedGeoYMMat, NumberOfScenarios }; const std::vector>> scenarioParameters = @@ -1116,31 +1002,24 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, std::vector>> scenarioParameters( NumberOfScenarios); scenarioParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}}; + scenarioParameters[YMAndGeo] = {{R, Theta, E}}; scenarioParameters[noYM] = {{A, I2, I3, J, R, Theta}}; scenarioParameters[SplittedMatGeo] = {{E, A, I2, I3, J}, {R, Theta}}; scenarioParameters[SplittedYM_before] = {{E}, {A, I2, I3, J, R, Theta}}; scenarioParameters[SplittedYM_after] = {{A, I2, I3, J, R, Theta}, {E}}; + scenarioParameters[SplittedGeoYMMat] = {{R, Theta}, {E}, {A, I2, I3, J}}; return scenarioParameters; }(); - constexpr OptimizationParameterComparisonScenarioIndex scenario = AllVar; + constexpr OptimizationParameterComparisonScenarioIndex scenario = noYM; const std::vector> scenarioParameterGroups = scenarioParameters[scenario]; //TODO:Ensure that the reduced pattern mesh has the initial parameter values before starting the optimization - dlib::function_evaluation optimalResult; + FunctionEvaluation optimalResult; for (const std::vector ¶meterGroup : scenarioParameterGroups) { - //Set min max values - dlib::matrix xMin(parameterGroup.size()); - dlib::matrix xMax(parameterGroup.size()); - 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; - } //Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup - function_updateReducedPattern = [&](const dlib::matrix &x, + function_updateReducedPattern = [&](const std::vector &x, std::shared_ptr &pMesh) { for (int xIndex = 0; xIndex < x.size(); xIndex++) { const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; @@ -1150,377 +1029,131 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, if (parameterIndex == R || parameterIndex == Theta) { //NOTE:Here I explore the geometry parameters linearly - return x(xIndex) + parameterInitialValue; + return x[xIndex] /*+ parameterInitialValue*/; } //and the material parameters exponentially(?).TODO: Check what happens if I make all linear - return x(xIndex) * parameterInitialValue; + return x[xIndex] * parameterInitialValue; }(); global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, pMesh); } + pMesh->reset(); //NOTE: I could put this code into each updateParameter function for avoiding unessecary calculations }; - optimalResult = dlib::find_min_global(objF, - xMin, - xMax, - dlib::max_function_calls( - settings.numberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), - settings.solverAccuracy); + std::vector xMin; + std::vector xMax; + xMin.resize(parameterGroup.size()); + xMax.resize(parameterGroup.size()); + 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; + } +#ifdef USE_ENSMALLEN + arma::mat x(parameterGroup.size(), 1); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + + x(xIndex, 0) = global.optimizationInitialValue[parameterIndex]; + } + + global.xMin = xMin; + global.xMax = xMax; + + // Create simulated annealing optimizer with default options. + // The ens::SA<> type can be replaced with any suitable ensmallen optimizer + // 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; + // arma::mat xMin_arma(global.xMin); + // arma::mat xMax_arma(global.xMax); + // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); + ens::LBestPSO optimizer; + const double minima = optimizer.Optimize(optimizationFunction, x); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + if (x[xIndex] > global.xMax[xIndex]) { + x[xIndex] = global.xMax[xIndex]; + } else if (x[xIndex] < global.xMin[xIndex]) { + x[xIndex] = global.xMin[xIndex]; + } + } + +#ifdef POLYSCOPE_DEFINED + std::cout << "optimal x:" + << "\n" + << x << std::endl; + std::cout << "Minima:" << minima << std::endl; +#endif + optimalResult.x.clear(); + optimalResult.x.resize(parameterGroup.size()); + std::copy(x.begin(), x.end(), optimalResult.x.begin()); + optimalResult.y = minima; + +#else + //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.parameterRanges[parameterIndex].min; + xMax_dlib(xIndex) = settings.parameterRanges[parameterIndex].max; + } + const dlib::function_evaluation optimalResult_dlib + = dlib::find_min_global(objF, + xMin_dlib, + xMax_dlib, + dlib::max_function_calls(settings.numberOfFunctionCalls), + std::chrono::hours(24 * 365 * 290), + settings.solverAccuracy); + + // constexpr bool useBOBYQA = false; + // if (useBOBYQA) { + // const size_t npt = 2 * global.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) = global.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); + // } + optimalResult.x.clear(); + optimalResult.x.resize(parameterGroup.size()); + std::copy(optimalResult_dlib.x.begin(), optimalResult_dlib.x.end(), optimalResult.x.begin()); + optimalResult.y = optimalResult_dlib.y; + +#endif function_updateReducedPattern( - optimalResult.x, + std::vector(optimalResult.x.begin(), optimalResult.x.end()), global.reducedPatternSimulationJobs[0] ->pMesh); //TODO: Check if its ok to update only the changed parameters + // std::cout << "GLOBAL MIN:" << global.minY << std::endl; + // std::cout << "opt res y:" << optimalResult.y << std::endl; } - - // optimalResult.x.set_size(global.numberOfOptimizationParameters); - // const auto hexAngleParameterIt = std::find_if(settings.parameterRanges.begin(), - // settings.parameterRanges.end(), - // [](const xRange &x) { - // return x.label == ParameterLabels::theta; - // }); - // const bool hasHexAngleParameter = hexAngleParameterIt != settings.parameterRanges.end(); - // const auto hexSizeParameterIt = std::find_if(settings.parameterRanges.begin(), - // settings.parameterRanges.end(), - // [](const xRange &x) { - // return x.label == ParameterLabels::R; - // }); - // const bool hasHexSizeParameter = hexSizeParameterIt != settings.parameterRanges.end(); - // const bool hasBothGeometricalParameters = hasHexAngleParameter && hasHexSizeParameter; - // assert(hasBothGeometricalParameters || (!hasHexAngleParameter && !hasHexSizeParameter)); - // const int numberOfGeometryOptimizationParameters = hasBothGeometricalParameters ? 2 : 0; - // const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters - // - numberOfGeometryOptimizationParameters; - // const bool hasGeometryAndMaterialParameters = hasBothGeometricalParameters - // && numberOfMaterialOptimizationParameters != 0; - // constexpr bool splitYMOptimization = false; - // const auto parameterIt_E = std::find_if(settings.parameterRanges.begin(), - // settings.parameterRanges.end(), - // [](const xRange &x) { - // return x.label == ParameterLabels::E; - // }); - // //E exists as a parameter - // if (splitYMOptimization && parameterIt_E == settings.parameterRanges.end()) { - // std::cerr << "YM is not a parameter." << std::endl; - // std::terminate(); - // } - - // std::unordered_map parameterLabelGlobalIndexMap; - // { - // for (int globalParameterIndex = 0; - // globalParameterIndex < global.numberOfOptimizationParameters; - // globalParameterIndex++) { - // parameterLabelGlobalIndexMap[settings.parameterRanges[globalParameterIndex].label] - // = globalParameterIndex; - // } - // } - // const int parameterIndexGlobal_E = parameterLabelGlobalIndexMap.contains(ParameterLabels::E) - // ? parameterLabelGlobalIndexMap[ParameterLabels::E] - // : -1; - // const int parameterIndexGlobal_A = parameterLabelGlobalIndexMap.contains(ParameterLabels::A) - // ? parameterLabelGlobalIndexMap[ParameterLabels::A] - // : -1; - // const int parameterIndexGlobal_I2 = parameterLabelGlobalIndexMap.contains(ParameterLabels::I2) - // ? parameterLabelGlobalIndexMap[ParameterLabels::I2] - // : -1; - // const int parameterIndexGlobal_I3 = parameterLabelGlobalIndexMap.contains(ParameterLabels::I3) - // ? parameterLabelGlobalIndexMap[ParameterLabels::I3] - // : -1; - // const int parameterIndexGlobal_J = parameterLabelGlobalIndexMap.contains(ParameterLabels::J) - // ? parameterLabelGlobalIndexMap[ParameterLabels::J] - // : -1; - // const int parameterIndexGlobal_R = parameterLabelGlobalIndexMap.contains(ParameterLabels::R) - // ? parameterLabelGlobalIndexMap[ParameterLabels::R] - // : -1; - // const int parameterIndexGlobal_theta = parameterLabelGlobalIndexMap.contains( - // ParameterLabels::theta) - // ? parameterLabelGlobalIndexMap[ParameterLabels::theta] - // : -1; - - // if (splitYMOptimization) { - // std::unordered_map parameterLabelLocalIndexMap; - // //Optimize all variables besides the youngs modulus - // const int parameterIndex_E = std::distance(settings.parameterRanges.begin(), parameterIt_E); - // dlib::matrix xMin(global.numberOfOptimizationParameters - 1); - // dlib::matrix xMax(global.numberOfOptimizationParameters - 1); - // { - // int xIndex = 0; - // for (int globalParameterIndex = 0; - // globalParameterIndex < global.numberOfOptimizationParameters; - // globalParameterIndex++) { - // if (globalParameterIndex == parameterIndexGlobal_E) { - // continue; - // } - // parameterLabelLocalIndexMap[settings.parameterRanges[globalParameterIndex].label] - // = xIndex; - // xMin(xIndex) = settings.parameterRanges[globalParameterIndex].min; - // xMax(xIndex) = settings.parameterRanges[globalParameterIndex].max; - // xIndex++; - // } - // } - // const int parameterLocalIndex_A = parameterLabelLocalIndexMap[ParameterLabels::A]; - // const int parameterLocalIndex_I2 = parameterLabelLocalIndexMap[ParameterLabels::I2]; - // const int parameterLocalIndex_I3 = parameterLabelLocalIndexMap[ParameterLabels::I3]; - // const int parameterLocalIndex_J = parameterLabelLocalIndexMap[ParameterLabels::J]; - // const int parameterLocalIndex_R = parameterLabelLocalIndexMap[ParameterLabels::R]; - // const int parameterLocalIndex_theta = parameterLabelLocalIndexMap[ParameterLabels::theta]; - // constexpr std::array<7, bool> optimizeAllBesidesE(false); - - // function_updateReducedPattern = - // [](const dlib::matrix &x, - // const std::vector &initialValues, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_material_A(pReducedPatternSimulationMesh, - // initialValues[parameterIndexGlobal_A] - // * x(parameterLocalIndex_A)); - // function_updateReducedPattern_material_I2(pReducedPatternSimulationMesh, - // initialValues[parameterIndexGlobal_I2] - // * x(parameterLocalIndex_I2)); - // function_updateReducedPattern_material_I3(pReducedPatternSimulationMesh, - // initialValues[parameterIndexGlobal_I3] - // * x(parameterLocalIndex_I3)); - // function_updateReducedPattern_material_J(pReducedPatternSimulationMesh, - // initialValues[parameterIndexGlobal_J] - // * x(parameterLocalIndex_J)); - // function_updateReducedPattern_geometry({x(parameterLocalIndex_R), - // x(parameterLocalIndex_theta)}, - // pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // }; - - // const dlib::function_evaluation result_dlib_noYM - // = dlib::find_min_global(objF, - // xMin, - // xMax, - // dlib::max_function_calls(settings.numberOfFunctionCalls), - // std::chrono::hours(24 * 365 * 290), - // settings.solverAccuracy); - // function_updateReducedPattern(result_dlib_noYM.x, - // global.reducedPatternSimulationJobs[0]->pMesh); - // //Optimize for the youngs modulus only - // dlib::matrix xEMin(1); - // dlib::matrix xEMax(1); - // xEMin(0) = parameterIt_E->min; - // xEMax(0) = parameterIt_E->max; - - // function_updateReducedPattern = [&](const dlib::matrix &x, - // const std::vector &initialValues, - // std::shared_ptr - // &pReducedPatternSimulationMesh) { - // // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - // function_updateReducedPattern_material_E(pReducedPatternSimulationMesh, - // x(0) * initialValues[parameterIndexGlobal_E]); - // function_updateReducedPattern_geometry( - // {result_dlib_noYM.x(parameterLocalIndex_R), - // result_dlib_noYM.x(parameterLocalIndex_theta)}, - // pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this - // pReducedPatternSimulationMesh->reset(); - // }; - // double optimalE = 1; - // double (*objective_singleVariable)(const double &) = &objective; - // const double optimalObjectiveValue = dlib::find_min_single_variable(objective_singleVariable, - // optimalE, - // parameterIt_E->min, - // parameterIt_E->max, - // 1e-2, - // 100); - // //Gather results of the two optimizations - // for (int parameterIndex = 0; parameterIndex < numberOfMaterialOptimizationParameters; - // parameterIndex++) { - // if (parameterIndex == parameterIndex_E) { - // optimalResult.x(parameterIndex) = optimalE; - // continue; - // } - // assert(parameterLabelLocalIndexMap.contains( - // settings.parameterRanges[parameterIndex].label)); - // const int firstOptimizationParameterIndex - // = parameterLabelLocalIndexMap[settings.parameterRanges[parameterIndex].label]; - // optimalResult.x(parameterIndex) = result_dlib_noYM.x(firstOptimizationParameterIndex); - // } - // function_updateReducedPattern = - // [&](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // }; - - // } else if (settings.splitGeometryMaterialOptimization && hasGeometryAndMaterialParameters) { - // //Geometry optimization of the reduced pattern - // dlib::matrix xGeometryMin(numberOfGeometryOptimizationParameters); - // dlib::matrix xGeometryMax(numberOfGeometryOptimizationParameters); - // // const int hexAngleParameterIndex = std::distance(settings.xRanges.begin(), - // // hexAngleParameterIt); - // // const int hexSizeParameterIndex = std::distance(settings.xRanges.begin(), - // // hexSizeParameterIt); - // xGeometryMin(0) = hexSizeParameterIt->min; - // xGeometryMax(0) = hexSizeParameterIt->max; - // xGeometryMin(1) = hexAngleParameterIt->min; - // xGeometryMax(1) = hexAngleParameterIt->max; - - // //Set reduced pattern update functions to be used during optimization - // function_updateReducedPattern = - // [&](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // }; - // const int optimizationFunctionCalls_geometry - // = global.numberOfOptimizationParameters == numberOfGeometryOptimizationParameters - // ? settings.numberOfFunctionCalls - // : static_cast(settings.numberOfFunctionCalls / 3.0); - // dlib::function_evaluation result_dlib_geometry - // = dlib::find_min_global(objF, - // xGeometryMin, - // xGeometryMax, - // dlib::max_function_calls(optimizationFunctionCalls_geometry), - // std::chrono::hours(24 * 365 * 290), - // settings.solverAccuracy); - // function_updateReducedPattern(result_dlib_geometry.x, - // global.reducedPatternSimulationJobs[0]->pMesh); - // //Material optimization of the reduced pattern - // // std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl; - // // std::cout << "opt size:" << global.minX[0] << std::endl; - // dlib::function_evaluation result_dlib_material; - // if (numberOfMaterialOptimizationParameters != 0) { - // dlib::matrix xMaterialMin(numberOfMaterialOptimizationParameters); - // dlib::matrix xMaterialMax(numberOfMaterialOptimizationParameters); - // for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) { - // xMaterialMin(i) = settings.parameterRanges[i].min; - // xMaterialMax(i) = settings.parameterRanges[i].max; - // } - // function_updateReducedPattern = - // [&](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - // function_updateReducedPattern_geometry( - // result_dlib_geometry.x, - // pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this - // pReducedPatternSimulationMesh->reset(); - // }; - // result_dlib_material = dlib::find_min_global(objF, - // xMaterialMin, - // xMaterialMax, - // dlib::max_function_calls(static_cast( - // 2.0 * settings.numberOfFunctionCalls - // / 3.0)), - // std::chrono::hours(24 * 365 * 290), - // settings.solverAccuracy); - // } - // // constexpr bool useBOBYQA = false; - // // if (useBOBYQA) { - // // const size_t npt = 2 * global.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 = 200; - // // dlib::find_min_bobyqa(objF, - // // result_dlib.x, - // // npt, - // // xMaterialMin, - // // xMaterialMax, - // // rhobeg, - // // rhoend, - // // bobyqa_maxFunctionCalls); - // // } - - // std::copy(result_dlib_material.x.begin(), - // result_dlib_material.x.end(), - // optimalResult.x.begin()); - // // debug_x.begin()); - // std::copy(result_dlib_geometry.x.begin(), - // result_dlib_geometry.x.end(), - // optimalResult.x.begin() + numberOfMaterialOptimizationParameters); - // // std::cout << "opt x:"; - // // for (const auto optx : optimalResult.x) { - // // std::cout << optx << " "; - // // } - // // std::cout << std::endl; - // function_updateReducedPattern = - // [&](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // }; - // // const auto meshLabel = global.reducedPatternSimulationJobs[0]->pMesh->getLabel(); - // // global.reducedPatternSimulationJobs[0]->pMesh->setLabel("pre"); - // // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial); - // // polyscope::show(); - // // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); - // } else { - // dlib::matrix xMin(global.numberOfOptimizationParameters); - // dlib::matrix xMax(global.numberOfOptimizationParameters); - // for (int i = 0; i < global.numberOfOptimizationParameters; i++) { - // xMin(i) = settings.parameterRanges[i].min; - // xMax(i) = settings.parameterRanges[i].max; - // } - // if (numberOfGeometryOptimizationParameters != 0 - // && numberOfMaterialOptimizationParameters != 0) { - // function_updateReducedPattern = [&](const dlib::matrix &x, - // std::shared_ptr - // &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // // function_updateReducedPattern_material_E(pReducedPatternSimulationMesh, - // // global.initialParameters( - // // parameterIndexGlobal_E) - // // * x(parameterIndexGlobal_E)); - // // function_updateReducedPattern_material_A(pReducedPatternSimulationMesh, - // // global.initialParameters( - // // parameterIndexGlobal_A) - // // * x(parameterIndexGlobal_A)); - // // function_updateReducedPattern_material_I2(pReducedPatternSimulationMesh, - // // global.initialParameters( - // // parameterIndexGlobal_I2) - // // * x(parameterIndexGlobal_I2)); - // // function_updateReducedPattern_material_I3(pReducedPatternSimulationMesh, - // // global.initialParameters( - // // parameterIndexGlobal_I3) - // // * x(parameterIndexGlobal_I3)); - // // function_updateReducedPattern_material_J(pReducedPatternSimulationMesh, - // // global.initialParameters( - // // parameterIndexGlobal_J) - // // * x(parameterIndexGlobal_J)); - // // function_updateReducedPattern_geometry({x(parameterIndexGlobal_R), - // // x(parameterIndexGlobal_theta)}, - // // pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // }; - // } /* else if (numberOfGeometryOptimizationParameters == 0) { - // function_updateReducedPattern = - // [&](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // }; - - // } else { - // function_updateReducedPattern = - // [&](const dlib::matrix &x, - // std::shared_ptr &pReducedPatternSimulationMesh) { - // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - // pReducedPatternSimulationMesh->reset(); - // }; - // }*/ - // optimalResult = dlib::find_min_global(objF, - // xMin, - // xMax, - // dlib::max_function_calls( - // settings.numberOfFunctionCalls), - // std::chrono::hours(24 * 365 * 290), - // settings.solverAccuracy); - // } - // function_updateReducedPattern(optimalResult.x, global.reducedPatternSimulationJobs[0]->pMesh); - // const auto debug_currObj = objective(optimalResult.x); - // if (std::abs(debug_currObj - optimalResult.y) > 1e-5) { - // std::cout << "wrong optimal result.Previous was:" << optimalResult.y - // << " .Current is:" << debug_currObj << std::endl; - // } - // optimalResult.y = debug_currObj; getResults(optimalResult, settings, results); } @@ -1720,6 +1353,34 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni return error; } +void search(const std::function objectiveFunction, + const double &targetY, + double &optimalX, + double &error, + const double &epsilon, + const double &maxIterations) +{ + error = std::numeric_limits::max(); + int iterationIndex = 0; + double low = optimalX / 1e-3, high = optimalX * 1e3; + while (error > epsilon && iterationIndex < maxIterations) { + const double y = objectiveFunction(optimalX); + error = std::abs(targetY - y); + if (y > targetY) { + high = optimalX; + } else { + low = optimalX; + } + optimalX = low + (high - low) / 2; + + iterationIndex++; + } + + if (iterationIndex == maxIterations) { + std::cerr << "Max iterations reached." << std::endl; + } +} + double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude) { SimulationJob job; @@ -1772,6 +1433,15 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa return error; } +#ifdef USE_ENSMALLEN +struct ForceMagnitudeOptimization +{ + std::function objectiveFunction; + ForceMagnitudeOptimization(std::function &f) : objectiveFunction(f) {} + double Evaluate(const arma::mat &x) { return objectiveFunction(x(0, 0)); } +}; +#endif + double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( const BaseSimulationScenario &scenario) { @@ -1813,12 +1483,34 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( break; } + constexpr int maxIterations = 1000; minimumError = dlib::find_min_single_variable(objectiveFunction, forceMagnitude, 1e-4, 1e4, forceMagnitudeEpsilon, - 1000); + 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 << ReducedPatternOptimization::baseSimulationScenarioNames[scenario] << ": " + // << optimalObjective << std::endl; + // forceMagnitude = forceMagnitude_mat(0, 0); + +// search(objectiveFunction, +// global.desiredMaxDisplacementValue, +// forceMagnitude, +// minimumError, +// objectiveEpsilon, +// maxIterations); +#endif + wasSuccessful = minimumError < objectiveEpsilon; #ifdef POLYSCOPE_DEFINED @@ -1855,6 +1547,12 @@ std::vector> ReducedModelOptimizer::createFullPat SimulationJob job; job.pMesh = pMesh; + global.scenarioEqualizationWeight.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 + global.scenarioUserWeights.resize(totalNumberOfSimulationScenarios); + const double userWeightsSum = std::accumulate(baseScenarioWeights.begin(), + baseScenarioWeights.end(), + 0); for (std::pair scenarioMaxForceMagnitudePair : scenarioMaxForceMagnitudePairs) { const BaseSimulationScenario scenario = scenarioMaxForceMagnitudePair.first; @@ -1871,6 +1569,11 @@ std::vector> ReducedModelOptimizer::createFullPat = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.begin() + scenario, 0); + const double baseScenarioEqualizationWeight = static_cast( + numberOfSimulationScenarios) + / totalNumberOfSimulationScenarios; + const double baseScenarioUserWeight = static_cast(baseScenarioWeights[scenario]) + / userWeightsSum; for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationScenarios; simulationScenarioIndex++) { job.nodalExternalForces.clear(); @@ -1887,6 +1590,9 @@ std::vector> ReducedModelOptimizer::createFullPat scenarios[baseSimulationScenarioIndexOffset + simulationScenarioIndex] = std::make_shared(job); + global.scenarioEqualizationWeight[simulationScenarioIndex] + = baseScenarioEqualizationWeight; + global.scenarioUserWeights[simulationScenarioIndex] = baseScenarioUserWeight; } } @@ -1898,6 +1604,8 @@ std::vector> ReducedModelOptimizer::createFullPat void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { + // m_pFullPatternSimulationMesh->registerForDrawing(); + // m_pFullPatternSimulationMesh->save(std::filesystem::current_path().append("initial.ply")); // Compute the sum of the displacement norms std::vector fullPatternTranslationalDisplacementNormSum( totalNumberOfSimulationScenarios); @@ -1920,8 +1628,16 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() const Vector6d &vertexDisplacement = global.fullPatternResults[simulationScenarioIndex] .displacements[fullPatternVi]; + // std::cout << "vi:" << fullPatternVi << std::endl; + // std::cout << "displacement:" << vertexDisplacement.getTranslation().norm() << std::endl; + translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm(); } + // global.fullPatternResults[simulationScenarioIndex].saveDeformedModel( + // std::filesystem::current_path()); + // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(); + // polyscope::show(); + // global.fullPatternResults[simulationScenarioIndex].unregister(); double angularDistanceSum = 0; for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { const int fullPatternVi = interfaceViPair.second; @@ -2029,11 +1745,33 @@ void ReducedModelOptimizer::optimize( const std::shared_ptr &pFullPatternSimulationJob = global.fullPatternSimulationJobs[simulationScenarioIndex]; - // SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, - // simulationSettings); + // std::filesystem::path patternMaxForceMagnitudesFilePath( + // std::filesystem::path(forceMagnitudesDirectoryPath) + // .append(pFullPatternSimulationJob->getLabel() + ".json")); + // const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists( + // patternMaxForceMagnitudesFilePath); + // if (fullPatternScenarioMagnitudesExist) { + // nlohmann::json json; + // std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); + // ifs >> json; + // fullPatternSimulationScenarioMaxMagnitudes + // = static_cast>>( + // json.at("maxMagn")); + // const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size() + // != desiredBaseSimulationScenarioIndices.size(); + // if (!shouldRecompute) { + // return fullPatternSimulationScenarioMaxMagnitudes; + // } + // } +#ifdef POLYSCOPE_DEFINED LinearSimulationModel linearSimulator; SimulationResults fullPatternResults = linearSimulator.executeSimulation( pFullPatternSimulationJob); + +#else + SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, + simulationSettings); +#endif // if (!fullPatternResults.converged) { // DRMSimulationModel::Settings simulationSettings_secondRound; // simulationSettings_secondRound.viscousDampingFactor = 2e-3; @@ -2048,9 +1786,9 @@ void ReducedModelOptimizer::optimize( if (!fullPatternResults.converged) { results.wasSuccessful = false; -#ifdef POLYSCOPE_DEFINED - std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel() + std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel() << " did not converge." << std::endl; +#ifdef POLYSCOPE_DEFINED // DRMSimulationModel::Settings debugSimulationSettings; // debugSimulationSettings.debugModeStep = 50; // // // debugSimulationSettings.maxDRMIterations = 100000; diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 8bcd6a2..638ef1e 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -9,8 +9,10 @@ #include "reducedmodeloptimizer_structs.hpp" #include "simulationmesh.hpp" #include +#ifdef DLIB_DEFINED #include #include +#endif #ifdef POLYSCOPE_DEFINED #include "polyscope/color_management.h" @@ -39,7 +41,16 @@ class ReducedModelOptimizer constexpr static double youngsModulus{1 * 1e9}; public: - constexpr static struct ParameterLabels + struct FunctionEvaluation + { + FunctionEvaluation() = default; + FunctionEvaluation(const std::vector &x, double y) : x(x), y(y) {} + + std::vector x; + double y = std::numeric_limits::quiet_NaN(); + }; + + struct ParameterLabels { inline const static std::string E = {"E"}; inline const static std::string A = {"A"}; @@ -48,11 +59,13 @@ public: inline const static std::string J = {"J"}; inline const static std::string theta = {"Theta"}; inline const static std::string R = {"R"}; - } parameterLabels; + }; + inline constexpr static ParameterLabels parameterLabels(); enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationParameters }; constexpr static std::array simulationScenariosResolution = {11, 11, 20, 20, 20}; + constexpr static std::array baseScenarioWeights = {1, 1, 5, 5, 5}; // constexpr static std::array simulationScenariosResolution = {3, 3, 3, 3, 3}; inline static int totalNumberOfSimulationScenarios = std::accumulate(simulationScenariosResolution.begin(), @@ -183,12 +196,9 @@ public: const std::vector> &oppositeInterfaceViPairs, SimulationJob &job); - static std::function &x, + static std::function &x, std::shared_ptr &pReducedPatternSimulationMesh)> function_updateReducedPattern; - static std::function &x, - std::shared_ptr &pReducedPatternSimulationMesh)> - function_updateReducedPattern_material; static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_E; @@ -207,9 +217,9 @@ public: static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_J; - static std::function &x, - std::shared_ptr &pReducedPatternSimulationMesh)> - function_updateReducedPattern_geometry; + static double objective(const std::vector &x); + static void initializeUpdateReducedPatternFunctions(); + static double objective(const double &xValue); private: static void computeDesiredReducedModelDisplacements( @@ -230,7 +240,8 @@ private: DRMSimulationModel simulator; void computeObjectiveValueNormalizationFactors(); - static void getResults(const dlib::function_evaluation &optimizationResult_dlib, + + static void getResults(const FunctionEvaluation &optimalObjective, const ReducedPatternOptimization::Settings &settings, ReducedPatternOptimization::Results &results); double computeFullPatternMaxSimulationForce( @@ -244,12 +255,11 @@ private: getFullPatternMaxSimulationForces( const std::vector &desiredBaseSimulationScenarioIndices); + +#ifdef DLIB_DEFINED static double objective(const dlib::matrix &x); - static void initializeUpdateReducedPatternFunctions(); - static double objective(const double &xValue); +#endif }; -inline std::function &x, std::shared_ptr &m)> - ReducedModelOptimizer::function_updateReducedPattern_material; inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_E; @@ -268,10 +278,7 @@ inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_J; -inline std::function &x, std::shared_ptr &m)> - ReducedModelOptimizer::function_updateReducedPattern_geometry; -inline std::function &x, - std::shared_ptr &m)> +inline std::function &x, std::shared_ptr &m)> ReducedModelOptimizer::function_updateReducedPattern; #endif // REDUCEDMODELOPTIMIZER_HPP