#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,
    S,
    NumberOfBaseSimulationScenarios
};
inline static std::vector<std::string> baseSimulationScenarioNames{"Axial",
                                                                   "Shear",
                                                                   "Bending",
                                                                   "Dome",
                                                                   "Saddle",
                                                                   "S"};

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::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0.590241,
    //                                                                                  0.888372,
    //                                                                                  0.368304,
    //                                                                                  0.0127508,
    //                                                                                  1.18079,
    //                                                                                  0}; //final
    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0.590241 / 3,
                                                                                  0.588372 / 3,
                                                                                  0.368304,
                                                                                  0.1,
                                                                                  1.18 / 2,
                                                                                  0}; //final
    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0, 0, 0, 0.1, 0};

    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMaxMagnitudes{20.85302947095844,
    //                                  1.8073431893126763,
    //                                  0.2864731720436702,
    //                                  0.14982243562639147,
    //                                  0.18514829631059054};//median
    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMaxMagnitudes{1.1725844893199244,
    //                                  0.3464275389927846,
    //                                  0.09527915004635197,
    //                                  0.06100757786262501,
    //                                  0.10631914784812076}; //5_1565 0.03axial
    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMaxMagnitudes{/*15*/ 0 /*1.711973658196369*/,
    //                                  1.878077115238504,
    //                                  0.8,
    //                                  0.15851675178327318,
    //                                  0.8,
    //                                  /*1.711973658196369*/ 0}; //custom
    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMaxMagnitudes{0,
    //                                  14.531854387244818,
    //                                  0.38321932238436796,
    //                                  0.21381267870193282,
    //                                  0.28901381608791094,
    //                                  1.711973658196369}; //9_14423
    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMaxMagnitudes{1.1725844893199244,
    //                                  0.3464275389927846,
    //                                  0.09527915004635197,
    //                                  0.06100757786262501,
    //                                  0.10631914784812076}; //5_1565 0.03axial
    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMagnitudes{9.82679, 0.138652, 0.247242, 0.739443, 0.00675865}; //Hyperparam opt
    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{
    //        30, 8, 0.4421382884449713, 0.22758433903942452, 0.3247935583883217};

    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMagnitudes{10 * 6.310485381644259,
    //                               10 * 1.7100142258819078,
    //                               10 * 0.18857048204421728,
    //                               10 * 0.10813697502645818,
    //                               10 * 0.11982893539207524}; //6_338
    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMagnitudes{7.72224,
    //                                                                               7.72224,
    //                                                                               0.89468,
    //                                                                               0.445912,
    //                                                                               0.625905};
    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMagnitudes{0.407714,
    //                                                                               22.3524,
    //                                                                               0.703164,
    //                                                                               0.0226138,
    //                                                                               0.161316};
    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMagnitudes{2, 1, 0.4, 0.2, 0.2}; //8_15444 magnitudes from randomBending0

    //        std::array<double, NumberOfBaseSimulationScenarios>
    //            baseScenarioMagnitudes{1.0600375226325425,
    //                                   0.6381040280710403,
    //                                   0.17201755995098306,
    //                                   0.0706601822149856,
    //                                   0.13578373479448072}; //8_15444 magnitudes from displacements
    //    std::array<double, NumberOfBaseSimulationScenarios>
    //        baseScenarioMagnitudes{10 * 1.0600375226325425,
    //                               10 * 0.6381040280710403,
    //                               10 * 0.17201755995098306,
    //                               10 * 0.0706601822149856,
    //                               10 * 0.13578373479448072}; //8_15444 magnitudes from displacements

    //    std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes;
    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", 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}};
    struct SettingsPSO
    {
        int numberOfParticles{200};
#ifdef USE_PSO
        inline static std::string optimizerName{"pso"};
#else
        inline static std::string optimizerName{"sa"};
#endif
    } pso;

    struct SettingsDlibGlobal
    {
        int numberOfFunctionCalls{100000};
    } dlib;

    double solverAccuracy{1e-2};
    double translationEpsilon{4e-3};
    //    double translationEpsilon{0};
    //    double angularDistanceEpsilon{vcg::math::ToRad(2.0)};
    double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
    double targetBaseTriangleSize{0.03};
    std::filesystem::path intermediateResultsDirectoryPath;
    struct ObjectiveWeights
    {
        double translational{1.2};
        double rotational{0.8};
        bool operator==(const ObjectiveWeights &other) const;
    };
    std::array<ObjectiveWeights, NumberOfBaseSimulationScenarios> 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;
    }

    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"};
    };

    nlohmann::json toJson() const
    {
        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(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);
        json[GET_VARIABLE_NAME(targetBaseTriangleSize)] = targetBaseTriangleSize;
        json[GET_VARIABLE_NAME(baseScenarioMaxMagnitudes)] = baseScenarioMaxMagnitudes;

#ifdef USE_ENSMALLEN
#ifdef USE_PSO
        json[GET_VARIABLE_NAME(pso.numberOfParticles)] = pso.numberOfParticles;
#endif
        json[GET_VARIABLE_NAME(pso.optimizerName)] = pso.optimizerName;
#else
        json[GET_VARIABLE_NAME(dlib.numberOfFunctionCalls)] = dlib.numberOfFunctionCalls;
#endif
        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;

        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)]);
        }

        //        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 << 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 << dlib.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.dlib.numberOfFunctionCalls == settings2.dlib.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;

      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] = reducedPatternSimulationJobs[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]};
          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(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 = 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());
              }

              //              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 << baseTriangleFullPattern.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 (!baseTriangleFullPattern.load(
                  std::filesystem::path(loadFromPath).append(fullPatternLabel + ".ply").string())) {
              !baseTriangleFullPattern.load(std::filesystem::path(loadFromPath)
                                                .append(loadFromPath.stem().string() + ".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"));

              const std::vector<std::filesystem::path> scenariosSortedByName = [&]() {
                  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);
                  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, fullPatternYoungsModulus);
                  fullPatternSimulationJobs.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);
                  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_reducedModel_nonFanned(
          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_reducedModel_nonFanned(optimalXVariables.at("HexSize"),
                                                              optimalXVariables.at("HexAngle"),
                                                              patternBaseTriangle,
                                                              reducedPattern);
              return;
          }
          applyOptimizationResults_reducedModel_nonFanned(optimalXVariables.at("R"),
                                                          optimalXVariables.at("Theta"),
                                                          patternBaseTriangle,
                                                          reducedPattern);
      }

      template<typename MeshType>
      static void applyOptimizationResults_reducedModel_nonFanned(
          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";
          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 < pReducedPattern_simulationMesh->EN(); ei++) {
                  Element &e = pReducedPattern_simulationMesh->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 < pReducedPattern_simulationMesh->EN(); ei++) {
                  Element &e = pReducedPattern_simulationMesh->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 < pReducedPattern_simulationMesh->EN(); ei++) {
                  Element &e = pReducedPattern_simulationMesh->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 < pReducedPattern_simulationMesh->EN(); ei++) {
                  Element &e = pReducedPattern_simulationMesh->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 < pReducedPattern_simulationMesh->EN(); ei++) {
                  Element &e = pReducedPattern_simulationMesh->elements[ei];
                  e.dimensions.inertia.I3 = I3;
              }
          }
          pReducedPattern_simulationMesh->reset();
      }

#if POLYSCOPE_DEFINED
      void draw(const std::vector<int> &desiredSimulationScenariosIndices = std::vector<int>()) 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();
          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> &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 };

  inline bool Settings::ObjectiveWeights::operator==(const ObjectiveWeights &other) const
  {
      return this->translational == other.translational && this->rotational == other.rotational;
  }
  }    // namespace ReducedModelOptimization
#endif // REDUCEDMODELOPTIMIZER_STRUCTS_HPP