diff --git a/CMakeLists.txt b/CMakeLists.txt index 2b725de..9467c15 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,16 +37,15 @@ download_project(PROJ MYSOURCES endif() - -#Polyscope -download_project(PROJ POLYSCOPE - GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git - GIT_TAG master - PREFIX ${CMAKE_CURRENT_SOURCE_DIR}/build/external/ - ${UPDATE_DISCONNECTED_IF_AVAILABLE} -) -add_subdirectory(${POLYSCOPE_SOURCE_DIR}) -add_compile_definitions(POLYSCOPE_DEFINED) +##Polyscope +#download_project(PROJ POLYSCOPE +# GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git +# GIT_TAG master +# PREFIX ${CMAKE_CURRENT_SOURCE_DIR}/build/external/ +# ${UPDATE_DISCONNECTED_IF_AVAILABLE} +#) +#add_subdirectory(${POLYSCOPE_SOURCE_DIR}) +#add_compile_definitions(POLYSCOPE_DEFINED) #dlib download_project(PROJ DLIB @@ -89,4 +88,5 @@ target_include_directories(${PROJECT_NAME} ) link_directories(${MYSOURCES_SOURCE_DIR}/boost_graph/libs) -target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen matplot dlib::dlib) +#target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen matplot dlib::dlib) +target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib) diff --git a/src/main.cpp b/src/main.cpp index eee3863..63bcf7e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -2,9 +2,6 @@ #include "csvfile.hpp" #include "edgemesh.hpp" #include "flatpattern.hpp" -#include "polyscope/curve_network.h" -#include "polyscope/point_cloud.h" -#include "polyscope/polyscope.h" #include "reducedmodeloptimizer.hpp" #include "simulationhistoryplotter.hpp" #include "trianglepattterntopology.hpp" @@ -16,6 +13,11 @@ #include #include +#ifdef POLYSCOPE_DEFINED +#include "polyscope/curve_network.h" +#include "polyscope/point_cloud.h" +#include "polyscope/polyscope.h" +#endif int main(int argc, char *argv[]) { if (argc < 3) { std::cerr << "Specify at least the two pattern filepaths to be " @@ -49,7 +51,10 @@ int main(int argc, char *argv[]) { const bool input_numberOfFunctionCallsDefined = argc >= 4; settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; - settings_optimization.shouldNormalizeObjectiveValue = true; + settings_optimization.normalizationStrategy = + ReducedModelOptimizer::Settings::NormalizationStrategy::Epsilon; + settings_optimization.normalizationParameter = 0.003; + settings_optimization.solutionAccuracy = 0.01; // Optimize pair const std::string pairName = @@ -65,21 +70,24 @@ int main(int argc, char *argv[]) { const bool input_resultDirectoryDefined = argc >= 5; std::string optimizationResultsDirectory = input_resultDirectoryDefined ? argv[4] : "OptimizationResults"; + std::string resultsOutputDir; if (optimizationResults.numberOfSimulationCrashes != 0) { const auto crashedJobsDirPath = std::filesystem::path(optimizationResultsDirectory) .append("CrashedJobs") .append(pairName); std::filesystem::create_directories(crashedJobsDirPath); - optimizationResults.save(crashedJobsDirPath.string()); + resultsOutputDir = crashedJobsDirPath.string(); } else { std::filesystem::path convergedJobsDirPath( std::filesystem::path(optimizationResultsDirectory) .append("ConvergedJobs") .append(pairName)); std::filesystem::create_directories(convergedJobsDirPath); - optimizationResults.save(convergedJobsDirPath.string()); + resultsOutputDir = convergedJobsDirPath.string(); } + optimizationResults.save(resultsOutputDir); + // Write results in csv csvFile csv_results({}, false); // csvFile csv_results(std::filesystem::path(dirPath_thisOptimization) // .append("results.csv") diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 0edbe87..6aeaeec 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -9,7 +9,7 @@ struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; std::vector> fullPatternDisplacements; - std::vector reducedPatternMaximumDisplacementNormSum; + std::vector objectiveNormalizationValues; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; std::unordered_map @@ -46,16 +46,10 @@ double ReducedModelOptimizer::computeError( const double rawError = computeRawError(reducedPatternDisplacements, fullPatternDisplacements, reducedToFullInterfaceViMap); - if (global.optimizationSettings.shouldNormalizeObjectiveValue) { - assert(rawError <= normalizationFactor); - if (rawError > normalizationFactor) { - std::cout << "Normalized error is bigger than one." << std::endl; - } - return rawError / - normalizationFactor; // std::max(interfaceDisplacementsNormSum, - // 0.00003); + if (global.optimizationSettings.normalizationStrategy != + Settings::NormalizationStrategy::NonNormalized) { + return rawError / normalizationFactor; } - return rawError; } @@ -66,24 +60,22 @@ double ReducedModelOptimizer::computeRawError( &reducedToFullInterfaceViMap) { double error = 0; for (const auto reducedFullViPair : reducedToFullInterfaceViMap) { - // const auto pos = - // g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP(); - // std::cout << "Interface vi " << reducedModelVi << " is at position " - // << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; + const VertexIndex reducedModelVi = reducedFullViPair.first; Eigen::Vector3d reducedVertexDisplacement( - reducedPatternDisplacements[reducedFullViPair.first][0], - reducedPatternDisplacements[reducedFullViPair.first][1], - reducedPatternDisplacements[reducedFullViPair.first][2]); + reducedPatternDisplacements[reducedModelVi][0], + reducedPatternDisplacements[reducedModelVi][1], + reducedPatternDisplacements[reducedModelVi][2]); if (!std::isfinite(reducedVertexDisplacement[0]) || !std::isfinite(reducedVertexDisplacement[1]) || !std::isfinite(reducedVertexDisplacement[2])) { std::cout << "Displacements are not finite" << std::endl; std::terminate(); } + const VertexIndex fullModelVi = reducedFullViPair.second; Eigen::Vector3d fullVertexDisplacement( - fullPatternDisplacements[reducedFullViPair.second][0], - fullPatternDisplacements[reducedFullViPair.second][1], - fullPatternDisplacements[reducedFullViPair.second][2]); + fullPatternDisplacements[fullModelVi][0], + fullPatternDisplacements[fullModelVi][1], + fullPatternDisplacements[fullModelVi][2]); Eigen::Vector3d errorVector = fullVertexDisplacement - reducedVertexDisplacement; // error += errorVector.squaredNorm(); @@ -179,7 +171,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // << std::endl; // run simulations - double error = 0; + double totalError = 0; FormFinder simulator; FormFinder::Settings simulationSettings; // simulationSettings.shouldDraw = true; @@ -192,75 +184,54 @@ double ReducedModelOptimizer::objective(long n, const double *x) { simulationSettings); // std::string filename; if (!reducedModelResults.converged) { - error += std::numeric_limits::max(); + totalError += std::numeric_limits::max(); global.numOfSimulationCrashes++; bool beVerbose = false; // TODO: use an extern or data member for that if (beVerbose) { std::cout << "Failed simulation" << std::endl; } } else { - const double thisSimulationScenarioError = - computeError(reducedModelResults.displacements, - global.fullPatternDisplacements[simulationScenarioIndex], - global.reducedToFullInterfaceViMap, - global.reducedPatternMaximumDisplacementNormSum - [simulationScenarioIndex]); - if (thisSimulationScenarioError > 1) { - std::cout << "Simulation scenario " - << simulationScenarioStrings[simulationScenarioIndex] - << " results in an error bigger than one." << std::endl; - for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { - std::cout << "x[" + std::to_string(parameterIndex) + "]=" - << x[parameterIndex] << std::endl; - } + double simulationScenarioError = computeError( + reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, + global.objectiveNormalizationValues[simulationScenarioIndex]); - //#ifdef POLYSCOPE_DEFINED - // ReducedModelOptimizer::visualizeResults( - // global.fullPatternSimulationJobs[simulationScenarioIndex], - // global.reducedPatternSimulationJobs[simulationScenarioIndex], - // global.reducedToFullInterfaceViMap, false); + // 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; + // } + // } - // ReducedModelOptimizer::visualizeResults( - // global.fullPatternSimulationJobs[simulationScenarioIndex], - // std::make_shared( - // reducedPatternMaximumDisplacementSimulationJobs - // [simulationScenarioIndex]), - // global.reducedToFullInterfaceViMap, true); - // polyscope::removeAllStructures(); - //#endif // POLYSCOPE_DEFINED - } - error += thisSimulationScenarioError; + //#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; } - // std::ofstream out(filename, std::ios_base::app); - // auto pMesh = - // global - // .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] - // ->pMesh; - // for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) - // { - // out << "x[" + std::to_string(parameterIndex) + "]=" << - // x[parameterIndex] - // << std::endl; - // } - // out << pMesh->elements[0].dimensions.toString() + "\n" + - // pMesh->elements[0].material.toString() + " \nA=" - // << pMesh->elements[0].A << " \nratio=" - // << pMesh->elements[0].dimensions.b / - // pMesh->elements[0].dimensions.h - // << " \naxialRig:" << pMesh->elements[0].rigidity.axial - // << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional - // << " \nfirstBendingRig:" << - // pMesh->elements[0].rigidity.firstBending - // << " \nsecondBendingRig:" << - // pMesh->elements[0].rigidity.secondBending - // << " \nscenario:" + - // simulationScenarioStrings[simulationScenarioIndex] + - // "\n\n"; - // out.close(); } // std::cout << error << std::endl; - if (error < global.minY) { - global.minY = error; + if (totalError < global.minY) { + global.minY = totalError; global.minX.assign(x, x + n); } // if (++global.numberOfFunctionCalls %100== 0) { @@ -269,7 +240,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { //} // compute error and return it - global.gObjectiveValueHistory.push_back(error); + // global.gObjectiveValueHistory.push_back(error); // auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), // gObjectiveValueHistory.size()); // std::vector colors(gObjectiveValueHistory.size(), 2); @@ -282,7 +253,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // value", // gObjectiveValueHistory); - return error; + return totalError; } void ReducedModelOptimizer::createSimulationMeshes( @@ -578,11 +549,73 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( #if POLYSCOPE_DEFINED void ReducedModelOptimizer::visualizeResults( + const std::vector> + &fullPatternSimulationJobs, + const std::vector> + &reducedPatternSimulationJobs, + const std::vector &simulationScenarios, + const std::unordered_map + &reducedToFullInterfaceViMap) { + FormFinder simulator; + std::shared_ptr pFullPatternSimulationMesh = + fullPatternSimulationJobs[0]->pMesh; + pFullPatternSimulationMesh->registerForDrawing(); + pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() + + "_undeformed.ply"); + reducedPatternSimulationJobs[0]->pMesh->savePly( + 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 != + Settings::NormalizationStrategy::NonNormalized) { + normalizationFactor = + global.objectiveNormalizationValues[simulationScenarioIndex]; + } + reducedModelResults.saveDeformedModel(); + fullModelResults.saveDeformedModel(); + double error = computeError( + reducedModelResults.displacements, fullModelResults.displacements, + reducedToFullInterfaceViMap, normalizationFactor); + std::cout << "Error of simulation scenario " + << simulationScenarioStrings[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() + "_" + + simulationScenarioStrings[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, - const bool &draw = true) { + &reducedToFullInterfaceViMap) { FormFinder simulator; std::shared_ptr pFullPatternSimulationMesh = pFullPatternSimulationJob->pMesh; @@ -609,19 +642,6 @@ void ReducedModelOptimizer::visualizeResults( global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); std::cout << "Error is " << error << std::endl; reducedModelResults.registerForDrawing(); - if (draw) { - polyscope::show(); - const std::string screenshotFilename = - "/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/RodModelOptimizationForPatterns/Results/" - "Images/" + - pFullPatternSimulationMesh->getLabel() + "_" + "noScenarioName"; - // simulationScenarioStrings[simulationScenarioIndex]; - polyscope::screenshot(screenshotFilename, false); - fullModelResults.unregister(); - reducedModelResults.unregister(); - pFullPatternSimulationMesh->unregister(); - } } #endif // POLYSCOPE_DEFINED @@ -651,46 +671,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::runOptimization(const Settings &settings) { global.gObjectiveValueHistory.clear(); - - // g_optimizeInnerHexagonSize ? 5: 4; - // 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 - - // const size_t initialGuess = 1; - // std::vector x(n, initialGuess); - // if (global.optimizeInnerHexagonSize) { - // x[n - 1] = global.g_innerHexagonInitialPos; - //} - /*if (!initialGuess.empty()) { - x = g_optimizationInitialGuess; - }*/ // {0.10000000000000 001, - // 2, 1.9999999971613847, 6.9560343643347764}; - // {1, 5.9277}; - // {0.0001, 2, 2.000000005047502, 1.3055270196964464}; - // {initialGuess(0), initialGuess(1), initialGuess(2), - // initialGuess(3)}; - // 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); }); - // const double rhobeg = std::min(0.95, 0.2 * maxX); - // 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); dlib::matrix xMin(global.numberOfOptimizationParameters); dlib::matrix xMax(global.numberOfOptimizationParameters); for (int i = 0; i < global.numberOfOptimizationParameters; i++) { @@ -720,17 +700,59 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) { results.numberOfSimulationCrashes = global.numOfSimulationCrashes; results.x = global.minX; results.objectiveValue = global.minY; - if (settings.shouldNormalizeObjectiveValue) { - global.optimizationSettings.shouldNormalizeObjectiveValue = false; + if (global.minY != result.y) { + std::cerr << "Global min objective is not equal to result objective" + << std::endl; + } + + // Compute raw objective value + if (global.optimizationSettings.normalizationStrategy != + Settings::NormalizationStrategy::NonNormalized) { + auto temp = global.optimizationSettings.normalizationStrategy; + global.optimizationSettings.normalizationStrategy = + Settings::NormalizationStrategy::NonNormalized; if (global.optimizeInnerHexagonSize) { - results.objectiveValue = objective(4, results.x.data()); + results.rawObjectiveValue = objective(4, results.x.data()); } else { - results.objectiveValue = + results.rawObjectiveValue = objective(results.x[0], results.x[1], results.x[2]); } - global.optimizationSettings.shouldNormalizeObjectiveValue = true; + global.optimizationSettings.normalizationStrategy = temp; + + } else { + results.rawObjectiveValue = results.objectiveValue; } + + // Compute obj value per simulation scenario + updateMesh(results.x.size(), results.x.data()); + results.objectiveValuePerSimulationScenario.resize( + NumberOfSimulationScenarios); + FormFinder::Settings simulationSettings; + FormFinder simulator; + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < NumberOfSimulationScenarios; + simulationScenarioIndex++) { + SimulationResults reducedModelResults = simulator.executeSimulation( + global.reducedPatternSimulationJobs[simulationScenarioIndex], + simulationSettings); + const double error = computeError( + reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap, + global.objectiveNormalizationValues[simulationScenarioIndex]); + + results.objectiveValuePerSimulationScenario[simulationScenarioIndex] = + error; + } + // if (result.y != + // std::accumulate(results.objectiveValuePerSimulationScenario.begin(), + // results.objectiveValuePerSimulationScenario.end(), 0)) + // { + // std::cerr + // << "Sum of per scenario objectives is not equal to result objective" + // << std::endl; + // } results.time = elapsed.count() / 1000.0; const bool printDebugInfo = false; if (printDebugInfo) { @@ -739,33 +761,6 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) { // std::cout << result.x << endl; std::cout << "Objective value:" << global.minY << 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 << "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; return results; } @@ -784,37 +779,7 @@ ReducedModelOptimizer::createScenarios( const double forceMagnitude = 1; // Assuming the patterns lays on the x-y plane const CoordType patternPlaneNormal(0, 0, 1); - // Make the first interface node lay on the x axis - // const size_t fullPatternFirstInterfaceNodeIndex = - // m_fullPatternOppositeInterfaceViMap.begin()->second; - // CoordType fullPatternFirstInterfaceNodePosition = - // m_pFullModelSimulationMesh->vert[fullPatternFirstInterfaceNodeIndex].cP(); - // CoordType centerOfMass(0, 0, 0); - // for (size_t vi = 0; vi < pMesh->VN(); vi++) { - // centerOfMass = centerOfMass + pMesh->vert[vi].P(); - // } - // centerOfMass /= pMesh->VN(); - // vcg::tri::UpdatePosition::Translate( - // *m_pFullModelSimulationMesh, -centerOfMass); - // vcg::tri::UpdatePosition::Translate( - // *m_pReducedPatternSimulationMesh, centerOfMass); - - // const vcg::Matrix33d R = vcg::RotationMatrix( - // fullPatternFirstInterfaceNodePosition, - // CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false); - // std::for_each(m_pFullModelSimulationMesh->vert.begin(), - // m_pFullModelSimulationMesh->vert.end(), [&](auto &v) { - // v.P() = R * v.P(); - // v.N() = R * v.N(); - // }); - // std::for_each(m_pReducedPatternSimulationMesh->vert.begin(), - // m_pReducedPatternSimulationMesh->vert.end(), [&](auto &v) { - // v.P() = R * v.P(); - // v.N() = R * v.N(); - // }); - - // m_pFullModelSimulationMesh->updateEigenEdgeAndVertices(); - // m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); + // Assumes that the first interface node lay on the x axis //// Axial SimulationScenario scenarioName = SimulationScenario::Axial; @@ -945,89 +910,87 @@ ReducedModelOptimizer::createScenarios( return scenarios; } -// void ReducedModelOptimizer::runBeamOptimization() { -// // load beams -// VCGEdgeMesh fullBeam; -// fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); -// VCGEdgeMesh reducedBeam; -// reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply"); -// fullBeam.registerForDrawing(); -// reducedBeam.registerForDrawing(); -// // polyscope::show(); -// // maps -// std::unordered_map displacementReducedToFullMap; -// displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2; -// g_reducedToFullViMap = displacementReducedToFullMap; -// std::unordered_map jobFullToReducedMap; -// jobFullToReducedMap[0] = 0; -// jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1; +void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { + assert(global.optimizationSettings.shouldNormalizeObjectiveValue); + if (global.optimizationSettings.normalizationStrategy == + Settings::NormalizationStrategy::Epsilon || + global.optimizationSettings.normalizationStrategy == + Settings::NormalizationStrategy::EqualDisplacements) { -// // full model simuilation job -// auto pFullPatternSimulationMesh = -// std::make_shared(fullBeam); -// pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, -// 0.02}); pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); -// std::unordered_map> fixedVertices; -// fixedVertices[0] = ::unordered_set({0, 1, 2, 3, 4, 5}); -// std::unordered_map forces; -// forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0}); -// const std::string fullBeamSimulationJobLabel = "Pull_Z"; -// std::shared_ptr pFullModelSimulationJob = -// make_shared(SimulationJob(pFullPatternSimulationMesh, -// fullBeamSimulationJobLabel, -// fixedVertices, forces)); -// auto fullModelResults = -// formFinder.executeSimulation(pFullModelSimulationJob); + // Compute the sum of the displacement norms + std::vector fullPatternDisplacementNormSum( + NumberOfSimulationScenarios); + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + double displacementNormSum = 0; + for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + const Vector6d &vertexDisplacement = + global.fullPatternDisplacements[simulationScenarioIndex] + [interfaceViPair.second]; + displacementNormSum += vertexDisplacement.getTranslation().norm(); + } + fullPatternDisplacementNormSum[simulationScenarioIndex] = + displacementNormSum; + } -// // Optimal reduced model displacements -// const size_t numberOfSimulationScenarios = 1; -// g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios); -// g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize( -// reducedBeam.VN(), 3); -// computeDesiredReducedModelDisplacements( -// fullModelResults, displacementReducedToFullMap, -// g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); - -// // reduced model simuilation job -// auto reducedSimulationMesh = std::make_shared(reducedBeam); -// reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02}); -// reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9); -// g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios); -// SimulationJob reducedSimJob; -// computeReducedModelSimulationJob(*pFullModelSimulationJob, -// jobFullToReducedMap, reducedSimJob); -// reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] = -// reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1; -// g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] = -// make_shared( -// reducedSimulationMesh, fullBeamSimulationJobLabel, -// reducedSimJob.constrainedVertices, -// reducedSimJob.nodalExternalForces); -// initializeOptimizationParameters(reducedSimulationMesh); - -// // const std::string simulationJobsPath = "SimulationJobs"; -// // std::filesystem::create_directory(simulationJobsPath); -// // g_reducedPatternSimulationJob[0].save(simulationJobsPath); -// // g_reducedPatternSimulationJob[0].load( -// // std::filesystem::path(simulationJobsPath) -// // .append(g_reducedPatternSimulationJob[0].mesh->getLabel() + -// // "_simScenario.json")); - -// runOptimization({}, &objective); - -// fullModelResults.registerForDrawing(); -// SimulationResults reducedModelResults = simulator.executeSimulation( -// g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]); -// double error = computeError( -// reducedModelResults, -// g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]); -// reducedModelResults.registerForDrawing(); -// std::cout << "Error between beams:" << error << endl; -// // registerWorldAxes(); -// polyscope::show(); -// fullModelResults.unregister(); -// reducedModelResults.unregister(); -//} + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + if (global.optimizationSettings.normalizationStrategy == + Settings::NormalizationStrategy::Epsilon) { + const double epsilon = + global.optimizationSettings.normalizationParameter; + if (epsilon > fullPatternDisplacementNormSum[simulationScenarioIndex]) { + // std::cout << "Epsilon used in " + // << + // simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; + } + global.objectiveNormalizationValues[simulationScenarioIndex] = std::max( + fullPatternDisplacementNormSum[simulationScenarioIndex], epsilon); + // displacementNormSum; + } else { // NormalizationStrategy::EqualDisplacements + global.objectiveNormalizationValues[simulationScenarioIndex] = + fullPatternDisplacementNormSum[simulationScenarioIndex] / + fullPatternDisplacementNormSum[SimulationScenario::Shear]; + } + } + } else if (global.optimizationSettings.normalizationStrategy == + Settings::NormalizationStrategy::MaxDisplacement) { + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + // Compute the quantities for normalizing the obj value + const double minB = global.optimizationSettings.xRanges[0].min; + const double maxRatio = global.optimizationSettings.xRanges[1].max; + const double minE = global.optimizationSettings.xRanges[2].min; + const double minHS = 0.3; + std::vector mostFlexibleOptimizationParameters{minB, maxRatio, + minE, minHS}; + if (global.optimizeInnerHexagonSize) { + updateMesh(4, mostFlexibleOptimizationParameters.data()); + } else { + updateMesh(3, mostFlexibleOptimizationParameters.data()); + } + reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] = + global.reducedPatternSimulationJobs[simulationScenarioIndex] + ->getCopy(); + reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] + .pMesh->setLabel("reduced_maxDisplacement"); + SimulationResults reducedModelResults = simulator.executeSimulation( + global.reducedPatternSimulationJobs[simulationScenarioIndex]); + const double errorOfMaxDisplacedReduced = computeRawError( + reducedModelResults.displacements, + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap); + const double errorOfNonDisplacedReduced = computeRawError( + std::vector(reducedModelResults.displacements.size(), + Vector6d(0)), + global.fullPatternDisplacements[simulationScenarioIndex], + global.reducedToFullInterfaceViMap); + const double displacementMultiplier = + global.optimizationSettings.normalizationParameter; + global.objectiveNormalizationValues[simulationScenarioIndex] = + displacementMultiplier * + std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced); + } + } +} ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( const Settings &optimizationSettings, @@ -1041,10 +1004,10 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( SimulationScenario::Saddle}; } - global.g_optimalReducedModelDisplacements.resize(6); - global.reducedPatternSimulationJobs.resize(6); - global.fullPatternDisplacements.resize(6); - global.reducedPatternMaximumDisplacementNormSum.resize(6); + global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios); + global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios); + global.fullPatternDisplacements.resize(NumberOfSimulationScenarios); + global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios); global.g_firstRoundIterationIndex = 0; global.minY = std::numeric_limits::max(); global.numOfSimulationCrashes = 0; @@ -1052,7 +1015,8 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.optimizationSettings = optimizationSettings; global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh); - reducedPatternMaximumDisplacementSimulationJobs.resize(6); + reducedPatternMaximumDisplacementSimulationJobs.resize( + NumberOfSimulationScenarios); // polyscope::removeAllStructures(); FormFinder::Settings simulationSettings; @@ -1071,40 +1035,11 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( reducedPatternSimulationJob); global.reducedPatternSimulationJobs[simulationScenarioIndex] = std::make_shared(reducedPatternSimulationJob); + } - if (optimizationSettings.shouldNormalizeObjectiveValue) { - // Compute the quantities for normalizing the obj value - const double minB = optimizationSettings.xRanges[0].min; - const double maxRatio = optimizationSettings.xRanges[1].max; - const double minE = optimizationSettings.xRanges[2].min; - const double minHS = 0.3; - std::vector mostFlexibleOptimizationParameters{minB, maxRatio, - minE, minHS}; - if (global.optimizeInnerHexagonSize) { - updateMesh(4, mostFlexibleOptimizationParameters.data()); - } else { - updateMesh(3, mostFlexibleOptimizationParameters.data()); - } - reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] = - global.reducedPatternSimulationJobs[simulationScenarioIndex] - ->getCopy(); - reducedPatternMaximumDisplacementSimulationJobs[simulationScenarioIndex] - .pMesh->setLabel("reduced_maxDisplacement"); - SimulationResults reducedModelResults = simulator.executeSimulation( - global.reducedPatternSimulationJobs[simulationScenarioIndex], - simulationSettings); - const double errorOfMaxDisplacedReduced = computeRawError( - reducedModelResults.displacements, fullModelResults.displacements, - global.reducedToFullInterfaceViMap); - const double errorOfNonDisplacedReduced = computeRawError( - std::vector(reducedModelResults.displacements.size(), - Vector6d(0)), - fullModelResults.displacements, global.reducedToFullInterfaceViMap); - const double displacementMultiplier = 2; - global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex] = - displacementMultiplier * - std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced); - } + if (global.optimizationSettings.normalizationStrategy != + Settings::NormalizationStrategy::NonNormalized) { + computeObjectiveValueNormalizationFactors(); } Results optResults = runOptimization(optimizationSettings); for (int simulationScenarioIndex : global.simulationScenarioIndices) { @@ -1113,12 +1048,9 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( optResults.reducedPatternSimulationJobs.push_back( global.reducedPatternSimulationJobs[simulationScenarioIndex]); } - optResults.draw(); - #ifdef POLYSCOPE_DEFINED -// visualizeResults( -// global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs, -// global.simulationScenarioIndices, global.reducedToFullInterfaceViMap); + optResults.draw(); #endif // POLYSCOPE_DEFINED + return optResults; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 2e93d4c..7317ce3 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -8,6 +8,9 @@ #include "matplot/matplot.h" #include +#ifdef POLYSCOPE_DEFINED +#include "polyscope/color_management.h" +#endif // POLYSCOPE_DEFINED using FullPatternVertexIndex = VertexIndex; using ReducedPatternVertexIndex = VertexIndex; @@ -21,6 +24,21 @@ class ReducedModelOptimizer { std::unordered_map nodeToSlot; std::unordered_map> slotToNode; std::vector initialGuess; +#ifdef POLYSCOPE_DEFINED + struct StaticColors { + glm::vec3 fullInitial; + glm::vec3 fullDeformed; + glm::vec3 reducedInitial; + glm::vec3 reducedDeformed; + StaticColors() { + fullInitial = {0.416666, 0.109804, 0.890196}; + fullDeformed = {0.583333, 0.890196, 0.109804}; + reducedInitial = {0.890196, 0.109804, 0.193138}; + reducedDeformed = {0.109804, 0.890196, 0.806863}; + } + }; + inline static StaticColors colors; +#endif // POLYSCOPE_DEFINED public: enum SimulationScenario { @@ -44,10 +62,19 @@ public: struct Results; struct Settings { + enum NormalizationStrategy { + NonNormalized, + Epsilon, + MaxDisplacement, + EqualDisplacements + }; + inline static vector normalizationStrategyStrings{ + "NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"}; std::vector xRanges; int numberOfFunctionCalls{100}; double solutionAccuracy{1e-2}; - bool shouldNormalizeObjectiveValue{true}; + NormalizationStrategy normalizationStrategy{Epsilon}; + double normalizationParameter{0.003}; std::string toString() const { std::string settingsString; @@ -61,7 +88,7 @@ public: settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) + " Accuracy=" + std::to_string(solutionAccuracy) + - " Norm=" + (shouldNormalizeObjectiveValue ? "yes" : "no"); + " Norm=" + normalizationStrategyStrings[normalizationStrategy]; return settingsString; } @@ -75,6 +102,7 @@ public: } os << "Function Calls"; os << "Solution Accuracy"; + os << "Normalization strategy"; // os << std::endl; } @@ -87,11 +115,211 @@ public: } os << numberOfFunctionCalls; os << solutionAccuracy; + os << normalizationStrategyStrings[normalizationStrategy] + "_" + + std::to_string(normalizationParameter); + } + }; + struct Results { + double time{-1}; + int numberOfSimulationCrashes{0}; + std::vector x; + double objectiveValue; + double rawObjectiveValue; + std::vector objectiveValuePerSimulationScenario; + std::vector> fullPatternSimulationJobs; + std::vector> reducedPatternSimulationJobs; + + void save(const string &saveToPath) const { + assert(std::filesystem::is_directory(saveToPath)); + + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + for (int simulationJobIndex = 0; + simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationJobIndex]; + std::filesystem::path simulationJobFolderPath( + std::filesystem::path(saveToPath) + .append(pFullPatternSimulationJob->label)); + std::filesystem::create_directory(simulationJobFolderPath); + const auto fullPatternDirectoryPath = + std::filesystem::path(simulationJobFolderPath).append("Full"); + std::filesystem::create_directory(fullPatternDirectoryPath); + pFullPatternSimulationJob->save(fullPatternDirectoryPath.string()); + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationJobIndex]; + const auto reducedPatternDirectoryPath = + std::filesystem::path(simulationJobFolderPath).append("Reduced"); + if (!std::filesystem::exists(reducedPatternDirectoryPath)) { + std::filesystem::create_directory(reducedPatternDirectoryPath); + } + pReducedPatternSimulationJob->save( + reducedPatternDirectoryPath.string()); + } + } + void load(const string &loadFromPath) { + assert(std::filesystem::is_directory(loadFromPath)); + + for (const auto &directoryEntry : + filesystem::directory_iterator(loadFromPath)) { + const auto simulationScenarioPath = directoryEntry.path(); + if (!std::filesystem::is_directory(simulationScenarioPath)) { + continue; + } + // Load reduced pattern files + for (const auto &fileEntry : filesystem::directory_iterator( + std::filesystem::path(simulationScenarioPath) + .append("Full"))) { + const auto filepath = fileEntry.path(); + if (filepath.extension() == ".json") { + SimulationJob job; + job.load(filepath.string()); + fullPatternSimulationJobs.push_back( + std::make_shared(job)); + } + } + + // Load full pattern files + for (const auto &fileEntry : filesystem::directory_iterator( + std::filesystem::path(simulationScenarioPath) + .append("Reduced"))) { + const auto filepath = fileEntry.path(); + if (filepath.extension() == ".json") { + SimulationJob job; + job.load(filepath.string()); + reducedPatternSimulationJobs.push_back( + std::make_shared(job)); + } + } + } + } +#if POLYSCOPE_DEFINED + void draw() const { + initPolyscope(); + FormFinder simulator; + assert(fullPatternSimulationJobs.size() == + reducedPatternSimulationJobs.size()); + fullPatternSimulationJobs[0]->pMesh->registerForDrawing( + colors.fullInitial); + reducedPatternSimulationJobs[0]->pMesh->registerForDrawing( + colors.reducedInitial); + + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + for (int simulationJobIndex = 0; + simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { + // Drawing of full pattern results + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationJobIndex]; + pFullPatternSimulationJob->registerForDrawing( + fullPatternSimulationJobs[0]->pMesh->getLabel()); + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob); + fullModelResults.registerForDrawing(colors.fullDeformed); + // Drawing of reduced pattern results + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationJobIndex]; + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob); + reducedModelResults.registerForDrawing(colors.reducedDeformed); + polyscope::options::programName = + fullPatternSimulationJobs[0]->pMesh->getLabel(); + polyscope::show(); + // Save a screensh + const std::string screenshotFilename = + "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/Results/Images/" + + pFullPatternSimulationJob->pMesh->getLabel() + "_" + + simulationScenarioStrings[simulationJobIndex]; + polyscope::screenshot(screenshotFilename, false); + fullModelResults.unregister(); + reducedModelResults.unregister(); + // double error = computeError( + // reducedModelResults.displacements,fullModelResults.displacements, + // ); + // std::cout << "Error of simulation scenario " + // << + // simulationScenarioStrings[simulationScenarioIndex] + // << " is " + // << error << std::endl; + } + } +#endif // POLYSCOPE_DEFINED + void saveMeshFiles() const { + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + assert(numberOfSimulationJobs != 0 && + fullPatternSimulationJobs.size() == + reducedPatternSimulationJobs.size()); + + fullPatternSimulationJobs[0]->pMesh->savePly( + "FullPattern_undeformed.ply"); + reducedPatternSimulationJobs[0]->pMesh->savePly( + "ReducedPattern_undeformed.ply"); + FormFinder simulator; + for (int simulationJobIndex = 0; + simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { + // Drawing of full pattern results + const std::shared_ptr &pFullPatternSimulationJob = + fullPatternSimulationJobs[simulationJobIndex]; + SimulationResults fullModelResults = + simulator.executeSimulation(pFullPatternSimulationJob); + fullModelResults.saveDeformedModel(); + + // Drawing of reduced pattern results + const std::shared_ptr &pReducedPatternSimulationJob = + reducedPatternSimulationJobs[simulationJobIndex]; + SimulationResults reducedModelResults = + simulator.executeSimulation(pReducedPatternSimulationJob); + reducedModelResults.saveDeformedModel(); + } + } + + void + writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization, + csvFile &os) { + os << "Total raw Obj value"; + os << "Total Obj value"; + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < + SimulationScenario::NumberOfSimulationScenarios; + simulationScenarioIndex++) { + os << "Obj value " + simulationScenarioStrings[simulationScenarioIndex]; + } + for (const ReducedModelOptimizer::xRange &range : + settings_optimization.xRanges) { + os << range.label; + } + os << "Time(s)"; + os << "#Crashes"; + } + + void + writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization, + csvFile &os) const { + os << rawObjectiveValue; + os << objectiveValue; + for (double scenarioObjValue : objectiveValuePerSimulationScenario) { + os << scenarioObjValue; + } + for (const double &optimalX : x) { + os << optimalX; + } + + for (int unusedXVarCounter = 0; + unusedXVarCounter < settings_optimization.xRanges.size() - x.size(); + unusedXVarCounter++) { + os << "-"; + } + + os << time; + if (numberOfSimulationCrashes == 0) { + os << "-"; + } else { + os << numberOfSimulationCrashes; + } } }; inline static const std::string simulationScenarioStrings[] = { - "Axial", "Shear", "Bending", "Double", "Saddle"}; + "Axial", "Shear", "Bending", "Dome", "Saddle"}; Results optimize(const Settings &xRanges, const std::vector &simulationScenarios = std::vector()); @@ -118,7 +346,7 @@ public: std::vector &x); static double objective(double x0, double x1, double x2, double x3); - static double objective(double b, double h, double E); + static double objective(double b, double r, double E); static std::vector> createScenarios(const std::shared_ptr &pMesh, @@ -138,20 +366,36 @@ public: &fullToReducedInterfaceViMap, std::unordered_map &fullPatternOppositeInterfaceViMap); - static void visualizeResults( - const std::shared_ptr &fullPatternSimulationJobs, - const std::shared_ptr &reducedPatternSimulationJobs, + static void + visualizeResults(const std::vector> + &fullPatternSimulationJobs, + const std::vector> + &reducedPatternSimulationJobs, + const std::vector &simulationScenarios, + const std::unordered_map + &reducedToFullInterfaceViMap); + static void registerResultsForDrawing( + const std::shared_ptr &pFullPatternSimulationJob, + const std::shared_ptr &pReducedPatternSimulationJob, const std::unordered_map - &reducedToFullInterfaceViMap, - const bool &shouldDraw); + &reducedToFullInterfaceViMap); - static double computeError(const std::vector &reducedPatternResults, - const std::vector &fullPatternResults, - const std::unordered_map - &reducedToFullInterfaceViMap, - const double &normalizationFactor); + static double + computeRawError(const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap); + + static double + computeError(const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor); private: static void computeDesiredReducedModelDisplacements( @@ -168,194 +412,8 @@ private: static void initializeOptimizationParameters(const std::shared_ptr &mesh); - static double - computeError(const SimulationResults &reducedPatternResults, - const Eigen::MatrixX3d &optimalReducedPatternDisplacements); static double objective(long n, const double *x); FormFinder simulator; - double static computeRawError( - const std::vector &reducedPatternDisplacements, - const std::vector &fullPatternDisplacements, - const std::unordered_map - &reducedToFullInterfaceViMap); + void computeObjectiveValueNormalizationFactors(); }; -struct ReducedModelOptimizer::Results { - double time{-1}; - int numberOfSimulationCrashes{0}; - std::vector x; - double objectiveValue; - std::vector> fullPatternSimulationJobs; - std::vector> reducedPatternSimulationJobs; - - void save(const string &saveToPath) const { - assert(std::filesystem::is_directory(saveToPath)); - - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - std::filesystem::path simulationJobFolderPath( - std::filesystem::path(saveToPath) - .append(pFullPatternSimulationJob->label)); - std::filesystem::create_directory(simulationJobFolderPath); - const auto fullPatternDirectoryPath = - std::filesystem::path(simulationJobFolderPath).append("Full"); - std::filesystem::create_directory(fullPatternDirectoryPath); - pFullPatternSimulationJob->save(fullPatternDirectoryPath.string()); - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - const auto reducedPatternDirectoryPath = - std::filesystem::path(simulationJobFolderPath).append("Reduced"); - if (!std::filesystem::exists(reducedPatternDirectoryPath)) { - std::filesystem::create_directory(reducedPatternDirectoryPath); - } - pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string()); - } - } - void load(const string &loadFromPath) { - assert(std::filesystem::is_directory(loadFromPath)); - - for (const auto &directoryEntry : - filesystem::directory_iterator(loadFromPath)) { - const auto simulationScenarioPath = directoryEntry.path(); - if (!std::filesystem::is_directory(simulationScenarioPath)) { - continue; - } - // Load reduced pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath).append("Full"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - job.load(filepath.string()); - fullPatternSimulationJobs.push_back( - std::make_shared(job)); - } - } - - // Load full pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath) - .append("Reduced"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - job.load(filepath.string()); - reducedPatternSimulationJobs.push_back( - std::make_shared(job)); - } - } - } - } -#if POLYSCOPE_DEFINED - void draw() const { - initPolyscope(); - FormFinder simulator; - assert(fullPatternSimulationJobs.size() == - reducedPatternSimulationJobs.size()); - fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); - reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); - - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - // Drawing of full pattern results - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - pFullPatternSimulationJob->registerForDrawing( - fullPatternSimulationJobs[0]->pMesh->getLabel()); - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.registerForDrawing(); - // Drawing of reduced pattern results - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - reducedModelResults.registerForDrawing(); - polyscope::show(); - // Save a screensh - // const std::string screenshotFilename = - // "/home/iason/Coding/Projects/Approximating shapes with flat " - // "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" - // + m_pFullPatternSimulationMesh->getLabel() + "_" + - // simulationScenarioStrings[simulationScenarioIndex]; - // polyscope::screenshot(screenshotFilename, false); - fullModelResults.unregister(); - reducedModelResults.unregister(); - // double error = computeError( - // reducedModelResults, - // global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); - // std::cout << "Error of simulation scenario " - // << simulationScenarioStrings[simulationScenarioIndex] << " - // is " - // << error << std::endl; - } - } -#endif // POLYSCOPE_DEFINED - void saveMeshFiles() const { - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - assert(numberOfSimulationJobs != 0 && - fullPatternSimulationJobs.size() == - reducedPatternSimulationJobs.size()); - - fullPatternSimulationJobs[0]->pMesh->savePly("FullPattern_undeformed.ply"); - reducedPatternSimulationJobs[0]->pMesh->savePly( - "ReducedPattern_undeformed.ply"); - FormFinder simulator; - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - // Drawing of full pattern results - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.saveDeformedModel(); - - // Drawing of reduced pattern results - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationJobIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - reducedModelResults.saveDeformedModel(); - } - } - - void - writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization, - csvFile &os) { - os << "Obj value"; - for (const ReducedModelOptimizer::xRange &range : - settings_optimization.xRanges) { - os << range.label; - } - os << "Time(s)"; - os << "#Crashes"; - } - - void - writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization, - csvFile &os) const { - os << objectiveValue; - for (const double &optimalX : x) { - os << optimalX; - } - - for (int unusedXVarCounter = 0; - unusedXVarCounter < settings_optimization.xRanges.size() - x.size(); - unusedXVarCounter++) { - os << "-"; - } - - os << time; - if (numberOfSimulationCrashes == 0) { - os << "-"; - } else { - os << numberOfSimulationCrashes; - } - } -}; - #endif // REDUCEDMODELOPTIMIZER_HPP