#ifndef REDUCEDMODELOPTIMIZER_HPP
#define REDUCEDMODELOPTIMIZER_HPP

#include "csvfile.hpp"
#include "drmsimulationmodel.hpp"
#include "edgemesh.hpp"
#include "linearsimulationmodel.hpp"
#ifdef POLYSCOPE_DEFINED
#include "matplot/matplot.h"
#endif
#include "reducedmodel.hpp"
#include "reducedmodeloptimizer_structs.hpp"
#include "simulationmesh.hpp"
#include <Eigen/Dense>
#ifdef DLIB_DEFINED
#include <dlib/global_optimization.h>
#include <dlib/optimization.h>
#endif

#ifdef POLYSCOPE_DEFINED
#include "polyscope/color_management.h"
#endif // POLYSCOPE_DEFINED
using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;

class ReducedModelOptimizer
{
public:
    struct OptimizationState
    {
        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>>
            fullPatternOppositeInterfaceViPairs;
        matplot::line_handle gPlotHandle;
        std::vector<size_t> objectiveValueHistory_iteration;
        std::vector<double> objectiveValueHistory;
        std::vector<double> plotColors;
        std::array<double,
                   ReducedModelOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
            parametersInitialValue;
        std::array<double,
                   ReducedModelOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
            optimizationInitialValue;
        std::vector<int> simulationScenarioIndices;
        double minY{DBL_MAX};
        std::vector<double> minX;
        int numOfSimulationCrashes{false};
        int numberOfFunctionCalls{0};
        //Variables for finding the full pattern simulation forces
        std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh;
        std::array<std::function<void(const double &newValue,
                                      std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>,
                   7>
            functions_updateReducedPatternParameter;
        std::vector<double> xMin;
        std::vector<double> xMax;
        std::vector<double> scenarioWeights;
        std::vector<ReducedModelOptimization::Settings::ObjectiveWeights> objectiveWeights;
    };

private:
    OptimizationState optimizationState;
    vcg::Triangle3<double> baseTriangle;
    std::function<void(const double &,
                       const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
                       SimulationJob &)>
        constructScenarioFunction;
    std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
    std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
        m_fullToReducedInterfaceViMap;
    std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
        m_fullPatternOppositeInterfaceViPairs;
    std::unordered_map<size_t, size_t> nodeToSlot;
    std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
    std::string optimizationNotes;
    std::array<std::function<void(
                   const double &,
                   const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
                   SimulationJob &)>,
               ReducedModelOptimization::NumberOfBaseSimulationScenarios>
        constructBaseScenarioFunctions;
    std::vector<bool> scenarioIsSymmetrical;
    int fullPatternNumberOfEdges;
    constexpr static double youngsModulus{1 * 1e9};
    std::string fullPatternLabel;
    //    ReducedModelOptimization::Settings optimizationSettings;

public:
    struct FunctionEvaluation
    {
        FunctionEvaluation() = default;
        FunctionEvaluation(const std::vector<double> &x, double y) : x(x), y(y) {}

        std::vector<double> x;
        double y = std::numeric_limits<double>::quiet_NaN();
    };

    //    struct ParameterLabels
    //    {
    //        inline const static std::string E = {"E"};
    //        inline const static std::string A = {"A"};
    //        inline const static std::string I2 ={"I2"};
    //        inline const static std::string I3 ={"I3"};
    //        inline const static std::string J = {"J"};
    //        inline const static std::string th= {"Theta"};
    //        inline const static std::string R = {"R"};
    //    };
    //    inline constexpr static ParameterLabels parameterLabels();

    inline static std::array<std::string, ReducedModelOptimization::NumberOfOptimizationVariables>
        parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"};
    constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
        simulationScenariosResolution = {12, 12, 22, 22, 22, 22};
    constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
        baseScenarioWeights = {1, 1, 2, 3, 2};
    inline static int totalNumberOfSimulationScenarios
        = std::accumulate(simulationScenariosResolution.begin(),
                          simulationScenariosResolution.end(),
                          0);
    inline static int fanSize{6};
    inline static double initialHexagonSize{0.3};
    int interfaceNodeIndex;
    double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;

    ReducedModelOptimizer();
    static void computeReducedModelSimulationJob(
        const SimulationJob &simulationJobOfFullModel,
        const std::unordered_map<size_t, size_t> &fullToReducedMap,
        SimulationJob &simulationJobOfReducedModel);

    SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);

    static void runSimulation(const std::string &filename, std::vector<double> &x);

    static std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
        const std::shared_ptr<SimulationMesh> &pMesh,
        const std::unordered_map<size_t, size_t> &fullPatternOppositeInterfaceViMap);

    static void createSimulationMeshes(
        PatternGeometry &fullModel,
        PatternGeometry &reducedModel,
        std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
        std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
    void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
                     PatternGeometry &fullPattern,
                     ReducedModel &reducedPattern,
                     std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
                         &reducedToFullInterfaceViMap,
                     std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
                         &fullToReducedInterfaceViMap,
                     std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
                         &fullPatternOppositeInterfaceViMap);
    static void visualizeResults(
        const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
        const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
        const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios,
        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
            &reducedToFullInterfaceViMap);
    static void registerResultsForDrawing(
        const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
        const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
            &reducedToFullInterfaceViMap);

    static double computeRawTranslationalError(
        const std::vector<Vector6d> &fullPatternDisplacements,
        const std::vector<Vector6d> &reducedPatternDisplacements,
        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
            &reducedToFullInterfaceViMap);

    static double computeDisplacementError(
        const std::vector<Vector6d> &fullPatternDisplacements,
        const std::vector<Vector6d> &reducedPatternDisplacements,
        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
            &reducedToFullInterfaceViMap,
        const double &normalizationFactor);

    static double computeRawRotationalError(
        const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
        const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
            &reducedToFullInterfaceViMap);

    static double computeRotationalError(const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
        const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
            &reducedToFullInterfaceViMap,
        const double &normalizationFactor);
    static double 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 &scenarioWeight,
        const ReducedModelOptimization::Settings::ObjectiveWeights &objectiveWeights);
    static void constructAxialSimulationScenario(
        const double &forceMagnitude,
        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
            &oppositeInterfaceViPairs,
        SimulationJob &job);

    static void constructShearSimulationScenario(
        const double &forceMagnitude,
        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
            &oppositeInterfaceViPairs,
        SimulationJob &job);

    static void constructBendingSimulationScenario(
        const double &forceMagnitude,
        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
            &oppositeInterfaceViPairs,
        SimulationJob &job);

    static void constructDomeSimulationScenario(
        const double &forceMagnitude,
        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
            &oppositeInterfaceViPairs,
        SimulationJob &job);

    static void constructSaddleSimulationScenario(
        const double &forceMagnitude,
        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
            &oppositeInterfaceViPairs,
        SimulationJob &job);
    static void constructSSimulationScenario(
        const double &forceMagnitude,
        const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
            &oppositeInterfaceViPairs,
        SimulationJob &job);
    static std::function<void(const std::vector<double> &x,
                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
        function_updateReducedPattern;
    static std::function<void(const double &newE,
                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
        function_updateReducedPattern_material_E;
    static std::function<void(const double &newA,
                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
        function_updateReducedPattern_material_A;
    static std::function<void(const double &newI,
                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
        function_updateReducedPattern_material_I;
    static std::function<void(const double &newI2,
                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
        function_updateReducedPattern_material_I2;
    static std::function<void(const double &newI3,
                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
        function_updateReducedPattern_material_I3;
    static std::function<void(const double &newJ,
                              std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
        function_updateReducedPattern_material_J;
    static double objective(const std::vector<double> &x);
    void initializeUpdateReducedPatternFunctions();
    //    static double objective(const double &xValue);

    ReducedModelOptimization::Results optimize(
        ConstPatternGeometry &fullPattern,
        const ReducedModelOptimization::Settings &optimizationSettings);

    void optimize(ConstPatternGeometry &fullPattern,
                  ReducedModelOptimization::Settings &optimizationSettings,
                  ReducedModelOptimization::Results &optimizationResults);
    static double objective(const std::vector<double> &x,
                            ReducedModelOptimizer::OptimizationState &optimizationState);

private:
    void optimize(
        ReducedModelOptimization::Settings &optimizationSettings,
        ReducedModelOptimization::Results &results,
        const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios
        = std::vector<ReducedModelOptimization::BaseSimulationScenario>(
            {ReducedModelOptimization::Axial,
             ReducedModelOptimization::Shear,
             ReducedModelOptimization::Bending,
             ReducedModelOptimization::Dome,
             ReducedModelOptimization::Saddle,
             ReducedModelOptimization::S}));

    void initializePatterns(PatternGeometry &fullPattern,
                            ReducedModel &reducedPattern,
                            const std::array<ReducedModelOptimization::xRange,
                                             ReducedModelOptimization::NumberOfOptimizationVariables>
                                &optimizationParameters);
    static void computeDesiredReducedModelDisplacements(
        const SimulationResults &fullModelResults,
        const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
        Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
    void runOptimization(const ReducedModelOptimization::Settings &settings,
                         ReducedModelOptimization::Results &results);
    void computeMaps(PatternGeometry &fullModel, ReducedModel &reducedModel);
    void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
    void initializeOptimizationParameters(
        const std::shared_ptr<SimulationMesh> &mesh,
        const std::array<ReducedModelOptimization::xRange,
                         ReducedModelOptimization::NumberOfOptimizationVariables>
            &optimizationParamters);

    DRMSimulationModel simulator;
    void computeObjectiveValueNormalizationFactors(
        const ReducedModelOptimization::Settings &optimizationSettings);

    void getResults(const FunctionEvaluation &optimalObjective,
                    const ReducedModelOptimization::Settings &settings,
                    ReducedModelOptimization::Results &results);
    double computeFullPatternMaxSimulationForce(
        const ReducedModelOptimization::BaseSimulationScenario &scenario) const;

#ifdef DLIB_DEFINED
    static double objective(const dlib::matrix<double, 0, 1> &x);
#endif
    std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
    computeFullPatternMaxSimulationForces(
        const std::vector<ReducedModelOptimization::BaseSimulationScenario>
            &desiredBaseSimulationScenario) const;
    std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
        const std::shared_ptr<SimulationMesh> &pMesh,
        const std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
            &baseScenarioMaxForceMagnitudes) const;
    std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
    getFullPatternMaxSimulationForces(
        const std::vector<ReducedModelOptimization::BaseSimulationScenario>
            &desiredBaseSimulationScenarioIndices,
        const std::filesystem::path &intermediateResultsDirectoryPath,
        const bool &recomputeForceMagnitudes);
    std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
    getFullPatternMaxSimulationForces();
    void computeScenarioWeights(
        const std::vector<ReducedModelOptimization::BaseSimulationScenario> &baseSimulationScenarios,
        const ReducedModelOptimization::Settings &optimizationSettings);
};
inline std::function<void(const double &newE,
                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
    ReducedModelOptimizer::function_updateReducedPattern_material_E;
inline std::function<void(const double &newA,
                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
    ReducedModelOptimizer::function_updateReducedPattern_material_A;
inline std::function<void(const double &newI,
                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
    ReducedModelOptimizer::function_updateReducedPattern_material_I;
inline std::function<void(const double &newI2,
                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
    ReducedModelOptimizer::function_updateReducedPattern_material_I2;
inline std::function<void(const double &newI3,
                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
    ReducedModelOptimizer::function_updateReducedPattern_material_I3;
inline std::function<void(const double &newJ,
                          std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
    ReducedModelOptimizer::function_updateReducedPattern_material_J;
inline std::function<void(const std::vector<double> &x, std::shared_ptr<SimulationMesh> &m)>
    ReducedModelOptimizer::function_updateReducedPattern;

extern ReducedModelOptimizer::OptimizationState global;

#endif // REDUCEDMODELOPTIMIZER_HPP