MySources/reducedmodeloptimizer_struc...

1214 lines
49 KiB
C++

#ifndef REDUCEDMODELOPTIMIZER_STRUCTS_HPP
#define REDUCEDMODELOPTIMIZER_STRUCTS_HPP
#include <string>
#include "chronoseulersimulationmodel.hpp"
#include "csvfile.hpp"
#include "drmsimulationmodel.hpp"
#include "linearsimulationmodel.hpp"
#include "matplot/matplot.h"
#include "simulation_structs.hpp"
#include "simulationmodelfactory.hpp"
#include "trianglepatterngeometry.hpp"
#include "unordered_map"
namespace ReducedModelOptimization {
enum BaseSimulationScenario {
Axial,
Shear,
Bending,
Dome,
Saddle,
S,
NumberOfBaseSimulationScenarios
};
inline static std::vector<std::string> baseSimulationScenarioNames{
"Axial", "Shear", "Bending", "Dome", "Saddle", "S"};
struct Colors {
using RGBColor = std::array<float, 3>;
inline static RGBColor patternInitial{0.518, 0.518, 0.518};
// inline static RGBColor fullDeformed{0.583333, 0.890196,
// 0.109804};
inline static RGBColor patternDeformed{0.094, 0.094, 0.094};
// inline static RGBColor reducedInitial{0.890196, 0.109804,
// 0.193138};
inline static RGBColor reducedInitial{0.518, 0.518, 0.518};
inline static RGBColor reducedDeformed{0.262, 0.627, 0.910};
};
struct xRange {
std::string label{};
double min{0};
double max{0};
NLOHMANN_DEFINE_TYPE_INTRUSIVE_WITH_DEFAULT(xRange, label, min, max)
inline bool operator<(const xRange& other) { return label < other.label; }
std::string toString() const {
return label + "=[" + std::to_string(min) + "," + std::to_string(max) + "]";
}
void fromString(const std::string& s) {
const std::size_t equalPos = s.find("=");
label = s.substr(0, equalPos);
const std::size_t commaPos = s.find(",");
const size_t minBeginPos = equalPos + 2;
min = std::stod(s.substr(minBeginPos, commaPos - minBeginPos));
const size_t maxBeginPos = commaPos + 1;
const std::size_t closingBracketPos = s.find("]");
max = std::stod(s.substr(maxBeginPos, closingBracketPos - maxBeginPos));
}
bool operator==(const xRange& xrange) const {
return label == xrange.label && min == xrange.min && max == xrange.max;
}
std::tuple<std::string, double, double> toTuple() const {
return std::make_tuple(label, min, max);
}
void set(const std::tuple<std::string, double, double>& inputTuple) {
if (std::get<1>(inputTuple) > std::get<2>(inputTuple)) {
std::cerr << "Invalid xRange tuple. Second argument(min) cant be smaller "
"than the third(max)"
<< std::endl;
std::terminate();
// return;
}
std::tie(label, min, max) = inputTuple;
}
};
enum OptimizationParameterIndex {
E,
A,
I2,
I3,
J,
R,
Theta,
NumberOfOptimizationVariables
};
inline int getParameterIndex(const std::string& s) {
if ("E" == s) {
return E;
} else if ("A" == s) {
return A;
} else if ("I2" == s) {
return I2;
} else if ("I3" == s) {
return I3;
} else if ("J" == s) {
return J;
} else if ("R" == s || "HexSize" == s) {
return R;
} else if ("Theta" == s || "HexAngle" == s) {
return Theta;
} else {
std::cerr
<< "Input is not recognized as a valid optimization variable index:"
<< s << std::endl;
return -1;
}
}
struct Settings {
inline static std::string defaultFilename{"OptimizationSettings.json"};
std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{
0.590241 / 3, 0.588372 / 3, 0.368304,
0.1, 1.18 / 2, 0}; // final b,h= 0.002
// std::array<double, NumberOfBaseSimulationScenarios>
// baseScenarioMaxMagnitudes{
// 0.590241, 0.888372, 0.368304, 0.0127508, 1.18079, 0 /*.3*/};
std::string getOptimizationSettingsLabel() const {
return simulationModelLabel_groundTruth + "_" +
simulationModelLabel_reducedModel;
}
enum OptimizationTool {
//#ifdef DLIB_DEFINED
dlib_globalMin = 0,
dlib_bobyqa,
//#else
ens_SA,
ens_DE,
ens_PSO,
//#endif
NumOptTools
};
OptimizationTool optimizationTool{ens_SA};
inline static const std::array<std::string, NumOptTools>
labels_optimizationTools{
//#ifdef DLIB_DEFINED
"dlib_globalMin", "dlib_bobyqa",
//#else
"ens_SA", "ens_DE", "ens_PSO"
//#endif
};
std::string getLabel_optimizationTool() const {
return labels_optimizationTools[optimizationTool];
}
std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
// {E, A, I2, I3, J, R, Theta}};
{A, I2, I3, J, R, Theta}};
// {A, I2, I3, J}};
// {R, Theta}};
// {R, Theta},
// {A, I2, I3, J}};
// {A, I2, I3, J},
// {R, Theta}};
std::optional<std::vector<double>>
optimizationVariablesGroupsWeights; // TODO:should be removed in the
// future if not a splitting
// strategy is used for the
// optimization
enum NormalizationStrategy { NonNormalized, Epsilon };
inline static std::vector<std::string> normalizationStrategyStrings{
"NonNormalized", "Epsilon"};
NormalizationStrategy normalizationStrategy{Epsilon};
std::array<xRange, NumberOfOptimizationVariables> variablesRanges{
xRange{"E", 1e-3, 1e3}, xRange{"A", 1e-3, 1e3}, xRange{"I2", 1e-3, 1e3},
xRange{"I3", 1e-3, 1e3}, xRange{"J", 1e-3, 1e3}, xRange{"R", 0.05, 0.95},
xRange{"Theta", -30, 30}};
inline static std::string simulationModelLabel_groundTruth{
// ChronosEulerNonLinearSimulationModel::label};
DRMSimulationModel::label};
inline static std::string simulationModelLabel_reducedModel{
LinearSimulationModel::label}; // namespace ReducedModelOptimization
// ChronosEulerLinearSimulationModel::label};
double targetBaseTriangleSize{0.03};
double translationEpsilon{4e-3};
// double translationEpsilon{targetBaseTriangleSize * 1e-2};
// double translationEpsilon{0};
double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
double solver_accuracy{1e-3};
double solver_maxIterations{100e3};
CrossSectionType beamDimensions_pattern{0.002, 0.002};
double youngsModulus_pattern{1e9};
std::filesystem::path intermediateResultsDirectoryPath;
bool beVerbose{false};
struct ObjectiveWeights {
double translational{1.2};
double rotational{0.8};
bool operator==(const ObjectiveWeights& other) const;
NLOHMANN_DEFINE_TYPE_INTRUSIVE(ObjectiveWeights, translational, rotational)
};
std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios>
perBaseScenarioObjectiveWeights;
NLOHMANN_DEFINE_TYPE_INTRUSIVE_WITH_DEFAULT(Settings,
optimizationTool,
baseScenarioMaxMagnitudes,
optimizationStrategy,
variablesRanges,
simulationModelLabel_groundTruth,
simulationModelLabel_reducedModel,
solver_accuracy,
solver_maxIterations,
translationEpsilon,
angularDistanceEpsilon,
targetBaseTriangleSize,
beamDimensions_pattern,
youngsModulus_pattern,
perBaseScenarioObjectiveWeights)
// std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios>
// perBaseScenarioObjectiveWeights{
// {{1.95, 0.05}, {0.87, 1.13}, {0.37, 1.63}, {0.01, 1.99},
// {0.94, 1.06}, {1.2, 0.8}}};
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
convertObjectiveWeightsToPairs() const {
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
objectiveWeightsPairs;
for (int baseScenario = Axial;
baseScenario != NumberOfBaseSimulationScenarios; baseScenario++) {
objectiveWeightsPairs[baseScenario] = std::make_pair(
perBaseScenarioObjectiveWeights[baseScenario].translational,
perBaseScenarioObjectiveWeights[baseScenario].rotational);
}
return objectiveWeightsPairs;
}
nlohmann::json toJson() const {
nlohmann::json json;
json = *this;
return json;
}
void save(const std::filesystem::path& saveToPath) {
assert(std::filesystem::is_directory(saveToPath));
nlohmann::json json = toJson();
std::filesystem::path jsonFilePath(
std::filesystem::path(saveToPath).append(defaultFilename));
std::ofstream jsonFile(jsonFilePath.string());
jsonFile << json;
jsonFile.close();
}
bool load(const std::filesystem::path& jsonFilePath) {
if (!std::filesystem::exists(jsonFilePath)) {
std::cerr << "Optimization settings could not be loaded because input "
"filepath does "
"not exist:"
<< jsonFilePath << std::endl;
assert(false);
return false;
}
std::ifstream ifs(jsonFilePath);
nlohmann::json json;
ifs >> json;
*this = json;
// if (json.contains(GET_VARIABLE_NAME(optimizationStrategy))) {
// optimizationStrategy = std::vector<
// std::vector<ReducedModelOptimization::OptimizationParameterIndex>>(
// (json.at(GET_VARIABLE_NAME(optimizationStrategy))));
// }
// if
// (json.contains(GET_VARIABLE_NAME(optimizationStrategyGroupWeights)))
// {
// optimizationVariablesGroupsWeights = std::vector<double>(
// json[GET_VARIABLE_NAME(optimizationStrategyGroupWeights)]);
// }
// // read x ranges
// if (json.contains(JsonKeys::OptimizationVariables)) {
// const std::array<std::tuple<std::string, double, double>,
// NumberOfOptimizationVariables>
// xRangesAsTuples = json.at(JsonKeys::OptimizationVariables);
// for (const auto& rangeTuple : xRangesAsTuples) {
// variablesRanges[getParameterIndex(std::get<0>(rangeTuple))].set(
// rangeTuple);
// }
// } else { // NOTE:legacy compatibility
// size_t xRangeIndex = 0;
// while (true) {
// const std::string jsonXRangeKey =
// JsonKeys::OptimizationVariables
// +
// "_" +
// std::to_string(xRangeIndex++);
// if (!json.contains(jsonXRangeKey)) {
// break;
// }
// xRange x;
// x.fromString(json.at(jsonXRangeKey));
// variablesRanges[getParameterIndex(x.label)] = x;
// }
// }
// if (json.contains(GET_VARIABLE_NAME(dlib.numberOfFunctionCalls))) {
// dlib.numberOfFunctionCalls =
// json.at(GET_VARIABLE_NAME(dlib.numberOfFunctionCalls));
// }
// if (json.contains(GET_VARIABLE_NAME(solverAccuracy))) {
// solverAccuracy = json.at(GET_VARIABLE_NAME(solverAccuracy));
// }
// // Objective weights
// if (json.contains(JsonKeys::ObjectiveWeights)) {
// std::array<std::pair<double, double>,
// NumberOfBaseSimulationScenarios>
// objectiveWeightsPairs = json.at(JsonKeys::ObjectiveWeights);
// std::transform(objectiveWeightsPairs.begin(),
// objectiveWeightsPairs.end(),
// perBaseScenarioObjectiveWeights.begin(),
// [](const std::pair<double, double>&
// objectiveWeightsPair) {
// return
// ObjectiveWeights{objectiveWeightsPair.first,
// objectiveWeightsPair.second};
// });
// }
// if
// (json.contains(GET_VARIABLE_NAME(translationalNormalizationEpsilon)))
// {
// translationEpsilon =
// json[GET_VARIABLE_NAME(translationalNormalizationEpsilon)];
// }
// if (json.contains(GET_VARIABLE_NAME(angularDistanceEpsilon))) {
// angularDistanceEpsilon = vcg::math::ToRad(
// static_cast<double>(json[GET_VARIABLE_NAME(angularDistanceEpsilon)]));
// }
// if (json.contains(GET_VARIABLE_NAME(targetBaseTriangleSize))) {
// targetBaseTriangleSize =
// static_cast<double>(json[GET_VARIABLE_NAME(targetBaseTriangleSize)]);
// }
// if (json.contains(GET_VARIABLE_NAME(pso.numberOfParticles))) {
// pso.numberOfParticles =
// static_cast<int>(json[GET_VARIABLE_NAME(pso.numberOfParticles)]);
// }
// if
// (json.contains(GET_VARIABLE_NAME(simulationModelLabel_reducedModel)))
// {
// simulationModelLabel_reducedModel = static_cast<std::string>(
// json[GET_VARIABLE_NAME(simulationModelLabel_reducedModel)]);
// }
// beamDimensions_pattern.from_json(json, beamDimensions_pattern);
// if (json.contains(GET_VARIABLE_NAME(youngsModulus_pattern))) {
// youngsModulus_pattern =
// static_cast<double>(json[GET_VARIABLE_NAME(youngsModulus_pattern)]);
// }
// perBaseScenarioObjectiveWeights =
// json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.translational =
// json.at(JsonKeys::ObjectiveWeights); objectiveWeights.rotational
// = 2 - objectiveWeights.translational;
return true;
}
std::string toString() const { return toJson().dump(); }
void writeHeaderTo(CSVFile& os) const {
if (!variablesRanges.empty()) {
for (const xRange& range : variablesRanges) {
os << range.label + " max";
os << range.label + " min";
}
}
os << "Function Calls";
os << "Solution Accuracy";
os << "Normalization trans epsilon";
os << "Normalization rot epsilon(deg)";
os << GET_VARIABLE_NAME(perBaseScenarioObjectiveWeights);
os << "Optimization parameters";
// os << std::endl;
}
void writeSettingsTo(CSVFile& os) const {
if (!variablesRanges.empty()) {
for (const xRange& range : variablesRanges) {
os << range.max;
os << range.min;
}
}
os << solver_accuracy;
os << std::to_string(translationEpsilon);
os << std::to_string(vcg::math::ToDeg(angularDistanceEpsilon));
std::string objectiveWeightsString;
objectiveWeightsString += "{";
for (int baseScenario = Axial;
baseScenario != NumberOfBaseSimulationScenarios; baseScenario++) {
objectiveWeightsString +=
"{" +
std::to_string(
perBaseScenarioObjectiveWeights[baseScenario].translational) +
"," +
std::to_string(
perBaseScenarioObjectiveWeights[baseScenario].rotational) +
"}";
}
objectiveWeightsString += "}";
os << objectiveWeightsString;
// export optimization parameters
std::vector<std::vector<int>> vv;
for (const std::vector<OptimizationParameterIndex>& v :
optimizationStrategy) {
std::vector<int> vi;
vi.reserve(v.size());
for (const OptimizationParameterIndex& parameter : v) {
vi.emplace_back(parameter);
}
vv.push_back(vi);
}
os << Utilities::toString(vv);
}
}; // namespace ReducedModelOptimization
inline bool operator==(const Settings& settings1,
const Settings& settings2) noexcept {
const bool haveTheSameObjectiveWeights =
std::mismatch(settings1.perBaseScenarioObjectiveWeights.begin(),
settings1.perBaseScenarioObjectiveWeights.end(),
settings2.perBaseScenarioObjectiveWeights.begin())
.first == settings1.perBaseScenarioObjectiveWeights.end();
return settings1.variablesRanges == settings2.variablesRanges &&
settings1.solver_accuracy == settings2.solver_accuracy &&
haveTheSameObjectiveWeights &&
settings1.translationEpsilon == settings2.translationEpsilon;
}
inline void updateMeshWithOptimalVariables(const std::vector<double>& x,
SimulationEdgeMesh& m) {
assert(CrossSectionType::name == RectangularBeamDimensions::name);
const double E = x[0];
const double A = x[1];
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
const double J = x[2];
const double I2 = x[3];
const double I3 = x[4];
for (EdgeIndex ei = 0; ei < m.EN(); ei++) {
Element& e = m.elements[ei];
e.setDimensions(CrossSectionType(beamWidth, beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
e.dimensions.inertia.J = J;
e.dimensions.inertia.I2 = I2;
e.dimensions.inertia.I3 = I3;
}
CoordType center_barycentric(1, 0, 0);
CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
CoordType movableVertex_barycentric(
(center_barycentric + interfaceEdgeMiddle_barycentric) * x[x.size() - 2]);
CoordType patternCoord0 = CoordType(0, 0, 0);
double bottomEdgeHalfSize = 0.03 / std::tan(M_PI / 3);
CoordType interfaceNodePosition(0, -0.03, 0);
CoordType patternBottomRight =
interfaceNodePosition + CoordType(bottomEdgeHalfSize, 0, 0);
CoordType patternBottomLeft =
interfaceNodePosition - CoordType(bottomEdgeHalfSize, 0, 0);
vcg::Triangle3<double> baseTriangle(patternCoord0, patternBottomLeft,
patternBottomRight);
CoordType baseTriangleMovableVertexPosition =
baseTriangle.cP(0) * movableVertex_barycentric[0] +
baseTriangle.cP(1) * movableVertex_barycentric[1] +
baseTriangle.cP(2) * movableVertex_barycentric[2];
VectorType patternPlaneNormal(0, 0, 1);
baseTriangleMovableVertexPosition =
vcg::RotationMatrix(patternPlaneNormal,
vcg::math::ToRad(x[x.size() - 1])) *
baseTriangleMovableVertexPosition;
const int fanSize = 6;
for (int rotationCounter = 0; rotationCounter < fanSize; rotationCounter++) {
m.vert[2 * rotationCounter + 1].P() =
vcg::RotationMatrix(patternPlaneNormal,
vcg::math::ToRad(60.0 * rotationCounter)) *
baseTriangleMovableVertexPosition;
}
m.reset();
}
struct Results {
std::string label{"EmptyLabel"};
double time{-1};
bool wasSuccessful{true};
// int numberOfSimulationCrashes{0};
Settings settings;
std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
std::vector<std::shared_ptr<SimulationJob>> pSimulationJobs_pattern;
std::vector<std::shared_ptr<SimulationJob>> pSimulationJobs_reducedModel;
// Full pattern
PatternGeometry baseTrianglePattern; // non-fanned,non-tiled full pattern
vcg::Triangle3<double> baseTriangle;
// std::string notes;
// Data gathered for csv exporting
struct ObjectiveValues {
double totalRaw;
double total;
std::vector<double> perSimulationScenario_rawTranslational;
std::vector<double> perSimulationScenario_rawRotational;
std::vector<double> perSimulationScenario_translational;
std::vector<double> perSimulationScenario_rotational;
std::vector<double> perSimulationScenario_total;
std::vector<double> perSimulationScenario_total_unweighted;
} objectiveValue;
std::vector<double> perScenario_fullPatternPotentialEnergy;
std::vector<double> objectiveValueHistory;
std::vector<size_t> objectiveValueHistory_iteration;
inline static std::string DefaultFileName{"OptimizationResults.json"};
struct CSVExportingSettings {
bool exportPE{false};
bool exportIterationOfMinima{false};
bool exportRawObjectiveValue{false};
CSVExportingSettings() {}
};
const CSVExportingSettings exportSettings;
struct JsonKeys {
inline static std::string optimizationVariables{"OptimizationVariables"};
inline static std::string baseTriangle{"BaseTriangle"};
inline static std::string Label{"Label"};
inline static std::string FullPatternLabel{"FullPatternLabel"};
inline static std::string Settings{"OptimizationSettings"};
inline static std::string FullPatternPotentialEnergies{"PE"};
inline static std::string fullPatternYoungsModulus{"youngsModulus"};
};
void saveObjectiveValuePlot(
const std::filesystem::path& outputImageDirPath) const {
std::vector<std::string> scenarioLabels(
objectiveValue.perSimulationScenario_total.size());
const double colorAxial = 1;
const double colorShear = 3;
const double colorBending = 5;
const double colorDome = 0.1;
const double colorSaddle = 0;
std::vector<double> colors(
objectiveValue.perSimulationScenario_total.size());
for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size();
scenarioIndex++) {
scenarioLabels[scenarioIndex] =
pSimulationJobs_reducedModel[scenarioIndex]->getLabel();
if (scenarioLabels[scenarioIndex].rfind("Axial", 0) == 0) {
colors[scenarioIndex] = colorAxial;
} else if (scenarioLabels[scenarioIndex].rfind("Shear", 0) == 0) {
colors[scenarioIndex] = colorShear;
} else if (scenarioLabels[scenarioIndex].rfind("Bending", 0) == 0) {
colors[scenarioIndex] = colorBending;
} else if (scenarioLabels[scenarioIndex].rfind("Dome", 0) == 0) {
colors[scenarioIndex] = colorDome;
} else if (scenarioLabels[scenarioIndex].rfind("Saddle", 0) == 0) {
colors[scenarioIndex] = colorSaddle;
} else {
std::cerr << "Label could not be identified" << std::endl;
}
}
std::vector<double> y(objectiveValue.perSimulationScenario_total.size());
for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size();
scenarioIndex++) {
y[scenarioIndex]
// =
// optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex]
// +
// optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex];
= objectiveValue
.perSimulationScenario_total_unweighted[scenarioIndex];
}
std::vector<double> x = matplot::linspace(0, y.size() - 1, y.size());
std::vector<double> markerSizes(y.size(), 5);
Utilities::createPlot("scenario index", "error", x, y, markerSizes, colors,
std::filesystem::path(outputImageDirPath)
.append("perScenarioObjectiveValues.svg"));
}
void save(const std::string& saveToPath,
const bool shouldExportDebugFiles = false) {
// clear directory
if (std::filesystem::exists(saveToPath)) {
for (const auto& entry :
std::filesystem::directory_iterator(saveToPath)) {
std::error_code ec;
std::filesystem::remove_all(entry.path(), ec);
}
}
std::filesystem::create_directories(saveToPath);
// Save optimal X
nlohmann::json json_optimizationResults;
json_optimizationResults[JsonKeys::Label] = label;
if (wasSuccessful) {
std::string jsonValue_optimizationVariables;
for (const auto& optimalXNameValuePair : optimalXNameValuePairs) {
jsonValue_optimizationVariables.append(optimalXNameValuePair.first +
",");
}
jsonValue_optimizationVariables
.pop_back(); // for deleting the last comma
json_optimizationResults[JsonKeys::optimizationVariables] =
jsonValue_optimizationVariables;
for (const auto& optimalXNameValuePair : optimalXNameValuePairs) {
json_optimizationResults[optimalXNameValuePair.first] =
optimalXNameValuePair.second;
}
}
// base triangle
json_optimizationResults[JsonKeys::baseTriangle] = std::vector<double>{
baseTriangle.cP0(0)[0], baseTriangle.cP0(0)[1], baseTriangle.cP0(0)[2],
baseTriangle.cP1(0)[0], baseTriangle.cP1(0)[1], baseTriangle.cP1(0)[2],
baseTriangle.cP2(0)[0], baseTriangle.cP2(0)[1], baseTriangle.cP2(0)[2]};
baseTrianglePattern.save(std::filesystem::path(saveToPath).string());
json_optimizationResults[JsonKeys::FullPatternLabel] =
baseTrianglePattern.getLabel();
// potential energies
// const int numberOfSimulationJobs =
// fullPatternSimulationJobs.size(); std::vector<double>
// fullPatternPE(numberOfSimulationJobs); for (int
// simulationScenarioIndex = 0; simulationScenarioIndex <
// numberOfSimulationJobs;
// simulationScenarioIndex++) {
// fullPatternPE[simulationScenarioIndex]
// =
// perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
// }
// json_optimizationResults[JsonKeys::FullPatternPotentialEnergies]
// = fullPatternPE;
////Save to json file
std::filesystem::path jsonFilePath(
std::filesystem::path(saveToPath).append(DefaultFileName));
std::ofstream jsonFile_optimizationResults(jsonFilePath.string());
jsonFile_optimizationResults << json_optimizationResults;
/*TODO: Refactor since the meshes are saved for each simulation scenario
* although they do not change*/
// Save jobs and meshes for each simulation scenario
if (shouldExportDebugFiles) {
// Save the reduced and full patterns
const std::filesystem::path simulationJobsPath(
std::filesystem::path(saveToPath).append("SimulationJobs"));
const int numberOfSimulationJobs = pSimulationJobs_pattern.size();
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < numberOfSimulationJobs;
simulationScenarioIndex++) {
const std::shared_ptr<SimulationJob>& pFullPatternSimulationJob =
pSimulationJobs_pattern[simulationScenarioIndex];
std::filesystem::path simulationJobFolderPath(
std::filesystem::path(simulationJobsPath)
.append(std::to_string(simulationScenarioIndex) + "_" +
pFullPatternSimulationJob->label));
std::filesystem::create_directories(simulationJobFolderPath);
const auto fullPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Full");
std::filesystem::create_directory(fullPatternDirectoryPath);
pFullPatternSimulationJob->save(fullPatternDirectoryPath.string());
const std::shared_ptr<SimulationJob>& pReducedPatternSimulationJob =
pSimulationJobs_reducedModel[simulationScenarioIndex];
const auto reducedPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Reduced");
if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
std::filesystem::create_directory(reducedPatternDirectoryPath);
}
pReducedPatternSimulationJob->save(
reducedPatternDirectoryPath.string());
}
// constexpr bool shouldSaveObjectiveValuePlot =
// shouldExportDebugFiles; if
// (shouldSaveObjectiveValuePlot)
// {
saveObjectiveValuePlot(saveToPath);
// }
}
CSVFile csv_resultsLocalFile(
std::filesystem::path(saveToPath).append("results.csv"), true);
csv_resultsLocalFile << "Name";
writeHeaderTo(csv_resultsLocalFile);
settings.writeHeaderTo(csv_resultsLocalFile);
csv_resultsLocalFile << endrow;
csv_resultsLocalFile << baseTrianglePattern.getLabel();
writeResultsTo(csv_resultsLocalFile);
settings.writeSettingsTo(csv_resultsLocalFile);
csv_resultsLocalFile << endrow;
// save minima info
// std::filesystem::path csvFilepathMinimaInfo =
// std::filesystem::path(saveToPath)
// .append("minimaInfo.csv");
// csvFile csv_minimaInfo(csvFilepathMinimaInfo, false);
// writeMinimaInfoTo(csv_minimaInfo);
settings.save(saveToPath);
#ifdef POLYSCOPE_DEFINED
std::cout << "Saved optimization results to:" << saveToPath << std::endl;
#endif
}
bool load(const std::filesystem::path& loadFromPath,
const bool& shouldLoadDebugFiles = false) {
assert(std::filesystem::is_directory(loadFromPath));
std::filesystem::path jsonFilepath(
std::filesystem::path(loadFromPath).append(DefaultFileName));
if (!std::filesystem::exists(jsonFilepath)) {
std::cerr << "Input path does not exist:" << loadFromPath << std::endl;
return false;
}
// Load optimal X
nlohmann::json json_optimizationResults;
std::ifstream ifs(
std::filesystem::path(loadFromPath).append(DefaultFileName));
ifs >> json_optimizationResults;
// std::cout << json_optimizationResults.dump() << std::endl;
label = json_optimizationResults.at(JsonKeys::Label);
std::string optimizationVariablesString =
*json_optimizationResults.find(JsonKeys::optimizationVariables);
std::string optimizationVariablesDelimeter = ",";
size_t pos = 0;
std::vector<std::string> optimizationVariablesNames;
while ((pos = optimizationVariablesString.find(
optimizationVariablesDelimeter)) != std::string::npos) {
const std::string variableName =
optimizationVariablesString.substr(0, pos);
optimizationVariablesNames.push_back(variableName);
optimizationVariablesString.erase(
0, pos + optimizationVariablesDelimeter.length());
}
optimizationVariablesNames.push_back(
optimizationVariablesString); // add last variable name
optimalXNameValuePairs.resize(optimizationVariablesNames.size());
for (int xVariable_index = 0;
xVariable_index < optimizationVariablesNames.size();
xVariable_index++) {
const std::string xVariable_name =
optimizationVariablesNames[xVariable_index];
const double xVariable_value =
*json_optimizationResults.find(xVariable_name);
optimalXNameValuePairs[xVariable_index] =
std::make_pair(xVariable_name, xVariable_value);
}
const std::string fullPatternLabel =
json_optimizationResults.at(JsonKeys::FullPatternLabel);
if (!baseTrianglePattern.load(std::filesystem::path(loadFromPath)
.append(fullPatternLabel + ".ply")
.string())) {
if (!baseTrianglePattern.load(
std::filesystem::path(loadFromPath)
.append(loadFromPath.stem().string() + ".ply")
.string())) {
if (!baseTrianglePattern.load(std::filesystem::path(loadFromPath)
.append(fullPatternLabel + ".obj")
.string())) {
baseTrianglePattern.load(
std::filesystem::path(loadFromPath)
.append(loadFromPath.stem().string() + ".obj")
.string());
}
}
}
std::vector<double> baseTriangleVertices =
json_optimizationResults.at(JsonKeys::baseTriangle);
baseTriangle.P0(0) =
CoordType(baseTriangleVertices[0], baseTriangleVertices[1],
baseTriangleVertices[2]);
baseTriangle.P1(0) =
CoordType(baseTriangleVertices[3], baseTriangleVertices[4],
baseTriangleVertices[5]);
baseTriangle.P2(0) =
CoordType(baseTriangleVertices[6], baseTriangleVertices[7],
baseTriangleVertices[8]);
if (json_optimizationResults.contains(JsonKeys::fullPatternYoungsModulus)) {
settings.youngsModulus_pattern =
json_optimizationResults.at(JsonKeys::fullPatternYoungsModulus);
} else {
settings.youngsModulus_pattern = 1 * 1e9;
}
const std::filesystem::path folderPath_simulationJobs(
std::filesystem::path(loadFromPath).append("SimulationJobs"));
if (shouldLoadDebugFiles &&
std::filesystem::exists(folderPath_simulationJobs)) {
const std::vector<std::filesystem::path> scenariosSortedByName = [&]() {
std::vector<std::filesystem::path> sortedByName;
for (auto& entry :
std::filesystem::directory_iterator(folderPath_simulationJobs))
sortedByName.push_back(entry.path());
std::sort(sortedByName.begin(), sortedByName.end(),
&Utilities::compareNat);
return sortedByName;
}();
for (const auto& simulationScenarioPath : scenariosSortedByName) {
if (!std::filesystem::is_directory(simulationScenarioPath)) {
continue;
}
// Load full job
const auto fullJobFilepath = Utilities::getFilepathWithExtension(
std::filesystem::path(simulationScenarioPath).append("Full"),
".json");
SimulationJob fullJob;
fullJob.load(fullJobFilepath.string());
fullJob.pMesh->setBeamMaterial(0.3, settings.youngsModulus_pattern);
pSimulationJobs_pattern.push_back(
std::make_shared<SimulationJob>(fullJob));
// Load reduced job
const auto reducedJobFilepath = Utilities::getFilepathWithExtension(
std::filesystem::path(simulationScenarioPath).append("Reduced"),
".json");
SimulationJob reducedJob;
reducedJob.load(reducedJobFilepath.string());
applyOptimizationResults_elements(*this, reducedJob.pMesh);
pSimulationJobs_reducedModel.push_back(
std::make_shared<SimulationJob>(reducedJob));
}
}
settings.load(
std::filesystem::path(loadFromPath).append(Settings::defaultFilename));
return true;
}
template <typename MeshType>
static void applyOptimizationResults_reducedModel_nonFanned(
const ReducedModelOptimization::Results&
reducedPattern_optimizationResults,
const vcg::Triangle3<double>& patternBaseTriangle,
MeshType& reducedModel) {
std::unordered_map<std::string, double> optimalXVariables(
reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
reducedPattern_optimizationResults.optimalXNameValuePairs.end());
assert((optimalXVariables.contains("R") &&
optimalXVariables.contains("Theta")) ||
(optimalXVariables.contains("HexSize") &&
optimalXVariables.contains("HexAngle")));
if (optimalXVariables.contains("HexSize")) {
applyOptimizationResults_reducedModel_nonFanned(
optimalXVariables.at("HexSize"), optimalXVariables.at("HexAngle"),
patternBaseTriangle, reducedModel);
return;
}
applyOptimizationResults_reducedModel_nonFanned(
optimalXVariables.at("R"), optimalXVariables.at("Theta"),
patternBaseTriangle, reducedModel);
}
template <typename MeshType>
static void applyOptimizationResults_reducedModel_nonFanned(
const double& hexSize,
const double& hexAngle,
const vcg::Triangle3<double>& patternBaseTriangle,
MeshType& reducedModel) {
// Set optimal geometrical params of the reduced pattern
// CoordType center_barycentric(1, 0, 0);
// CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
// CoordType movableVertex_barycentric(
// (center_barycentric * (1 - hexSize) +
// interfaceEdgeMiddle_barycentric));
CoordType movableVertex_barycentric(1 - hexSize, hexSize / 2, hexSize / 2);
reducedModel.vert[0].P() =
patternBaseTriangle.cP(0) * movableVertex_barycentric[0] +
patternBaseTriangle.cP(1) * movableVertex_barycentric[1] +
patternBaseTriangle.cP(2) * movableVertex_barycentric[2];
if (hexAngle != 0) {
reducedModel.vert[0].P() =
vcg::RotationMatrix(CoordType{0, 0, 1}, vcg::math::ToRad(hexAngle)) *
reducedModel.vert[0].cP();
}
// for (int rotationCounter = 0;
// rotationCounter < ReducedModelOptimizer::fanSize;
// rotationCounter++) {
// pReducedPatternSimulationEdgeMesh->vert[2 * rotationCounter].P()
// =
// vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
// vcg::math::ToRad(60.0 *
// rotationCounter))
// * baseTriangleMovableVertexPosition;
// }
// reducedPattern.registerForDrawing();
// polyscope::show();
// CoordType p0 = reducedPattern.vert[0].P();
// CoordType p1 = reducedPattern.vert[1].P();
// int i = 0;
// i++;
}
static void applyOptimizationResults_elements(
const ReducedModelOptimization::Results&
reducedPattern_optimizationResults,
const std::shared_ptr<SimulationEdgeMesh>&
pReducedModel_SimulationEdgeMesh) {
assert(CrossSectionType::name == RectangularBeamDimensions::name);
// Set reduced parameters fron the optimization results
std::unordered_map<std::string, double> optimalXVariables(
reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
reducedPattern_optimizationResults.optimalXNameValuePairs.end());
const std::string ALabel = "A";
if (optimalXVariables.contains(ALabel)) {
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
for (int ei = 0; ei < pReducedModel_SimulationEdgeMesh->EN(); ei++) {
Element& e = pReducedModel_SimulationEdgeMesh->elements[ei];
e.setDimensions(elementDimensions);
}
}
const double poissonsRatio = 0.3;
const std::string ymLabel = "E";
if (optimalXVariables.contains(ymLabel)) {
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
for (int ei = 0; ei < pReducedModel_SimulationEdgeMesh->EN(); ei++) {
Element& e = pReducedModel_SimulationEdgeMesh->elements[ei];
e.setMaterial(elementMaterial);
}
}
const std::string JLabel = "J";
if (optimalXVariables.contains(JLabel)) {
const double J = optimalXVariables.at(JLabel);
for (int ei = 0; ei < pReducedModel_SimulationEdgeMesh->EN(); ei++) {
Element& e = pReducedModel_SimulationEdgeMesh->elements[ei];
e.dimensions.inertia.J = J;
}
}
const std::string I2Label = "I2";
if (optimalXVariables.contains(I2Label)) {
const double I2 = optimalXVariables.at(I2Label);
for (int ei = 0; ei < pReducedModel_SimulationEdgeMesh->EN(); ei++) {
Element& e = pReducedModel_SimulationEdgeMesh->elements[ei];
e.dimensions.inertia.I2 = I2;
}
}
const std::string I3Label = "I3";
if (optimalXVariables.contains(I3Label)) {
const double I3 = optimalXVariables.at(I3Label);
for (int ei = 0; ei < pReducedModel_SimulationEdgeMesh->EN(); ei++) {
Element& e = pReducedModel_SimulationEdgeMesh->elements[ei];
e.dimensions.inertia.I3 = I3;
}
}
pReducedModel_SimulationEdgeMesh->reset();
}
#if POLYSCOPE_DEFINED
void draw(const std::vector<int>& desiredSimulationScenariosIndices =
std::vector<int>()) const {
PolyscopeInterface::init();
assert(pSimulationJobs_pattern.size() ==
pSimulationJobs_reducedModel.size());
pSimulationJobs_pattern[0]->pMesh->registerForDrawing(
Colors::patternInitial);
pSimulationJobs_reducedModel[0]->pMesh->registerForDrawing(
Colors::reducedInitial, false);
const int numberOfSimulationJobs = pSimulationJobs_pattern.size();
const std::vector<int> scenariosToDraw = [&]() {
if (desiredSimulationScenariosIndices.empty()) {
std::vector<int> v(numberOfSimulationJobs);
std::iota(v.begin(), v.end(), 0); // draw all
return v;
} else {
return desiredSimulationScenariosIndices;
}
}();
for (const int& simulationJobIndex : scenariosToDraw) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob>& pSimulationJob_pattern =
pSimulationJobs_pattern[simulationJobIndex];
pSimulationJob_pattern->registerForDrawing(
pSimulationJobs_pattern[0]->pMesh->getLabel());
std::unique_ptr<SimulationModel> pPatternSimulator =
SimulationModelFactory::create(
settings.simulationModelLabel_groundTruth);
pPatternSimulator->setStructure(pSimulationJob_pattern->pMesh);
SimulationResults simulationResults_pattern =
pPatternSimulator->executeSimulation(pSimulationJob_pattern);
// ChronosEulerSimulationModel simulator_pattern;
// simulator_pattern.settings.analysisType =
// ChronosEulerSimulationModel::Settings::AnalysisType::NonLinear;
// SimulationResults simulationResults_pattern =
// simulator_pattern.executeSimulation(pSimulationJob_pattern);
simulationResults_pattern.registerForDrawing(Colors::patternDeformed,
true);
// SimulationResults fullModelLinearResults =
// linearSimulator.executeSimulation(pFullPatternSimulationJob);
// fullModelLinearResults.setLabelPrefix("linear");
// fullModelLinearResults.registerForDrawing(Colors::fullDeformed,false);
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob>& pSimulationJob_reducedModel =
pSimulationJobs_reducedModel[simulationJobIndex];
// pReducedPatternSimulationJob->pMesh->registerForDrawing();
// polyscope::show();
// SimulationResults reducedModelResults =
// drmSimulator.executeSimulation(
// pReducedPatternSimulationJob);
// reducedModelResults.registerForDrawing(Colors::reducedDeformed,
// false);
// SimulationResults reducedModelResults
// =
// drmSimulator.executeSimulation(pReducedPatternSimulationJob,
// DRMSimulationModel::Settings());
std::unique_ptr<SimulationModel> pSimulator_reducedModel =
SimulationModelFactory::create(
settings.simulationModelLabel_reducedModel);
pSimulator_reducedModel->setStructure(pSimulationJob_reducedModel->pMesh);
// ChronosEulerSimulationModel simulator_reducedModel;
// simulator_reducedModel.settings.analysisType =
// ChronosEulerSimulationModel::Settings::AnalysisType::Linear;
SimulationResults simulationResults_reducedModel =
pSimulator_reducedModel->executeSimulation(
pSimulationJob_reducedModel);
simulationResults_reducedModel.registerForDrawing(Colors::reducedDeformed,
true);
polyscope::options::programName =
pSimulationJobs_pattern[0]->pMesh->getLabel();
// polyscope::view::resetCameraToHomeView();
polyscope::show();
// Save a screensh
const std::string screenshotFilename =
"/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/RodModelOptimizationForPatterns/Results/Images/" +
pSimulationJobs_pattern[0]->pMesh->getLabel() + "_" +
pSimulationJob_pattern->getLabel();
polyscope::screenshot(screenshotFilename, false);
pSimulationJob_pattern->unregister(
pSimulationJobs_pattern[0]->pMesh->getLabel());
simulationResults_reducedModel.unregister();
simulationResults_pattern.unregister();
// double error =
// computeError(reducedModelLinearResults.displacements,
// fullModelResults.displacements);
// std::cout << "Error of simulation scenario "
// <<
// simulationScenarioStrings[simulationScenarioIndex]
// << " is " << error
// << std::endl;
}
}
#endif // POLYSCOPE_DEFINED
void saveMeshFiles() const {
const int numberOfSimulationJobs = pSimulationJobs_pattern.size();
assert(numberOfSimulationJobs != 0 &&
pSimulationJobs_pattern.size() ==
pSimulationJobs_reducedModel.size());
pSimulationJobs_pattern[0]->pMesh->save(
"undeformed " + pSimulationJobs_pattern[0]->pMesh->getLabel() + ".ply");
pSimulationJobs_reducedModel[0]->pMesh->save(
"undeformed " + pSimulationJobs_reducedModel[0]->pMesh->getLabel() +
".ply");
DRMSimulationModel simulator_drm;
LinearSimulationModel simulator_linear;
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob>& pFullPatternSimulationJob =
pSimulationJobs_pattern[simulationJobIndex];
DRMSimulationModel::Settings drmSettings;
SimulationResults fullModelResults = simulator_drm.executeSimulation(
pFullPatternSimulationJob, drmSettings);
fullModelResults.saveDeformedModel();
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob>& pReducedPatternSimulationJob =
pSimulationJobs_reducedModel[simulationJobIndex];
SimulationResults reducedModelResults =
simulator_linear.executeSimulation(pReducedPatternSimulationJob);
reducedModelResults.saveDeformedModel();
}
}
void writeHeaderTo(CSVFile& os) const {
if (exportSettings.exportRawObjectiveValue) {
os << "Total raw Obj value";
}
os << "Total Obj value";
if (exportSettings.exportIterationOfMinima) {
os << "Iteration of minima";
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < pSimulationJobs_pattern.size();
simulationScenarioIndex++) {
const std::string simulationScenarioName =
pSimulationJobs_pattern[simulationScenarioIndex]->getLabel();
os << "Obj value Trans " + simulationScenarioName;
os << "Obj value Rot " + simulationScenarioName;
os << "Obj value Total " + simulationScenarioName;
}
if (exportSettings.exportPE) {
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < pSimulationJobs_pattern.size();
simulationScenarioIndex++) {
const std::string simulationScenarioName =
pSimulationJobs_pattern[simulationScenarioIndex]->getLabel();
os << "PE " + simulationScenarioName;
}
}
for (const auto& nameValuePair : optimalXNameValuePairs) {
os << nameValuePair.first;
}
os << "Time(s)";
// os << "#Crashes";
// os << "notes";
}
void writeHeaderTo(
std::vector<CSVFile*>& vectorOfPointersToOutputStreams) const {
for (CSVFile* outputStream : vectorOfPointersToOutputStreams) {
writeHeaderTo(*outputStream);
}
}
void writeResultsTo(CSVFile& os) const {
if (exportSettings.exportRawObjectiveValue) {
os << objectiveValue.totalRaw;
}
os << objectiveValue.total;
if (exportSettings.exportIterationOfMinima &&
!objectiveValueHistory_iteration.empty()) {
os << objectiveValueHistory_iteration.back();
}
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < pSimulationJobs_pattern.size();
simulationScenarioIndex++) {
os << objectiveValue
.perSimulationScenario_translational[simulationScenarioIndex];
os << objectiveValue
.perSimulationScenario_rotational[simulationScenarioIndex];
os << objectiveValue.perSimulationScenario_total[simulationScenarioIndex];
}
if (exportSettings.exportPE) {
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < pSimulationJobs_pattern.size();
simulationScenarioIndex++) {
os << perScenario_fullPatternPotentialEnergy[simulationScenarioIndex];
}
}
for (const auto& optimalXNameValuePair : optimalXNameValuePairs) {
os << optimalXNameValuePair.second;
}
os << time;
// if (numberOfSimulationCrashes == 0) {
// os << "-";
// } else {
// os << numberOfSimulationCrashes;
// }
// os << notes;
}
void writeResultsTo(
std::vector<CSVFile*>& vectorOfPointersToOutputStreams) const {
for (CSVFile*& outputStream : vectorOfPointersToOutputStreams) {
writeResultsTo(*outputStream);
}
}
void writeMinimaInfoTo(CSVFile& outputCsv) {
outputCsv << "Iteration";
outputCsv << "Objective value";
for (int objectiveValueIndex = 0;
objectiveValueIndex < objectiveValueHistory.size();
objectiveValueIndex++) {
outputCsv.endrow();
outputCsv << objectiveValueHistory_iteration[objectiveValueIndex];
outputCsv << objectiveValueHistory[objectiveValueIndex];
}
}
};
enum SimulationModelTag { DRM, Linear };
inline bool Settings::ObjectiveWeights::operator==(
const ObjectiveWeights& other) const {
return this->translational == other.translational &&
this->rotational == other.rotational;
}
} // namespace ReducedModelOptimization
#endif // REDUCEDMODELOPTIMIZER_STRUCTS_HPP