#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/matchers/catch_matchers_range_equals.hpp>

#include "mini-odeint.hpp"

unsigned int Factorial(unsigned int number) {
  return number <= 1 ? number : Factorial(number - 1) * number;
}

TEST_CASE("Factorials are computed", "[factorial]") {
  REQUIRE(Factorial(1) == 1);
  REQUIRE(Factorial(2) == 2);
  REQUIRE(Factorial(3) == 6);
  REQUIRE(Factorial(10) == 3628800);
}

TEST_CASE("Mini-odeint works", "[mini-odeint]") {
  using namespace mini_odeint;

  SECTION("float") {
    std::vector<float> times;
    times.reserve(1001);
    for (int i = 0; i <= 1000; ++i) {
      times.push_back(i / 1000.0f);
    }

    SECTION("scalar") {
      std::vector<float> ys(times.size());

      explicitRungeKutta(std::span(ys), std::span<const float>(times), 1.0f,
                         1e-6f, [](auto y, auto t) { return -y; });

      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
                     return std::abs(a - std::exp(-b)) < 1e-6;
                   }));
    }

    SECTION("2 vector") {
      std::vector<std::array<float, 2>> ys(times.size());

      explicitRungeKutta(
          std::span(ys), std::span<const float>(times), std::array{1.0f, 2.0f},
          1e-6f, [](auto y, auto t) { return std::array{-y[0], -y[1]}; });

      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
                     const auto s = std::exp(-b);

                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s}) <
                            1e-6;
                   }));
    }

    SECTION("4 vector") {
      std::vector<std::array<float, 4>> ys(times.size());

      explicitRungeKutta(std::span(ys), std::span<const float>(times),
                         std::array{1.0f, 2.0f, 3.0f, 4.0f}, 1e-6f,
                         [](auto y, auto t) {
                           return std::array{-y[0], -y[1], -y[2], -y[3]};
                         });

      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
                     const auto s = std::exp(-b);

                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s,
                                                      a[2] - 3 * s,
                                                      a[3] - 4 * s}) < 1e-6;
                   }));
    }
  }

  SECTION("float") {
    std::vector<double> times;
    times.reserve(1001);
    for (int i = 0; i <= 1000; ++i) {
      times.push_back(i / 1000.0);
    }

    SECTION("scalar") {
      std::vector<double> ys(times.size());

      explicitRungeKutta(std::span(ys), std::span<const double>(times), 1.0,
                         1e-6, [](auto y, auto t) { return -y; });

      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
                     return std::abs(a - std::exp(-b)) < 1e-6;
                   }));
    }

    SECTION("2 vector") {
      std::vector<std::array<double, 2>> ys(times.size());

      explicitRungeKutta(
          std::span(ys), std::span<const double>(times), std::array{1.0, 2.0},
          1e-6, [](auto y, auto t) { return std::array{-y[0], -y[1]}; });

      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
                     const auto s = std::exp(-b);

                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s}) <
                            1e-6;
                   }));
    }

    SECTION("4 vector") {
      std::vector<std::array<double, 4>> ys(times.size());

      explicitRungeKutta(std::span(ys), std::span<const double>(times),
                         std::array{1.0, 2.0, 3.0, 4.0}, 1e-6,
                         [](auto y, auto t) {
                           return std::array{-y[0], -y[1], -y[2], -y[3]};
                         });

      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
                     const auto s = std::exp(-b);

                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s,
                                                      a[2] - 3 * s,
                                                      a[3] - 4 * s}) < 1e-6;
                   }));
    }
  }
}