-
Notifications
You must be signed in to change notification settings - Fork 23
Expand file tree
/
Copy pathode_secir_parameter_sampling.cpp
More file actions
111 lines (98 loc) · 4.72 KB
/
ode_secir_parameter_sampling.cpp
File metadata and controls
111 lines (98 loc) · 4.72 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Daniel Abele, Martin J. Kuehn
*
* Contact: Martin J. Kuehn <[email protected]>
*
* 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/utils/parameter_distributions.h"
#include "memilio/utils/random_number_generator.h"
#include "ode_secir/parameter_space.h"
#include "ode_secir/model.h"
#include <stdio.h>
// verify manually that the correct distribution is chosen
int main()
{
/*
* Real valued variable element
*/
ScalarType mean = 5;
ScalarType stddev = 1.5;
ScalarType min = 1;
ScalarType max = 10;
// check if constructor works correctly
mio::ParameterDistributionNormal some_parameter(min, max, mean, stddev);
// "some parameter",
// std::make_unique<mio::ParameterDistributionNormal>(mio::ParameterDistributionNormal(min, max, mean, stddev))};
printf("\n N(%.0f,%.0f)-distribution with sampling only in [%.0f,%.0f]", mean, stddev, min, max);
int counter[10] = {0};
for (int i = 0; i < 1000; i++) {
int rounded = (int)(some_parameter.get_sample(mio::thread_local_rng()) - 1);
if (rounded >= 0 && rounded < 10) {
counter[rounded]++;
}
}
ScalarType acc = 0;
for (int i = 0; i < 9; i++) {
acc += (ScalarType)counter[i] / 1000.0;
printf("\n [%d-%d): %.2f %.2f ", i + 1, i + 2, (ScalarType)counter[i] / 1000.0, acc);
}
printf("\n");
// check if constructor works correctly
printf("\n U(%.0f,%.0f)-distribution", min, max);
mio::ParameterDistributionUniform some_other_parameter(1.0, 10.0);
ScalarType counter_unif[10] = {0};
for (int i = 0; i < 1000; i++) {
int rounded = (int)(some_other_parameter.get_sample(mio::thread_local_rng()) - 1);
if (rounded >= 0 && rounded < 10) {
counter_unif[rounded]++;
}
}
acc = 0;
for (int i = 0; i < 9; i++) {
acc += (ScalarType)counter_unif[i] / 1000.0;
printf("\n [%d-%d): %.2f %.2f ", i + 1, i + 2, (ScalarType)counter_unif[i] / 1000.0, acc);
}
/*
* Contact frequency and dampings variable element
*/
mio::osecir::Model<ScalarType> model(3);
auto& params = model.parameters;
mio::AgeGroup nb_groups = params.get_num_groups();
mio::ContactMatrixGroup<ScalarType> cm_group{mio::ContactMatrix<ScalarType>(
Eigen::MatrixX<ScalarType>::Constant((size_t)nb_groups, (size_t)nb_groups, 0.5))};
params.get<mio::osecir::ContactPatterns<ScalarType>>() = cm_group;
params.get<mio::osecir::ContactPatterns<ScalarType>>().get_dampings().push_back(mio::DampingSampling<ScalarType>(
mio::UncertainValue<ScalarType>(0.5), mio::DampingLevel(0), mio::DampingType(0),
mio::SimulationTime<ScalarType>(30.), std::vector<size_t>(1, size_t(0)),
Eigen::VectorX<ScalarType>::Constant(Eigen::Index(nb_groups.get()), 1.0)));
params.get<mio::osecir::ContactPatterns<ScalarType>>().get_dampings()[0].get_value().set_distribution(
mio::ParameterDistributionNormal(0.0, 1.0, 0.5, 0.2));
params.get<mio::osecir::ContactPatterns<ScalarType>>().get_dampings().push_back(mio::DampingSampling<ScalarType>(
mio::UncertainValue<ScalarType>(0.3), mio::DampingLevel(1), mio::DampingType(0),
mio::SimulationTime<ScalarType>(10.), std::vector<size_t>(1, size_t(0)),
Eigen::VectorX<ScalarType>::Constant(Eigen::Index(nb_groups.get()), 1.0)));
params.get<mio::osecir::ContactPatterns<ScalarType>>().get_dampings()[0].get_value().set_distribution(
mio::ParameterDistributionNormal(0.0, 1.0, 0.4, 0.05));
draw_sample(model);
auto& cfmat_sample = params.get<mio::osecir::ContactPatterns<ScalarType>>().get_cont_freq_mat();
printf("\n\n Number of dampings: %zu\n", cfmat_sample[0].get_dampings().size());
//Dampings are sorted automatically by time, therefore the second DampingSamping is at the first position
printf("\n First damping at %.2f with factor %.2f\n", ScalarType(cfmat_sample[0].get_dampings()[0].get_time()),
cfmat_sample[0].get_dampings()[0].get_coeffs()(0, 0));
// printout the second damping
printf("\n Damping at day %.2f\n\t", ScalarType(cfmat_sample[0].get_dampings()[1].get_time()));
std::cout << cfmat_sample[0].get_dampings()[1].get_coeffs() << std::endl;
}