Add Street::meanSpeed and Dynamics::saveStreetSpeeds functions

benchmark/Bench_Street.cpp

@@ -25,7 +25,7 @@ static void BM_Street_AddAgent(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (auto _ : state) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
   }
 }
 
@@ -35,7 +35,7 @@ static void BM_Street_Enqueue(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
   }
   size_t queueId = 0;
   for (auto _ : state) {
@@ -52,13 +52,13 @@ static void BM_Street_Dequeue(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
     street.enqueue(0);
   }
   size_t index = 0;
   for (auto _ : state) {
     if (!street.queue(index).empty()) {
-      auto agent = street.dequeue(index);
+      auto agent = street.dequeue(index, spawnTime++);
       benchmark::DoNotOptimize(agent);
     }
   }
@@ -70,7 +70,7 @@ static void BM_Street_nAgents(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
     if (i % 2 == 0)
       street.enqueue(0);
   }
@@ -86,7 +86,7 @@ static void BM_Street_Density(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
     if (i % 2 == 0)
       street.enqueue(0);
   }
@@ -102,7 +102,7 @@ static void BM_Street_nMovingAgents(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
   }
   for (auto _ : state) {
     int n = street.nMovingAgents();
@@ -116,7 +116,7 @@ static void BM_Street_nExitingAgents(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
     street.enqueue(0);
   }
   for (auto _ : state) {
@@ -142,7 +142,7 @@ static void BM_CoilStreet_AddAgent(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (auto _ : state) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
   }
 }
 
@@ -153,10 +153,10 @@ static void BM_CoilStreet_MeanFlow(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
     street.enqueue(0);
     if (i % 2 == 0) {
-      auto dequeued = street.dequeue(0);
+      auto dequeued = street.dequeue(0, spawnTime++);
     }
   }
   for (auto _ : state) {
@@ -172,13 +172,13 @@ static void BM_CoilStreet_Dequeue(benchmark::State& state) {
   auto pItinerary = std::make_shared<dsf::mobility::Itinerary>(1, 1);
   for (int i = 0; i < 50; ++i) {
     auto agent = std::make_unique<dsf::mobility::Agent>(0, spawnTime++, pItinerary, 0);
-    street.addAgent(std::move(agent));
+    street.addAgent(std::move(agent), spawnTime);
     street.enqueue(0);
   }
   size_t index = 0;
   for (auto _ : state) {
     if (!street.queue(index).empty()) {
-      auto agent = street.dequeue(index);
+      auto agent = street.dequeue(index, spawnTime++);
       benchmark::DoNotOptimize(agent);
     }
   }

src/dsf/base/Dynamics.hpp

@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Network.hpp"
+#include "../utility/Measurement.hpp"
 #include "../utility/Typedef.hpp"
 
 #include <algorithm>
@@ -36,33 +37,6 @@
 #include <tbb/tbb.h>
 
 namespace dsf {
-  /// @brief The Measurement struct represents the mean of a quantity and its standard deviation
-  /// @tparam T The type of the quantity
-  /// @param mean The mean
-  /// @param std The standard deviation of the sample
-  template <typename T>
-  struct Measurement {
-    T mean;
-    T std;
-
-    Measurement(T mean, T std) : mean{mean}, std{std} {}
-    Measurement(std::span<T> data) {
-      auto x_mean = static_cast<T>(0), x2_mean = static_cast<T>(0);
-      if (data.empty()) {
-        mean = static_cast<T>(0);
-        std = static_cast<T>(0);
-        return;
-      }
-
-      std::for_each(data.begin(), data.end(), [&x_mean, &x2_mean](auto value) -> void {
-        x_mean += value;
-        x2_mean += value * value;
-      });
-      mean = x_mean / data.size();
-      std = std::sqrt(x2_mean / data.size() - mean * mean);
-    }
-  };
-
   /// @brief The Dynamics class represents the dynamics of the network.
   /// @tparam network_t The type of the network
   template <typename network_t>

src/dsf/bindings.cpp

@@ -681,9 +681,15 @@ PYBIND11_MODULE(dsf_cpp, m) {
           "saveStreetDensities",
           &dsf::mobility::FirstOrderDynamics::saveStreetDensities,
           pybind11::arg("filename"),
-          pybind11::arg("normalized") = true,
           pybind11::arg("separator") = ';',
+          pybind11::arg("normalized") = true,
           dsf::g_docstrings.at("dsf::mobility::RoadDynamics::saveStreetDensities").c_str())
+      .def("saveStreetSpeeds",
+           &dsf::mobility::FirstOrderDynamics::saveStreetSpeeds,
+           pybind11::arg("filename"),
+           pybind11::arg("separator") = ';',
+           pybind11::arg("normalized") = false,
+           dsf::g_docstrings.at("dsf::mobility::RoadDynamics::saveStreetSpeeds").c_str())
       .def("saveCoilCounts",
            &dsf::mobility::FirstOrderDynamics::saveCoilCounts,
            pybind11::arg("filename"),

src/dsf/dsf.hpp

@@ -6,7 +6,7 @@
 
 static constexpr uint8_t DSF_VERSION_MAJOR = 4;
 static constexpr uint8_t DSF_VERSION_MINOR = 7;
-static constexpr uint8_t DSF_VERSION_PATCH = 7;
+static constexpr uint8_t DSF_VERSION_PATCH = 8;
 
 static auto const DSF_VERSION =
     std::format("{}.{}.{}", DSF_VERSION_MAJOR, DSF_VERSION_MINOR, DSF_VERSION_PATCH);

src/dsf/mobility/RoadDynamics.hpp

@@ -358,11 +358,18 @@ namespace dsf::mobility {
 
     /// @brief Save the street densities in csv format
     /// @param filename The name of the file (default is "{datetime}_{simulation_name}_street_densities.csv")
-    /// @param normalized If true, the densities are normalized in [0, 1]
     /// @param separator The separator character (default is ';')
+    /// @param normalized If true, the densities are normalized in [0, 1] dividing by the street capacity attribute
     void saveStreetDensities(std::string filename = std::string(),
-                             bool normalized = true,
-                             char const separator = ';') const;
+                             char const separator = ';',
+                             bool const normalized = true) const;
+    /// @brief Save the street speeds in csv format
+    /// @param filename The name of the file (default is "{datetime}_{simulation_name}_street_speeds.csv")
+    /// @param separator The separator character (default is ';')
+    /// @param bNormalized If true, the speeds are normalized in [0, 1] dividing by the street maxSpeed attribute
+    void saveStreetSpeeds(std::string filename = std::string(),
+                          char const separator = ';',
+                          bool bNormalized = false) const;
     /// @brief Save the street input counts in csv format
     /// @param filename The name of the file
     /// @param reset If true, the input counts are cleared after the computation
@@ -742,7 +749,7 @@ namespace dsf::mobility {
                 timeTolerance,
                 timeDiff);
             // Kill the agent
-            this->m_killAgent(pStreet->dequeue(queueIndex));
+            this->m_killAgent(pStreet->dequeue(queueIndex, this->time_step()));
             ++m_nKilledAgents;
             continue;
           }
@@ -898,7 +905,8 @@ namespace dsf::mobility {
           }
         }
         if (bArrived) {
-          auto pAgent = this->m_killAgent(pStreet->dequeue(queueIndex));
+          auto pAgent =
+              this->m_killAgent(pStreet->dequeue(queueIndex, this->time_step()));
           ++m_nArrivedAgents;
           if (reinsert_agents) {
             // reset Agent's values
@@ -920,7 +928,7 @@ namespace dsf::mobility {
               *nextStreet);
           continue;
         }
-        auto pAgent{pStreet->dequeue(queueIndex)};
+        auto pAgent{pStreet->dequeue(queueIndex, this->time_step())};
         spdlog::debug(
             "{} at time {} has been dequeued from street {} and enqueued on street {} "
             "with free time {}.",
@@ -990,7 +998,7 @@ namespace dsf::mobility {
               pNode->id(),
               nextStreet->id(),
               pAgent->freeTime());
-          nextStreet->addAgent(std::move(pAgent));
+          nextStreet->addAgent(std::move(pAgent), this->time_step());
           it = intersection.agents().erase(it);
           break;
         }
@@ -1018,7 +1026,7 @@ namespace dsf::mobility {
               nextStreet->id(),
               pAgent->freeTime(),
               *pAgent);
-          nextStreet->addAgent(std::move(pAgent));
+          nextStreet->addAgent(std::move(pAgent), this->time_step());
         } else {
           return;
         }
@@ -1351,7 +1359,7 @@ namespace dsf::mobility {
       this->setAgentSpeed(pAgent);
       pAgent->setFreeTime(this->time_step() +
                           std::ceil(street->length() / pAgent->speed()));
-      street->addAgent(std::move(pAgent));
+      street->addAgent(std::move(pAgent), this->time_step());
       this->m_agents.pop_back();
     }
   }
@@ -2302,8 +2310,8 @@ namespace dsf::mobility {
   template <typename delay_t>
     requires(is_numeric_v<delay_t>)
   void RoadDynamics<delay_t>::saveStreetDensities(std::string filename,
-                                                  bool normalized,
-                                                  char const separator) const {
+                                                  char const separator,
+                                                  bool const normalized) const {
     if (filename.empty()) {
       filename =
           this->m_safeDateTime() + '_' + this->m_safeName() + "_street_densities.csv";
@@ -2333,6 +2341,48 @@ namespace dsf::mobility {
     file << std::endl;
     file.close();
   }
+  template <typename delay_t>
+    requires(is_numeric_v<delay_t>)
+  void RoadDynamics<delay_t>::saveStreetSpeeds(std::string filename,
+                                               char const separator,
+                                               bool const bNormalized) const {
+    if (filename.empty()) {
+      filename = this->m_safeDateTime() + '_' + this->m_safeName() + "_street_speeds.csv";
+    }
+    bool bEmptyFile{false};
+    {
+      std::ifstream file(filename);
+      bEmptyFile = file.peek() == std::ifstream::traits_type::eof();
+    }
+    std::ofstream file(filename, std::ios::app);
+    if (!file.is_open()) {
+      throw std::runtime_error("Error opening file \"" + filename + "\" for writing.");
+    }
+    if (bEmptyFile) {
+      file << "datetime" << separator << "time_step";
+      for (auto const& [streetId, pStreet] : this->graph().edges()) {
+        file << separator << streetId;
+      }
+      file << std::endl;
+    }
+    file << this->strDateTime() << separator << this->time_step();
+    for (auto const& [streetId, pStreet] : this->graph().edges()) {
+      auto const measure = pStreet->meanSpeed(true);
+      file << separator;
+      // If not valid, write empty value (less space w.r.t. NaN)
+      if (!measure.is_valid) {
+        continue;
+      }
+
+      double speed{measure.mean};
+      if (bNormalized) {
+        speed /= pStreet->maxSpeed();
+      }
+      file << std::fixed << std::setprecision(2) << speed;
+    }
+    file << std::endl;
+    file.close();
+  }
   template <typename delay_t>
     requires(is_numeric_v<delay_t>)
   void RoadDynamics<delay_t>::saveCoilCounts(const std::string& filename,

src/dsf/mobility/Street.cpp

@@ -126,9 +126,10 @@ namespace dsf::mobility {
     m_counter->reset();
   }
 
-  void Street::addAgent(std::unique_ptr<Agent> pAgent) {
+  void Street::addAgent(std::unique_ptr<Agent> pAgent, std::time_t const currentTime) {
     assert(!isFull());
-    spdlog::debug("Adding {} on {}", *pAgent, *this);
+    spdlog::trace("Adding {} on {}", *pAgent, *this);
+    m_agentsInsertionTimes[pAgent->id()] = currentTime;
     m_movingAgents.push(std::move(pAgent));
     if (m_counter.has_value() && m_counterPosition == CounterPosition::ENTRY) {
       ++(*m_counter);
@@ -144,9 +145,15 @@ namespace dsf::mobility {
       ++(*m_counter);
     }
   }
-  std::unique_ptr<Agent> Street::dequeue(std::size_t const& index) {
+  std::unique_ptr<Agent> Street::dequeue(std::size_t const& index,
+                                         std::time_t const currentTime) {
     assert(!m_exitQueues[index].empty());
     auto pAgent{std::move(m_exitQueues[index].front())};
+    // Keep track of average speed
+    m_avgSpeeds.push_back(m_length /
+                          (currentTime - m_agentsInsertionTimes[pAgent->id()]));
+    m_agentsInsertionTimes.erase(pAgent->id());
+
     m_exitQueues[index].pop();
     if (m_counter.has_value() && m_counterPosition == CounterPosition::EXIT) {
       ++(*m_counter);
@@ -202,5 +209,12 @@ namespace dsf::mobility {
     }
     return nAgents;
   }
+  Measurement<double> Street::meanSpeed(bool const bReset) {
+    auto const avgSpeed{Measurement<double>(m_avgSpeeds)};
+    if (bReset) {
+      m_avgSpeeds.clear();
+    }
+    return avgSpeed;
+  }
 
 };  // namespace dsf::mobility