Merge branch 'elliot/float' into 'main'

Baseline Flop Count Update, New Float Baseline

See merge request course-asl-fastcode26/team70!28

Makefile

@@ -49,7 +49,13 @@ VARIANT ?= baseline
 
 VARIANT_SRC := ./variants/variant_$(VARIANT).cpp
 
-SRC := main.cpp $(VARIANT_SRC)
+ifneq ($(filter float_%,$(VARIANT)),)
+MAIN_SRC := main-float.cpp
+else
+MAIN_SRC := main.cpp
+endif
+
+SRC := $(MAIN_SRC) $(VARIANT_SRC)
 HEADERS := ./headers/board.hpp ./headers/pedestrian.hpp ./headers/vector2d.hpp ./headers/constants.hpp ./headers/util.hpp
 
 TEST_SRC := ./tests/test_vector2d.cpp ./tests/test_board.cpp ./tests/test_forces.cpp $(VARIANT_SRC)
@@ -57,7 +63,6 @@ TEST_TARGET := run_tests
 TEST_LIBS := -lgtest -lgtest_main -pthread
 
 TARGET := main_$(VARIANT)
-SIZES_TARGET := sizes
 
 
 # Default simulation parameters and output files.
@@ -94,12 +99,18 @@ BENCH_SIZES ?= 50 100 200 500 1000 2000 5000 10000
 BENCH_ITER ?= 1000
 
 # Verification parameters
+
+ifneq ($(filter float_%,$(VARIANT)),)
+REFERENCE_CSV ?= reference/baseline_float_n60_i4000.csv
+else
 REFERENCE_CSV ?= reference/baseline_n60_i4000.csv
+endif
+
 REFERENCE_N_P ?= 60
 REFERENCE_N_ITER ?= 4000
 VERIFY_CSV := verify_$(VARIANT).csv
 
-.PHONY: all main run bench run-bench bench-suite video sizes asm test clean clean-all verify
+.PHONY: all main run bench run-bench bench-suite video asm test clean clean-all verify
 
 all: main
 
@@ -147,10 +158,6 @@ bench-suite: bench
 video: run
 	python ./scripts/visualization.py $(CSV) $(MP4)
 
-# make sizes : compile the size-check helper.
-sizes: clean
-	$(CXX) $(CXXFLAGS) sizes.cpp -o $(SIZES_TARGET)
-
 # make test : compile and run all unit tests.
 test: clean
 	$(CXX) $(CXXFLAGS) $(TEST_SRC) $(TEST_LIBS) -o $(TEST_TARGET)
@@ -158,7 +165,7 @@ test: clean
 
 # make asm : compile and generate assembly files.
 
-ASM_CPP_SRCS := main.cpp $(VARIANT_SRC) sizes.cpp
+ASM_CPP_SRCS := main.cpp main-float.cpp $(VARIANT_SRC)
 ASM_HPP_SRCS := ./headers/util.hpp
 ASM_SRCS := $(ASM_CPP_SRCS) $(ASM_HPP_SRCS)
 
@@ -174,7 +181,7 @@ asm: clean
 
 # make clean : delete generated files from previous builds.
 clean:
-	rm -f main_* $(SIZES_TARGET) $(TEST_TARGET) *.o *.s
+	rm -f main_* $(TEST_TARGET) *.o *.s
 
 # make clean-all : delete all generated files from previous builds and runs.
 # WARNING: this deletes all *.csv and *.mp4 files in this directory.

headers/board_float.hpp

@@ -0,0 +1,97 @@
+#pragma once
+#include "vector2d_float.hpp"
+
+struct Obstacle {
+  Vector2d start;       // First endpoint of the obstacle segment.
+  Vector2d end;         // Second endpoint of the obstacle segment.
+  Vector2d obstacle_v;  // Direction vector of the segment: end - start.
+  float length_sq_inv; // Inverse squared length of the segment, used for projection.
+
+  // Default constructor: creates a degenerate obstacle at the origin.
+  // This is mainly used for array allocation before assigning real obstacles.
+  Obstacle() : start(0, 0), end(0, 0), obstacle_v(0, 0), length_sq_inv(0.0) {}
+
+  // Constructs an obstacle as a finite line segment from s to e.
+  //
+  // The segment direction and inverse squared length are precomputed and cached because
+  // they are reused when projecting pedestrian positions onto the obstacle.
+  //
+  // @param s The first endpoint of the obstacle segment.
+  // @param e The second endpoint of the obstacle segment.
+  //
+  Obstacle(Vector2d s, Vector2d e) : start(s), end(e) {
+    obstacle_v = end - start;
+    float length_sq = obstacle_v.dot(obstacle_v);
+    length_sq_inv = length_sq != 0.0 ? 1.0 / length_sq : 0.0;
+  }
+
+  // Flop count: 5 ADDs, 5 MULTs
+  //
+  // Finds the point on this finite obstacle segment that is closest to p.
+  //
+  // The method projects the vector (p - start) onto the obstacle direction
+  // vector, then clamps the projection factor to [0, 1] so that the returned
+  // point lies on the finite segment.
+  //
+  // Note : If the obstacle is a single point, meaning start == end, length_sq_inv is zero
+  // and the method returns start.
+  //
+  // @param p The point for which the closest obstacle point is computed.
+  // @return The closest point on the obstacle segment to p.
+  //
+  Vector2d find_nearest_obstacle_point(const Vector2d &p) const {
+
+    // # 3 ADDs, 3 MULTs
+    float projection = (p - start).dot(obstacle_v) * length_sq_inv;
+
+    if (projection < 0.0)
+      projection = 0.0;
+    else if (projection > 1.0)
+      projection = 1.0;
+
+    // # 2 ADDs, 2 MULTs
+    return start + obstacle_v * projection;
+  }
+};
+
+struct Board {
+  float length;       // Board length.
+  float width;        // Board width.
+  Obstacle *obstacles; // Array of obstacles contained in the board.
+  int n_obstacles;     // Total number of obstacles
+
+  // Constructs a rectangular board of size length x width.
+  //
+  // The input array o contains n user-defined obstacles. The constructor copies
+  // those obstacles and then appends the two horizontal borders of the board (bottom and top) as additional
+  // obstacle segments:
+  //
+  //   bottom: (0, 0) -> (length, 0)
+  //   top:    (0, width) -> (length, width)
+  //
+  // Note: only the horizontal borders are modeled as obstacles. The vertical
+  // borders are left open so pedestrians can enter and leave the walkway without
+  // being affected by repulsive forces from those borders. This is consistent
+  // with Fig. 2 in the paper, where the horizontal walkway boundaries are
+  // emphasized but the vertical ends are not shown as walls.
+  //
+  // @param l The board length.
+  // @param w The board width.
+  // @param n The number of user-defined obstacles in o.
+  // @param o Array of user-defined obstacles to copy into the board obstacles array.
+  //
+  Board(float l, float w, int n, Obstacle *o) : length(l), width(w), n_obstacles(n + 2) {
+    obstacles = new Obstacle[n + 2];
+
+    for (int i = 0; i < n; i++) {
+      obstacles[i] = o[i];
+    }
+
+    // Board borders: bottom, top.
+    obstacles[n] = Obstacle(Vector2d(0, 0), Vector2d(l, 0));
+    obstacles[n + 1] = Obstacle(Vector2d(0, w), Vector2d(l, w));
+  }
+
+  // Releases the dynamically allocated obstacle array.
+  ~Board() { delete[] obstacles; }
+};

headers/board_reduced.hpp

@@ -0,0 +1,44 @@
+// # Replaced Vector2d with double
+// # Predecessor: board.hpp
+
+#pragma once
+
+struct Obstacle {
+  double start_x;
+  double start_y;
+  double end_x;
+  double end_y;
+  double obstacle_v_x;
+  double obstacle_v_y;
+  double length_sq_inv;
+
+  Obstacle() : start_x(0), start_y(0), end_x(0), end_y(0), obstacle_v_x(0), obstacle_v_y(0), length_sq_inv(0.0) {}
+
+  Obstacle(double s_x, double s_y, double e_x, double e_y) : start_x(s_x), start_y(s_y), end_x(e_x), end_y(e_y) {
+    obstacle_v_x = end_x - start_x;
+    obstacle_v_y = end_y - start_y;
+
+    double length_sq = obstacle_v_x * obstacle_v_x + obstacle_v_y * obstacle_v_y;
+    length_sq_inv = length_sq != 0.0 ? 1.0 / length_sq : 0.0;
+  }
+};
+
+struct Board {
+  double length;
+  double width;
+  Obstacle *obstacles;
+  int n_obstacles;
+
+  Board(double l, double w, int n, Obstacle *o) : length(l), width(w), n_obstacles(n + 2) {
+    obstacles = new Obstacle[n + 2];
+
+    for (int i = 0; i < n; i++) {
+      obstacles[i] = o[i];
+    }
+
+    obstacles[n] = Obstacle(0, 0, l, 0);
+    obstacles[n + 1] = Obstacle(0, w, l, w);
+  }
+
+  ~Board() { delete[] obstacles; }
+};

headers/board_reduced_float.hpp

@@ -0,0 +1,44 @@
+// # Replaced double with float
+// # Predecessor: board_reduced.hpp
+
+#pragma once
+
+struct Obstacle {
+  float start_x;
+  float start_y;
+  float end_x;
+  float end_y;
+  float obstacle_v_x;
+  float obstacle_v_y;
+  float length_sq_inv;
+
+  Obstacle() : start_x(0), start_y(0), end_x(0), end_y(0), obstacle_v_x(0), obstacle_v_y(0), length_sq_inv(0.0) {}
+
+  Obstacle(float s_x, float s_y, float e_x, float e_y) : start_x(s_x), start_y(s_y), end_x(e_x), end_y(e_y) {
+    obstacle_v_x = end_x - start_x;
+    obstacle_v_y = end_y - start_y;
+
+    float length_sq = obstacle_v_x * obstacle_v_x + obstacle_v_y * obstacle_v_y;
+    length_sq_inv = length_sq != 0.0 ? 1.0 / length_sq : 0.0;
+  }
+};
+
+struct Board {
+  float length;
+  float width;
+  Obstacle *obstacles;
+  int n_obstacles;
+
+  Board(float l, float w, int n, Obstacle *o) : length(l), width(w), n_obstacles(n + 2) {
+    obstacles = new Obstacle[n + 2];
+
+    for (int i = 0; i < n; i++) {
+      obstacles[i] = o[i];
+    }
+
+    obstacles[n] = Obstacle(0, 0, l, 0);
+    obstacles[n + 1] = Obstacle(0, w, l, w);
+  }
+
+  ~Board() { delete[] obstacles; }
+};

headers/constants_float.hpp

@@ -0,0 +1,105 @@
+#pragma once
+#include <cmath>
+#include <new>
+
+//
+// $\frac{1}{\tau_{\alpha}}$; $\tau_{\alpha} = 0.5s$
+//
+const float tau_alpha_inv = 2.0;
+
+// $V^0_{\alpha \beta}$
+const float V0_AB = 2.1;
+
+// $\sigma$
+const float SIGMA = 0.3;
+
+//
+// $\frac{-1}{2\sigma}$
+//
+const float NEG_TWO_SIGMA_INV = -1.0 / (2 * SIGMA);
+
+//
+// $\frac{-V^0_{\alpha\beta}}{2\sigma}$
+//
+const float NEG_V0_AB_div_TWO_SIGMA = V0_AB * NEG_TWO_SIGMA_INV;
+
+// $U^0_{\alpha B}$
+const float U0_AB = 10.0;
+
+// $R$
+const float R = 0.2;
+
+//
+// $\frac{-1}{R}$
+//
+const float NEG_R_INV = -1.0 / R;
+
+//
+// $\frac{-U^0_{\alpha B}}{R}$
+//
+const float NEG_U0_AB_div_R = -U0_AB / R;
+
+// $-c$
+const float c_neg = -0.5;
+
+// $\varphi$
+const float phi = 100.0 / 180.0 * M_PI;
+
+// $cos(\varphi)$
+const float cos_phi = cos(phi);
+
+//
+// $\frac{1}{cos(\varphi)}$
+//
+const float cos_phi_inv = 1 / cos(phi);
+
+// $<v^0>$
+const float v0_mean = 1.34;
+
+// $\sqrt{\theta}$
+const float v0_std_dev = 0.26;
+
+// $\varDelta_{t}$
+const float delta_t = 2;
+
+// $\varDelta_{t}^2$
+const float delta_t_sq = delta_t * delta_t;
+
+// $dt$
+const float dt = 0.02;
+
+// $-dt$
+const float dt_neg = -dt;
+
+// $-c \cdot dt$
+const float c_neg_times_dt = c_neg * dt;
+
+//
+// $\frac{1}{\tau_{\alpha}} \cdot dt$
+//
+const float dt_times_tau_alpha_inv = dt * tau_alpha_inv;
+
+const float cost_add = 1;
+const float cost_mult = 1;
+const float 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 float cost_exp = 30;
+
+// -O3 maps each sqrt() call to an instruction
+const float cost_sqrt = 1;
+
+// Numerical tolerance used to avoid division by zero and unstable operations
+// when distances or norms are extremely small.
+const float EPS = 1e-9;
+
+// Squared numerical tolerance
+const float EPS_SQ = EPS * EPS;
+
+// Cutoff for pedestrian-pedestrian repulsion
+const float CUTOFF = 12.0;
+const float CUTOFF_SQ = CUTOFF * CUTOFF;
+
+const float cos_phi_sq = cos_phi * cos_phi;
+
+const size_t CACHE_BLOCK_SIZE = std::hardware_destructive_interference_size;

headers/pedestrian_float.hpp

@@ -0,0 +1,65 @@
+#pragma once
+#include "constants_float.hpp"
+#include "vector2d_float.hpp"
+
+// Needed for visual simulation, in order to distinguish two groups of pedestrians:
+//    - White pedestrians move from left to right
+//    - Black pedestrians move from right to left
+enum Color { WHITE, BLACK};
+
+struct Pedestrian {
+
+  Color color;
+
+  // $\vec{v}_{\alpha}(t)$
+  Vector2d curr_vel; // v_alpha  (actual, ≤ v_max)
+
+  // $\vec{w}_{\alpha}$
+  Vector2d w; // w_alpha  (preferred, integrates force)
+
+  // $\vec{r}_\alpha(t)$
+  Vector2d curr_pos;
+
+  // $\vec{e}_{\alpha}(t)$
+  Vector2d e; // desired direction vector
+
+  // $v_{\alpha}^{max} = 1.3 v^{0}_{\alpha}$
+  float v_max;
+
+  // $v_\alpha^0$
+  float desired_v0;
+
+  // $\vec{r}^0_\alpha$
+  Vector2d destination;
+  // $s_\beta^2 = (v_\beta \Delta t)^2$
+  float s_b_sq;
+  // $s_\beta \vec{e}_\beta$
+  Vector2d s_b_e;
+
+  // Default constructor: initializes a pedestrian with color WHITE and all vector/scalar fields set to zero.
+  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), s_b_sq(0.0), s_b_e(0.0, 0.0) {}
+
+  // Constructs a pedestrian with the given initial state and precomputes the
+  // quantities used in the pedestrian-pedestrian repulsion term.
+  //
+  // The preferred velocity w is initialized to the current velocity c_v. The
+  // cached quantities s_b_sq and s_b_e correspond to the predicted step of this
+  // pedestrian over $\Delta t$:
+  //
+  //   s_b     = ||c_v|| * dt
+  //   s_b_sq  = s_b^2
+  //   s_b_e   = e * s_b
+  //
+  // These values are used when another pedestrian computes the elliptical
+  // interaction term involving this pedestrian.
+  //
+  // @param clr The display/color category of the pedestrian.
+  // @param c_v The initial current velocity.
+  // @param c_p The initial current position.
+  // @param e The initial desired direction of motion.
+  // @param v_max The maximum speed that the pedestrian can reach.
+  // @param desired_v0 The desired walking speed.
+  // @param d The pedestrian's target destination point.
+  //
+  Pedestrian(Color clr, Vector2d c_v, Vector2d c_p, Vector2d e, float v_max, float desired_v0, Vector2d d) : color(clr), curr_vel(c_v), w(c_v), curr_pos(c_p), e(e), v_max(v_max), desired_v0(desired_v0), destination(d), s_b_sq(c_v.dot(c_v) * delta_t_sq), s_b_e(e * c_v.norm() * delta_t) {}
+};

headers/pedestrian_reduced_float.hpp

@@ -0,0 +1,38 @@
+// # Replaced double with float
+// # Predecessor: pedestrian_reduced.hpp
+
+#pragma once
+#include "constants_float.hpp"
+
+struct Pedestrian {
+
+  int color;
+
+  float curr_vel_x;
+  float curr_vel_y;
+
+  float w_x;
+  float w_y;
+
+  float curr_pos_x;
+  float curr_pos_y;
+
+  float e_x;
+  float e_y;
+
+  float v_max;
+
+  float desired_v0;
+
+  float destination_x;
+  float destination_y;
+
+  float s_b_sq;
+  float s_b_e_x;
+  float s_b_e_y;
+
+  Pedestrian() : color(0), curr_vel_x(0.0), curr_vel_y(0), w_x(0.0), w_y(0.0), curr_pos_x(0.0), curr_pos_y(0.0), e_x(0.0), e_y(0.0), v_max(0.0), desired_v0(0.0), destination_x(0.0), destination_y(0.0), s_b_sq(0.0), s_b_e_x(0.0), s_b_e_y(0.0) {}
+
+  Pedestrian(int clr, float c_v_x, float c_v_y, float c_p_x, float c_p_y, float e_x, float e_y, float v_max, float desired_v0, float d_x, float d_y)
+      : color(clr), curr_vel_x(c_v_x), curr_vel_y(c_v_y), w_x(c_v_x), w_y(c_v_y), curr_pos_x(c_p_x), curr_pos_y(c_p_y), e_x(e_x), e_y(e_y), v_max(v_max), desired_v0(desired_v0), destination_x(d_x), destination_y(d_y), s_b_sq((c_v_x * c_v_x + c_v_y * c_v_y) * delta_t_sq), s_b_e_x(e_x * (sqrt(c_v_x * c_v_x + c_v_y * c_v_y)) * delta_t), s_b_e_y(e_y * (sqrt(c_v_x * c_v_x + c_v_y * c_v_y)) * delta_t) {}
+};

headers/social_force_float.hpp

@@ -0,0 +1,34 @@
+#pragma once
+
+#include "board_float.hpp"
+#include "pedestrian_float.hpp"
+
+// Initializes the pedestrians used in the simulation.
+//
+// The returned array contains n_p pedestrians split into two groups:
+// one group walking from left to right and one group walking from right to left.
+// Pedestrian positions, desired speeds, destinations, and initial directions are
+// initialized according to the simulation setup.
+//
+// @param n_p : Number of pedestrians to initialize.
+// @return Pointer to a dynamically allocated array of pedestrians.
+Pedestrian *init_pedestrians(int n_p);
+
+// Initializes a simulation board with no obstacles inside.
+//
+// @return Initialized Board instance.
+Board init_board_without_obstacles();
+
+// Runs the social force model with n_p pedestrians for n_iter iterations,
+// i.e., n_iter * dt simulated seconds.
+// If dump_path is not nullptr, the simulation states are written to that CSV
+// file, which can be used to generate a visualization video. If dump_path is
+// nullptr, no output file is generated.
+//
+// @param pedestrians : Array of pedestrians to simulate. The array is updated in place.
+// @param board : Board containing the simulation dimensions and obstacles.
+// @param n_p : Number of pedestrians.
+// @param n_iter : Number of simulation iterations.
+// @param dump_path : Optional path to the CSV output file.
+// @return 0 on success.
+int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path);

headers/vector2d_float.hpp

@@ -0,0 +1,34 @@
+#pragma once
+#include <math.h>
+
+struct Vector2d {
+  float x, y;
+  Vector2d(float a, float b) : x(a), y(b) {}
+  Vector2d() : x(.0), y(.0) {}
+  Vector2d operator+(Vector2d v) const { return Vector2d(v.x + x, v.y + y); }
+  Vector2d operator-(Vector2d v) const { return Vector2d(x - v.x, y - v.y); }
+  Vector2d operator-(float a) const { return Vector2d(x - a, y - a); }
+  Vector2d operator+(float a) const { return Vector2d(x + a, y + a); }
+  Vector2d operator/(float a) const { return Vector2d(x / a, y / a); }
+  Vector2d operator*(float a) const { return Vector2d(x * a, y * a); }
+  Vector2d operator+=(const Vector2d &v) {
+    *this = *this + v;
+    return *this;
+  }
+  Vector2d operator-=(const Vector2d &v) {
+    *this = *this - v;
+    return *this;
+  }
+
+  Vector2d reverse_sign() { return Vector2d(-x, -y); }
+
+  // # 1 ADD, 2 MULTs, 1 SQRT
+  float norm() { return sqrt(x * x + y * y); }
+
+  // # 1 ADD, 2 MULTs
+  float dot(Vector2d v) { return x * v.x + y * v.y; }
+  void reset() {
+    x = 0;
+    y = 0;
+  }
+};

main-float.cpp

@@ -0,0 +1,282 @@
+#include "./headers/board_float.hpp"
+#include "./headers/constants_float.hpp"
+#include "./headers/pedestrian_float.hpp"
+#include "./headers/social_force_float.hpp"
+
+#include <algorithm>
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+
+#ifndef BENCHMARK
+#define BENCHMARK 0
+#endif
+
+#ifndef ROOFLINE
+#define ROOFLINE 0
+#endif
+
+#ifdef __x86_64__
+#include "./headers/tsc_x86.h"
+#include <immintrin.h>
+#define TRASH_BUFFER_SIZE (16ULL * 1024 * 1024 * 2) // 2x L3 (~16MB) on the x86 machine we use for measurements
+#define PAGE_SIZE 4096
+#endif
+
+#ifdef __aarch64__
+#include "./headers/kperf.h"
+#include <arm_acle.h>
+#define TRASH_BUFFER_SIZE (1502838784ULL * 2) // 2x SLC (~1.5GB) on M3 Pro
+#define PAGE_SIZE 16384
+#endif
+
+#define N_OBSTACLES 2
+
+// Benchmark parameters (mostly copied from mvm.c in hw1)
+#define WARMUP_ITERS 3
+#define NUM_REPS 15
+#define NUM_RUNS 1
+#define CYCLES_REQUIRED 1e8
+#define MAX_NUM_RUNS (1 << 14)
+
+// Flushes caches by streaming reads through a large persistent buffer.
+long long flush_cache() {
+  static volatile char *trash_buffer = nullptr;
+  static bool initialized = false;
+
+  if (!initialized) {
+    trash_buffer = new volatile char[TRASH_BUFFER_SIZE];
+
+    // Touch one byte per page once so pages are physically mapped. And we avoid paying the cost of page faults during cache trashing below.
+    for (size_t i = 0; i < TRASH_BUFFER_SIZE; i += PAGE_SIZE) {
+      trash_buffer[i] = 1;
+    }
+
+    initialized = true;
+  }
+
+  long long acc = 0;
+
+  for (size_t i = 0; i < TRASH_BUFFER_SIZE; i += CACHE_BLOCK_SIZE) {
+    // Access one byte per cache block is enough to populate the entire block
+    acc += trash_buffer[i];
+  }
+
+#ifdef __x86_64__
+  _mm_lfence();
+#endif
+
+#ifdef __aarch64__
+  __isb(0xF);
+#endif
+
+  asm volatile("" ::: "memory");
+
+  return acc;
+}
+
+int main(int argc, const char *argv[]) {
+  // Check that the correct number of arguments has been passed.
+  //
+  // Expected usage:
+  //
+  //   ./main <number_of_pedestrians> <number_of_iterations> [csv_output_path]
+  //
+  // The optional third argument enables CSV dumping for visualization.
+  if (argc < 3) {
+    printf("Usage: ./main <number_of_pedestrians> <number_of_iterations> [csv_output_path]\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'\n", argv[1]);
+    return 1;
+  }
+  if (n_iter <= 0) {
+    printf("Error: <number_of_iterations> must be a positive integer, got '%s'\n", argv[2]);
+    return 1;
+  }
+
+  // If provided, dump_path is used to write simulation states to CSV.
+  // In benchmark mode, this path is intentionally ignored to avoid measuring
+  // file I/O overhead.
+  const char *dump_path = (argc >= 4) ? argv[3] : nullptr;
+
+#if BENCHMARK
+
+#if ROOFLINE
+  volatile long long trash_value;
+#endif
+
+#ifdef __aarch64__
+  // Initialize ARM performance counters.
+  kperf_init();
+  struct performance_counters startperf, endperf;
+#endif
+
+  // Warmup runs, to bring caches, predictors, frequency to steady state. These are not measured.
+  // Each run gets a fresh initialized state because run_simulation mutates the
+  // pedestrians and board state.
+  for (int i = 0; i < WARMUP_ITERS; i++) {
+    Pedestrian *P = init_pedestrians(n_p);
+    Board board = init_board_without_obstacles();
+    run_simulation(P, board, n_p, n_iter, nullptr);
+    delete[] P;
+  }
+
+  // Calibration: finds num_runs such that one batch runs for at least
+  // CYCLES_REQUIRED cycles. A long batch keeps makes the timer and other
+  // noise overhead, negligible relative to the kernel execution cost we're
+  // actually trying to measure.
+  //
+  // Skipped when ROOFLINE is enabled: cold-cache measurements need
+  // exactly one timed simulation per batch, otherwise only the first
+  // call starts cold and subsequent calls measure warm-cache behavior.
+  // We can use large enough inputs so a single call already exceeds CYCLES_REQUIRED,
+  // when measuring for roofline plot.
+  // Note: we implemented the calibration exactly as it is done in mvm.c from hw1, we also
+  // used the same values for NUM_RUNS, etc
+  int num_runs = NUM_RUNS;
+#if !ROOFLINE
+  while (num_runs < MAX_NUM_RUNS) {
+    Pedestrian *P = init_pedestrians(n_p);
+    Board board = init_board_without_obstacles();
+
+    float cycles;
+
+#ifdef __x86_64__
+    myInt64 start = start_tsc();
+#endif
+
+#ifdef __aarch64__
+    startperf = kperf_get_counters();
+#endif
+
+    for (int j = 0; j < num_runs; j++) {
+      run_simulation(P, board, n_p, n_iter, nullptr);
+    }
+
+#ifdef __x86_64__
+    cycles = (float)stop_tsc(start);
+#endif
+
+#ifdef __aarch64__
+    endperf = kperf_get_counters();
+    cycles = (float)(endperf.cycles - startperf.cycles);
+#endif
+
+    delete[] P;
+
+    if (cycles >= CYCLES_REQUIRED)
+      break;
+    num_runs *= 2;
+  }
+#endif
+
+  // Store the cycle count of each measured repetition.
+  // The median is then used as the final timing result.
+  std::vector<float> cycles_v;
+  cycles_v.reserve(NUM_REPS);
+
+  // Measurement phase
+  for (int i = 0; i < NUM_REPS; i++) {
+    Pedestrian *P = init_pedestrians(n_p);
+    Board board = init_board_without_obstacles();
+
+#if ROOFLINE
+    trash_value = flush_cache();
+#endif
+
+#ifdef __x86_64__
+    myInt64 start = start_tsc();
+#endif
+#ifdef __aarch64__
+    startperf = kperf_get_counters();
+#endif
+
+    for (int j = 0; j < num_runs; j++) {
+      // Disable CSV dumping during benchmarking to avoid measuring file I/O
+      // overhead. Initialization is also excluded from the measured region.
+
+      run_simulation(P, board, n_p, n_iter, nullptr);
+    }
+
+#ifdef __x86_64__
+    float cycles = (float)stop_tsc(start) / num_runs;
+#endif
+#ifdef __aarch64__
+    endperf = kperf_get_counters();
+    float cycles = (float)(endperf.cycles - startperf.cycles) / num_runs;
+#endif
+
+    delete[] P;
+    cycles_v.push_back(cycles);
+  }
+
+  std::sort(cycles_v.begin(), cycles_v.end());
+  // Use the median cycle count to reduce the effect of noisy measurements.
+  float T = cycles_v[NUM_REPS / 2];
+
+#else
+
+  // Normal execution mode. If dump_path is provided, the simulation writes a CSV
+  // file that can later be used to generate a video of the simulation.
+  Pedestrian *P = init_pedestrians(n_p);
+  Board board = init_board_without_obstacles();
+  run_simulation(P, board, n_p, n_iter, dump_path);
+
+  delete[] P;
+
+#endif
+
+#if BENCHMARK
+
+  // Flop-count per simulation iteration.
+  //
+  // The total work W(n) is estimated from the number of additions,
+  // multiplications, divisions, square roots, and exponentials performed by the
+  // simulation kernels.
+
+  long long add_count_per_iter = n_p * (2 + (n_p - 1) * 14 + N_OBSTACLES * 10 + 2) + n_p * 7;
+  long long mult_count_per_iter = n_p * (4 + (n_p - 1) * 24 + N_OBSTACLES * 13) + n_p * 16;
+  long long div_count_per_iter = n_p * ((n_p - 1) * 1 + N_OBSTACLES * 1) + n_p * 2;
+  long long sqrt_count_per_iter = n_p * ((n_p - 1) * 3 + N_OBSTACLES * 1) + n_p * 3;
+  long long exp_count_per_iter = n_p * ((n_p - 1) * 1 + N_OBSTACLES * 1);
+
+  float 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;
+
+  // Estimated memory traffic per iteration.
+  // TODO: refine this model to account for actual memory accesses, cache reuse,
+  // and the full interaction loop.
+
+  long long Q_per_iter = sizeof(Pedestrian) * n_p + sizeof(Obstacle) * N_OBSTACLES;
+
+  // Operational intensity:
+  //
+  //   I(n) = W(n) / Q(n)
+  //
+  // where W is estimated work in flops and Q is estimated memory traffic in
+  // bytes.
+
+  float I = W_per_iter / Q_per_iter;
+  float W = W_per_iter * n_iter;
+
+  // Performance in estimated flops per cycle.
+  float P = W / T;
+
+  printf("W(n) per iteration: %.2f flops\n", W_per_iter);
+  printf("Q(n) per iteration: %lld bytes\n", Q_per_iter);
+  printf("I(n): %.2f flops/byte\n", I);
+  printf("W(n): %.2f flops\n", W);
+  printf("T(n): %.2f cycles\n", T);
+  printf("P(n): %.2f flops/cycle\n", P);
+
+#endif
+
+  return 0;
+}

main.cpp

@@ -242,10 +242,10 @@ int main(int argc, const char *argv[]) {
   // multiplications, divisions, square roots, and exponentials performed by the
   // simulation kernels.
 
-  long long add_count_per_iter = n_p * (4 + (n_p - 1) * 14 + N_OBSTACLES * 12) + n_p * 7;
+  long long add_count_per_iter = n_p * (2 + (n_p - 1) * 14 + N_OBSTACLES * 10 + 2) + n_p * 7;
-  long long mult_count_per_iter = n_p * (4 + (n_p - 1) * 23 + N_OBSTACLES * 18) + n_p * 16;
+  long long mult_count_per_iter = n_p * (4 + (n_p - 1) * 24 + N_OBSTACLES * 13) + n_p * 16;
   long long div_count_per_iter = n_p * ((n_p - 1) * 1 + N_OBSTACLES * 1) + n_p * 2;
-  long long sqrt_count_per_iter = n_p * ((n_p - 1) * 4 + N_OBSTACLES * 2) + n_p * 3;
+  long long sqrt_count_per_iter = n_p * ((n_p - 1) * 3 + N_OBSTACLES * 1) + n_p * 3;
   long long exp_count_per_iter = n_p * ((n_p - 1) * 1 + 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;

scripts/verify.py

@@ -3,33 +3,37 @@ import csv
 import math
 
 if len(sys.argv) < 3:
-    sys.exit("Use it like: python verify.py reference.csv generated.csv [tol]")
+  sys.exit('Use it like: python verify.py reference.csv generated.csv [tol]')
-# We can use this later in case that we wanna compare something with a tolerance
-tol = float(sys.argv[3]) if len(sys.argv) > 3 else 0.00001
 
-with open(sys.argv[1], newline="") as f:
+TOL = float(sys.argv[3]) if len(sys.argv) > 3 else 1
+DIFF_COL_RATIO_THRESHOLD = 0.025
+
+with open(sys.argv[1], newline='') as f:
   ref = list(csv.reader(f))
-with open(sys.argv[2], newline="") as f:
+with open(sys.argv[2], newline='') as f:
   gen = list(csv.reader(f))
 
 if len(ref) != len(gen):
-    sys.exit(f"row count differs: ref={len(ref)} got={len(gen)}")
+  sys.exit(f'row count differs: ref={len(ref)} got={len(gen)}')
 
 header = ref[0]
+diff_col_count = 0
 for i in range(len(ref)):
   if len(ref[i]) != len(gen[i]):
-        sys.exit(f"line {i+1}: column count differs")
+    sys.exit(f'line {i+1}: column count differs')
   for j in range(len(ref[i])):
     a = ref[i][j].strip()
     b = gen[i][j].strip()
     if a == b:
       continue
-        # in case that we wanna use tolerance
+    if TOL > 0:
-        if tol > 0:
       try:
-                if math.isclose(float(a), float(b), abs_tol=tol):
+        if math.isclose(float(a), float(b), abs_tol=TOL):
           continue
+        else:
+          diff_col_count += 1
       except ValueError:
         pass
-        col = header[j] if j < len(header) else f"col{j}"
+
-        sys.exit(f"mismatch at line {i+1}, col {j} ('{col}'): ref='{a}' got='{b}'")
+if diff_col_count / len(ref) > DIFF_COL_RATIO_THRESHOLD:
+  sys.exit(f'Error Threshold Exceeded!')

sizes.cpp

@@ -1,18 +0,0 @@
-#include "./headers/pedestrian.hpp"
-#include "./headers/board.hpp"
-#include <iostream>
-#include <math.h>
-
-int main() {
-  std::cout << sizeof(Color);
-  std::cout << "\n";
-  std::cout << sizeof(Vector2d);
-  std::cout << "\n";
-  std::cout << sizeof(Pedestrian);
-  std::cout << "\n";
-  std::cout << sizeof(Obstacle);
-  std::cout << "\n";
-  std::cout << sizeof(Board);
-  std::cout << "\n";
-  return 0;
-}

variants/variant_baseline.cpp

@@ -35,9 +35,21 @@ using namespace std;
 
 // forces.cpp
 
-// Flop count: 2 ADDs, 5 MULTs, 1 SQRT
+// Flop count: 2 ADDs, 6 MULTs
-double w(Vector2d e, Vector2d f) {
+double w_ped(Vector2d e, Vector2d f) {
-  if (e.dot(f) >= f.norm() * cos_phi)
+  double f_norm_sq = f.dot(f);
+  double rhs = f_norm_sq * cos_phi_sq;
+  double e_dot_f = e.dot(f);
+  double e_dot_f_sq = e_dot_f * e_dot_f;
+  if (e_dot_f >= 0.0 || e_dot_f_sq <= rhs)
+    return dt_neg;
+  else
+    return c_neg_times_dt;
+}
+
+double w_obstacle(Vector2d e, Vector2d f) {
+  double e_dot_f = (f.y == 0) ? 0 : (f.y < 0 ? -e.y : e.y);
+  if (e_dot_f >= cos_phi)
     return dt_neg;
   else
     return c_neg_times_dt;
@@ -57,7 +69,7 @@ Vector2d acceleration_term_dt(const Pedestrian &alpha) {
   return f_alpha_0_dt;
 }
 
-// Flop count: 12 ADDs, 23 MULTs, 1 DIV, 4 SQRTs, 1 EXP
+// Flop count: 12 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP
 Vector2d pedestrian_repulsion_dt(const Pedestrian &alpha, const Pedestrian &beta) {
 
   // # 2 ADDs
@@ -129,16 +141,16 @@ Vector2d pedestrian_repulsion_dt(const Pedestrian &alpha, const Pedestrian &beta
   //
   Vector2d f_alpha_beta_neg = r_ab * scalar_mult_r_ab + r_ab2 * scalar_mult_r_ab2;
 
-  // # 2 ADDs, 7 MULTs, 1 SQRT
+  // # 2 ADDs, 8 MULTs
   //
   // $\vec{F}_{\alpha\beta}(\vec{e}_{\alpha}, \vec{r}_{\alpha \beta}).dt = w(\vec{e}_{\alpha}, -\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta})).\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta}).dt$
   //
-  Vector2d f_alpha_beta_dt = f_alpha_beta_neg * w(alpha.e, f_alpha_beta_neg);
+  Vector2d f_alpha_beta_dt = f_alpha_beta_neg * w_ped(alpha.e, f_alpha_beta_neg);
 
   return f_alpha_beta_dt;
 }
 
-// Flop count: 10 ADDs, 18 MULTs, 1 DIVs, 2 SQRTs, 1 EXP
+// Flop count: 8 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP
 Vector2d obstacle_repulsion_dt(const Pedestrian &alpha, const Obstacle &B) {
 
   // # 7 ADDs, 5 MULTs
@@ -162,11 +174,11 @@ Vector2d obstacle_repulsion_dt(const Pedestrian &alpha, const Obstacle &B) {
   //
   Vector2d f_alpha_B_neg = r_aB * scalar_mult;
 
-  // # 2 ADDs, 7 MULTs, 1 SQRT
+  // # 2 MULTs
   //
   // $\vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha B}).dt = w(\vec{e}_{\alpha}, -\vec{F}_{\alpha B}(\vec{r}_{\alpha B})).\vec{F}_{\alpha B}(\vec{r}_{\alpha B}).dt$
   //
-  Vector2d f_alpha_B_dt = f_alpha_B_neg * w(alpha.e, f_alpha_B_neg);
+  Vector2d f_alpha_B_dt = f_alpha_B_neg * w_obstacle(alpha.e, f_alpha_B_neg);
 
   return f_alpha_B_dt;
 }
@@ -280,8 +292,7 @@ int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, c
       //
       Vector2d f_alpha_dt = acceleration_term_dt(p);
 
-      // # (N-1) * (12 ADDs, 23 MULTs, 1 DIV, 4 SQRTs, 1 EXP)
+      // # (N-1) * (14 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP)
-      //   + (N-1) * 2 ADDs
       //
       // $\sum_{\beta} \vec{F}_{\alpha\beta}(\vec{e}_{\alpha},\vec{r}_{\alpha} - \vec{r}_{\beta})$
       //
@@ -292,8 +303,7 @@ int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, c
         f_alpha_dt += pedestrian_repulsion_dt(p, pedestrians[j]);
       }
 
-      // # board.n_obstacles * (10 ADDs, 18 MULTs, 1 DIV, 2 SQRTs, 1 EXP)
+      // # board.n_obstacles * (10 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP)
-      //   + board.n_obstacles * 2 ADDs
       //
       // $\sum_{B} \vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{B})$
       //

variants/variant_float_baseline.cpp

@@ -0,0 +1,409 @@
+// This file is merged from the forces.cpp and social_forces_baseline.cpp
+// The goal is that we from now on only have one file as our baseline
+// This should make it more easy in the future to implement optimizations
+// Our naming convention from now on is: variant_<variant_name>.cpp (Please make sure to keep that consistent, otherwise our Makefile wont work anymore)
+
+#include "../headers/board_float.hpp"
+#include "../headers/constants_float.hpp"
+#include "../headers/pedestrian_float.hpp"
+#include "../headers/vector2d_float.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <math.h> // keep — forces code uses unqualified sqrt/exp
+#include <random>
+#include <tuple>
+#include <vector>
+
+#define SEED 42
+#define BOARD_WIDTH 10.0f
+#define BOARD_LENGTH 50.0f
+#define WALL_EPS 1e-3f
+#define BOARD_WIDTH_MINUS_WALL_EPS 9.999f
+
+#ifndef RANDOM_SEED
+#define RANDOM_SEED 0
+#endif
+
+#ifndef BENCHMARK
+#define BENCHMARK 0
+#endif
+
+using namespace std;
+
+// forces.cpp
+
+// Flop count: 2 ADDs, 6 MULTs
+float w_ped(Vector2d e, Vector2d f) {
+  float f_norm_sq = f.dot(f);
+  float rhs = f_norm_sq * cos_phi_sq;
+  float e_dot_f = e.dot(f);
+  float e_dot_f_sq = e_dot_f * e_dot_f;
+  if (e_dot_f >= 0.0 || e_dot_f_sq <= rhs)
+    return dt_neg;
+  else
+    return c_neg_times_dt;
+}
+
+float w_obstacle(Vector2d e, Vector2d f) {
+  float e_dot_f = (f.y == 0) ? 0 : (f.y < 0 ? -e.y : e.y);
+  if (e_dot_f >= cos_phi)
+    return dt_neg;
+  else
+    return c_neg_times_dt;
+}
+
+// Flop count: 2 ADDs, 4 MULTs
+Vector2d acceleration_term_dt(const Pedestrian &alpha) {
+
+  // # 2 MULTs
+  Vector2d desired_velocity = alpha.e * alpha.desired_v0;
+
+  Vector2d actual_velocity = alpha.curr_vel;
+
+  // # 2 ADDs, 2 MULTs
+  Vector2d f_alpha_0_dt = (desired_velocity - actual_velocity) * dt_times_tau_alpha_inv;
+
+  return f_alpha_0_dt;
+}
+
+// Flop count: 12 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP
+Vector2d pedestrian_repulsion_dt(const Pedestrian &alpha, const Pedestrian &beta) {
+
+  // # 2 ADDs
+  // $\vec{r}_{\alpha\beta}$
+  Vector2d r_ab = alpha.curr_pos - beta.curr_pos;
+
+  // # 1 ADD, 2 MULTs, 1 SQRT
+  // $\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert$
+  float r_ab_norm = r_ab.norm();
+  r_ab_norm = r_ab_norm > EPS ? r_ab_norm : EPS; // Avoid division by zero
+
+  // $v_\beta \Delta t \vec{e}_\beta$
+  Vector2d s_b_e = beta.s_b_e;
+
+  // # 2 ADDs
+  // $\vec{r}_{\alpha\beta} - \vec{v}_{\beta}\varDelta_{t}\vec{e}_{\beta}$
+  Vector2d r_ab2 = r_ab - s_b_e;
+
+  // # 1 ADD, 2 MULTs, 1 SQRT
+  // $\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert$
+  float r_ab2_norm = r_ab2.norm();
+  r_ab2_norm = r_ab2_norm > EPS ? r_ab2_norm : EPS; // Avoid division by zero
+
+  // # 1 ADD
+  // $\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert$
+  float r_ab_norm_sum = r_ab_norm + r_ab2_norm;
+
+  // # 1 ADD, 1 MULTs, 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}$
+  //
+  float two_b_sq = r_ab_norm_sum * r_ab_norm_sum - beta.s_b_sq;
+  two_b_sq = two_b_sq > EPS_SQ ? two_b_sq : EPS_SQ; // Avoid sqrt of negative and avoid division by zero
+  float two_b = sqrt(two_b_sq);
+
+  // # 5 MULTs, 1 DIV, 1 EXP
+  //
+  // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}$
+  //
+  float scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV));
+
+  // # 1 MULT
+  //
+  // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}$
+  //
+  float scalar_mult_r_ab = scalar_mult_common * r_ab2_norm;
+
+  // # 1 MULT
+  //
+  // $\frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}$
+  //
+  float scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm;
+
+  // # 2 ADDs, 4 MULTs
+  //
+  // $\vec{-f}_{\alpha \beta}(\vec{r}_{\alpha \beta}) = \frac{-V^0_{\alpha \beta}}{2\sigma}e^{-b / \sigma} 2\grad_{\vec{r}_{\alpha \beta}}b $
+  //
+  //
+  //          $= \frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}\,\vec{r}_{\alpha \beta} + \frac{-V^0_{\alpha\beta}}{2\sigma}\,e^{-b/\sigma}\,\frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} -
+  //          v_\beta \Delta t \vec{e}_\beta \Arrowvert}\,(\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta)$
+  //
+  //
+  //
+  // $2\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})}{2b}$
+  //
+  //
+  //       $= \frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert}\,\vec{r}_{\alpha \beta} + \frac{\Arrowvert \vec{r}_{\alpha \beta} \Arrowvert + \Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}{2b\,\Arrowvert \vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta \Arrowvert}\,(\vec{r}_{\alpha \beta} - v_\beta \Delta t \vec{e}_\beta)$
+  //
+  //
+  Vector2d f_alpha_beta_neg = r_ab * scalar_mult_r_ab + r_ab2 * scalar_mult_r_ab2;
+
+  // # 2 ADDs, 8 MULTs
+  //
+  // $\vec{F}_{\alpha\beta}(\vec{e}_{\alpha}, \vec{r}_{\alpha \beta}).dt = w(\vec{e}_{\alpha}, -\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta})).\vec{f}_{\alpha\beta}(\vec{r}_{\alpha \beta}).dt$
+  //
+  Vector2d f_alpha_beta_dt = f_alpha_beta_neg * w_ped(alpha.e, f_alpha_beta_neg);
+
+  return f_alpha_beta_dt;
+}
+
+// Flop count: 8 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP
+Vector2d obstacle_repulsion_dt(const Pedestrian &alpha, const Obstacle &B) {
+
+  // # 7 ADDs, 5 MULTs
+  // $\vec{r}_{\alpha B}$
+  Vector2d r_aB = alpha.curr_pos - B.find_nearest_obstacle_point(alpha.curr_pos);
+
+  // # 1 ADD, 2 MULTs, 1 SQRT
+  // $\Arrowvert \vec{r}_{\alpha B} \Arrowvert$
+  float r_aB_norm = r_aB.norm();
+  r_aB_norm = r_aB_norm > EPS ? r_aB_norm : EPS; // Avoid division by zero
+
+  // # 2 MULTs, 1 DIV, 1 EXP
+  //
+  // $\frac{-U^0_{\alpha B}}{R} e^{- \Arrowvert \vec{r}_{\alpha B} \Arrowvert / R}$
+  //
+  float scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm;
+
+  // # 2 MULTs
+  //
+  // $-\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_alpha_B_neg = r_aB * scalar_mult;
+
+  // # 2 MULTs
+  //
+  // $\vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha B}).dt = w(\vec{e}_{\alpha}, -\vec{F}_{\alpha B}(\vec{r}_{\alpha B})).\vec{F}_{\alpha B}(\vec{r}_{\alpha B}).dt$
+  //
+  Vector2d f_alpha_B_dt = f_alpha_B_neg * w_obstacle(alpha.e, f_alpha_B_neg);
+
+  return f_alpha_B_dt;
+}
+
+// social_force_baseline.cpp
+
+Pedestrian *init_pedestrians(int n_p) {
+#if RANDOM_SEED
+  random_device rd;  // Platform-provided random seed source.
+  mt19937 gen(rd()); // Non-reproducible initialization.
+#else
+  mt19937 gen(SEED); // Reproducible initialization.
+#endif
+
+  normal_distribution<float> velocity_gen(v0_mean, v0_std_dev);
+  uniform_real_distribution<float> y_gen(0.3, BOARD_WIDTH - 0.3);
+  uniform_real_distribution<float> x_gen_left(2.0, BOARD_LENGTH / 2.0);
+  uniform_real_distribution<float> x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0);
+
+  Pedestrian *P = new Pedestrian[n_p];
+
+  int n_p_half = n_p / 2;
+
+  // Group of pedestrians walking from left to right.
+  for (int i = 0; i < n_p_half; i++) {
+    float x = x_gen_left(gen);
+    float y = y_gen(gen);
+
+    Vector2d start_pos = Vector2d(x, y);
+    Vector2d destination = Vector2d(BOARD_LENGTH + 30.0, y);
+
+    Vector2d init_e = destination - start_pos;
+    float init_e_norm = init_e.norm();
+    Vector2d init_e_unit = init_e / init_e_norm;
+
+    // Reject near-zero desired speeds to avoid pedestrians that barely move,
+    // and reject negative speeds because they would reverse the intended
+    // walking direction.
+    float desired_v0;
+    do {
+      desired_v0 = velocity_gen(gen);
+    } while (desired_v0 <= 0.1);
+
+    float v_max = 1.3 * desired_v0;
+
+    P[i] = Pedestrian(WHITE, Vector2d(), start_pos, init_e_unit, v_max, desired_v0, destination);
+  }
+
+  // Group of pedestrians walking from right to left.
+  for (int i = n_p_half; i < n_p; i++) {
+    float x = x_gen_right(gen);
+    float y = y_gen(gen);
+
+    Vector2d start_pos = Vector2d(x, y);
+    Vector2d destination = Vector2d(-30.0, y);
+
+    Vector2d init_e = destination - start_pos;
+    float init_e_norm = init_e.norm();
+    Vector2d init_e_unit = init_e / init_e_norm;
+
+    // Reject near-zero desired speeds to avoid pedestrians that barely move,
+    // and reject negative speeds because they would reverse the intended
+    // walking direction.
+    float desired_v0;
+    do {
+      desired_v0 = velocity_gen(gen);
+    } while (desired_v0 <= 0.1);
+
+    float v_max = 1.3 * desired_v0;
+
+    P[i] = Pedestrian(BLACK, Vector2d(), start_pos, init_e_unit, v_max, desired_v0, destination);
+  }
+
+  return P;
+}
+
+Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); }
+
+int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) {
+
+#if BENCHMARK
+  // Avoid compiler unused parameter warning
+  (void)dump_path;
+#endif
+
+#if !BENCHMARK
+  const int dump_interval = 1;
+
+  // Each sample stores: time step, pedestrian id, x, y, e.x, e.y, color.
+  vector<tuple<int, int, float, float, float, float, int>> samples;
+
+  if (dump_path) {
+    samples.reserve(static_cast<size_t>(n_iter / dump_interval + 1) * n_p);
+
+    // Store the initial state before the first simulation step.
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+      samples.emplace_back(0, i, p.curr_pos.x, p.curr_pos.y, p.e.x, p.e.y, p.color);
+    }
+  }
+#endif
+
+  for (int t = 0; t < n_iter; t++) {
+    // Pass 1: Compute all forces from the current state.
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+
+      // # 2 ADDs, 4 MULTs
+      //
+      // $\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)$
+      //
+      Vector2d f_alpha_dt = acceleration_term_dt(p);
+
+      // # (N-1) * (14 ADDs, 24 MULTs, 1 DIV, 3 SQRTs, 1 EXP)
+      //
+      // $\sum_{\beta} \vec{F}_{\alpha\beta}(\vec{e}_{\alpha},\vec{r}_{\alpha} - \vec{r}_{\beta})$
+      //
+      for (int j = 0; j < n_p; j++) {
+        if (i == j)
+          continue;
+
+        f_alpha_dt += pedestrian_repulsion_dt(p, pedestrians[j]);
+      }
+
+      // # board.n_obstacles * (10 ADDs, 13 MULTs, 1 DIV, 1 SQRT, 1 EXP)
+      //
+      // $\sum_{B} \vec{F}_{\alpha B}(\vec{e}_{\alpha}, \vec{r}_{\alpha} - \vec{r}_{B})$
+      //
+      for (int k = 0; k < board.n_obstacles; k++) {
+        f_alpha_dt += obstacle_repulsion_dt(p, board.obstacles[k]);
+      }
+
+      // # 2 ADDs
+      //
+      //
+      // $\frac{d\vec{w}_\alpha}{dt} = \vec{F}_{\alpha}(t)$
+      //
+      p.w += f_alpha_dt;
+    }
+
+    // Pass 2: Update pedestrian states after all forces have been computed.
+    // # (7 ADDs, 16 MULTs, 2 DIVs, 3 SQRTs) * n_p
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+
+      // # 1 ADD, 2 MULTs, 1 SQRT
+      float w_norm = p.w.norm();
+
+      float w_norm2 = w_norm > EPS ? w_norm : EPS; // Avoid division by zero.
+
+      // # 1 DIV
+      float correction_factor = std::min(1.0f, p.v_max / w_norm2);
+
+      // # 2 MULTs
+      p.curr_vel = p.w * correction_factor;
+
+      // # 2 ADDs, 2 MULTs
+      //
+      // $\vec{r}_i^{t+1} = \vec{r}_i^t + \vec{v}_i^t dt$
+      //
+      p.curr_pos = p.curr_pos + p.curr_vel * dt;
+
+      // When a pedestrian reaches the end of the corridor, teleport them back
+      // to the opposite entrance while keeping the same y-coordinate. This
+      // maintains a constant pedestrian flow and makes lane formation easier to
+      // observe.
+      if (p.color == WHITE && p.curr_pos.x > (BOARD_LENGTH - 1.5)) {
+        p.curr_pos.x = 1.0;
+        p.destination.y = p.curr_pos.y;
+        p.w = Vector2d(p.desired_v0, 0.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.destination.y = p.curr_pos.y;
+        p.w = Vector2d(-p.desired_v0, 0.0);
+        p.curr_vel = p.w;
+      }
+
+      // Keep pedestrians inside the board boundaries.
+      p.curr_pos.y = clamp(p.curr_pos.y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS);
+
+      // # 2 ADDs
+      Vector2d e = p.destination - p.curr_pos;
+
+      // # 1 ADD, 2 MULTs, 1 DIV, 1 SQRT
+      float e_norm = e.norm();
+      float e_norm2 = e_norm > EPS ? e_norm : EPS; // Avoid division by zero.
+      float e_norm_inv = 1.0 / e_norm2;
+
+      // # 2 MULTs
+      p.e = e * e_norm_inv;
+
+      // # 1 ADD, 2 MULTs
+      float velocity_norm_sq = p.curr_vel.dot(p.curr_vel);
+
+      // # 1 SQRT
+      float velocity_norm = std::sqrt(velocity_norm_sq);
+
+      // # 1 MULT
+      p.s_b_sq = velocity_norm_sq * delta_t_sq;
+
+      // # 3 MULTs
+      p.s_b_e = p.e * (velocity_norm * delta_t);
+
+#if !BENCHMARK
+      // Store the current pedestrian state for visualization.
+      if (dump_path && (t % dump_interval == 0)) {
+        samples.emplace_back(t + 1, i, p.curr_pos.x, p.curr_pos.y, p.e.x, p.e.y, p.color);
+      }
+#endif
+    }
+  }
+
+#if !BENCHMARK
+  // Export recorded simulation states to CSV.
+  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';
+    }
+  }
+#endif
+
+  return 0;
+}

variants/variant_float_inlining.cpp

@@ -0,0 +1,296 @@
+// # Changes: Replaced «double» with «float» in variable declarations
+// # Predecessor: variant_float_baseline.cpp
+
+#include "../headers/board_reduced_float.hpp"
+#include "../headers/constants_float.hpp"
+#include "../headers/pedestrian_reduced_float.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <math.h>
+#include <random>
+#include <tuple>
+#include <vector>
+
+#define SEED 42
+#define BOARD_WIDTH 10.0f
+#define BOARD_LENGTH 50.0f
+#define WALL_EPS 1e-3f
+#define BOARD_WIDTH_MINUS_WALL_EPS 9.999f
+
+#ifndef RANDOM_SEED
+#define RANDOM_SEED 0
+#endif
+
+#ifndef BENCHMARK
+#define BENCHMARK 0
+#endif
+
+using namespace std;
+
+Pedestrian *init_pedestrians(int n_p) {
+#if RANDOM_SEED
+  random_device rd;
+  mt19937 gen(rd());
+#else
+  mt19937 gen(SEED);
+#endif
+
+  normal_distribution<float> velocity_gen(v0_mean, v0_std_dev);
+  uniform_real_distribution<float> y_gen(0.3, BOARD_WIDTH - 0.3);
+  uniform_real_distribution<float> x_gen_left(2.0, BOARD_LENGTH / 2.0);
+  uniform_real_distribution<float> x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0);
+
+  Pedestrian *P = new Pedestrian[n_p];
+
+  int n_p_half = n_p / 2;
+
+  for (int i = 0; i < n_p_half; i++) {
+    float x = x_gen_left(gen);
+    float y = y_gen(gen);
+
+    float start_pos_x = x;
+    float start_pos_y = y;
+
+    float destination_x = BOARD_LENGTH + 30.0;
+    float destination_y = y;
+
+    float init_e_x = destination_x - start_pos_x;
+    float init_e_y = destination_y - start_pos_y;
+    float init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y);
+    float init_e_unit_x = init_e_x / init_e_norm;
+    float init_e_unit_y = init_e_y / init_e_norm;
+
+    float desired_v0;
+    do {
+      desired_v0 = velocity_gen(gen);
+    } while (desired_v0 <= 0.1);
+
+    float v_max = 1.3 * desired_v0;
+
+    P[i] = Pedestrian(0, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y);
+  }
+
+  for (int i = n_p_half; i < n_p; i++) {
+    float x = x_gen_right(gen);
+    float y = y_gen(gen);
+
+    float start_pos_x = x;
+    float start_pos_y = y;
+    float destination_x = -30;
+    float destination_y = y;
+
+    float init_e_x = destination_x - start_pos_x;
+    float init_e_y = destination_y - start_pos_y;
+    float init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y);
+    float init_e_unit_x = init_e_x / init_e_norm;
+    float init_e_unit_y = init_e_y / init_e_norm;
+
+    float desired_v0;
+    do {
+      desired_v0 = velocity_gen(gen);
+    } while (desired_v0 <= 0.1);
+
+    float v_max = 1.3 * desired_v0;
+
+    P[i] = Pedestrian(1, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y);
+  }
+
+  return P;
+}
+
+Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); }
+
+int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) {
+
+#if BENCHMARK
+  (void)dump_path;
+#endif
+
+#if !BENCHMARK
+  const int dump_interval = 1;
+
+  vector<tuple<int, int, float, float, float, float, int>> samples;
+
+  if (dump_path) {
+    samples.reserve(static_cast<size_t>(n_iter / dump_interval + 1) * n_p);
+
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+      samples.emplace_back(0, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color);
+    }
+  }
+#endif
+
+  for (int t = 0; t < n_iter; t++) {
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+
+      float desired_velocity_x = p.e_x * p.desired_v0;
+      float desired_velocity_y = p.e_y * p.desired_v0;
+      float f_alpha_dt_x = (desired_velocity_x - p.curr_vel_x) * dt_times_tau_alpha_inv;
+      float f_alpha_dt_y = (desired_velocity_y - p.curr_vel_y) * dt_times_tau_alpha_inv;
+
+      for (int j = 0; j < n_p; j++) {
+        if (i == j)
+          continue;
+
+        float r_ab_x = p.curr_pos_x - pedestrians[j].curr_pos_x;
+        float r_ab_y = p.curr_pos_y - pedestrians[j].curr_pos_y;
+
+        float r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y);
+        r_ab_norm = r_ab_norm > EPS ? r_ab_norm : EPS;
+
+        float r_ab2_x = r_ab_x - pedestrians[j].s_b_e_x;
+        float r_ab2_y = r_ab_y - pedestrians[j].s_b_e_y;
+
+        float r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y);
+        r_ab2_norm = r_ab2_norm > EPS ? r_ab2_norm : EPS;
+
+        float r_ab_norm_sum = r_ab_norm + r_ab2_norm;
+
+        float two_b_sq = r_ab_norm_sum * r_ab_norm_sum - pedestrians[j].s_b_sq;
+        two_b_sq = two_b_sq > EPS_SQ ? two_b_sq : EPS_SQ;
+        float two_b = sqrt(two_b_sq);
+        float scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV));
+        float scalar_mult_r_ab = scalar_mult_common * r_ab2_norm;
+        float scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm;
+
+        float f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2;
+        float f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2;
+
+        float w_val;
+
+        float f_norm_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y;
+        float rhs = f_norm_sq * cos_phi_sq;
+        float e_dot_f = p.e_x * f_alpha_beta_neg_x + p.e_y * f_alpha_beta_neg_y;
+        float e_dot_f_sq = e_dot_f * e_dot_f;
+        if (e_dot_f >= 0.0 || e_dot_f_sq <= rhs)
+          w_val = dt_neg;
+        else
+          w_val = c_neg_times_dt;
+
+        f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val;
+        f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val;
+      }
+
+      for (int k = 0; k < board.n_obstacles; k++) {
+
+        float v1_x = p.curr_pos_x - board.obstacles[k].start_x;
+        float v1_y = p.curr_pos_y - board.obstacles[k].start_y;
+        float v2_x = board.obstacles[k].obstacle_v_x;
+        float v2_y = board.obstacles[k].obstacle_v_y;
+
+        float projection = (v1_x * v2_x + v1_y * v2_y) * board.obstacles[k].length_sq_inv;
+
+        if (projection < 0.0)
+          projection = 0.0;
+        else if (projection > 1.0)
+          projection = 1.0;
+
+        float r_aB_x = p.curr_pos_x - (board.obstacles[k].start_x + board.obstacles[k].obstacle_v_x * projection);
+        float r_aB_y = p.curr_pos_y - (board.obstacles[k].start_y + board.obstacles[k].obstacle_v_y * projection);
+
+        float r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y);
+        r_aB_norm = r_aB_norm > EPS ? r_aB_norm : EPS;
+
+        float scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm;
+        float f_alpha_B_neg_x = r_aB_x * scalar_mult;
+        float f_alpha_B_neg_y = r_aB_y * scalar_mult;
+
+        float w_val;
+
+        float e_dot_f = (f_alpha_B_neg_y == 0) ? 0 : (f_alpha_B_neg_y < 0 ? -p.e_y : p.e_y);
+
+        if (e_dot_f >= cos_phi)
+          w_val = dt_neg;
+        else
+          w_val = c_neg_times_dt;
+
+        float f_alpha_B_dt_x = f_alpha_B_neg_x * w_val;
+        float f_alpha_B_dt_y = f_alpha_B_neg_y * w_val;
+
+        f_alpha_dt_x += f_alpha_B_dt_x;
+        f_alpha_dt_y += f_alpha_B_dt_y;
+      }
+
+      p.w_x += f_alpha_dt_x;
+      p.w_y += f_alpha_dt_y;
+    }
+
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+
+      float w_norm = sqrt(p.w_x * p.w_x + p.w_y * p.w_y);
+
+      float w_norm2 = w_norm > EPS ? w_norm : EPS;
+
+      float correction_factor = min(1.0f, p.v_max / w_norm2);
+
+      p.curr_vel_x = p.w_x * correction_factor;
+      p.curr_vel_y = p.w_y * correction_factor;
+
+      p.curr_pos_x = p.curr_pos_x + p.curr_vel_x * dt;
+      p.curr_pos_y = p.curr_pos_y + p.curr_vel_y * dt;
+
+      if (p.color == 0 && p.curr_pos_x > (BOARD_LENGTH - 1.5)) {
+        p.curr_pos_x = 1.0;
+        p.destination_y = p.curr_pos_y;
+        p.w_x = p.desired_v0;
+        p.w_y = 0.0;
+        p.curr_vel_x = p.w_x;
+        p.curr_vel_y = p.w_y;
+      }
+      if (p.color == 1 && p.curr_pos_x < 1.5) {
+        p.curr_pos_x = BOARD_LENGTH - 1.0;
+        p.destination_y = p.curr_pos_y;
+        p.w_x = -p.desired_v0;
+        p.w_y = 0.0;
+        p.curr_vel_x = p.w_x;
+        p.curr_vel_y = p.w_y;
+      }
+
+      p.curr_pos_y = clamp(p.curr_pos_y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS);
+
+      float e_x = p.destination_x - p.curr_pos_x;
+      float e_y = p.destination_y - p.curr_pos_y;
+
+      float e_norm = sqrt(e_x * e_x + e_y * e_y);
+      float e_norm2 = e_norm > EPS ? e_norm : EPS;
+      float e_norm_inv = 1.0 / e_norm2;
+
+      p.e_x = e_x * e_norm_inv;
+      p.e_y = e_y * e_norm_inv;
+
+      float velocity_norm_sq = p.curr_vel_x * p.curr_vel_x + p.curr_vel_y * p.curr_vel_y;
+
+      float velocity_norm = sqrt(velocity_norm_sq);
+
+      p.s_b_sq = velocity_norm_sq * delta_t_sq;
+
+      p.s_b_e_x = p.e_x * (velocity_norm * delta_t);
+      p.s_b_e_y = p.e_y * (velocity_norm * delta_t);
+
+#if !BENCHMARK
+      if (dump_path && (t % dump_interval == 0)) {
+        samples.emplace_back(t + 1, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color);
+      }
+#endif
+    }
+  }
+
+#if !BENCHMARK
+  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';
+    }
+  }
+#endif
+
+  return 0;
+}

variants/variant_inlining_8.cpp

@@ -0,0 +1,280 @@
+// # This file performs in-lining for run_simulation.
+// # Changes: Replaced references to Vector2d structs inside each Obstacle of Board.obstacles with references to doubles
+// # Predecessor: variant_inlining_7.cpp
+
+#include "../headers/board_reduced.hpp"
+#include "../headers/constants.hpp"
+#include "../headers/pedestrian_reduced.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <math.h>
+#include <random>
+#include <tuple>
+#include <vector>
+
+#define SEED 42
+#define BOARD_WIDTH 10.0
+#define BOARD_LENGTH 50.0
+#define WALL_EPS 1e-3
+#define BOARD_WIDTH_MINUS_WALL_EPS 9.999
+
+#ifndef RANDOM_SEED
+#define RANDOM_SEED 0
+#endif
+
+#ifndef BENCHMARK
+#define BENCHMARK 0
+#endif
+
+using namespace std;
+
+Pedestrian *init_pedestrians(int n_p) {
+#if RANDOM_SEED
+  random_device rd;
+  mt19937 gen(rd());
+#else
+  mt19937 gen(SEED);
+#endif
+
+  normal_distribution<double> velocity_gen(v0_mean, v0_std_dev);
+  uniform_real_distribution<double> y_gen(0.3, BOARD_WIDTH - 0.3);
+  uniform_real_distribution<double> x_gen_left(2.0, BOARD_LENGTH / 2.0);
+  uniform_real_distribution<double> x_gen_right(BOARD_LENGTH / 2.0, BOARD_LENGTH - 2.0);
+
+  Pedestrian *P = new Pedestrian[n_p];
+
+  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);
+
+    double start_pos_x = x;
+    double start_pos_y = y;
+
+    double destination_x = BOARD_LENGTH + 30.0;
+    double destination_y = y;
+
+    double init_e_x = destination_x - start_pos_x;
+    double init_e_y = destination_y - start_pos_y;
+    double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y);
+    double init_e_unit_x = init_e_x / init_e_norm;
+    double init_e_unit_y = init_e_y / 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(0, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y);
+  }
+
+  for (int i = n_p_half; i < n_p; i++) {
+    double x = x_gen_right(gen);
+    double y = y_gen(gen);
+
+    double start_pos_x = x;
+    double start_pos_y = y;
+    double destination_x = -30;
+    double destination_y = y;
+
+    double init_e_x = destination_x - start_pos_x;
+    double init_e_y = destination_y - start_pos_y;
+    double init_e_norm = sqrt(init_e_x * init_e_x + init_e_y * init_e_y);
+    double init_e_unit_x = init_e_x / init_e_norm;
+    double init_e_unit_y = init_e_y / 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(1, 0, 0, start_pos_x, start_pos_y, init_e_unit_x, init_e_unit_y, v_max, desired_v0, destination_x, destination_y);
+  }
+
+  return P;
+}
+
+Board init_board_without_obstacles() { return Board(BOARD_LENGTH, BOARD_WIDTH, 0, nullptr); }
+
+int run_simulation(Pedestrian *pedestrians, Board &board, int n_p, int n_iter, const char *dump_path) {
+
+#if BENCHMARK
+  (void)dump_path;
+#endif
+
+#if !BENCHMARK
+  const int dump_interval = 1;
+
+  vector<tuple<int, int, double, double, double, double, int>> samples;
+
+  if (dump_path) {
+    samples.reserve(static_cast<size_t>(n_iter / dump_interval + 1) * n_p);
+
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+      samples.emplace_back(0, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color);
+    }
+  }
+#endif
+
+  for (int t = 0; t < n_iter; t++) {
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+
+      double desired_velocity_x = p.e_x * p.desired_v0;
+      double desired_velocity_y = p.e_y * p.desired_v0;
+
+      double f_alpha_dt_x = (desired_velocity_x - p.curr_vel_x) * dt_times_tau_alpha_inv;
+      double f_alpha_dt_y = (desired_velocity_y - p.curr_vel_y) * dt_times_tau_alpha_inv;
+
+      for (int j = 0; j < n_p; j++) {
+        if (i == j)
+          continue;
+
+        double r_ab_x = p.curr_pos_x - pedestrians[j].curr_pos_x;
+        double r_ab_y = p.curr_pos_y - pedestrians[j].curr_pos_y;
+
+        double r_ab_norm = sqrt(r_ab_x * r_ab_x + r_ab_y * r_ab_y);
+
+        double r_ab2_x = r_ab_x - pedestrians[j].s_b_e_x;
+        double r_ab2_y = r_ab_y - pedestrians[j].s_b_e_y;
+
+        double r_ab2_norm = sqrt(r_ab2_x * r_ab2_x + r_ab2_y * r_ab2_y);
+        double r_ab_norm_sum = r_ab_norm + r_ab2_norm;
+        double two_b_sq = r_ab_norm_sum * r_ab_norm_sum - pedestrians[j].s_b_sq;
+        double two_b = sqrt(two_b_sq);
+        double scalar_mult_common = (r_ab_norm_sum / (two_b * r_ab_norm * r_ab2_norm)) * (NEG_V0_AB_div_TWO_SIGMA * exp(two_b * NEG_TWO_SIGMA_INV));
+        double scalar_mult_r_ab = scalar_mult_common * r_ab2_norm;
+        double scalar_mult_r_ab2 = scalar_mult_common * r_ab_norm;
+
+        double f_alpha_beta_neg_x = r_ab_x * scalar_mult_r_ab + r_ab2_x * scalar_mult_r_ab2;
+        double f_alpha_beta_neg_y = r_ab_y * scalar_mult_r_ab + r_ab2_y * scalar_mult_r_ab2;
+
+        double w_val = dt_neg;
+        double a = (p.e_x * f_alpha_beta_neg_x + p.e_y * f_alpha_beta_neg_y) * cos_phi_inv;
+        double b_sq = f_alpha_beta_neg_x * f_alpha_beta_neg_x + f_alpha_beta_neg_y * f_alpha_beta_neg_y;
+        
+        if (a >= 0 && (a * a > b_sq)) {
+          w_val = c_neg_times_dt;
+        }
+
+        f_alpha_dt_x = f_alpha_dt_x + f_alpha_beta_neg_x * w_val;
+        f_alpha_dt_y = f_alpha_dt_y + f_alpha_beta_neg_y * w_val;
+      }
+
+      for (int k = 0; k < board.n_obstacles; k++) {
+
+        double v1_x = p.curr_pos_x - board.obstacles[k].start_x;
+        double v1_y = p.curr_pos_y - board.obstacles[k].start_y;
+
+        double projection = (v1_x * board.obstacles[k].obstacle_v_x + v1_y * board.obstacles[k].obstacle_v_y) * board.obstacles[k].length_sq_inv;
+
+        if (projection < 0.0)
+          projection = 0.0;
+        else if (projection > 1.0)
+          projection = 1.0;
+
+        double r_aB_x = p.curr_pos_x - (board.obstacles[k].start_x + board.obstacles[k].obstacle_v_x * projection);
+        double r_aB_y = p.curr_pos_y - (board.obstacles[k].start_y + board.obstacles[k].obstacle_v_y * projection);
+
+        double r_aB_norm = sqrt(r_aB_x * r_aB_x + r_aB_y * r_aB_y);
+        double scalar_mult = (NEG_U0_AB_div_R * exp(r_aB_norm * NEG_R_INV)) / r_aB_norm;
+        double f_alpha_B_neg_x = r_aB_x * scalar_mult;
+        double f_alpha_B_neg_y = r_aB_y * scalar_mult;
+
+        double w_val;
+        double e_dot_f = (f_alpha_B_neg_y == 0) ? 0 : (f_alpha_B_neg_y < 0 ? -p.e_y : p.e_y);
+        if (e_dot_f >= cos_phi)
+          w_val = dt_neg;
+        else
+          w_val = c_neg_times_dt;
+
+        double f_alpha_B_dt_x = f_alpha_B_neg_x * w_val;
+        double f_alpha_B_dt_y = f_alpha_B_neg_y * w_val;
+
+        f_alpha_dt_x += f_alpha_B_dt_x;
+        f_alpha_dt_y += f_alpha_B_dt_y;
+      }
+
+      p.w_x += f_alpha_dt_x;
+      p.w_y += f_alpha_dt_y;
+    }
+
+    for (int i = 0; i < n_p; i++) {
+      Pedestrian &p = pedestrians[i];
+
+      double w_norm = sqrt(p.w_x * p.w_x + p.w_y * p.w_y);
+
+      double correction_factor = min(1.0, p.v_max / w_norm);
+
+      p.curr_vel_x = p.w_x * correction_factor;
+      p.curr_vel_y = p.w_y * correction_factor;
+
+      p.curr_pos_x = p.curr_pos_x + p.curr_vel_x * dt;
+      p.curr_pos_y = p.curr_pos_y + p.curr_vel_y * dt;
+
+      if (p.color == 0 && p.curr_pos_x > (BOARD_LENGTH - 1.5)) {
+        p.curr_pos_x = 1.0;
+        p.destination_y = p.curr_pos_y;
+        p.w_x = p.desired_v0;
+        p.w_y = 0.0;
+        p.curr_vel_x = p.w_x;
+        p.curr_vel_y = p.w_y;
+      }
+      if (p.color == 1 && p.curr_pos_x < 1.5) {
+        p.curr_pos_x = BOARD_LENGTH - 1.0;
+        p.destination_y = p.curr_pos_y;
+        p.w_x = -p.desired_v0;
+        p.w_y = 0.0;
+        p.curr_vel_x = p.w_x;
+        p.curr_vel_y = p.w_y;
+      }
+
+      p.curr_pos_y = clamp(p.curr_pos_y, WALL_EPS, BOARD_WIDTH_MINUS_WALL_EPS);
+
+      double e_x = p.destination_x - p.curr_pos_x;
+      double e_y = p.destination_y - p.curr_pos_y;
+
+      double e_norm_inv = 1.0 / sqrt(e_x * e_x + e_y * e_y);
+
+      p.e_x = e_x * e_norm_inv;
+      p.e_y = e_y * e_norm_inv;
+
+      double velocity_norm_sq = p.curr_vel_x * p.curr_vel_x + p.curr_vel_y * p.curr_vel_y;
+
+      double velocity_norm = sqrt(velocity_norm_sq);
+
+      p.s_b_sq = velocity_norm_sq * delta_t_sq;
+
+      p.s_b_e_x = p.e_x * (velocity_norm * delta_t);
+      p.s_b_e_y = p.e_y * (velocity_norm * delta_t);
+
+#if !BENCHMARK
+      if (dump_path && (t % dump_interval == 0)) {
+        samples.emplace_back(t + 1, i, p.curr_pos_x, p.curr_pos_y, p.e_x, p.e_y, p.color);
+      }
+#endif
+    }
+  }
+
+#if !BENCHMARK
+  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';
+    }
+  }
+#endif
+
+  return 0;
+}