/// @file
/// @brief Task A
/// Debug: g++ -O0 -g -std=c++20 taskB.cpp -Ieigen -Imodules -o taskB && ./taskB
/// Release: g++ -DNDEBUG -O3 -std=c++20 taskB.cpp -Ieigen -Imodules -o taskB && ./taskB

// http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#title4
#include <Eigen/Dense> // MatrixXd, VectorXd

#include <cassert>
#include <chrono>
#include <cmath>
#include <iostream>
#include <random>

/// @brief Funtionality equivalent function 'poly_fit' in taskA.py
std::vector<double> poly_fit(std::vector<double> x_coords, std::vector<double> y_coords, size_t order) {

  assert(x_coords.size() == y_coords.size());

  using Eigen::MatrixXd;
  using Eigen::VectorXd;

  size_t m = y_coords.size();
  size_t n = order + 1;

  auto A = MatrixXd(m, n);
  auto b = VectorXd(m);

  for (size_t i = 0; i != m; ++i) {
    auto row = VectorXd(n);
    for (size_t j = 0; j != n; ++j) {
      row(j) = std::pow(x_coords[i], j);
    }
    A.row(i) = row;
    b(i) = y_coords[i];
  }

  using std::chrono::duration;
  using std::chrono::duration_cast;
  using std::chrono::high_resolution_clock;
  using std::chrono::milliseconds;

  auto t1 = high_resolution_clock::now();
  MatrixXd x = (A.transpose() * A).ldlt().solve(A.transpose() * b);
  auto t2 = high_resolution_clock::now();

  duration<double, std::milli> ms_double = t2 - t1;
  std::cout << "N =" << x_coords.size() << " runtime: " << ms_double.count() << "ms\n";
  std::vector<double> res(&x(0), x.data() + x.cols() * x.rows());
  return res;
}

int main(int argc, char* argv[]) {

  /// @todo Change implementation s.t. rhe value for N can be supplied via a command line argument
  size_t N = 30000;

  {
    std::random_device rd;
    std::mt19937 gen(1);
    std::normal_distribution<double> dis(0, 0.2);

    auto x = std::vector<double>(N, 0);
    auto y = std::vector<double>(N, 0);
    for (size_t i = 0; i != N; ++i) {
      x[i] = 0 + i * (5.0 - 0) / (N - 1);
      //y = 2*(1 - exp(-x))
      y[i] = 2.0 * (1.0 - std::exp(-x[i])); 
    }

    {
      auto coeffs = poly_fit(x, y, 3);
      std::cout << "[ ";
      for (const auto& coeff : coeffs) {
        std::cout << coeff << " ";
      }
      std::cout << "]" << std::endl;
    }
    for (size_t i = 0; i != N; ++i) {
      y[i] = y[i] + dis(gen);
    }

    {
      auto coeffs = poly_fit(x, y, 3);
      std::cout << "[ ";
      for (const auto& coeff : coeffs) {
        std::cout << coeff << " ";
      }
      std::cout << "]" << std::endl;
    }
    {
      auto coeffs = poly_fit(x, y, 4);
      std::cout << "[ ";
      for (const auto& coeff : coeffs) {
        std::cout << coeff << " ";
      }
      std::cout << "]" << std::endl;
    }
  }
  return 0;
}