1034 ide model simulation does not work under some conditions (#1035)

compute_compartment_from_flows() now works when the TransitionDistribution SusceptibleToExposed is set to a ConstantFunction (or any other function with an invalid maximum support).

Author: lenaploetzke

cpp/memilio/epidemiology/state_age_function.h

@@ -658,8 +658,8 @@ struct ConstantFunction : public StateAgeFunction {
         unused(tol);
         m_support_max = -2.;
 
-        log_error("This function is not suited to be a TransitionDistribution. Do not call in case of StateAgeFunctions"
-                  "of type b); see documentation of StateAgeFunction Base class.");
+        log_error("This function is not suited to be a TransitionDistribution. Do not call in case of "
+                  "StateAgeFunctions of type b); see documentation of StateAgeFunction Base class.");
 
         return m_support_max;
     }

cpp/models/ide_secir/model.cpp

@@ -60,12 +60,17 @@ void Model::compute_compartment_from_flows(ScalarType dt, Eigen::Index idx_Infec
                                            Eigen::Index idx_IncomingFlow, int idx_TransitionDistribution1,
                                            int idx_TransitionDistribution2)
 {
-    ScalarType sum = 0;
-
+    ScalarType sum       = 0;
+    ScalarType calc_time = 0;
     // Determine relevant calculation area and corresponding index.
-    ScalarType calc_time =
-        std::max(parameters.get<TransitionDistributions>()[idx_TransitionDistribution1].get_support_max(dt, m_tol),
-                 parameters.get<TransitionDistributions>()[idx_TransitionDistribution2].get_support_max(dt, m_tol));
+    if ((1 - parameters.get<TransitionProbabilities>()[idx_TransitionDistribution1]) > 0) {
+        calc_time =
+            std::max(parameters.get<TransitionDistributions>()[idx_TransitionDistribution1].get_support_max(dt, m_tol),
+                     parameters.get<TransitionDistributions>()[idx_TransitionDistribution2].get_support_max(dt, m_tol));
+    }
+    else {
+        calc_time = parameters.get<TransitionDistributions>()[idx_TransitionDistribution1].get_support_max(dt, m_tol);
+    }
 
     Eigen::Index calc_time_index = (Eigen::Index)std::ceil(calc_time / dt) - 1;