#include <catch2/catch.hpp>
#include <iostream>

#include "../utils.hpp"
#include "../TOMS748.hpp"
#include "utils_test.hpp"

#define print(x) PRINT_VAR(x);
#define printvec(x) PRINT_VEC(x);
#define printstr(x) PRINT_STR(x);

template<typename T> T f(T E) { return E - 0.5*sin(E) - 0.3; }
template<typename T> T f2(T x) { return exp(x)-T(2); }
template<typename T> T f3(T x) { return (x*x)-T(1); }

using namespace TOMS748;
using namespace TOMS748::internal;

TEST_CASE( "sureAbs", "[TOMS748]" )
{
    REQUIRE(TOMS748::internal::sureAbs(-1) == 1);
    REQUIRE(TOMS748::internal::sureAbs(-1.) == 1.);
}

TEST_CASE( "checkTwoValuesIdentical", "[TOMS748]" )
{
    std::cout.precision(16);

    SECTION( "checkTwoValuesIdentical" ) {
        REQUIRE(checkTwoValuesIdentical(1,2,3,4) == false);
		REQUIRE(checkTwoValuesIdentical(1,2,3,1) == true);
		REQUIRE(checkTwoValuesIdentical(1,1,3,4) == true);
		REQUIRE(checkTwoValuesIdentical(1,2,2,4) == true);
		REQUIRE(checkTwoValuesIdentical(1,2,3,3) == true);
		REQUIRE(checkTwoValuesIdentical(1,2,1,4) == true);
		REQUIRE(checkTwoValuesIdentical(1,2,3,2) == true);
    }
}

TEST_CASE( "swap_ab_fafb_if_not_ascending", "[TOMS748]" )
{
    std::cout.precision(16);

    SECTION( "ascending (do nothing)" ) {
        double  a = 1., b = 2.,
                fa = f(a), fb = f(b);

        swap_ab_fafb_if_not_ascending(a,b,fa,fb);

        REQUIRE(a == 1.);
        REQUIRE(b == 2.);
        REQUIRE(fa == f(1.));
        REQUIRE(fb == f(2.));
    }

    SECTION( "descending (swap)" ) {
        double  a = 2., b = 1.,
                fa = f(a), fb = f(b);

        swap_ab_fafb_if_not_ascending(a,b,fa,fb);

        REQUIRE(a == 1.);
        REQUIRE(b == 2.);
        REQUIRE(fa == f(1.));
        REQUIRE(fb == f(2.));
    }
}

TEST_CASE( "Bracket", "[TOMS748]" )
{
    std::cout.precision(16);

    double a = 0, b = 1, c = 0.3, d = 1.2;
    double fa = f(a), fb = f(b), fc = f(c), fd = f(d);
    double tol = 1e-12;
    bool ok;

    SECTION( "checkTwoValuesIdentical" ) {
        REQUIRE(check_almost_equal(fbracket1(a, b, fa, fb), 0.5792645075960517, tol));
        REQUIRE(check_almost_equal(fbracket2(a, b, d, fa, fb, fd), 0.1619293662546865, tol));
    }

    SECTION( "bracket normal operation [c b]" ) {
        std::tie(a, b, d, fa, fb, fd) = bracket(a, b, c, fa, fb, fc);
        REQUIRE(check_almost_equal(a, 0.3, tol));
        REQUIRE(check_almost_equal(b, 1., tol));
        REQUIRE(check_almost_equal(d, 0., tol));
        REQUIRE(check_almost_equal(fa, f(0.3), tol));
        REQUIRE(check_almost_equal(fb, f(1.), tol));
        REQUIRE(check_almost_equal(fd, f(0.), tol));
    }

    SECTION( "bracket normal operation [a c]" ) {
        std::tie(a, b, d, fa, fb, fd) = bracket(a, b, 0.7, fa, fb, f(0.7));
        REQUIRE(check_almost_equal(a, 0., tol));
        REQUIRE(check_almost_equal(b, 0.7, tol));
        REQUIRE(check_almost_equal(d, 1., tol));
        REQUIRE(check_almost_equal(fa, f(0.), tol));
        REQUIRE(check_almost_equal(fb, f(0.7), tol));
        REQUIRE(check_almost_equal(fd, f(1.), tol));
    }

    SECTION( "bracket and check normal operation [c b]" ) {
        std::tie(a, b, d, fa, fb, fd, ok) = bracketAndCheckConvergence(a, b, c, fa, fb, fc, tol);
        REQUIRE(check_almost_equal(a, 0.3, tol));
        REQUIRE(check_almost_equal(b, 1., tol));
        REQUIRE(check_almost_equal(d, 0., tol));
        REQUIRE(check_almost_equal(fa, f(0.3), tol));
        REQUIRE(check_almost_equal(fb, f(1.), tol));
        REQUIRE(check_almost_equal(fd, f(0.), tol));
        REQUIRE_FALSE(ok);
    }

    SECTION( "bracket and check normal operation [a c]" ) {
        std::tie(a, b, d, fa, fb, fd, ok) = bracketAndCheckConvergence(a, b, 0.7, fa, fb, f(0.7), tol);
        REQUIRE(check_almost_equal(a, 0., tol));
        REQUIRE(check_almost_equal(b, 0.7, tol));
        REQUIRE(check_almost_equal(d, 1., tol));
        REQUIRE(check_almost_equal(fa, f(0.), tol));
        REQUIRE(check_almost_equal(fb, f(0.7), tol));
        REQUIRE(check_almost_equal(fd, f(1.), tol));
        REQUIRE_FALSE(ok);
    }

    SECTION( "bracket and check normal operation c is root" ) {
        std::tie(a, b, d, fa, fb, fd, ok) = bracketAndCheckConvergence(a, b, 0.569682256443945, fa, fb, f(0.569682256443945), tol);
        REQUIRE(check_almost_equal(a, 0., tol));
        REQUIRE(check_almost_equal(b, 0.569682256443945, tol));
        REQUIRE(check_almost_equal(d, 1., tol));
        REQUIRE(check_almost_equal(fa, f(0.), tol));
        REQUIRE(check_almost_equal(fb, f(0.569682256443945), tol));
        REQUIRE(check_almost_equal(fd, f(1.), tol));
        REQUIRE(ok);
    }
}

TEST_CASE( "TOMS748_solve1", "[TOMS748]")
{
    double a = 0, b = 1, tol = 1e-12;

    SECTION( "TOMS748_solve1 good" ) {
        auto [x, fx, ok] = TOMS748_solve1(f<double>, a, b, tol);
        CHECK(check_almost_equal(fx, 0., tol));
        CHECK(ok == true);
    }

    SECTION( "TOMS748_solve1 wrong initial bracket" ) {
        auto [x, fx, ok] = TOMS748_solve1(f<double>, 1., 2., tol);// still converges !
        CHECK(check_almost_equal(fx, 0., tol));
        CHECK(ok == true);
    }
}

TEST_CASE( "findEnclosingBracket", "[TOMS748]")
{
    double tol = 1e-12, factor = 2;

    SECTION( "findEnclosingBracket start at 0" ) {
        double x = 0;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
        CHECK(a == x);
        CHECK(a == b);
        CHECK_FALSE(fa*fb <= 0.);
        CHECK_FALSE(ok);
    }

    SECTION( "findEnclosingBracket sol below" ) {
        double x = 2.;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
        CHECK(check_almost_equal(b, x, tol));
        CHECK(a < b);
        CHECK(fa*fb <= 0.);
        CHECK(ok);
    }

    SECTION( "findEnclosingBracket sol above" ) {
        double x = .1;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
        CHECK(check_almost_equal(a, x, tol));
        CHECK(a < b);
        CHECK(fa*fb <= 0.);
        CHECK(ok);
    }

    SECTION( "findEnclosingBracket initial bracket good" ) {
        double x = .5;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
        CHECK(check_almost_equal(a, x, tol));
        CHECK(a < b);
        CHECK(fa*fb <= 0.);
        CHECK(ok);
    }

    SECTION( "findEnclosingBracket sol above far far away" ) {
        double x = 1e-10;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
        CHECK(check_almost_equal(a, x, tol));
        CHECK(a < b);
        CHECK(fa*fb <= 0.);
        CHECK(ok);
    }

    SECTION( "findEnclosingBracket sol on the other side of 0 (does not converge)" ) {
        double x = -1.;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
        CHECK(a < b);
        CHECK_FALSE(fa*fb <= 0.);
        CHECK_FALSE(ok);
    }

    SECTION( "findEnclosingBracket sol below (negative)" ) {
        double x = -.1;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f3<double>, x, factor);
        CHECK(check_almost_equal(b, x, tol));
        CHECK(a < b);
        CHECK(fa*fb <= 0.);
        CHECK(ok);
    }

    SECTION( "findEnclosingBracket sol above (negative)" ) {
        double x = -2.;
        auto [a, b, fa, fb, ok] = findEnclosingBracket(f3<double>, x, factor);
        CHECK(check_almost_equal(a, x, tol));
        CHECK(a < b);
        CHECK(fa*fb <= 0.);
        CHECK(ok);
    }
}

TEST_CASE( "bracket_and_solve", "[TOMS748]" )
{
    double tol = 1e-12;
    SECTION( "initial guess close to solution" ) {
        auto [x, fx, ok] = bracket_and_solve(f<double>, 0.5, tol);
        CHECK(check_almost_equal(fx, 0., tol));
        CHECK(ok == true);
    }

    SECTION( "initial guess far from solution" ) {
        auto [x, fx, ok] = bracket_and_solve(f<double>, 100., tol);
        CHECK(check_almost_equal(fx, 0., tol));
        CHECK(ok == true);
    }

    SECTION( "initial guess on other side of 0 compared to solution" ) {
        auto [x, fx, ok] = bracket_and_solve(f<double>, -1., tol);
        CHECK_FALSE(check_almost_equal(fx, 0., tol));
        CHECK_FALSE(ok);
    }
}