Skip to content

1370-2 Add support for reading templated SECIR models from JSON data #1375

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 7 commits into
base: 1370-move-source-to-header
Choose a base branch
from
Open

1370-2 Add support for reading templated SECIR models from JSON data #1375

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
8e9f700
Add support for reading templated SECIR models from JSON
julianlitz Sep 15, 2025
ccddcd1
Safer FP to int conversion
julianlitz Sep 16, 2025
5641fbd
std::round and static_cast
julianlitz Sep 16, 2025
052e31f
Fix FP metaprogramming madness.
julianlitz Sep 26, 2025
31e299a
FP fixes in Graph_Simulation.h
julianlitz Sep 26, 2025
c20d181
Resolve Issue 1382, 1383, 1384
julianlitz Sep 29, 2025
11bb56c
Implement suggestion to merge round and trunc namespace region
julianlitz Oct 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions cpp/examples/ode_secir_parameter_study.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,8 @@ write_single_run_result(const size_t run,
return node.id;
});
auto num_groups = (int)(size_t)graph.nodes()[0].property.get_simulation().get_model().parameters.get_num_groups();
BOOST_OUTCOME_TRY(mio::save_result(all_results, ids, num_groups,
mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5"))));
BOOST_OUTCOME_TRY(mio::save_result<ScalarType>(
all_results, ids, num_groups, mio::path_join(abs_path, ("Results_run" + std::to_string(run) + ".h5"))));

return mio::success();
}
Expand Down
9 changes: 5 additions & 4 deletions cpp/examples/ode_secir_parameter_study_graph.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -322,12 +322,13 @@ int main()
mio::log_info("Directory '{:s}' was created.", results_dir.string());
}

auto county_ids = std::vector<int>{1001, 1002, 1003};
auto save_results_status = save_results(ensemble_results, ensemble_params, county_ids, results_dir, false);
auto pairs_edges = std::vector<std::pair<int, int>>{};
auto county_ids = std::vector<int>{1001, 1002, 1003};
auto save_results_status =
save_results<ScalarType>(ensemble_results, ensemble_params, county_ids, results_dir, false);
auto pairs_edges = std::vector<std::pair<int, int>>{};
for (auto& edge : params_graph.edges()) {
pairs_edges.push_back({county_ids[edge.start_node_idx], county_ids[edge.end_node_idx]});
}
auto save_edges_status = save_edges(ensemble_edges, pairs_edges, "test_results", false, true);
auto save_edges_status = save_edges<ScalarType>(ensemble_edges, pairs_edges, "test_results", false, true);
}
}
4 changes: 2 additions & 2 deletions cpp/examples/ode_secir_read_graph.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ int main(int argc, char** argv)
mio::osecir::set_params_distributions_normal(model, t0, tmax, 0.2);

std::cout << "Reading Mobility File..." << std::flush;
auto read_mobility_result = mio::read_mobility_plain(filename);
auto read_mobility_result = mio::read_mobility_plain<ScalarType>(filename);
if (!read_mobility_result) {
std::cout << read_mobility_result.error().formatted_message() << '\n';
std::cout << "Create the mobility file with MEmilio Epidata's getCommuterMobility.py file." << '\n';
Expand All @@ -135,7 +135,7 @@ int main(int argc, char** argv)
std::cout << "Done" << std::endl;

std::cout << "Writing Json Files..." << std::flush;
auto write_status = mio::write_graph(graph, "graph_parameters");
auto write_status = mio::write_graph<ScalarType>(graph, "graph_parameters");
if (!write_status) {
std::cout << "Error: " << write_status.error().formatted_message();
}
Expand Down
3 changes: 2 additions & 1 deletion cpp/examples/ode_secir_save_results.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,5 +83,6 @@ int main()
std::vector<mio::TimeSeries<ScalarType>> results_from_sim = {result_from_sim, result_from_sim};
std::vector<int> ids = {1, 2};

auto save_result_status = mio::save_result(results_from_sim, ids, (int)(size_t)nb_groups, "test_result.h5");
auto save_result_status =
mio::save_result<ScalarType>(results_from_sim, ids, (int)(size_t)nb_groups, "test_result.h5");
}
68 changes: 68 additions & 0 deletions cpp/memilio/ad/ad.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,74 @@
#include <cmath>
#include <limits>

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

@julianlitz julianlitz Sep 16, 2025 •
edited
Loading

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.

namespace ad
{
namespace internal
{
using std::round;
template <class AD_TAPE_REAL, class DATA_HANDLER_1>
static inline double round(const ad::internal::active_type<AD_TAPE_REAL, DATA_HANDLER_1>& x)
{
return round(x._value());
}
template <class AD_TAPE_REAL, class A1_T1, class A1_T2, class A1_OP>
static inline double round(const ad::internal::binary_intermediate_aa<AD_TAPE_REAL, A1_T1, A1_T2, A1_OP>& x)
{
return round(x._value());
}
template <class AD_TAPE_REAL, class A1_T1, class A1_OP>
static inline double round(const ad::internal::binary_intermediate_ap<AD_TAPE_REAL, A1_T1, A1_OP>& x)
{
return round(x._value());
}
template <class AD_TAPE_REAL, class A1_T2, class A1_OP>
static inline double round(const ad::internal::binary_intermediate_pa<AD_TAPE_REAL, A1_T2, A1_OP>& x)
{
return round(x._value());
}
template <class AD_TAPE_REAL, class A1_T, class A1_OP>
static inline double round(const ad::internal::unary_intermediate<AD_TAPE_REAL, A1_T, A1_OP>& x)
{
return round(x._value());
}
} // namespace internal
} // namespace ad

// Extend automatic differentiation (AD) library to support std::trunc.
namespace ad
{
namespace internal
{
using std::trunc;
template <class AD_TAPE_REAL, class DATA_HANDLER_1>
static inline double trunc(const ad::internal::active_type<AD_TAPE_REAL, DATA_HANDLER_1>& x)
{
return trunc(x._value());
}
template <class AD_TAPE_REAL, class A1_T1, class A1_T2, class A1_OP>
static inline double trunc(const ad::internal::binary_intermediate_aa<AD_TAPE_REAL, A1_T1, A1_T2, A1_OP>& x)
{
return trunc(x._value());
}
template <class AD_TAPE_REAL, class A1_T1, class A1_OP>
static inline double trunc(const ad::internal::binary_intermediate_ap<AD_TAPE_REAL, A1_T1, A1_OP>& x)
{
return trunc(x._value());
}
template <class AD_TAPE_REAL, class A1_T2, class A1_OP>
static inline double trunc(const ad::internal::binary_intermediate_pa<AD_TAPE_REAL, A1_T2, A1_OP>& x)
{
return trunc(x._value());
}
template <class AD_TAPE_REAL, class A1_T, class A1_OP>
static inline double trunc(const ad::internal::unary_intermediate<AD_TAPE_REAL, A1_T, A1_OP>& x)
{
return trunc(x._value());
}
} // namespace internal
} // namespace ad

// Allow std::numeric_limits to work with AD types.
template <class FP, class DataHandler>
struct std::numeric_limits<ad::internal::active_type<FP, DataHandler>> : public numeric_limits<FP> {
Expand Down
4 changes: 3 additions & 1 deletion cpp/memilio/data/analyze_result.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -309,6 +309,8 @@ std::vector<TimeSeries<FP>> ensemble_mean(const std::vector<std::vector<TimeSeri
template <typename FP>
std::vector<TimeSeries<FP>> ensemble_percentile(const std::vector<std::vector<TimeSeries<FP>>>& ensemble_result, FP p)
{
using std::floor;

assert(p > 0.0 && p < 1.0 && "Invalid percentile value.");

auto num_runs = ensemble_result.size();
Expand All @@ -328,7 +330,7 @@ std::vector<TimeSeries<FP>> ensemble_percentile(const std::vector<std::vector<Ti
return run[node][time][elem];
});
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

@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.

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

}
}
}
Expand Down
72 changes: 36 additions & 36 deletions cpp/memilio/io/mobility_io.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,13 +60,13 @@ IOResult<int> count_lines(const std::string& filename);
* where N is the number of regions
* @param filename name of file to be read
*/
template <typename FP = ScalarType>
IOResult<Eigen::MatrixXd> read_mobility_formatted(const std::string& filename)
template <typename FP>
IOResult<Eigen::MatrixX<FP>> read_mobility_formatted(const std::string& filename)
{
BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename));

if (num_lines == 0) {
return success(Eigen::MatrixXd(0, 0));
return success(Eigen::MatrixX<FP>(0, 0));
}

std::fstream file;
Expand All @@ -75,7 +75,7 @@ IOResult<Eigen::MatrixXd> read_mobility_formatted(const std::string& filename)
return failure(StatusCode::FileNotFound, filename);
}

Eigen::MatrixXd txt_matrix(num_lines - 1, 3);
Eigen::MatrixX<FP> txt_matrix(num_lines - 1, 3);
std::vector<int> ids;

try {
Expand All @@ -89,9 +89,9 @@ IOResult<Eigen::MatrixXd> read_mobility_formatted(const std::string& filename)
filename + ":" + std::to_string(linenumber) + ": Not enough entries in line.");
}
ids.push_back(std::stoi(line[2]));
txt_matrix(linenumber - 1, 0) = std::stoi(line[2]);
txt_matrix(linenumber - 1, 1) = std::stoi(line[3]);
txt_matrix(linenumber - 1, 2) = std::stod(line[4]);
txt_matrix(linenumber - 1, 0) = static_cast<FP>(std::stoi(line[2]));
txt_matrix(linenumber - 1, 1) = static_cast<FP>(std::stoi(line[3]));
txt_matrix(linenumber - 1, 2) = static_cast<FP>(std::stod(line[4]));
}
linenumber++;
}
Expand All @@ -104,7 +104,7 @@ IOResult<Eigen::MatrixXd> read_mobility_formatted(const std::string& filename)
std::vector<int>::iterator iter = std::unique(ids.begin(), ids.end());
ids.resize(std::distance(ids.begin(), iter));

Eigen::MatrixXd mobility = Eigen::MatrixXd::Zero(ids.size(), ids.size());
Eigen::MatrixX<FP> mobility = Eigen::MatrixX<FP>::Zero(ids.size(), ids.size());

for (int k = 0; k < num_lines - 1; k++) {
int row_ind = 0;
Expand All @@ -127,13 +127,13 @@ IOResult<Eigen::MatrixXd> read_mobility_formatted(const std::string& filename)
* Matrix, where N is the number of regions
* @param filename name of file to be read
*/
template <typename FP = ScalarType>
IOResult<Eigen::MatrixXd> read_mobility_plain(const std::string& filename)
template <typename FP>
IOResult<Eigen::MatrixX<FP>> read_mobility_plain(const std::string& filename)
{
BOOST_OUTCOME_TRY(auto&& num_lines, count_lines(filename));

if (num_lines == 0) {
return success(Eigen::MatrixXd(0, 0));
return success(Eigen::MatrixX<FP>(0, 0));
}

std::fstream file;
Expand All @@ -142,7 +142,7 @@ IOResult<Eigen::MatrixXd> read_mobility_plain(const std::string& filename)
return failure(StatusCode::FileNotFound, filename);
}

Eigen::MatrixXd mobility(num_lines, num_lines);
Eigen::MatrixX<FP> mobility(num_lines, num_lines);

try {
std::string tp;
Expand All @@ -154,7 +154,7 @@ IOResult<Eigen::MatrixXd> read_mobility_plain(const std::string& filename)
}
Eigen::Index i = static_cast<Eigen::Index>(linenumber);
for (Eigen::Index j = 0; j < static_cast<Eigen::Index>(line.size()); j++) {
mobility(i, j) = std::stod(line[j]);
mobility(i, j) = static_cast<FP>(std::stod(line[j]));
}
linenumber++;
}
Expand Down Expand Up @@ -302,9 +302,9 @@ IOResult<Graph<Model, MobilityParameters<FP>>> read_graph(const std::string& dir
* @param[in] filename Name of the file where the results will be saved.
* @return Any io errors that occur during writing of the files.
*/
template <typename FP = ScalarType>
IOResult<void> save_edges(const std::vector<TimeSeries<ScalarType>>& results,
const std::vector<std::pair<int, int>>& ids, const std::string& filename)
template <typename FP>
IOResult<void> save_edges(const std::vector<TimeSeries<FP>>& results, const std::vector<std::pair<int, int>>& ids,
const std::string& filename)
{
const auto num_edges = results.size();
size_t edge_indx = 0;
Expand Down Expand Up @@ -337,18 +337,18 @@ IOResult<void> save_edges(const std::vector<TimeSeries<ScalarType>>& results,
"Failed to create the 'Time' DataSet in group " + h5group_name +
" in the file: " + filename);

auto values_t = std::vector<ScalarType>(result.get_times().begin(), result.get_times().end());
auto values_t = std::vector<FP>(result.get_times().begin(), result.get_times().end());
MEMILIO_H5_CHECK(H5Dwrite(dset_t.id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, values_t.data()),
StatusCode::UnknownError,
"Failed to write 'Time' data in group " + h5group_name + " in the file: " + filename);

int start_id = ids[edge_indx].first;
auto total = Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>::Zero(
num_timepoints, num_elements)
.eval();
auto total =
Eigen::Matrix<FP, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>::Zero(num_timepoints, num_elements)
.eval();
while (edge_indx < num_edges && ids[edge_indx].first == start_id) {
const auto& result_edge = results[edge_indx];
auto edge_result = Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>::Zero(
auto edge_result = Eigen::Matrix<FP, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>::Zero(
num_timepoints, num_elements)
.eval();
for (Eigen::Index t_idx = 0; t_idx < result_edge.get_num_time_points(); ++t_idx) {
Expand Down Expand Up @@ -412,17 +412,17 @@ IOResult<void> save_edges(const std::vector<TimeSeries<ScalarType>>& results,
* @param[in] save_percentiles [Default: true] Defines if percentiles are written.
* @return Any io errors that occur during writing of the files.
*/
template <typename FP = ScalarType>
IOResult<void> save_edges(const std::vector<std::vector<TimeSeries<ScalarType>>>& ensemble_edges,
template <typename FP>
IOResult<void> save_edges(const std::vector<std::vector<TimeSeries<FP>>>& ensemble_edges,
const std::vector<std::pair<int, int>>& pairs_edges, const fs::path& result_dir,
bool save_single_runs = true, bool save_percentiles = true)
{
//save results and sum of results over nodes
auto ensemble_edges_sum = sum_nodes(ensemble_edges);
if (save_single_runs) {
for (size_t i = 0; i < ensemble_edges_sum.size(); ++i) {
BOOST_OUTCOME_TRY(save_edges(ensemble_edges[i], pairs_edges,
(result_dir / ("Edges_run" + std::to_string(i) + ".h5")).string()));
BOOST_OUTCOME_TRY(save_edges<FP>(ensemble_edges[i], pairs_edges,
(result_dir / ("Edges_run" + std::to_string(i) + ".h5")).string()));
}
}

Expand All @@ -441,17 +441,17 @@ IOResult<void> save_edges(const std::vector<std::vector<TimeSeries<ScalarType>>>

// save percentiles of Edges
{
auto ensemble_edges_p05 = ensemble_percentile(ensemble_edges, 0.05);
auto ensemble_edges_p25 = ensemble_percentile(ensemble_edges, 0.25);
auto ensemble_edges_p50 = ensemble_percentile(ensemble_edges, 0.50);
auto ensemble_edges_p75 = ensemble_percentile(ensemble_edges, 0.75);
auto ensemble_edges_p95 = ensemble_percentile(ensemble_edges, 0.95);

BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p05, pairs_edges, (result_dir_p05 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p25, pairs_edges, (result_dir_p25 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p50, pairs_edges, (result_dir_p50 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p75, pairs_edges, (result_dir_p75 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges(ensemble_edges_p95, pairs_edges, (result_dir_p95 / "Edges.h5").string()));
auto ensemble_edges_p05 = ensemble_percentile<FP>(ensemble_edges, 0.05);
auto ensemble_edges_p25 = ensemble_percentile<FP>(ensemble_edges, 0.25);
auto ensemble_edges_p50 = ensemble_percentile<FP>(ensemble_edges, 0.50);
auto ensemble_edges_p75 = ensemble_percentile<FP>(ensemble_edges, 0.75);
auto ensemble_edges_p95 = ensemble_percentile<FP>(ensemble_edges, 0.95);

BOOST_OUTCOME_TRY(save_edges<FP>(ensemble_edges_p05, pairs_edges, (result_dir_p05 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges<FP>(ensemble_edges_p25, pairs_edges, (result_dir_p25 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges<FP>(ensemble_edges_p50, pairs_edges, (result_dir_p50 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges<FP>(ensemble_edges_p75, pairs_edges, (result_dir_p75 / "Edges.h5").string()));
BOOST_OUTCOME_TRY(save_edges<FP>(ensemble_edges_p95, pairs_edges, (result_dir_p95 / "Edges.h5").string()));
}
}
return success();
Expand Down
23 changes: 11 additions & 12 deletions cpp/memilio/io/parameters_io.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ int get_region_id(const EpiDataEntry& data_entry)
*
* @return An IOResult indicating success or failure.
*/
template <typename FP = ScalarType>
template <typename FP>
IOResult<void> compute_divi_data(const std::vector<DiviEntry>& divi_data, const std::vector<int>& vregion, Date date,
std::vector<FP>& vnum_icu)
{
Expand Down Expand Up @@ -106,12 +106,12 @@ IOResult<void> compute_divi_data(const std::vector<DiviEntry>& divi_data, const
*
* @return An IOResult indicating success or failure.
*/
template <typename FP = ScalarType>
template <typename FP>
IOResult<void> read_divi_data(const std::string& path, const std::vector<int>& vregion, Date date,
std::vector<FP>& vnum_icu)
{
BOOST_OUTCOME_TRY(auto&& divi_data, mio::read_divi_data(path));
return compute_divi_data(divi_data, vregion, date, vnum_icu);
return compute_divi_data<FP>(divi_data, vregion, date, vnum_icu);
}

/**
Expand All @@ -122,12 +122,12 @@ IOResult<void> read_divi_data(const std::string& path, const std::vector<int>& v
* @return An IOResult containing a vector of vectors, where each inner vector represents the population
* distribution across age groups for a specific region, or an error if the function fails.
*/
template <typename FP = ScalarType>
IOResult<std::vector<std::vector<ScalarType>>>
read_population_data(const std::vector<PopulationDataEntry>& population_data, const std::vector<int>& vregion)
template <typename FP>
IOResult<std::vector<std::vector<FP>>> read_population_data(const std::vector<PopulationDataEntry>& population_data,
const std::vector<int>& vregion)
{
std::vector<std::vector<ScalarType>> vnum_population(
vregion.size(), std::vector<ScalarType>(ConfirmedCasesDataEntry::age_group_names.size(), 0.0));
std::vector<std::vector<FP>> vnum_population(vregion.size(),
std::vector<FP>(ConfirmedCasesDataEntry::age_group_names.size(), 0.0));

for (auto&& county_entry : population_data) {
//accumulate population of states or country from population of counties
Expand Down Expand Up @@ -163,12 +163,11 @@ read_population_data(const std::vector<PopulationDataEntry>& population_data, co
* @return An IOResult containing a vector of vectors, where each inner vector represents the population
* distribution across age groups for a specific region, or an error if the function fails.
*/
template <typename FP = ScalarType>
IOResult<std::vector<std::vector<ScalarType>>> read_population_data(const std::string& path,
const std::vector<int>& vregion)
template <typename FP>
IOResult<std::vector<std::vector<FP>>> read_population_data(const std::string& path, const std::vector<int>& vregion)
{
BOOST_OUTCOME_TRY(auto&& population_data, mio::read_population_data(path));
return read_population_data(population_data, vregion);
return read_population_data<FP>(population_data, vregion);
}

} // namespace mio
Expand Down
Loading