//
// Created by Anna  Proost on 11/12/2023.
//

#include <fstream>
#include <string>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>

#include "../inputs/headers.hpp"
#include "../methods/headers.hpp"
#include "../outputs/Solution.hpp"
#include "Solver.hpp"


double Solver::EvaluateFunction(const std::string& eq, double x) {

    MonoFunction Eq(eq, "", "");
    return Eq(x);
}


double Solver::EvaluateMultiFunction(const std::string& eq, Eigen::VectorXd x, const std::vector<std::string>& var) {
    MultiFunction Eq(eq, std::vector<std::string>(var.size(), ""), var);
    return Eq(x);
}

Solution Solver::NewtonMethod(const std::string& eq, const std::string& der, double start_point, double stop_tol, int max_iter) {

    MonoFunction Eq(eq, "", der);

    Newton NewtonMethod;
    return NewtonMethod(Eq, start_point, stop_tol, max_iter) ;

}

Solution Solver::NewtonMethodA(const std::string& eq, const std::string& der, double start_point, double stop_tol, int max_iter) {

    MonoFunction Eq(eq, "", der);

    Newton NewtonMethod;
    Aitken AitkenMethod(&NewtonMethod);
    return AitkenMethod(Eq, start_point, stop_tol, max_iter);

}

Solution Solver::SystemNewtonMethod(const std::vector<std::string>& eqs, const std::vector<std::vector<std::string>>& der, const std::vector<std::string>& var, const std::vector<double>& start_points, double stop_tol, int max_iter) {

    std::vector<MultiFunction> functions;

    functions.reserve((int)eqs.size());
    for (int i = 0; i < (int)var.size(); i++) { functions.push_back( MultiFunction(eqs[i], der[i], var) ); }

    Eigen::VectorXd startVector = Eigen::VectorXd::Map(start_points.data(), (int)start_points.size());

    System system(functions);

    Newton NewtonMethod;
    return NewtonMethod(system, startVector, stop_tol, max_iter) ;

}

Solution Solver::FixedPointMethod(const std::string& eq, const std::string& fp, double start_point, double stop_tol, int max_iter) {

    MonoFunction Eq(eq, fp, "");

    FixedPoint FixedPointMethod;
    return FixedPointMethod(Eq, start_point, stop_tol, max_iter) ;

}

Solution Solver::FixedPointMethodA(const std::string& eq, const std::string& fp, double start_point, double stop_tol, int max_iter) {

    MonoFunction Eq(eq, fp, "");

    FixedPoint FixedPointMethod;
    Aitken AitkenMethod(&FixedPointMethod);
    return AitkenMethod(Eq, start_point, stop_tol, max_iter);

}

Solution Solver::ChordMethod(const std::string& eq, std::pair<double, double> start_point, double stop_tol, int max_iter) {

    MonoFunction Eq(eq, "", "");

    Chord ChordMethod;
    return ChordMethod(Eq, start_point, stop_tol, max_iter) ;

}

Solution Solver::BisectionMethod(const std::string& eq, std::pair<double, double> start_point, double stop_tol, int max_iter) {

    MonoFunction Eq(eq, "", "");

    Bisection BisectionMethod;
    return BisectionMethod(Eq, start_point, stop_tol, max_iter) ;

}

void Solver::PointOutput(Solution root, const std::string& method, bool Aitken, const std::string& function){

    if (Aitken) std::cout << "\n\033[1mMethod was run with Aitken acceleration\033[0m\n";
    if (root.hasConverged()) std::cout << "The root \033[1m"<< std::fixed << std::setprecision(15) << root.getRoot() << "\033[0m of the function f(x) = " << function <<" was found by the " << method << " method in " <<  root.getNIterations() <<" iterations. The final function evaluation is " << EvaluateFunction(function, root.getRoot()) <<  "\n";
    else std::cout << "The " << method << " method does not converge for the initial point or within the max. number of iterations";

    std::cout << "\nThis is an overview of the results:";
    root.printIterations();

}

void Solver::IntervalOutput(Solution root, const std::string& method, std::string function){

    if (root.hasConverged()) std::cout << "The root \033[1m"<< std::fixed << std::setprecision(15) << root.getRoot() << "\033[0m of the function f(x) = " << function << " in interval [" << root.getIntervalIterations()[root.getNIterations()].first << ", " << root.getIntervalIterations()[root.getNIterations()].second << "] was found by the " << method << " method in " <<  root.getNIterations() <<" iterations. The final function evaluation is " << EvaluateFunction(function, root.getRoot()) <<  "\n";
    else std::cout << "The " << method << " method does not converge for the initial point or within the max. number of iterations";

    std::cout << "\nThis is an overview of the results:";
    root.printIterations();

}

void Solver::SystemOutput(Solution root, std::vector<std::string> functions){

    std::string combinedFunctions = std::accumulate(
            functions.begin() + 1, functions.end(),
            functions[0],
            [](const std::string& acc, const std::string& str) {
                return acc + ", " + str;
            }
    );

    if (root.hasConverged()) {
        std::cout << "The root\033[1m (" << std::fixed << std::setprecision(15) << root.getVectorRoot()[0] << ", "
                  << root.getVectorRoot()[1] << ")\033[0m of the functions {" << combinedFunctions <<
                                                "} was found by the Newton method for Systems in "
                  << root.getNIterations() << " iterations. \n";
    }

    else std::cout << "The Newton method for Systems does not converge for the initial points or within the max. number of iterations";

    std::cout << "\nThis is an overview of the results:";
    root.printIterations();

}