1486 Modernizations and minor fixes in core C++ library

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.
@@ -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.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);
 }
 

cpp/memilio/epidemiology/adoption_rate.h

@@ -20,8 +20,7 @@
 #ifndef MIO_EPI_ADOPTIONRATE_H
 #define MIO_EPI_ADOPTIONRATE_H
 
-#include "memilio/utils/index.h"
-#include "memilio/config.h"
+#include "memilio/config.h" // IWYU pragma: keep
 #include "memilio/geography/regions.h"
 
 namespace mio

cpp/memilio/epidemiology/contact_matrix.h

@@ -17,14 +17,12 @@
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
-#ifndef EPI_ODE_CONTACT_FREQUENCY_MATRIX_H
-#define EPI_ODE_CONTACT_FREQUENCY_MATRIX_H
+#ifndef MIO_EPI_CONTACT_MATRIX_H
+#define MIO_EPI_CONTACT_MATRIX_H
 
 #include "memilio/epidemiology/damping.h"
-#include "memilio/math/matrix_shape.h"
 #include "memilio/math/math_utils.h"
 #include "memilio/utils/stl_util.h"
-#include "memilio/utils/logging.h"
 
 #include <numeric>
 #include <ostream>
@@ -84,7 +82,8 @@ class DampingMatrixExpression
      * @param shape_args shape arguments.
      * @tparam T shape arguments.
      */
-    template <class... T, class = std::enable_if_t<std::is_constructible<Shape, T...>::value, void>>
+    template <class... T>
+        requires std::is_constructible_v<Shape, T...>
     explicit DampingMatrixExpression(T... shape_args)
         : DampingMatrixExpression(Matrix::Zero(Shape(shape_args...).rows(), Shape(shape_args...).cols()))
     {
@@ -287,7 +286,8 @@ class DampingMatrixExpressionGroup
      * @param num_groups number of groups.
      * @param num_matrices number of matrices.
      */
-    template <class... T, class = std::enable_if_t<std::is_constructible<Shape, T...>::value, int>>
+    template <class... T>
+        requires(std::is_constructible_v<Shape, T...>)
     explicit DampingMatrixExpressionGroup(size_t num_matrices, T... shape_args)
         : m_matrices(num_matrices, value_type{shape_args...})
     {
@@ -559,4 +559,4 @@ class ContactMatrixGroup : public DampingMatrixExpressionGroup<FP, ContactMatrix
 
 } // namespace mio
 
-#endif //EPI_ODE_CONTACT_FREQUENCY_MATRIX_H
+#endif // MIO_EPI_CONTACT_MATRIX_H

cpp/memilio/epidemiology/damping.cpp

@@ -18,13 +18,6 @@
 * limitations under the License.
 */
 #include "memilio/epidemiology/damping.h"
-#include "memilio/utils/stl_util.h"
-#include "memilio/math/smoother.h"
-
-#include <algorithm>
-#include <cassert>
-#include <stdio.h>
-#include <cmath>
 
 namespace mio
 {