1038 odeintegrator can run indefinitely #1046

reneSchm · 2024-05-30T18:03:36Z

Changes and Information

Please briefly list the changes (main added features, changed items, or corrected bugs) made:

Fix a bug where OdeIntegrator would run indefinitely.
Allow integrators to force a given step size.

If need be, add additional information and what the reviewer should look out for in particular:

Please look through the newly added tests. They impose assumptions on the IntegratorCore::step method we may want to discuss.
- e.g., a forced step currently always counts as successful
- is there maybe something missing?
Please proofread the documentation, especially for IntegratorCore::step
Should we document the unused argument force_step_size of EulerIntegratorCore::step? How?

Merge Request - Guideline Checklist

Please check our git workflow. Use the draft feature if the Pull Request is not yet ready to review.

Checks by code author

Every addressed issue is linked (use the "Closes #ISSUE" keyword below)
New code adheres to coding guidelines
No large data files have been added (files should in sum not exceed 100 KB, avoid PDFs, Word docs, etc.)
Tests are added for new functionality and a local test run was successful (with and without OpenMP)
Appropriate documentation for new functionality has been added (Doxygen in the code and Markdown files if necessary)
Proper attention to licenses, especially no new third-party software with conflicting license has been added
(For ABM development) Checked benchmark results and ran and posted a local test above from before and after development to ensure performance is monitored.

Checks by code reviewer(s)

Corresponding issue(s) is/are linked and addressed
Code is clean of development artifacts (no deactivated or commented code lines, no debugging printouts, etc.)
Appropriate unit tests have been added, CI passes, code coverage and performance is acceptable (did not decrease)
No large data files added in the whole history of commits(files should in sum not exceed 100 KB, avoid PDFs, Word docs, etc.)
On merge, add 2-5 lines with the changes (main added features, changed items, or corrected bugs) to the merge-commit-message. This can be taken from the briefly-list-the-changes above (best case) or the separate commit messages (worst case).

Closes #1038

cpp/memilio/math/euler.h

lenaploetzke · 2024-05-31T06:26:39Z

cpp/tests/test_graph_simulation.cpp

+    auto expected_values_n0 = std::vector<double>{692.0, 43.647065696491225, 95.758901343443526, 159.59403296006528};
    auto actual_values_n0   = std::vector<double>{result_n0[0], result_n0[1], result_n0[2], result_n0[3]};
-    auto expected_values_n1 = std::vector<double>{708.0, 44.063147085799322, 96.485223892060375, 160.45162902214025};
+    auto expected_values_n1 = std::vector<double>{708.0, 44.04685402850226, 96.476851270558342, 160.47629470093941};


Can you explain to me why these values have to be changed?

Yes, it is due to the newly forced steps. Disabling forced step sizing (i.e. overriding the bool with true) does recover the old results.
The graph simulation makes small time steps in which it calls each node simulation, which usually calls advance with dt < t_max - t, hence it makes only one forced time step. There are at least two steps smaller than dt_min, hence the discrepancy.

Actually, there was a problem with the OdeIntegrator. In its advance function, I forced a step whenever making the last step, i.e. dt < t_max - t. However, we do not care whether this is the last step or not, we want to know whether the last step would be smaller than dt_min (so whenever the step with dt = t_max - t would fail to adapt). Only then we actually need to force a smaller step, otherwise we may loose accuracy.

The new logic in advance re-takes a step, which is not optimal, but currently there is no way to get dt_min from advance. Maybe we need another virtual method?

codecov · 2024-05-31T07:23:19Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 96.24%. Comparing base (2994f76) to head (0716e48).

❗ Current head 0716e48 differs from pull request most recent head 59f4f1a

Please upload reports for the commit 59f4f1a to get more accurate results.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1046      +/-   ##
==========================================
+ Coverage   96.23%   96.24%   +0.01%     
==========================================
  Files         128      128              
  Lines       10852    10861       +9     
==========================================
+ Hits        10443    10453      +10     
+ Misses        409      408       -1

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

lenaploetzke · 2024-06-10T09:04:03Z

With these changes, the application where the problem occurred works as expected. Will review the code soon.

lenaploetzke

Thank you for the work! I still need to review the test.

You stated in the description: a forced step currently always counts as successful.
Nevertheless, I still get the warning "Adaptive step sizing failed. Forced at least one integration step of size dt_min." if the last step is forced.
Is this a problem?

lenaploetzke · 2024-06-05T08:58:47Z

cpp/memilio/math/adapt_rk.h

+            }

-        dt = std::min(dt, m_dt_max);
+            dt = std::min(dt, m_dt_max);


Why is this line not part of the if statement? It will only change dt if dt>m_dt_max.

dont we need dt = std::max(dt, m_dt_min); (or usestd::clamp(dt, m_dt_min, m_dt_max)as in the stepper_wrapper) additionally?

Why is this line not part of the if statement?

I thought about it, and the code behavior is not affected by the placement of the min. Personally, I like it more outside such that the if only affects the log message, and no further code behavior. I can also move the min down above the while loop so it is closer to the "max", i.e. the if statement right at the beginning of the while. Maybe this also avoids confusion about the missing max/clamp?

cpp/memilio/math/euler.h

lenaploetzke · 2024-06-10T09:24:58Z

cpp/memilio/math/integrator.h

     *     (If adaptive, the given dt is used as the maximum step size, and must be within [dt_min, dt_max].
     *      During integration, dt is adjusted in [dt_min, dt_max] to have an optimal size for the next step.)
     * @param[out] ytp1 Set to the approximated value of y at time t + dt.
     *     (If adaptive, this time may be smaller, but it is at least t + dt_min, at most t + dt_max.


Should this be adjusted?

I decided against it to not duplicate information, but I think I could just add "If adaptive and force_step_size==false" or something similar to hint at the other parameter.

lenaploetzke · 2024-06-10T09:31:14Z

cpp/memilio/math/integrator.h

+                // recalculate the erroneous last step and overwrite its results
+                step_okay = m_core->step(f, results[i], t, dt, results[i + 1], true);
+            }
+            is_step_sizing_okay &= step_okay;


I think it would be enough to just use one bool.

With only one bool and a forced last step it is difficult to report what the state before the forced step was. I assume that if we discard the step that went beyond tmax, we also want to discard its "step_okay" value. And if that step was not okay (i.e. false), the &= will overwrite the previous state, making it unrecoverable.

Do you agree with this assumption?

lenaploetzke · 2024-06-10T09:35:33Z

cpp/memilio/math/stepper_wrapper.h

+     * @param[in,out] t Current time step t_{k} for some k. Will be set to t_{k+1} in [t_{k} + dt_min, t_{k} + dt].
+     * @param[in,out] dt Current time step size h=dt. Overwritten by an estimated optimal step size for the next step.
+     * @param[out] ytp1 The approximated value of y(t_{k+1}).
+     * @param[in] force_step_size Forces a step with exactly size dt.


Maybe a little more detailed (as in the integrator)?

reneSchm · 2024-06-10T11:04:11Z

Nevertheless, I still get the warning "Adaptive step sizing failed. Forced at least one integration step of size dt_min." if the last step is forced.

The message should only appear if a non-forced step failed to adapt, i.e. the stepper wanted to make a step smaller than dt_min. Are you sure this message is due to the last step?

cpp/tests/test_numericalIntegration.cpp

Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com>

HenrZu

Just a few short comments. Otherwise really great, especially the tests :)

HenrZu · 2024-06-11T07:09:54Z

cpp/memilio/math/adapt_rk.h

        assert(0 <= m_dt_min);
        assert(m_dt_min <= m_dt_max);


It is important that we check this. But only in the debug build. I think it makes sense to issue at least a warning to the user here.

IMO this is a perfect use case for for an assert, since we cannot continue if these conditions fail, but checking this every iteration is needlessly expensive.
And since the asserts are printed in plain text on failure, I don't see the need to add text to it.

HenrZu · 2024-06-11T07:14:03Z

cpp/memilio/math/adapt_rk.h

+            }

-        dt = std::min(dt, m_dt_max);
+            dt = std::min(dt, m_dt_max);


dont we need dt = std::max(dt, m_dt_min); (or usestd::clamp(dt, m_dt_min, m_dt_max)as in the stepper_wrapper) additionally?

HenrZu · 2024-06-11T07:32:26Z

cpp/memilio/math/integrator.h

            m_is_adaptive |= !floating_point_equal(dt, dt_copy);
+
+            // catch case where dt_min of an adaptive stepper is larger than tmax - t
+            if (t > tmax) {


Maybe add a tolerance here? However, it would be hard to set since we have the template FP type

We do have floating_point compare functions, which work at least with those ad types, and can use both relative and absolute tolerances. To be clear, that tolerance is then used as "t >= tmax, and w.r.t. tolerances t != tmax" (see here)

Do you think that tolerance is needed? Should we use the same tolerance as in the for loop?

cpp/memilio/math/euler.h

HenrZu · 2024-06-11T07:39:36Z

cpp/memilio/math/stepper_wrapper.h

+        if ((dt < m_dt_min || dt > m_dt_max) && !force_step_size) {
            mio::log_warning("IntegratorCore: Restricting given step size dt = {} to [{}, {}].", dt, m_dt_min,
                             m_dt_max);
+            dt = std::min(dt, m_dt_max);


same as in adapt_rk.h

I do not use clamp here nor in adapt_rk as a small optimization, as the main loop of both integrators already check and bound to dt_min.

Though the std::min should trip up some AD simulation, don't the tests cover simulations?

HenrZu · 2024-06-11T07:49:56Z

cpp/memilio/math/stepper_wrapper.h

+        // clamp dt (again)
+        // the last adaptive step (successful or not) may have calculated a new step size smaller than m_dt_min,
+        // and when forcing step size, dt may have been larger than m_dt_max
+        dt = std::clamp(dt, m_dt_min, m_dt_max);


isnt this only necessary if !is_dt_valid && !force_step_size ?

No, as mentioned in the comment. Though it is now technically possible for a successful step without forcing to grow larger than dt_max, I will update that.

I also cannot disable the adaptive step sizing (without some additional code or run time cost).

reneSchm · 2024-06-11T15:39:58Z

@HenrZu @lenaploetzke I think I don't like the force_step_size approach anymore, it makes both the OdeIntegrator as well as the IntegratorCore weirdly complicated. I made another version in #1049, could you please check that one out instead?

reneSchm · 2024-06-20T13:56:48Z

Closed in favor of v2.

reneSchm added 2 commits May 23, 2024 13:45

add an option to force a step size with adaptive integrators

2173956

add test for forced step sizing and update integrators

bb413a3

reneSchm added class::bug Bugs found in the software loc::backend This issue concerns the C++ backend implementation. labels May 30, 2024

reneSchm self-assigned this May 30, 2024

reneSchm linked an issue May 30, 2024 that may be closed by this pull request

OdeIntegrator can run indefinitely #1038

Closed

2 tasks

reneSchm requested review from HenrZu and lenaploetzke May 30, 2024 18:03

lenaploetzke reviewed May 31, 2024

View reviewed changes

cpp/memilio/math/euler.h Show resolved Hide resolved

lenaploetzke reviewed May 31, 2024

View reviewed changes

reneSchm added 2 commits May 31, 2024 08:46

try to fix issue with "previous versions" of msvc

b0cba3b

Merge branch 'main' into 1038-odeintegrator-can-run-indefinitely

802550f

reneSchm added 5 commits May 31, 2024 12:38

change when a step is forced, undo changes to consistency test

c1f3bdc

fix condition when to force the step size in OdeIntegrator

782d81c

add documentation for unused force_step_size in euler integrator

9c73129

test whether odeintegrator forces steps correctly

a46a18a

use correct assert_eq in new test

0716e48

lenaploetzke reviewed Jun 10, 2024

View reviewed changes

cpp/tests/test_numericalIntegration.cpp Outdated Show resolved Hide resolved

Update cpp/tests/test_numericalIntegration.cpp

59f4f1a

Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com>

HenrZu requested changes Jun 11, 2024

View reviewed changes

reneSchm closed this Jun 20, 2024

reneSchm deleted the 1038-odeintegrator-can-run-indefinitely branch June 21, 2024 08:34

1038 odeintegrator can run indefinitely #1046

1038 odeintegrator can run indefinitely #1046

Uh oh!

Conversation

reneSchm commented May 30, 2024

Changes and Information

Merge Request - Guideline Checklist

Checks by code author

Checks by code reviewer(s)

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

codecov bot commented May 31, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

lenaploetzke commented Jun 10, 2024

Uh oh!

lenaploetzke left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

reneSchm commented Jun 10, 2024

Uh oh!

Uh oh!

HenrZu left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

reneSchm commented Jun 11, 2024

Uh oh!

reneSchm commented Jun 20, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

codecov bot commented May 31, 2024 •

edited

Loading

lenaploetzke left a comment •

edited

Loading