diff --git a/reducedmodeloptimizer.cpp b/reducedmodeloptimizer.cpp index 1ad60c3..7a3ac04 100644 --- a/reducedmodeloptimizer.cpp +++ b/reducedmodeloptimizer.cpp @@ -11,58 +11,28 @@ #include #include #include - //#define USE_SCENARIO_WEIGHTS using namespace ReducedModelOptimization; -struct GlobalOptimizationVariables +struct BaseScenarioMaxDisplacementHelperData { - std::vector fullPatternResults; - std::vector translationalDisplacementNormalizationValues; - std::vector rotationalDisplacementNormalizationValues; - std::vector> fullPatternSimulationJobs; - std::vector> reducedPatternSimulationJobs; - std::unordered_map reducedToFullInterfaceViMap; - std::vector> - fullPatternInterfaceViPairs; - matplot::line_handle gPlotHandle; - std::vector objectiveValueHistory_iteration; - std::vector objectiveValueHistory; - std::vector plotColors; - std::array - parametersInitialValue; - std::array - optimizationInitialValue; - std::vector simulationScenarioIndices; - double minY{DBL_MAX}; - std::vector minX; - int numOfSimulationCrashes{false}; - int numberOfFunctionCalls{0}; - ReducedModelOptimization::Settings optimizationSettings; - vcg::Triangle3 baseTriangle; - //Variables for finding the full pattern simulation forces - std::shared_ptr pFullPatternSimulationMesh; std::function> &, SimulationJob &)> constructScenarioFunction; - FullPatternVertexIndex interfaceViForComputingScenarioError; double desiredMaxDisplacementValue; double desiredMaxRotationAngle; + FullPatternVertexIndex interfaceViForComputingScenarioError; std::string currentScenarioName; - std::array &pReducedPatternSimulationMesh)>, - 7> - functions_updateReducedPatternParameter; - std::vector xMin; - std::vector xMax; - std::vector scenarioWeights; - std::vector objectiveWeights; -}; -GlobalOptimizationVariables global; + std::shared_ptr pFullPatternSimulationMesh; + + std::vector> + fullPatternOppositeInterfaceViPairs; +} g_baseScenarioMaxDisplacementHelperData; + +ReducedModelOptimizer::OptimizationState + g_optimizationState; //TODO: instead of having this global I could encapsulate it into a struct as I did for ensmallen double ReducedModelOptimizer::computeDisplacementError( const std::vector &fullPatternDisplacements, @@ -177,29 +147,34 @@ double ReducedModelOptimizer::computeError( //} #ifdef USE_ENSMALLEN -struct EnsmallenOptimizationObjective +struct EnsmallenReducedModelOptimizationObjective { - static double Evaluate(const arma::mat &x_arma) + EnsmallenReducedModelOptimizationObjective( + ReducedModelOptimizer::OptimizationState &optimizationState) + : optimizationState(optimizationState) + {} + ReducedModelOptimizer::OptimizationState &optimizationState; + 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]) { + for (int xi = 0; xi < x_arma.size(); xi++) { + if (x_arma[xi] > optimizationState.xMax[xi]) { //#ifdef POLYSCOPE_DEFINED // std::cout << "Out of range" << std::endl; //#endif return std::numeric_limits::max(); // return 1e10; - // x[xi] = global.xMax[xi]; - } else if (x[xi] < global.xMin[xi]) { + // x[xi] = optimizationState.xMax[xi]; + } else if (x_arma[xi] < optimizationState.xMin[xi]) { //#ifdef POLYSCOPE_DEFINED // std::cout << "Out of range" << std::endl; //#endif return std::numeric_limits::max(); // return 1e10; - // x[xi] = global.xMin[xi]; + // x[xi] = optimizationState.xMin[xi]; } } - return ReducedModelOptimizer::objective(x); + std::vector x(x_arma.begin(), x_arma.end()); + return ReducedModelOptimizer::objective(x, optimizationState); } }; #endif @@ -207,16 +182,17 @@ struct EnsmallenOptimizationObjective #ifdef DLIB_DEFINED double ReducedModelOptimizer::objective(const dlib::matrix &x) { - return objective(std::vector(x.begin(), x.end())); + return objective(std::vector(x.begin(), x.end()), g_optimizationState); } #endif -double ReducedModelOptimizer::objective(const double &xValue) -{ - return objective(std::vector({xValue})); -} +//double ReducedModelOptimizer::objective(const double &xValue) +//{ +// return objective(std::vector({xValue})); +//} -double ReducedModelOptimizer::objective(const std::vector &x) +double ReducedModelOptimizer::objective(const std::vector &x, + ReducedModelOptimizer::OptimizationState &optimizationState) { // std::cout.precision(17); // for (int i = 0; i < x.size(); i++) { @@ -226,19 +202,21 @@ double ReducedModelOptimizer::objective(const std::vector &x) // std::cout << x(x.size() - 2) << " " << x(x.size() - 1) << std::endl; // const Element &e = - // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << + // optimizationState.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << // e.axialConstFactor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor // << std::endl; // const int n = x.size(); std::shared_ptr &pReducedPatternSimulationMesh - = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; + = optimizationState + .reducedPatternSimulationJobs[optimizationState.simulationScenarioIndices[0]] + ->pMesh; function_updateReducedPattern(x, pReducedPatternSimulationMesh); - // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); - // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); + // optimizationState.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); + // optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); // polyscope::show(); - // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); + // optimizationState.reducedPatternSimulationJobs[0]->pMesh->unregister(); // run simulations std::vector simulationErrorsPerScenario(totalNumberOfSimulationScenarios, 0); @@ -249,12 +227,12 @@ double ReducedModelOptimizer::objective(const std::vector &x) std::for_each( std::execution::par_unseq, - global.simulationScenarioIndices.begin(), - global.simulationScenarioIndices.end(), + optimizationState.simulationScenarioIndices.begin(), + optimizationState.simulationScenarioIndices.end(), [&](const int &simulationScenarioIndex) { - // for (const int simulationScenarioIndex : global.simulationScenarioIndices) { + // for (const int simulationScenarioIndex : optimizationState.simulationScenarioIndices) { const std::shared_ptr &reducedJob - = global.reducedPatternSimulationJobs[simulationScenarioIndex]; + = optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]; //#ifdef POLYSCOPE_DEFINED // std::cout << reducedJob->getLabel() << ":" << std::endl; //#endif @@ -263,19 +241,20 @@ double ReducedModelOptimizer::objective(const std::vector &x) if (!reducedModelResults.converged) { simulationErrorsPerScenario[simulationScenarioIndex] = std::numeric_limits::max(); - global.numOfSimulationCrashes++; + optimizationState.numOfSimulationCrashes++; #ifdef POLYSCOPE_DEFINED std::cout << "Failed simulation" << std::endl; #endif } else { // const double simulationScenarioError_PE = std::abs( // (reducedModelResults.internalPotentialEnergy -// - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy) -// / global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); +// - optimizationState.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy) +// / optimizationState.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); // else { #ifdef POLYSCOPE_DEFINED #ifdef USE_SCENARIO_WEIGHTS - const double scenarioWeight = global.scenarioWeights[simulationScenarioIndex]; + const double scenarioWeight = optimizationState + .scenarioWeights[simulationScenarioIndex]; #else const double scenarioWeight = 1; #endif @@ -283,13 +262,15 @@ double ReducedModelOptimizer::objective(const std::vector &x) const double scenarioWeight = 1; #endif const double simulationScenarioError_displacements = computeError( - global.fullPatternResults[simulationScenarioIndex], + optimizationState.fullPatternResults[simulationScenarioIndex], reducedModelResults, - global.reducedToFullInterfaceViMap, - global.translationalDisplacementNormalizationValues[simulationScenarioIndex], - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + optimizationState.reducedToFullInterfaceViMap, + optimizationState + .translationalDisplacementNormalizationValues[simulationScenarioIndex], + optimizationState + .rotationalDisplacementNormalizationValues[simulationScenarioIndex], scenarioWeight, - global.objectiveWeights[simulationScenarioIndex]); + optimizationState.objectiveWeights[simulationScenarioIndex]); simulationErrorsPerScenario[simulationScenarioIndex] = simulationScenarioError_displacements /*+ simulationScenarioError_PE*/; @@ -297,60 +278,68 @@ double ReducedModelOptimizer::objective(const std::vector &x) // #ifdef POLYSCOPE_DEFINED // reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); // reducedModelResults.registerForDrawing(Colors::reducedDeformed); - // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); - // global.fullPatternResults[simulationScenarioIndex].registerForDrawing( + // optimizationState.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); + // optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing( // Colors::fullDeformed); // polyscope::show(); // reducedModelResults.unregister(); - // global.pFullPatternSimulationMesh->unregister(); - // global.fullPatternResults[simulationScenarioIndex].unregister(); + // optimizationState.pFullPatternSimulationMesh->unregister(); + // optimizationState.fullPatternResults[simulationScenarioIndex].unregister(); // #endif } }); - // double totalError = std::accumulate(simulationErrorsPerScenario.begin(), - // simulationErrorsPerScenario.end(), - // 0.0); const double totalError = std::reduce(std::execution::par_unseq, simulationErrorsPerScenario.begin(), simulationErrorsPerScenario.end()); // std::cout << totalError << std::endl; - // global.objectiveValueHistory.push_back(totalError); - // global.plotColors.push_back(10); - ++global.numberOfFunctionCalls; - if (totalError < global.minY) { - global.minY = totalError; - // global.objectiveValueHistory.push_back(totalError); - // global.objectiveValueHistory_iteration.push_back(global.numberOfFunctionCalls); + // optimizationState.objectiveValueHistory.push_back(totalError); + // optimizationState.plotColors.push_back(10); + ++optimizationState.numberOfFunctionCalls; + if (totalError < optimizationState.minY) { + optimizationState.minY = totalError; + // optimizationState.objectiveValueHistory.push_back(totalError); + // optimizationState.objectiveValueHistory_iteration.push_back(optimizationState.numberOfFunctionCalls); #ifdef POLYSCOPE_DEFINED std::cout << "New best:" << totalError << std::endl; - // // global.minX.assign(x.begin(), x.begin() + n); + // // optimizationState.minX.assign(x.begin(), x.begin() + n); std::cout.precision(17); for (int i = 0; i < x.size(); i++) { std::cout << x[i] << " "; } std::cout << std::endl; #endif - // global.objectiveValueHistoryY.push_back(std::log(totalError)); - // global.objectiveValueHistoryX.push_back(global.numberOfFunctionCalls); - // global.plotColors.push_back(0.1); + // optimizationState.objectiveValueHistoryY.push_back(std::log(totalError)); + // optimizationState.objectiveValueHistoryX.push_back(optimizationState.numberOfFunctionCalls); + // optimizationState.plotColors.push_back(0.1); // auto xPlot = matplot::linspace(0, - // global.objectiveValueHistoryY.size(), - // global.objectiveValueHistoryY.size()); - // global.gPlotHandle = matplot::scatter(global.objectiveValueHistoryX, - // global.objectiveValueHistoryY, + // optimizationState.objectiveValueHistoryY.size(), + // optimizationState.objectiveValueHistoryY.size()); + // optimizationState.gPlotHandle = matplot::scatter(optimizationState.objectiveValueHistoryX, + // optimizationState.objectiveValueHistoryY, // 4, - // global.plotColors); + // optimizationState.plotColors); // matplot::show(); // SimulationResultsReporter::createPlot("Number of Steps", // "Objective value", - // global.objectiveValueHistoryY); + // optimizationState.objectiveValueHistoryY); } #ifdef POLYSCOPE_DEFINED - if (global.optimizationSettings.numberOfFunctionCalls >= 100 - && global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100) - == 0) { - std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl; + +#ifdef USE_ENSMALLEN + if (optimizationState.numberOfFunctionCalls % 1000 == 0) { + std::cout << "Number of function calls:" << optimizationState.numberOfFunctionCalls + << std::endl; } +#else + if (optimizationState.optimizationSettings.dlib.numberOfFunctionCalls >= 100 + && optimizationState.numberOfFunctionCalls + % (optimizationState.optimizationSettings.dlib.numberOfFunctionCalls / 100) + == 0) { + std::cout << "Number of function calls:" << optimizationState.numberOfFunctionCalls + << std::endl; + } + +#endif #endif // compute error and return it return totalError; @@ -420,13 +409,13 @@ void ReducedModelOptimizer::computeMaps( // std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair << // std::endl; - // Save excluded edges TODO:this changes the global object. Should this be a + // Save excluded edges TODO:this changes the optimizationState object. Should this be a // function parameter? - // global.reducedPatternExludedEdges.clear(); + // optimizationState.reducedPatternExludedEdges.clear(); // const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN(); // for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { // for (const size_t ei : reducedModelExcludedEdges) { - // global.reducedPatternExludedEdges.insert( + // optimizationState.reducedPatternExludedEdges.insert( // fanIndex * reducedBaseTriangleNumberOfEdges + ei); // } // } @@ -440,7 +429,7 @@ void ReducedModelOptimizer::computeMaps( = reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/; Results::applyOptimizationResults_innerHexagon(initialHexagonSize, 0, - global.baseTriangle, + baseTriangle, reducedPattern); reducedPattern.createFan({0}); //TODO: should be an input parameter from main for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { @@ -462,7 +451,6 @@ void ReducedModelOptimizer::computeMaps( assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); fullPatternOppositeInterfaceViPairs.emplace_back(std::make_pair(vi0, vi1)); } - global.fullPatternInterfaceViPairs = fullPatternOppositeInterfaceViPairs; #if POLYSCOPE_DEFINED const bool debugMapping = false; @@ -471,7 +459,7 @@ void ReducedModelOptimizer::computeMaps( // std::vector colors_reducedPatternExcludedEdges( // reducedPattern.EN(), glm::vec3(0, 0, 0)); - // for (const size_t ei : global.reducedPatternExludedEdges) { + // for (const size_t ei : optimizationState.reducedPatternExludedEdges) { // colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); // } // const std::string label = reducedPattern.getLabel(); @@ -497,10 +485,11 @@ void ReducedModelOptimizer::computeMaps( glm::vec3(0, 0, 0)); std::vector nodeColorsReducedToFull_full(fullPattern.VN(), glm::vec3(0, 0, 0)); for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { - if (global.reducedToFullInterfaceViMap.contains(vi)) { + if (optimizationState.reducedToFullInterfaceViMap.contains(vi)) { auto color = polyscope::getNextUniqueColor(); nodeColorsReducedToFull_reduced[vi] = color; - nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] = color; + nodeColorsReducedToFull_full[optimizationState.reducedToFullInterfaceViMap[vi]] + = color; } } polyscope::getCurveNetwork(reducedPattern.getLabel()) @@ -520,7 +509,7 @@ void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern, ReducedModelOptimizer::computeMaps(slotToNode, fullPattern, reducedPattern, - global.reducedToFullInterfaceViMap, + optimizationState.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap, m_fullPatternOppositeInterfaceViPairs); } @@ -531,20 +520,23 @@ ReducedModelOptimizer::ReducedModelOptimizer() FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); interfaceNodeIndex = numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]; - constructBaseScenarioFunctions.resize(BaseSimulationScenario::NumberOfBaseSimulationScenarios); + // constructBaseScenarioFunctions.resize(BaseSimulationScenario::NumberOfBaseSimulationScenarios); scenarioIsSymmetrical.resize(BaseSimulationScenario::NumberOfBaseSimulationScenarios); constructBaseScenarioFunctions[BaseSimulationScenario::Axial] = &constructAxialSimulationScenario; + // constructBaseScenarioFunctions[BaseSimulationScenario::Axial] = &constructSSimulationScenario; scenarioIsSymmetrical[BaseSimulationScenario::Axial] = false; constructBaseScenarioFunctions[BaseSimulationScenario::Shear] = &constructShearSimulationScenario; scenarioIsSymmetrical[BaseSimulationScenario::Shear] = false; constructBaseScenarioFunctions[BaseSimulationScenario::Bending] = &constructBendingSimulationScenario; - scenarioIsSymmetrical[BaseSimulationScenario::Bending] = true; + scenarioIsSymmetrical[BaseSimulationScenario::Bending] = false; constructBaseScenarioFunctions[BaseSimulationScenario::Dome] = &constructDomeSimulationScenario; - scenarioIsSymmetrical[BaseSimulationScenario::Dome] = true; + scenarioIsSymmetrical[BaseSimulationScenario::Dome] = false; constructBaseScenarioFunctions[BaseSimulationScenario::Saddle] = &constructSaddleSimulationScenario; - scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = true; + scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = false; + constructBaseScenarioFunctions[BaseSimulationScenario::S] = &constructSSimulationScenario; + scenarioIsSymmetrical[BaseSimulationScenario::S] = false; } void ReducedModelOptimizer::initializePatterns( @@ -562,16 +554,19 @@ void ReducedModelOptimizer::initializePatterns( PatternGeometry copyReducedPattern; copyFullPattern.copy(fullPattern); copyReducedPattern.copy(reducedPattern); - global.baseTriangle = copyReducedPattern.getBaseTriangle(); + copyFullPattern.updateEigenEdgeAndVertices(); + copyReducedPattern.updateEigenEdgeAndVertices(); + baseTriangle = copyReducedPattern.getBaseTriangle(); computeMaps(copyFullPattern, copyReducedPattern); + optimizationState.fullPatternOppositeInterfaceViPairs = m_fullPatternOppositeInterfaceViPairs; createSimulationMeshes(copyFullPattern, copyReducedPattern); initializeUpdateReducedPatternFunctions(); initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters); } void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern, - const ReducedModelOptimization::Settings &optimizationSettings, + ReducedModelOptimization::Settings &optimizationSettings, ReducedModelOptimization::Results &optimizationResults) { //scale pattern and reduced model @@ -581,9 +576,10 @@ void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern, auto start = std::chrono::system_clock::now(); optimizationResults.label = fullPattern.getLabel(); optimizationResults.baseTriangleFullPattern.copy(fullPattern); - optimizationResults.settings = optimizationSettings; initializePatterns(fullPattern, reducedModel, optimizationSettings.variablesRanges); optimize(optimizationSettings, optimizationResults); + optimizationResults.settings + = optimizationSettings; //NOTE: currently I set the max force base magn in optimize thus this must be called after the optimize function auto end = std::chrono::system_clock::now(); auto elapsed = std::chrono::duration_cast(end - start); optimizationResults.time = elapsed.count() / 1000.0; @@ -591,15 +587,16 @@ void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern, void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() { - global.functions_updateReducedPatternParameter[R] = - [](const double &newR, std::shared_ptr &pReducedPatternSimulationMesh) { + const vcg::Triangle3 &baseTriangle = this->baseTriangle; + optimizationState.functions_updateReducedPatternParameter[R] = + [=](const double &newR, std::shared_ptr &pReducedPatternSimulationMesh) { const CoordType barycentricCoordinates_hexagonBaseTriangleVertex(1 - newR, newR / 2, newR / 2); const CoordType hexagonBaseTriangleVertexPosition - = global.baseTriangle.cP(0) * barycentricCoordinates_hexagonBaseTriangleVertex[0] - + global.baseTriangle.cP(1) * barycentricCoordinates_hexagonBaseTriangleVertex[1] - + global.baseTriangle.cP(2) * barycentricCoordinates_hexagonBaseTriangleVertex[2]; + = baseTriangle.cP(0) * barycentricCoordinates_hexagonBaseTriangleVertex[0] + + baseTriangle.cP(1) * barycentricCoordinates_hexagonBaseTriangleVertex[1] + + baseTriangle.cP(2) * barycentricCoordinates_hexagonBaseTriangleVertex[2]; for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { @@ -610,7 +607,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() } }; - global.functions_updateReducedPatternParameter[Theta] = + optimizationState.functions_updateReducedPatternParameter[Theta] = [](const double &newTheta, std::shared_ptr &pReducedPatternSimulationMesh) { const CoordType baseTriangleHexagonVertexPosition = pReducedPatternSimulationMesh->vert[0].cP(); @@ -628,14 +625,14 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() } }; - global.functions_updateReducedPatternParameter[E] = + optimizationState.functions_updateReducedPatternParameter[E] = [](const double &newE, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE)); } }; - global.functions_updateReducedPatternParameter[A] = + optimizationState.functions_updateReducedPatternParameter[A] = [](const double &newA, std::shared_ptr &pReducedPatternSimulationMesh) { assert(pReducedPatternSimulationMesh != nullptr); const double beamWidth = std::sqrt(newA); @@ -646,7 +643,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() } }; - global.functions_updateReducedPatternParameter[I2] = + optimizationState.functions_updateReducedPatternParameter[I2] = [](const double &newI2, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; @@ -654,7 +651,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() e.updateRigidity(); } }; - global.functions_updateReducedPatternParameter[I3] = + optimizationState.functions_updateReducedPatternParameter[I3] = [](const double &newI3, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; @@ -662,7 +659,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() e.updateRigidity(); } }; - global.functions_updateReducedPatternParameter[J] = + optimizationState.functions_updateReducedPatternParameter[J] = [](const double &newJ, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; @@ -686,37 +683,37 @@ void ReducedModelOptimizer::initializeOptimizationParameters( // const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex]; switch (optimizationParameterIndex) { case E: - global.parametersInitialValue[optimizationParameterIndex] + optimizationState.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].material.youngsModulus; - global.optimizationInitialValue[optimizationParameterIndex] = 1; + optimizationState.optimizationInitialValue[optimizationParameterIndex] = 1; break; case A: - global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] - .dimensions.A; - global.optimizationInitialValue[optimizationParameterIndex] = 1; + optimizationState.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].dimensions.A; + optimizationState.optimizationInitialValue[optimizationParameterIndex] = 1; break; case I2: - global.parametersInitialValue[optimizationParameterIndex] + optimizationState.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].dimensions.inertia.I2; - global.optimizationInitialValue[optimizationParameterIndex] = 1; + optimizationState.optimizationInitialValue[optimizationParameterIndex] = 1; break; case I3: - global.parametersInitialValue[optimizationParameterIndex] + optimizationState.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].dimensions.inertia.I3; - global.optimizationInitialValue[optimizationParameterIndex] = 1; + optimizationState.optimizationInitialValue[optimizationParameterIndex] = 1; break; case J: - global.parametersInitialValue[optimizationParameterIndex] + optimizationState.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].dimensions.inertia.J; - global.optimizationInitialValue[optimizationParameterIndex] = 1; + optimizationState.optimizationInitialValue[optimizationParameterIndex] = 1; break; case R: - global.parametersInitialValue[optimizationParameterIndex] = 0; - global.optimizationInitialValue[optimizationParameterIndex] = 0.5; + optimizationState.parametersInitialValue[optimizationParameterIndex] = 0; + optimizationState.optimizationInitialValue[optimizationParameterIndex] = 0.5; break; case Theta: - global.parametersInitialValue[optimizationParameterIndex] = 0; - global.optimizationInitialValue[optimizationParameterIndex] = 0; + optimizationState.parametersInitialValue[optimizationParameterIndex] = 0; + optimizationState.optimizationInitialValue[optimizationParameterIndex] = 0; break; } } @@ -777,15 +774,16 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv const Settings &settings, ReducedModelOptimization::Results &results) { - std::shared_ptr &pReducedPatternSimulationMesh - = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; + // std::shared_ptr &pReducedPatternSimulationMesh + // = g_optimizationState.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]] + // ->pMesh; //Number of crashes - // results.numberOfSimulationCrashes = global.numOfSimulationCrashes; + // results.numberOfSimulationCrashes = optimizationState.numOfSimulationCrashes; //Value of optimal objective Y results.objectiveValue.total = optimalObjective.y; // results.objectiveValue.total = 0; - if (optimalObjective.y != global.minY) { - std::cout << "Different optimal y:" << optimalObjective.y << " " << global.minY + if (std::abs(optimalObjective.y - optimizationState.minY) > 1e-1) { + std::cout << "Different optimal y:" << optimalObjective.y << " " << optimizationState.minY << std::endl; } //Optimal X values @@ -795,10 +793,10 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv optimizationParameterIndex++) { if (optimizationParameterIndex != OptimizationParameterIndex::R && optimizationParameterIndex != OptimizationParameterIndex::Theta) { - results.optimalXNameValuePairs[optimizationParameterIndex] - = std::make_pair(settings.variablesRanges[optimizationParameterIndex].label, - optimalObjective.x[optimizationParameterIndex] - * global.parametersInitialValue[optimizationParameterIndex]); + results.optimalXNameValuePairs[optimizationParameterIndex] = std::make_pair( + settings.variablesRanges[optimizationParameterIndex].label, + optimalObjective.x[optimizationParameterIndex] + * optimizationState.parametersInitialValue[optimizationParameterIndex]); } else { //Hex size and angle are pure values (not multipliers of the initial values) results.optimalXNameValuePairs[optimizationParameterIndex] @@ -806,9 +804,9 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv optimalObjective.x[optimizationParameterIndex]); } - global.functions_updateReducedPatternParameter[optimizationParameterIndex]( + optimizationState.functions_updateReducedPatternParameter[optimizationParameterIndex]( results.optimalXNameValuePairs[optimizationParameterIndex].second, - pReducedPatternSimulationMesh); + m_pReducedPatternSimulationMesh); //NOTE:due to this call this function is not const #ifdef POLYSCOPE_DEFINED std::cout << results.optimalXNameValuePairs[optimizationParameterIndex].first << ":" << results.optimalXNameValuePairs[optimizationParameterIndex].second << " "; @@ -817,7 +815,7 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv #ifdef POLYSCOPE_DEFINED std::cout << std::endl; #endif - pReducedPatternSimulationMesh->reset(); + m_pReducedPatternSimulationMesh->reset(); results.fullPatternYoungsModulus = youngsModulus; // Compute obj value per simulation scenario and the raw objective value @@ -827,30 +825,31 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv LinearSimulationModel simulator; results.objectiveValue.totalRaw = 0; results.objectiveValue.perSimulationScenario_translational.resize( - global.simulationScenarioIndices.size()); + optimizationState.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_rawTranslational.resize( - global.simulationScenarioIndices.size()); + optimizationState.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_rotational.resize( - global.simulationScenarioIndices.size()); + optimizationState.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_rawRotational.resize( - global.simulationScenarioIndices.size()); + optimizationState.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_total.resize( - global.simulationScenarioIndices.size()); + optimizationState.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_total_unweighted.resize( - global.simulationScenarioIndices.size()); + optimizationState.simulationScenarioIndices.size()); //#ifdef POLYSCOPE_DEFINED - // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); + // optimizationState.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); //#endif - results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size()); - for (int i = 0; i < global.simulationScenarioIndices.size(); i++) { - const int simulationScenarioIndex = global.simulationScenarioIndices[i]; + results.perScenario_fullPatternPotentialEnergy.resize( + optimizationState.simulationScenarioIndices.size()); + for (int i = 0; i < optimizationState.simulationScenarioIndices.size(); i++) { + const int simulationScenarioIndex = optimizationState.simulationScenarioIndices[i]; SimulationResults reducedModelResults = simulator.executeSimulation( - global.reducedPatternSimulationJobs[simulationScenarioIndex]); + optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]); #ifdef POLYSCOPE_DEFINED #ifdef USE_SCENARIO_WEIGHTS - const double scenarioWeight = global.scenarioWeights[simulationScenarioIndex]; + const double scenarioWeight = optimizationState.scenarioWeights[simulationScenarioIndex]; #else const double scenarioWeight = 1; #endif @@ -858,77 +857,83 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv const double scenarioWeight = 1; #endif results.objectiveValue.perSimulationScenario_total[i] = computeError( - global.fullPatternResults[simulationScenarioIndex], + optimizationState.fullPatternResults[simulationScenarioIndex], reducedModelResults, - global.reducedToFullInterfaceViMap, - global.translationalDisplacementNormalizationValues[simulationScenarioIndex], - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + optimizationState.reducedToFullInterfaceViMap, + optimizationState.translationalDisplacementNormalizationValues[simulationScenarioIndex], + optimizationState.rotationalDisplacementNormalizationValues[simulationScenarioIndex], scenarioWeight, - global.objectiveWeights[simulationScenarioIndex]); + optimizationState.objectiveWeights[simulationScenarioIndex]); results.objectiveValue.perSimulationScenario_total_unweighted[i] = computeError( - global.fullPatternResults[simulationScenarioIndex], + optimizationState.fullPatternResults[simulationScenarioIndex], reducedModelResults, - global.reducedToFullInterfaceViMap, - global.translationalDisplacementNormalizationValues[simulationScenarioIndex], - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + optimizationState.reducedToFullInterfaceViMap, + optimizationState.translationalDisplacementNormalizationValues[simulationScenarioIndex], + optimizationState.rotationalDisplacementNormalizationValues[simulationScenarioIndex], 1, - global.objectiveWeights[simulationScenarioIndex]); + optimizationState.objectiveWeights[simulationScenarioIndex]); // results.objectiveValue.total += results.objectiveValue.perSimulationScenario_total[i]; //Raw translational const double rawTranslationalError = computeRawTranslationalError( - global.fullPatternResults[simulationScenarioIndex].displacements, + optimizationState.fullPatternResults[simulationScenarioIndex].displacements, reducedModelResults.displacements, - global.reducedToFullInterfaceViMap); + optimizationState.reducedToFullInterfaceViMap); results.objectiveValue.perSimulationScenario_rawTranslational[i] = rawTranslationalError; //Raw rotational - const double rawRotationalError = computeRawRotationalError( - global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, - reducedModelResults.rotationalDisplacementQuaternion, - global.reducedToFullInterfaceViMap); + const double rawRotationalError + = computeRawRotationalError(optimizationState + .fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion, + reducedModelResults.rotationalDisplacementQuaternion, + optimizationState.reducedToFullInterfaceViMap); results.objectiveValue.perSimulationScenario_rawRotational[i] = rawRotationalError; //Normalized translational const double normalizedTranslationalError = computeDisplacementError( - global.fullPatternResults[simulationScenarioIndex].displacements, + optimizationState.fullPatternResults[simulationScenarioIndex].displacements, reducedModelResults.displacements, - global.reducedToFullInterfaceViMap, - global.translationalDisplacementNormalizationValues[simulationScenarioIndex]); + optimizationState.reducedToFullInterfaceViMap, + optimizationState.translationalDisplacementNormalizationValues[simulationScenarioIndex]); results.objectiveValue.perSimulationScenario_translational[i] = normalizedTranslationalError; // const double test_normalizedTranslationError = computeDisplacementError( - // global.fullPatternResults[simulationScenarioIndex].displacements, + // optimizationState.fullPatternResults[simulationScenarioIndex].displacements, // reducedModelResults.displacements, - // global.reducedToFullInterfaceViMap, - // global.translationalDisplacementNormalizationValues[simulationScenarioIndex]); + // optimizationState.reducedToFullInterfaceViMap, + // optimizationState.translationalDisplacementNormalizationValues[simulationScenarioIndex]); //Normalized rotational const double normalizedRotationalError = computeRotationalError( - global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, + optimizationState.fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion, reducedModelResults.rotationalDisplacementQuaternion, - global.reducedToFullInterfaceViMap, - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + optimizationState.reducedToFullInterfaceViMap, + optimizationState.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); results.objectiveValue.perSimulationScenario_rotational[i] = normalizedRotationalError; // const double test_normalizedRotationalError = computeRotationalError( - // global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, + // optimizationState.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, // reducedModelResults.rotationalDisplacementQuaternion, - // global.reducedToFullInterfaceViMap, - // global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + // optimizationState.reducedToFullInterfaceViMap, + // optimizationState.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); // assert(test_normalizedTranslationError == normalizedTranslationalError); // assert(test_normalizedRotationalError == normalizedRotationalError); results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError; results.perScenario_fullPatternPotentialEnergy[i] - = global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy; + = optimizationState.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy; #ifdef POLYSCOPE_DEFINED std::cout << "Simulation scenario:" - << global.reducedPatternSimulationJobs[simulationScenarioIndex]->getLabel() + << optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex] + ->getLabel() << std::endl; std::cout << "raw translational error:" << rawTranslationalError << std::endl; std::cout << "translation normalization value:" - << global.translationalDisplacementNormalizationValues[simulationScenarioIndex] + << optimizationState + .translationalDisplacementNormalizationValues[simulationScenarioIndex] << std::endl; std::cout << "raw Rotational error:" << rawRotationalError << std::endl; std::cout << "rotational normalization value:" - << global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] + << optimizationState + .rotationalDisplacementNormalizationValues[simulationScenarioIndex] << std::endl; std::cout << "Translational error:" << normalizedTranslationalError << std::endl; std::cout << "Rotational error:" << normalizedRotationalError << std::endl; @@ -939,27 +944,27 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv std::cout << std::endl; // reducedModelResults.registerForDrawing(Colors::reducedDeformed); - // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed); + // optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed); // polyscope::show(); // reducedModelResults.unregister(); - // global.fullPatternResults[simulationScenarioIndex].unregister(); + // optimizationState.fullPatternResults[simulationScenarioIndex].unregister(); #endif } - for (int simulationScenarioIndex : global.simulationScenarioIndices) { + for (int simulationScenarioIndex : optimizationState.simulationScenarioIndices) { results.fullPatternSimulationJobs.push_back( - global.fullPatternSimulationJobs[simulationScenarioIndex]); + optimizationState.fullPatternSimulationJobs[simulationScenarioIndex]); results.reducedPatternSimulationJobs.push_back( - global.reducedPatternSimulationJobs[simulationScenarioIndex]); - // const std::string temp = global.reducedPatternSimulationJobs[simulationScenarioIndex] + optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]); + // const std::string temp = optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex] // ->pMesh->getLabel(); - // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel("temp"); - // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(); - // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp); + // optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel("temp"); + // optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(); + // optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp); } - results.objectiveValueHistory = global.objectiveValueHistory; - results.objectiveValueHistory_iteration = global.objectiveValueHistory_iteration; + results.objectiveValueHistory = optimizationState.objectiveValueHistory; + results.objectiveValueHistory_iteration = optimizationState.objectiveValueHistory_iteration; // results.draw(); } @@ -967,7 +972,7 @@ std::array ReducedModelOptimizer::computeFullPatternMaxSimulationForces( const std::vector &desiredBaseSimulationScenario) const { - std::array fullPatternMaxSimulationForces; + std::array fullPatternMaxSimulationForces{0}; for (const BaseSimulationScenario &scenario : desiredBaseSimulationScenario) { const double maxForce = computeFullPatternMaxSimulationForce(scenario); fullPatternMaxSimulationForces[scenario] = maxForce; @@ -1027,13 +1032,18 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces( void ReducedModelOptimizer::runOptimization(const Settings &settings, ReducedModelOptimization::Results &results) { - global.objectiveValueHistory.clear(); - global.objectiveValueHistory_iteration.clear(); - global.objectiveValueHistory.reserve(settings.numberOfFunctionCalls / 100); - global.objectiveValueHistory_iteration.reserve(settings.numberOfFunctionCalls / 100); + optimizationState.objectiveValueHistory.clear(); + optimizationState.objectiveValueHistory_iteration.clear(); + optimizationState.objectiveValueHistory.reserve(settings.dlib.numberOfFunctionCalls / 100); + optimizationState.objectiveValueHistory_iteration.reserve(settings.dlib.numberOfFunctionCalls + / 100); + optimizationState.minY = std::numeric_limits::max(); + optimizationState.numOfSimulationCrashes = 0; + optimizationState.numberOfFunctionCalls = 0; + optimizationState.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh; #if POLYSCOPE_DEFINED -// global.plotColors.reserve(settings.numberOfFunctionCalls); +// optimizationState.plotColors.reserve(settings.numberOfFunctionCalls); #endif assert(!settings.optimizationStrategy.empty()); const std::vector> &optimizationParametersGroups @@ -1056,13 +1066,13 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, optimizationParameterIndex != NumberOfOptimizationVariables; optimizationParameterIndex++) { optimization_optimalResult.x[optimizationParameterIndex] - = global.optimizationInitialValue[optimizationParameterIndex]; + = optimizationState.optimizationInitialValue[optimizationParameterIndex]; } for (int parameterGroupIndex=0;parameterGroupIndex ¶meterGroup = optimizationParametersGroups[parameterGroupIndex]; FunctionEvaluation parameterGroup_optimalResult; - //Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup + //Set update function. TODO: Make this function immutable by defining it once and using the optimizationState variable to set parameterGroup function_updateReducedPattern = [&](const std::vector &x, std::shared_ptr &pMesh) { for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { @@ -1074,13 +1084,14 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, } //and the material parameters exponentially(?).TODO: Check what happens if I make all linear const double parameterInitialValue - = global.parametersInitialValue[parameterIndex]; + = optimizationState.parametersInitialValue[parameterIndex]; return x[xIndex] * parameterInitialValue; }(); // std::cout << "Optimization parameter:" << parameterIndex << std::endl; // std::cout << "New value:" << parameterNewValue << std::endl; - global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, - pMesh); + optimizationState + .functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, + pMesh); } pMesh->reset(); //NOTE: I could put this code into each updateParameter function for avoiding unessecary calculations }; @@ -1096,43 +1107,47 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, xMax[xIndex] = settings.variablesRanges[parameterIndex].max; } #ifdef USE_ENSMALLEN - arma::mat x(parameterGroup.size(), 1); - for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { - const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; +#ifdef POLYSCOPE_DEFINED + std::cout << "Optimizing using ensmallen" << std::endl; +#endif + 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]; - } + x(xIndex, 0) = optimizationState.optimizationInitialValue[parameterIndex]; + } - global.xMin = xMin; - global.xMax = xMax; + optimizationState.xMin = xMin; + optimizationState.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(100, + EnsmallenReducedModelOptimizationObjective optimizationFunction(optimizationState); +// optimizationFunction.optimizationState = optimizationState; +//Set min max values +// ens::CNE optimizer; +// ens::DE optimizer; +// ens::SPSA optimizer; +#ifdef USE_PSO + arma::mat xMin_arma(optimizationState.xMin); + arma::mat xMax_arma(optimizationState.xMax); + // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); + ens::LBestPSO optimizer(settings.pso.numberOfParticles, xMin_arma, xMax_arma, - settings.numberOfFunctionCalls, + 1e5, 500, settings.solverAccuracy, 2.05, 2.35); - // ens::LBestPSO optimizer; - const double minima = optimizer.Optimize(optimizationFunction, x); +#else + ens::SA optimizer; +#endif + // ens::LBestPSO optimizer; + parameterGroup_optimalResult.y = 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]; + // if (x[xIndex] > optimizationState.xMax[xIndex]) { + // x[xIndex] = optimizationState.xMax[xIndex]; + // } else if (x[xIndex] < optimizationState.xMin[xIndex]) { + // x[xIndex] = optimizationState.xMin[xIndex]; // } // } @@ -1146,102 +1161,108 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, parameterGroup_optimalResult.x.clear(); parameterGroup_optimalResult.x.resize(parameterGroup.size()); std::copy(x.begin(), x.end(), parameterGroup_optimalResult.x.begin()); - parameterGroup_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.variablesRanges[parameterIndex].min; - xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max; - } - const double portionOfTotalFunctionCallsBudget = [&]() { - if (settings.optimizationVariablesGroupsWeights.has_value()) { - assert(settings.optimizationVariablesGroupsWeights->size() - == settings.optimizationStrategy.size()); - return settings.optimizationVariablesGroupsWeights.value()[parameterGroupIndex]; - - } else { - return static_cast(parameterGroup.size()) - / totalNumberOfOptimizationParameters; - } - }(); - const int optimizationNumberOfFunctionCalls = portionOfTotalFunctionCallsBudget - * settings.numberOfFunctionCalls; - const dlib::function_evaluation optimalResult_dlib = [&]() { - if (parameterGroup.size() == 1) { - dlib::function_evaluation result; - double optimalX = global.optimizationInitialValue[parameterGroup[0]]; - double (*singleVariableObjective)(const double &) = &objective; - try { - result.y = dlib::find_min_single_variable(singleVariableObjective, - optimalX, - xMin_dlib(0), - xMax_dlib(0), - 1e-5, - optimizationNumberOfFunctionCalls); - } catch (const std::exception &e) { // reference to the base of a polymorphic object + g_optimizationState = optimizationState; #ifdef POLYSCOPE_DEFINED - std::cout << e.what() << std::endl; + std::cout << "Optimizing using dlib" << std::endl; #endif - } + //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]; - result.x.set_size(1); - result.x(0) = optimalX; - return result; + xMin_dlib(xIndex) = settings.variablesRanges[parameterIndex].min; + xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max; } - double (*objF)(const dlib::matrix &) = &objective; - return dlib::find_min_global(objF, - xMin_dlib, - xMax_dlib, - dlib::max_function_calls(optimizationNumberOfFunctionCalls), - 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]; + const double portionOfTotalFunctionCallsBudget = [&]() { + if (settings.optimizationVariablesGroupsWeights.has_value()) { + assert(settings.optimizationVariablesGroupsWeights->size() + == settings.optimizationStrategy.size()); + return settings.optimizationVariablesGroupsWeights.value()[parameterGroupIndex]; - // std::cout << "xIndex:" << xIndex << std::endl; - // std::cout << "xInit:" << optimalResult_dlib.x(xIndex) << std::endl; - // } + } else { + return static_cast(parameterGroup.size()) + / totalNumberOfOptimizationParameters; + } + }(); + const int optimizationNumberOfFunctionCalls = portionOfTotalFunctionCallsBudget + * settings.dlib.numberOfFunctionCalls; + const dlib::function_evaluation optimalResult_dlib = [&]() { + if (parameterGroup.size() == 1) { + dlib::function_evaluation result; + double optimalX = optimizationState.optimizationInitialValue[parameterGroup[0]]; + double (*singleVariableObjective)(const double &) = &objective; + try { + result.y = dlib::find_min_single_variable(singleVariableObjective, + optimalX, + xMin_dlib(0), + xMax_dlib(0), + 1e-5, + optimizationNumberOfFunctionCalls); + } catch ( + const std::exception &e) { // reference to the base of a polymorphic object +#ifdef POLYSCOPE_DEFINED + std::cout << e.what() << std::endl; +#endif + } - // optimalResult_dlib.y = dlib::find_min_bobyqa(objF, - // optimalResult_dlib.x, - // npt, - // xMin, - // xMax, - // rhobeg, - // rhoend, - // bobyqa_maxFunctionCalls); - // } + result.x.set_size(1); + result.x(0) = optimalX; + return result; + } + double (*objF)(const dlib::matrix &) = &objective; + return dlib::find_min_global(objF, + xMin_dlib, + xMax_dlib, + dlib::max_function_calls( + optimizationNumberOfFunctionCalls), + std::chrono::hours(24 * 365 * 290), + settings.solverAccuracy); + }(); + // constexpr bool useBOBYQA = false; + // if (useBOBYQA) { + // const size_t npt = 2 * optimizationState.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) = optimizationState.optimizationInitialValue[parameterIndex]; - parameterGroup_optimalResult.x.clear(); - parameterGroup_optimalResult.x.resize(parameterGroup.size()); - std::copy(optimalResult_dlib.x.begin(), - optimalResult_dlib.x.end(), - parameterGroup_optimalResult.x.begin()); - parameterGroup_optimalResult.y = optimalResult_dlib.y; + // 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); + // } + + parameterGroup_optimalResult.x.clear(); + parameterGroup_optimalResult.x.resize(parameterGroup.size()); + std::copy(optimalResult_dlib.x.begin(), + optimalResult_dlib.x.end(), + parameterGroup_optimalResult.x.begin()); + parameterGroup_optimalResult.y = optimalResult_dlib.y; + + optimizationState = g_optimizationState; #endif const auto firstNonNullReducedSimulationJobIt - = std::find_if(global.reducedPatternSimulationJobs.begin(), - global.reducedPatternSimulationJobs.end(), + = std::find_if(optimizationState.reducedPatternSimulationJobs.begin(), + optimizationState.reducedPatternSimulationJobs.end(), [](const std::shared_ptr &pJob) { return pJob != nullptr; }); @@ -1256,13 +1277,14 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, for (const auto &x : parameterGroup_optimalResult.x) { std::cout << x << std::endl; } - std::cout << "optimal objective:" << objective(parameterGroup_optimalResult.x) << std::endl; + std::cout << "optimal objective:" + << objective(parameterGroup_optimalResult.x, optimizationState) << std::endl; // std::cout << "Minima:" << minima << std::endl; - std::cout << "GLOBAL MIN:" << global.minY << std::endl; + std::cout << "optimizationState MIN:" << optimizationState.minY << std::endl; std::cout << "opt res y:" << parameterGroup_optimalResult.y << std::endl; #endif -optimization_optimalResult.y=parameterGroup_optimalResult.y; + optimization_optimalResult.y = parameterGroup_optimalResult.y; for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; optimization_optimalResult.x[parameterIndex] = parameterGroup_optimalResult.x[xIndex]; @@ -1290,6 +1312,28 @@ void ReducedModelOptimizer::constructAxialSimulationScenario( } } +void ReducedModelOptimizer::constructSSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job) +{ + for (auto viPairIt = oppositeInterfaceViPairs.begin(); + viPairIt != oppositeInterfaceViPairs.end(); + viPairIt++) { + if (viPairIt == oppositeInterfaceViPairs.begin()) { + job.constrainedVertices[viPairIt->first] = std::unordered_set{0, 1, 2, 3, 4, 5}; + job.constrainedVertices[viPairIt->second] = std::unordered_set{0, 1, 2, 3, 4, 5}; + } else { + CoordType forceDirection(0, 0, 1); + job.nodalExternalForces[viPairIt->first] = Vector6d({0, 0, 1, 0, 0, 0}) + * forceMagnitude; + job.nodalExternalForces[viPairIt->second] = Vector6d({0, 0, -1, 0, 0, 0}) + * forceMagnitude; + } + } +} + void ReducedModelOptimizer::constructShearSimulationScenario( const double &forceMagnitude, const std::vector> @@ -1403,10 +1447,13 @@ void ReducedModelOptimizer::constructSaddleSimulationScenario( double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude) { SimulationJob job; - job.pMesh = global.pFullPatternSimulationMesh; - global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job); + job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; + g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( + forceMagnitude, + g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, + job); // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, - // global.fullPatternInterfaceViPairs, + // optimizationState.fullPatternInterfaceViPairs, // job); DRMSimulationModel simulator; @@ -1421,7 +1468,9 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni settings.maxDRMIterations = 200000; SimulationResults results = simulator.executeSimulation(std::make_shared(job), settings); - const double &desiredRotationAngle = global.desiredMaxRotationAngle; + const double &desiredRotationAngle + = g_baseScenarioMaxDisplacementHelperData + .desiredMaxRotationAngle; //TODO: move from OptimizationState to a new struct double error; // job.pMesh->setLabel("initial"); // job.pMesh->registerForDrawing(); @@ -1449,14 +1498,16 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni saveJobToPath = "../nonConvergingJobs"; } else { error = std::abs( - results.rotationalDisplacementQuaternion[global.interfaceViForComputingScenarioError] + results + .rotationalDisplacementQuaternion[g_baseScenarioMaxDisplacementHelperData + .interfaceViForComputingScenarioError] .angularDistance(Eigen::Quaterniond::Identity()) - desiredRotationAngle); saveJobToPath = "../convergingJobs"; } // std::filesystem::path outputPath(std::filesystem::path(saveJobToPath) // .append(job.pMesh->getLabel()) - // .append("mag_" + global.currentScenarioName)); + // .append("mag_" + optimizationState.currentScenarioName)); // std::filesystem::create_directories(outputPath); // job.save(outputPath); // settings.save(outputPath); @@ -1498,10 +1549,13 @@ void search(const std::function objectiveFunction, double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude) { SimulationJob job; - job.pMesh = global.pFullPatternSimulationMesh; - global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job); + job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; + g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( + forceMagnitude, + g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, + job); // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, - // global.fullPatternInterfaceViPairs, + // optimizationState.fullPatternInterfaceViPairs, // job); DRMSimulationModel simulator; @@ -1523,21 +1577,27 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa settings.maxDRMIterations = 200000; SimulationResults results = simulator.executeSimulation(std::make_shared(job), settings); - const double &desiredDisplacementValue = global.desiredMaxDisplacementValue; + const double &desiredDisplacementValue + = g_baseScenarioMaxDisplacementHelperData + .desiredMaxDisplacementValue; //TODO: move from OptimizationState to a new struct double error; if (!results.converged) { error = std::numeric_limits::max(); - std::filesystem::path outputPath(std::filesystem::path("../nonConvergingJobs") - .append(job.pMesh->getLabel()) - .append("mag_" + global.currentScenarioName)); + std::filesystem::path outputPath( + std::filesystem::path("../nonConvergingJobs") + .append(job.pMesh->getLabel()) + .append("mag_" + g_baseScenarioMaxDisplacementHelperData.currentScenarioName)); std::filesystem::create_directories(outputPath); job.save(outputPath); settings.save(outputPath); // std::terminate(); } else { - error = std::abs( - results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm() - - desiredDisplacementValue); + error = std::abs(results + .displacements[g_baseScenarioMaxDisplacementHelperData + .interfaceViForComputingScenarioError] + .getTranslation() + .norm() + - desiredDisplacementValue); } #ifdef POLYSCOPE_DEFINED @@ -1562,41 +1622,58 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( double forceMagnitude = 1; double minimumError; bool wasSuccessful = false; - global.constructScenarioFunction = constructBaseScenarioFunctions[scenario]; - const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0), - (global.baseTriangle.cP(1) - + global.baseTriangle.cP(2)) - / 2); + g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction + = constructBaseScenarioFunctions[scenario]; + const double baseTriangleHeight = vcg::Distance(baseTriangle.cP(0), + (baseTriangle.cP(1) + baseTriangle.cP(2)) / 2); std::function objectiveFunction; double translationalOptimizationEpsilon{baseTriangleHeight * 0.001}; double objectiveEpsilon = translationalOptimizationEpsilon; objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective; - global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; - global.desiredMaxDisplacementValue = 0.1 * baseTriangleHeight; - global.currentScenarioName = baseSimulationScenarioNames[scenario]; + g_baseScenarioMaxDisplacementHelperData.interfaceViForComputingScenarioError + = m_fullPatternOppositeInterfaceViPairs[1].first; + g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = 0.1 * baseTriangleHeight; + g_baseScenarioMaxDisplacementHelperData.currentScenarioName + = baseSimulationScenarioNames[scenario]; + g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh; double forceMagnitudeEpsilon = 1e-4; switch (scenario) { case Axial: - global.desiredMaxDisplacementValue = 0.04 * baseTriangleHeight; + g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = 0.03 + * baseTriangleHeight; break; case Shear: - global.desiredMaxDisplacementValue = 0.04 * baseTriangleHeight; + g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = 0.04 + * baseTriangleHeight; break; case Bending: - global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; - global.desiredMaxDisplacementValue = 0.05 * baseTriangleHeight; + g_baseScenarioMaxDisplacementHelperData.interfaceViForComputingScenarioError + = m_fullPatternOppositeInterfaceViPairs[0].first; + g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = 0.05 + * baseTriangleHeight; break; case Dome: - global.desiredMaxRotationAngle = vcg::math::ToRad(20.0); + g_baseScenarioMaxDisplacementHelperData.desiredMaxRotationAngle = vcg::math::ToRad(20.0); objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective; forceMagnitudeEpsilon *= 1e-2; objectiveEpsilon = vcg::math::ToRad(1.0); forceMagnitude = 0.6; break; case Saddle: - global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; - global.desiredMaxDisplacementValue = 0.05 * baseTriangleHeight; + g_baseScenarioMaxDisplacementHelperData.interfaceViForComputingScenarioError + = m_fullPatternOppositeInterfaceViPairs[0].first; + g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = 0.05 + * baseTriangleHeight; + break; + case S: + g_baseScenarioMaxDisplacementHelperData.interfaceViForComputingScenarioError + = m_fullPatternOppositeInterfaceViPairs[1].first; + g_baseScenarioMaxDisplacementHelperData.desiredMaxDisplacementValue = 0.05 + * baseTriangleHeight; + break; + default: + std::cerr << "Unrecognized base scenario" << std::endl; break; } @@ -1621,7 +1698,7 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( // forceMagnitude = forceMagnitude_mat(0, 0); // search(objectiveFunction, -// global.desiredMaxDisplacementValue, +// optimizationState.desiredMaxDisplacementValue, // forceMagnitude, // minimumError, // objectiveEpsilon, @@ -1635,8 +1712,11 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( << " magnitude:" << forceMagnitude << std::endl; if (!wasSuccessful) { SimulationJob job; - job.pMesh = global.pFullPatternSimulationMesh; - global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job); + job.pMesh = optimizationState.pFullPatternSimulationMesh; + g_baseScenarioMaxDisplacementHelperData + .constructScenarioFunction(forceMagnitude, + optimizationState.fullPatternOppositeInterfaceViPairs, + job); std::cerr << ReducedModelOptimization::baseSimulationScenarioNames[scenario] + " max scenario magnitude was not succefully determined." << std::endl; @@ -1655,7 +1735,8 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( } void ReducedModelOptimizer::computeScenarioWeights( - const std::vector &baseSimulationScenarios) + const std::vector &baseSimulationScenarios, + const ReducedModelOptimization::Settings &optimizationSettings) { std::array baseScenario_equalizationWeights; baseScenario_equalizationWeights.fill(0); @@ -1692,8 +1773,8 @@ void ReducedModelOptimizer::computeScenarioWeights( } Utilities::normalize(baseScenarios_weights.begin(), baseScenarios_weights.end()); - global.objectiveWeights.resize(totalNumberOfSimulationScenarios); - global.scenarioWeights.resize( + optimizationState.objectiveWeights.resize(totalNumberOfSimulationScenarios); + optimizationState.scenarioWeights.resize( totalNumberOfSimulationScenarios, 0); //This essentially holds the base scenario weights but I use totalNumberOfSimulationScenarios elements instead in order to have O(1) access time during the optimization for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { @@ -1706,15 +1787,16 @@ void ReducedModelOptimizer::computeScenarioWeights( simulationScenarioIndex++) { const int globalScenarioIndex = baseSimulationScenarioIndexOffset + simulationScenarioIndex; - global.scenarioWeights[globalScenarioIndex] = baseScenarios_weights[baseScenario]; + optimizationState.scenarioWeights[globalScenarioIndex] + = baseScenarios_weights[baseScenario]; - global.objectiveWeights[globalScenarioIndex] - = global.optimizationSettings.perBaseScenarioObjectiveWeights[baseScenario]; + optimizationState.objectiveWeights[globalScenarioIndex] + = optimizationSettings.perBaseScenarioObjectiveWeights[baseScenario]; } } - //Dont normalize since we want the base scenarios to be normalized not the "global scenarios" - // Utilities::normalize(global.scenarioWeights.begin(), global.scenarioWeights.end()); + //Dont normalize since we want the base scenarios to be normalized not the "optimizationState scenarios" + // Utilities::normalize(optimizationState.scenarioWeights.begin(), optimizationState.scenarioWeights.end()); } std::vector> ReducedModelOptimizer::createFullPatternSimulationJobs( @@ -1765,12 +1847,13 @@ std::vector> ReducedModelOptimizer::createFullPat } #ifdef POLYSCOPE_DEFINED - std::cout << "Computed full pattern scenario magnitudes" << std::endl; + std::cout << "Created pattern simulation jobs" << std::endl; #endif return scenarios; } -void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() +void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors( + const ReducedModelOptimization::Settings &optimizationSettings) { // m_pFullPatternSimulationMesh->registerForDrawing(); // m_pFullPatternSimulationMesh->save(std::filesystem::current_path().append("initial.ply")); @@ -1778,15 +1861,15 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() std::vector fullPatternTranslationalDisplacementNormSum( totalNumberOfSimulationScenarios); std::vector fullPatternAngularDistance(totalNumberOfSimulationScenarios); - for (int simulationScenarioIndex : global.simulationScenarioIndices) { + for (int simulationScenarioIndex : optimizationState.simulationScenarioIndices) { double translationalDisplacementNormSum = 0; - for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + for (auto interfaceViPair : optimizationState.reducedToFullInterfaceViMap) { const int fullPatternVi = interfaceViPair.second; //If the full pattern vertex is translationally constrained dont take it into account - if (global.fullPatternSimulationJobs[simulationScenarioIndex] + if (optimizationState.fullPatternSimulationJobs[simulationScenarioIndex] ->constrainedVertices.contains(fullPatternVi)) { const std::unordered_set constrainedDof - = global.fullPatternSimulationJobs[simulationScenarioIndex] + = optimizationState.fullPatternSimulationJobs[simulationScenarioIndex] ->constrainedVertices.at(fullPatternVi); if (constrainedDof.contains(0) && constrainedDof.contains(1) && constrainedDof.contains(2)) { @@ -1794,35 +1877,37 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() } } - const Vector6d &vertexDisplacement = global.fullPatternResults[simulationScenarioIndex] + const Vector6d &vertexDisplacement = optimizationState + .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( + // optimizationState.fullPatternResults[simulationScenarioIndex].saveDeformedModel( // std::filesystem::current_path()); - // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(); + // optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing(); // polyscope::show(); - // global.fullPatternResults[simulationScenarioIndex].unregister(); + // optimizationState.fullPatternResults[simulationScenarioIndex].unregister(); double angularDistanceSum = 0; - for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + for (auto interfaceViPair : optimizationState.reducedToFullInterfaceViMap) { const int fullPatternVi = interfaceViPair.second; //If the full pattern vertex is rotationally constrained dont take it into account - if (global.fullPatternSimulationJobs[simulationScenarioIndex] + if (optimizationState.fullPatternSimulationJobs[simulationScenarioIndex] ->constrainedVertices.contains(fullPatternVi)) { const std::unordered_set constrainedDof - = global.fullPatternSimulationJobs[simulationScenarioIndex] + = optimizationState.fullPatternSimulationJobs[simulationScenarioIndex] ->constrainedVertices.at(fullPatternVi); if (constrainedDof.contains(3) && constrainedDof.contains(5) && constrainedDof.contains(4)) { continue; } } - angularDistanceSum += std::abs(global.fullPatternResults[simulationScenarioIndex] - .rotationalDisplacementQuaternion[fullPatternVi] - .angularDistance(Eigen::Quaterniond::Identity())); + angularDistanceSum += std::abs( + optimizationState.fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion[fullPatternVi] + .angularDistance(Eigen::Quaterniond::Identity())); } fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] @@ -1830,31 +1915,47 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum; } - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - if (global.optimizationSettings.normalizationStrategy - == Settings::NormalizationStrategy::Epsilon) { - const double epsilon_translationalDisplacement = global.optimizationSettings - .translationEpsilon; - global.translationalDisplacementNormalizationValues[simulationScenarioIndex] + optimizationState.translationalDisplacementNormalizationValues.resize( + totalNumberOfSimulationScenarios); + optimizationState.rotationalDisplacementNormalizationValues.resize( + totalNumberOfSimulationScenarios); + for (int simulationScenarioIndex : optimizationState.simulationScenarioIndices) { + if (optimizationSettings.normalizationStrategy == Settings::NormalizationStrategy::Epsilon) { + const double epsilon_translationalDisplacement = optimizationSettings.translationEpsilon; + optimizationState.translationalDisplacementNormalizationValues[simulationScenarioIndex] = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], epsilon_translationalDisplacement); - const double epsilon_rotationalDisplacement = global.optimizationSettings - .angularDistanceEpsilon; - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] + const double epsilon_rotationalDisplacement = optimizationSettings.angularDistanceEpsilon; + optimizationState.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = std::max(fullPatternAngularDistance[simulationScenarioIndex], epsilon_rotationalDisplacement); } else { - global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; + optimizationState.translationalDisplacementNormalizationValues[simulationScenarioIndex] + = 1; + optimizationState.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; } } } void ReducedModelOptimizer::optimize( - const Settings &optimizationSettings, + Settings &optimizationSettings, ReducedModelOptimization::Results &results, - const std::vector &desiredBaseSimulationScenarioIndices) + const std::vector &desiredBaseSimulationScenarios) { + std::unordered_set zeroMagnitudeScenarios; + for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios; baseScenario++) { + if (optimizationSettings.baseScenarioMaxMagnitudes[baseScenario] == 0) { + zeroMagnitudeScenarios.insert(static_cast(baseScenario)); + continue; + } + } + std::vector usedBaseSimulationScenarios; + std::set_difference(desiredBaseSimulationScenarios.begin(), + desiredBaseSimulationScenarios.end(), + zeroMagnitudeScenarios.begin(), + zeroMagnitudeScenarios.end(), + std::back_inserter(usedBaseSimulationScenarios)); + assert(!optimizationSettings.optimizationStrategy.empty()); assert(!optimizationSettings.optimizationVariablesGroupsWeights.has_value() || (optimizationSettings.optimizationStrategy.size() @@ -1863,58 +1964,47 @@ void ReducedModelOptimizer::optimize( optimizationSettings.optimizationVariablesGroupsWeights->end(), 0))); - for (int baseSimulationScenarioIndex : desiredBaseSimulationScenarioIndices) { + for (int baseSimulationScenarioIndex : usedBaseSimulationScenarios) { //Increase the size of the vector holding the simulation scenario indices - global.simulationScenarioIndices.resize( - global.simulationScenarioIndices.size() + optimizationState.simulationScenarioIndices.resize( + optimizationState.simulationScenarioIndices.size() + simulationScenariosResolution[baseSimulationScenarioIndex]); //Add the simulation scenarios indices that correspond to this base simulation scenario - std::iota(global.simulationScenarioIndices.end() + std::iota(optimizationState.simulationScenarioIndices.end() - simulationScenariosResolution[baseSimulationScenarioIndex], - global.simulationScenarioIndices.end(), + optimizationState.simulationScenarioIndices.end(), std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.begin() + baseSimulationScenarioIndex, 0)); } - // global.totalNumberOfSimulationScenarios = totalNumberOfSimulationScenarios; - global.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios); - global.fullPatternResults.resize(totalNumberOfSimulationScenarios); - global.translationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios); - global.rotationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios); - global.minY = std::numeric_limits::max(); - global.numOfSimulationCrashes = 0; - global.numberOfFunctionCalls = 0; - global.optimizationSettings = optimizationSettings; - global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh; - -#ifdef POLYSCOPE_DEFINED - constexpr bool recomputeMagnitudes = false; -#else - constexpr bool recomputeMagnitudes = true; -#endif - // std::array fullPatternSimulationScenarioMaxMagnitudes - // = getFullPatternMaxSimulationForces(desiredBaseSimulationScenarioIndices, - // intermediateResultsDirectoryPath, - // recomputeMagnitudes); + // constexpr bool recomputeForceMagnitudes = true; + // std::array fullPatternSimulationScenarioMaxMagnitudes + // = getFullPatternMaxSimulationForces(usedBaseSimulationScenarios, + // optimizationSettings.intermediateResultsDirectoryPath, + // recomputeForceMagnitudes); + // optimizationSettings.baseScenarioMaxMagnitudes = fullPatternSimulationScenarioMaxMagnitudes; std::array fullPatternSimulationScenarioMaxMagnitudes - = optimizationSettings.baseScenarioMagnitudes; - global.fullPatternSimulationJobs + = optimizationSettings.baseScenarioMaxMagnitudes; + optimizationState.fullPatternSimulationJobs = createFullPatternSimulationJobs(m_pFullPatternSimulationMesh, fullPatternSimulationScenarioMaxMagnitudes); // polyscope::removeAllStructures(); // ComputeScenarioWeights({Axial, Shear, Bending, Dome, Saddle}, - computeScenarioWeights(desiredBaseSimulationScenarioIndices); + computeScenarioWeights(usedBaseSimulationScenarios, optimizationSettings); - results.baseTriangle = global.baseTriangle; + results.baseTriangle = baseTriangle; DRMSimulationModel::Settings simulationSettings; simulationSettings.totalExternalForcesNormPercentageTermination = 1e-5; -// simulationSettings.maxDRMIterations = 200000; -// simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-10; -// simulationSettings.viscousDampingFactor = 5e-3; + simulationSettings.maxDRMIterations = 200000; + simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-15; +// simulationSettings.viscousDampingFactor = 5e-3; +// simulationSettings.debugModeStep = 100000; +// simulationSettings.shouldDraw = true; +// simulationSettings.beVerbose = true; // simulationSettings.useKineticDamping = true; // simulationSettings.save(std::filesystem::path(std::string(__FILE__)) // .parent_path() @@ -1932,31 +2022,40 @@ void ReducedModelOptimizer::optimize( #ifdef POLYSCOPE_DEFINED constexpr bool drawFullPatternSimulationResults = false; if (drawFullPatternSimulationResults) { - global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( + optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( ReducedModelOptimization::Colors::fullInitial); } #endif results.wasSuccessful = true; - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - const std::shared_ptr &pFullPatternSimulationJob - = global.fullPatternSimulationJobs[simulationScenarioIndex]; + optimizationState.fullPatternResults.resize(totalNumberOfSimulationScenarios); + optimizationState.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios); - std::filesystem::path jobResultsDirectoryPath( - std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath) - .append("FullPatternResults") - .append(m_pFullPatternSimulationMesh->getLabel()) - .append(pFullPatternSimulationJob->getLabel())); - // .append(pFullPatternSimulationJob->getLabel() + ".json") -#ifdef POLYSCOPE_DEFINED - constexpr bool recomputeFullPatternResults = false; -#else - constexpr bool recomputeFullPatternResults = true; -#endif - SimulationResults fullPatternResults; - if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) { - fullPatternResults.load(std::filesystem::path(jobResultsDirectoryPath).append("Results"), - pFullPatternSimulationJob); - } else { + std::for_each( + std::execution::par_unseq, + optimizationState.simulationScenarioIndices.begin(), + optimizationState.simulationScenarioIndices.end(), + [&](const int &simulationScenarioIndex) { + // for (int simulationScenarioIndex : optimizationState.simulationScenarioIndices) { + const std::shared_ptr &pFullPatternSimulationJob + = optimizationState.fullPatternSimulationJobs[simulationScenarioIndex]; + + std::filesystem::path jobResultsDirectoryPath( + std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath) + .append("FullPatternResults") + .append(m_pFullPatternSimulationMesh->getLabel()) + .append(pFullPatternSimulationJob->getLabel())); + // .append(pFullPatternSimulationJob->getLabel() + ".json") + //#ifdef POLYSCOPE_DEFINED + // constexpr bool recomputeFullPatternResults = true; + //#else + // constexpr bool recomputeFullPatternResults = true; + //#endif + SimulationResults fullPatternResults; + // if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) { + // fullPatternResults + // .load(std::filesystem::path(jobResultsDirectoryPath).append("Results"), + // pFullPatternSimulationJob); + // } else { // const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists( // patternMaxForceMagnitudesFilePath); // if (fullPatternScenarioMagnitudesExist) { @@ -1978,10 +2077,11 @@ void ReducedModelOptimizer::optimize( // pFullPatternSimulationJob); //#else + DRMSimulationModel simulator; fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, simulationSettings); - fullPatternResults.save(jobResultsDirectoryPath); - } +// fullPatternResults.save(jobResultsDirectoryPath); +// } //#endif // if (!fullPatternResults.converged) { // DRMSimulationModel::Settings simulationSettings_secondRound; @@ -1995,77 +2095,77 @@ void ReducedModelOptimizer::optimize( // << " used viscous damping." << std::endl; #endif - if (!fullPatternResults.converged) { - results.wasSuccessful = false; - std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel() - << " did not converge." << std::endl; + if (!fullPatternResults.converged) { + results.wasSuccessful = false; + std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel() + << " did not converge." << std::endl; #ifdef POLYSCOPE_DEFINED - // DRMSimulationModel::Settings debugSimulationSettings; - // debugSimulationSettings.debugModeStep = 50; - // // // debugSimulationSettings.maxDRMIterations = 100000; - // debugSimulationSettings.shouldDraw = true; - // // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; - // debugSimulationSettings.shouldCreatePlots = true; - // // // debugSimulationSettings.Dtini = 0.06; - // debugSimulationSettings.beVerbose = true; - // debugSimulationSettings.averageResidualForcesCriterionThreshold = 1e-5; - // debugSimulationSettings.maxDRMIterations = 100000; - // debugSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; - // debugSimulationSettings.viscousDampingFactor = 1e-2; - // // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; - // // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; - // auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob, - // debugSimulationSettings); - // debugResults.setLabelPrefix("debugResults"); - // debugResults.registerForDrawing(); - // polyscope::show(); - // debugResults.unregister(); - std::filesystem::path outputPath( - std::filesystem::path("../nonConvergingJobs") - .append(m_pFullPatternSimulationMesh->getLabel()) - .append("final_" + pFullPatternSimulationJob->getLabel())); - std::filesystem::create_directories(outputPath); - pFullPatternSimulationJob->save(outputPath); - simulationSettings.save(outputPath); - std::terminate(); - return; + // DRMSimulationModel::Settings debugSimulationSettings; + // debugSimulationSettings.debugModeStep = 50; + // // // debugSimulationSettings.maxDRMIterations = 100000; + // debugSimulationSettings.shouldDraw = true; + // // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; + // debugSimulationSettings.shouldCreatePlots = true; + // // // debugSimulationSettings.Dtini = 0.06; + // debugSimulationSettings.beVerbose = true; + // debugSimulationSettings.averageResidualForcesCriterionThreshold = 1e-5; + // debugSimulationSettings.maxDRMIterations = 100000; + // debugSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; + // debugSimulationSettings.viscousDampingFactor = 1e-2; + // // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; + // // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; + // auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob, + // debugSimulationSettings); + // debugResults.setLabelPrefix("debugResults"); + // debugResults.registerForDrawing(); + // polyscope::show(); + // debugResults.unregister(); + std::filesystem::path outputPath( + std::filesystem::path("../nonConvergingJobs") + .append(m_pFullPatternSimulationMesh->getLabel()) + .append("final_" + pFullPatternSimulationJob->getLabel())); + std::filesystem::create_directories(outputPath); + pFullPatternSimulationJob->save(outputPath); + simulationSettings.save(outputPath); + std::terminate(); + return; #endif - // } - } + // } + } #ifdef POLYSCOPE_DEFINED - if (drawFullPatternSimulationResults) { - // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( - // pFullPatternSimulationJob); - fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, - true); - // fullPatternResults_linear.labelPrefix += "_linear"; - // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, - // true, - // true); - polyscope::show(); - fullPatternResults.unregister(); - // fullPatternResults_linear.unregister(); - } + if (drawFullPatternSimulationResults) { + // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( + // pFullPatternSimulationJob); + fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + true); + // fullPatternResults_linear.labelPrefix += "_linear"; + // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + // true, + // true); + polyscope::show(); + fullPatternResults.unregister(); + // fullPatternResults_linear.unregister(); + } #endif - global.fullPatternResults[simulationScenarioIndex] = fullPatternResults; - SimulationJob reducedPatternSimulationJob; - reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; - computeReducedModelSimulationJob(*pFullPatternSimulationJob, - m_fullToReducedInterfaceViMap, - reducedPatternSimulationJob); - global.reducedPatternSimulationJobs[simulationScenarioIndex] - = std::make_shared(reducedPatternSimulationJob); - // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; - } + optimizationState.fullPatternResults[simulationScenarioIndex] = fullPatternResults; + SimulationJob reducedPatternSimulationJob; + reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; + computeReducedModelSimulationJob(*pFullPatternSimulationJob, + m_fullToReducedInterfaceViMap, + reducedPatternSimulationJob); + optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex] + = std::make_shared(reducedPatternSimulationJob); + // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; + }); #ifdef POLYSCOPE_DEFINED if (drawFullPatternSimulationResults) { - global.fullPatternSimulationJobs[0]->pMesh->unregister(); + optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); } #endif - // if (global.optimizationSettings.normalizationStrategy + // if (optimizationState.optimizationSettings.normalizationStrategy // == Settings::NormalizationStrategy::Epsilon) { - computeObjectiveValueNormalizationFactors(); + computeObjectiveValueNormalizationFactors(optimizationSettings); // } #ifdef POLYSCOPE_DEFINED std::cout << "Running reduced model optimization" << std::endl; diff --git a/reducedmodeloptimizer.hpp b/reducedmodeloptimizer.hpp index c9a2648..26644ca 100644 --- a/reducedmodeloptimizer.hpp +++ b/reducedmodeloptimizer.hpp @@ -25,6 +25,52 @@ using ReducedPatternVertexIndex = VertexIndex; class ReducedModelOptimizer { +public: + struct OptimizationState + { + std::vector fullPatternResults; + std::vector translationalDisplacementNormalizationValues; + std::vector rotationalDisplacementNormalizationValues; + std::vector> fullPatternSimulationJobs; + std::vector> reducedPatternSimulationJobs; + std::unordered_map + reducedToFullInterfaceViMap; + std::vector> + fullPatternOppositeInterfaceViPairs; + matplot::line_handle gPlotHandle; + std::vector objectiveValueHistory_iteration; + std::vector objectiveValueHistory; + std::vector plotColors; + std::array + parametersInitialValue; + std::array + optimizationInitialValue; + std::vector simulationScenarioIndices; + double minY{DBL_MAX}; + std::vector minX; + int numOfSimulationCrashes{false}; + int numberOfFunctionCalls{0}; + //Variables for finding the full pattern simulation forces + std::shared_ptr pFullPatternSimulationMesh; + std::array &pReducedPatternSimulationMesh)>, + 7> + functions_updateReducedPatternParameter; + std::vector xMin; + std::vector xMax; + std::vector scenarioWeights; + std::vector objectiveWeights; + }; + +private: + OptimizationState optimizationState; + vcg::Triangle3 baseTriangle; + std::function> &, + SimulationJob &)> + constructScenarioFunction; std::shared_ptr m_pReducedPatternSimulationMesh; std::shared_ptr m_pFullPatternSimulationMesh; std::unordered_map @@ -34,15 +80,17 @@ class ReducedModelOptimizer std::unordered_map nodeToSlot; std::unordered_map> slotToNode; std::string optimizationNotes; - std::vector> &, - SimulationJob &)>> + std::array> &, + SimulationJob &)>, + ReducedModelOptimization::NumberOfBaseSimulationScenarios> constructBaseScenarioFunctions; std::vector scenarioIsSymmetrical; int fullPatternNumberOfEdges; constexpr static double youngsModulus{1 * 1e9}; std::string fullPatternLabel; + // ReducedModelOptimization::Settings optimizationSettings; public: struct FunctionEvaluation @@ -69,7 +117,7 @@ public: inline static std::array parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"}; constexpr static std::array - simulationScenariosResolution = {11, 11, 20, 20, 20}; + simulationScenariosResolution = {11, 11, 21, 21, 21, 21}; constexpr static std::array baseScenarioWeights = {1, 1, 2, 3, 2}; inline static int totalNumberOfSimulationScenarios @@ -100,15 +148,15 @@ public: PatternGeometry &reducedModel, std::shared_ptr &pFullPatternSimulationMesh, std::shared_ptr &pReducedPatternSimulationMesh); - static void computeMaps(const std::unordered_map> &slotToNode, - PatternGeometry &fullPattern, - PatternGeometry &reducedPattern, - std::unordered_map - &reducedToFullInterfaceViMap, - std::unordered_map - &fullToReducedInterfaceViMap, - std::vector> - &fullPatternOppositeInterfaceViMap); + void computeMaps(const std::unordered_map> &slotToNode, + PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, + std::unordered_map + &reducedToFullInterfaceViMap, + std::unordered_map + &fullToReducedInterfaceViMap, + std::vector> + &fullPatternOppositeInterfaceViMap); static void visualizeResults( const std::vector> &fullPatternSimulationJobs, const std::vector> &reducedPatternSimulationJobs, @@ -134,10 +182,6 @@ public: &reducedToFullInterfaceViMap, const double &normalizationFactor); - static double objective(double E, - double A, - double innerHexagonSize, - double innerHexagonRotationAngle); static double computeRawRotationalError( const std::vector &rotatedQuaternion_fullPattern, const std::vector &rotatedQuaternion_reducedPattern, @@ -187,6 +231,11 @@ public: const std::vector> &oppositeInterfaceViPairs, SimulationJob &job); + static void constructSSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job); static std::function &x, std::shared_ptr &pReducedPatternSimulationMesh)> function_updateReducedPattern; @@ -209,20 +258,22 @@ public: std::shared_ptr &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_J; static double objective(const std::vector &x); - static void initializeUpdateReducedPatternFunctions(); - static double objective(const double &xValue); + void initializeUpdateReducedPatternFunctions(); + // static double objective(const double &xValue); ReducedModelOptimization::Results optimize( ConstPatternGeometry &fullPattern, - const ReducedModelOptimization::Settings &optimizationSettings); + const ReducedModelOptimization::Settings &optimizatifullyonSettings); void optimize(ConstPatternGeometry &fullPattern, - const ReducedModelOptimization::Settings &optimizationSettings, + ReducedModelOptimization::Settings &optimizationSettings, ReducedModelOptimization::Results &optimizationResults); + static double objective(const std::vector &x, + ReducedModelOptimizer::OptimizationState &optimizationState); private: void optimize( - const ReducedModelOptimization::Settings &xRanges, + ReducedModelOptimization::Settings &optimizationSettings, ReducedModelOptimization::Results &results, const std::vector &simulationScenarios = std::vector( @@ -230,7 +281,8 @@ private: ReducedModelOptimization::Shear, ReducedModelOptimization::Bending, ReducedModelOptimization::Dome, - ReducedModelOptimization::Saddle})); + ReducedModelOptimization::Saddle, + ReducedModelOptimization::S})); void initializePatterns(PatternGeometry &fullPattern, PatternGeometry &reducedPattern, @@ -241,30 +293,29 @@ private: const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - static void runOptimization(const ReducedModelOptimization::Settings &settings, - ReducedModelOptimization::Results &results); + void runOptimization(const ReducedModelOptimization::Settings &settings, + ReducedModelOptimization::Results &results); void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); - static void initializeOptimizationParameters( + void initializeOptimizationParameters( const std::shared_ptr &mesh, const std::array &optimizationParamters); DRMSimulationModel simulator; - void computeObjectiveValueNormalizationFactors(); + void computeObjectiveValueNormalizationFactors( + const ReducedModelOptimization::Settings &optimizationSettings); - static void getResults(const FunctionEvaluation &optimalObjective, - const ReducedModelOptimization::Settings &settings, - ReducedModelOptimization::Results &results); + void getResults(const FunctionEvaluation &optimalObjective, + const ReducedModelOptimization::Settings &settings, + ReducedModelOptimization::Results &results); double computeFullPatternMaxSimulationForce( const ReducedModelOptimization::BaseSimulationScenario &scenario) const; #ifdef DLIB_DEFINED static double objective(const dlib::matrix &x); #endif - void computeScenarioWeights(const std::vector - &baseSimulationScenarios); std::array computeFullPatternMaxSimulationForces( const std::vector @@ -281,6 +332,9 @@ private: const bool &recomputeForceMagnitudes); std::array getFullPatternMaxSimulationForces(); + void computeScenarioWeights( + const std::vector &baseSimulationScenarios, + const ReducedModelOptimization::Settings &optimizationSettings); }; inline std::function &pReducedPatternSimulationMesh)> @@ -303,7 +357,6 @@ inline std::function &x, std::shared_ptr &m)> ReducedModelOptimizer::function_updateReducedPattern; -struct GlobalOptimizationVariables; -extern GlobalOptimizationVariables global; +extern ReducedModelOptimizer::OptimizationState global; #endif // REDUCEDMODELOPTIMIZER_HPP