From 50c15751a3b10b5d7ebb69052760c0135ba4a408 Mon Sep 17 00:00:00 2001 From: Iason Date: Wed, 24 Feb 2021 22:21:09 +0200 Subject: [PATCH] Normalizing using an estimate of the maximum displacement. Removed the R and D command line argument.Refactoring --- src/reducedmodeloptimizer.cpp | 327 +++++++++++++++++++++------------- 1 file changed, 200 insertions(+), 127 deletions(-) diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index d8623e6..0b70a03 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -9,7 +9,7 @@ struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; std::vector> fullPatternDisplacements; - std::vector fullPatternDisplacementNormSum; + std::vector reducedPatternMaximumDisplacementNormSum; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; std::unordered_map @@ -35,26 +35,49 @@ struct GlobalOptimizationVariables { ReducedModelOptimizer::Settings optimizationSettings; } global; +std::vector reducedPatternMaximumDisplacementSimulationJobs; + double ReducedModelOptimizer::computeError( const std::vector &reducedPatternDisplacements, const std::vector &fullPatternDisplacements, - const double &interfaceDisplacementsNormSum, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor) { + 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); + } + + return rawError; +} + +double ReducedModelOptimizer::computeRawError( + const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, const std::unordered_map &reducedToFullInterfaceViMap) { double error = 0; for (const auto reducedFullViPair : reducedToFullInterfaceViMap) { - VertexIndex reducedModelVi = reducedFullViPair.first; // 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; Eigen::Vector3d reducedVertexDisplacement( - reducedPatternDisplacements[reducedModelVi][0], - reducedPatternDisplacements[reducedModelVi][1], - reducedPatternDisplacements[reducedModelVi][2]); + reducedPatternDisplacements[reducedFullViPair.first][0], + reducedPatternDisplacements[reducedFullViPair.first][1], + reducedPatternDisplacements[reducedFullViPair.first][2]); if (!std::isfinite(reducedVertexDisplacement[0]) || !std::isfinite(reducedVertexDisplacement[1]) || !std::isfinite(reducedVertexDisplacement[2])) { + std::cout << "Displacements are not finite" << std::endl; std::terminate(); } Eigen::Vector3d fullVertexDisplacement( @@ -67,9 +90,6 @@ double ReducedModelOptimizer::computeError( error += errorVector.norm(); } - if (global.optimizationSettings.normalizeObjectiveValue) { - return error / std::max(interfaceDisplacementsNormSum, 0.00003); - } return error; } @@ -132,8 +152,8 @@ void updateMesh(long n, const double *x) { pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); } -double ReducedModelOptimizer::objective(double b, double h, double E) { - std::vector x{b, h, E}; +double ReducedModelOptimizer::objective(double b, double r, double E) { + std::vector x{b, r, E}; return ReducedModelOptimizer::objective(x.size(), x.data()); } @@ -146,12 +166,9 @@ double ReducedModelOptimizer::objective(double b, double h, double E, double ReducedModelOptimizer::objective(long n, const double *x) { // std::cout.precision(17); - // for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { - // std::cout << "x[" + std::to_string(parameterIndex) + "]=" - // << x[parameterIndex] << std::endl; - // } - const Element &e = global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; - // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // const Element &e = + // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << + // e.axialConstFactor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor // << std::endl; @@ -167,51 +184,79 @@ double ReducedModelOptimizer::objective(long n, const double *x) { FormFinder::Settings simulationSettings; // simulationSettings.shouldDraw = true; for (const int simulationScenarioIndex : global.simulationScenarioIndices) { + // std::cout << "Executing " + // << simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; SimulationResults reducedModelResults = simulator.executeSimulation( global.reducedPatternSimulationJobs[simulationScenarioIndex], simulationSettings); - std::string filename; + // std::string filename; if (!reducedModelResults.converged) { - error += std::numeric_limits::max(); global.numOfSimulationCrashes++; -if(beVerbose){ - std::cout << "Failed simulation" << std::endl; - } }else { - error += computeError( - reducedModelResults.displacements, - global.fullPatternDisplacements[simulationScenarioIndex], - global.fullPatternDisplacementNormSum[simulationScenarioIndex], - global.reducedToFullInterfaceViMap); - } - std::ofstream out(filename, std::ios_base::app); - auto pMesh = - global - .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] - ->pMesh; + 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; + } - 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(); #ifdef POLYSCOPE_DEFINED - ReducedModelOptimizer::visualizeResults( - global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs, - std::vector{ - static_cast(simulationScenarioIndex)}, - global.reducedToFullInterfaceViMap); + 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 + } + error += thisSimulationScenarioError; + } + // 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) { @@ -533,69 +578,50 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( #if POLYSCOPE_DEFINED void ReducedModelOptimizer::visualizeResults( - const std::vector> - &fullPatternSimulationJobs, - const std::vector> - &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, + const std::shared_ptr &pFullPatternSimulationJob, + const std::shared_ptr &pReducedPatternSimulationJob, const std::unordered_map - &reducedToFullInterfaceViMap) { + &reducedToFullInterfaceViMap, + const bool draw = true) { FormFinder simulator; std::shared_ptr pFullPatternSimulationMesh = - fullPatternSimulationJobs[0]->pMesh; + pFullPatternSimulationJob->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 interfaceDisplacementNormSum = 0; - for (const auto &interfaceViPair : reducedToFullInterfaceViMap) { - const int fullPatternInterfaceIndex = interfaceViPair.second; - Eigen::Vector3d fullPatternDisplacementVector( - fullModelResults.displacements[fullPatternInterfaceIndex][0], - fullModelResults.displacements[fullPatternInterfaceIndex][1], - fullModelResults.displacements[fullPatternInterfaceIndex][2]); - interfaceDisplacementNormSum += fullPatternDisplacementVector.norm(); - } - reducedModelResults.saveDeformedModel(); - fullModelResults.saveDeformedModel(); - double error = computeError( - reducedModelResults.displacements, fullModelResults.displacements, - interfaceDisplacementNormSum, reducedToFullInterfaceViMap); - std::cout << "Error of simulation scenario " - << simulationScenarioStrings[simulationScenarioIndex] << " is " - << error << std::endl; - totalError += error; - reducedModelResults.registerForDrawing(); - // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); - // registerWorldAxes(); + // 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 = computeRawError( + reducedModelResults.displacements, fullModelResults.displacements, + reducedToFullInterfaceViMap/*, + 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/build/OptimizationResults/" "Images/" + - pFullPatternSimulationMesh->getLabel() + "_" + - simulationScenarioStrings[simulationScenarioIndex]; - polyscope::show(); + pFullPatternSimulationMesh->getLabel() + "_" + "noScenarioName"; + // simulationScenarioStrings[simulationScenarioIndex]; polyscope::screenshot(screenshotFilename, false); fullModelResults.unregister(); reducedModelResults.unregister(); - // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); + pFullPatternSimulationMesh->unregister(); } - std::cout << "Total error:" << totalError << std::endl; } #endif // POLYSCOPE_DEFINED @@ -694,7 +720,19 @@ 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.optimizeInnerHexagonSize) { + results.objectiveValue = objective(4, results.x.data()); + + } else { + results.objectiveValue = + objective(results.x[0], results.x[1], results.x[2]); + } + global.optimizationSettings.shouldNormalizeObjectiveValue = true; + } results.time = elapsed.count() / 1000.0; + const bool printDebugInfo = false; if (printDebugInfo) { std::cout << "Finished optimizing." << endl; // std::cout << "Solution x:" << endl; @@ -1006,7 +1044,7 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.g_optimalReducedModelDisplacements.resize(6); global.reducedPatternSimulationJobs.resize(6); global.fullPatternDisplacements.resize(6); - global.fullPatternDisplacementNormSum.resize(6); + global.reducedPatternMaximumDisplacementNormSum.resize(6); global.g_firstRoundIterationIndex = 0; global.minY = std::numeric_limits::max(); global.numOfSimulationCrashes = 0; @@ -1014,34 +1052,18 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( global.optimizationSettings = optimizationSettings; global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh); + reducedPatternMaximumDisplacementSimulationJobs.resize(6); // polyscope::removeAllStructures(); - FormFinder::Settings settings; + FormFinder::Settings simulationSettings; // settings.shouldDraw = true; for (int simulationScenarioIndex : global.simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = global.fullPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob, settings); + SimulationResults fullModelResults = simulator.executeSimulation( + pFullPatternSimulationJob, simulationSettings); global.fullPatternDisplacements[simulationScenarioIndex] = fullModelResults.displacements; - double interfaceDisplacementNormSum = 0; - for (const auto &interfaceViPair : global.reducedToFullInterfaceViMap) { - const int fullPatternInterfaceIndex = interfaceViPair.second; - Eigen::Vector3d fullPatternDisplacementVector( - fullModelResults.displacements[fullPatternInterfaceIndex][0], - fullModelResults.displacements[fullPatternInterfaceIndex][1], - fullModelResults.displacements[fullPatternInterfaceIndex][2]); - interfaceDisplacementNormSum += fullPatternDisplacementVector.norm(); - } - global.fullPatternDisplacementNormSum[simulationScenarioIndex] = - interfaceDisplacementNormSum; - // global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize( - // m_pReducedPatternSimulationMesh->VN(), 3); - // computeDesiredReducedModelDisplacements( - // fullModelResults, global.reducedToFullInterfaceViMap, - // global.g_optimalReducedModelDisplacements[simulationScenarioIndex]); - SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; computeReducedModelSimulationJob(*pFullPatternSimulationJob, @@ -1049,6 +1071,54 @@ 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; + // if (global.optimizeInnerHexagonSize) { + // displacementMultiplier = 2; + // } else { + // displacementMultiplier = 2; + // } + global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex] = + displacementMultiplier * + std::max(errorOfMaxDisplacedReduced, errorOfNonDisplacedReduced); + // if (errorOfMaxDisplacedReduced > errorOfNonDisplacedReduced) { + // std::cout << "Max disp results in a bigger error for scenario " + // << simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; + // } else { + // std::cout << "Zero disp results in a bigger error for scenario " + // << simulationScenarioStrings[simulationScenarioIndex] + // << std::endl; + // } + } } Results optResults = runOptimization(optimizationSettings); for (int simulationScenarioIndex : global.simulationScenarioIndices) { @@ -1061,7 +1131,10 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( // optResults.draw(); // visualizeResults(simulationJobs, global.simulationScenarioIndices); - // visualizeResults(simulationJobs, global.reducedPatternSimulationJobs, - // global.simulationScenarioIndices,global.reducedToFullInterfaceViMap); +#ifdef POLYSCOPE_DEFINED +// visualizeResults( +// global.fullPatternSimulationJobs, global.reducedPatternSimulationJobs, +// global.simulationScenarioIndices, global.reducedToFullInterfaceViMap); +#endif // POLYSCOPE_DEFINED return optResults; }