Framework for simulations

main.cpp

@@ -4,7 +4,10 @@
 #include "pedestrian.hpp"
 #include "tsc_x86.h"
 #include "vector2d.hpp"
+#include <fstream>
 #include <iostream>
+#include <tuple>
+#include <vector>
 
 #define BOARD_WIDTH 10.0
 #define BOARD_LENGTH 50.0
@@ -22,22 +25,33 @@ using namespace std;
 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
+  // std::mt19937 gen(rd());                                    // use rd instead of 42 as seed for real randomness
+  std::mt19937 gen(42);                                      // use fixed seed for reproducibility
   std::normal_distribution<double> velocity_gen(1.34, 0.26); // 1.34 and 0.26 are the mean and stdrd deviation used in paper
   std::uniform_real_distribution<double> x_gen(2, BOARD_LENGTH - 2);
   std::uniform_real_distribution<double> y_gen(2, BOARD_WIDTH - 2);
   for (int i = 0; i < n_p / 2; i++) {
-    Vector2d start_pos = Vector2d(1, y_gen(gen));
+    double y = y_gen(gen);
-    Vector2d destination = Vector2d(BOARD_LENGTH - 1, y_gen(gen));
+    Vector2d start_pos = Vector2d(x_gen(gen), y);
-    double desired_v0 = velocity_gen(gen);
+    Vector2d destination = Vector2d(BOARD_LENGTH + 300, BOARD_WIDTH / 2);
+    // double desired_v0 = velocity_gen(gen);
+    double desired_v0;
+    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)
     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());
   }
 
   for (int i = n_p / 2; i < n_p; i++) {
-    Vector2d start_pos = Vector2d(BOARD_LENGTH - 1, y_gen(gen));
+    double y = y_gen(gen);
-    Vector2d destination = Vector2d(1, y_gen(gen));
+    Vector2d start_pos = Vector2d(x_gen(gen), y);
-    double desired_v0 = velocity_gen(gen);
+    Vector2d destination = Vector2d(-300, BOARD_WIDTH / 2);
+    // double desired_v0 = velocity_gen(gen);
+    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));
   }
@@ -51,6 +65,15 @@ int main(int argc, const char *argv[]) {
 
   int n_p = atoi(argv[1]);
   int n_iter = atoi(argv[2]);
+
+  // Logic to get/dump the simulation data
+  const int dump_interval = 5; // dump every 5 iterations -> 20 fps (since dt = 0.01)
+  const char *dump_path = (argc >= 4) ? argv[3] : nullptr;
+
+  std::vector<std::tuple<int, int, double, double, int>> samples; // t, id, x, y, color
+  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
@@ -88,8 +111,26 @@ int main(int argc, const char *argv[]) {
       // # 4 flops = 2 ADDs, 2 MULTs
       p.curr_pos = p.curr_pos + p.curr_vel * dt;
 
+      // how to handle if a pedestrian reached the end of our simulation/end of the hallway
+      // 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;
+      }
+      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;
+      }
+
       // # 10 flops = 5 ADDs, 2 DIVs, 1 SQRT, 2 MULTs
       p.e = (p.destination - p.curr_pos) / (p.destination - p.curr_pos).norm();
+
+      // if we have a dump path/file we wanna save the data for visualisation
+      if (dump_path && (t % dump_interval == 0)) {
+        samples.emplace_back(t, i, p.curr_pos.x, p.curr_pos.y, p.color);
+      }
     }
   }
   long long W = (long long)(8 + 55 * n_p - 55 + 192 + 4 + 8 + 4 + 10) * n_p * n_iter;
@@ -108,6 +149,15 @@ int main(int argc, const char *argv[]) {
   std::cout << P;
   std::cout << " flops/cycle\n";
 
+  // now we wanna save the data to a csv
+  if (dump_path) {
+    std::ofstream out(dump_path);
+    out << "time,id,x,y,color\n";
+    for (auto &[t, id, x, y, c] : samples) {
+      out << t << ',' << id << ',' << x << ',' << y << ',' << c << '\n';
+    }
+  }
+
   delete[] pedestrians;
   return 0;
 }

visualization.py

@@ -0,0 +1,52 @@
+import sys
+import pandas as pd
+import matplotlib.pyplot as plt
+
+from matplotlib.animation import FuncAnimation
+
+
+
+
+if len(sys.argv)<2:
+    sys.exit("Use it like the following: python visualization.py positions.csv (out.gif)")
+
+
+df = pd.read_csv(sys.argv[1])
+frames = sorted(df["time"].unique())
+
+fig, ax = plt.subplots(figsize = (10,2.2))
+ax.set_xlim(0, 50); ax.set_ylim(0, 10)
+ax.set_aspect("equal"); ax.set_facecolor("white")
+ax.set_xticks([]); ax.set_yticks([])
+for s in ("top", "right", "left", "bottom"):
+    ax.spines[s].set_color("gray")
+
+
+# Lets define the 2 groups
+# WHITE = 1 (left -> right, drawn as filled black circles, like paper)
+# BLACK = 0 (right -> left, drawn as empty circles)
+s_white = ax.scatter([],[], s=30, facecolors = "black", edgecolor = "black", zorder = 2)
+s_black = ax.scatter([],[], s=30, facecolors = "white", edgecolor = "black", zorder = 2)
+
+
+
+# update our graphic from t to t+1
+title = ax.set_title("")
+def update(t):
+    sub = df[df["time"] == t]
+    w = sub[sub["color"] == 1]
+    b = sub[sub["color"] == 0]
+    s_white.set_offsets(w[["x", "y"]].values)
+    s_black.set_offsets(b[["x", "y"]].values)
+    title.set_text(f"t = {t}")
+    return s_white, s_black, title
+
+anim = FuncAnimation(fig,update, frames = frames, interval=50, blit=False)
+
+
+
+if len(sys.argv) >= 3:
+    anim.save(sys.argv[2], writer="ffmpeg", fps=20)
+    print(f"saved {sys.argv[2]}")
+else:
+    plt.show()