#include "board.hpp"
#include "forces.hpp"
#include "math.h"
#include "pedestrian.hpp"
#include "tsc_x86.h"
#include "vector2d.hpp"
#include <fstream>
#include <iostream>
#include <random>
#include <tuple>
#include <vector>

#define SEED 42

#define BOARD_WIDTH 10.0
#define BOARD_LENGTH 50.0

#define BOARD_LENGTH_MINUS_ONE BOARD_LENGTH - 1
#define BOARD_LENGTH_MINUS_ONE_POINT_FIVE BOARD_LENGTH - 1.5
#define WALL_EPS 1e-3
#define BOARD_WIDTH_MINUS_WALL_EPS BOARD_WIDTH - WALL_EPS
#define BOARD_LENGTH_HALF BOARD_LENGTH / 2
#define BOARD_LENGTH_MINUS_TWO BOARD_LENGTH - 2
#define BOARD_WIDTH_MINUS_ZERO_POINT_THREE BOARD_WIDTH - 0.3
#define BOARD_LENGTH_PLUS_THIRTY BOARD_LENGTH + 30

using namespace std;

Pedestrian *init_pedestrians(int n_p) {

  Pedestrian *P = new Pedestrian[n_p];

  // Platform-specific non-deterministic random number generator (if available)
  // random_device rd;
  // mt19937 gen(rd()); // Use rd for true randomness

  mt19937 gen(SEED); // Use fixed seed for reproducibility

  normal_distribution<double> velocity_gen(v0_mean, v0_std_dev);
  uniform_real_distribution<double> x_gen_left(2, BOARD_LENGTH_HALF);
  uniform_real_distribution<double> x_gen_right(BOARD_LENGTH_HALF, BOARD_LENGTH_MINUS_TWO);
  uniform_real_distribution<double> y_gen(0.3, BOARD_WIDTH_MINUS_ZERO_POINT_THREE);

  int n_p_half = n_p / 2;

  for (int i = 0; i < n_p_half; i++) {

    double x = x_gen_left(gen);
    double y = y_gen(gen);

    Vector2d start_pos = Vector2d(x, y);
    Vector2d destination = Vector2d(BOARD_LENGTH_PLUS_THIRTY, y);

    Vector2d init_e = destination - start_pos;
    double init_e_norm = init_e.norm();
    Vector2d init_e_unit = init_e / init_e_norm;

    double desired_v0;

    // Avoid values are too small since they would make the pedestrians walk backwards
    do {
      desired_v0 = velocity_gen(gen);
    } while (desired_v0 <= 0.1);

    double v_max = 1.3 * desired_v0;

    P[i] = Pedestrian(WHITE, init_e_unit * desired_v0, start_pos, init_e_unit, v_max, desired_v0, destination, Vector2d());
  }

  for (int i = n_p_half; i < n_p; i++) {

    double x = x_gen_right(gen);
    double y = y_gen(gen);

    Vector2d start_pos = Vector2d(x, y);
    Vector2d destination = Vector2d(-30, y);

    Vector2d init_e = destination - start_pos;
    double init_e_norm = init_e.norm();
    Vector2d init_e_unit = init_e / init_e_norm;

    double desired_v0;
    do {
      desired_v0 = velocity_gen(gen);
    } while (desired_v0 <= 0.1);

    double v_max = 1.3 * desired_v0;

    P[i] = Pedestrian(BLACK, init_e_unit * desired_v0, start_pos, init_e_unit, v_max, desired_v0, destination, Vector2d(0, 0));
  }

  return P;
}

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

  // Ensure that arguments are passed
  if (argc < 3) {
    printf("Usage: ./sim <number_of_pedestrians> <number_of_iterations>\n");
    printf("Each iteration takes %.2f seconds in the simulation.\n", dt);
    return 1;
  }

  // Parse command-line arguments
  int n_p = atoi(argv[1]);
  int n_iter = atoi(argv[2]);

  if (n_p <= 0) {
    printf("Error: <number_of_pedestrians> must be a positive integer, got '%s' instead\n", argv[1]);
    return 1;
  }

  if (n_iter <= 0) {
    printf("Error: <number_of_iterations> must be a positive integer, got '%s' instead\n", argv[2]);
    return 1;
  }

  Obstacle obs[] = {};
  Board board = Board(BOARD_LENGTH, BOARD_WIDTH, 0, obs);

  // Logic for importing/exporting the simulation data
  const int dump_interval = 1;
  const char *dump_path = (argc >= 4) ? argv[3] : nullptr;

  vector<tuple<int, int, double, double, double, double, int>> samples; // t, id, x, y, e.x, e.y, color
  if (dump_path) {
    samples.reserve((size_t)n_iter / dump_interval * n_p);
  }

  Pedestrian *pedestrians = init_pedestrians(n_p);

  uint64_t start = start_tsc();

  for (int t = 0; t < n_iter; t++) {

    // Pass 1: Conmpute all the forces given the current state
    for (int i = 0; i < n_p; i++) {

      Pedestrian &p = pedestrians[i];

      p.f_i.reset();

      // # 6 ADDs, 2 DIVs
      //
      // $\vec{F}_\alpha(\vec{v}_\alpha, v^0_\alpha \vec{e}_\alpha) = \frac{1}{\tau_\alpha}(v^0_\alpha \vec{e}_\alpha - \vec{v}_\alpha)$
      //
      add_acceleration_term(&p);

      // # 15 ADDs, 24 MULTs, 8 DIVs, 5 SQRTs, 1 EXP
      //
      // $\sum_{\beta} \vec{F}_{\alpha\beta}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{\beta})$
      //
      add_repulsive_effects_from_pedestrians(pedestrians, &p, n_p);

      // # board.n_obstacles * (13 ADDs, 19 MULTs, 4 DIVs, 2 SQRTs, 1 EXP)
      //
      // $\sum_{B} \vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{B})$
      //
      add_repulsive_effects_from_borders(board, &p);
    }

    // Pass 2: Integrate the forces and update the current state
    for (int i = 0; i < n_p; i++) {

      Pedestrian &p = pedestrians[i];

      // # 2 ADDs, 2 MULTs
      //
      // $\frac{d\vec{w_\alpha}}{dt} = \vec{F}_{\alpha}(t)$
      //
      //
      p.w = p.w + p.f_i * dt;

      // # 1 ADD, 2 MULTs, 1 SQRT
      double w_norm = p.w.norm();

      // # Assumption for calculating W(n): The "if" statement always evaluates to true.
      if (w_norm > p.v_max) {
        // # 2 MULTs, 1 DIV
        p.w = p.w * (p.v_max / w_norm);
      }

      // # 1 ADD, 4 MULTs, 1 DIV, 1 SQRT
      p.curr_vel = p.w * g(p.v_max, p.w.norm());

      // # 2 ADDs, 2 MULTs
      //
      // $\vec{r}_i^{t+1} = \vec{r}^t_i + v^t_idt$
      //
      p.curr_pos = p.curr_pos + p.curr_vel * dt;

      // When a pedestrian reaches the end of the hallway we simply "teleport" the person back to the beginning.
      // This way, we have a constant flow of pedestrians, and it is easier to see lanes forming.
      if (p.color == WHITE && p.curr_pos.x > BOARD_LENGTH_MINUS_ONE_POINT_FIVE) {

        // keep y, reset x
        p.curr_pos.x = 1.0;
        p.destination.y = p.curr_pos.y;

        p.w = Vector2d(p.desired_v0, 0);
        p.curr_vel = p.w;
      }

      if (p.color == BLACK && p.curr_pos.x < 1.5) {
        p.curr_pos.x = BOARD_LENGTH_MINUS_ONE;
        p.destination.y = p.curr_pos.y;
        p.w = Vector2d(-p.desired_v0, 0);
        p.curr_vel = p.w;
      }

      // # 2 ADDs
      Vector2d e = p.destination - p.curr_pos;

      // # 1 ADD, 2 MULTs, 1 SQRT
      double e_norm = e.norm();

      // # 2 DIVs
      p.e = e / e_norm;

      if (p.curr_pos.y < 0.0) {
        p.curr_pos.y = WALL_EPS;
        if (p.curr_vel.y < 0)
          p.curr_vel.y = 0;
        if (p.w.y < 0)
          p.w.y = 0;
      }
      if (p.curr_pos.y > BOARD_WIDTH) {
        p.curr_pos.y = BOARD_WIDTH_MINUS_WALL_EPS;
        if (p.curr_vel.y > 0)
          p.curr_vel.y = 0;
        if (p.w.y > 0)
          p.w.y = 0;
      }

      // Save the data for visualization
      if (dump_path && (t % dump_interval == 0)) {
        samples.emplace_back(t, i, p.curr_pos.x, p.curr_pos.y, p.e.x, p.e.y, p.color);
      }
    }
  }

  long double T = stop_tsc(start);

  long long add_count_per_iter = n_p * (21 + board.n_obstacles * 13) + n_p * 9;
  long long mult_count_per_iter = n_p * (24 + board.n_obstacles * 19) + n_p * 14;
  long long div_count_per_iter = n_p * (10 + board.n_obstacles * 4) + n_p * 4;
  long long sqrt_count_per_iter = n_p * (5 + board.n_obstacles * 2) + n_p * 3;
  long long exp_count_per_iter = n_p * (1 + board.n_obstacles * 1);

  long double W_per_iter = add_count_per_iter * cost_add + mult_count_per_iter * cost_mult + div_count_per_iter * cost_div + sqrt_count_per_iter * cost_sqrt + exp_count_per_iter * cost_exp;
  long long Q_per_iter = sizeof(Pedestrian) * n_p + sizeof(Obstacle) * board.n_obstacles; // TODO: Refine this
  long double I = W_per_iter / Q_per_iter;
  long double W = W_per_iter * n_iter;
  long double P = W / T;

  printf("W(n) per iteration: %.2Lf flops\n", W_per_iter);
  printf("Q(n) per iteration: %Ld bytes\n", Q_per_iter);
  printf("I(n): %.2Lf flops/byte\n", I);
  printf("W(n): %.2Lf flops\n", W);
  printf("T(n): %.2Lf cycles\n", T);
  printf("P(n): %.2Lf flops/cycle\n", P);

  // Export the data to a .csv file
  if (dump_path) {
    ofstream out(dump_path);
    out << "time,id,x,y,p.e.x,p.e.y,color\n";
    for (auto &[t, id, x, y, e_x, e_y, c] : samples) {
      out << t << ',' << id << ',' << x << ',' << y << ',' << e_x << ',' << e_y << ',' << c << '\n';
    }
  }

  delete[] pedestrians;
  return 0;
}