MySources/reducedmodeloptimizer_struc...

1062 lines
50 KiB
C++

#ifndef REDUCEDMODELOPTIMIZER_STRUCTS_HPP
#define REDUCEDMODELOPTIMIZER_STRUCTS_HPP
#include "csvfile.hpp"
#include "drmsimulationmodel.hpp"
#include "linearsimulationmodel.hpp"
#include "simulation_structs.hpp"
#include "unordered_map"
#include <string>
namespace ReducedModelOptimization {
enum BaseSimulationScenario {
Axial,
Shear,
Bending,
Dome,
Saddle,
NumberOfBaseSimulationScenarios
};
inline static std::vector<std::string> baseSimulationScenarioNames{"Axial",
"Shear",
"Bending",
"Dome",
"Saddle"};
struct Colors
{
inline static std::array<double, 3> fullInitial{0.094, 0.094, 0.094};
// inline static std::array<double, 3> fullDeformed{0.583333, 0.890196, 0.109804};
inline static std::array<double, 3> fullDeformed{0.094, 0.094, 0.094};
inline static std::array<double, 3> reducedInitial{0.890196, 0.109804, 0.193138};
inline static std::array<double, 3> reducedDeformed{0.262, 0.627, 0.910};
};
struct xRange
{
std::string label{};
double min{0};
double max{0};
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::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
{E, A, I2, I3, J, R, Theta}};
// {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", 0.001, 1000},
xRange{"A", 0.001, 1000},
xRange{"I2", 0.001, 1000},
xRange{"I3", 0.001, 1000},
xRange{"J", 0.001, 1000},
xRange{"R", 0.05, 0.95},
xRange{"Theta", -30, 30}};
int numberOfFunctionCalls{100000};
double solverAccuracy{1e-3};
double translationEpsilon{3e-3};
double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
struct ObjectiveWeights
{
double translational{1.2};
double rotational{0.8};
bool operator==(const ObjectiveWeights &other) const
{
return this->translational == other.translational
&& this->rotational == other.rotational;
}
};
std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios> perBaseScenarioObjectiveWeights;
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;
}
struct JsonKeys
{
inline static std::string OptimizationStrategy{"OptimizationStrategy"};
inline static std::string OptimizationStrategyGroupWeights{
"OptimizationStrategyGroupWeights"};
inline static std::string OptimizationVariables{"OptimizationVariables"};
inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"};
inline static std::string SolverAccuracy{"SolverAccuracy"};
inline static std::string ObjectiveWeights{"ObjectiveWeight"};
};
void setDefault()
{
// ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
// ReducedPatternOptimization::xRange beamA{"A", 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{"R", 0.05, 0.95};
// ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
// variablesRanges = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
// numberOfFunctionCalls = 100000;
enum OptimizationParameterComparisonScenarioIndex {
AllVar,
GeoYM,
MatGeo,
YMMat_Geo,
YM_MatGeo,
MatGeo_YM,
Geo_YM_Mat,
YM_Geo_Mat,
Geo_Mat,
YMGeo_Mat,
NumberOfScenarios
};
const std::vector<
std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>>
optimizationParameters = [&]() {
std::vector<
std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>>
optimizationParameters(NumberOfScenarios);
using namespace ReducedModelOptimization;
optimizationParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}};
optimizationParameters[GeoYM] = {{R, Theta, E}};
optimizationParameters[MatGeo] = {{A, I2, I3, J, R, Theta}};
optimizationParameters[YMMat_Geo] = {{E, A, I2, I3, J}, {R, Theta}};
optimizationParameters[YM_MatGeo] = {{E}, {A, I2, I3, J, R, Theta}};
optimizationParameters[MatGeo_YM] = {{A, I2, I3, J, R, Theta}, {E}};
optimizationParameters[Geo_YM_Mat] = {{R, Theta}, {E}, {A, I2, I3, J}};
optimizationParameters[YM_Geo_Mat] = {{E}, {R, Theta}, {A, I2, I3, J}};
optimizationParameters[Geo_Mat] = {{R, Theta}, {A, I2, I3, J}};
optimizationParameters[YMGeo_Mat] = {{E, R, Theta}, {A, I2, I3, J}};
return optimizationParameters;
}();
constexpr OptimizationParameterComparisonScenarioIndex scenario = AllVar;
optimizationStrategy = optimizationParameters[scenario];
if (scenario == YMGeo_Mat) {
optimizationVariablesGroupsWeights = {0.15, 0.85};
}
// normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
// translationNormalizationEpsilon = 0;
// rotationNormalizationEpsilon = vcg::math::ToRad(0.0);
// solverAccuracy = 1e-2;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational = 1.174011;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational = 1.690460;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational = 1.812089;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational = 0.127392;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational = 0.503043;
// for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
// baseScenario++) {
// perBaseScenarioObjectiveWeights[baseScenario].rotational
// = 2 - perBaseScenarioObjectiveWeights[baseScenario].translational;
// }
}
void save(const std::filesystem::path &saveToPath)
{
assert(std::filesystem::is_directory(saveToPath));
nlohmann::json json;
json[GET_VARIABLE_NAME(optimizationStrategy)] = optimizationStrategy;
if (optimizationVariablesGroupsWeights.has_value()) {
json[GET_VARIABLE_NAME(ptimizationStrategyGroupWeights)]
= optimizationVariablesGroupsWeights.value();
}
//write x ranges
const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
xRangesAsTuples = [=]() {
std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
xRangesAsTuples;
for (int optimizationParameterIndex = E;
optimizationParameterIndex != NumberOfOptimizationVariables;
optimizationParameterIndex++) {
xRangesAsTuples[optimizationParameterIndex]
= variablesRanges[optimizationParameterIndex].toTuple();
}
return xRangesAsTuples;
}();
json[JsonKeys::OptimizationVariables] = xRangesAsTuples;
// for (size_t xRangeIndex = 0; xRangeIndex < variablesRanges.size(); xRangeIndex++) {
// const auto &xRange = variablesRanges[xRangeIndex];
// json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)]
// = xRange.toString();
// }
json[GET_VARIABLE_NAME(numberOfFunctionCalls)] = numberOfFunctionCalls;
json[GET_VARIABLE_NAME(solverAccuracy)] = solverAccuracy;
//Objective weights
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
std::transform(perBaseScenarioObjectiveWeights.begin(),
perBaseScenarioObjectiveWeights.end(),
objectiveWeightsPairs.begin(),
[](const ObjectiveWeights &objectiveWeights) {
return std::make_pair(objectiveWeights.translational,
objectiveWeights.rotational);
});
json[JsonKeys::ObjectiveWeights] = objectiveWeightsPairs;
json[GET_VARIABLE_NAME(translationEpsilon)] = translationEpsilon;
json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(
angularDistanceEpsilon);
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)
{
nlohmann::json json;
if (!std::filesystem::exists(jsonFilepath)) {
std::cerr << "Optimization settings could not be loaded because input filepath does "
"not exist:"
<< jsonFilepath << std::endl;
return false;
}
std::ifstream ifs(jsonFilepath);
ifs >> 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(numberOfFunctionCalls))) {
numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(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)]));
}
// perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.rotational = 2 - objectiveWeights.translational;
return true;
}
std::string toString() const
{
std::string settingsString;
if (!variablesRanges.empty()) {
std::string xRangesString;
for (const xRange &range : variablesRanges) {
xRangesString += range.toString() + " ";
}
settingsString += xRangesString;
}
settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls)
+ " Accuracy=" + std::to_string(solverAccuracy);
for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
baseScenario++) {
+" W_t=" + std::to_string(perBaseScenarioObjectiveWeights[baseScenario].translational)
+ " W_r="
+ std::to_string(perBaseScenarioObjectiveWeights[baseScenario].rotational);
}
return settingsString;
}
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 << JsonKeys::ObjectiveWeights;
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 << numberOfFunctionCalls;
os << solverAccuracy;
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);
}
};
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.numberOfFunctionCalls == settings2.numberOfFunctionCalls
&& settings1.variablesRanges == settings2.variablesRanges
&& settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
&& settings1.translationEpsilon
== settings2.translationEpsilon;
}
inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &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>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
double fullPatternYoungsModulus{0};
PatternGeometry baseTriangleFullPattern; //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;
struct CSVExportingSettings
{
bool exportPE{false};
bool exportIterationOfMinima{false};
bool exportRawObjectiveValue{true};
CSVExportingSettings() {}
};
const CSVExportingSettings exportSettings;
struct JsonKeys
{
inline static std::string filename{"OptimizationResults.json"};
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 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]};
baseTriangleFullPattern.save(std::filesystem::path(saveToPath).string());
json_optimizationResults[JsonKeys::FullPatternLabel] = baseTriangleFullPattern.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;
json_optimizationResults[JsonKeys::fullPatternYoungsModulus] = fullPatternYoungsModulus;
////Save to json file
std::filesystem::path jsonFilePath(
std::filesystem::path(saveToPath).append(JsonKeys::filename));
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 = fullPatternSimulationJobs.size();
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < numberOfSimulationJobs;
simulationScenarioIndex++) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= fullPatternSimulationJobs[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
= reducedPatternSimulationJobs[simulationScenarioIndex];
const auto reducedPatternDirectoryPath
= std::filesystem::path(simulationJobFolderPath).append("Reduced");
if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
std::filesystem::create_directory(reducedPatternDirectoryPath);
}
pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string());
}
}
//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::string &loadFromPath, const bool &shouldLoadDebugFiles = false)
{
assert(std::filesystem::is_directory(loadFromPath));
std::filesystem::path jsonFilepath(
std::filesystem::path(loadFromPath).append(JsonKeys::filename));
if (!std::filesystem::exists(jsonFilepath)) {
return false;
}
//Load optimal X
nlohmann::json json_optimizationResults;
std::ifstream ifs(std::filesystem::path(loadFromPath).append(JsonKeys::filename));
ifs >> json_optimizationResults;
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);
baseTriangleFullPattern.load(
std::filesystem::path(loadFromPath).append(fullPatternLabel + ".ply").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)) {
fullPatternYoungsModulus = json_optimizationResults.at(
JsonKeys::fullPatternYoungsModulus);
} else {
fullPatternYoungsModulus = 1 * 1e9;
}
if (shouldLoadDebugFiles) {
const std::filesystem::path simulationJobsFolderPath(
std::filesystem::path(loadFromPath).append("SimulationJobs"));
std::vector<std::filesystem::path> sortedByName;
for (auto &entry : std::filesystem::directory_iterator(simulationJobsFolderPath))
sortedByName.push_back(entry.path());
std::sort(sortedByName.begin(), sortedByName.end(), &Utilities::compareNat);
for (const auto &simulationScenarioPath : sortedByName) {
if (!std::filesystem::is_directory(simulationScenarioPath)) {
continue;
}
const auto fullJobFilepath = Utilities::getFilepathWithExtension(
std::filesystem::path(simulationScenarioPath).append("Full"), ".json");
SimulationJob fullJob;
fullJob.load(fullJobFilepath.string());
fullJob.pMesh->setBeamMaterial(0.3, fullPatternYoungsModulus);
fullPatternSimulationJobs.push_back(std::make_shared<SimulationJob>(fullJob));
const auto reducedJobFilepath = Utilities::getFilepathWithExtension(
std::filesystem::path(simulationScenarioPath).append("Reduced"), ".json");
SimulationJob reducedJob;
reducedJob.load(reducedJobFilepath.string());
// applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh);
applyOptimizationResults_elements(*this, reducedJob.pMesh);
reducedPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(reducedJob));
}
}
settings.load(std::filesystem::path(loadFromPath).append(Settings::defaultFilename));
return true;
}
template<typename MeshType>
static void applyOptimizationResults_innerHexagon(
const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
const vcg::Triangle3<double> &patternBaseTriangle,
MeshType &reducedPattern)
{
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_innerHexagon(optimalXVariables.at("HexSize"),
optimalXVariables.at("HexAngle"),
patternBaseTriangle,
reducedPattern);
return;
}
applyOptimizationResults_innerHexagon(optimalXVariables.at("R"),
optimalXVariables.at("Theta"),
patternBaseTriangle,
reducedPattern);
}
template<typename MeshType>
static void applyOptimizationResults_innerHexagon(
const double &hexSize,
const double &hexAngle,
const vcg::Triangle3<double> &patternBaseTriangle,
MeshType &reducedPattern)
{
//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);
reducedPattern.vert[0].P() = patternBaseTriangle.cP(0) * movableVertex_barycentric[0]
+ patternBaseTriangle.cP(1) * movableVertex_barycentric[1]
+ patternBaseTriangle.cP(2) * movableVertex_barycentric[2];
if (hexAngle != 0) {
reducedPattern.vert[0].P() = vcg::RotationMatrix(CoordType{0, 0, 1},
vcg::math::ToRad(hexAngle))
* reducedPattern.vert[0].cP();
}
// for (int rotationCounter = 0;
// rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
// pReducedPatternSimulationMesh->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<SimulationMesh> &pReducedPattern_simulationMesh)
{
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";
assert(optimalXVariables.contains(ALabel));
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
const double poissonsRatio = 0.3;
const std::string ymLabel = "E";
assert(optimalXVariables.contains(ymLabel));
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
const std::string JLabel = "J";
assert(optimalXVariables.contains(JLabel));
const double J = optimalXVariables.at(JLabel);
const std::string I2Label = "I2";
assert(optimalXVariables.contains(I2Label));
const double I2 = optimalXVariables.at(I2Label);
const std::string I3Label = "I3";
assert(optimalXVariables.contains(I3Label));
const double I3 = optimalXVariables.at(I3Label);
for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pReducedPattern_simulationMesh->elements[ei];
e.setDimensions(elementDimensions);
e.setMaterial(elementMaterial);
e.dimensions.inertia.J = J;
e.dimensions.inertia.I2 = I2;
e.dimensions.inertia.I3 = I3;
}
}
#if POLYSCOPE_DEFINED
void draw() const {
PolyscopeInterface::init();
DRMSimulationModel drmSimulator;
LinearSimulationModel linearSimulator;
assert(fullPatternSimulationJobs.size() == reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::fullInitial);
reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial, false);
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
for (int simulationJobIndex = 0; simulationJobIndex < numberOfSimulationJobs;
simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= fullPatternSimulationJobs[simulationJobIndex];
pFullPatternSimulationJob->registerForDrawing(
fullPatternSimulationJobs[0]->pMesh->getLabel());
DRMSimulationModel::Settings drmSettings;
SimulationResults fullModelResults
= drmSimulator.executeSimulation(pFullPatternSimulationJob, drmSettings);
fullModelResults.registerForDrawing(Colors::fullDeformed, true);
// SimulationResults fullModelLinearResults =
// linearSimulator.executeSimulation(pFullPatternSimulationJob);
// fullModelLinearResults.setLabelPrefix("linear");
// fullModelLinearResults.registerForDrawing(Colors::fullDeformed,false);
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob
= reducedPatternSimulationJobs[simulationJobIndex];
// pReducedPatternSimulationJob->pMesh->registerForDrawing();
// polyscope::show();
// SimulationResults reducedModelResults = drmSimulator.executeSimulation(
// pReducedPatternSimulationJob);
// reducedModelResults.registerForDrawing(Colors::reducedDeformed, false);
// SimulationResults reducedModelResults
// = drmSimulator.executeSimulation(pReducedPatternSimulationJob,
// DRMSimulationModel::Settings());
SimulationResults reducedModelResults = linearSimulator.executeSimulation(
pReducedPatternSimulationJob);
reducedModelResults.setLabelPrefix("linear");
reducedModelResults.registerForDrawing(Colors::reducedDeformed, true);
polyscope::options::programName = fullPatternSimulationJobs[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/"
+ fullPatternSimulationJobs[0]->pMesh->getLabel() + "_"
+ pFullPatternSimulationJob->getLabel();
polyscope::screenshot(screenshotFilename, false);
pFullPatternSimulationJob->unregister(fullPatternSimulationJobs[0]->pMesh->getLabel());
fullModelResults.unregister();
// reducedModelResults.unregister();
reducedModelResults.unregister();
// fullModelLinearResults.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 = fullPatternSimulationJobs.size();
assert(numberOfSimulationJobs != 0 &&
fullPatternSimulationJobs.size() ==
reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->save(
"undeformed " + fullPatternSimulationJobs[0]->pMesh->getLabel() + ".ply");
reducedPatternSimulationJobs[0]->pMesh->save(
"undeformed " + reducedPatternSimulationJobs[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
= fullPatternSimulationJobs[simulationJobIndex];
DRMSimulationModel::Settings drmSettings;
SimulationResults fullModelResults
= simulator_drm.executeSimulation(pFullPatternSimulationJob, drmSettings);
fullModelResults.saveDeformedModel();
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob
= reducedPatternSimulationJobs[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 < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
const std::string simulationScenarioName
= fullPatternSimulationJobs[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 < fullPatternSimulationJobs.size();
simulationScenarioIndex++) {
const std::string simulationScenarioName
= fullPatternSimulationJobs[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 < fullPatternSimulationJobs.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 < fullPatternSimulationJobs.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 };
} // namespace ReducedPatternOptimization
#endif // REDUCEDMODELOPTIMIZER_STRUCTS_HPP