Add support for reading in templated SECIR/TS/VVS models

cpp/examples/ad_odeint_example.cpp

@@ -18,7 +18,7 @@
 * limitations under the License.
 */
 
-#include "ad/ad.hpp"
+#include "memilio/ad/ad.h"
 #include "boost/numeric/odeint.hpp"
 #include <array>
 
@@ -72,10 +72,10 @@ int main()
 
     // We want to compare AD derivatives with difference quotient
     // To this end, we simulate again with a perturbation of the initial value of x[0]
-    const double h        = 1e-3; // pertubation for finite differences
+    const double h                  = 1e-3; // pertubation for finite differences
     std::array<double, dimension> y = {ad::value(x[0]), ad::value(x[1])};
-    x[0]                  = 1.0 + h; // add perturbation to initial value of x[0]
-    x[1]                  = 0.0;
+    x[0]                            = 1.0 + h; // add perturbation to initial value of x[0]
+    x[1]                            = 0.0;
 
     // integrate perturbed system
     boost::numeric::odeint::integrate_adaptive(make_controlled<error_stepper_type>(abs_tol, rel_tol),

cpp/examples/ad_square_example.cpp

@@ -18,7 +18,7 @@
 * limitations under the License.
 */
 
-#include "ad/ad.hpp"
+#include "memilio/ad/ad.h"
 #include <iostream>
 
 /*  This example computes the derivative of f(x) = x^2 in two different ways:

cpp/examples/ode_seair.cpp

@@ -23,7 +23,7 @@
 // A detailed description of the model can be found in the publication
 // Tsay et al. (2020), Modeling, state estimation, and optimal control for the US COVID-19 outbreak.
 
-#include "ad/ad.hpp"
+#include "memilio/ad/ad.h"
 
 #include "ode_seair/model.h"
 #include "ode_seair/infection_state.h"

cpp/examples/ode_seair_optimization.cpp

@@ -21,7 +21,7 @@
 * is extracted from the Ipopt documentation
 */
 
-#include "ad/ad.hpp"
+#include "memilio/ad/ad.h"
 
 #include "memilio/utils/compiler_diagnostics.h"
 #include "ode_seair/model.h"
@@ -47,13 +47,13 @@
 class Seair_NLP : public Ipopt::TNLP
 {
 public:
-    static constexpr double N   = 327167434; // total US population
-    Seair_NLP()                 = default;
-    Seair_NLP(const Seair_NLP&) = delete;
-    Seair_NLP(Seair_NLP&&)      = delete;
+    static constexpr double N              = 327167434; // total US population
+    Seair_NLP()                            = default;
+    Seair_NLP(const Seair_NLP&)            = delete;
+    Seair_NLP(Seair_NLP&&)                 = delete;
     Seair_NLP& operator=(const Seair_NLP&) = delete;
-    Seair_NLP& operator=(Seair_NLP&&) = delete;
-    ~Seair_NLP()                      = default;
+    Seair_NLP& operator=(Seair_NLP&&)      = delete;
+    ~Seair_NLP()                           = default;
 
     /** Method to request the initial information about the problem.
     *

cpp/memilio/CMakeLists.txt

@@ -115,6 +115,7 @@ add_library(memilio
     utils/string_literal.h
     utils/type_list.h
     utils/base_dir.h
+    ad/ad.h
 )
 
 target_include_directories(memilio PUBLIC

cpp/memilio/ad/ad.h

@@ -0,0 +1,126 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Julian Litz
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+#ifndef MIO_AD_H
+#define MIO_AD_H
+
+#include "ad/ad.hpp"
+#include "memilio/math/eigen.h"
+
+#include <cmath>
+#include <limits>
+
+// Extend automatic differentiation (AD) library to support std::round.
+namespace ad
+{
+namespace internal
+{
+using std::round;
+template <class AD_TAPE_REAL, class DATA_HANDLER_1>
+static inline double round(const ad::internal::active_type<AD_TAPE_REAL, DATA_HANDLER_1>& x)
+{
+    return round(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T1, class A1_T2, class A1_OP>
+static inline double round(const ad::internal::binary_intermediate_aa<AD_TAPE_REAL, A1_T1, A1_T2, A1_OP>& x)
+{
+    return round(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T1, class A1_OP>
+static inline double round(const ad::internal::binary_intermediate_ap<AD_TAPE_REAL, A1_T1, A1_OP>& x)
+{
+    return round(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T2, class A1_OP>
+static inline double round(const ad::internal::binary_intermediate_pa<AD_TAPE_REAL, A1_T2, A1_OP>& x)
+{
+    return round(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T, class A1_OP>
+static inline double round(const ad::internal::unary_intermediate<AD_TAPE_REAL, A1_T, A1_OP>& x)
+{
+    return round(x._value());
+}
+} // namespace internal
+} // namespace ad
+
+// Extend automatic differentiation (AD) library to support std::trunc.
+namespace ad
+{
+namespace internal
+{
+using std::trunc;
+template <class AD_TAPE_REAL, class DATA_HANDLER_1>
+static inline double trunc(const ad::internal::active_type<AD_TAPE_REAL, DATA_HANDLER_1>& x)
+{
+    return trunc(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T1, class A1_T2, class A1_OP>
+static inline double trunc(const ad::internal::binary_intermediate_aa<AD_TAPE_REAL, A1_T1, A1_T2, A1_OP>& x)
+{
+    return trunc(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T1, class A1_OP>
+static inline double trunc(const ad::internal::binary_intermediate_ap<AD_TAPE_REAL, A1_T1, A1_OP>& x)
+{
+    return trunc(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T2, class A1_OP>
+static inline double trunc(const ad::internal::binary_intermediate_pa<AD_TAPE_REAL, A1_T2, A1_OP>& x)
+{
+    return trunc(x._value());
+}
+template <class AD_TAPE_REAL, class A1_T, class A1_OP>
+static inline double trunc(const ad::internal::unary_intermediate<AD_TAPE_REAL, A1_T, A1_OP>& x)
+{
+    return trunc(x._value());
+}
+} // 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> {
+};
+
+// Ensures that Eigen recognizes your AD types as valid scalars.
+namespace Eigen
+{
+template <class FP, class DataHandler>
+struct NumTraits<ad::internal::active_type<FP, DataHandler>>
+    : GenericNumTraits<ad::internal::active_type<FP, DataHandler>> {
+    using Scalar     = ad::internal::active_type<FP, DataHandler>;
+    using Real       = Scalar;
+    using NonInteger = Scalar;
+    using Nested     = Scalar;
+    enum
+    {
+        IsComplex             = 0,
+        IsInteger             = 0,
+        IsSigned              = 1,
+        RequireInitialization = 1,
+        ReadCost              = 1,
+        AddCost               = 3,
+        MulCost               = 3
+    };
+};
+} // namespace Eigen
+
+#endif // MIO_AD_H

cpp/memilio/ad/include/ad/ad.hpp

@@ -4162,4 +4162,3 @@ namespace ad {
 }
 
 #endif 
-

cpp/memilio/data/analyze_result.cpp

@@ -1,4 +1,4 @@
-/* 
+/*
 * Copyright (C) 2020-2025 MEmilio
 *
 * Authors: Wadim Koslow, Daniel Abele, David Kerkmann, Sascha Korf
@@ -25,177 +25,4 @@
 
 namespace mio
 {
-TimeSeries<double> interpolate_simulation_result(const TimeSeries<double>& simulation_result, const double abs_tol)
-{
-    const auto t0    = simulation_result.get_time(0);
-    const auto t_max = simulation_result.get_last_time();
-    // add another day if the first time point is equal to day_0 up to absolute tolerance tol
-    const auto day0 = (t0 - abs_tol < std::ceil(t0) - 1) ? std::floor(t0) : std::ceil(t0);
-    // add another day if the last time point is equal to day_max up to absolute tolerance tol
-    const auto day_max = (t_max + abs_tol > std::floor(t_max) + 1) ? std::ceil(t_max) : std::floor(t_max);
-
-    // create interpolation_times vector with all days between day0 and day_max
-    std::vector<double> tps(static_cast<int>(day_max) - static_cast<int>(day0) + 1);
-    std::iota(tps.begin(), tps.end(), day0);
-
-    return interpolate_simulation_result(simulation_result, tps);
-}
-
-TimeSeries<double> interpolate_simulation_result(const TimeSeries<double>& simulation_result,
-                                                 const std::vector<double>& interpolation_times)
-{
-    assert(simulation_result.get_num_time_points() > 0 && "TimeSeries must not be empty.");
-
-    assert(std::is_sorted(interpolation_times.begin(), interpolation_times.end()) &&
-           "Time points for interpolation have to be sorted in non-descending order.");
-
-    if (interpolation_times.size() >= 2) {
-        assert((interpolation_times[1] > simulation_result.get_time(0) &&
-                interpolation_times.rbegin()[1] <= simulation_result.get_last_time()) &&
-               "All but the first and the last time point of interpolation have lie between simulation times (strictly "
-               "for lower boundary).");
-    }
-
-    TimeSeries<double> interpolated(simulation_result.get_num_elements());
-
-    if (interpolation_times.size() == 0) {
-        return interpolated;
-    }
-
-    size_t interp_idx = 0;
-    // add first time point of interpolation times in case it is smaller than the first time point of simulation_result
-    // this is used for the case that it equals the first time point of simulation up to tolerance
-    // this is necessary even if the tolerance is 0 due to the way the comparison in the loop is implemented (< and >=)
-    if (simulation_result.get_time(0) >= interpolation_times[0]) {
-        interpolated.add_time_point(interpolation_times[0], simulation_result[0]);
-        ++interp_idx;
-    }
-
-    //interpolate between pair of time points that lie on either side of each interpolation point
-    for (Eigen::Index sim_idx = 0;
-         sim_idx < simulation_result.get_num_time_points() - 1 && interp_idx < interpolation_times.size();) {
-        //only go to next pair of time points if no time point is added.
-        //otherwise check the same time points again
-        //in case there is more than one interpolation point between the two time points
-        if (simulation_result.get_time(sim_idx) < interpolation_times[interp_idx] &&
-            simulation_result.get_time(sim_idx + 1) >= interpolation_times[interp_idx]) {
-            interpolated.add_time_point(
-                interpolation_times[interp_idx],
-                linear_interpolation(interpolation_times[interp_idx], simulation_result.get_time(sim_idx),
-                                     simulation_result.get_time(sim_idx + 1), simulation_result[sim_idx],
-                                     simulation_result[sim_idx + 1]));
-            ++interp_idx;
-        }
-        else {
-            ++sim_idx;
-        }
-    }
-
-    // add last time point of interpolation times in case it is larger than the last time point of simulation_result
-    // this is used for the case that it equals the last time point of simulation up to tolerance
-    if (interp_idx < interpolation_times.size() &&
-        simulation_result.get_last_time() < interpolation_times[interp_idx]) {
-        interpolated.add_time_point(interpolation_times[interp_idx], simulation_result.get_last_value());
-    }
-
-    return interpolated;
-}
-
-std::vector<std::vector<TimeSeries<double>>>
-sum_nodes(const std::vector<std::vector<TimeSeries<double>>>& ensemble_result)
-{
-    auto num_runs        = ensemble_result.size();
-    auto num_nodes       = ensemble_result[0].size();
-    auto num_time_points = ensemble_result[0][0].get_num_time_points();
-    auto num_elements    = ensemble_result[0][0].get_num_elements();
-
-    std::vector<std::vector<TimeSeries<double>>> sum_result(
-        num_runs, std::vector<TimeSeries<double>>(1, TimeSeries<double>::zero(num_time_points, num_elements)));
-
-    for (size_t run = 0; run < num_runs; run++) {
-        for (Eigen::Index time = 0; time < num_time_points; time++) {
-            sum_result[run][0].get_time(time) = ensemble_result[run][0].get_time(time);
-            for (size_t node = 0; node < num_nodes; node++) {
-                sum_result[run][0][time] += ensemble_result[run][node][time];
-            }
-        }
-    }
-    return sum_result;
-}
-
-std::vector<TimeSeries<double>> ensemble_mean(const std::vector<std::vector<TimeSeries<double>>>& ensemble_result)
-{
-    auto num_runs        = ensemble_result.size();
-    auto num_nodes       = ensemble_result[0].size();
-    auto num_time_points = ensemble_result[0][0].get_num_time_points();
-    auto num_elements    = ensemble_result[0][0].get_num_elements();
-
-    std::vector<TimeSeries<double>> mean(num_nodes, TimeSeries<double>::zero(num_time_points, num_elements));
-
-    for (size_t run = 0; run < num_runs; run++) {
-        assert(ensemble_result[run].size() == num_nodes && "ensemble results not uniform.");
-        for (size_t node = 0; node < num_nodes; node++) {
-            assert(ensemble_result[run][node].get_num_time_points() == num_time_points &&
-                   "ensemble results not uniform.");
-            for (Eigen::Index time = 0; time < num_time_points; time++) {
-                assert(ensemble_result[run][node].get_num_elements() == num_elements &&
-                       "ensemble results not uniform.");
-                mean[node].get_time(time) = ensemble_result[run][node].get_time(time);
-                mean[node][time] += ensemble_result[run][node][time] / num_runs;
-            }
-        }
-    }
-
-    return mean;
-}
-
-std::vector<TimeSeries<double>> ensemble_percentile(const std::vector<std::vector<TimeSeries<double>>>& ensemble_result,
-                                                    double p)
-{
-    assert(p > 0.0 && p < 1.0 && "Invalid percentile value.");
-
-    auto num_runs        = ensemble_result.size();
-    auto num_nodes       = ensemble_result[0].size();
-    auto num_time_points = ensemble_result[0][0].get_num_time_points();
-    auto num_elements    = ensemble_result[0][0].get_num_elements();
-
-    std::vector<TimeSeries<double>> percentile(num_nodes, TimeSeries<double>::zero(num_time_points, num_elements));
-
-    std::vector<double> single_element_ensemble(num_runs); //reused for each element
-    for (size_t node = 0; node < num_nodes; node++) {
-        for (Eigen::Index time = 0; time < num_time_points; time++) {
-            percentile[node].get_time(time) = ensemble_result[0][node].get_time(time);
-            for (Eigen::Index elem = 0; elem < num_elements; elem++) {
-                std::transform(ensemble_result.begin(), ensemble_result.end(), single_element_ensemble.begin(),
-                               [=](auto& run) {
-                                   return run[node][time][elem];
-                               });
-                std::sort(single_element_ensemble.begin(), single_element_ensemble.end());
-                percentile[node][time][elem] = single_element_ensemble[static_cast<size_t>(num_runs * p)];
-            }
-        }
-    }
-    return percentile;
-}
-
-double result_distance_2norm(const std::vector<mio::TimeSeries<double>>& result1,
-                             const std::vector<mio::TimeSeries<double>>& result2)
-{
-    assert(result1.size() == result2.size());
-    assert(result1.size() > 0);
-    assert(result1[0].get_num_time_points() > 0);
-    assert(result1[0].get_num_elements() > 0);
-
-    auto norm_sqr = 0.0;
-    for (auto iter_node1 = result1.begin(), iter_node2 = result2.begin(); iter_node1 < result1.end();
-         ++iter_node1, ++iter_node2) {
-        for (Eigen::Index time_idx = 0; time_idx < iter_node1->get_num_time_points(); ++time_idx) {
-            auto v1 = (*iter_node1)[time_idx];
-            auto v2 = (*iter_node2)[time_idx];
-            norm_sqr += ((v1 - v2).array() * (v1 - v2).array()).sum();
-        }
-    }
-    return std::sqrt(norm_sqr);
-}
-
 } // namespace mio