#ifndef REDUCEDMODELOPTIMIZER_HPP
#define REDUCEDMODELOPTIMIZER_HPP

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

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

class ReducedModelOptimizer
{
    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::vector<std::function<
        void(const double &,
             const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
             SimulationJob &)>>
        constructBaseScenarioFunctions;
    std::vector<bool> scenarioIsSymmetrical;
    int fullPatternNumberOfEdges;

public:
    constexpr static std::array<int, 5> simulationScenariosResolution = {11, 11, 20, 20, 20};
    //    constexpr static std::array<int, 5> simulationScenariosResolution = {3, 3, 3, 3, 3};
    inline static int totalNumberOfSimulationScenarios
        = std::accumulate(simulationScenariosResolution.begin(),
                          simulationScenariosResolution.end(),
                          0);
    inline static int fanSize{6};
    inline static double initialHexagonSize{0.3};
    inline static VectorType patternPlaneNormal{0, 0, 1};
    void optimize(
        const ReducedPatternOptimization::Settings &xRanges,
        ReducedPatternOptimization::Results &results,
        const std::vector<ReducedPatternOptimization::BaseSimulationScenario> &simulationScenarios
        = std::vector<ReducedPatternOptimization::BaseSimulationScenario>(
            {ReducedPatternOptimization::Axial,
             ReducedPatternOptimization::Shear,
             ReducedPatternOptimization::Bending,
             ReducedPatternOptimization::Dome,
             ReducedPatternOptimization::Saddle}));
    double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;

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

    SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);

    void initializePatterns(
        PatternGeometry &fullPattern,
        PatternGeometry &reducedPatterm,
        const std::vector<ReducedPatternOptimization::xRange> &optimizationParameters);

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

    static double objective(double x2,
                            double A,
                            double J,
                            double I2,
                            double I3,
                            double x3,
                            double innerHexagonRotationAngle);
    static double objective(double b, double r, double E);

    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);
    static void 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>>
                                &fullPatternOppositeInterfaceViMap);
    static void visualizeResults(
        const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
        const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
        const std::vector<ReducedPatternOptimization::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 objective(double E,
                            double A,
                            double innerHexagonSize,
                            double innerHexagonRotationAngle);
    static double computeRawRotationalError(
        const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
        const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
        const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
            &reducedToFullInterfaceViMap);

    static double 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);
    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);
    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 std::function<void(const dlib::matrix<double, 0, 1> &x,
                              std::shared_ptr<SimulationMesh> &m)>
        function_updateReducedPattern;
    static std::function<void(const dlib::matrix<double, 0, 1> &x,
                              std::shared_ptr<SimulationMesh> &m)>
        function_updateReducedPattern_material;
    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                              const double &newE)>
        function_updateReducedPattern_material_E;
    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                              const double &newA)>
        function_updateReducedPattern_material_A;
    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                              const double &newI)>
        function_updateReducedPattern_material_I;
    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                              const double &newI2)>
        function_updateReducedPattern_material_I2;
    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                              const double &newI3)>
        function_updateReducedPattern_material_I3;
    static std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                              const double &newJ)>
        function_updateReducedPattern_material_J;
    static std::function<void(const dlib::matrix<double, 0, 1> &x,
                              std::shared_ptr<SimulationMesh> &m)>
        function_updateReducedPattern_geometry;

private:
    static void computeDesiredReducedModelDisplacements(
        const SimulationResults &fullModelResults,
        const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
        Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
    static void runOptimization(const ReducedPatternOptimization::Settings &settings,
                                ReducedPatternOptimization::Results &results);
    std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
        const std::shared_ptr<SimulationMesh> &pMesh,
        const std::vector<std::pair<ReducedPatternOptimization::BaseSimulationScenario, double>>
            &maxForceMagnitudes);
    void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern);
    void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
    static void initializeOptimizationParameters(
        const std::shared_ptr<SimulationMesh> &mesh,
        const std::vector<ReducedPatternOptimization::xRange> &optimizationParamters);

    static double objective(long n, const double *x);
    DRMSimulationModel simulator;
    void computeObjectiveValueNormalizationFactors();
    static void getResults(const dlib::function_evaluation &optimizationResult_dlib,
                           const ReducedPatternOptimization::Settings &settings,
                           ReducedPatternOptimization::Results &results);
    double computeFullPatternMaxSimulationForce(
        const ReducedPatternOptimization::BaseSimulationScenario &scenario);
    std::vector<std::pair<ReducedPatternOptimization::BaseSimulationScenario, double>>
    computeFullPatternMaxSimulationForces(
        const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
            &desiredBaseSimulationScenarioIndices);

    std::vector<std::pair<ReducedPatternOptimization::BaseSimulationScenario, double>>
    getFullPatternMaxSimulationForces(
        const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
            &desiredBaseSimulationScenarioIndices);
    static double objective(const dlib::matrix<double, 0, 1> &x);
    static void initializeUpdateReducedPatternFunctions();
};
inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
    ReducedModelOptimizer::function_updateReducedPattern;
inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
    ReducedModelOptimizer::function_updateReducedPattern_material;
inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                          const double &newE)>
    ReducedModelOptimizer::function_updateReducedPattern_material_E;
inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                          const double &newA)>
    ReducedModelOptimizer::function_updateReducedPattern_material_A;
inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                          const double &newI)>
    ReducedModelOptimizer::function_updateReducedPattern_material_I;
inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                          const double &newI2)>
    ReducedModelOptimizer::function_updateReducedPattern_material_I2;
inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                          const double &newI3)>
    ReducedModelOptimizer::function_updateReducedPattern_material_I3;
inline std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
                          const double &newJ)>
    ReducedModelOptimizer::function_updateReducedPattern_material_J;
inline std::function<void(const dlib::matrix<double, 0, 1> &x, std::shared_ptr<SimulationMesh> &m)>
    ReducedModelOptimizer::function_updateReducedPattern_geometry;

#endif // REDUCEDMODELOPTIMIZER_HPP