diff --git a/src/main.cpp b/src/main.cpp index beead98..241fac9 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -38,20 +38,21 @@ int main(int argc, char *argv[]) { reducedPattern.scale(0.03, interfaceNodeIndex); // Set the optization settings - ReducedPatternOptimization::xRange beamE{"E", 0.001, 100}; + ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000}; ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000}; + ReducedPatternOptimization::xRange beamI{"I", 0.001, 1000}; ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000}; ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000}; ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000}; ReducedPatternOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95}; ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0}; ReducedPatternOptimization::Settings settings_optimization; - // settings_optimization.xRanges - // = {beamE, beamA, beamJ, beamI2, beamI3, innerHexagonSize, innerHexagonAngle}; - settings_optimization.xRanges = {beamE, beamA, beamI2, innerHexagonSize, innerHexagonAngle}; + settings_optimization.xRanges + = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle}; const bool input_numberOfFunctionCallsDefined = argc >= 4; - settings_optimization.numberOfFunctionCalls = - input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; + settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined + ? std::atoi(argv[3]) + : 100; settings_optimization.normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon; @@ -116,9 +117,7 @@ int main(int argc, char *argv[]) { const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]); ReducedModelOptimizer optimizer(numberOfNodesPerSlot); - optimizer.initializePatterns(fullPattern, - reducedPattern, - settings_optimization.xRanges.size()); + optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges); optimizer.optimize(settings_optimization, optimizationResults); optimizationResults.label = optimizationName; optimizationResults.baseTriangleFullPattern.copy(fullPattern); diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index d08078b..5fd16f4 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -20,7 +20,9 @@ struct GlobalOptimizationVariables { std::vector> fullPatternInterfaceViPairs; matplot::line_handle gPlotHandle; - std::vector objectiveValueHistory; + std::vector objectiveValueHistoryY; + std::vector objectiveValueHistoryX; + std::vector plotColors; Eigen::VectorXd initialParameters; std::vector simulationScenarioIndices; double minY{DBL_MAX}; @@ -40,6 +42,10 @@ struct GlobalOptimizationVariables { double desiredMaxDisplacementValue; double desiredMaxRotationAngle; std::string currentScenarioName; + std::vector &pReducedPatternSimulationMesh, + const double &newValue)>> + updateReducedPatternFunctions_material; + } global; double ReducedModelOptimizer::computeDisplacementError( @@ -165,7 +171,7 @@ double ReducedModelOptimizer::objective(const dlib::matrix &x) // polyscope::show(); // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); - // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " + // std::cout << e.axialConstFact|A,I2,I3,J,r,thetaor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor // << std::endl; @@ -249,28 +255,36 @@ double ReducedModelOptimizer::objective(const dlib::matrix &x) // std::cout << error << std::endl; if (totalError < global.minY) { global.minY = totalError; - global.minX.assign(x.begin(), x.begin() + n); + std::cout << "New best:" << totalError << std::endl; + // global.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; + global.objectiveValueHistoryY.push_back(std::log(totalError)); + global.objectiveValueHistoryX.push_back(global.numberOfFunctionCalls); + global.plotColors.push_back(0.1); + // auto xPlot = matplot::linspace(0, + // global.objectiveValueHistoryY.size(), + // global.objectiveValueHistoryY.size()); + global.gPlotHandle = matplot::scatter(global.objectiveValueHistoryX, + global.objectiveValueHistoryY, + 4, + global.plotColors); + // matplot::show(); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Objective value", + // global.objectiveValueHistoryY); } +// global.objectiveValueHistory.push_back(totalError); +// global.plotColors.push_back(10); #ifdef POLYSCOPE_DEFINED ++global.numberOfFunctionCalls; - global.objectiveValueHistory.push_back( - std::sqrt(totalError / global.simulationScenarioIndices.size())); if (global.optimizationSettings.numberOfFunctionCalls >= 100 && global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100) == 0) { std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl; - auto xPlot = matplot::linspace(0, - global.objectiveValueHistory.size(), - global.objectiveValueHistory.size()); - // std::vector colors(global.gObjectiveValueHistory.size(), 2); - // if (global.g_firstRoundIterationIndex != 0) { - // for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), - // [](double &c) { c = 0.7; }); - // } - global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory); - SimulationResultsReporter::createPlot("Number of Steps", - "Objective value", - global.objectiveValueHistory); } #endif // compute error and return it @@ -476,7 +490,7 @@ ReducedModelOptimizer::ReducedModelOptimizer(const std::vector &numberOf void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern, PatternGeometry &reducedPattern, - const int &optimizationParameters) + const std::vector &optimizationParameters) { assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN()); fullPatternNumberOfEdges = fullPattern.EN(); @@ -492,113 +506,12 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern, computeMaps(copyFullPattern, copyReducedPattern); createSimulationMeshes(copyFullPattern, copyReducedPattern); + initializeUpdateReducedPatternFunctions(); initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters); } -void ReducedModelOptimizer::initializeOptimizationParameters( - const std::shared_ptr &mesh, const int &optimizationParamters) +void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() { - global.numberOfOptimizationParameters = optimizationParamters; - global.initialParameters.resize(global.numberOfOptimizationParameters); - - switch (optimizationParamters) { - case 2: - function_updateReducedPattern_material = - [](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) {}; - break; - - case 3: - global.initialParameters(0) = mesh->elements[0].material.youngsModulus; - function_updateReducedPattern_material = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - const double E = global.initialParameters(0) * x(0); - for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); - } - }; - break; - case 4: - global.initialParameters(0) = mesh->elements[0].material.youngsModulus; - global.initialParameters(1) = mesh->elements[0].A; - function_updateReducedPattern_material = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - const double E = global.initialParameters(0) * x(0); - const double A = global.initialParameters(1) * x(1); - const double beamWidth = std::sqrt(A); - const double beamHeight = beamWidth; - for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); - e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); - } - }; - break; - case 5: - global.initialParameters(0) = mesh->elements[0].material.youngsModulus; - global.initialParameters(1) = mesh->elements[0].A; - global.initialParameters(2) = mesh->elements[0].inertia.I2; - function_updateReducedPattern_material = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - const double E = global.initialParameters(0) * x(0); - const double A = global.initialParameters(1) * x(1); - const double beamWidth = std::sqrt(A); - const double beamHeight = beamWidth; - const double I = global.initialParameters(2) * x(2); - for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); - e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); - e.inertia.I2 = I; - e.inertia.I3 = I; - } - }; - break; - case 7: - global.initialParameters(0) = mesh->elements[0].material.youngsModulus; - global.initialParameters(1) = mesh->elements[0].A; - global.initialParameters(2) = mesh->elements[0].inertia.J; - global.initialParameters(3) = mesh->elements[0].inertia.I2; - global.initialParameters(4) = mesh->elements[0].inertia.I3; - function_updateReducedPattern_material = [&](const dlib::matrix &x, - std::shared_ptr - &pReducedPatternSimulationMesh) { - const double E = global.initialParameters(0) * x(0); - const double A = global.initialParameters(1) * x(1); - const double beamWidth = std::sqrt(A); - const double beamHeight = beamWidth; - - const double J = global.initialParameters(2) * x(2); - const double I2 = global.initialParameters(3) * x(3); - const double I3 = global.initialParameters(4) * x(4); - for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); - e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); - e.inertia.J = J; - e.inertia.I2 = I2; - e.inertia.I3 = I3; - } - //#ifdef POLYSCOPE_DEFINED - // pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); - // pReducedPatternSimulationMesh->registerForDrawing(); - // std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2]) - // << std::endl; - // std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl; - // polyscope::show(); - //#endif - }; - break; - default: - std::cerr << "Wrong number of parameters:" << optimizationParamters << std::endl; - std::terminate(); - break; - } - function_updateReducedPattern_geometry = [&](const dlib::matrix &x, std::shared_ptr &pReducedPatternSimulationMesh) { @@ -630,6 +543,120 @@ void ReducedModelOptimizer::initializeOptimizationParameters( * baseTriangleMovableVertexPosition; } }; + + function_updateReducedPattern_material_E = + [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newE) { + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE)); + } + }; + function_updateReducedPattern_material_A = + [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newA) { + const double beamWidth = std::sqrt(newA); + const double beamHeight = beamWidth; + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); + } + }; + function_updateReducedPattern_material_I = + [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI) { + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.inertia.I2 = newI; + e.inertia.I3 = newI; + } + }; + function_updateReducedPattern_material_I2 = + [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI2) { + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.inertia.I2 = newI2; + } + }; + function_updateReducedPattern_material_I3 = + [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI3) { + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.inertia.I3 = newI3; + } + }; + function_updateReducedPattern_material_J = + [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newJ) { + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.inertia.J = newJ; + } + }; +} + +void ReducedModelOptimizer::initializeOptimizationParameters( + const std::shared_ptr &mesh, const std::vector &optimizationParameters) +{ + global.numberOfOptimizationParameters = optimizationParameters.size(); + global.initialParameters.resize(global.numberOfOptimizationParameters); + + for (int optimizationParameterIndex = 0; + optimizationParameterIndex < optimizationParameters.size(); + optimizationParameterIndex++) { + const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex]; + if (xOptimizationParameter.label == "E") { + global.initialParameters(optimizationParameterIndex) = mesh->elements[0] + .material.youngsModulus; + global.updateReducedPatternFunctions_material.push_back( + function_updateReducedPattern_material_E); + } else if (xOptimizationParameter.label == "A") { + global.initialParameters(optimizationParameterIndex) = mesh->elements[0].A; + global.updateReducedPatternFunctions_material.push_back( + function_updateReducedPattern_material_A); + } else if (xOptimizationParameter.label == "I") { + global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2; + global.updateReducedPatternFunctions_material.push_back( + function_updateReducedPattern_material_I); + } else if (xOptimizationParameter.label == "I2") { + global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2; + global.updateReducedPatternFunctions_material.push_back( + function_updateReducedPattern_material_I2); + } else if (xOptimizationParameter.label == "I3") { + global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I3; + global.updateReducedPatternFunctions_material.push_back( + function_updateReducedPattern_material_I3); + } else if (xOptimizationParameter.label == "J") { + global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.J; + global.updateReducedPatternFunctions_material.push_back( + function_updateReducedPattern_material_J); + } + //NOTE:Assuming that I2,I3 and J are be passed after E and A because the update of E and A changes I2,I3,J. + // case "HexSize": + // updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_E); + // break; + // case "HexAngle": + // updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_E); + // break; + } + + if (global.updateReducedPatternFunctions_material.size() != 0) { + function_updateReducedPattern_material = + [&](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) { + for (int optimizationParameterIndex = 0; + optimizationParameterIndex + < global.updateReducedPatternFunctions_material.size(); + optimizationParameterIndex++) { + global.updateReducedPatternFunctions_material[optimizationParameterIndex]( + pReducedPatternSimulationMesh, + global.initialParameters(optimizationParameterIndex) + * x(optimizationParameterIndex)); + } + }; + } else { + function_updateReducedPattern_material = + [](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) {}; + } } void ReducedModelOptimizer::computeReducedModelSimulationJob( @@ -795,16 +822,15 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza if (xVariableIndex < settings.xRanges.size() - 2) { results.optimalXNameValuePairs[xVariableIndex] = std::make_pair(settings.xRanges[xVariableIndex].label, - global.minX[xVariableIndex] + optimizationResult_dlib.x(xVariableIndex) * global.initialParameters(xVariableIndex)); } else { //Hex size and angle are pure values (not multipliers of the initial values) results.optimalXNameValuePairs[xVariableIndex] = std::make_pair(settings.xRanges[xVariableIndex].label, - global.minX[xVariableIndex]); + optimizationResult_dlib.x(xVariableIndex)); } - assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex)); optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex); #ifdef POLYSCOPE_DEFINED std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":" @@ -834,9 +860,9 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza global.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_total.resize( global.simulationScenarioIndices.size()); -#ifdef POLYSCOPE_DEFINED - global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); -#endif + //#ifdef POLYSCOPE_DEFINED + // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); + //#endif results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size()); for (int i = 0; i < global.simulationScenarioIndices.size(); i++) { @@ -1000,98 +1026,180 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces( void ReducedModelOptimizer::runOptimization(const Settings &settings, ReducedPatternOptimization::Results &results) { - global.objectiveValueHistory.clear(); +#if POLYSCOPE_DEFINED + global.objectiveValueHistoryY.clear(); + global.objectiveValueHistoryY.reserve(settings.numberOfFunctionCalls); + global.objectiveValueHistoryX.reserve(settings.numberOfFunctionCalls); + global.plotColors.reserve(settings.numberOfFunctionCalls); +#endif double (*objF)(const dlib::matrix &) = &objective; - //Geometry optimization of the reduced pattern - constexpr int numberOfGeometryOptimizationParameters = 2; - dlib::matrix xGeometryMin(numberOfGeometryOptimizationParameters); - dlib::matrix xGeometryMax(numberOfGeometryOptimizationParameters); - const int geometryParametersOffset - = global.numberOfOptimizationParameters - - numberOfGeometryOptimizationParameters; //I assume the geometry parameters come at the end - int paramIndex = 0; - for (int i = geometryParametersOffset; i < global.numberOfOptimizationParameters; - i++, paramIndex++) { - xGeometryMin(paramIndex) = settings.xRanges[i].min; - xGeometryMax(paramIndex) = settings.xRanges[i].max; - } - - //Set reduced pattern update functions to be used during optimization - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - pReducedPatternSimulationMesh->reset(); - }; - - dlib::function_evaluation result_dlib_geometry - = dlib::find_min_global(objF, - xGeometryMin, - xGeometryMax, - dlib::max_function_calls(settings.numberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), - settings.solverAccuracy); - //Material optimization of the reduced pattern + dlib::function_evaluation optimalResult; + constexpr bool shouldJointlyOptimizeGeometryAndMaterial = false; + const auto hexAngleParameterIt = std::find_if(settings.xRanges.begin(), + settings.xRanges.end(), + [](const xRange &x) { + return x.label == "HexAngle"; + }); + const bool hasHexAngleParameter = hexAngleParameterIt != settings.xRanges.end(); + const auto hexSizeParameterIt = std::find_if(settings.xRanges.begin(), + settings.xRanges.end(), + [](const xRange &x) { + return x.label == "HexSize"; + }); + const bool hasHexSizeParameter = hexSizeParameterIt != settings.xRanges.end(); + const bool hasBothGeometricalParameters = hasHexAngleParameter && hasHexSizeParameter; + assert(hasBothGeometricalParameters || (!hasHexAngleParameter && !hasHexSizeParameter)); + const int numberOfGeometryOptimizationParameters = hasBothGeometricalParameters ? 2 : 0; const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters - numberOfGeometryOptimizationParameters; - std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl; - std::cout << "opt size:" << global.minX[0] << std::endl; - dlib::function_evaluation result_dlib_material; - if (numberOfMaterialOptimizationParameters != 0) { - dlib::matrix xMaterialMin(numberOfMaterialOptimizationParameters); - dlib::matrix xMaterialMax(numberOfMaterialOptimizationParameters); - for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) { - xMaterialMin(i) = settings.xRanges[i].min; - xMaterialMax(i) = settings.xRanges[i].max; + const bool hasGeometryAndMaterialParameters = hasBothGeometricalParameters + && numberOfMaterialOptimizationParameters != 0; + if (!shouldJointlyOptimizeGeometryAndMaterial && hasGeometryAndMaterialParameters) { + //Geometry optimization of the reduced pattern + dlib::matrix xGeometryMin(numberOfGeometryOptimizationParameters); + dlib::matrix xGeometryMax(numberOfGeometryOptimizationParameters); + // const int hexAngleParameterIndex = std::distance(settings.xRanges.begin(), + // hexAngleParameterIt); + // const int hexSizeParameterIndex = std::distance(settings.xRanges.begin(), + // hexSizeParameterIt); + xGeometryMin(0) = hexSizeParameterIt->min; + xGeometryMax(0) = hexSizeParameterIt->max; + xGeometryMin(1) = hexAngleParameterIt->min; + xGeometryMax(1) = hexAngleParameterIt->max; + + //Set reduced pattern update functions to be used during optimization + function_updateReducedPattern = + [&](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) { + function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + pReducedPatternSimulationMesh->reset(); + }; + const int optimizationFunctionCalls_geometry + = global.numberOfOptimizationParameters == numberOfGeometryOptimizationParameters + ? settings.numberOfFunctionCalls + : static_cast(settings.numberOfFunctionCalls / 3.0); + dlib::function_evaluation result_dlib_geometry + = dlib::find_min_global(objF, + xGeometryMin, + xGeometryMax, + dlib::max_function_calls(optimizationFunctionCalls_geometry), + std::chrono::hours(24 * 365 * 290), + settings.solverAccuracy); + //Material optimization of the reduced pattern + // std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl; + // std::cout << "opt size:" << global.minX[0] << std::endl; + dlib::function_evaluation result_dlib_material; + if (numberOfMaterialOptimizationParameters != 0) { + dlib::matrix xMaterialMin(numberOfMaterialOptimizationParameters); + dlib::matrix xMaterialMax(numberOfMaterialOptimizationParameters); + for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) { + xMaterialMin(i) = settings.xRanges[i].min; + xMaterialMax(i) = settings.xRanges[i].max; + } + function_updateReducedPattern = + [&](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) { + function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + function_updateReducedPattern_geometry( + result_dlib_geometry.x, + pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this + pReducedPatternSimulationMesh->reset(); + }; + result_dlib_material = dlib::find_min_global(objF, + xMaterialMin, + xMaterialMax, + dlib::max_function_calls(static_cast( + 2.0 * settings.numberOfFunctionCalls + / 3.0)), + std::chrono::hours(24 * 365 * 290), + settings.solverAccuracy); } + // constexpr bool useBOBYQA = false; + // if (useBOBYQA) { + // const size_t npt = 2 * global.numberOfOptimizationParameters; + // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; + // const double rhobeg = 0.1; + // // const double rhobeg = 10; + // const double rhoend = rhobeg * 1e-6; + // // const size_t maxFun = 10 * (x.size() ^ 2); + // const size_t bobyqa_maxFunctionCalls = 200; + // dlib::find_min_bobyqa(objF, + // result_dlib.x, + // npt, + // xMaterialMin, + // xMaterialMax, + // rhobeg, + // rhoend, + // bobyqa_maxFunctionCalls); + // } + + optimalResult.x.set_size(global.numberOfOptimizationParameters); + std::copy(result_dlib_material.x.begin(), + result_dlib_material.x.end(), + optimalResult.x.begin()); + // debug_x.begin()); + std::copy(result_dlib_geometry.x.begin(), + result_dlib_geometry.x.end(), + optimalResult.x.begin() + numberOfMaterialOptimizationParameters); + std::cout << "opt x:"; + for (const auto optx : optimalResult.x) { + std::cout << optx << " "; + } + std::cout << std::endl; function_updateReducedPattern = [&](const dlib::matrix &x, std::shared_ptr &pReducedPatternSimulationMesh) { function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); pReducedPatternSimulationMesh->reset(); }; - result_dlib_material = dlib::find_min_global(objF, - xMaterialMin, - xMaterialMax, - dlib::max_function_calls( - settings.numberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), - settings.solverAccuracy); - } - // constexpr bool useBOBYQA = false; - // if (useBOBYQA) { - // const size_t npt = 2 * global.numberOfOptimizationParameters; - // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; - // const double rhobeg = 0.1; - // // const double rhobeg = 10; - // const double rhoend = rhobeg * 1e-6; - // // const size_t maxFun = 10 * (x.size() ^ 2); - // const size_t bobyqa_maxFunctionCalls = 200; - // dlib::find_min_bobyqa(objF, - // result_dlib.x, - // npt, - // xMaterialMin, - // xMaterialMax, - // rhobeg, - // rhoend, - // bobyqa_maxFunctionCalls); - // } + // const auto meshLabel = global.reducedPatternSimulationJobs[0]->pMesh->getLabel(); + // global.reducedPatternSimulationJobs[0]->pMesh->setLabel("pre"); + // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial); + // polyscope::show(); + // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); + optimalResult.y = objective(optimalResult.x); + } else { + dlib::matrix xMin(global.numberOfOptimizationParameters); + dlib::matrix xMax(global.numberOfOptimizationParameters); + for (int i = 0; i < global.numberOfOptimizationParameters; i++) { + xMin(i) = settings.xRanges[i].min; + xMax(i) = settings.xRanges[i].max; + } + if (numberOfGeometryOptimizationParameters != 0 + && numberOfMaterialOptimizationParameters != 0) { + function_updateReducedPattern = + [&](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) { + function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + pReducedPatternSimulationMesh->reset(); + }; + } else if (numberOfGeometryOptimizationParameters == 0) { + function_updateReducedPattern = + [&](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) { + function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + pReducedPatternSimulationMesh->reset(); + }; - dlib::function_evaluation optimalResult; - optimalResult.x.set_size(global.numberOfOptimizationParameters); - std::copy(result_dlib_material.x.begin(), result_dlib_material.x.end(), optimalResult.x.begin()); - std::copy(result_dlib_geometry.x.begin(), - result_dlib_geometry.x.end(), - optimalResult.x.begin() + geometryParametersOffset); - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - pReducedPatternSimulationMesh->reset(); - }; - optimalResult.y = objective(optimalResult.x); + } else { + function_updateReducedPattern = + [&](const dlib::matrix &x, + std::shared_ptr &pReducedPatternSimulationMesh) { + function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + pReducedPatternSimulationMesh->reset(); + }; + } + optimalResult = dlib::find_min_global(objF, + xMin, + xMax, + dlib::max_function_calls( + settings.numberOfFunctionCalls), + std::chrono::hours(24 * 365 * 290), + settings.solverAccuracy); + } getResults(optimalResult, settings, results); } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index b5d08fb..454d59f 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -67,9 +67,10 @@ public: SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); - void initializePatterns(PatternGeometry &fullPattern, - PatternGeometry &reducedPatterm, - const int &optimizationParameters); + void initializePatterns( + PatternGeometry &fullPattern, + PatternGeometry &reducedPatterm, + const std::vector &optimizationParameters); static void runSimulation(const std::string &filename, std::vector &x); @@ -183,9 +184,24 @@ public: static std::function &x, std::shared_ptr &m)> function_updateReducedPattern_material; - static std::function &x, - std::shared_ptr &m)> + static std::function &pReducedPatternSimulationMesh, + const double &newE)> function_updateReducedPattern_material_E; + static std::function &pReducedPatternSimulationMesh, + const double &newA)> + function_updateReducedPattern_material_A; + static std::function &pReducedPatternSimulationMesh, + const double &newI)> + function_updateReducedPattern_material_I; + static std::function &pReducedPatternSimulationMesh, + const double &newI2)> + function_updateReducedPattern_material_I2; + static std::function &pReducedPatternSimulationMesh, + const double &newI3)> + function_updateReducedPattern_material_I3; + static std::function &pReducedPatternSimulationMesh, + const double &newJ)> + function_updateReducedPattern_material_J; static std::function &x, std::shared_ptr &m)> function_updateReducedPattern_geometry; @@ -203,8 +219,9 @@ private: &maxForceMagnitudes); void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); - static void initializeOptimizationParameters(const std::shared_ptr &mesh, - const int &optimizationParamters); + static void initializeOptimizationParameters( + const std::shared_ptr &mesh, + const std::vector &optimizationParamters); static double objective(long n, const double *x); DRMSimulationModel simulator; @@ -224,11 +241,30 @@ private: const std::vector &desiredBaseSimulationScenarioIndices); static double objective(const dlib::matrix &x); + static void initializeUpdateReducedPatternFunctions(); }; inline std::function &x, std::shared_ptr &m)> ReducedModelOptimizer::function_updateReducedPattern; inline std::function &x, std::shared_ptr &m)> ReducedModelOptimizer::function_updateReducedPattern_material; +inline std::function &pReducedPatternSimulationMesh, + const double &newE)> + ReducedModelOptimizer::function_updateReducedPattern_material_E; +inline std::function &pReducedPatternSimulationMesh, + const double &newA)> + ReducedModelOptimizer::function_updateReducedPattern_material_A; +inline std::function &pReducedPatternSimulationMesh, + const double &newI)> + ReducedModelOptimizer::function_updateReducedPattern_material_I; +inline std::function &pReducedPatternSimulationMesh, + const double &newI2)> + ReducedModelOptimizer::function_updateReducedPattern_material_I2; +inline std::function &pReducedPatternSimulationMesh, + const double &newI3)> + ReducedModelOptimizer::function_updateReducedPattern_material_I3; +inline std::function &pReducedPatternSimulationMesh, + const double &newJ)> + ReducedModelOptimizer::function_updateReducedPattern_material_J; inline std::function &x, std::shared_ptr &m)> ReducedModelOptimizer::function_updateReducedPattern_geometry;