Draft

examples/simulate_grid.py

@@ -0,0 +1,110 @@
+import argparse
+from datetime import datetime
+import logging
+
+from dsf.cartography import create_manhattan_cartography
+from dsf.cartography.cartography import get_cartography
+from dsf.mobility import (
+    RoadNetwork,
+    Dynamics,
+    AgentInsertionMethod,
+    PathWeight,
+    SpeedFunction,
+)
+
+from tqdm import trange
+from numba import cfunc, float64
+import numpy as np
+import networkx as nx
+
+
+@cfunc(float64(float64, float64), nopython=True, cache=True)
+def custom_speed(max_speed, density):
+    if density < 0.35:
+        return max_speed * (0.9 - 0.1 * density)
+    return max_speed * (1.2 - 0.7 * density)
+
+
+logging.basicConfig(
+    level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s"
+)
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--seed", type=int, default=69, help="Random seed for reproducibility"
+    )
+    parser.add_argument(
+        "--dim", type=str, default="12x12", help="Dimensions of the grid (e.g., 10x10)"
+    )
+    parser.add_argument("--amp", type=int, help="Amplitude of the vehicle input")
+    args = parser.parse_args()
+    np.random.seed(args.seed)
+
+    # Parse the grid dimensions
+    try:
+        rows, cols = map(int, args.dim.split("x"))
+    except ValueError:
+        raise ValueError(
+            "Invalid grid dimensions. Please use the format 'rowsxcols' (e.g., 10x10)."
+        )
+
+    logging.info(f"Creating manhattan cartography for {rows}x{cols} grid...")
+    # Get the cartography of the specified city
+    df_edges, df_nodes = create_manhattan_cartography(rows, cols)
+
+    df_nodes["type"] = (
+        "traffic_signals"  # Set all nodes as traffic lights for simplicity
+    )
+
+    df_edges.to_csv(f"grid_{args.dim}_edges.csv", sep=";", index=False)
+    df_nodes.to_csv(f"grid_{args.dim}_nodes.csv", sep=";", index=False)
+
+    del df_edges, df_nodes
+
+    logging.info("Creating road network and dynamics model...")
+
+    # Create a road network from the cartography
+    road_network = RoadNetwork()
+    road_network.importEdges(f"grid_{args.dim}_edges.csv", ";")
+    road_network.importNodeProperties(f"grid_{args.dim}_nodes.csv", ";")
+    # Adjust network parameters
+    road_network.adjustNodeCapacities()
+    road_network.autoMapStreetLanes()
+    road_network.autoAssignRoadPriorities()
+    road_network.autoInitTrafficLights()
+    road_network.describe()
+
+    # Generaate a random vector of integer values for vehicle input
+    # We want values to have a 10s entry for a whole day
+    vehicle_input = np.random.normal(args.amp, args.amp * 0.1, size=8640)
+    vehicle_input = np.clip(vehicle_input, 0, None).astype(int)
+
+    # Create a dynamics model for the road network
+    dynamics = Dynamics(road_network, seed=args.seed)
+    # To use a custom speed function, you must pass the pointer to the compiled function using the address attribute
+    # dynamics.setSpeedFunction(SpeedFunction.CUSTOM, custom_speed.address)
+    # Get epoch time of today at midnight
+    epoch_time = int(
+        datetime.combine(datetime.today(), datetime.min.time()).timestamp()
+    )
+
+    dynamics.setMeanTravelDistance(10e3)  # Set mean travel distance to 10 km
+    dynamics.killStagnantAgents(40.0)
+    dynamics.setInitTime(epoch_time)
+    dynamics.connectDataBase(f"grid_{args.dim}.db")
+    dynamics.saveData(300, True, True, True)
+
+    # Simulate traffic for 24 hours with a time step of 10 seconds
+    for time_step in trange(86400):
+        # Update paths every 5 minutes (300 seconds)
+        if time_step % 300 == 0:
+            dynamics.updatePaths()
+        # Add agents every 10 seconds
+        if time_step % 10 == 0:
+            dynamics.addAgents(
+                vehicle_input[time_step // 10], AgentInsertionMethod.RANDOM
+            )
+        dynamics.evolve(False)
+
+    dynamics.summary()

src/dsf/cartography/cartography.py

@@ -342,6 +342,7 @@ def create_manhattan_cartography(
         n_x (int): Number of nodes in the x-direction (longitude). Defaults to 10.
         n_y (int): Number of nodes in the y-direction (latitude). Defaults to 10.
         spacing (float): Distance between nodes in meters. Defaults to 2000.0.
+        maxspeed (float): Maximum speed for all edges in km/h. Defaults to 50.0.
         center_lat (float): Latitude of the network center. Defaults to 0.0.
         center_lon (float): Longitude of the network center. Defaults to 0.0.
 

src/dsf/mobility/FirstOrderDynamics.cpp

@@ -870,7 +870,9 @@ namespace dsf::mobility {
         "mean_speed_kph REAL, "
         "std_speed_kph REAL, "
         "mean_density_vpk REAL NOT NULL, "
-        "std_density_vpk REAL NOT NULL)");
+        "std_density_vpk REAL NOT NULL, "
+        "mean_travel_time_s REAL, "
+        "mean_queue_length REAL NOT NULL)");
 
     spdlog::info("Initialized avg_stats table in the database.");
   }
@@ -1412,7 +1414,8 @@ namespace dsf::mobility {
 
   void FirstOrderDynamics::evolve(bool const reinsert_agents) {
     auto const n_threads{std::max<std::size_t>(1, this->concurrency())};
-    std::atomic<double> mean_speed{0.}, mean_density{0.};
+    std::atomic<double> mean_speed{0.}, mean_density{0.}, mean_traveltime{0.},
+        mean_queue_length{0.};
     std::atomic<double> std_speed{0.}, std_density{0.};
     std::atomic<std::size_t> nValidEdges{0};
     bool const bComputeStats = this->database() != nullptr &&
@@ -1426,6 +1429,7 @@ namespace dsf::mobility {
       std::optional<std::string> coilName;
       double density;
       std::optional<double> avgSpeed;
+      std::optional<double> avgTravelTime;
       std::optional<double> stdSpeed;
       std::optional<std::size_t> nObservations;
       std::optional<std::size_t> counts;
@@ -1466,6 +1470,7 @@ namespace dsf::mobility {
                 m_evolveStreet(pStreet, reinsert_agents);
                 if (bComputeStats) {
                   auto const& density{pStreet->density() * 1e3};
+                  auto const& queueLength{pStreet->nExitingAgents()};
 
                   auto const speedMeasure = pStreet->meanSpeed(true);
                   if (speedMeasure.is_valid) {
@@ -1475,13 +1480,15 @@ namespace dsf::mobility {
                       mean_speed.fetch_add(speed, std::memory_order_relaxed);
                       std_speed.fetch_add(speed * speed + speed_std * speed_std,
                                           std::memory_order_relaxed);
-
+                      mean_traveltime.fetch_add(pStreet->length() / speedMeasure.mean,
+                                                std::memory_order_relaxed);
                       ++nValidEdges;
                     }
                   }
                   if (m_bSaveAverageStats) {
                     mean_density.fetch_add(density, std::memory_order_relaxed);
                     std_density.fetch_add(density * density, std::memory_order_relaxed);
+                    mean_queue_length.fetch_add(queueLength, std::memory_order_relaxed);
                   }
 
                   if (m_bSaveStreetData) {
@@ -1499,7 +1506,7 @@ namespace dsf::mobility {
                       record.stdSpeed = speedMeasure.std * 3.6;
                       record.nObservations = speedMeasure.n;
                     }
-                    record.queueLength = pStreet->nExitingAgents();
+                    record.queueLength = queueLength;
                     streetDataRecords.push_back(record);
                   }
                 }
@@ -1589,6 +1596,8 @@ namespace dsf::mobility {
       if (m_bSaveAverageStats) {  // Average Stats Table
         mean_speed.store(mean_speed.load() / nValidEdges.load());
         mean_density.store(mean_density.load() / numEdges);
+        mean_traveltime.store(mean_traveltime.load() / nValidEdges.load());
+        mean_queue_length.store(mean_queue_length.load() / numEdges);
         {
           double std_speed_val = std_speed.load();
           double mean_speed_val = mean_speed.load();
@@ -1605,8 +1614,9 @@ namespace dsf::mobility {
             *this->database(),
             "INSERT INTO avg_stats ("
             "simulation_id, datetime, time_step, n_ghost_agents, n_agents, "
-            "mean_speed_kph, std_speed_kph, mean_density_vpk, std_density_vpk) "
-            "VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?)");
+            "mean_speed_kph, std_speed_kph, mean_density_vpk, std_density_vpk, "
+            "mean_travel_time_s, mean_queue_length) "
+            "VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)");
         insertStmt.bind(1, static_cast<std::int64_t>(this->id()));
         insertStmt.bind(2, this->strDateTime());
         insertStmt.bind(3, static_cast<std::int64_t>(this->time_step()));
@@ -1621,6 +1631,8 @@ namespace dsf::mobility {
         }
         insertStmt.bind(8, mean_density);
         insertStmt.bind(9, std_density);
+        insertStmt.bind(10, mean_traveltime);
+        insertStmt.bind(11, mean_queue_length);
         insertStmt.exec();
       }
       // Special case: if m_savingInterval == 0, it was a triggered saveData() call, so we need to reset all flags

webapp/index.html

@@ -76,6 +76,15 @@ <h2>Load SQLite Database</h2>
 
     <!-- Legend container -->
     <div class="legend-container">
+      <div class="legend-controls">
+        <label for="edgeColorObservableSelector">Edge color by</label>
+        <select id="edgeColorObservableSelector">
+          <option value="density">density</option>
+          <option value="speed">speed</option>
+          <option value="traveltime">traveltime</option>
+          <option value="queue_length">queue_length</option>
+        </select>
+      </div>
       <div class="legend-title">Density</div>
       <div class="legend-bar"></div>
       <div class="legend-labels">