988 info on sucessful adaptive step sizing is also shown for fixed step size integrators

cpp/memilio/math/floating_point.h

@@ -72,7 +72,7 @@ template <class T>
 bool floating_point_less(T v1, T v2, T abs_tol = 0, T rel_tol = std::numeric_limits<T>::min())
 {
     auto diff = v1 - v2;
-    return diff < -abs_tol || diff < -abs_max(v1, v2) * rel_tol;
+    return diff < -(abs_tol + rel_tol * abs_max(v1, v2));
 }
 
 /**
@@ -110,7 +110,7 @@ bool floating_point_greater(T v1, T v2, T abs_tol = 0, T rel_tol = std::numeric_
 template <class T>
 bool floating_point_less_equal(T v1, T v2, T abs_tol = 0, T rel_tol = std::numeric_limits<T>::min())
 {
-    return !floating_point_less(v2, v1, abs_tol, rel_tol);
+    return floating_point_less(v1, v2, abs_tol, rel_tol) || floating_point_equal(v1, v2, abs_tol, rel_tol);
 }
 
 /**
@@ -129,7 +129,7 @@ bool floating_point_less_equal(T v1, T v2, T abs_tol = 0, T rel_tol = std::numer
 template <class T>
 bool floating_point_greater_equal(T v1, T v2, T abs_tol = 0, T rel_tol = std::numeric_limits<T>::min())
 {
-    return !floating_point_less(v1, v2, abs_tol, rel_tol);
+    return floating_point_greater(v1, v2, abs_tol, rel_tol) || floating_point_equal(v1, v2, abs_tol, rel_tol);
 }
 } // namespace mio
 

cpp/memilio/math/integrator.cpp

@@ -38,10 +38,11 @@ Eigen::Ref<Eigen::VectorXd> OdeIntegrator::advance(const DerivFunction& f, const
 
     bool step_okay = true;
 
+    double dt_cp; // used to check whether step sizing is adaptive
     double dt_restore = 0; // used to restore dt if dt was decreased to reach tmax
     double t          = t0;
-    size_t i          = results.get_num_time_points() - 1;
-    while (std::abs((tmax - t) / (tmax - t0)) > 1e-10) {
+
+    for (size_t i = results.get_num_time_points() - 1; std::abs((tmax - t) / (tmax - t0)) > 1e-10; ++i) {
         //we don't make timesteps too small as the error estimator of an adaptive integrator
         //may not be able to handle it. this is very conservative and maybe unnecessary,
         //but also unlikely to happen. may need to be reevaluated
@@ -50,24 +51,29 @@ Eigen::Ref<Eigen::VectorXd> OdeIntegrator::advance(const DerivFunction& f, const
             dt_restore = dt;
             dt         = tmax - t;
         }
+        dt_cp = dt;
+
         results.add_time_point();
         step_okay &= m_core->step(f, results[i], t, dt, results[i + 1]);
         results.get_last_time() = t;
 
-        ++i;
+        // if dt has been changed (even slighly) by step, register the current m_core as adaptive
+        m_is_adaptive |= !floating_point_equal(dt, dt_cp);
     }
     // if dt was decreased to reach tmax in the last time iteration,
     // we restore it as it is now probably smaller than required for tolerances
     dt = std::max(dt, dt_restore);
 
-    if (!step_okay) {
-        log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
-    }
-    else if (std::abs((tmax - t) / (tmax - t0)) > 1e-14) {
-        log_warning("Last time step too small. Could not reach tmax exactly.");
-    }
-    else {
-        log_info("Adaptive step sizing successful to tolerances.");
+    if (m_is_adaptive) {
+        if (!step_okay) {
+            log_warning("Adaptive step sizing failed. Forcing an integration step of size dt_min.");
+        }
+        else if (std::abs((tmax - t) / (tmax - t0)) > 1e-14) {
+            log_warning("Last time step too small. Could not reach tmax exactly.");
+        }
+        else {
+            log_info("Adaptive step sizing successful to tolerances.");
+        }
     }
 
     return results.get_last_value();

cpp/memilio/math/integrator.h

@@ -78,6 +78,7 @@ class OdeIntegrator
      */
     OdeIntegrator(std::shared_ptr<IntegratorCore> core)
         : m_core(core)
+        , m_is_adaptive(false)
     {
     }
 
@@ -95,11 +96,13 @@ class OdeIntegrator
 
     void set_integrator(std::shared_ptr<IntegratorCore> integrator)
     {
-        m_core = integrator;
+        m_core        = integrator;
+        m_is_adaptive = false;
     }
 
 private:
     std::shared_ptr<IntegratorCore> m_core;
+    bool m_is_adaptive;
 };
 
 } // namespace mio

cpp/tests/CMakeLists.txt

@@ -37,6 +37,7 @@ set(TESTSOURCES
     test_abm_simulation.cpp
     test_abm_testing_strategy.cpp
     test_abm_world.cpp
+    test_math_floating_point.cpp
     test_analyze_result.cpp
     test_contact_matrix.cpp
     test_type_safe.cpp

cpp/tests/test_math_floating_point.cpp

@@ -0,0 +1,169 @@
+/* 
+* Copyright (C) 2020-2024 MEmilio
+*
+* Authors: Rene Schmieding
+*
+* 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.
+*/
+#include "memilio/math/floating_point.h"
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <string>
+
+template <class FP>
+class TestMathFloatingPoint : public ::testing::Test
+{
+public:
+    using ParamType      = std::array<FP, 5>;
+    using TruthTableType = std::array<std::array<bool, 3>, 4>;
+
+    FP a = FP(1.1), b = FP(3.2); // arbitrary values s.t. a < b
+    const ParamType params = {
+        b - a, // abs_tol
+        (b - a) / b, // rel_tol
+        FP(0.0), // zero
+        FP(10.0), // increased
+        FP(0.1) // reduced
+    };
+};
+
+using FpTypes = ::testing::Types<float, double>;
+
+TYPED_TEST_SUITE(TestMathFloatingPoint, FpTypes);
+
+TYPED_TEST(TestMathFloatingPoint, abs_max)
+{
+    TypeParam v1 = this->a, v2 = this->b;
+    // check for maximum
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(v1, v2));
+    // check symmetry and signs
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(v2, v1));
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(v1, -v2));
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(-v2, v1));
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(-v1, v2));
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(v2, -v1));
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(-v1, -v2));
+    EXPECT_DOUBLE_EQ(v2, mio::abs_max(-v2, -v1));
+}
+
+/**
+ * Check floating point comparisions of v1 and v2 against a truth table.
+ * The table contents are described in the tests below.
+ */
+template <class FP, class Func>
+void test_fp_compare(FP v1, FP v2, typename TestMathFloatingPoint<FP>::ParamType params, Func fp_compare,
+                     typename TestMathFloatingPoint<FP>::TruthTableType truth_table)
+{
+    ASSERT_LT(v1, v2) << "This test is not set up correctly!";
+    const auto [abs_tol, rel_tol, zero, increased, reduced] = params;
+    // on error, print out the name of the test, which is the same as the tested function
+    auto info = std::string("  With fp_compare as ") + testing::UnitTest::GetInstance()->current_test_info()->name();
+    // check basics
+    EXPECT_EQ(fp_compare(v1, v1, zero, zero), truth_table[0][0]) << info;
+    EXPECT_EQ(fp_compare(v1, v2, zero, zero), truth_table[0][1]) << info;
+    EXPECT_EQ(fp_compare(v2, v1, zero, zero), truth_table[0][2]) << info;
+    // check equalities with absolute tolerances
+    EXPECT_EQ(fp_compare(v1, v1, abs_tol, zero), truth_table[1][0]) << info;
+    EXPECT_EQ(fp_compare(v1, v1, reduced * abs_tol, zero), truth_table[1][1]) << info;
+    EXPECT_EQ(fp_compare(v1, v1, increased * abs_tol, zero), truth_table[1][2]) << info;
+    // check equalities with relative tolerances
+    EXPECT_EQ(fp_compare(v1, v1, zero, rel_tol), truth_table[1][0]) << info;
+    EXPECT_EQ(fp_compare(v1, v1, zero, reduced * rel_tol), truth_table[1][1]) << info;
+    EXPECT_EQ(fp_compare(v1, v1, zero, increased * rel_tol), truth_table[1][2]) << info;
+    // check absolute tolerances
+    EXPECT_EQ(fp_compare(v1, v2, abs_tol, zero), truth_table[2][0]) << info;
+    EXPECT_EQ(fp_compare(v1, v2, reduced * abs_tol, zero), truth_table[2][1]) << info;
+    EXPECT_EQ(fp_compare(v1, v2, increased * abs_tol, zero), truth_table[2][2]) << info;
+    // check relative tolerances
+    EXPECT_EQ(fp_compare(v1, v2, zero, rel_tol), truth_table[2][0]) << info;
+    EXPECT_EQ(fp_compare(v1, v2, zero, reduced * rel_tol), truth_table[2][1]) << info;
+    EXPECT_EQ(fp_compare(v1, v2, zero, increased * rel_tol), truth_table[2][2]) << info;
+    // check absolute tolerances, with switched FPs
+    EXPECT_EQ(fp_compare(v2, v1, abs_tol, zero), truth_table[3][0]) << info;
+    EXPECT_EQ(fp_compare(v2, v1, reduced * abs_tol, zero), truth_table[3][1]) << info;
+    EXPECT_EQ(fp_compare(v2, v1, increased * abs_tol, zero), truth_table[3][2]) << info;
+    // check relative tolerances, with switched FPs
+    EXPECT_EQ(fp_compare(v2, v1, zero, rel_tol), truth_table[3][0]) << info;
+    EXPECT_EQ(fp_compare(v2, v1, zero, reduced * rel_tol), truth_table[3][1]) << info;
+    EXPECT_EQ(fp_compare(v2, v1, zero, increased * rel_tol), truth_table[3][2]) << info;
+}
+
+TYPED_TEST(TestMathFloatingPoint, floating_point_equal)
+{
+    typename TestMathFloatingPoint<TypeParam>::TruthTableType truth_table = {{
+        // {equal, strict less, strict greater}, both tolerances zero
+        {true, false, false}, // basic fp operations without tolerances
+        // {exact tolerance, reduced tolerance, increased tolerance} using either abs_tol or rel_tol
+        {true, true, true}, // equality
+        {true, false, true}, // less (or equal)
+        {true, false, true} // greater (or equal)
+    }};
+
+    test_fp_compare(this->a, this->b, this->params, &mio::floating_point_equal<TypeParam>, truth_table);
+}
+
+TYPED_TEST(TestMathFloatingPoint, floating_point_less)
+{
+    typename TestMathFloatingPoint<TypeParam>::TruthTableType truth_table = {{
+        // {equal, strict less, strict greater}, both tolerances zero
+        {false, true, false}, // basic fp operations without tolerances
+        // {exact tolerance, reduced tolerance, increased tolerance} using either abs_tol or rel_tol
+        {false, false, false}, // equality
+        {false, true, false}, // less (or equal)
+        {false, false, false} // greater (or equal)
+    }};
+    test_fp_compare(this->a, this->b, this->params, &mio::floating_point_less<TypeParam>, truth_table);
+}
+
+TYPED_TEST(TestMathFloatingPoint, floating_point_less_equal)
+{
+    typename TestMathFloatingPoint<TypeParam>::TruthTableType truth_table = {{
+        // {equal, strict less, strict greater}, both tolerances zero
+        {true, true, false}, // basic fp operations without tolerances
+        // {exact tolerance, reduced tolerance, increased tolerance} using either abs_tol or rel_tol
+        {true, true, true}, // equality
+        {true, true, true}, // less (or equal)
+        {true, false, true} // greater (or equal)
+    }};
+    test_fp_compare(this->a, this->b, this->params, &mio::floating_point_less_equal<TypeParam>, truth_table);
+}
+
+TYPED_TEST(TestMathFloatingPoint, floating_point_greater)
+{
+    typename TestMathFloatingPoint<TypeParam>::TruthTableType truth_table = {{
+        // {equal, strict less, strict greater}, both tolerances zero
+        {false, false, true}, // basic fp operations without tolerances
+        // {exact tolerance, reduced tolerance, increased tolerance} using either abs_tol or rel_tol
+        {false, false, false}, // equality
+        {false, false, false}, // less (or equal)
+        {false, true, false} // greater (or equal)
+    }};
+    test_fp_compare(this->a, this->b, this->params, &mio::floating_point_greater<TypeParam>, truth_table);
+}
+
+TYPED_TEST(TestMathFloatingPoint, floating_point_greater_equal)
+{
+    typename TestMathFloatingPoint<TypeParam>::TruthTableType truth_table = {{
+        // {equal, strict less, strict greater}, both tolerances zero
+        {true, false, true}, // basic fp operations without tolerances
+        // {exact tolerance, reduced tolerance, increased tolerance} using either abs_tol or rel_tol
+        {true, true, true}, // equality
+        {true, false, true}, // less (or equal)
+        {true, true, true} // greater (or equal)
+    }};
+    test_fp_compare(this->a, this->b, this->params, &mio::floating_point_greater_equal<TypeParam>, truth_table);
+}