894 Implement dynamic optimization

[hannemann-tamas]

Changes and Information

Please briefly list the changes made, additional Information and what the Reviewer should look out for:

• Implemented ODE-SEAIR model to reproduce optimization results from literature.
• Integrating automatic differentiation (AD) for ODE-SEAIR model.
• Extending template structure for AD data type.
• Implemented dynamic optimization example using IPopt.

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

• [x] Every addressed issue is linked (use the "Closes #ISSUE" keyword below)
• [x] New code adheres to coding guidelines
• [x] No large data files have been added (files should in sum not exceed 100 KB, avoid PDFs, Word docs, etc.)
• [x] Tests are added for new functionality and a local test run was successful
• [x] Appropriate documentation for new functionality has been added (Doxygen in the code and Markdown files if necessary)
• [x] Proper attention to licenses, especially no new third-party software with conflicting license has been added
• [x] (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)
• [x] No large data files added in the whole history of commits(files should in sum not exceed 100 KB, avoid PDFs, Word docs, etc.)

Closes #894

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 98.48649% with 28 lines in your changes are missing coverage. Please review.

Project coverage is 96.51%. Comparing base (4badae5) to head (e37d8e1).

:exclamation: Current head e37d8e1 differs from pull request most recent head e2bfbe6

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

Files	Patch %	Lines
cpp/models/ode_secirvvs/parameter_space.h	95.09%	5 Missing :warning:
cpp/memilio/utils/logging.h	0.00%	4 Missing :warning:
cpp/memilio/epidemiology/uncertain_matrix.h	93.75%	3 Missing :warning:
cpp/models/ode_secir/parameters_io.h	96.93%	3 Missing :warning:
cpp/memilio/compartments/compartmentalmodel.h	50.00%	2 Missing :warning:
cpp/memilio/math/integrator.h	93.33%	2 Missing :warning:
cpp/models/ide_seir/model.h	96.29%	2 Missing :warning:
cpp/models/ode_secirvvs/model.h	99.17%	2 Missing :warning:
cpp/memilio/epidemiology/populations.h	80.00%	1 Missing :warning:
cpp/memilio/utils/uncertain_value.h	94.73%	1 Missing :warning:
... and 3 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #895      +/-   ##
==========================================
+ Coverage   96.42%   96.51%   +0.08%     
==========================================
  Files         130      128       -2     
  Lines       10364    10537     +173     
==========================================
+ Hits         9994    10170     +176     
+ Misses        370      367       -3     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[dabele]

@mknaranja asked me to take a look here. maybe the templates aren't so terrible, since we are already making heavy use of templates. has anyone compared compile times before and after, especially of the unit tests since those are compiled most often by developers. (we haven't really watched compile times very much until now to be fair, but we shouldn't make it explode either).

But a proposal anyway to avoid making everything a template: We had before the configurable ScalarType. It wasn't used consistently, but maybe that's not so hard to fix. What if we gave users the choice (using CMake config) between a memilio library with ScalarType = double and a library with ScalarType = ad::.... The second version could be used for AD. This has it's own draw backs, e.g. if you want to use a different AD type, you need to recompile it. Not sure how much variability there is, how often a user would want to change the type.

If we stick with templates, the best we can hope for is consistency. Since almost everything is going to use FP as a template argument, I think it should always come first. There is not much gained from having a default value, but I guess it's ok to have it if the first is also the last template argument. Also think about where it makes sense to inherit the FP type, e.g.

template<class S> struct A { using Scalar = S; };
template<class T> struct B { using Scalar = typename T::Scalar; }; 

This is mostly if B is already a template that requires a compatible argument, i.e. something with the same FP type. In some cases, this can also be used if B wasn't a template, i.e. template<class A> struct B { A a; } instead of template<class FP> struct B { A<FP> a; }. This may give a cleaner result in some cases. All to make sure that FP is really the same everywhere.

I don't know how to handle the ParameterSet nicely, there's a lot of repetition.

[hannemann-tamas]

But it would be visible as third party and this is explicitly not wanted by the author. I know it is not the clearest solution, however I would stick to it.

Am Fr., 9. Feb. 2024 um 17:58 Uhr schrieb dabele @.***>:

@.**** commented on this pull request.

In cpp/CMakeLists.txt https://github.com/SciCompMod/memilio/pull/895#discussion_r1484580719:

@@ -113,6 +115,7 @@ endif()

add parts of the project

include(thirdparty/CMakeLists.txt) +add_subdirectory(memilio/ad)

to be clear, it would still be in the repository, just a different folder.

— Reply to this email directly, view it on GitHub https://github.com/SciCompMod/memilio/pull/895#discussion_r1484580719, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADM4BCDOFBORRQACCI5QCYTYSZISPAVCNFSM6AAAAABB3IBBL6VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMYTQNZSHAZTGNJZGA . You are receiving this because you authored the thread.Message ID: @.***>

[hannemann-tamas]

@mknaranja asked me to take a look here. maybe the templates aren't so terrible, since we are already making heavy use of templates. has anyone compared compile times before and after, especially of the unit tests since those are compiled most often by developers. (we haven't really watched compile times very much until now to be fair, but we shouldn't make it explode either).

But a proposal anyway to avoid making everything a template: We had before the configurable ScalarType. It wasn't used consistently, but maybe that's not so hard to fix. What if we gave users the choice (using CMake config) between a memilio library with ScalarType = double and a library with ScalarType = ad::.... The second version could be used for AD. This has it's own draw backs, e.g. if you want to use a different AD type, you need to recompile it. Not sure how much variability there is, how often a user would want to change the type.

Using the ScalarType-mechanism would limit the possible application a performance of AD significantly. For the dynamic optimization, I already use the double and ad::... data types in the same program. Restricting only to ad::.... would hurt runtime performance of the optimization. Since optimization requires more resources than simulation, I would like to avoid this. Another application scenario would be the use of an implicit or linear-implicit solver for the ODEs (that requires the system's Jacobian) if the ODEs appear to be stiff. Then. one would certainly also want to have double and ad::... data types in the same program.

[hannemann-tamas]

Unfortunately, I was a bit premature about the AD possibility of the ABM model. It is indeed a problem that at the end of the simulation the result is a vector of natural numbers. This means that derivatives in the conventional sense do not exist. I had imagined a stochastic process with real-valued states, which would then already be differentiable.

For discrete probability distributions, it is possible to calculate a so-called "stochastic derivative". Apparently there are also more extensions of the classical derivative concept to discrete random variables. However, I am not able to implement this in the short term.

It probably makes sense to limit the templatization to the ODE and IDE models. Nevertheless, I would archive my changes to the ABM models in a suitable location. Then one would have a starting point if one wants to implement "stochastic derivatives" after all.

[dabele]

Using the ScalarType-mechanism would limit the possible application a performance of AD significantly. For the dynamic optimization, I already use the double and ad::... data types in the same program. Restricting only to ad::.... would hurt runtime performance of the optimization. Since optimization requires more resources than simulation, I would like to avoid this. Another application scenario would be the use of an implicit or linear-implicit solver for the ODEs (that requires the system's Jacobian) if the ODEs appear to be stiff. Then. one would certainly also want to have double and ad::... data types in the same program

The idea would be to have both available, in the same library or probably better two libraries. The AD enabled library would need it's own namespace. So the normal version would still be available. So it would only be a problem if different AD data types are necessary, which I believe is needed at least for forward and reverse.

But I also think the added templates are much less severe changes if the ABM is excluded as discussed above, so maybe it's not worth coming up with alternative setups.

It probably makes sense to limit the templatization to the ODE and IDE models. Nevertheless, I would archive my changes to the ABM models in a suitable location. Then one would have a starting point if one wants to implement "stochastic derivatives" after all

Makes sense. As far as I can tell, the changes to ABM are almost exclusively to add a template and move the code to the header file. The ABM is constantly changing (and there are big changes imminent), so having the code available is not that useful, we would have to start from scratch. Storing the code in a branch is no effort, of course. But if there were any additional complications (e.g., have you done anything about the openmp parts) maybe it's better to write that down somewhere explicitly.

[reneSchm]

We use the somewhat complicated type Eigen::Matrix<FP, Eigen::Dynamic, 1> regularly, even for e.g. get_flows, which has to be implemented by users. Should we alias it to make it easier to use?

I think something like

namespace mio {

template <class FP = ScalarType>
using Vector = Eigen::Matrix<FP, Eigen::Dynamic, 1>;

}

would work, maybe in utils/eigen.h. Of course, we could also add a template for the vector size (defaulting to dynamic) or, similarly, a Matrix alias. What do you think?

[mknaranja]

@reneSchm I think this is a good idea! Maybe we should only check if Vector is not something that hides another structure behind and take a naming that is less common?

[reneSchm]

@reneSchm I think this is a good idea! Maybe we should only check if Vector is not something that hides another structure behind and take a naming that is less common?

With our current dependencies there is only std::vector or e.g. Eigen::Vector3f with some postfix, but there is no "Eigen::Vector". So the name mio::Vector would be unique, while Matrix would be in both the mio and the Eigen namespace.

Using a common name like Vector does make it easier to confuse types from different namespaces, but since we always use the Eigen and std namespaces explicitly, I don't think that will be a problem. Also, it makes it more clear what the type does, since it is just a specific "Eigen::Vector" (technically a Eigen::Matrix). Naming it something like "Array" would be misleading. But if you are more creative than me, we could of course use something different.

[mknaranja]

No let's go with it then.

[reneSchm]

While aliasing Vector, I looked for uses of ScalarType or double in the files I changed (to check if I used the correct Vector), and noticed that some of the changed files use FP inconsistently:

• populations.h
• adapt_rk.h
• integrator.h
• stepper_wrapper.h
• model.h in ode_secir, ode_secirvvs, ode_seir

We need to decide whether to use ScalarType or FP here. At least I think we don't want double anywhere explicitly? Are there any reasons we want to use something other than FP as a floating point type inside a class or function that has the FP template?

[mknaranja]

@reneSchm I think there are some things like Eigen structures that always explicitly use double while our custom stuff could be built with float or double for instance. This could maybe be a reason. But you have a more holistic view on all the different appearances / usages right now. Please decide on your own. Just let us bring this PR to a merge soon :)

[mknaranja]

@lenaploetzke Can you check if #983 resolves the pygeneration issue here? @reneSchm Can you merge the new main and check that the integrator gets resolved correctly?

[lenaploetzke]

I think #983 has nothing to do with the bug here

[mknaranja]

@MaxBetzDLR @reneSchm I cannot merge this. The CI is still failing. Can you address this tomorrow?

[mknaranja]

@hannemann-tamas Can you check the boxes in the pull request above and provide a bullet point list with the changes you did?

[MaxBetzDLR]

Changed reference file for the memilio-generation test to make the test succeed. Python package needs an update in a followup Issue (#1019)

[reneSchm]

@mknaranja Do we really want FP templates in the SEIR/SIR models? Since there is no working optimization for them, it would be much easier to e.g. not have every model parameter use its own FP template, and instead just use ScalarType directly.

We could also think about setting a model specific type ScalarType/FP instead.

[reneSchm]

No large data files added in the whole history of commits(files should in sum not exceed 100 KB, avoid PDFs, Word docs, etc.)

This is ticked, but IIRC there was an update to the boost archive, that was almost certainly larger than 1MB, at some point before its removal - is this a problem? @mknaranja

[reneSchm]

All right, all tests passed, no open conversations. But there are some things missing from test coverage, most not originating from this PR. Here is a list of the missed lines:

• cpp/models/ode_secirvvs/parameter_space.h
  • 64f: untested log_error and assert
  • 182: unused value for delta_fac
• cpp/memilio/utils/logging.h
  • 136f: The struct fmt::formatter<ad::gt1s<double>::type> is not covered by tests, the other struct appears unused. We should cover this, probably only using some fmt function or the structs directly.
• cpp/memilio/epidemiology/uncertain_matrix.h
  • 82f: operator=, maybe should be covered.
• cpp/models/ode_secir/parameters_io.h
  • 156f: A warning in set_confirmed_cases_data; Reading the warning suggests that we may want to return a failure there instead, which then should definitely be covered.
  • 249: Another warning; similar to above.
• cpp/models/ide_seir/model.h
  • 85f: log_warning
• cpp/memilio/math/integrator.h
  • 142 & 145: log_warnings on step sizing. I think they don't have to be covered.
• cpp/memilio/compartments/compartmentalmodel.h
  • 98f: get_derivatives: If someone knows how to test a virtual undefined function, have at it.
• cpp/memilio/epidemiology/populations.h
  • 257: log_error
• cpp/models/ode_secirvvs/parameters_io.h
  • 451: log_warning; should maybe return a failure?
• cpp/models/sde_sirs/simulation.h
  • 78f: simulate function, in my opinion should not be covered.
• cpp/models/sde_sir/simulation.h
  • 78f: simulate function, in my opinion should not be covered.
• cpp/memilio/utils/uncertain_value.h
  • 233f: PrintTo function used for testing. Maybe we should move it into the test directory?

Other than maybe changes for the tests, I think we are done here @HenrZu @jubicker @lenaploetzke @mknaranja ?

[lenaploetzke]

I think we are done here @HenrZu @jubicker @lenaploetzke @mknaranja ?

Yes, I think so too, the things you mentioned are probably something for a new issue.

[mknaranja]

I hope to have a brief look at this tomorrow. Is there only someone who could resolve the last merge conflicts?

Thank you all for this tremendous work!

[reneSchm]

Is there only someone who could resolve the last merge conflicts?

Hey, those weren't there before! Anyways, I'm on it...

[reneSchm]

Merge is done, thanks @MaxBetzDLR for your help

[mknaranja]

Fantastic! Thanks again for the tremendous work to finally merge this @hannemann-tamas @reneSchm @HenrZu @lenaploetzke @jubicker @MaxBetzDLR