1486 Modernizations and minor fixes in core C++ library

cpp/benchmarks/secir_ageres_setups.h

@@ -30,7 +30,7 @@ namespace mio
 {
 namespace benchmark
 {
-namespace detail
+namespace details
 {
 /**
  * @brief Helper function to create a secir model with consistent setup for use in benchmarking.
@@ -89,7 +89,7 @@ mio::osecir::Model<ScalarType> make_model(int num)
 
     return model;
 }
-} // namespace detail
+} // namespace details
 
 namespace model
 {
@@ -98,7 +98,7 @@ namespace model
  */
 mio::osecir::Model<ScalarType> SecirAgeres(size_t num_agegroups)
 {
-    mio::osecir::Model<ScalarType> model = mio::benchmark::detail::make_model(num_agegroups);
+    mio::osecir::Model<ScalarType> model = mio::benchmark::details::make_model(num_agegroups);
 
     auto nb_groups       = model.parameters.get_num_groups();
     ScalarType cont_freq = 10, fact = 1.0 / (ScalarType)(size_t)nb_groups;
@@ -117,7 +117,7 @@ mio::osecir::Model<ScalarType> SecirAgeres(size_t num_agegroups)
  */
 mio::osecir::Model<ScalarType> SecirAgeresDampings(size_t num_agegroups)
 {
-    mio::osecir::Model<ScalarType> model = mio::benchmark::detail::make_model(num_agegroups);
+    mio::osecir::Model<ScalarType> model = mio::benchmark::details::make_model(num_agegroups);
 
     auto nb_groups       = model.parameters.get_num_groups();
     ScalarType cont_freq = 10, fact = 1.0 / (ScalarType)(size_t)nb_groups;
@@ -145,7 +145,7 @@ mio::osecir::Model<ScalarType> SecirAgeresDampings(size_t num_agegroups)
  */
 mio::osecir::Model<ScalarType> SecirAgeresAbsurdDampings(size_t num_agegroups)
 {
-    mio::osecir::Model<ScalarType> model = mio::benchmark::detail::make_model(num_agegroups);
+    mio::osecir::Model<ScalarType> model = mio::benchmark::details::make_model(num_agegroups);
 
     auto nb_groups       = model.parameters.get_num_groups();
     ScalarType cont_freq = 10, fact = 1.0 / (ScalarType)(size_t)nb_groups;

cpp/memilio/ad/ad.h

@@ -27,6 +27,26 @@
 #include <cmath>
 #include <limits>
 
+namespace ad
+{
+namespace internal
+{
+
+/**
+ * @brief Format AD types (like ad::gt1s<double>::type) using their value for logging with spdlog.
+ *
+ * If derivative information is needed as well, use `ad::derivative(...)` or define a `fmt::formatter<...>`.
+ */
+template <class FP, class DataHandler>
+const FP& format_as(const active_type<FP, DataHandler>& ad_type)
+{
+    // Note: the format_as function needs to be in the same namespace as the value it takes
+    return value(ad_type);
+}
+
+} // namespace internal
+} // namespace ad
+
 // Allow std::numeric_limits to work with AD types.
 template <class FP, class DataHandler>
 struct std::numeric_limits<ad::internal::active_type<FP, DataHandler>> : public numeric_limits<FP> {

cpp/memilio/compartments/compartmental_model.h

@@ -20,8 +20,9 @@
 #ifndef MIO_COMPARTMENTS_COMPARTMENTAL_MODEL_H
 #define MIO_COMPARTMENTS_COMPARTMENTAL_MODEL_H
 
-#include "memilio/config.h"
-#include "memilio/math/eigen.h"
+#include "memilio/config.h" // IWYU pragma: keep
+#include "memilio/math/eigen.h" // IWYU pragma: keep
+
 #include <concepts>
 
 namespace mio

cpp/memilio/compartments/feedback_simulation.h

@@ -20,13 +20,11 @@
 #ifndef MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
 #define MIO_COMPARTMENTS_FEEDBACK_SIMULATION_H
 
-#include <cassert>
-#include "memilio/compartments/simulation.h"
-#include "memilio/utils/time_series.h"
-#include "memilio/utils/parameter_set.h"
 #include "memilio/epidemiology/age_group.h"
-#include "memilio/utils/uncertain_value.h"
 #include "memilio/epidemiology/damping_sampling.h"
+#include "memilio/utils/parameter_set.h"
+#include "memilio/utils/time_series.h"
+#include "memilio/utils/uncertain_value.h"
 
 namespace mio
 {

cpp/memilio/compartments/flow_model.h

@@ -24,7 +24,7 @@
 #include "memilio/compartments/compartmental_model.h"
 #include "memilio/utils/index_range.h"
 #include "memilio/utils/flow.h"
-#include "memilio/utils/type_list.h"
+#include "memilio/utils/type_list.h" // IWYU pragma: keep for easier flow definitions
 
 namespace mio
 {
@@ -37,8 +37,8 @@ namespace details
 // First a list of tuples is generated for each Tag in Tags, where the tuple is either of type tuple<Tag>, or if
 // Tag == OmittedTag, of type tuple<>. This list is then concatenated, effectively removing OmittedTag.
 template <class OmittedTag, class... Tags>
-decltype(std::tuple_cat(std::declval<typename std::conditional<std::is_same<OmittedTag, Tags>::value, std::tuple<>,
-                                                               std::tuple<Tags>>::type>()...))
+decltype(std::tuple_cat(
+    std::declval<std::conditional_t<std::is_same_v<OmittedTag, Tags>, std::tuple<>, std::tuple<Tags>>>()...))
     filter_tuple(std::tuple<Tags...>);
 
 // Function declaration used to replace type T by std::tuple.
@@ -117,7 +117,7 @@ class FlowModel : public CompartmentalModel<FP, Comp, Pop, Params>
             get_rhs_impl(flows, dydt, Index<>{});
         }
         else {
-            for (FlowIndex I : make_index_range(reduce_index<FlowIndex>(this->populations.size()))) {
+            for (FlowIndex I : reduce_index<FlowIndex>(this->populations.size())) {
                 get_rhs_impl(flows, dydt, I);
             }
         }
@@ -135,7 +135,7 @@ class FlowModel : public CompartmentalModel<FP, Comp, Pop, Params>
      * @param[out] dydt A reference to the calculated output.
      */
     void get_derivatives(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
-                         Eigen::Ref<Eigen::VectorX<FP>> dydt) const override final
+                         Eigen::Ref<Eigen::VectorX<FP>> dydt) const final
     {
         m_flow_values.setZero();
         get_flows(pop, y, t, m_flow_values);
@@ -199,7 +199,7 @@ class FlowModel : public CompartmentalModel<FP, Comp, Pop, Params>
     template <Comp Source, Comp Target>
     constexpr size_t get_flat_flow_index() const
     {
-        static_assert(std::is_same<FlowIndex, Index<>>::value, "Other indices must be specified");
+        static_assert(std::is_same_v<FlowIndex, Index<>>, "Other indices must be specified");
         return index_of_type_v<Flow<Source, Target>, Flows>;
     }
 

cpp/memilio/compartments/parameter_studies.h

@@ -22,9 +22,7 @@
 
 #include "memilio/io/binary_serializer.h"
 #include "memilio/io/io.h"
-#include "memilio/mobility/graph_simulation.h"
 #include "memilio/utils/logging.h"
-#include "memilio/utils/metaprogramming.h"
 #include "memilio/utils/miompi.h"
 #include "memilio/utils/random_number_generator.h"
 #include "memilio/mobility/metapopulation_mobility_instant.h"
@@ -160,14 +158,11 @@ class ParameterStudy
     run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result)
     {
         using ResultT = EnsembleResultT<CreateSimulationFunction, ProcessSimulationResultFunction>;
-        int num_procs, rank;
+        int num_procs = 1, rank = 0;
 
 #ifdef MEMILIO_ENABLE_MPI
         MPI_Comm_size(mpi::get_world(), &num_procs);
         MPI_Comm_rank(mpi::get_world(), &rank);
-#else
-        num_procs = 1;
-        rank      = 0;
 #endif
 
         //The ParameterDistributions used for sampling parameters use thread_local_rng()
@@ -176,9 +171,8 @@ class ParameterStudy
         m_rng.synchronize();
 
         std::vector<size_t> run_distribution = distribute_runs(m_num_runs, num_procs);
-        size_t start_run_idx =
-            std::accumulate(run_distribution.begin(), run_distribution.begin() + size_t(rank), size_t(0));
-        size_t end_run_idx = start_run_idx + run_distribution[size_t(rank)];
+        size_t start_run_idx = std::accumulate(run_distribution.begin(), run_distribution.begin() + rank, size_t(0));
+        size_t end_run_idx   = start_run_idx + run_distribution[rank];
 
         if constexpr (std::is_void_v<ResultT>) {
             // if the processor returns nothing, there is nothing to synchronize

cpp/memilio/compartments/simulation.h

@@ -23,7 +23,6 @@
 #include "memilio/compartments/compartmental_model.h"
 #include "memilio/compartments/simulation_base.h"
 #include "memilio/math/integrator.h"
-#include "memilio/utils/metaprogramming.h"
 #include "memilio/math/stepper_wrapper.h"
 #include "memilio/utils/time_series.h"
 
@@ -75,27 +74,15 @@ class Simulation : public details::SimulationBase<FP, M, OdeIntegrator<FP>>
 };
 
 /**
- * Defines the return type of the `advance` member function of a type.
- * Template is invalid if this member function does not exist.
- *
- * @tparam FP floating point type, e.g., double
- * @tparam Sim a compartment model simulation type.
- */
-template <typename FP, class Sim>
-using advance_expr_t = decltype(std::declval<Sim>().advance(std::declval<FP>()));
-
-/**
- * Template meta function to check if a type is a compartment model simulation.
- * Defines a static constant of name `value`.
- * The constant `value` will be equal to true if Sim is a valid compartment simulation type.
- * Otherwise, `value` will be equal to false.
- * @tparam FP floating point type, e.g., double
- * @tparam Sim a type that may or may not be a compartment model simulation.
+ * Concept to check if a type is a simulation for a compartmental model.
+ * @tparam Simulation A type that may or may not be a compartmental model simulation.
+ * @tparam FP A floating point type, e.g., double.
  */
-template <typename FP, class Sim>
-using is_compartment_model_simulation =
-    std::integral_constant<bool, (is_expression_valid<advance_expr_t, FP, Sim>::value &&
-                                  IsCompartmentalModel<typename Sim::Model, FP>)>;
+template <class Simulation, typename FP>
+concept IsCompartmentalModelSimulation = requires(Simulation simulation, FP t) {
+    requires IsCompartmentalModel<typename Simulation::Model, FP>;
+    simulation.advance(t);
+};
 
 /**
  * @brief Run a Simulation of a CompartmentalModel.

cpp/memilio/compartments/stochastic_model.h

@@ -22,7 +22,6 @@
 
 #include "memilio/compartments/compartmental_model.h"
 #include "memilio/compartments/flow_model.h"
-#include "memilio/utils/metaprogramming.h"
 #include "memilio/utils/random_number_generator.h"
 
 namespace mio

cpp/memilio/data/analyze_result.cpp

@@ -18,10 +18,6 @@
 * limitations under the License.
 */
 #include "memilio/data/analyze_result.h"
-#include "memilio/math/interpolation.h"
-
-#include <algorithm>
-#include <cassert>
 
 namespace mio
 {

cpp/memilio/data/analyze_result.h

@@ -21,6 +21,7 @@
 #define MEMILIO_DATA_ANALYZE_RESULT_H
 
 #include "memilio/config.h"
+#include "memilio/utils/logging.h"
 #include "memilio/utils/time_series.h"
 #include "memilio/mobility/metapopulation_mobility_instant.h"
 #include "memilio/math/interpolation.h"
@@ -376,81 +377,51 @@ template <class FP>
 IOResult<TimeSeries<FP>> merge_time_series(const TimeSeries<FP>& ts1, const TimeSeries<FP>& ts2,
                                            bool add_values = false)
 {
-    TimeSeries<FP> merged_ts(ts1.get_num_elements());
     if (ts1.get_num_elements() != ts2.get_num_elements()) {
         log_error("TimeSeries have a different number of elements.");
         return failure(mio::StatusCode::InvalidValue);
     }
-    else {
-        Eigen::Index t1_iterator = 0;
-        Eigen::Index t2_iterator = 0;
-        bool t1_finished         = false;
-        bool t2_finished         = false;
-        while (!t1_finished || !t2_finished) {
-            if (!t1_finished) {
-                if (ts1.get_time(t1_iterator) < ts2.get_time(t2_iterator) ||
-                    t2_finished) { // Current time point of first TimeSeries is smaller than current time point of second TimeSeries or second TimeSeries has already been copied entirely
-                    merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
-                    t1_iterator += 1;
-                }
-                else if (!t2_finished && ts1.get_time(t1_iterator) ==
-                                             ts2.get_time(t2_iterator)) { // Both TimeSeries have the current time point
-                    if (add_values) {
-                        merged_ts.add_time_point(ts1.get_time(t1_iterator),
-                                                 ts1.get_value(t1_iterator) + ts2.get_value(t2_iterator));
-                    }
-                    else {
-                        merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
-                        log_warning("Both TimeSeries have values for t={}. The value of the first TimeSeries is used",
-                                    ts1.get_time(t1_iterator));
-                    }
-                    t1_iterator += 1;
-                    t2_iterator += 1;
-                    if (t2_iterator >=
-                        ts2.get_num_time_points()) { // Check if all values of second TimeSeries have been copied
-                        t2_finished = true;
-                        t2_iterator = ts2.get_num_time_points() - 1;
-                    }
-                }
-                if (t1_iterator >=
-                    ts1.get_num_time_points()) { // Check if all values of first TimeSeries have been copied
-                    t1_finished = true;
-                    t1_iterator = ts1.get_num_time_points() - 1;
-                }
+    if (!ts1.is_strictly_monotonic() || !ts2.is_strictly_monotonic()) {
+        log_error("TimeSeries need to have strictly monotonic time points to be merged.");
+        return failure(mio::StatusCode::InvalidValue);
+    }
+    Eigen::Index t1_iterator   = 0;
+    Eigen::Index t2_iterator   = 0;
+    const Eigen::Index t1_size = ts1.get_num_time_points();
+    const Eigen::Index t2_size = ts2.get_num_time_points();
+    TimeSeries<FP> merged_ts(ts1.get_num_elements());
+    merged_ts.reserve(t1_size + t2_size);
+    // merge entries of both time series until one finishes
+    while (t1_iterator < t1_size && t2_iterator < t2_size) {
+        // check which ts has the smaller time at the current iterator, and merge it
+        if (ts1.get_time(t1_iterator) < ts2.get_time(t2_iterator)) {
+            merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
+            ++t1_iterator;
+        }
+        else if (ts1.get_time(t1_iterator) == ts2.get_time(t2_iterator)) {
+            merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
+            if (add_values) {
+                merged_ts.get_last_value() += ts2.get_value(t2_iterator);
             }
-            if (!t2_finished) {
-                if (ts2.get_time(t2_iterator) < ts1.get_time(t1_iterator) ||
-                    t1_finished) { // Current time point of second TimeSeries is smaller than current time point of first TimeSeries or first TimeSeries has already been copied entirely
-                    merged_ts.add_time_point(ts2.get_time(t2_iterator), ts2.get_value(t2_iterator));
-                    t2_iterator += 1;
-                }
-                else if (!t1_finished && ts2.get_time(t2_iterator) ==
-                                             ts1.get_time(t1_iterator)) { // Both TimeSeries have the current time point
-                    if (add_values) {
-                        merged_ts.add_time_point(ts1.get_time(t1_iterator),
-                                                 ts1.get_value(t1_iterator) + ts2.get_value(t2_iterator));
-                    }
-                    else {
-                        merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
-                        log_warning("Both TimeSeries have values for t={}. The value of the first TimeSeries is used",
-                                    ts1.get_time(t1_iterator));
-                    }
-                    t1_iterator += 1;
-                    t2_iterator += 1;
-                    if (t1_iterator >=
-                        ts1.get_num_time_points()) { // Check if all values of first TimeSeries have been copied
-                        t1_finished = true;
-                        t1_iterator = ts1.get_num_time_points() - 1;
-                    }
-                }
-                if (t2_iterator >=
-                    ts2.get_num_time_points()) { // Check if all values of second TimeSeries have been copied
-                    t2_finished = true;
-                    t2_iterator = ts2.get_num_time_points() - 1;
-                }
+            else {
+                log_warning("Both TimeSeries have values for t={}. The value of the first TimeSeries is used",
+                            ts1.get_time(t1_iterator));
             }
+            ++t1_iterator;
+            ++t2_iterator;
+        }
+        else { // " > "
+            merged_ts.add_time_point(ts2.get_time(t2_iterator), ts2.get_value(t2_iterator));
+            ++t2_iterator;
         }
     }
+    // append remaining entries. at most one of the following for loops will be executed
+    for (; t1_iterator < t1_size; ++t1_iterator) {
+        merged_ts.add_time_point(ts1.get_time(t1_iterator), ts1.get_value(t1_iterator));
+    }
+    for (; t2_iterator < t2_size; ++t2_iterator) {
+        merged_ts.add_time_point(ts2.get_time(t2_iterator), ts2.get_value(t2_iterator));
+    }
     return success(merged_ts);
 }