diff --git a/CMakeLists.txt b/CMakeLists.txt index 118d6bc..fa0d48a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -44,6 +44,17 @@ add_compile_definitions(BOBYQA_DEBUG) add_subdirectory(${BOBYQA_SOURCE_DIR}) message(STATUS "BOBYQA bin dir:${BOBYQA_BINARY_DIR}") +#dlib +download_project(PROJ DLIB + GIT_REPOSITORY https://github.com/davisking/dlib.git + GIT_TAG master + PREFIX ${CMAKE_CURRENT_SOURCE_DIR}/build/external/ + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${DLIB_SOURCE_DIR}) +#message(STATUS "BOBYQA bin dir:${BOBYQA_BINARY_DIR}") + + ##vcglib devel branch download_project(PROJ vcglib_devel GIT_REPOSITORY https://github.com/cnr-isti-vclab/vcglib.git @@ -78,4 +89,4 @@ else(MSVC) target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) endif(MSVC) -target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen OpenMP::OpenMP_CXX matplot bobyqa_shared) +target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen OpenMP::OpenMP_CXX matplot bobyqa_shared dlib::dlib) diff --git a/src/main.cpp b/src/main.cpp index 2f8cf23..ec5c1ee 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -45,7 +45,7 @@ int main(int argc, char *argv[]) { // SimulationJob job; // job.load("../" - // "ProblematicSimulationJobs/12:24_4_1_2021/" + // "ProblematicSimulationJobs/17:16_4_1_2021/" // "reduced_pattern_Single bar reduced model_simScenario.json"); // FormFinder ff; // ff.executeSimulation(std::make_shared(job), false, false, @@ -86,27 +86,13 @@ int main(int argc, char *argv[]) { std::string fullPatternsTestSetDirectory = // "/home/iason/Models/TestSet_validPatterns"; - "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - "1v_0v_2e_1e_1c_6fan/3/Valid"; + "/home/iason/Documents/PhD/Research/Approximating shapes with flat " + "patterns/Pattern_enumerator/Results/1v_0v_2e_1e_1c_6fan/3/Valid"; for (const auto &entry : filesystem::directory_iterator(fullPatternsTestSetDirectory)) { const auto filepath = - // std::filesystem::path("/home/iason/Models/valid_6777.ply"); - // std::filesystem::path( - // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - // "1v_0v_2e_1e_1c_6fan/3/Valid/222.ply"); - std::filesystem::path( - "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - "1v_0v_2e_1e_1c_6fan/3/Valid/624.ply"); - // std::filesystem::path( - // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - // "1v_0v_2e_1e_1c_6fan/2/Valid/33.ply"); - // std::filesystem::path( - // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - // "1v_0v_2e_1e_1c_6fan/3/Valid/431.ply"); - // std::filesystem::path( - // "/home/iason/Models/TestSet_validPatterns/59.ply"); - // entry.path(); + // std::filesystem::path(fullPatternsTestSetDirectory).append("305.ply"); + entry.path(); const auto filepathString = filepath.string(); // Use only the base triangle version const std::string tiledSuffix = "_tiled.ply"; @@ -114,6 +100,7 @@ int main(int argc, char *argv[]) { tiledSuffix.size(), tiledSuffix) == 0) { continue; } + std::cout << "Full pattern:" << filepathString << std::endl; FlatPattern pattern(filepathString); pattern.setLabel(filepath.stem().string()); pattern.scale(0.03); @@ -130,11 +117,11 @@ int main(int argc, char *argv[]) { // // first in the reducedModels vector // optimizationExcludedEi.insert(0); // } - std::cout << "Full pattern:" << filepathString << std::endl; - FlatPattern cp; - cp.copy(*reducedModels[0]); - optimizer.initializePatterns(cp, *reducedModels[0], optimizationExcludedEi); - optimizer.optimize({ReducedModelOptimizer::SimulationScenario::Bending}); + // FlatPattern cp; + // cp.copy(*reducedModels[0]); + optimizer.initializePatterns(pattern, *reducedModels[0], + optimizationExcludedEi); + optimizer.optimize(); // } } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 3495fe1..bac17b0 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -4,6 +4,9 @@ #include "gradientDescent.h" #include "simulationhistoryplotter.hpp" #include "trianglepattterntopology.hpp" +#include +#include +#include const bool gShouldDraw = true; @@ -28,6 +31,8 @@ double g_innerHexagonInitialPos = 0; bool g_optimizeInnerHexagonSize{false}; std::vector firstOptimizationRoundResults; int g_firstRoundIterationIndex{0}; +double minY{std::numeric_limits::max()}; +std::vector minX; // struct OptimizationCallback { // double operator()(const size_t &iterations, const Eigen::VectorXd &x, @@ -120,6 +125,11 @@ double ReducedModelOptimizer::computeError( reducedPatternResults.displacements[reducedModelVi][0], reducedPatternResults.displacements[reducedModelVi][1], reducedPatternResults.displacements[reducedModelVi][2]); + if (!std::isfinite(vertexDisplacement[0]) || + !std::isfinite(vertexDisplacement[1]) || + !std::isfinite(vertexDisplacement[2])) { + return std::numeric_limits::max(); + } Eigen::Vector3d errorVector = Eigen::Vector3d( optimalReducedPatternDisplacements.row(reducedModelVi)) - @@ -184,6 +194,12 @@ void updateMesh(long n, const double *x) { } } +double ReducedModelOptimizer::objective(double x0, double x1, double x2, + double x3) { + std::vector x{x0, x1, x2, x3}; + return ReducedModelOptimizer::objective(4, x.data()); +} + double ReducedModelOptimizer::objective(long n, const double *x) { std::cout.precision(17); @@ -205,23 +221,29 @@ double ReducedModelOptimizer::objective(long n, const double *x) { double error = 0; for (const int simulationScenarioIndex : g_simulationScenarioIndices) { SimulationResults reducedModelResults = simulator.executeSimulation( - g_reducedPatternSimulationJob[simulationScenarioIndex], false, false); + g_reducedPatternSimulationJob[simulationScenarioIndex]); error += computeError( reducedModelResults, g_optimalReducedModelDisplacements[simulationScenarioIndex]); } + if (error < minY) { + minY = error; + minX = std::vector({x[0], x[1], x[2], x[3]}); + } // compute error and return it gObjectiveValueHistory.push_back(error); - auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), - gObjectiveValueHistory.size()); - std::vector colors(gObjectiveValueHistory.size(), 2); - if (g_firstRoundIterationIndex != 0) { - for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), - [](double &c) { c = 0.7; }); - } - gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); + // auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), + // gObjectiveValueHistory.size()); + // std::vector colors(gObjectiveValueHistory.size(), 2); + // if (g_firstRoundIterationIndex != 0) { + // for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), + // [](double &c) { c = 0.7; }); + // } + // gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); std::cout << std::endl; + SimulationResultsReporter::createPlot("Number of Steps", "Objective value", + gObjectiveValueHistory); return error; } @@ -503,14 +525,15 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( gObjectiveValueHistory.clear(); const size_t n = g_optimizeInnerHexagonSize ? 5 : 4; - const size_t npt = 2 * n; - // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; - assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); - assert(npt <= 2 * n + 1 && "The choice of the number of interpolation " - "conditions is not recommended."); + // const size_t npt = (n + 1) * (n + 2) / 2; + // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; + // assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); + // assert(npt <= 2 * n + 1 && "The choice of the number of interpolation " + // "conditions is not recommended."); // Set initial guess of solution - std::vector x(n, 50); + const size_t initialGuess = 1; + std::vector x(n, initialGuess); if (g_optimizeInnerHexagonSize) { x[n - 1] = g_innerHexagonInitialPos; } @@ -523,50 +546,76 @@ Eigen::VectorXd ReducedModelOptimizer::runOptimization( // {initialGuess(0), initialGuess(1), initialGuess(2), // initialGuess(3)}; const double xMin = 1e-2; - const double xMax = 500; + const double xMax = 1e2; // assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) { // return d >= xMax || d <= xMin; // })); - std::vector xLow(x.size(), xMin); - std::vector xUpper(x.size(), xMax); - if (g_optimizeInnerHexagonSize) { - xLow[n - 1] = 0.1; - xUpper[n - 1] = 0.9; - } - const double maxX = *std::max_element( - x.begin(), x.end(), - [](const double &a, const double &b) { return abs(a) < abs(b); }); + // std::vector xLow(x.size(), xMin); + // std::vector xUpper(x.size(), xMax); + // if (g_optimizeInnerHexagonSize) { + // xLow[n - 1] = 0.1; + // xUpper[n - 1] = 0.9; + // } + // const double maxX = *std::max_element( + // x.begin(), x.end(), + // [](const double &a, const double &b) { return abs(a) < abs(b); }); // const double rhobeg = std::min(0.95, 0.2 * maxX); - const double rhobeg = 1; - // const double rhobeg = 10; - const double rhoend = rhobeg * 1e-6; - const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2; - std::vector w(wSize); - // const size_t maxFun = 10 * (x.size() ^ 2); - const size_t maxFun = 200; - bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), - rhobeg, rhoend, maxFun, w.data()); - std::cout << "Finished first optimization round" << std::endl; - firstOptimizationRoundResults.resize(6); - for (int simulationScenarioIndex = SimulationScenario::Axial; - simulationScenarioIndex != - SimulationScenario::NumberOfSimulationScenarios; - simulationScenarioIndex++) { - SimulationResults reducedModelResults = simulator.executeSimulation( - g_reducedPatternSimulationJob[simulationScenarioIndex], false, false); - reducedModelResults.setLabelPrefix("FirstRound"); - firstOptimizationRoundResults[simulationScenarioIndex] = - std::move(reducedModelResults); - } - g_firstRoundIterationIndex = gObjectiveValueHistory.size(); - bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), - rhobeg * 1e-3, rhoend, maxFun, w.data()); - std::cout << "Finished second optimization round" << std::endl; + // double rhobeg = 1; + // double rhoend = rhobeg * 1e-8; + // const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2; + // std::vector w(wSize); + // const size_t maxFun = std::min(100.0 * (x.size() + 1), 1000.0); + const size_t maxFun = 150; + + double (*objF)(double, double, double, double) = &objective; + auto start = std::chrono::system_clock::now(); + dlib::function_evaluation result = dlib::find_min_global( + objF, {xMin, xMin, xMin, xMin}, {xMax, xMax, xMax, xMax}, + dlib::max_function_calls(maxFun)); + auto end = std::chrono::system_clock::now(); + auto elapsed = std::chrono::duration_cast(end - start); + std::cout << "Finished optimization with dlib" << endl; + std::cout << "Solution x:" << endl; + std::cout << result.x << endl; + std::cout << "Solution y:" << endl; + std::cout << result.y << endl; + std::cout << minY << endl; + std::cout << "Time(sec):" << elapsed.count() << std::endl; + std::cout << "Max function evaluations:" << maxFun << std::endl; + std::cout << "[" + std::to_string(xMin) + "," + std::to_string(xMax) + "]" + << std::endl; + std::cout << "Initial guess:" << initialGuess << std::endl; + // const size_t maxFun = 200; + // bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), + // rhobeg, rhoend, maxFun, w.data()); + // std::cout << "Finished first optimization round" << std::endl; + // firstOptimizationRoundResults.resize(6); + // for (int simulationScenarioIndex = SimulationScenario::Axial; + // simulationScenarioIndex != + // SimulationScenario::NumberOfSimulationScenarios; + // simulationScenarioIndex++) { + // SimulationResults reducedModelResults = simulator.executeSimulation( + // g_reducedPatternSimulationJob[simulationScenarioIndex], false, + // false); + // reducedModelResults.setLabelPrefix("FirstRound"); + // firstOptimizationRoundResults[simulationScenarioIndex] = + // std::move(reducedModelResults); + // } + // g_firstRoundIterationIndex = gObjectiveValueHistory.size(); + // rhobeg *= 1e1; + // // rhoend *= 1e2; + // bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(), + // rhobeg, rhoend, maxFun, w.data()); + // std::cout << "Finished second optimization round" << std::endl; Eigen::VectorXd eigenX(x.size(), 1); for (size_t xi = 0; xi < x.size(); xi++) { - eigenX(xi) = x[xi]; + eigenX(xi) = minX[xi]; + // eigenX(xi) = x[xi]; + // eigenX(xi) = result.x(xi); } + // for (size_t xi = 0; xi < x.size(); xi++) { + // } return eigenX; } @@ -815,8 +864,7 @@ void ReducedModelOptimizer::runBeamOptimization() { fullModelResults.registerForDrawing(); SimulationResults reducedModelResults = simulator.executeSimulation( - g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1], false, - false); + g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]); double error = computeError( reducedModelResults, g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); @@ -840,13 +888,13 @@ void ReducedModelOptimizer::visualizeResults( pFullPatternSimulationJob->registerForDrawing( m_pFullPatternSimulationMesh->getLabel()); SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob, false); + simulator.executeSimulation(pFullPatternSimulationJob); fullModelResults.registerForDrawing(); fullModelResults.saveDeformedModel(); const std::shared_ptr &pReducedPatternSimulationJob = g_reducedPatternSimulationJob[simulationScenarioIndex]; SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob, false, false); + simulator.executeSimulation(pReducedPatternSimulationJob); error += computeError( reducedModelResults, g_optimalReducedModelDisplacements[simulationScenarioIndex]); @@ -854,39 +902,46 @@ void ReducedModelOptimizer::visualizeResults( << simulationScenarioStrings[simulationScenarioIndex] << " is " << error << std::endl; reducedModelResults.registerForDrawing(); - firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); - reducedModelResults.saveDeformedModel(); + // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); + // reducedModelResults.saveDeformedModel(); // registerWorldAxes(); + const std::string screenshotFilename = + "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" + + m_pFullPatternSimulationMesh->getLabel() + "_" + + simulationScenarioStrings[simulationScenarioIndex]; polyscope::show(); + polyscope::screenshot(screenshotFilename, false); fullModelResults.unregister(); reducedModelResults.unregister(); - firstOptimizationRoundResults[simulationScenarioIndex].unregister(); + // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); } } void ReducedModelOptimizer::optimize( - const std::vector &simulationScenariosIndices) { + const std::vector &simulationScenarios) { + + g_simulationScenarioIndices = simulationScenarios; + if (g_simulationScenarioIndices.empty()) { + g_simulationScenarioIndices = { + SimulationScenario::Axial, SimulationScenario::Shear, + SimulationScenario::Bending, SimulationScenario::Dome, + SimulationScenario::Saddle}; + } std::vector> simulationJobs = createScenarios(m_pFullPatternSimulationMesh); g_optimalReducedModelDisplacements.resize(6); g_reducedPatternSimulationJob.resize(6); g_firstRoundIterationIndex = 0; + minY = std::numeric_limits::max(); // polyscope::removeAllStructures(); - for (int simulationScenarioIndex = SimulationScenario::Axial; - simulationScenarioIndex != - SimulationScenario::NumberOfSimulationScenarios; - simulationScenarioIndex++) { + for (int simulationScenarioIndex : g_simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = simulationJobs[simulationScenarioIndex]; - // fullPatternSimulationJob.mesh->savePly( - // "Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply"); SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob); - // fullModelResults.label = - // "fullModel_" + SimulationScenarioStrings[simulationScenarioIndex]; - // fullModelResults.registerForDrawing(fullPatternSimulationJob); g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( m_pReducedPatternSimulationMesh->VN(), 3); computeDesiredReducedModelDisplacements( @@ -901,14 +956,7 @@ void ReducedModelOptimizer::optimize( std::make_shared(reducedPatternSimulationJob); } - g_simulationScenarioIndices = simulationScenariosIndices; - if (simulationScenariosIndices.empty()) { - g_simulationScenarioIndices = { - SimulationScenario::Axial, SimulationScenario::Shear, - SimulationScenario::Bending, SimulationScenario::Dome, - SimulationScenario::Saddle}; - } - - runOptimization(&objective); + Eigen::VectorXd optimalParameters = runOptimization(&objective); + updateMesh(4, optimalParameters.data()); visualizeResults(simulationJobs, g_simulationScenarioIndices); } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index eb7a52c..78b0561 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -32,7 +32,8 @@ public: }; inline static const std::string simulationScenarioStrings[] = { "Axial", "Shear", "Bending", "Double", "Saddle"}; - void optimize(const std::vector &simulationScenarios); + void optimize(const std::vector &simulationScenarios = + std::vector()); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); @@ -55,6 +56,8 @@ public: static void runSimulation(const std::string &filename, std::vector &x); + static double objective(double x0, double x1, double x2, double x3); + private: void visualizeResults(const std::vector>