#include "reducedmodeloptimizer.hpp"
#include "flatpattern.hpp"
#include "simulationhistoryplotter.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>
#include <dlib/global_optimization.h>
#include <dlib/optimization.h>

struct GlobalOptimizationVariables {
  std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
  std::vector<std::vector<Vector6d>> fullPatternDisplacements;
  std::vector<double> fullPatternDisplacementNormSum;
  std::vector<SimulationJob> g_fullPatternSimulationJob;
  std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
  std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
      reducedToFullInterfaceViMap;
  matplot::line_handle gPlotHandle;
  std::vector<double> gObjectiveValueHistory;
  Eigen::Vector2d g_initialX;
  std::unordered_set<size_t> reducedPatternExludedEdges;
  Eigen::VectorXd g_initialParameters;
  std::vector<ReducedModelOptimizer::SimulationScenario>
      simulationScenarioIndices;
  std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
  double g_innerHexagonInitialPos{0};
  bool optimizeInnerHexagonSize{false};
  std::vector<SimulationResults> firstOptimizationRoundResults;
  int g_firstRoundIterationIndex{0};
  double minY{DBL_MAX};
  std::vector<double> minX;
  std::vector<std::vector<double>> failedSimulationsXRatio;
  int numOfSimulationCrashes{false};
  int numberOfFunctionCalls{0};
  int numberOfOptimizationParameters{3};
  ReducedModelOptimizer::Settings optimizationSettings;
};

// static GlobalOptimizationVariables global;

const static int MAX_THREAD = 64;
#if defined(_MSC_VER)
__declspec(align(64)) GlobalOptimizationVariables tls[MAX_THREAD];
#elif defined(__GNUC__)
GlobalOptimizationVariables tls[MAX_THREAD] __attribute__((aligned(64)));
#endif
//#pragma omp threadprivate(global)

// struct OptimizationCallback {
//  double operator()(const size_t &iterations, const Eigen::VectorXd &x,
//                    const double &fval, Eigen::VectorXd &gradient) const {
//    // run simulation
//    //    SimulationResults reducedModelResults =
//    //        simulator.executeSimulation(reducedModelSimulationJob);
//    //    reducedModelResults.draw(reducedModelSimulationJob);
//    gObjectiveValueHistory.push_back(fval);
//    auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
//                                   gObjectiveValueHistory.size());
//    gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
//    //    const std::string plotImageFilename = "objectivePlot.png";
//    //    matplot::save(plotImageFilename);
//    //    if (numberOfOptimizationRounds % 30 == 0) {
//    //      std::filesystem::copy_file(
//    //          std::filesystem::path(plotImageFilename),
//    //          std::filesystem::path("objectivePlot_copy.png"));
//    //    }
//    //    std::stringstream ss;
//    //    ss << x;
//    //    reducedModelResults.simulationLabel = ss.str();
//    //    SimulationResultsReporter resultsReporter;
//    //    resultsReporter.reportResults(
//    //        {reducedModelResults},
//    //        std::filesystem::current_path().append("Results"));

//    return true;
//  }
//};

// struct Objective {
//  double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const {
//    assert(x.rows() == 4);

//    //      drawSimulationJob(simulationJob);
//    // Set mesh from x
//    std::shared_ptr<SimulationMesh> reducedModel =
//        g_reducedPatternSimulationJob.mesh;
//    for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) {
//      if (g_reducedPatternExludedEdges.contains(ei)) {
//        continue;
//      }
//      Element &e = reducedModel->elements[ei];
//      e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * x(0);
//      e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * x(1);
//      e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * x(2);
//      e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * x(3);
//    }
//    // run simulation
//    SimulationResults reducedModelResults =
//        simulator.executeSimulation(g_reducedPatternSimulationJob);
//    //    std::stringstream ss;
//    //    ss << x;
//    //    reducedModelResults.simulationLabel = ss.str();
//    //    SimulationResultsReporter resultsReporter;
//    //    resultsReporter.reportResults(
//    //        {reducedModelResults},
//    //        std::filesystem::current_path().append("Results"));
//    // compute error and return it
//    double error = 0;
//    for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
//      VertexIndex reducedModelVi = reducedFullViPair.first;
//      Eigen::Vector3d vertexDisplacement(
//          reducedModelResults.displacements[reducedModelVi][0],
//          reducedModelResults.displacements[reducedModelVi][1],
//          reducedModelResults.displacements[reducedModelVi][2]);
//      Eigen::Vector3d errorVector =
//          Eigen::Vector3d(
//              g_optimalReducedModelDisplacements.row(reducedModelVi)) -
//          vertexDisplacement;
//      error += errorVector.norm();
//    }

//    return error;
//  }
//};

double ReducedModelOptimizer::computeError(
    const std::vector<Vector6d> &reducedPatternDisplacements,
    const std::vector<Vector6d> &fullPatternDisplacements,
    const double &interfaceDisplacementsNormSum,
    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap) {
  auto &global = tls[omp_get_thread_num()];
  double error = 0;
  for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
    VertexIndex reducedModelVi = reducedFullViPair.first;
    //    const auto pos =
    //        g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
    //    std::cout << "Interface vi " << reducedModelVi << " is at position "
    //              << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
    Eigen::Vector3d reducedVertexDisplacement(
        reducedPatternDisplacements[reducedModelVi][0],
        reducedPatternDisplacements[reducedModelVi][1],
        reducedPatternDisplacements[reducedModelVi][2]);
    if (!std::isfinite(reducedVertexDisplacement[0]) ||
        !std::isfinite(reducedVertexDisplacement[1]) ||
        !std::isfinite(reducedVertexDisplacement[2])) {
      std::terminate();
    }
    Eigen::Vector3d fullVertexDisplacement(
        fullPatternDisplacements[reducedFullViPair.second][0],
        fullPatternDisplacements[reducedFullViPair.second][1],
        fullPatternDisplacements[reducedFullViPair.second][2]);
    Eigen::Vector3d errorVector =
        fullVertexDisplacement - reducedVertexDisplacement;
    //    error += errorVector.squaredNorm();
    error += errorVector.norm();
  }

  if (global.optimizationSettings.normalizeObjectiveValue) {
    return error / std::max(interfaceDisplacementsNormSum, 0.00003);
  }
  return error;
}

void updateMesh(long n, const double *x) {
  auto &global = tls[omp_get_thread_num()];
  std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
      global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
          ->pMesh;
  //  const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0];
  //  std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
  //  "
  //            << elem.firstBendingConstFactor << " "
  //            << elem.secondBendingConstFactor << std::endl;
  for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
    Element &e = pReducedPatternSimulationMesh->elements[ei];
    //    if (g_reducedPatternExludedEdges.contains(ei)) {
    //      continue;
    //    }
    //    e.properties.E = g_initialParameters * x[ei];
    //    e.properties.E = g_initialParameters(0) * x[0];
    //    e.properties.G = g_initialParameters(1) * x[1];
    e.setDimensions(
        RectangularBeamDimensions(global.g_initialParameters(0) * x[0],
                                  global.g_initialParameters(0) * x[0] /
                                      (global.g_initialParameters(1) * x[1])));
    e.setMaterial(ElementMaterial(e.material.poissonsRatio,
                                  global.g_initialParameters(2) * x[2]));
    //    e.properties.A = g_initialParameters(0) * x[0];
    //    e.properties.J = g_initialParameters(1) * x[1];
    //    e.properties.I2 = g_initialParameters(2) * x[2];
    //    e.properties.I3 = g_initialParameters(3) * x[3];
    //    e.properties.G = e.properties.E / (2 * (1 + 0.3));
    //    e.axialConstFactor = e.properties.E * e.properties.A /
    //    e.initialLength; e.torsionConstFactor = e.properties.G *
    //    e.properties.J / e.initialLength; e.firstBendingConstFactor =
    //        2 * e.properties.E * e.properties.I2 / e.initialLength;
    //    e.secondBendingConstFactor =
    //        2 * e.properties.E * e.properties.I3 / e.initialLength;
  }

  //  std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
  //  "
  //            << elem.firstBendingConstFactor << " "
  //            << elem.secondBendingConstFactor << std::endl;
  //  const Element &e = pReducedPatternSimulationMesh->elements[0];
  //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
  //            << e.firstBendingConstFactor << " " <<
  //            e.secondBendingConstFactor
  //            << std::endl;

  if (global.optimizeInnerHexagonSize) {
    assert(pReducedPatternSimulationMesh->EN() == 12);
    for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
         vi += 2) {
      pReducedPatternSimulationMesh->vert[vi].P() =
          global.g_innerHexagonVectors[vi / 2] * x[n - 1];
    }

    pReducedPatternSimulationMesh->reset();
    pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
    //    pReducedPatternSimulationMesh->registerForDrawing("Optimized
    //    hexagon"); polyscope::show();
  }
}

double ReducedModelOptimizer::objective(double b, double h, double E) {
  std::vector<double> x{b, h, E};
  return ReducedModelOptimizer::objective(x.size(), x.data());
}

double ReducedModelOptimizer::objective(double b, double h, double E,
                                        double innerHexagonSize) {
  std::vector<double> x{b, h, E, innerHexagonSize};
  return ReducedModelOptimizer::objective(x.size(), x.data());
}

double ReducedModelOptimizer::objective(long n, const double *x) {
  auto &global = tls[omp_get_thread_num()];
  //  std::cout.precision(17);

  //  for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
  //    std::cout << "x[" + std::to_string(parameterIndex) + "]="
  //              << x[parameterIndex] << 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);
  //  std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
  //            << e.firstBendingConstFactor << " " <<
  //            e.secondBendingConstFactor
  //            << std::endl;

  // run simulations
  double error = 0;
  FormFinder simulator;
  FormFinder::Settings simulationSettings;
  //    simulationSettings.shouldDraw = true;
  for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
    SimulationResults reducedModelResults = simulator.executeSimulation(
        global.reducedPatternSimulationJobs[simulationScenarioIndex],
        simulationSettings);
    std::string filename;
    if (!reducedModelResults.converged /*&&
        g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]
                    ->pMesh->elements[0]
                    .A > 1e-8 &
            x[0] / x[1] < 60*/) {
      std::cout << "Failed simulation" << std::endl;
      error += std::numeric_limits<double>::max();
      filename = "/home/iason/Coding/Projects/Approximating shapes with flat "
                 "patterns/RodModelOptimizationForPatterns/build/"
                 "ProblematicSimulationJobs/nonConv_dimensions.txt";
      //      if (failedSimulationsXRatio.empty()) {
      //        failedSimulationsXRatio.resize(2);
      //      }
      //      failedSimulationsXRatio[0].push_back(std::log(x[0] / x[1]));
      //      failedSimulationsXRatio[1].push_back(
      //          std::log(g_reducedPatternSimulationJob[g_simulationScenarioIndices[0]]
      //                       ->pMesh->elements[0]
      //                       .A));

      //      SimulationResultsReporter::createPlot(
      //          "log(b/h)", "log(A)", failedSimulationsXRatio[0],
      //          failedSimulationsXRatio[1], "ratioToAPlot.png");
      //      std::cout << "Failed simulation" << std::endl;
      //      simulationSettings.shouldDraw = true;
      //      simulationSettings.debugMessages = true;
      //      simulator.executeSimulation(
      //          g_reducedPatternSimulationJob[simulationScenarioIndex],
      //          simulationSettings);
      global.numOfSimulationCrashes++;
    } else {
      error += computeError(
          reducedModelResults.displacements,
          global.fullPatternDisplacements[simulationScenarioIndex],
          global.fullPatternDisplacementNormSum[simulationScenarioIndex],
          global.reducedToFullInterfaceViMap);
      filename = "/home/iason/Coding/Projects/Approximating shapes with flat "
                 "patterns/RodModelOptimizationForPatterns/build/"
                 "ProblematicSimulationJobs/conv_dimensions.txt";
    }
    std::ofstream out(filename, std::ios_base::app);
    auto pMesh =
        global
            .reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
            ->pMesh;

    for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
      out << "x[" + std::to_string(parameterIndex) + "]=" << x[parameterIndex]
          << std::endl;
    }
    out << pMesh->elements[0].dimensions.toString() + "\n" +
               pMesh->elements[0].material.toString() + " \nA="
        << pMesh->elements[0].A << " \nratio="
        << pMesh->elements[0].dimensions.b / pMesh->elements[0].dimensions.h
        << " \naxialRig:" << pMesh->elements[0].rigidity.axial
        << " \ntorsionalRig:" << pMesh->elements[0].rigidity.torsional
        << " \nfirstBendingRig:" << pMesh->elements[0].rigidity.firstBending
        << " \nsecondBendingRig:" << pMesh->elements[0].rigidity.secondBending
        << " \nscenario:" + simulationScenarioStrings[simulationScenarioIndex] +
               "\n\n";
    out.close();
  }
  //  std::cout << error << std::endl;
  if (error < global.minY) {
    global.minY = error;
    global.minX.assign(x, x + n);
  }
  // if (++global.numberOfFunctionCalls %100== 0) {
  // std::cout << "Number of function calls:" << global.numberOfFunctionCalls
  //          << std::endl;
  //}

  // compute error and return it
  global.gObjectiveValueHistory.push_back(error);
  //  auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
  //                                 gObjectiveValueHistory.size());
  //  std::vector<double> colors(gObjectiveValueHistory.size(), 2);
  //  if (g_firstRoundIterationIndex != 0) {
  //    for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
  //             [](double &c) { c = 0.7; });
  //  }
  //  gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors);
  //  SimulationResultsReporter::createPlot("Number of Steps", "Objective
  //  value",
  //                                        gObjectiveValueHistory);

  return error;
}

void ReducedModelOptimizer::createSimulationMeshes(
    FlatPattern &fullModel, FlatPattern &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(FlatPattern &fullModel,
                                                   FlatPattern &reducedModel) {
  ReducedModelOptimizer::createSimulationMeshes(
      fullModel, reducedModel, m_pFullPatternSimulationMesh,
      m_pReducedPatternSimulationMesh);
}

void ReducedModelOptimizer::computeMaps(
    const std::unordered_set<size_t> &reducedModelExcludedEdges,
    const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
    FlatPattern &fullPattern, FlatPattern &reducedPattern,
    std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap,
    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
        &fullToReducedInterfaceViMap,
    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
        &fullPatternOppositeInterfaceViMap) {

  auto &global = tls[omp_get_thread_num()];
  // 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<FlatPattern>::PointerUpdater<FlatPattern::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 fanSize = 6;
  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<FlatPattern>::PointerUpdater<FlatPattern::VertexPointer>
      pu_reducedModel;
  reducedPattern.deleteDanglingVertices(pu_reducedModel);
  const size_t reducedModelBaseTriangleInterfaceVi =
      pu_reducedModel.remap[baseTriangleInterfaceVi];
  const size_t reducedModelInterfaceVertexOffset =
      reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/;
  reducedPattern.createFan();
  for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
    reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
                                reducedModelBaseTriangleInterfaceVi] =
        fullModelBaseTriangleInterfaceVi +
        fanIndex * fullPatternInterfaceVertexOffset;
  }
  fullToReducedInterfaceViMap.clear();
  constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap);

  // Create opposite vertex map
  fullPatternOppositeInterfaceViMap.clear();
  for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
    const size_t vi0 = fullModelBaseTriangleInterfaceVi +
                       fanIndex * fullPatternInterfaceVertexOffset;
    const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
    assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
    fullPatternOppositeInterfaceViMap[vi0] = vi1;
  }

  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 :
         fullPatternOppositeInterfaceViMap) {
      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();
  }
}

void ReducedModelOptimizer::computeMaps(
    FlatPattern &fullPattern, FlatPattern &reducedPattern,
    const std::unordered_set<size_t> &reducedModelExcludedEdges) {
  auto &global = tls[omp_get_thread_num()];
  ReducedModelOptimizer::computeMaps(
      reducedModelExcludedEdges, slotToNode, fullPattern, reducedPattern,
      global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap,
      m_fullPatternOppositeInterfaceViMap);
}

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

void ReducedModelOptimizer::initializePatterns(
    FlatPattern &fullPattern, FlatPattern &reducedPattern,
    const std::unordered_set<size_t> &reducedModelExcludedEdges) {
  // fullPattern.setLabel("full_pattern_" + fullPattern.getLabel());
  // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
  assert(fullPattern.VN() == reducedPattern.VN() &&
         fullPattern.EN() >= reducedPattern.EN());
  polyscope::removeAllStructures();
  // Create copies of the input models
  FlatPattern copyFullPattern;
  FlatPattern copyReducedPattern;
  copyFullPattern.copy(fullPattern);
  copyReducedPattern.copy(reducedPattern);
  auto &global = tls[omp_get_thread_num()];
  global.optimizeInnerHexagonSize = copyReducedPattern.EN() == 2;
  if (global.optimizeInnerHexagonSize) {
    const double h = copyReducedPattern.getBaseTriangleHeight();
    double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
    VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
                                           0);

    const int fanSize = 6;
    const CoordType rotationAxis(0, 0, 1);
    for (int rotationCounter = 0; rotationCounter < fanSize;
         rotationCounter++) {
      VectorType rotatedVector =
          vcg::RotationMatrix(rotationAxis,
                              vcg::math::ToRad(rotationCounter * 60.0)) *
          baseTriangle_leftBottomNode;
      global.g_innerHexagonVectors[rotationCounter] = rotatedVector;
    }
    const double innerHexagonInitialPos_x =
        copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0];
    const double innerHexagonInitialPos_y =
        copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
    global.g_innerHexagonInitialPos = innerHexagonInitialPos_x;
  }
  computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
  createSimulationMeshes(copyFullPattern, copyReducedPattern);
  initializeOptimizationParameters(m_pReducedPatternSimulationMesh);
}

void ReducedModelOptimizer::initializeOptimizationParameters(
    const std::shared_ptr<SimulationMesh> &mesh) {
  auto &global = tls[omp_get_thread_num()];
  global.numberOfOptimizationParameters = 3;
  global.g_initialParameters.resize(
      global.optimizeInnerHexagonSize ? ++global.numberOfOptimizationParameters
                                      : global.numberOfOptimizationParameters);
  // Save save the beam stiffnesses
  //  for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
  //    Element &e = pReducedModelElementalMesh->elements[ei];
  //    if (g_reducedPatternExludedEdges.contains(ei)) {
  //      const double stiffnessFactor = 5;
  //      e.axialConstFactor *= stiffnessFactor;
  //      e.torsionConstFactor *= stiffnessFactor;
  //      e.firstBendingConstFactor *= stiffnessFactor;
  //      e.secondBendingConstFactor *= stiffnessFactor;
  //    }
  const double initialB = std::sqrt(mesh->elements[0].A);
  const double initialRatio = 1;
  ;
  global.g_initialParameters(0) = initialB;
  global.g_initialParameters(1) = initialRatio;
  global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus;
  if (global.optimizeInnerHexagonSize) {
    global.g_initialParameters(3) = global.g_innerHexagonInitialPos;
  }
  //  g_initialParameters =
  //      m_pReducedPatternSimulationMesh->elements[0].properties.E;
  //  for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
  //  }
  //  g_initialParameters(0) = mesh->elements[0].properties.E;
  //  g_initialParameters(1) = mesh->elements[0].properties.G;
  //  g_initialParameters(0) = mesh->elements[0].properties.A;
  //  g_initialParameters(1) = mesh->elements[0].properties.J;
  //  g_initialParameters(2) = mesh->elements[0].properties.I2;
  //  g_initialParameters(3) = mesh->elements[0].properties.I3;
  //  }
}

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

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

  //  std::unordered_map<VertexIndex, VectorType>
  //  reducedModelNodalForcedNormals; for (auto fullModelNodalForcedRotation :
  //       simulationJobOfFullModel.nodalForcedNormals) {
  //    reducedModelNodalForcedNormals[simulationJobFullToReducedMap.at(
  //        fullModelNodalForcedRotation.first)] =
  //        fullModelNodalForcedRotation.second;
  //  }
  simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices;
  simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces;
  simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel();
  //  simulationJobOfReducedModel.nodalForcedNormals =
  //      reducedModelNodalForcedNormals;
}

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

ReducedModelOptimizer::Results
ReducedModelOptimizer::runOptimization(const Settings &settings) {
  auto &global = tls[omp_get_thread_num()];

  global.gObjectiveValueHistory.clear();

  // g_optimizeInnerHexagonSize ? 5: 4;
  //  const size_t npt = (n + 1) * (n + 2) / 2;
  //  //      ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
  //  assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2);
  //  assert(npt <= 2 * n + 1 && "The choice of the number of interpolation "
  //                             "conditions is not recommended.");
  // Set initial guess of solution

  // const size_t initialGuess = 1;
  // std::vector<double> x(n, initialGuess);
  // if (global.optimizeInnerHexagonSize) {
  //  x[n - 1] = global.g_innerHexagonInitialPos;
  //}
  /*if (!initialGuess.empty()) {
    x = g_optimizationInitialGuess;
  }*/ //      {0.10000000000000	001,
  //  2, 1.9999999971613847, 6.9560343643347764};
  //      {1, 5.9277};
  //      {0.0001, 2, 2.000000005047502, 1.3055270196964464};
  //      {initialGuess(0), initialGuess(1), initialGuess(2),
  //      initialGuess(3)};
  //  assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
  //           return d >= xMax || d <= xMin;
  //         }));
  //  std::vector<double> xLow(x.size(), xMin);
  //  std::vector<double> xUpper(x.size(), xMax);
  //  if (g_optimizeInnerHexagonSize) {
  //    xLow[n - 1] = 0.1;
  //    xUpper[n - 1] = 0.9;
  //  }
  //  const double maxX = *std::max_element(
  //      x.begin(), x.end(),
  //      [](const double &a, const double &b) { return abs(a) < abs(b); });
  //  const double rhobeg = std::min(0.95, 0.2 * maxX);
  //  double rhobeg = 1;
  //  double rhoend = rhobeg * 1e-8;
  //  const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
  //  std::vector<double> w(wSize);
  //  const size_t maxFun = std::min(100.0 * (x.size() + 1), 1000.0);
  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;
  }

  auto start = std::chrono::system_clock::now();
  dlib::function_evaluation result;
  if (global.optimizeInnerHexagonSize) {
    double (*objF)(double, double, double, double) = &objective;
    result = dlib::find_min_global(
        objF, xMin, xMax,
        dlib::max_function_calls(settings.numberOfFunctionCalls),
        std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
  } else {
    double (*objF)(double, double, double) = &objective;
    result = dlib::find_min_global(
        objF, xMin, xMax,
        dlib::max_function_calls(settings.numberOfFunctionCalls),
        std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
  }
  auto end = std::chrono::system_clock::now();
  auto elapsed =
      std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
  Results results;
  results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
  results.x = global.minX;
  results.objectiveValue = global.minY;
  results.time = elapsed.count() / 1000.0;
  if (printDebugInfo) {
    std::cout << "Finished optimizing." << endl;
    //  std::cout << "Solution x:" << endl;
    //  std::cout << result.x << endl;
    std::cout << "Objective value:" << global.minY << endl;
  }
  //  std::cout << result.y << endl;
  //  std::cout << minY << endl;
  //  std::cout << "Time(sec):" << elapsed.count() << std::endl;
  //  std::cout << "Max function evaluations:" << maxFun << std::endl;
  //  std::cout << "Initial guess:" << initialGuess << std::endl;
  //  const size_t maxFun = 200;
  //  bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
  //         rhobeg, rhoend, maxFun, w.data());
  //  std::cout << "Finished first optimization round" << std::endl;
  //  firstOptimizationRoundResults.resize(6);
  //  for (int simulationScenarioIndex = SimulationScenario::Axial;
  //       simulationScenarioIndex !=
  //       SimulationScenario::NumberOfSimulationScenarios;
  //       simulationScenarioIndex++) {
  //    SimulationResults reducedModelResults = simulator.executeSimulation(
  //        g_reducedPatternSimulationJob[simulationScenarioIndex], false,
  //        false);
  //    reducedModelResults.setLabelPrefix("FirstRound");
  //    firstOptimizationRoundResults[simulationScenarioIndex] =
  //        std::move(reducedModelResults);
  //  }
  //  g_firstRoundIterationIndex = gObjectiveValueHistory.size();
  //  rhobeg *= 1e1;
  //  //  rhoend *= 1e2;
  //  bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
  //         rhobeg, rhoend, maxFun, w.data());
  //  std::cout << "Finished second optimization round" << std::endl;

  return results;
}

void ReducedModelOptimizer::setInitialGuess(std::vector<double> v) {
  initialGuess = v;
}

std::vector<std::shared_ptr<SimulationJob>>
ReducedModelOptimizer::createScenarios(
    const std::shared_ptr<SimulationMesh> &pMesh) {
  std::vector<std::shared_ptr<SimulationJob>> scenarios;
  scenarios.resize(SimulationScenario::NumberOfSimulationScenarios);
  std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
  std::unordered_map<VertexIndex, Vector6d> nodalForces;
  const double forceMagnitude = 1;
  // Assuming the patterns lays on the x-y plane
  const CoordType patternPlaneNormal(0, 0, 1);
  // Make the first interface node lay on the x axis
  //  const size_t fullPatternFirstInterfaceNodeIndex =
  //      m_fullPatternOppositeInterfaceViMap.begin()->second;
  //  CoordType fullPatternFirstInterfaceNodePosition =
  //      m_pFullModelSimulationMesh->vert[fullPatternFirstInterfaceNodeIndex].cP();
  //  CoordType centerOfMass(0, 0, 0);
  //  for (size_t vi = 0; vi < pMesh->VN(); vi++) {
  //    centerOfMass = centerOfMass + pMesh->vert[vi].P();
  //  }
  //  centerOfMass /= pMesh->VN();
  //  vcg::tri::UpdatePosition<SimulationMesh>::Translate(
  //      *m_pFullModelSimulationMesh, -centerOfMass);
  //  vcg::tri::UpdatePosition<SimulationMesh>::Translate(
  //      *m_pReducedPatternSimulationMesh, centerOfMass);

  //  const vcg::Matrix33d R = vcg::RotationMatrix(
  //      fullPatternFirstInterfaceNodePosition,
  //      CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false);
  //  std::for_each(m_pFullModelSimulationMesh->vert.begin(),
  //                m_pFullModelSimulationMesh->vert.end(), [&](auto &v) {
  //                  v.P() = R * v.P();
  //                  v.N() = R * v.N();
  //                });
  //  std::for_each(m_pReducedPatternSimulationMesh->vert.begin(),
  //                m_pReducedPatternSimulationMesh->vert.end(), [&](auto &v) {
  //                  v.P() = R * v.P();
  //                  v.N() = R * v.N();
  //                });

  //  m_pFullModelSimulationMesh->updateEigenEdgeAndVertices();
  //  m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();

  //// Axial
  SimulationScenario scenarioName = SimulationScenario::Axial;
  for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
    CoordType forceDirection =
        (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
            .Normalize();
    nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1],
                                          forceDirection[2], 0, 0, 0}) *
                                forceMagnitude * 10;
    fixedVertices[viPair.second] =
        std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
  }
  scenarios[scenarioName] = std::make_shared<SimulationJob>(
      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
                    fixedVertices, nodalForces, {}));

  //// Shear
  scenarioName = SimulationScenario::Shear;
  fixedVertices.clear();
  nodalForces.clear();
  for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
    CoordType v =
        (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
            .Normalize();
    CoordType forceDirection = (v ^ patternPlaneNormal).Normalize();
    nodalForces[viPair.first] = Vector6d({forceDirection[0], forceDirection[1],
                                          forceDirection[2], 0, 0, 0}) *
                                0.40 * forceMagnitude;
    fixedVertices[viPair.second] =
        std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
  }
  scenarios[scenarioName] = std::make_shared<SimulationJob>(
      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
                    fixedVertices, nodalForces, {}));

  //  //// Torsion
  //  fixedVertices.clear();
  //  nodalForces.clear();
  //  for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
  //       viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
  //    const auto &viPair = *viPairIt;
  //    if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
  //      CoordType v =
  //          (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
  //              .Normalize();
  //      CoordType normalVDerivativeDir = (v ^ patternPlaneNormal).Normalize();
  //      nodalForces[viPair.first] = Vector6d{
  //          0, 0, 0, normalVDerivativeDir[0], normalVDerivativeDir[1], 0};
  //      fixedVertices[viPair.second] =
  //          std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
  //      fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
  //    } else {
  //      fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
  //      fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
  //    }
  //  }
  //  scenarios.push_back({pMesh, fixedVertices, nodalForces});

  ////  Bending
  scenarioName = SimulationScenario::Bending;
  fixedVertices.clear();
  nodalForces.clear();
  for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
    nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 1;
    fixedVertices[viPair.second] =
        std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
  }
  scenarios[scenarioName] = std::make_shared<SimulationJob>(
      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
                    fixedVertices, nodalForces, {}));

  //// Double using moments
  scenarioName = SimulationScenario::Dome;
  fixedVertices.clear();
  nodalForces.clear();
  for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
       viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
    const auto viPair = *viPairIt;
    if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
      fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
      fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 2};
    } else {
      fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
      fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
    }
    CoordType v =
        (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
            .Normalize();
    nodalForces[viPair.first] =
        Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 0.1;
    nodalForces[viPair.second] =
        Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1;
  }
  scenarios[scenarioName] = std::make_shared<SimulationJob>(
      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
                    fixedVertices, nodalForces, {}));

  //// Saddle
  scenarioName = SimulationScenario::Saddle;
  fixedVertices.clear();
  nodalForces.clear();
  for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
       viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
    const auto &viPair = *viPairIt;
    CoordType v =
        (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
            .Normalize();
    if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
      nodalForces[viPair.first] =
          Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.02 * forceMagnitude;
      nodalForces[viPair.second] =
          Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.02 * forceMagnitude;
    } else {
      fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
      fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};

      nodalForces[viPair.first] =
          Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.01 * forceMagnitude;
      nodalForces[viPair.second] =
          Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.01 * forceMagnitude;
    }
  }
  scenarios[scenarioName] = std::make_shared<SimulationJob>(
      SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
                    fixedVertices, nodalForces, {}));

  return scenarios;
}

// void ReducedModelOptimizer::runBeamOptimization() {
//  // load beams
//  VCGEdgeMesh fullBeam;
//  fullBeam.loadPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
//  VCGEdgeMesh reducedBeam;
//  reducedBeam.loadPly("/home/iason/Models/simple_beam_model_4elem_1m.ply");
//  fullBeam.registerForDrawing();
//  reducedBeam.registerForDrawing();
//  //  polyscope::show();
//  // maps
//  std::unordered_map<size_t, size_t> displacementReducedToFullMap;
//  displacementReducedToFullMap[reducedBeam.VN() / 2] = fullBeam.VN() / 2;
//  g_reducedToFullViMap = displacementReducedToFullMap;
//  std::unordered_map<size_t, size_t> jobFullToReducedMap;
//  jobFullToReducedMap[0] = 0;
//  jobFullToReducedMap[fullBeam.VN() - 1] = reducedBeam.VN() - 1;

//  // full model simuilation job
//  auto pFullPatternSimulationMesh =
//  std::make_shared<SimulationMesh>(fullBeam);
//  pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.02,
//  0.02}); pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
//  std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices;
//  fixedVertices[0] = ::unordered_set<int>({0, 1, 2, 3, 4, 5});
//  std::unordered_map<VertexIndex, Vector6d> forces;
//  forces[fullBeam.VN() - 1] = Vector6d({0, 0, 10, 0, 0, 0});
//  const std::string fullBeamSimulationJobLabel = "Pull_Z";
//  std::shared_ptr<SimulationJob> pFullModelSimulationJob =
//      make_shared<SimulationJob>(SimulationJob(pFullPatternSimulationMesh,
//                                               fullBeamSimulationJobLabel,
//                                               fixedVertices, forces));
//  auto fullModelResults =
//  formFinder.executeSimulation(pFullModelSimulationJob);

//  // Optimal reduced model displacements
//  const size_t numberOfSimulationScenarios = 1;
//  g_optimalReducedModelDisplacements.resize(numberOfSimulationScenarios);
//  g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1].resize(
//      reducedBeam.VN(), 3);
//  computeDesiredReducedModelDisplacements(
//      fullModelResults, displacementReducedToFullMap,
//      g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);

//  // reduced model simuilation job
//  auto reducedSimulationMesh = std::make_shared<SimulationMesh>(reducedBeam);
//  reducedSimulationMesh->setBeamCrossSection(CrossSectionType{0.02, 0.02});
//  reducedSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
//  g_reducedPatternSimulationJob.resize(numberOfSimulationScenarios);
//  SimulationJob reducedSimJob;
//  computeReducedModelSimulationJob(*pFullModelSimulationJob,
//                                   jobFullToReducedMap, reducedSimJob);
//  reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] =
//      reducedSimJob.nodalExternalForces[reducedBeam.VN() - 1] * 0.1;
//  g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1] =
//      make_shared<SimulationJob>(
//          reducedSimulationMesh, fullBeamSimulationJobLabel,
//          reducedSimJob.constrainedVertices,
//          reducedSimJob.nodalExternalForces);
//  initializeOptimizationParameters(reducedSimulationMesh);

//  //  const std::string simulationJobsPath = "SimulationJobs";
//  //  std::filesystem::create_directory(simulationJobsPath);
//  //  g_reducedPatternSimulationJob[0].save(simulationJobsPath);
//  //  g_reducedPatternSimulationJob[0].load(
//  //      std::filesystem::path(simulationJobsPath)
//  //          .append(g_reducedPatternSimulationJob[0].mesh->getLabel() +
//  //                  "_simScenario.json"));

//  runOptimization({}, &objective);

//  fullModelResults.registerForDrawing();
//  SimulationResults reducedModelResults = simulator.executeSimulation(
//      g_reducedPatternSimulationJob[numberOfSimulationScenarios - 1]);
//  double error = computeError(
//      reducedModelResults,
//      g_optimalReducedModelDisplacements[numberOfSimulationScenarios - 1]);
//  reducedModelResults.registerForDrawing();
//  std::cout << "Error between beams:" << error << endl;
//  //    registerWorldAxes();
//  polyscope::show();
//  fullModelResults.unregister();
//  reducedModelResults.unregister();
//}

void ReducedModelOptimizer::visualizeResults(
    const std::vector<std::shared_ptr<SimulationJob>>
        &fullPatternSimulationJobs,
    const std::vector<std::shared_ptr<SimulationJob>>
        &reducedPatternSimulationJobs,
    const std::vector<SimulationScenario> &simulationScenarios,
    const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        &reducedToFullInterfaceViMap) {
  FormFinder simulator;
  std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
      fullPatternSimulationJobs[0]->pMesh;
  pFullPatternSimulationMesh->registerForDrawing();
  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 interfaceDisplacementNormSum = 0;
    for (const auto &interfaceViPair : reducedToFullInterfaceViMap) {
      const int fullPatternInterfaceIndex = interfaceViPair.second;
      Eigen::Vector3d fullPatternDisplacementVector(
          fullModelResults.displacements[fullPatternInterfaceIndex][0],
          fullModelResults.displacements[fullPatternInterfaceIndex][1],
          fullModelResults.displacements[fullPatternInterfaceIndex][2]);
      interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
    }
    double error = computeError(
        reducedModelResults.displacements, fullModelResults.displacements,
        interfaceDisplacementNormSum, reducedToFullInterfaceViMap);
    std::cout << "Error of simulation scenario "
              << simulationScenarioStrings[simulationScenarioIndex] << " is "
              << error << std::endl;
    totalError += error;
    reducedModelResults.registerForDrawing();
    //    firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
    //    reducedModelResults.saveDeformedModel();
    //    registerWorldAxes();
    const std::string screenshotFilename =
        "/home/iason/Coding/Projects/Approximating shapes with flat "
        "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
        "Images/" +
        pFullPatternSimulationMesh->getLabel() + "_" +
        simulationScenarioStrings[simulationScenarioIndex];
    polyscope::show();
    polyscope::screenshot(screenshotFilename, false);
    fullModelResults.unregister();
    reducedModelResults.unregister();
    //    firstOptimizationRoundResults[simulationScenarioIndex].unregister();
  }
  std::cout << "Total error:" << totalError << std::endl;
}

ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
    const Settings &optimizationSettings,
    const std::vector<SimulationScenario> &simulationScenarios) {
  auto &global = tls[omp_get_thread_num()];

  global.simulationScenarioIndices = simulationScenarios;
  if (global.simulationScenarioIndices.empty()) {
    global.simulationScenarioIndices = {
        SimulationScenario::Axial, SimulationScenario::Shear,
        SimulationScenario::Bending, SimulationScenario::Dome,
        SimulationScenario::Saddle};
  }

  std::vector<std::shared_ptr<SimulationJob>> simulationJobs =
      createScenarios(m_pFullPatternSimulationMesh);
  global.g_optimalReducedModelDisplacements.resize(6);
  global.reducedPatternSimulationJobs.resize(6);
  global.fullPatternDisplacements.resize(6);
  global.fullPatternDisplacementNormSum.resize(6);
  global.g_firstRoundIterationIndex = 0;
  global.minY = std::numeric_limits<double>::max();
  global.numOfSimulationCrashes = 0;
  global.numberOfFunctionCalls = 0;
  global.optimizationSettings = optimizationSettings;
  //  polyscope::removeAllStructures();

  FormFinder::Settings settings;
  //  settings.shouldDraw = true;
  for (int simulationScenarioIndex : global.simulationScenarioIndices) {
    const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
        simulationJobs[simulationScenarioIndex];
    SimulationResults fullModelResults =
        simulator.executeSimulation(pFullPatternSimulationJob, settings);
    global.fullPatternDisplacements[simulationScenarioIndex] =
        fullModelResults.displacements;
    double interfaceDisplacementNormSum = 0;
    for (const auto &interfaceViPair : global.reducedToFullInterfaceViMap) {
      const int fullPatternInterfaceIndex = interfaceViPair.second;
      Eigen::Vector3d fullPatternDisplacementVector(
          fullModelResults.displacements[fullPatternInterfaceIndex][0],
          fullModelResults.displacements[fullPatternInterfaceIndex][1],
          fullModelResults.displacements[fullPatternInterfaceIndex][2]);
      interfaceDisplacementNormSum += fullPatternDisplacementVector.norm();
    }
    global.fullPatternDisplacementNormSum[simulationScenarioIndex] =
        interfaceDisplacementNormSum;
    // global.g_optimalReducedModelDisplacements[simulationScenarioIndex].resize(
    //    m_pReducedPatternSimulationMesh->VN(), 3);
    // computeDesiredReducedModelDisplacements(
    //    fullModelResults, global.reducedToFullInterfaceViMap,
    //    global.g_optimalReducedModelDisplacements[simulationScenarioIndex]);

    SimulationJob reducedPatternSimulationJob;
    reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
    computeReducedModelSimulationJob(*pFullPatternSimulationJob,
                                     m_fullToReducedInterfaceViMap,
                                     reducedPatternSimulationJob);
    global.reducedPatternSimulationJobs[simulationScenarioIndex] =
        std::make_shared<SimulationJob>(reducedPatternSimulationJob);
  }
  Results optResults = runOptimization(optimizationSettings);
  for (int simulationScenarioIndex : global.simulationScenarioIndices) {
    optResults.fullPatternSimulationJobs.push_back(
        simulationJobs[simulationScenarioIndex]);
    optResults.reducedPatternSimulationJobs.push_back(
        global.reducedPatternSimulationJobs[simulationScenarioIndex]);
  }
  //  updateMesh(optResults.x.size(), optResults.x.data());
  //  optResults.draw();

  //  visualizeResults(simulationJobs, global.simulationScenarioIndices);
  // visualizeResults(simulationJobs, global.reducedPatternSimulationJobs,
  //                 global.simulationScenarioIndices,global.reducedToFullInterfaceViMap);
  return optResults;
}