#include "reducedmodelevaluator.hpp"
#include "hexagonremesher.hpp"
#include "reducedmodel.hpp"
#include "trianglepatterngeometry.hpp"
#include <filesystem>

using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;

ReducedModelEvaluator::ReducedModelEvaluator()
{

}
ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
    ReducedModelOptimization::Results &optimizationResult)
{
    const std::filesystem::path tileInto_triMesh_filePath
        = "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/"
          "instantMeshes_plane_34.ply";
    const std::filesystem::path scenariosDirectoryPath
        = "/home/iason/Coding/Projects/Approximating shapes with flat "
          "patterns/ReducedModelEvaluator/Scenarios";
    //    const std::filesystem::path reducedPatternFilePath
    //        = "/home/iason/Coding/Projects/Approximating shapes with flat "
    //          "patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply";
    const std::filesystem::path fullPatternTessellatedResultsDirectoryPath
        = "/home/iason/Coding/Projects/Approximating shapes with flat "
          "patterns/ReducedModelEvaluator/TessellatedResults";

    return evaluateReducedModel(optimizationResult,
                                tileInto_triMesh_filePath,
                                scenariosDirectoryPath,
                                //                                reducedPatternFilePath,
                                fullPatternTessellatedResultsDirectoryPath);
}
//void ReducedModelEvaluator::printResults(const Results &evaluationResults,
//                                         const std::string &resultsLabel,
//                                         const std::filesystem::path &resultsOutputPath)
//{
//    const bool outputPathIsDirectory = !resultsOutputPath.empty()
//                                       && !resultsOutputPath.has_extension();
//    const bool outputPathIsFile = !resultsOutputPath.empty() && resultsOutputPath.has_extension();
//    assert(outputPathIsDirectory && outputPathIsFile);
//    if (outputPathIsDirectory) {
//        std::filesystem::create_directories(resultsOutputPath);
//    }

//#else
//    std::filesystem::path csvOutputFilePath;
//    bool shouldOverwrite = false;
//    if (outputPathIsDirectory) {
//        csvOutputFilePath = std::filesystem::path(resultsOutputPath)
//                                .append("distances_" + resultsLabel + ".csv")
//                                .string();
//        shouldOverwrite = true;
//    } else if (outputPathIsFile) {
//        csvOutputFilePath = resultsOutputPath;
//    }

//    csvFile csvOutput(csvOutputFilePath, shouldOverwrite);
//    printResults(evaluationResults, resultsLabel, csvOutput);
//}

void ReducedModelEvaluator::printResults(const Results &evaluationResults,
                                         const std::string &resultsLabel)
{
    csvFile csvOutputToCout({}, true);
    Settings exportSettings;
    exportSettings.exportingDirection = Vertical;
    exportSettings.shouldWriteHeader = false;
    exportSettings.resultsLabel = resultsLabel;
    printResults(evaluationResults, exportSettings, csvOutputToCout);
}

void ReducedModelEvaluator::printResults(const Results &evaluationResults,
                                         const Settings &settings,
                                         csvFile &csvOutput)
{
    if (settings.shouldWriteHeader) {
        printHeader(settings, csvOutput);
        //    csvOutput << "Average error";
        //    csvOutput<<"Cumulative error";
        csvOutput << endrow;
    }
    if (!settings.resultsLabel.empty()) {
        csvOutput << settings.resultsLabel;
    }
    if (settings.exportingDirection == Vertical) {
        printResultsVertically(evaluationResults, csvOutput);
    } else {
        printResultsHorizontally(evaluationResults, csvOutput);
    }
}

void ReducedModelEvaluator::printHeader(const Settings &settings, csvFile &csvOutput)
{
    if (settings.exportingDirection == Horizontal) {
        csvOutput << csvExportingDataStrings[settings.exportingData];
        csvOutput.endrow();
        //        csvOutput << "Job label";
        for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
             jobIndex++) {
            const std::string &jobLabel = ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
            csvOutput << jobLabel;
        }
    } else {
        std::cout << "Vertical header not defined" << std::endl;
        assert(false);
        std::terminate();
    }
}

void ReducedModelEvaluator::printResultsHorizontally(const Results &evaluationResults,
                                                     csvFile &csvOutput)
{
    //print header
    //print raw error
    constexpr bool shouldPrintRawError = false;
    if (shouldPrintRawError) {
        //        csvOutput << "drm2Reduced";
        double sumOfFull2Reduced = 0;
        int numOfNonEvaluatedScenarios = 0;
        for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
             jobIndex++) {
            const double &distance_fullDrmToReduced = evaluationResults
                                                          .distances_drm2reduced[jobIndex];
            if (distance_fullDrmToReduced == -1) {
                csvOutput << "notEvaluated";
                numOfNonEvaluatedScenarios++;
            } else {
                csvOutput << distance_fullDrmToReduced;
                sumOfFull2Reduced += distance_fullDrmToReduced;
            }
        }
        //    const int numOfEvaluatedScenarios = ReducedModelEvaluator::scenariosTestSetLabels.size()
        //                                        - numOfNonEvaluatedScenarios;
        //    const double averageDistance_full2Reduced = sumOfFull2Reduced / numOfEvaluatedScenarios;
        //    csvOutput << averageDistance_full2Reduced;
        //        csvOutput << endrow;
    }
    //print normalized error
    //    csvOutput << "norm_drm2Reduced";
    double sumOfNormalizedFull2Reduced = 0;
    for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
         jobIndex++) {
        const double &distance_normalizedFullDrmToReduced
            = evaluationResults.distances_normalizedDrm2reduced[jobIndex];
        if (distance_normalizedFullDrmToReduced == -1) {
            csvOutput << "notEvaluated";
        } else {
            csvOutput << distance_normalizedFullDrmToReduced;
            sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
        }
    }
    //    const double averageDistance_normalizedFull2Reduced = sumOfNormalizedFull2Reduced
    //                                                          / numOfEvaluatedScenarios;
    //    csvOutput << averageDistance_normalizedFull2Reduced;
    //    csvOutput << endrow;
}

void ReducedModelEvaluator::printResultsVertically(const Results &evaluationResults,
                                                   csvFile &csvOutput)
{
#ifdef POLYSCOPE_DEFINED
    csvOutput << "drm2Reduced"
              << "norm_drm2Reduced";
#else
    csvOutput << "Job Label"
              << "drm2Reduced"
              << "norm_drm2Reduced";

#endif
    csvOutput << endrow;
    double sumOfFull2Reduced = 0;
    double sumOfNormalizedFull2Reduced = 0;
    int numOfNonEvaluatedScenarios = 0;
    for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
         jobIndex++) {
        const double &distance_fullDrmToReduced = evaluationResults.distances_drm2reduced[jobIndex];
        const double &distance_normalizedFullDrmToReduced
            = evaluationResults.distances_normalizedDrm2reduced[jobIndex];
#ifndef POLYSCOPE_DEFINED
        const std::string &jobLabel = ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
        csvOutput << jobLabel;
#endif
        if (distance_fullDrmToReduced == -1) {
            csvOutput << "notEvaluated"
                      << "notEvaluated";
            numOfNonEvaluatedScenarios++;
        } else {
            csvOutput << distance_fullDrmToReduced << distance_normalizedFullDrmToReduced;
            sumOfFull2Reduced += distance_fullDrmToReduced;
            sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
        }
        csvOutput << endrow;
        //        sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
    }
    const int numOfEvaluatedScenarios = ReducedModelEvaluator::scenariosTestSetLabels.size()
                                        - numOfNonEvaluatedScenarios;
    const double averageDistance_full2Reduced = sumOfFull2Reduced / numOfEvaluatedScenarios;
    const double averageDistance_normalizedFull2Reduced = sumOfNormalizedFull2Reduced
                                                          / numOfEvaluatedScenarios;
#ifndef POLYSCOPE_DEFINED
    csvOutput << "Average error";
#endif
    csvOutput << averageDistance_full2Reduced << averageDistance_normalizedFull2Reduced;
    csvOutput << endrow;
#ifndef POLYSCOPE_DEFINED
    csvOutput << "Cumulative error";
#endif
    csvOutput << sumOfFull2Reduced << sumOfNormalizedFull2Reduced;
    csvOutput << endrow;
    csvOutput << endrow;
}

ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
    ReducedModelOptimization::Results &optimizationResult,
    const std::filesystem::path &tileInto_triMesh_filePath,
    const std::filesystem::path &scenariosDirectoryPath,
    //    const std::filesystem::path &reducedPatternFilePath,
    const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath)
{
    //    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = std::make_shared<VCGPolyMesh>();
    //    std::filesystem::path tileIntoSurfaceFilePath{
    //        "/home/iason/Coding/Projects/Approximating shapes with flat "
    //        "patterns/ReducedModelOptimization/TestSet/TileIntoSurface/plane_34Polygons.ply"};
    //    assert(std::filesystem::exists(tileIntoSurfaceFilePath));
    //    pTileIntoSurface->load(tileIntoSurfaceFilePath);

    //Set required file paths
    //    const std::filesystem::path drmSettingsFilePath
    //        = "/home/iason/Coding/Projects/Approximating shapes with flat "
    //          "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/"
    //          "defaultDRMSimulationSettings.json";
    DRMSimulationModel::Settings drmSimulationSettings;
    drmSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
    //    drmSimulationSettings.load(drmSettingsFilePath);
    drmSimulationSettings.linearGuessForceScaleFactor = 0.5;
    drmSimulationSettings.debugModeStep = 10000;
    drmSimulationSettings.beVerbose = true;
    drmSimulationSettings.maxDRMIterations = 1e6;
    //    drmSimulationSettings.translationalKineticEnergyThreshold = 1e-15;
    //    drmSimulationSettings.viscousDampingFactor = 1e-2;
    constexpr bool shouldRerunFullPatternSimulation = false;

    //Load surface
    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
        VCGTriMesh tileInto_triMesh;
        const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filePath);
        assert(surfaceLoadSuccessfull);
        return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
    }();
    const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResult.baseTriangle.cP(0),
                                                             (optimizationResult.baseTriangle.cP(1)
                                                              + optimizationResult.baseTriangle.cP(
                                                                  2))
                                                                 / 2);
    pTileIntoSurface->moveToCenter();
    const double scaleFactor = optimizedBaseTriangleHeight
                               / pTileIntoSurface->getAverageFaceRadius();
    vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
    //    tileIntoSurface.registerForDrawing();
    //    polyscope::show();

    //Tile full pattern into surface
    std::vector<PatternGeometry> fullPatterns(1);
    fullPatterns[0].copy(optimizationResult.baseTriangleFullPattern);
    ////        Base triangle pattern might contain dangling vertices.Remove those
    fullPatterns[0].interfaceNodeIndex = 3;
    fullPatterns[0].deleteDanglingVertices();
    std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN(), 0);
    std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
    std::vector<size_t> tileIntoEdgeToTiledFullVi;
    std::shared_ptr<PatternGeometry> pTiledFullPattern
        = PatternGeometry::tilePattern(fullPatterns,
                                       {},
                                       *pTileIntoSurface,
                                       perSurfaceFacePatternIndices,
                                       tileIntoEdgeToTiledFullVi,
                                       perPatternIndexTiledFullPatternEdgeIndices);
    pTiledFullPattern->setLabel("Tiled_full_patterns");
    //    pTiledFullPattern->registerForDrawing();
    //Tile reduced pattern into surface
    //    PatternGeometry reducedPattern;
    ReducedModel reducedModel;
    //    reducedModel.registerForDrawing();
    //    polyscope::show();
    reducedModel.deleteDanglingVertices();
    //    reducedPattern.interfaceNodeIndex = 1;
    //    assert(reducedPattern.interfaceNodeIndex == 1);
    std::vector<PatternGeometry> reducedPatterns(1);
    reducedPatterns[0].copy(reducedModel);
    const auto reducedPatternBaseTriangle = reducedModel.computeBaseTriangle();
    ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon(
        optimizationResult, reducedPatternBaseTriangle, reducedPatterns[0]);

    std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
    std::vector<size_t> tileIntoEdgeToTiledReducedVi;
    std::shared_ptr<PatternGeometry> pTiledReducedPattern
        = PatternGeometry::tilePattern(reducedPatterns,
                                       {0},
                                       *pTileIntoSurface,
                                       perSurfaceFacePatternIndices,
                                       tileIntoEdgeToTiledReducedVi,
                                       perPatternIndexTiledReducedPatternEdgeIndices);
    pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
    //    pTiledReducedPattern->registerForDrawing();
    //    polyscope::show();

    std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
        fullToReducedViMap; //of only the common vertices
    std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
        reducedToFullViMap; //of only the common vertices
    for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
        fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
    }
    constructInverseMap(fullToReducedViMap, reducedToFullViMap);

    std::vector<size_t> tilledFullPatternInterfaceVi;
    tilledFullPatternInterfaceVi.clear();
    tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
    size_t viIndex = 0;
    for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
        tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
    }

    //Create simulation meshes
    ////Tessellated full pattern simulation mesh
    std::shared_ptr<SimulationMesh> pTiledFullPattern_simulationMesh;
    pTiledFullPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledFullPattern);
    //NOTE: Those should be derived from the optimization results instead of hardcoding them
    pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
    if (optimizationResult.fullPatternYoungsModulus == 0) {
        std::cerr << "Full pattern's young modulus not found." << std::endl;
        std::terminate();
    }
    pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
                                                      optimizationResult.fullPatternYoungsModulus);
    pTiledFullPattern_simulationMesh->reset();

    ////Tessellated reduced pattern simulation mesh
    std::shared_ptr<SimulationMesh> pTiledReducedPattern_simulationMesh;
    pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
    const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
        = perPatternIndexTiledReducedPatternEdgeIndices[0];
    ReducedModelOptimization::Results::applyOptimizationResults_elements(
        optimizationResult, pTiledReducedPattern_simulationMesh);
    pTiledReducedPattern_simulationMesh->reset();
    Results evaluationResults;
    evaluationResults.distances_drm2reduced.fill(-1);
    evaluationResults.distances_normalizedDrm2reduced.fill(-1);
    for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
        const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
        const std::filesystem::path tiledReducedPatternJobFilePath
            = std::filesystem::path(scenariosDirectoryPath)
                  .append(pTileIntoSurface->getLabel())
                  .append(jobLabel)
                  .append("ReducedJob")
                  .append(SimulationJob::jsonDefaultFileName);
        if (!std::filesystem::exists(tiledReducedPatternJobFilePath)) {
            std::cerr << "Scenario " << jobLabel
                      << " not found in:" << tiledReducedPatternJobFilePath << std::endl;
            continue;
        }
        //set jobs
        std::shared_ptr<SimulationJob> pJob_tiledReducedPattern;
        pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
        pJob_tiledReducedPattern->load(tiledReducedPatternJobFilePath, false);
        pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
        std::shared_ptr<SimulationJob> pJob_tiledFullPattern;
        pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
        pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
        pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
        //        pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
        //        pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
        //    polyscope::show();
        //Save reduced job

        const std::filesystem::path surfaceFolderPath
            = std::filesystem::path(fullPatternTessellatedResultsDirectoryPath)
                  .append(pTileIntoSurface->getLabel());
        const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
        const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath)
                                                                .append(scenarioLabel);
        const std::filesystem::path reducedJobDirectoryPath
            = std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
        std::filesystem::create_directories(reducedJobDirectoryPath);
        pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
        //Run scenario
        ////Full
        //        const std::string patternLabel = std::to_string(
        //                                             optimizationResult.baseTriangleFullPattern.EN())
        //                                         + "_"
        //                                         + optimizationResult.baseTriangleFullPattern.getLabel();
        const std::string &patternLabel = optimizationResult.baseTriangleFullPattern.getLabel();
        const auto fullResultsFolderPath
            = std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results");
        if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
            std::filesystem::remove_all(fullResultsFolderPath);
        }

        const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
                                                                   scenarioDirectoryPath)
                                                                   .append(patternLabel)
                                                                   .append("SimulationJob");
        SimulationResults simulationResults_tiledFullPattern_drm;
        if (std::filesystem::exists(fullResultsFolderPath)) {
            //Load full pattern results
            assert(std::filesystem::exists(fullPatternJobFolderPath));
            simulationResults_tiledFullPattern_drm.load(fullResultsFolderPath,
                                                        fullPatternJobFolderPath);
            //#ifdef POLYSCOPE_DEFINED
            //            std::array<double, 3> resultsColor({28.0, 99.0, 227.0});
            //            simulationResults_tiledFullPattern_drm.registerForDrawing(resultsColor);
            //            std::ifstream ifs("CameraSettings.json");
            //            nlohmann::json json;
            //            ifs >> json;
            //            polyscope::view::setCameraFromJson(json.dump(), false);
            //            //            polyscope::show();
            //            const std::string cameraJson = polyscope::view::getCameraJson();
            //            std::filesystem::path jsonFilePath("CameraSettings.json");
            //            std::ofstream jsonFile_cameraSettings(jsonFilePath.string());
            //            jsonFile_cameraSettings << cameraJson;
            //            jsonFile_cameraSettings.close();
            //            std::filesystem::create_directories("screenshots");
            //            const std::string screenshotOutputFilePath
            //                = (std::filesystem::current_path()
            //                       .append("screenshots")
            //                       .append(optimizationResult.label + "_" + pJob_tiledFullPattern->getLabel()))
            //                      .string()
            //                  + ".png";
            //            //            std::cout << "Saving image to:" << screenshotOutputFilePath << std::endl;
            //            polyscope::screenshot(screenshotOutputFilePath, false);
            //            simulationResults_tiledFullPattern_drm.unregister();
            //#endif

            simulationResults_tiledFullPattern_drm.converged = true;
        } else {
            std::cout << "Drm results not found in:" << fullResultsFolderPath << std::endl;
            //Full
            std::cout << "Executing:" << jobLabel << std::endl;
            DRMSimulationModel drmSimulationModel;
            simulationResults_tiledFullPattern_drm
                = drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
            simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
        }
        std::filesystem::create_directories(fullResultsFolderPath);
        const std::filesystem::path drmResultsOutputPath
            = std::filesystem::path(scenarioDirectoryPath).append(patternLabel);
        simulationResults_tiledFullPattern_drm.save(drmResultsOutputPath);
        if (!simulationResults_tiledFullPattern_drm.converged) {
            std::cerr << "Full pattern simulation failed." << std::endl;
            std::cerr << "Saved failed simulation to:" << drmResultsOutputPath << std::endl;
            continue;
        }

        LinearSimulationModel linearSimulationModel;
        SimulationResults simulationResults_tiledReducedPattern
            = linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);

        //        simulationResults_tiledReducedPattern.registerForDrawing();
        //        simulationResults_tiledFullPattern_drm.registerForDrawing();
        //        polyscope::show();
        //measure distance
        const double distance_fullDrmToReduced
            = simulationResults_tiledFullPattern_drm
                  .computeDistance(simulationResults_tiledReducedPattern, fullToReducedViMap);
        double distance_fullSumOfAllVerts = 0;
        for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
            distance_fullSumOfAllVerts += simulationResults_tiledFullPattern_drm
                                              .displacements[fullToReducedPair.first]
                                              .getTranslation()
                                              .norm();
        }
        const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced
                                                           / distance_fullSumOfAllVerts;
        evaluationResults.distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
        evaluationResults.distances_normalizedDrm2reduced[jobIndex]
            = distance_normalizedFullDrmToReduced;
    }

    return evaluationResults;
}