#include "reducedmodeloptimizer.hpp"
#include "linearsimulationmodel.hpp"
#include "simulationhistoryplotter.hpp"
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>

using namespace ReducedPatternOptimization;

struct GlobalOptimizationVariables {
    //    std::vector<std::vector<Vector6d>> fullPatternDisplacements;
    std::vector<SimulationResults> fullPatternResults;
    std::vector<double> translationalDisplacementNormalizationValues;
    std::vector<double> rotationalDisplacementNormalizationValues;
    std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
    std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
    std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        reducedToFullInterfaceViMap;
    std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
        fullPatternInterfaceViPairs;
    matplot::line_handle gPlotHandle;
    std::vector<double> objectiveValueHistory;
    Eigen::VectorXd initialParameters;
    std::vector<int> simulationScenarioIndices;
    double minY{DBL_MAX};
    std::vector<double> minX;
    int numOfSimulationCrashes{false};
    int numberOfFunctionCalls{0};
    int numberOfOptimizationParameters{5};
    ReducedPatternOptimization::Settings optimizationSettings;
    vcg::Triangle3<double> baseTriangle;
    //Variables for finding the full pattern simulation forces
    std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh;
    std::function<void(const double &,
                       const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
                       SimulationJob &)>
        constructScenarioFunction;
    FullPatternVertexIndex interfaceViForComputingScenarioError;
    double desiredMaxDisplacementValue;
    double desiredMaxRotationAngle;
} global;

double ReducedModelOptimizer::computeDisplacementError(
    const std::vector<Vector6d> &fullPatternDisplacements,
    const std::vector<Vector6d> &reducedPatternDisplacements,
    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap,
    const double &normalizationFactor)
{
    const double rawError = computeRawTranslationalError(fullPatternDisplacements,
                                                         reducedPatternDisplacements,
                                                         reducedToFullInterfaceViMap);
    return rawError / normalizationFactor;
}

double ReducedModelOptimizer::computeRawTranslationalError(
    const std::vector<Vector6d> &fullPatternDisplacements,
    const std::vector<Vector6d> &reducedPatternDisplacements,
    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap)
{
    double error = 0;
    for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
        const VertexIndex reducedModelVi = reducedFullViPair.first;
        const VertexIndex fullModelVi = reducedFullViPair.second;
        const Eigen::Vector3d fullPatternVertexDiplacement = fullPatternDisplacements[fullModelVi]
                                                                 .getTranslation();
        const Eigen::Vector3d reducedPatternVertexDiplacement
            = reducedPatternDisplacements[reducedModelVi].getTranslation();
        const double vertexError = (fullPatternVertexDiplacement - reducedPatternVertexDiplacement)
                                       .norm();
        error += vertexError;
    }

    return error;
}

double ReducedModelOptimizer::computeRawRotationalError(
    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap)
{
    double rawRotationalError = 0;
    for (const auto &reducedToFullInterfaceViPair : reducedToFullInterfaceViMap) {
        const double vertexRotationalError
            = rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance(
                rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first]);
        rawRotationalError += vertexRotationalError;
    }

    return rawRotationalError;
}

double ReducedModelOptimizer::computeRotationalError(
    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
    const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap,
    const double &normalizationFactor)
{
    const double rawRotationalError = computeRawRotationalError(rotatedQuaternion_fullPattern,
                                                                rotatedQuaternion_reducedPattern,
                                                                reducedToFullInterfaceViMap);
    return rawRotationalError / normalizationFactor;
}

double ReducedModelOptimizer::computeError(
    const SimulationResults &simulationResults_fullPattern,
    const SimulationResults &simulationResults_reducedPattern,
    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap,
    const double &normalizationFactor_translationalDisplacement,
    const double &normalizationFactor_rotationalDisplacement)
{
    const double translationalError
        = computeDisplacementError(simulationResults_fullPattern.displacements,
                                   simulationResults_reducedPattern.displacements,
                                   reducedToFullInterfaceViMap,
                                   normalizationFactor_translationalDisplacement);
    const double rotationalError
        = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
                                 simulationResults_reducedPattern.rotationalDisplacementQuaternion,
                                 reducedToFullInterfaceViMap,
                                 normalizationFactor_rotationalDisplacement);
    return global.optimizationSettings.objectiveWeights.translational * translationalError
           + global.optimizationSettings.objectiveWeights.rotational * rotationalError;
}

double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
                                        double innerHexagonSize,
                                        double innerHexagonRotationAngle) {
    std::vector<double> x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle};
    return ReducedModelOptimizer::objective(x.size(), x.data());
}

double ReducedModelOptimizer::objective(long n, const double *x) {
    //  std::cout.precision(17);
    //        for (int i = 0; i < n; i++) {
    //    std::cout << x[i] << " ";
    //        }
    //    std::cout << std::endl;

    //    std::cout << x[n - 2] << " " << x[n - 1] << std::endl;
    //  const Element &e =
    //  global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout <<
    //  e.axialConstFactor << " " << e.torsionConstFactor << " "
    //            << e.firstBendingConstFactor << " " <<
    //            e.secondBendingConstFactor
    //            << std::endl;
    updateMesh(n, x);
    //    global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
    //    global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
    //    polyscope::show();
    //    global.reducedPatternSimulationJobs[0]->pMesh->unregister();

    //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
    //            << e.firstBendingConstFactor << " " <<
    //            e.secondBendingConstFactor
    //            << std::endl;

    // run simulations
    double totalError = 0;
    LinearSimulationModel simulator;
    //    FormFinder simulator;
    for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
        const std::shared_ptr<SimulationJob> &reducedJob
            = global.reducedPatternSimulationJobs[simulationScenarioIndex];
        SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob);
        //    std::string filename;
        if (!reducedModelResults.converged) {
            totalError += std::numeric_limits<double>::max();
            global.numOfSimulationCrashes++;
#ifdef POLYSCOPE_DEFINED
            std::cout << "Failed simulation" << std::endl;
#endif
        } else {
            const bool usePotentialEnergy = false;
            double simulationScenarioError;
            if (usePotentialEnergy) {
                simulationScenarioError = std::abs(
                    reducedModelResults.internalPotentialEnergy
                    - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
            } else {
                simulationScenarioError = computeError(
                    global.fullPatternResults[simulationScenarioIndex],
                    reducedModelResults,
                    global.reducedToFullInterfaceViMap,
                    global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
                    global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
            }
            //#ifdef POLYSCOPE_DEFINED
            //            std::cout << "sim error:" << simulationScenarioError << std::endl;
            //            global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(
            //                ReducedModelOptimization::Colors::reducedInitial);
            //            reducedModelResults.registerForDrawing(
            //                ReducedModelOptimization::Colors::reducedDeformed);
            //            polyscope::show();
            //            reducedModelResults.unregister();
            //#endif
            //        if (global.optimizationSettings.normalizationStrategy !=
            //        NormalizationStrategy::Epsilon &&
            //            simulationScenarioError > 1) {
            //          std::cout << "Simulation scenario "
            //                    <<
            //                    simulationScenarioStrings[simulationScenarioIndex]
            //                    << " results in an error bigger than one." <<
            //                    std::endl;
            //          for (size_t parameterIndex = 0; parameterIndex < n;
            //               parameterIndex++) {
            //            std::cout << "x[" + std::to_string(parameterIndex) + "]="
            //                      << x[parameterIndex] << std::endl;
            //          }
            //        }

            //#ifdef POLYSCOPE_DEFINED
            //        ReducedModelOptimizer::visualizeResults(
            //            global.fullPatternSimulationJobs[simulationScenarioIndex],
            //            global.reducedPatternSimulationJobs[simulationScenarioIndex],
            //            global.reducedToFullInterfaceViMap, false);

            //        ReducedModelOptimizer::visualizeResults(
            //            global.fullPatternSimulationJobs[simulationScenarioIndex],
            //            std::make_shared<SimulationJob>(
            //                reducedPatternMaximumDisplacementSimulationJobs
            //                    [simulationScenarioIndex]),
            //            global.reducedToFullInterfaceViMap, true);
            //        polyscope::removeAllStructures();
            //#endif // POLYSCOPE_DEFINED
            totalError += simulationScenarioError;
        }
    }
    //  std::cout << error << std::endl;
    if (totalError < global.minY) {
        global.minY = totalError;
        global.minX.assign(x, x + n);
    }
#ifdef POLYSCOPE_DEFINED
    ++global.numberOfFunctionCalls;
    if (global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100)
        == 0) {
        std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl;
    }
#endif
//     compute error and return it
//      global.objectiveValueHistory.push_back(totalError);
//      auto xPlot = matplot::linspace(0, global.objectiveValueHistory.size(),
//                                     global.objectiveValueHistory.size());
//      std::vector<double> colors(global.gObjectiveValueHistory.size(), 2);
//      if (global.g_firstRoundIterationIndex != 0) {
//        for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
//                 [](double &c) { c = 0.7; });
//      }
//      global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory);
//      SimulationResultsReporter::createPlot("Number of Steps", "Objective value",
//                                            global.objectiveValueHistory);

    return totalError;
}

void ReducedModelOptimizer::createSimulationMeshes(
    PatternGeometry &fullModel, PatternGeometry &reducedModel,
    std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
    if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
        std::cerr << "Error: A rectangular cross section is expected." << std::endl;
        terminate();
    }
    // Full pattern
    pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
    pFullPatternSimulationMesh->setBeamCrossSection(
        CrossSectionType{0.002, 0.002});
    pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);

    // Reduced pattern
    pReducedPatternSimulationMesh =
        std::make_shared<SimulationMesh>(reducedModel);
    pReducedPatternSimulationMesh->setBeamCrossSection(
        CrossSectionType{0.002, 0.002});
    pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
}

void ReducedModelOptimizer::createSimulationMeshes(
    PatternGeometry &fullModel, PatternGeometry &reducedModel) {
    ReducedModelOptimizer::createSimulationMeshes(
        fullModel, reducedModel, m_pFullPatternSimulationMesh,
        m_pReducedPatternSimulationMesh);
}

void ReducedModelOptimizer::computeMaps(
    const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
    PatternGeometry &fullPattern,
    PatternGeometry &reducedPattern,
    std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap,
    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
        &fullToReducedInterfaceViMap,
    std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
        &fullPatternOppositeInterfaceViPairs)
{
    // Compute the offset between the interface nodes
    const size_t interfaceSlotIndex = 4; // bottom edge
    assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() &&
           slotToNode.find(interfaceSlotIndex)->second.size() == 1);
    // Assuming that in the bottom edge there is only one vertex which is also the
    // interface
    const size_t baseTriangleInterfaceVi =
        *(slotToNode.find(interfaceSlotIndex)->second.begin());

    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<PatternGeometry::VertexPointer> pu_fullModel;
    fullPattern.deleteDanglingVertices(pu_fullModel);
    const size_t fullModelBaseTriangleInterfaceVi =
        pu_fullModel.remap.empty() ? baseTriangleInterfaceVi
                                   : pu_fullModel.remap[baseTriangleInterfaceVi];
    const size_t fullModelBaseTriangleVN = fullPattern.VN();
    fullPattern.createFan();
    const size_t duplicateVerticesPerFanPair =
        fullModelBaseTriangleVN - fullPattern.VN() / 6;
    const size_t fullPatternInterfaceVertexOffset =
        fullModelBaseTriangleVN - duplicateVerticesPerFanPair;
    //  std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair <<
    //  std::endl;

    // Save excluded edges TODO:this changes the global object. Should this be a
    // function parameter?
    //  global.reducedPatternExludedEdges.clear();
    //  const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
    //  for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
    //    for (const size_t ei : reducedModelExcludedEdges) {
    //      global.reducedPatternExludedEdges.insert(
    //          fanIndex * reducedBaseTriangleNumberOfEdges + ei);
    //    }
    //  }

    // Construct reduced->full and full->reduced interface vi map
    reducedToFullInterfaceViMap.clear();
    vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<PatternGeometry::VertexPointer> pu_reducedModel;
    reducedPattern.deleteDanglingVertices(pu_reducedModel);
    const size_t reducedModelBaseTriangleInterfaceVi =
        pu_reducedModel.remap[baseTriangleInterfaceVi];
    const size_t reducedModelInterfaceVertexOffset
        = reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/;
    Results::applyOptimizationResults_innerHexagon(initialHexagonSize,
                                                   0,
                                                   global.baseTriangle,
                                                   reducedPattern);
    reducedPattern.createFan({0}); //TODO: should be an input parameter from main
    for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
        reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
                                    reducedModelBaseTriangleInterfaceVi] =
                fullModelBaseTriangleInterfaceVi +
                fanIndex * fullPatternInterfaceVertexOffset;
    }
    fullToReducedInterfaceViMap.clear();
    constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap);

    // Create opposite vertex map
    fullPatternOppositeInterfaceViPairs.clear();
    fullPatternOppositeInterfaceViPairs.reserve(fanSize / 2);
    //    for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
    for (int fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) {
        const size_t vi0 = fullModelBaseTriangleInterfaceVi
                           + fanIndex * fullPatternInterfaceVertexOffset;
        const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
        assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
        fullPatternOppositeInterfaceViPairs.emplace_back(std::make_pair(vi0, vi1));
    }
    global.fullPatternInterfaceViPairs = fullPatternOppositeInterfaceViPairs;

#if POLYSCOPE_DEFINED
    const bool debugMapping = false;
    if (debugMapping) {
        reducedPattern.registerForDrawing();

        //    std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
        //        reducedPattern.EN(), glm::vec3(0, 0, 0));
        //    for (const size_t ei : global.reducedPatternExludedEdges) {
        //      colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
        //    }
        //    const std::string label = reducedPattern.getLabel();
        //    polyscope::getCurveNetwork(label)
        //        ->addEdgeColorQuantity("Excluded edges",
        //                               colors_reducedPatternExcludedEdges)
        //        ->setEnabled(true);
        //    polyscope::show();

        std::vector<glm::vec3> nodeColorsOpposite(fullPattern.VN(),
                                                  glm::vec3(0, 0, 0));
        for (const std::pair<size_t, size_t> oppositeVerts : fullPatternOppositeInterfaceViPairs) {
            auto color = polyscope::getNextUniqueColor();
            nodeColorsOpposite[oppositeVerts.first] = color;
            nodeColorsOpposite[oppositeVerts.second] = color;
        }
        fullPattern.registerForDrawing();
        polyscope::getCurveNetwork(fullPattern.getLabel())
            ->addNodeColorQuantity("oppositeMap", nodeColorsOpposite)
            ->setEnabled(true);
        polyscope::show();

        std::vector<glm::vec3> nodeColorsReducedToFull_reduced(reducedPattern.VN(),
                                                               glm::vec3(0, 0, 0));
        std::vector<glm::vec3> nodeColorsReducedToFull_full(fullPattern.VN(),
                                                            glm::vec3(0, 0, 0));
        for (size_t vi = 0; vi < reducedPattern.VN(); vi++) {
            if (global.reducedToFullInterfaceViMap.contains(vi)) {

                auto color = polyscope::getNextUniqueColor();
                nodeColorsReducedToFull_reduced[vi] = color;
                nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] =
                    color;
            }
        }
        polyscope::getCurveNetwork(reducedPattern.getLabel())
            ->addNodeColorQuantity("reducedToFull_reduced",
                                   nodeColorsReducedToFull_reduced)
            ->setEnabled(true);
        polyscope::getCurveNetwork(fullPattern.getLabel())
            ->addNodeColorQuantity("reducedToFull_full",
                                   nodeColorsReducedToFull_full)
            ->setEnabled(true);
        polyscope::show();
    }
#endif
}

void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern,
                                        PatternGeometry &reducedPattern)
{
    ReducedModelOptimizer::computeMaps(slotToNode,
                                       fullPattern,
                                       reducedPattern,
                                       global.reducedToFullInterfaceViMap,
                                       m_fullToReducedInterfaceViMap,
                                       m_fullPatternOppositeInterfaceViPairs);
}

ReducedModelOptimizer::ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot)
{
    FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot);
    FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode);
}

void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
                                               PatternGeometry &reducedPattern,
                                               const int &optimizationParameters)
{
    assert(fullPattern.VN() == reducedPattern.VN() &&
           fullPattern.EN() >= reducedPattern.EN());
#if POLYSCOPE_DEFINED
    polyscope::removeAllStructures();
#endif
    // Create copies of the input models
    PatternGeometry copyFullPattern;
    PatternGeometry copyReducedPattern;
    copyFullPattern.copy(fullPattern);
    copyReducedPattern.copy(reducedPattern);
    global.baseTriangle = copyReducedPattern.getBaseTriangle();

    computeMaps(copyFullPattern, copyReducedPattern);
    createSimulationMeshes(copyFullPattern, copyReducedPattern);
    initializeOptimizationParameters(m_pReducedPatternSimulationMesh,optimizationParameters);
}

void updateMesh(long n, const double *x) {
    std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
        global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
            ->pMesh;

    const double E=global.initialParameters(0)*x[0];
    const double A=global.initialParameters(1) * x[1];
    const double beamWidth=std::sqrt(A);
        const double beamHeight=beamWidth;

const double J=global.initialParameters(2) * x[2];
const double I2=global.initialParameters(3) * x[3];
const double I3=global.initialParameters(4) * x[4];
    for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
        Element &e = pReducedPatternSimulationMesh->elements[ei];
        e.setDimensions(
            RectangularBeamDimensions(beamWidth,
                                      beamHeight));
        e.setMaterial(ElementMaterial(e.material.poissonsRatio,
                                      E));
            e.J = J;
            e.I2 = I2;
            e.I3 = I3;
    }
    assert(pReducedPatternSimulationMesh->EN() == 12);
    assert(n >= 2);

    //    CoordType center_barycentric(1, 0, 0);
    //    CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
    //    CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
    //                                        * x[n - 2]);

    CoordType movableVertex_barycentric(1 - x[n - 2], x[n - 2] / 2, x[n - 2] / 2);
    CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
                                                      * movableVertex_barycentric[0]
                                                  + global.baseTriangle.cP(1)
                                                        * movableVertex_barycentric[1]
                                                  + global.baseTriangle.cP(2)
                                                        * movableVertex_barycentric[2];
    baseTriangleMovableVertexPosition
        = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, vcg::math::ToRad(x[n - 1]))
          * baseTriangleMovableVertexPosition;

    for (int rotationCounter = 0;
         rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
        pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
            = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
                                  vcg::math::ToRad(60.0 * rotationCounter))
              * baseTriangleMovableVertexPosition;
    }

    pReducedPatternSimulationMesh->reset();
    //#ifdef POLYSCOPE_DEFINED
    //    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
    //    pReducedPatternSimulationMesh->registerForDrawing();
    //    std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2])
    //              << std::endl;
    //    std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl;
    //    polyscope::show();
    //#endif
}


void ReducedModelOptimizer::initializeOptimizationParameters(
    const std::shared_ptr<SimulationMesh> &mesh,const int& optimizationParamters) {
    global.numberOfOptimizationParameters = optimizationParamters;
    global.initialParameters.resize(global.numberOfOptimizationParameters);

      global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
      global.initialParameters(1) = mesh->elements[0].A;
      global.initialParameters(2) = mesh->elements[0].J;
      global.initialParameters(3) = mesh->elements[0].I2;
      global.initialParameters(4) = mesh->elements[0].I3;
}

void ReducedModelOptimizer::computeReducedModelSimulationJob(
    const SimulationJob &simulationJobOfFullModel,
    const std::unordered_map<size_t, size_t> &fullToReducedMap,
    SimulationJob &simulationJobOfReducedModel)
{
    assert(simulationJobOfReducedModel.pMesh->VN() != 0);
    std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
        reducedModelFixedVertices;
    for (auto fullModelFixedVertex :
         simulationJobOfFullModel.constrainedVertices) {
        reducedModelFixedVertices[fullToReducedMap.at(fullModelFixedVertex.first)]
            = fullModelFixedVertex.second;
    }

    std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
    for (auto fullModelNodalForce :
         simulationJobOfFullModel.nodalExternalForces) {
        reducedModelNodalForces[fullToReducedMap.at(fullModelNodalForce.first)]
            = fullModelNodalForce.second;
    }

    std::unordered_map<VertexIndex, Eigen::Vector3d> reducedNodalForcedDisplacements;
    for (auto fullForcedDisplacement : simulationJobOfFullModel.nodalForcedDisplacements) {
        reducedNodalForcedDisplacements[fullToReducedMap.at(fullForcedDisplacement.first)]
            = fullForcedDisplacement.second;
    }

    simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices;
    simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces;
    simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel();
    simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements;
}

//#if POLYSCOPE_DEFINED
//void ReducedModelOptimizer::visualizeResults(
//    const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
//    const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
//    const std::vector<BaseSimulationScenario> &simulationScenarios,
//    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
//        &reducedToFullInterfaceViMap)
//{
//    DRMSimulationModel simulator;
//    std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
//        fullPatternSimulationJobs[0]->pMesh;
//    pFullPatternSimulationMesh->registerForDrawing();
//    pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply");
//    reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel()
//                                                 + "_undeformed.ply");
//    double totalError = 0;
//    for (const int simulationScenarioIndex : simulationScenarios) {
//        const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
//            fullPatternSimulationJobs[simulationScenarioIndex];
//        pFullPatternSimulationJob->registerForDrawing(
//            pFullPatternSimulationMesh->getLabel());
//        SimulationResults fullModelResults =
//            simulator.executeSimulation(pFullPatternSimulationJob);
//        fullModelResults.registerForDrawing();
//        // fullModelResults.saveDeformedModel();
//        const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
//            reducedPatternSimulationJobs[simulationScenarioIndex];
//        SimulationResults reducedModelResults =
//            simulator.executeSimulation(pReducedPatternSimulationJob);
//        double normalizationFactor = 1;
//        if (global.optimizationSettings.normalizationStrategy !=
//            ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) {
//            normalizationFactor
//                = global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
//        }
//        reducedModelResults.saveDeformedModel();
//        fullModelResults.saveDeformedModel();
//        double error = computeDisplacementError(reducedModelResults.displacements,
//                                                fullModelResults.displacements,
//                                                reducedToFullInterfaceViMap,
//                                                normalizationFactor);
//        std::cout << "Error of simulation scenario "
//                  <baseSimulationScenarioNames[simulationScenarioIndex]
//                  << " is " << error << std::endl;
//        totalError += error;
//        reducedModelResults.registerForDrawing();
//        //    firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
//        //    registerWorldAxes();
//        const std::string screenshotFilename
//            = "/home/iason/Coding/Projects/Approximating shapes with flat "
//              "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
//              "Images/"
//              + pFullPatternSimulationMesh->getLabel() + "_"
//              + baseSimulationScenarioNames[simulationScenarioIndex];
//        polyscope::show();
//        polyscope::screenshot(screenshotFilename, false);
//        fullModelResults.unregister();
//        reducedModelResults.unregister();
//        //    firstOptimizationRoundResults[simulationScenarioIndex].unregister();
//    }
//    std::cout << "Total error:" << totalError << std::endl;
//}

//void ReducedModelOptimizer::registerResultsForDrawing(
//    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
//    const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
//    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
//        &reducedToFullInterfaceViMap) {
//    DRMSimulationModel simulator;
//    std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
//        pFullPatternSimulationJob->pMesh;
//    pFullPatternSimulationMesh->registerForDrawing();
//    //  pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel()
//    //  +
//    //                                      "_undeformed.ply");
//    //  reducedPatternSimulationJobs[0]->pMesh->savePly(
//    //      reducedPatternSimulationJobs[0]->pMesh->getLabel() +
//    //      "_undeformed.ply");
//    pFullPatternSimulationJob->registerForDrawing(
//        pFullPatternSimulationMesh->getLabel());
//    SimulationResults fullModelResults =
//        simulator.executeSimulation(pFullPatternSimulationJob);
//    fullModelResults.registerForDrawing();
//    // fullModelResults.saveDeformedModel();
//    SimulationResults reducedModelResults =
//        simulator.executeSimulation(pReducedPatternSimulationJob);
//    //    reducedModelResults.saveDeformedModel();
//    //    fullModelResults.saveDeformedModel();
//    double error = computeRawDisplacementError(
//        reducedModelResults.displacements, fullModelResults.displacements,
//        reducedToFullInterfaceViMap/*,
//      global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/);
//    std::cout << "Error is " << error << std::endl;
//    reducedModelResults.registerForDrawing();
//}
//#endif // POLYSCOPE_DEFINED

void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
    const SimulationResults &fullModelResults,
    const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
    Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) {

    assert(optimalDisplacementsOfReducedModel.rows() != 0 &&
           optimalDisplacementsOfReducedModel.cols() == 3);
    optimalDisplacementsOfReducedModel.setZero(
        optimalDisplacementsOfReducedModel.rows(),
        optimalDisplacementsOfReducedModel.cols());

    for (auto reducedFullViPair : displacementsReducedToFullMap) {
        const VertexIndex fullModelVi = reducedFullViPair.second;
        const Vector6d fullModelViDisplacements =
            fullModelResults.displacements[fullModelVi];
        optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) =
            Eigen::Vector3d(fullModelViDisplacements[0],
                            fullModelViDisplacements[1],
                            fullModelViDisplacements[2]);
    }
}

void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimizationResult_dlib,
                                       const Settings &settings,
                                       ReducedPatternOptimization::Results &results)
{
    //Number of crashes
    //    results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
    //Value of optimal objective Y
    results.objectiveValue.total = optimizationResult_dlib.y;
    //Optimal X values
    results.optimalXNameValuePairs.resize(settings.xRanges.size());
        std::vector<double> optimalX(settings.xRanges.size());
    for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) {
        if (xVariableIndex < 5) {
            results.optimalXNameValuePairs[xVariableIndex]
                = std::make_pair(settings.xRanges[xVariableIndex].label,
                                 global.minX[xVariableIndex]
                                     * global.initialParameters(xVariableIndex));
        } else {
            //Hex size and angle are pure values (not multipliers of the initial values)
            results.optimalXNameValuePairs[xVariableIndex]
                = std::make_pair(settings.xRanges[xVariableIndex].label,
                                 global.minX[xVariableIndex]);
        }

        assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex));
        optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex);
#ifdef POLYSCOPE_DEFINED
        std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":"
                  << optimalX[xVariableIndex] << " ";
#endif
    }
#ifdef POLYSCOPE_DEFINED
    std::cout << std::endl;
#endif

    // Compute obj value per simulation scenario and the raw objective value
    updateMesh(optimalX.size(), optimalX.data());
    //    results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios);
    //TODO:use push_back it will make the code more readable
    LinearSimulationModel simulator;
    results.objectiveValue.totalRaw = 0;
    results.objectiveValue.perSimulationScenario_translational.resize(
        global.simulationScenarioIndices.size());
    results.objectiveValue.perSimulationScenario_rawTranslational.resize(
        global.simulationScenarioIndices.size());
    results.objectiveValue.perSimulationScenario_rotational.resize(
        global.simulationScenarioIndices.size());
    results.objectiveValue.perSimulationScenario_rawRotational.resize(
        global.simulationScenarioIndices.size());
    results.objectiveValue.perSimulationScenario_total.resize(
        global.simulationScenarioIndices.size());
    for (int i = 0; i < global.simulationScenarioIndices.size(); i++) {
        const int simulationScenarioIndex = global.simulationScenarioIndices[i];
        SimulationResults reducedModelResults = simulator.executeSimulation(
            global.reducedPatternSimulationJobs[simulationScenarioIndex]);

        results.objectiveValue.perSimulationScenario_total[i] = computeError(
            global.fullPatternResults[simulationScenarioIndex],
            reducedModelResults,
            global.reducedToFullInterfaceViMap,
            global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
            global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);

        //Raw translational
        const double rawTranslationalError = computeRawTranslationalError(
            global.fullPatternResults[simulationScenarioIndex].displacements,
            reducedModelResults.displacements,
            global.reducedToFullInterfaceViMap);
        results.objectiveValue.perSimulationScenario_rawTranslational[i] = rawTranslationalError;
        //Raw rotational
        const double rawRotationalError = computeRawRotationalError(
            global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
            reducedModelResults.rotationalDisplacementQuaternion,
            global.reducedToFullInterfaceViMap);
        results.objectiveValue.perSimulationScenario_rawRotational[i] = rawRotationalError;

        //Normalized translational
        const double normalizedTranslationalError = computeDisplacementError(
            global.fullPatternResults[simulationScenarioIndex].displacements,
            reducedModelResults.displacements,
            global.reducedToFullInterfaceViMap,
            global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
        results.objectiveValue.perSimulationScenario_translational[i] = normalizedTranslationalError;
        //            const double test_normalizedTranslationError = computeDisplacementError(
        //                global.fullPatternResults[simulationScenarioIndex].displacements,
        //                reducedModelResults.displacements,
        //                global.reducedToFullInterfaceViMap,
        //                global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
        //Normalized rotational
        const double normalizedRotationalError = computeRotationalError(
            global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
            reducedModelResults.rotationalDisplacementQuaternion,
            global.reducedToFullInterfaceViMap,
            global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
        results.objectiveValue.perSimulationScenario_rotational[i] = normalizedRotationalError;
        //            const double test_normalizedRotationalError = computeRotationalError(
        //                global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
        //                reducedModelResults.rotationalDisplacementQuaternion,
        //                global.reducedToFullInterfaceViMap,
        //                global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
        //            assert(test_normalizedTranslationError == normalizedTranslationalError);
        //            assert(test_normalizedRotationalError == normalizedRotationalError);
        results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError;
#ifdef POLYSCOPE_DEFINED
        std::cout << "Simulation scenario:"
                  << global.reducedPatternSimulationJobs[simulationScenarioIndex]->getLabel()
                  << std::endl;
        std::cout << "raw translational error:" << rawTranslationalError << std::endl;
        std::cout << "translation normalization value:"
                  << global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
                  << std::endl;
        std::cout << "raw Rotational error:" << rawRotationalError << std::endl;
        std::cout << "rotational normalization value:"
                  << global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
                  << std::endl;
        std::cout << "Translational error:" << normalizedTranslationalError << std::endl;
        std::cout << "Rotational error:" << normalizedRotationalError << std::endl;
        //        results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
        //            = normalizedTranslationalError + normalizedRotationalError;
        std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i]
                  << std::endl;
        std::cout << std::endl;
#endif
    }

    const bool printDebugInfo = false;
    if (printDebugInfo) {
        std::cout << "Finished optimizing." << endl;
    std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl;
    assert(global.minY == optimizationResult_dlib.y);
    if (global.minY != optimizationResult_dlib.y) {
        std::cerr << "ERROR in objective value" << std::endl;
    }
    }

    for (int simulationScenarioIndex : global.simulationScenarioIndices) {
        results.fullPatternSimulationJobs.push_back(
            global.fullPatternSimulationJobs[simulationScenarioIndex]);
        results.reducedPatternSimulationJobs.push_back(
            global.reducedPatternSimulationJobs[simulationScenarioIndex]);
        //        const std::string temp = global.reducedPatternSimulationJobs[simulationScenarioIndex]
        //                                     ->pMesh->getLabel();
        //        global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel("temp");
        //        global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing();
        //        global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp);
    }
}

void ReducedModelOptimizer::runOptimization(const Settings &settings,
                                            ReducedPatternOptimization::Results &results)
{
    global.objectiveValueHistory.clear();
    dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
    dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
    for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
        xMin(i) = settings.xRanges[i].min;
        xMax(i) = settings.xRanges[i].max;
    }

    dlib::function_evaluation result_dlib;
    double (*objF)(double, double, double, double, double,double,double) = &objective;
    result_dlib = dlib::find_min_global(objF,
                                        xMin,
                                        xMax,
                                        dlib::max_function_calls(settings.numberOfFunctionCalls),
                                        std::chrono::hours(24 * 365 * 290),
                                        settings.solverAccuracy);
    getResults(result_dlib, settings, results);
}

void ReducedModelOptimizer::constructAxialSimulationScenario(
    const double &forceMagnitude,
    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
        &oppositeInterfaceViPairs,
    SimulationJob &job)
{
    for (auto viPairIt = oppositeInterfaceViPairs.begin();
         viPairIt != oppositeInterfaceViPairs.end();
         viPairIt++) {
        if (viPairIt != oppositeInterfaceViPairs.begin()) {
            CoordType forceDirection(1, 0, 0);
            job.nodalExternalForces[viPairIt->first]
                = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
                  * forceMagnitude;
            job.constrainedVertices[viPairIt->second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
        }
    }

}

void ReducedModelOptimizer::constructShearSimulationScenario(
    const double &forceMagnitude,
    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
        &oppositeInterfaceViPairs,
    SimulationJob &job)
{
    for (auto viPairIt = oppositeInterfaceViPairs.begin();
         viPairIt != oppositeInterfaceViPairs.end();
         viPairIt++) {
        if (viPairIt != oppositeInterfaceViPairs.begin()) {
            CoordType forceDirection(0, 1, 0);
            const auto viPair = *viPairIt;
            job.nodalExternalForces[viPair.first]
                = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
                  * forceMagnitude;
            job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
        }
    }
}

void ReducedModelOptimizer::constructBendingSimulationScenario(
    const double &forceMagnitude,
    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
        &oppositeInterfaceViPairs,
    SimulationJob &job)
{
    for (const auto &viPair : oppositeInterfaceViPairs) {
        job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude;
        job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
    }
}

/*NOTE: From the results it seems as if the forced displacements are different in the linear and in the drm
 * */
void ReducedModelOptimizer::constructDomeSimulationScenario(
    const double &forceMagnitude,
    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
        &oppositeInterfaceViPairs,
    SimulationJob &job)
{
    for (auto viPairIt = oppositeInterfaceViPairs.begin();
         viPairIt != oppositeInterfaceViPairs.end();
         viPairIt++) {
        const auto viPair = *viPairIt;
        CoordType interfaceVector = (job.pMesh->vert[viPair.first].cP()
                                     - job.pMesh->vert[viPair.second].cP());
        VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0))
                                * interfaceVector.Normalize();
        if (viPairIt == oppositeInterfaceViPairs.begin()) {
            job.nodalForcedDisplacements[viPair.first]
                = Eigen::Vector3d(-interfaceVector[0],
                                  -interfaceVector[1],
                                  0)
                  * 0.01
                  * std::abs(
                      forceMagnitude); //NOTE:Should the forced displacement change relatively to the magnitude?
            //                  * std::abs(forceMagnitude / maxForceMagnitude_dome);
            job.nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0],
                                                                          interfaceVector[1],
                                                                          0)
                                                          * 0.001 * std::abs(forceMagnitude);
            //                  * std::abs(forceMagnitude / maxForceMagnitude_dome);
            //            CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
            //                          ^ CoordType(0, 0, -1).Normalize();
            //            nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
            //                                        * 0.0001;
            //            nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
            //                                         * 0.0001;
        } else {
            job.nodalExternalForces[viPair.first]
                = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude
                  / 3;
            job.nodalExternalForces[viPair.second]
                = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]})
                  * forceMagnitude / 3;
            job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
            job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{2};
        }
    }
}

void ReducedModelOptimizer::constructSaddleSimulationScenario(
    const double &forceMagnitude,
    const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
        &oppositeInterfaceViPairs,
    SimulationJob &job)
{
    for (auto viPairIt = oppositeInterfaceViPairs.begin();
         viPairIt != oppositeInterfaceViPairs.end();
         viPairIt++) {
        const auto &viPair = *viPairIt;
        CoordType v = (job.pMesh->vert[viPair.first].cP() - job.pMesh->vert[viPair.second].cP())
                      ^ CoordType(0, 0, -1).Normalize();
        if (viPairIt == oppositeInterfaceViPairs.begin()) {
            job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0})
                                                    * forceMagnitude;
            job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0})
                                                     * forceMagnitude;
        } else {
            job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
            job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};

            job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0})
                                                    * forceMagnitude / 2;
            job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0})
                                                     * forceMagnitude / 2;
        }
    }
}

double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude)
{
    SimulationJob job;
    job.pMesh = global.pFullPatternSimulationMesh;
    global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job);
    //    ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
    //                                                   global.fullPatternInterfaceViPairs,
    //                                                   job);

    DRMSimulationModel simulator;
    DRMSimulationModel::Settings settings;
    //    settings.totalResidualForcesNormThreshold = 1e-3;
    //        settings.totalTranslationalKineticEnergyThreshold = 1e-6;
    //        settings.shouldUseTranslationalKineticEnergyThreshold = true;
    //    settings.shouldDraw = true;
    //    settings.isDebugMode = true;
    //    settings.drawingStep = 500;
    //    settings.beVerbose = true;
    //    settings.debugModeStep = 200000;
    //    settings.shouldCreatePlots = true;
    settings.maxDRMIterations = 100000;
    SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
                                                            settings);
    const double &desiredRotationAngle = global.desiredMaxRotationAngle;
    if (!results.converged) {
        return std::numeric_limits<double>::max();
    }
    const double error = std::abs(
        //        results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm()
        results.rotationalDisplacementQuaternion[global.interfaceViForComputingScenarioError]
            .angularDistance(Eigen::Quaterniond::Identity())
        - desiredRotationAngle);

#ifdef POLYSCOPE_DEFINED
    std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl;
#endif
    return error;
}

double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude)
{
    SimulationJob job;
    job.pMesh = global.pFullPatternSimulationMesh;
    global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job);
    //    ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
    //                                                   global.fullPatternInterfaceViPairs,
    //                                                   job);

    DRMSimulationModel simulator;
    DRMSimulationModel::Settings settings;
    settings.totalResidualForcesNormThreshold = 1e-3;
    //    settings.totalTranslationalKineticEnergyThreshold = 1e-10;
    //    settings.shouldUseTranslationalKineticEnergyThreshold = true;
    //    settings.shouldDraw = true;
    //    settings.isDebugMode = true;
    //    settings.drawingStep = 200000;
    //    settings.beVerbose = true;
    //    settings.debugModeStep = 200000;
    settings.maxDRMIterations = 50000;
#ifdef POLYSCOPE_DEFINED
//    settings.shouldDraw = true;
//    settings.isDebugMode = true;
//    settings.drawingStep = 100000;
//    settings.debugModeStep = 10000;
//    settings.shouldCreatePlots = true;
#endif
    SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job));
    const double &desiredDisplacementValue = global.desiredMaxDisplacementValue;
    if (!results.converged) {
        return std::numeric_limits<double>::max();
    }
    const double error = std::abs(
        //        results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm()
        results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm()
        - desiredDisplacementValue);

#ifdef POLYSCOPE_DEFINED
    std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl;
#endif

    return error;
}

double ReducedModelOptimizer::getFullPatternMaxSimulationForce(
    const BaseSimulationScenario &scenario, bool &wasSuccessful)
{
    double forceMagnitude = 1;
    const double forceMagnitudeEpsilon = 1e-2;
    double minimumError;
    double translationalOptimizationEpsilon;
    dlib::function_evaluation result;
    switch (scenario) {
    case Axial:
        global.desiredMaxDisplacementValue = 0.03
                                             * vcg::Distance(global.baseTriangle.cP(0),
                                                             (global.baseTriangle.cP(1)
                                                              + global.baseTriangle.cP(2))
                                                                 / 2);
        global.constructScenarioFunction = &ReducedModelOptimizer::constructAxialSimulationScenario;
        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
        minimumError
            = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
                                             forceMagnitude,
                                             1e-2,
                                             1e2,
                                             forceMagnitudeEpsilon);
        translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
        wasSuccessful = minimumError < translationalOptimizationEpsilon;
        break;
    case Shear:
        global.desiredMaxDisplacementValue = 0.04
                                             * vcg::Distance(global.baseTriangle.cP(0),
                                                             (global.baseTriangle.cP(1)
                                                              + global.baseTriangle.cP(2))
                                                                 / 2);
        global.constructScenarioFunction = &ReducedModelOptimizer::constructShearSimulationScenario;
        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
        minimumError
            = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
                                             forceMagnitude,
                                             1e-2,
                                             1e2,
                                             forceMagnitudeEpsilon);
        translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
        wasSuccessful = minimumError < translationalOptimizationEpsilon;
        break;
    case Bending:
        global.desiredMaxDisplacementValue = 0.1
                                             * vcg::Distance(global.baseTriangle.cP(0),
                                                             (global.baseTriangle.cP(1)
                                                              + global.baseTriangle.cP(2))
                                                                 / 2);
        global.constructScenarioFunction = &ReducedModelOptimizer::constructBendingSimulationScenario;
        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
        minimumError
            = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
                                             forceMagnitude,
                                             1e-2,
                                             1e2,
                                             forceMagnitudeEpsilon);
        translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
        wasSuccessful = minimumError < translationalOptimizationEpsilon;
        break;
    case Dome:
        global.desiredMaxRotationAngle = vcg::math::ToRad(35.0);
        global.constructScenarioFunction = &ReducedModelOptimizer::constructDomeSimulationScenario;
        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
        minimumError
            = dlib::find_min_single_variable(&fullPatternMaxSimulationForceRotationalObjective,
                                             forceMagnitude,
                                             1e-2,
                                             1e2,
                                             forceMagnitudeEpsilon);
        wasSuccessful = minimumError < vcg::math::ToRad(3.0);
        //        result = dlib::find_min_global(&fullPatternMaxSimulationForceRotationalObjective,
        //                                       1e-2,
        //                                       1,
        //                                       dlib::max_function_calls(50));
        //        wasSuccessful = result.y < vcg::math::ToRad(3.0);
        break;
    case Saddle:
        global.desiredMaxDisplacementValue = 0.1
                                             * vcg::Distance(global.baseTriangle.cP(0),
                                                             (global.baseTriangle.cP(1)
                                                              + global.baseTriangle.cP(2))
                                                                 / 2);
        global.constructScenarioFunction = &ReducedModelOptimizer::constructSaddleSimulationScenario;
        global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
        minimumError
            = dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
                                             forceMagnitude,
                                             1e-2,
                                             1e2,
                                             forceMagnitudeEpsilon);
        translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
        wasSuccessful = minimumError < translationalOptimizationEpsilon;
        break;
    }

#ifdef POLYSCOPE_DEFINED
    std::cout << "Max " << ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
              << " magnitude:" << forceMagnitude << std::endl;
    if (!wasSuccessful)
        std::cout << "Was not successfull" << std::endl;
#endif
    if (!wasSuccessful) {
        const std::string debugMessage
            = ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
              + " max scenario magnitude was not succefully determined.";
        optimizationNotes.append(debugMessage);
    }

    return forceMagnitude;
}

//TODO: Make more compact
std::vector<std::shared_ptr<SimulationJob>>
ReducedModelOptimizer::createFullPatternSimulationScenarios(
    const std::shared_ptr<SimulationMesh> &pMesh)
{
    std::vector<std::shared_ptr<SimulationJob>> scenarios;
    scenarios.resize(totalNumberOfSimulationScenarios);

    SimulationJob job;
    job.pMesh = pMesh;

    //// Axial
    bool wasSuccessful;
    const double maxForceMagnitude_axial
        = getFullPatternMaxSimulationForce(BaseSimulationScenario::Axial, wasSuccessful);
    const double minForceMagnitude_axial = -maxForceMagnitude_axial;
    const int numberOfSimulationScenarios_axial
        = simulationScenariosResolution[BaseSimulationScenario::Axial];
    const double forceMagnitudeStep_axial = numberOfSimulationScenarios_axial == 1
                                                ? maxForceMagnitude_axial
                                                : (maxForceMagnitude_axial - minForceMagnitude_axial)
                                                      / (numberOfSimulationScenarios_axial - 1);
    const int baseSimulationScenarioIndexOffset_axial
        = std::accumulate(simulationScenariosResolution.begin(),
                          simulationScenariosResolution.begin() + BaseSimulationScenario::Axial,
                          0);
    for (int axialSimulationScenarioIndex = 0;
         axialSimulationScenarioIndex < numberOfSimulationScenarios_axial;
         axialSimulationScenarioIndex++) {
        job.nodalExternalForces.clear();
        job.constrainedVertices.clear();
        job.nodalForcedDisplacements.clear();
        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_"
                    + std::to_string(axialSimulationScenarioIndex);

        const double forceMagnitude = (forceMagnitudeStep_axial * axialSimulationScenarioIndex
                                       + minForceMagnitude_axial);
        constructAxialSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);

        scenarios[baseSimulationScenarioIndexOffset_axial + axialSimulationScenarioIndex]
            = std::make_shared<SimulationJob>(job);
    }

    //// Shear
    const double maxForceMagnitude_shear
        = getFullPatternMaxSimulationForce(BaseSimulationScenario::Shear, wasSuccessful);
    const double minForceMagnitude_shear = -maxForceMagnitude_shear;
    const int numberOfSimulationScenarios_shear
        = simulationScenariosResolution[BaseSimulationScenario::Shear];
    const double forceMagnitudeStep_shear = numberOfSimulationScenarios_shear == 1
                                                ? maxForceMagnitude_shear
                                                : (maxForceMagnitude_shear - minForceMagnitude_shear)
                                                      / (numberOfSimulationScenarios_shear - 1);
    const int baseSimulationScenarioIndexOffset_shear
        = std::accumulate(simulationScenariosResolution.begin(),
                          simulationScenariosResolution.begin() + BaseSimulationScenario::Shear,
                          0);
    for (int shearSimulationScenarioIndex = 0;
         shearSimulationScenarioIndex < numberOfSimulationScenarios_shear;
         shearSimulationScenarioIndex++) {
        job.constrainedVertices.clear();
        job.nodalExternalForces.clear();
        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_"
                    + std::to_string(shearSimulationScenarioIndex);
        const double forceMagnitude = (forceMagnitudeStep_shear * shearSimulationScenarioIndex
                                       + minForceMagnitude_shear);
        constructShearSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
        scenarios[baseSimulationScenarioIndexOffset_shear + shearSimulationScenarioIndex]
            = std::make_shared<SimulationJob>(job);
    }

    ////  Bending
    const double maxForceMagnitude_bending
        = getFullPatternMaxSimulationForce(BaseSimulationScenario::Bending, wasSuccessful);
    const double minForceMagnitude_bending = 0;
    const int numberOfSimulationScenarios_bending
        = simulationScenariosResolution[BaseSimulationScenario::Bending];
    const double forceMagnitudeStep_bending = (maxForceMagnitude_bending - minForceMagnitude_bending)
                                              / (numberOfSimulationScenarios_bending);
    const int baseSimulationScenarioIndexOffset_bending
        = std::accumulate(simulationScenariosResolution.begin(),
                          simulationScenariosResolution.begin() + BaseSimulationScenario::Bending,
                          0);
    for (int bendingSimulationScenarioIndex = 0;
         bendingSimulationScenarioIndex < numberOfSimulationScenarios_bending;
         bendingSimulationScenarioIndex++) {
        job.nodalExternalForces.clear();
        job.constrainedVertices.clear();
        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_"
                    + std::to_string(bendingSimulationScenarioIndex);
        const double forceMagnitude = (forceMagnitudeStep_bending
                                           * (bendingSimulationScenarioIndex + 1)
                                       + minForceMagnitude_bending);
        constructBendingSimulationScenario(forceMagnitude,
                                           m_fullPatternOppositeInterfaceViPairs,
                                           job);
        scenarios[baseSimulationScenarioIndexOffset_bending + bendingSimulationScenarioIndex]
            = std::make_shared<SimulationJob>(job);
    }

    //// Dome
    const double maxForceMagnitude_dome
        = getFullPatternMaxSimulationForce(BaseSimulationScenario::Dome, wasSuccessful);
    const double minForceMagnitude_dome = 0;
    const int numberOfSimulationScenarios_dome
        = simulationScenariosResolution[BaseSimulationScenario::Dome];
    const double forceMagnitudeStep_dome = (maxForceMagnitude_dome - minForceMagnitude_dome)
                                           / (numberOfSimulationScenarios_dome);
    const int baseSimulationScenarioIndexOffset_dome
        = std::accumulate(simulationScenariosResolution.begin(),
                          simulationScenariosResolution.begin() + BaseSimulationScenario::Dome,
                          0);
    for (int domeSimulationScenarioIndex = 0;
         domeSimulationScenarioIndex < numberOfSimulationScenarios_dome;
         domeSimulationScenarioIndex++) {
        job.constrainedVertices.clear();
        job.nodalExternalForces.clear();
        job.nodalForcedDisplacements.clear();
        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_"
                    + std::to_string(domeSimulationScenarioIndex);
        const double forceMagnitude = (forceMagnitudeStep_dome * (domeSimulationScenarioIndex + 1)
                                       + minForceMagnitude_dome);
        constructDomeSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
        scenarios[baseSimulationScenarioIndexOffset_dome + domeSimulationScenarioIndex]
            = std::make_shared<SimulationJob>(job);
    }

    //// Saddle
    const double maxForceMagnitude_saddle
        = getFullPatternMaxSimulationForce(BaseSimulationScenario::Saddle, wasSuccessful);
    const double minForceMagnitude_saddle = 0;
    const int numberOfSimulationScenarios_saddle
        = simulationScenariosResolution[BaseSimulationScenario::Saddle];
    const double forceMagnitudeStep_saddle = (maxForceMagnitude_saddle - minForceMagnitude_saddle)
                                             / numberOfSimulationScenarios_saddle;
    const int baseSimulationScenarioIndexOffset_saddle
        = std::accumulate(simulationScenariosResolution.begin(),
                          simulationScenariosResolution.begin() + BaseSimulationScenario::Saddle,
                          0);
    for (int saddleSimulationScenarioIndex = 0;
         saddleSimulationScenarioIndex < numberOfSimulationScenarios_saddle;
         saddleSimulationScenarioIndex++) {
        job.constrainedVertices.clear();
        job.nodalExternalForces.clear();
        job.nodalForcedDisplacements.clear();
        job.label = baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_"
                    + std::to_string(saddleSimulationScenarioIndex);
        const double forceMagnitude = (forceMagnitudeStep_saddle
                                           * (saddleSimulationScenarioIndex + 1)
                                       + minForceMagnitude_saddle);
        constructSaddleSimulationScenario(forceMagnitude,
                                          m_fullPatternOppositeInterfaceViPairs,
                                          job);
        scenarios[baseSimulationScenarioIndexOffset_saddle + saddleSimulationScenarioIndex]
            = std::make_shared<SimulationJob>(job);
    }

#ifdef POLYSCOPE_DEFINED
    std::cout << "Computed full pattern scenario magnitudes" << std::endl;
#endif
    return scenarios;
}

void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {

        // Compute the sum of the displacement norms
        std::vector<double> fullPatternTranslationalDisplacementNormSum(
            totalNumberOfSimulationScenarios);
        std::vector<double> fullPatternAngularDistance(totalNumberOfSimulationScenarios);
        for (int simulationScenarioIndex : global.simulationScenarioIndices) {
            double translationalDisplacementNormSum = 0;
            for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
                const int fullPatternVi = interfaceViPair.second;
                //If the full pattern vertex is translationally constrained dont take it into account
                if (global.fullPatternSimulationJobs[simulationScenarioIndex]
                        ->constrainedVertices.contains(fullPatternVi)) {
                    const std::unordered_set<int> constrainedDof
                        = global.fullPatternSimulationJobs[simulationScenarioIndex]
                              ->constrainedVertices.at(fullPatternVi);
                    if (constrainedDof.contains(0) && constrainedDof.contains(1)
                        && constrainedDof.contains(2)) {
                        continue;
                    }
                }

                const Vector6d &vertexDisplacement = global
                                                         .fullPatternResults[simulationScenarioIndex]
                                                         .displacements[fullPatternVi];
                translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm();
            }
            double angularDistanceSum = 0;
            for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
                const int fullPatternVi = interfaceViPair.second;
                //If the full pattern vertex is rotationally constrained dont take it into account
                if (global.fullPatternSimulationJobs[simulationScenarioIndex]
                        ->constrainedVertices.contains(fullPatternVi)) {
                    const std::unordered_set<int> constrainedDof
                        = global.fullPatternSimulationJobs[simulationScenarioIndex]
                              ->constrainedVertices.at(fullPatternVi);
                    if (constrainedDof.contains(3) && constrainedDof.contains(5)
                        && constrainedDof.contains(4)) {
                        continue;
                    }
                }
                angularDistanceSum += global.fullPatternResults[simulationScenarioIndex]
                                          .rotationalDisplacementQuaternion[fullPatternVi]
                                          .angularDistance(Eigen::Quaterniond::Identity());
            }

            fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex]
                = translationalDisplacementNormSum;
            fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum;
        }

        for (int simulationScenarioIndex : global.simulationScenarioIndices) {
            if (global.optimizationSettings.normalizationStrategy ==
                Settings::NormalizationStrategy::Epsilon) {
                const double epsilon_translationalDisplacement = global.optimizationSettings
                                                                     .normalizationParameter;
                global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
                    = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
                               epsilon_translationalDisplacement);
                //                const double epsilon_rotationalDisplacement = vcg::math::ToRad(10.0);
                global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
                    = /*std::max(*/ fullPatternAngularDistance[simulationScenarioIndex] /*,
                               epsilon_rotationalDisplacement)*/
                    ;
            } else {
                global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
                global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
            }
        }
}

void ReducedModelOptimizer::optimize(
    const Settings &optimizationSettings,
    ReducedPatternOptimization::Results &results,
    const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices)
{
    for (int baseSimulationScenarioIndex : desiredBaseSimulationScenarioIndices) {
        //Increase the size of the vector holding the simulation scenario indices
        global.simulationScenarioIndices.resize(
            global.simulationScenarioIndices.size()
            + simulationScenariosResolution[baseSimulationScenarioIndex]);
        //Add the simulation scenarios indices that correspond to this base simulation scenario
        std::iota(global.simulationScenarioIndices.end()
                      - simulationScenariosResolution[baseSimulationScenarioIndex],
                  global.simulationScenarioIndices.end(),
                  std::accumulate(simulationScenariosResolution.begin(),
                                  simulationScenariosResolution.begin()
                                      + baseSimulationScenarioIndex,
                                  0));
    }
    if (desiredBaseSimulationScenarioIndices.empty()) {
        global.simulationScenarioIndices.resize(totalNumberOfSimulationScenarios);
        std::iota(global.simulationScenarioIndices.begin(),
                  global.simulationScenarioIndices.end(),
                  0);
    }

    global.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios);
    global.fullPatternResults.resize(totalNumberOfSimulationScenarios);
    global.translationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios);
    global.rotationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios);
    global.minY = std::numeric_limits<double>::max();
    global.numOfSimulationCrashes = 0;
    global.numberOfFunctionCalls = 0;
    global.optimizationSettings = optimizationSettings;
    global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh;
    global.fullPatternSimulationJobs = createFullPatternSimulationScenarios(
        m_pFullPatternSimulationMesh);
    //  polyscope::removeAllStructures();

    results.baseTriangle = global.baseTriangle;

    DRMSimulationModel::Settings simulationSettings;
    simulationSettings.maxDRMIterations = 200000;
//    simulationSettings.shouldDraw = true;
//    simulationSettings.isDebugMode = true;
//    simulationSettings.debugModeStep = 100000;
#ifdef POLYSCOPE_DEFINED
    const bool drawFullPatternSimulationResults = false;
    if (drawFullPatternSimulationResults) {
        global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
            ReducedPatternOptimization::Colors::fullInitial);
    }
#endif
    //    LinearSimulationModel linearSimulator;
    for (int simulationScenarioIndex : global.simulationScenarioIndices) {
        const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
            = global.fullPatternSimulationJobs[simulationScenarioIndex];

        SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
                                                                           simulationSettings);
        if (!fullPatternResults.converged) {
            results.wasSuccessful = false;
            return;
        }
        results.wasSuccessful = true;
#ifdef POLYSCOPE_DEFINED
        if (drawFullPatternSimulationResults) {
            //                SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
            //                    pFullPatternSimulationJob);
            fullPatternResults.registerForDrawing(ReducedPatternOptimization::Colors::fullDeformed,
                                                  true,
                                                  true);
            //                fullPatternResults_linear.labelPrefix += "_linear";
            //                fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
            //                                                             true,
            //                                                             true);
            polyscope::show();
            fullPatternResults.unregister();
            //                fullPatternResults_linear.unregister();
        }
#endif
        global.fullPatternResults[simulationScenarioIndex] = fullPatternResults;
        SimulationJob reducedPatternSimulationJob;
        reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
        computeReducedModelSimulationJob(*pFullPatternSimulationJob,
                                         m_fullToReducedInterfaceViMap,
                                         reducedPatternSimulationJob);
        global.reducedPatternSimulationJobs[simulationScenarioIndex]
            = std::make_shared<SimulationJob>(reducedPatternSimulationJob);
        //        std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl;
    }

#ifdef POLYSCOPE_DEFINED
    if (drawFullPatternSimulationResults) {
        global.fullPatternSimulationJobs[0]->pMesh->unregister();
    }
#endif
    //    if (global.optimizationSettings.normalizationStrategy
    //        != Settings::NormalizationStrategy::NonNormalized) {
    computeObjectiveValueNormalizationFactors();
    //    }
#ifdef POLYSCOPE_DEFINED
    std::cout << "Running reduced model optimization" << std::endl;
#endif
    runOptimization(optimizationSettings, results);
    results.notes = optimizationNotes;
}