diff --git a/examples/simulate_grid.py b/examples/simulate_grid.py new file mode 100644 index 00000000..e8a0921a --- /dev/null +++ b/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() diff --git a/src/dsf/cartography/cartography.py b/src/dsf/cartography/cartography.py index 392ea465..d00b86d9 100644 --- a/src/dsf/cartography/cartography.py +++ b/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. diff --git a/src/dsf/mobility/FirstOrderDynamics.cpp b/src/dsf/mobility/FirstOrderDynamics.cpp index 5e3efde5..966aefe5 100644 --- a/src/dsf/mobility/FirstOrderDynamics.cpp +++ b/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(1, this->concurrency())}; - std::atomic mean_speed{0.}, mean_density{0.}; + std::atomic mean_speed{0.}, mean_density{0.}, mean_traveltime{0.}, + mean_queue_length{0.}; std::atomic std_speed{0.}, std_density{0.}; std::atomic nValidEdges{0}; bool const bComputeStats = this->database() != nullptr && @@ -1426,6 +1429,7 @@ namespace dsf::mobility { std::optional coilName; double density; std::optional avgSpeed; + std::optional avgTravelTime; std::optional stdSpeed; std::optional nObservations; std::optional 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(this->id())); insertStmt.bind(2, this->strDateTime()); insertStmt.bind(3, static_cast(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 diff --git a/webapp/index.html b/webapp/index.html index 96edb272..35520dc9 100644 --- a/webapp/index.html +++ b/webapp/index.html @@ -76,6 +76,15 @@

Load SQLite Database

+
+ + +
Density
diff --git a/webapp/script.js b/webapp/script.js index 36904be5..d3a8a7d7 100644 --- a/webapp/script.js +++ b/webapp/script.js @@ -822,6 +822,26 @@ let highlightedNode = null; let chart; let db = null; let selectedSimulationId = null; +let edgeObservableData = { + density: [], + speed: [], + traveltime: [], + queue_length: [] +}; +let edgeObservableDomains = { + density: [0, MAX_DENSITY], + speed: [0, 1], + traveltime: [0, 1], + queue_length: [0, 1] +}; +let selectedEdgeColorObservable = 'density'; + +const EDGE_OBSERVABLE_CONFIG = { + density: { label: 'Density' }, + speed: { label: 'Speed', reverseColorScale: true }, + traveltime: { label: 'Travel Time' }, + queue_length: { label: 'Queue Length' } +}; function formatTime(date) { const year = date.getFullYear(); @@ -832,11 +852,6 @@ function formatTime(date) { return `${year}-${month}-${day} ${hours}:${minutes}`; } -// Create a color scale for density values using three color stops -const colorScale = d3.scaleLinear() - .domain([0, MAX_DENSITY / 2, MAX_DENSITY]) - .range(["green", "yellow", "red"]); - // Update node highlight position function updateNodeHighlight() { g.selectAll(".node-highlight").remove(); @@ -906,73 +921,239 @@ function loadEdgesFromDB() { }); } -// Load road_data from SQLite for selected simulation and transform to density format +function getRoadDataColumns() { + const result = db.exec('PRAGMA table_info(road_data)'); + if (result.length === 0) return []; + return result[0].values.map(row => row[1]); +} + +function getTravelTimeExpression() { + const columns = getRoadDataColumns(); + + if (columns.includes('traveltime')) return 'r.traveltime'; + if (columns.includes('travel_time')) return 'r.travel_time'; + if (columns.includes('travel_time_s')) return 'r.travel_time_s'; + + // Fallback: estimate travel time [s] from edge length [m] and avg speed [km/h]. + return 'CASE WHEN r.avg_speed_kph > 0 THEN e.length / (r.avg_speed_kph / 3.6) ELSE 0 END'; +} + +function computeObservableDomain(observableRows) { + const values = observableRows + .flatMap(row => row.values) + .map(v => +v) + .filter(v => Number.isFinite(v)); + + if (values.length === 0) return [0, 1]; + + let minValue = values[0]; + let maxValue = values[0]; + for (let i = 1; i < values.length; i++) { + if (values[i] < minValue) minValue = values[i]; + if (values[i] > maxValue) maxValue = values[i]; + } + + if (minValue === maxValue) { + maxValue = minValue + 1; + } + + return [minValue, maxValue]; +} + +// Load road_data from SQLite for selected simulation and transform to edge-wise time series function loadRoadDataFromDB() { - // Get edge IDs in order const edgeIds = edges.map(e => e.id); - - // Single query to get all data ordered by datetime and street_id + const travelTimeExpression = getTravelTimeExpression(); + const result = db.exec( - `SELECT datetime, street_id, density_vpk FROM road_data WHERE simulation_id = ${selectedSimulationId} ORDER BY datetime, street_id` + `SELECT r.datetime, + r.street_id, + r.density_vpk, + r.avg_speed_kph, + ${travelTimeExpression} AS traveltime, + r.queue_length + FROM road_data r + LEFT JOIN edges e ON e.id = r.street_id + WHERE r.simulation_id = ${selectedSimulationId} + ORDER BY r.datetime, r.street_id` ); - if (result.length === 0) return []; - + + if (result.length === 0) { + return { + densities: [], + observables: { + density: [], + speed: [], + traveltime: [], + queue_length: [] + }, + domains: { + density: [0, 1], + speed: [0, 1], + traveltime: [0, 1], + queue_length: [0, 1] + } + }; + } + const densityData = []; + const speedData = []; + const travelTimeData = []; + const queueLengthData = []; + let currentTs = null; - let currentMap = {}; - + let currentDensityMap = {}; + let currentSpeedMap = {}; + let currentTravelTimeMap = {}; + let currentQueueLengthMap = {}; + for (const row of result[0].values) { - const [ts, streetId, density] = row; + const [ts, streetId, density, speed, travelTime, queueLength] = row; + if (ts !== currentTs) { if (currentTs !== null) { - // Build densities array in same order as edges for previous timestamp - const densityArray = edgeIds.map(id => currentMap[id] || 0); + const densityArray = edgeIds.map(id => +currentDensityMap[id] || 0); + const speedArray = edgeIds.map(id => +currentSpeedMap[id] || 0); + const travelTimeArray = edgeIds.map(id => +currentTravelTimeMap[id] || 0); + const queueLengthArray = edgeIds.map(id => +currentQueueLengthMap[id] || 0); + densityData.push({ datetime: new Date(currentTs), densities: densityArray }); + speedData.push({ + datetime: new Date(currentTs), + values: speedArray + }); + travelTimeData.push({ + datetime: new Date(currentTs), + values: travelTimeArray + }); + queueLengthData.push({ + datetime: new Date(currentTs), + values: queueLengthArray + }); } + currentTs = ts; - currentMap = {}; + currentDensityMap = {}; + currentSpeedMap = {}; + currentTravelTimeMap = {}; + currentQueueLengthMap = {}; } - currentMap[streetId] = density; + + currentDensityMap[streetId] = density; + currentSpeedMap[streetId] = speed; + currentTravelTimeMap[streetId] = travelTime; + currentQueueLengthMap[streetId] = queueLength; } - - // Handle the last timestamp + if (currentTs !== null) { - const densityArray = edgeIds.map(id => currentMap[id] || 0); + const densityArray = edgeIds.map(id => +currentDensityMap[id] || 0); + const speedArray = edgeIds.map(id => +currentSpeedMap[id] || 0); + const travelTimeArray = edgeIds.map(id => +currentTravelTimeMap[id] || 0); + const queueLengthArray = edgeIds.map(id => +currentQueueLengthMap[id] || 0); + densityData.push({ datetime: new Date(currentTs), densities: densityArray }); + speedData.push({ + datetime: new Date(currentTs), + values: speedArray + }); + travelTimeData.push({ + datetime: new Date(currentTs), + values: travelTimeArray + }); + queueLengthData.push({ + datetime: new Date(currentTs), + values: queueLengthArray + }); } - - return densityData; + + const observables = { + density: densityData.map(row => ({ datetime: row.datetime, values: row.densities })), + speed: speedData, + traveltime: travelTimeData, + queue_length: queueLengthData + }; + + const domains = { + density: computeObservableDomain(observables.density), + speed: computeObservableDomain(observables.speed), + traveltime: computeObservableDomain(observables.traveltime), + queue_length: computeObservableDomain(observables.queue_length) + }; + + return { + densities: densityData, + observables, + domains + }; } -// Load global data (aggregated statistics per timestamp) +// Load global data (aggregated statistics per timestamp) from avg-stats table. function loadGlobalDataFromDB() { - // Calculate mean density, avg_speed, etc. per timestamp for selected simulation + const tablesResult = db.exec("SELECT name FROM sqlite_master WHERE type='table'"); + const tableNames = tablesResult.length > 0 ? tablesResult[0].values.map(r => r[0]) : []; + + // Support both historical names. + const avgStatsTable = tableNames.includes('avg_stats') + ? 'avg_stats' + : (tableNames.includes('avgstats') ? 'avgstats' : null); + + if (!avgStatsTable) { + // Fallback for legacy DBs without avg-stats table. + const result = db.exec(` + SELECT datetime, + AVG(density_vpk) as mean_density_vpk, + AVG(avg_speed_kph) as mean_speed_kph, + SUM(counts) as total_counts + FROM road_data + WHERE simulation_id = ${selectedSimulationId} + GROUP BY datetime + ORDER BY datetime + `); + + if (result.length === 0) return []; + + const columns = result[0].columns; + const values = result[0].values; + + return values.map(row => { + const data = { datetime: new Date(row[0]) }; + for (let i = 1; i < columns.length; i++) { + data[columns[i]] = +row[i] || 0; + } + return data; + }); + } + + const colsResult = db.exec(`PRAGMA table_info(${avgStatsTable})`); + if (colsResult.length === 0) return []; + + const allColumns = colsResult[0].values.map(row => row[1]); + const metricColumns = allColumns.filter( + col => !['id', 'simulation_id', 'datetime', 'time_step'].includes(col) + ); + + if (metricColumns.length === 0) return []; + const result = db.exec(` - SELECT datetime, - AVG(density_vpk) as mean_density_vpk, - AVG(avg_speed_kph) as mean_speed_kph, - SUM(counts) as total_counts - FROM road_data + SELECT datetime, ${metricColumns.join(', ')} + FROM ${avgStatsTable} WHERE simulation_id = ${selectedSimulationId} - GROUP BY datetime ORDER BY datetime `); - + if (result.length === 0) return []; - - const columns = result[0].columns; + const values = result[0].values; - return values.map(row => { const data = { datetime: new Date(row[0]) }; - for (let i = 1; i < columns.length; i++) { - data[columns[i]] = +row[i] || 0; + for (let i = 0; i < metricColumns.length; i++) { + data[metricColumns[i]] = +row[i + 1] || 0; } return data; }); @@ -993,7 +1174,10 @@ function getSimulations() { function initializeApp() { // Load data from database edges = loadEdgesFromDB(); - densities = loadRoadDataFromDB(); + const roadDataBundle = loadRoadDataFromDB(); + densities = roadDataBundle.densities; + edgeObservableData = roadDataBundle.observables; + edgeObservableDomains = roadDataBundle.domains; globalData = loadGlobalDataFromDB(); console.log("Loaded edges:", edges.length); @@ -1035,6 +1219,44 @@ function initializeApp() { const canvasEdges = new L.CanvasEdges(edges); canvasEdges.addTo(map); + const edgeColorObservableSelector = document.getElementById('edgeColorObservableSelector'); + selectedEdgeColorObservable = edgeColorObservableSelector.value || 'density'; + + function formatLegendValue(value) { + const numeric = +value; + if (!Number.isFinite(numeric)) return 'N/A'; + if (Math.abs(numeric) >= 100) return numeric.toFixed(0); + if (Math.abs(numeric) >= 10) return numeric.toFixed(1); + return numeric.toFixed(2); + } + + function updateLegend() { + const config = EDGE_OBSERVABLE_CONFIG[selectedEdgeColorObservable] || EDGE_OBSERVABLE_CONFIG.density; + const title = document.querySelector('.legend-title'); + const labels = document.querySelectorAll('.legend-labels span'); + const legendBar = document.querySelector('.legend-bar'); + const domain = edgeObservableDomains[selectedEdgeColorObservable] || [0, 1]; + const middleValue = (domain[0] + domain[1]) / 2; + + title.textContent = config.label; + if (legendBar) { + legendBar.style.background = config.reverseColorScale + ? 'linear-gradient(to right, red, yellow, green)' + : 'linear-gradient(to right, green, yellow, red)'; + } + if (labels.length >= 3) { + labels[0].textContent = formatLegendValue(domain[0]); + labels[1].textContent = formatLegendValue(middleValue); + labels[2].textContent = formatLegendValue(domain[1]); + } + } + + edgeColorObservableSelector.onchange = () => { + selectedEdgeColorObservable = edgeColorObservableSelector.value; + updateLegend(); + updateDensityVisualization(); + }; + let currentChartColumn = 'mean_density_vpk'; // Initialize Chart @@ -1199,23 +1421,33 @@ function initializeApp() { } map.on("zoomend", update); - update(); // Initial render - // Update edge colors based on the current time step density data + // Update edge colors based on the selected time step and observable function updateDensityVisualization() { - const currentDensityRow = densities.find(d => d.datetime.getTime() === timeStamp.getTime()); - if (!currentDensityRow) { - console.error("No density data for time step:", timeStamp); + const currentIndex = densities.findIndex(d => d.datetime.getTime() === timeStamp.getTime()); + if (currentIndex < 0) { + console.error("No road data for time step:", timeStamp); return; } - const currentDensities = currentDensityRow.densities; + + const currentDensities = densities[currentIndex].densities; + const observableRows = edgeObservableData[selectedEdgeColorObservable] || edgeObservableData.density; + const currentObservableValues = observableRows[currentIndex]?.values || currentDensities; + const domain = edgeObservableDomains[selectedEdgeColorObservable] || [0, MAX_DENSITY]; + const colorConfig = EDGE_OBSERVABLE_CONFIG[selectedEdgeColorObservable] || EDGE_OBSERVABLE_CONFIG.density; + const colorRange = colorConfig.reverseColorScale + ? ['red', 'yellow', 'green'] + : ['green', 'yellow', 'red']; + const colorScale = d3.scaleLinear() + .domain([domain[0], (domain[0] + domain[1]) / 2, domain[1]]) + .range(colorRange); const colors = edges.map((edge, index) => { - let density = currentDensities[index]; - if (density === undefined || isNaN(density)) { - density = 0; + let value = currentObservableValues[index]; + if (value === undefined || isNaN(value)) { + value = 0; } - const rgb = d3.rgb(colorScale(density)); + const rgb = d3.rgb(colorScale(value)); return `rgba(${rgb.r}, ${rgb.g}, ${rgb.b}, 0.69)`; }); @@ -1223,6 +1455,9 @@ function initializeApp() { canvasEdges.setDensities(currentDensities); } + updateLegend(); + update(); // Initial render + // Set up the time slider based on the density data's maximum time value const timeSlider = document.getElementById('timeSlider'); const timeLabel = document.getElementById('timeLabel'); diff --git a/webapp/styles.css b/webapp/styles.css index 218de5f2..fea271a5 100644 --- a/webapp/styles.css +++ b/webapp/styles.css @@ -290,6 +290,25 @@ body, html { font-family: sans-serif; font-size: 12px; } + + .legend-controls { + display: flex; + flex-direction: column; + gap: 4px; + margin-bottom: 8px; + } + + .legend-controls label { + font-weight: bold; + } + + .legend-controls select { + padding: 4px; + border: 1px solid #ccc; + border-radius: 3px; + background: white; + font-size: 12px; + } .legend-title { font-weight: bold;