reduced flop count

.clang-format

@@ -0,0 +1 @@
+ColumnLimit: 500

.gitignore

@@ -1,14 +1,17 @@
 # Build artifacts
 *.o
 *.out
+*.s
+
+# Runtime data
 *.csv
-*.mp4
 *.gif
-*.log
+*.mp4
-sim
+
-main*
+# Compiled binaries
+main
+sizes
 
 # Editor files
-.clang*
 .vscode
 .DS_Store

Makefile

@@ -0,0 +1,7 @@
+main: board.hpp forces.cpp forces.hpp main.cpp pedestrian.hpp vector2d.hpp
+	g++ -O3 -march=native -std=c++23 -S main.cpp -o main.s
+	g++ -O3 -march=native -std=c++23 -S forces.cpp -o force.s
+	g++ -O3 -march=native -std=c++23 main.cpp forces.cpp -o main
+
+sizes: sizes.cpp
+	g++ -O3 -march=native -std=c++23 sizes.cpp -o sizes

board.hpp

@@ -8,14 +8,14 @@ struct Obstacle {
   Obstacle() : start(0, 0), end(0, 0) {}
   Obstacle(Vector2d s, Vector2d e) : start(s), end(e) {}
 
-  // Find the point of the obstacle that is nearest to p
+  // Find the point of the obstacle that is closest to p
-  // # 15 flops
+  // # 8 ADDs, 6 MULTs, 1 DIV
   Vector2d find_nearest_obstacle_point(Vector2d p) {
 
     // # 2 ADDs
     Vector2d obstacle_v = end - start;
 
-    // # 9 flops =  2 ADDs, 2 MULTs, 1 ADD, 1 DIV, 2 MULTs, 1 ADD
+    // # 4 ADDs, 4 MULTs, 1 DIV
     double projection = (p - start).dot(obstacle_v) / obstacle_v.dot(obstacle_v);
 
     if (projection < 0.0)
@@ -23,17 +23,15 @@ struct Obstacle {
     else if (projection > 1.0)
       projection = 1.0;
 
-    // # 4 flops = 2 ADDs, 2 MULTs
+    // # 2 ADDs, 2 MULTs
     return start + obstacle_v * projection;
-
   }
-
 };
 
 struct Board {
   double length, width;
-  Obstacle *obstacles; // array of obstacles (represented as lines) that are present on the board
+  Obstacle *obstacles; // Array of obstacles, represented as lines on the board
-  int n_obstacles;     // number of obstacles
+  int n_obstacles;     // Number of obstacles
 
   Board(double l, double w, int n, Obstacle *o) : length(l), width(w), n_obstacles(n + 4) {
     obstacles = new Obstacle[n + 4];

forces.cpp

@@ -1,7 +1,7 @@
 #include "forces.hpp"
 #include <math.h>
 
-// # 8 flops = 5 MULTs, 2 ADD, 1 SQRT
+// # 2 ADDs, 5 MULTs, 1 SQRT
 double w(Vector2d e, Vector2d f) {
   if (e.dot(f) >= f.norm() * cos_phi)
     return 1.0;
@@ -9,95 +9,97 @@ double w(Vector2d e, Vector2d f) {
     return c;
 }
 
-// # 2 DIVs
+// # 1 DIV
 double g(double v_max_alpha, double w_alpha_norm) {
+
+  // # Assumption for calculating W(n): The "if" statement always evaluates to false.
   if (w_alpha_norm <= v_max_alpha) {
     return 1;
   }
+
   return v_max_alpha / w_alpha_norm;
 }
 
-// # 8 flops
+// # 6 ADDs, 2 DIVs
-void compute_acceleration_term(Pedestrian *p) {
+void add_acceleration_term(Pedestrian *p) {
-  // # 2 flops
+
+  // # 2 ADDs
   Vector2d desired_velocity = p->e * p->desired_v0;
+
   Vector2d actual_velocity = p->curr_vel;
-  // # 6 flops = 4 ADDs, 2 DIVs
+
+  // # 4 ADDs, 2 DIVs
   p->f_i += (desired_velocity - actual_velocity) / tau_alpha;
 }
 
-// # 55(N - 1) = 55N - 55 flops
+// # 15 ADDs, 24 MULTs, 8 DIVs, 5 SQRTs, 1 EXP
 //
 // $\vec{F}_{\alpha\beta}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{\beta})$
 //
-void compute_repulsive_effects_from_pedestrians(Pedestrian *pedestrians, Pedestrian *p, int N) {
+void add_repulsive_effects_from_pedestrians(Pedestrian *pedestrians, Pedestrian *p, int N) {
 
   for (int i = 0; i < N; i++) {
 
     if (pedestrians + i == p)
       continue;
 
-    // # Optional: Do not consider far away pedestrians from current pedestrian
-    // # Might get better speed up, remove if not
-    // if ((pedestrians[i].curr_pos - p->curr_pos).norm() > DISTANCE) {
-    //   continue;
-    // }
-
     // # 2 ADDs
     // $\vec{r}_{\alpha\beta}$
     Vector2d r_ab = p->curr_pos - pedestrians[i].curr_pos;
 
-    // # 4 flops = 1 SQRT, 2 MULTs, 1 ADD
+    // # 1 ADD, 2 MULTs, 1 SQRT
     // $\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert$
     double r_ab_norm = r_ab.norm();
 
     // $\vec{e}_{\beta}$
     Vector2d e_b = pedestrians[i].e;
 
-    // # 2 MULTs
+    // # 1 ADD, 3 MULTs, 1 SQRT
     // $v_\beta \Delta t$
     double s_b = pedestrians[i].curr_vel.norm() * delta_t;
 
-    // # 4 flops = 2 ADDs, 2 MULTs
+    // # 2 ADDs, 2 MULTs
     // $\vec{r}_{\alpha\beta} - \vec{v}_{\beta}\varDelta_{t}\vec{e}_{\beta}$
     Vector2d r_ab2 = r_ab - (e_b * s_b);
 
-    // # 4 flops = 1 SQRT, 2 MULTs, 1 ADD
+    // # 1 ADD, 2 MULTs, 1 SQRT
     // $\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert$
     double r_ab2_norm = r_ab2.norm();
 
-    // # 7 flops = 1 SQRT, 3 ADDs, 2 MULTs, 1 DIV
+    // # 1 ADD
+    double r_ab_norm_sum = r_ab_norm + r_ab2_norm;
+
+    // # 1 ADD, 2 MULTs, 1 DIV, 1 SQRT
     //
     // $2b = \sqrt{(\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert)^2 - (v_\beta \Delta t)^2}$
     //
-    double b = sqrt((r_ab_norm + r_ab2_norm) * (r_ab_norm + r_ab2_norm) - s_b * s_b) / 2;
+    double b = sqrt(r_ab_norm_sum * r_ab_norm_sum - s_b * s_b) / 2;
 
-    // # 12 flops = 6 DIVs, 3 ADDs, 3 MULTs
+    // # 2 ADDs, 3 MULTs, 6 DIVs
     //
     // $\grad_{\vec{r}_{\alpha \beta}}b = \frac{(\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert) (\frac{\vec{r}_{\alpha \beta}}{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert} + \frac{\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta}{\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert})}{4b}$
     //
-    Vector2d b_grad = ((r_ab / r_ab_norm + r_ab2 / r_ab2_norm) * (r_ab_norm + r_ab2_norm)) / (4 * b);
+    Vector2d b_grad = ((r_ab / r_ab_norm + r_ab2 / r_ab2_norm) * r_ab_norm_sum) / (4 * b);
 
-    // # 6 flops = 3 MULTs, 2 DIVs, 1 EXP
+    // # 3 MULTs, 1 DIV, 1 EXP
     //
     // $\vec{f}_{\alpha \beta}(\vec{r}_{\alpha \beta}) = \frac{V^0_{\alpha \beta}}{\sigma}e^{-b / \sigma} \grad_{\vec{r}_{\alpha \beta}}b$
     //
-    Vector2d f_ped = b_grad * ((V0_AB / SIGMA) * (exp(-b / SIGMA)));
+    Vector2d f_ped = b_grad * (V0_AB_div_SIGMA * exp(-b / SIGMA));
 
-    // # 14 flops: 9 MULTs, 4 ADDs, 1 SQRT
+    // # 4 ADDs, 7 MULTs, 1 SQRT
     //
     // $\vec{F}_{\alpha\beta}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{\beta}) = w(\vec{e}_{\alpha}, -\vec{f}_{\alpha\beta})\vec{f}_{\alpha\beta}(\vec{r}_{\alpha} - \vec{r}_{\beta})$
     //
-    p->f_i += f_ped * w(p->e, f_ped * -1);
+    p->f_i += f_ped * w(p->e, f_ped.reverse_sign());
   }
 }
 
-// # Assume that the four borders are also 4 obstacles
+// # board.n_obstacles * (13 ADDs, 19 MULTs, 4 DIVs, 2 SQRTs, 1 EXP)
-// # 48 n_obstacles flops = 192 flops
 //
 // $\vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{B})$
 //
-void compute_repulsive_effects_from_borders(Board &board, Pedestrian *p) {
+void add_repulsive_effects_from_borders(Board &board, Pedestrian *p) {
 
   Obstacle *obstacles = board.obstacles;
   int n_borders = board.n_obstacles;
@@ -106,26 +108,23 @@ void compute_repulsive_effects_from_borders(Board &board, Pedestrian *p) {
 
     Obstacle obstacle = obstacles[i];
 
-    // # 17 flops
+    // # 8 ADDs, 6 MULTs, 1 DIV
     // $\vec{r}_{\alpha B}$
     Vector2d r_aB = p->curr_pos - obstacle.find_nearest_obstacle_point(p->curr_pos);
 
-    // # 17 flops = 2 DIVs, 1 SQRT, 2 MULTs, 1 ADD, 2 MULTs, 1 DIV, 2 MULTs, 1 EXP, 1 SQRT, 2 MULTs, 1 ADD, 1 DIV
+    // # 1 ADD, 2 MULTs, 1 SQRT
+    double r_aB_norm = r_aB.norm();
+
+    // # 4 MULTs, 3 DIVs, 1 EXP
     //
     // $\vec{F}_{\alpha B}(\vec{r}_{\alpha B}) = \frac{U^0_{\alpha B}}{R} e^{- \Arrowvert \vec{r}_{\alpha B} \Arrowvert / R} \frac{\vec{r}_{\alpha B}}{\Arrowvert \vec{r}_{\alpha B} \Arrowvert}$
     //
-    Vector2d f_obst = (r_aB / r_aB.norm()) * (U0_AB / R) * std::exp(-r_aB.norm() / R);
+    Vector2d f_obst = (r_aB / r_aB_norm) * U0_AB_div_R * exp(-r_aB_norm / R);
 
-    // # 2 ADDs
+    // # 4 ADDS, 7 MULTs, 1 SQRT
-    // p->f_i += f_obst; //commented since this version dont takes the field of view / angle to the border into account
-
-    // I have read the paper again and on page 6 they say
-    // However, the formulas above for attractive and repulsive effects only hold for situations that are perceived in the desired direction e_alpha(t) of motion.
-    // I assumed that with formulas above its also meant for equation 8 (even though not explicitly mentioned their)
-    // # 14 flops: 9 MULTs, 4 ADDS, 1 SQRT
     //
-    // $(12):\vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{ B}) = w(\vec{e}_{\alpha}, -\vec{F}_{\alpha B})\vec{F}_{\alpha B}(\vec{r}_{\alpha} - \vec{r}_{ B})$
+    // $\vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{ B}) = w(\vec{e}_{\alpha}, -\vec{F}_{\alpha B})\vec{F}_{\alpha B}(\vec{r}_{\alpha} - \vec{r}_{ B})$
     //
-    p->f_i += f_obst * w(p->e, f_obst * -1);
+    p->f_i += f_obst * w(p->e, f_obst.reverse_sign());
   }
 }

forces.hpp

@@ -12,35 +12,59 @@ const double V0_AB = 2.1;
 // $\sigma$
 const double SIGMA = 0.3;
 
+//
+// $\frac{V^0_{\alpha\beta}}{\sigma}$
+// 
+const double V0_AB_div_SIGMA =  V0_AB / SIGMA;
+
 // $U^0_{\alpha B}$
 const double U0_AB = 10.0;
 
 // $R$
 const double R = 0.2;
 
+// 
+// $\frac{U^0_{\alpha B}}{R}$
+// 
+const double U0_AB_div_R = U0_AB / R;
+
 // $c$
 const double c = 0.5;
 
 // $\varphi$
 const double phi = 45.0 / 180.0 * M_PI;
 
+// $cos(\varphi)$
 const double cos_phi = cos(phi);
 
+// $<v^0>$
+const double v0_mean = 1.34;
+
+// $\sqrt{\theta}$
+const double v0_std_dev = 0.26;
+
 // $\varDelta_{t}$
-// step size = 0.2 seconds
 const double delta_t = 2;
 
+// $dt$
 const double dt = 0.2;
 
-// # Optional (forces.cpp: compute_repulsive_effects_from_pedestrians)
+const double cost_add = 1;
-// const double DISTANCE = 2.0;
+const double cost_mult = 1;
+const double cost_div = 1;
+
+// https://acl.inf.ethz.ch/teaching/fastcode/knowledgebase/#:~:text=How%20do%20exp(x)%2C%20sin(x)%2C%20etc.%20count%20in%20the%20flopcount%20of%20an%20algorithm%3F
+const double cost_exp = 30;
+
+// -O3 maps each sqrt() call to an instruction
+const double cost_sqrt = 1;
 
 double w(Vector2d e, Vector2d f);
 
 double g(double v_max_alpha, double w_alpha_norm);
 
-void compute_acceleration_term(Pedestrian *p);
+void add_acceleration_term(Pedestrian *p);
 
-void compute_repulsive_effects_from_pedestrians(Pedestrian *P, Pedestrian *p, int N);
+void add_repulsive_effects_from_pedestrians(Pedestrian *P, Pedestrian *p, int N);
 
-void compute_repulsive_effects_from_borders(Board &B, Pedestrian *p);
+void add_repulsive_effects_from_borders(Board &B, Pedestrian *p);

main.cpp

@@ -6,172 +6,217 @@
 #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;
 
-// Numer of pedestrians
-// $n_p$
-// const int n_p = 100;
-
-// Number of iterations
-// $n_{iter}$
-// const int n_iter = 100;
-
 Pedestrian *init_pedestrians(int n_p) {
+
   Pedestrian *P = new Pedestrian[n_p];
-  std::random_device rd; // use for truly random simulation
+
-  // std::mt19937 gen(rd());                                    // use rd instead of 42 as seed for real randomness
+  // Platform-specific non-deterministic random number generator (if available)
-  std::mt19937 gen(42);                                      // use fixed seed for reproducibility
+  // random_device rd;
-  std::normal_distribution<double> velocity_gen(1.34, 0.26); // 1.34 and 0.26 are the mean and stdrd deviation used in paper
+  // mt19937 gen(rd()); // Use rd for true randomness
-  std::uniform_real_distribution<double> x_gen_left(2, BOARD_LENGTH / 2);
+
-  std::uniform_real_distribution<double> x_gen_right(BOARD_LENGTH / 2, BOARD_LENGTH - 2);
+  mt19937 gen(SEED); // Use fixed seed for reproducibility
-  std::uniform_real_distribution<double> y_gen(0.3, BOARD_WIDTH - 0.3);
+
-  for (int i = 0; i < n_p / 2; i++) {
+  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_gen_left(gen), y);
+
-    Vector2d destination = Vector2d(BOARD_LENGTH + 30, y);
+    Vector2d start_pos = Vector2d(x, y);
-    // double desired_v0 = velocity_gen(gen);
+    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); // to avoid small values (Way below the social force -> that would make the walk backwards) (we discovered some weird stuff in the simulation)
+    } while (desired_v0 <= 0.1);
+
     double v_max = 1.3 * desired_v0;
-    P[i] = Pedestrian(WHITE, (destination - start_pos) / (destination - start_pos).norm() * desired_v0, start_pos, (destination - start_pos) / (destination - start_pos).norm(), v_max, desired_v0, destination, Vector2d());
+
+    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 / 2; i < n_p; i++) {
+  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_gen_right(gen), y);
+
+    Vector2d start_pos = Vector2d(x, y);
     Vector2d destination = Vector2d(-30, y);
-    // double desired_v0 = velocity_gen(gen);
+
+    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, (destination - start_pos) / (destination - start_pos).norm() * desired_v0, start_pos, (destination - start_pos) / (destination - start_pos).norm(), v_max, desired_v0, destination, Vector2d(0, 0));
+
+    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[]) {
 
-  // Check command line arguments
+  // Ensure that arguments are passed
   if (argc < 3) {
     printf("Usage: ./sim <number_of_pedestrians> <number_of_iterations>\n");
-    printf("  number_of_iterations: each iteration = %.2f seconds of simulation time\n", dt);
+    printf("Each iteration takes %.2f seconds in the simulation.\n", dt);
     return 1;
   }
 
-  // Parse command line arguments
+  // 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'\n", argv[1]);
+    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'\n", argv[2]);
+    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 to get/dump the simulation data
+  // Logic for importing/exporting the simulation data
-  const int dump_interval = 1; // dump every 5 iterations -> 20 fps (since dt = 0.01)
+  const int dump_interval = 1;
   const char *dump_path = (argc >= 4) ? argv[3] : nullptr;
 
-  std::vector<std::tuple<int, int, double, double, double, double, int>> samples; // t, id, x, y, e.x, e.y, color
+  vector<tuple<int, int, double, double, double, double, int>> samples; // t, id, x, y, e.x, e.y, color
-  if (dump_path)
+  if (dump_path) {
     samples.reserve((size_t)n_iter / dump_interval * n_p);
+  }
 
   Pedestrian *pedestrians = init_pedestrians(n_p);
 
-  // # Start here or right after function signature
   uint64_t start = start_tsc();
 
-  // # (8 + 55N - 55 + 192 + 4 + 8 + 4 + 10)*n_p*n_iter
   for (int t = 0; t < n_iter; t++) {
-    // Pass1: conmpute all forces from the current state
+
+    // 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();
 
-      // # 8 flops
+      // # 6 ADDs, 2 DIVs
       //
-      // $\vec{F}_\alpha(\vec{v}_\alpha, v^0_\alpha \vec{e}_\alpha) \vec{e}_\alpha(t) = \frac{1}{\tau_\alpha}(v^0_\alpha \vec{e}_\alpha - \vec{v}_\alpha)$
+      // $\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)$
       //
-      compute_acceleration_term(&p);
+      add_acceleration_term(&p);
 
-      // # 55N - 55 flops
+      // # 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})$
       //
-      compute_repulsive_effects_from_pedestrians(pedestrians, &p, n_p);
+      add_repulsive_effects_from_pedestrians(pedestrians, &p, n_p);
 
-      // # 192 flops
+      // # 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})$
       //
-      compute_repulsive_effects_from_borders(board, &p);
+      add_repulsive_effects_from_borders(board, &p);
     }
-    // Pass2: integrate the forces and update the current state
+
+    // Pass 2: Integrate the forces and update the current state
     for (int i = 0; i < n_p; i++) {
+
       Pedestrian &p = pedestrians[i];
-      // # 4 flops = 2 ADDs, 2 MULTs
+
-      // integrate preferred velocity:  dw/dt = F           (Eq. 10)
+      // # 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);
       }
-      // Vector2d w_alpha = p.curr_vel + p.f_i * dt; // commented since this was using p.curr_vel instead of p.w (different from paper)
 
-      // # 8 flops = 2 MULTs, 2 DIVs, 1 SQRT, 2 MULTs, 1 ADD
+      // # 1 ADD, 4 MULTs, 1 DIV, 1 SQRT
       p.curr_vel = p.w * g(p.v_max, p.w.norm());
 
-      // # 4 flops = 2 ADDs, 2 MULTs
+      // # 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;
 
-      // # 10 flops = 5 ADDs, 2 DIVs, 1 SQRT, 2 MULTs
+      // When a pedestrian reaches the end of the hallway we simply "teleport" the person back to the beginning.
-      // p.e = (p.destination - p.curr_pos) / (p.destination - p.curr_pos).norm();
+      // 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) {
 
-      // how to handle if a pedestrian reached the end of our simulation/end of the hallway
+        // keep y, reset x
-      // The option we do here is to just "teleport" the person back to the beginning
-      // the good thing is we will have a constant flow of pedestrians and in the simulation we should see how lines build up
-      if (p.color == WHITE && p.curr_pos.x > (BOARD_LENGTH - 1.5)) {
-        // basically just keep the y position and reset the 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 - 1.0;
+        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;
       }
 
-      // # 10 flops = 5 ADDs, 2 DIVs, 1 SQRT, 2 MULTs
+      // # 2 ADDs
-      p.e = (p.destination - p.curr_pos) / (p.destination - p.curr_pos).norm();
+      Vector2d e = p.destination - p.curr_pos;
-      // p.e = (p.curr_vel - p.destination) / (p.curr_vel - p.destination).norm();
+
-      // p.e = (p.destination + (p.curr_vel / p.curr_vel.norm())) / (p.destination + (p.curr_vel / p.curr_vel.norm())).norm();
+      // # 1 ADD, 2 MULTs, 1 SQRT
-      const double WALL_EPS = 1e-3;
+      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)
@@ -180,53 +225,44 @@ int main(int argc, const char *argv[]) {
           p.w.y = 0;
       }
       if (p.curr_pos.y > BOARD_WIDTH) {
-        p.curr_pos.y = BOARD_WIDTH - WALL_EPS;
+        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;
       }
 
-      // if we have a dump path/file we wanna save the data for visualisation
+      // 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 long W = (long long)(8 + 55 * n_p - 55 + 192 + 4 + 8 + 4 + 10) * n_p * n_iter;
-  std::cout << "W(n): ";
-  std::cout << W;
-  std::cout << " flops\n";
 
-  double T = stop_tsc(start);
+  long double T = stop_tsc(start);
-  std::cout << "T(n): ";
-  std::cout << T;
-  std::cout << " cycles\n";
 
-  double P = W / T;
+  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);
 
-  std::cout << "P(n): ";
+  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;
-  std::cout << P;
+  long long Q_per_iter = sizeof(Pedestrian) * n_p + sizeof(Obstacle) * board.n_obstacles; // TODO: Refine this
-  std::cout << " flops/cycle\n";
+  long double I = W_per_iter / Q_per_iter;
+  long double W = W_per_iter * n_iter;
+  long double P = W / T;
 
-  long long W_per_iter = (long long)(8 + 55 * n_p - 55 + 192 + 4 + 8 + 4 + 10) * n_p;
+  printf("W(n) per iteration: %.2Lf flops\n", W_per_iter);
-  std::cout << "W(n) per iteration: ";
+  printf("Q(n) per iteration: %Ld bytes\n", Q_per_iter);
-  std::cout << W_per_iter;
+  printf("I(n): %.2Lf flops/byte\n", I);
-  std::cout << " flops\n";
+  printf("W(n): %.2Lf flops\n", W);
+  printf("T(n): %.2Lf cycles\n", T);
+  printf("P(n): %.2Lf flops/cycle\n", P);
 
-  long long Q = 104 * n_p + 32 * 4 + 80;
+  // Export the data to a .csv file
-  std::cout << "Q(n) per iteration: ";
-  std::cout << Q;
-  std::cout << " bytes\n";
-
-  long long I = (long long)(8 + 55 * n_p - 55 + 192 + 4 + 8 + 4 + 10) * n_p / Q;
-  std::cout << "I(n): ";
-  std::cout << I;
-  std::cout << " flops/byte\n";
-
-  // now we wanna save the data to a csv
   if (dump_path) {
-    std::ofstream out(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';

pedestrian.hpp

@@ -1,5 +1,4 @@
 #pragma once
-#include "random"
 #include "vector2d.hpp"
 
 // For visual simulation
@@ -14,7 +13,7 @@ struct Pedestrian {
   // $\vec{r}_\alpha(t)$
   Vector2d curr_pos;
   // $\vec{e}_{\alpha}(t)$
-  Vector2d e;  // desired direction vector
+  Vector2d e; // desired direction vector
   // $v_{\alpha}^{max} = 1.3 v^{0}_{\alpha}$
   double v_max;
   // $v_\alpha^0$
@@ -24,9 +23,14 @@ struct Pedestrian {
   // social force
   Vector2d f_i;
 
-  Pedestrian() : color(WHITE), curr_vel(0.0, 0.0), w(0.0, 0.0), curr_pos(0.0, 0.0), e(0.0, 0.0), v_max(0.0), desired_v0(0.0), destination(0.0, 0.0), f_i(0.0, 0.0) {}
+  Pedestrian()
+      : color(WHITE), curr_vel(0.0, 0.0), w(0.0, 0.0), curr_pos(0.0, 0.0),
+        e(0.0, 0.0), v_max(0.0), desired_v0(0.0), destination(0.0, 0.0),
+        f_i(0.0, 0.0) {}
 
-  Pedestrian(Color clr, Vector2d c_v, Vector2d c_p, Vector2d e, double v_max, double desired_v0, Vector2d d, Vector2d f)
+  Pedestrian(Color clr, Vector2d c_v, Vector2d c_p, Vector2d e, double v_max,
+             double desired_v0, Vector2d d, Vector2d f)
       : color(clr), curr_vel(c_v), w(c_v), // initialize w = v at t=0
-        curr_pos(c_p), e(e), v_max(v_max), desired_v0(desired_v0), destination(d), f_i(f) {}
+        curr_pos(c_p), e(e), v_max(v_max), desired_v0(desired_v0),
+        destination(d), f_i(f) {}
 };

questions.md

@@ -1,6 +0,0 @@
-- [ ] SQRT == 1 flop? EXP == 1 flop?
-- [ ] Where to start measurements, do we include measurements?
-- [ ] Can Vector2d be reused?
-- [ ] How precise do we have to be? How do we define "correctness"?
-- [ ] Do borders count as obstacles?
-- [ ] Do we use a fixed machine, frequency scaling, etc.

vector2d.hpp

@@ -15,9 +15,13 @@ struct Vector2d {
     *this = *this + v;
     return *this;
   }
-  // # 4 flops = 1 SQRT, 2 MULTs, 1 ADD
+
+  Vector2d reverse_sign() { return Vector2d(-x, -y); }
+
+  // # 1 ADD, 2 MULTs, 1 SQRT
   double norm() { return sqrt(x * x + y * y); }
-  // # 3 flops = 2 MULTs, 1 ADD
+
+  // # 1 ADD, 2 MULTs
   double dot(Vector2d v) { return x * v.x + y * v.y; }
   void reset() {
     x = 0;