Skip to content

Conversation

@julianlitz
Copy link
Contributor

@julianlitz julianlitz commented Sep 15, 2025

Changes and Information

This is part 2 of adding support for reading templated SECIR models from JSON data.

In PR #1374, we moved the IO functions from the source file to the header to enable their implementation. This branch builds on that work, so its initial commits are identical to the changes introduced there. Subsequent commits focus on addressing AD compatibility and directly modifying the implementation.

Please briefly list the changes (main added features, modified items, or fixed issues):

  • Use FP instead of double in IO functions
  • Ensure to correctly convert AD types to integers
  • Added ad::value() to convert a FP to ScalarType for the ParameterDistribution
  • Data is read using H5T_NATIVE_DOUBLE and then converted to FP later

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

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 within the code (Doxygen) for new functionality has been added in the code
  • Appropriate external documentation (ReadTheDocs) for new functionality has been added to the online documentation
  • 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 #1370, #1382, #1383, #1384

@codecov
Copy link

codecov bot commented Sep 15, 2025

Codecov Report

❌ Patch coverage is 97.94521% with 9 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.29%. Comparing base (71833e4) to head (11bb56c).

Files with missing lines Patch % Lines
cpp/memilio/io/mobility_io.h 92.59% 2 Missing ⚠️
cpp/models/ode_secir/parameters_io.h 96.00% 2 Missing ⚠️
cpp/models/ode_secirts/parameters_io.h 96.15% 2 Missing ⚠️
cpp/models/ode_secirvvs/parameters_io.h 98.95% 2 Missing ⚠️
...memilio/mobility/metapopulation_mobility_instant.h 83.33% 1 Missing ⚠️
Additional details and impacted files
@@                      Coverage Diff                       @@
##           1370-move-source-to-header    #1375      +/-   ##
==============================================================
- Coverage                       97.31%   97.29%   -0.03%     
==============================================================
  Files                             171      171              
  Lines                           15213    15194      -19     
==============================================================
- Hits                            14805    14783      -22     
- Misses                            408      411       +3     

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

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

#include <cmath>
#include <limits>

// Extend automatic differentiation (AD) library to support std::round.
Copy link
Contributor Author

@julianlitz julianlitz Sep 16, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We extend AD support for std::round and std::trunc. We used the same structure as from std::floor, provided in ad/include/ad.hpp.
Note that we don't modify the original ad/include/ad.hpp library.

});
std::sort(single_element_ensemble.begin(), single_element_ensemble.end());
percentile[node][time][elem] = single_element_ensemble[static_cast<size_t>(num_runs * p)];
percentile[node][time][elem] = single_element_ensemble[static_cast<size_t>(floor(num_runs * p))];
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To cast a FP to size_t, we need to call either floor, ceil, round or trunc.

Comment on lines +289 to +290
tnt_value.set_distribution(
mio::ParameterDistributionUniform(0.8 * ad::value(tnt_capacity), 1.2 * ad::value(tnt_capacity)));
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ParameterDistributionUniform takes in ScalarType as an argument.
To cast FP to ScalarType we need to call ad::value.

Comment on lines +316 to +317
compartment_value.set_distribution(mio::ParameterDistributionUniform(
0.9 * ad::value(compartment_value.value()), 1.1 * ad::value(compartment_value.value())));
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ParameterDistributionUniform takes in ScalarType as an argument.
To cast FP to ScalarType we need to call ad::value.

if (it != vregion.end()) {
auto region_idx = size_t(it - vregion.begin());
if (rki_entry.date == offset_date_by_days(date, int(-delay))) {
if (rki_entry.date == offset_date_by_days(date, int(trunc(-delay)))) {
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

int(-x) seems to be equivalent to int(std::trunc(-x)).
This is needed for AD to convert a FP to int correctly.

round(static_cast<FP>(model[county].parameters.template get<TimeExposed<FP>>()[(AgeGroup)group]))));
t_InfectedNoSymptoms[county].push_back(static_cast<int>(round(static_cast<FP>(
model[county].parameters.template get<TimeInfectedNoSymptoms<FP>>()[(AgeGroup)group] * reduc_t))));
t_InfectedSymptoms[county].push_back(static_cast<int>(round(static_cast<FP>(
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cast UncertainValue<FP> to FP.
Then call round to be able to cast to an integer.
We can't call round directly on an Uncertain Value.

@julianlitz julianlitz self-assigned this Sep 16, 2025
@julianlitz julianlitz added loc::backend This issue concerns the C++ backend implementation. class::improvement Cleanup that doesn't affect functionality labels Sep 16, 2025
@reneSchm reneSchm changed the base branch from main to 1370-move-source-to-header September 16, 2025 09:45
Comment on lines +544 to +545
using test_commuters_expr_t = decltype(test_commuters<FP>(
std::declval<Sim&>(), std::declval<Eigen::Ref<Eigen::VectorX<FP>>>(), std::declval<FP>()));
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Metaprogramming uses incorrect function declaration.
const& must be removed.

}
template <typename FP, class Sim,
std::enable_if_t<is_expression_valid<get_mobility_factors_expr_t, Sim>::value, void*> = nullptr>
std::enable_if_t<is_expression_valid<get_mobility_factors_expr_t, FP, Sim>::value, void*> = nullptr>
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Important to add a FP in the metaprogramming as well. Otherwise the function can't be found.

auto get_mobility_factors(const SimulationNode<FP, Sim>& node, FP t, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
{
return get_mobility_factors<FP, Sim>(node.get_simulation(), t, y);
return get_mobility_factors<FP>(node.get_simulation(), t, y);
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using Sim here seems to cause issues. Better remove it and let the compiler choose it for us.

Copy link
Member

@MaxBetzDLR MaxBetzDLR Sep 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Max: I dont understand how using Sim causes problems at this point. Can u elaborate on that.

Julian: I am not sure why, but using <FP, Sim> instead of in some places connected with the metaprograming, it wont be able to find the overloading function of the specialized model. It is really strange. I think in the apply_mobility function, it is fine to use <FP, Sim> but I wouldn't risk it in the other functions above which transfer to the overloaded function, as it doesn't seem to work properly then.
Better not risk it here and let the compiler decide which template to use as it can be deduced from the arguments anyway.

Copy link
Contributor Author

@julianlitz julianlitz Sep 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure why, but using <FP, Sim> instead of <FP> in some places connected with the metaprograming, it wont be able to find the overloading function of the specialized model. It is really strange. I think in the apply_mobility function, it is fine to use <FP, Sim> but I wouldn't risk it in the other functions above which transfer to the overloaded function, as it doesn't seem to work properly then.
Better not risk it here and let the compiler decide which template to use as it can be deduced from the arguments anyway.

Copy link
Member

@MaxBetzDLR MaxBetzDLR left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changes requested

auto get_mobility_factors(const SimulationNode<FP, Sim>& node, FP t, const Eigen::Ref<const Eigen::VectorX<FP>>& y)
{
return get_mobility_factors<FP, Sim>(node.get_simulation(), t, y);
return get_mobility_factors<FP>(node.get_simulation(), t, y);
Copy link
Member

@MaxBetzDLR MaxBetzDLR Sep 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Max: I dont understand how using Sim causes problems at this point. Can u elaborate on that.

Julian: I am not sure why, but using <FP, Sim> instead of in some places connected with the metaprograming, it wont be able to find the overloading function of the specialized model. It is really strange. I think in the apply_mobility function, it is fine to use <FP, Sim> but I wouldn't risk it in the other functions above which transfer to the overloaded function, as it doesn't seem to work properly then.
Better not risk it here and let the compiler decide which template to use as it can be deduced from the arguments anyway.

std::string path = mio::path_join(pydata_dir_Germany, "county_current_population.json");
const std::vector<int> region{1002};
auto population_data = mio::read_population_data(path, region).value();
auto population_data = mio::read_population_data<double>(path, region).value();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use a consistent type in all tests. Until now it seems, that ScalarType wasnt used and instead double was picked as FP type.

Copy link
Contributor Author

@julianlitz julianlitz Sep 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the moment the tests use a wild mix of doubles and ScalarType.

The plan is to progress in three steps:

  1. Template missing functionalities (Done)
    Here we had a large PR which also removed doubles from the memilio/ and models/ folder.
  2. Template IO Functions (This PR)
  3. Switch from doubles to ScalarType everywhere. These mostly effect the tests and the python binding. (Future PR)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dont think tests are suppose to use ScalarType. The type should be fixed to produce consistent results

Copy link
Contributor Author

@julianlitz julianlitz Sep 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The MEmilio software is currently compiled with a single, fixed floating-point type, which is defined as:

using ScalarType = double;

However, our goal is to support alternative types, such as:

using ScalarType = float;

If we keep double in the tests while compiling the library with a different ScalarType, the build will fail because of type mismatches. Therefore, the tests must also consistently use ScalarType to ensure compatibility across different configurations.

Only the occurrences of double that interact with MEmilio’s numeric types or APIs should be replaced with ScalarType. Plain double usage that is unrelated to MEmilio internals does not necessarily need to be changed.

I think there are very few instances where we use thresholds e.g. 1e-12 which are tuned for doubles. You are right, that these would have to be addressed.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right for variables that expect ScalarType, but all these instances, where you added ScalarType in the tests, are as a template instatiation of FP. And in those cases they SHOULD also work with dissimilar types.
So please change that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don’t think this PR is the right place to fully sort out the floating-point type policy. Let’s wrap this one up focusing only on adding the missing templates in the IO functions.

About the tests: I’m not sure yet if keeping double everywhere will even compile once ScalarType changes (e.g. to float). My gut feeling is that using whatever ScalarType the library was built with is the safer option — but I’d rather test it properly before we decide.

So I suggest:
Finish this PR first (IO templating only and some bug fixes).
Afterwards I’ll check which tests actually need ScalarType and where double still works.
Then we can discuss it again (maybe with René) and agree on a consistent rule.

Does that sound okay?

std::vector<std::vector<double>>& vnum_InfectedSymptoms, std::vector<std::vector<double>>& vnum_InfectedSevere,
std::vector<std::vector<double>>& vnum_icu, std::vector<std::vector<double>>& vnum_death,
std::vector<std::vector<double>>& vnum_rec, const std::vector<std::vector<int>>& vt_Exposed,
std::vector<std::vector<FP>>& vnum_Exposed, std::vector<std::vector<FP>>& vnum_InfectedNoSymptoms,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should those values really be FP or ScalarType? The IO could work a lot without FP and just cast to that type when writing into the model. Similar with input parameters like scaling_factor_inf

Copy link
Contributor Author

@julianlitz julianlitz Sep 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would argue that templating the read_confirmed_cases_data<FP> is the straight forward and an easy way to do it.
Manually casting lots of std::vectors to FP in some intermediate step where I don't know where this should even happen doesn't sound great.

I like the current approach where we read in something from a file and directly convert that value to float, double or AD.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

class::improvement Cleanup that doesn't affect functionality loc::backend This issue concerns the C++ backend implementation.

Projects

3 participants