Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "MorrisReporter.h"
11 :
12 : #include "MorrisSampler.h"
13 :
14 : registerMooseObject("StochasticToolsApp", MorrisReporter);
15 :
16 : InputParameters
17 128 : MorrisReporter::validParams()
18 : {
19 128 : InputParameters params = GeneralReporter::validParams();
20 128 : params.addClassDescription("Compute global sensitivities using the Morris method.");
21 :
22 256 : params.addRequiredParam<SamplerName>("sampler", "Morris sampler used to generate samples.");
23 256 : params.addParam<std::vector<VectorPostprocessorName>>(
24 : "vectorpostprocessors",
25 : "List of VectorPostprocessor(s) to utilized for statistic computations.");
26 256 : params.addParam<std::vector<ReporterName>>(
27 : "reporters", {}, "List of Reporter values to utilized for statistic computations.");
28 :
29 : // Confidence Levels
30 256 : params.addParam<std::vector<Real>>(
31 : "ci_levels",
32 128 : std::vector<Real>(),
33 : "A vector of confidence levels to consider, values must be in (0, 1).");
34 256 : params.addParam<unsigned int>("ci_replicates",
35 256 : 10000,
36 : "The number of replicates to use when computing confidence level "
37 : "intervals. This is basically the number of times the statistics "
38 : "are recomputed with a random selection of indices.");
39 256 : params.addParam<unsigned int>("ci_seed",
40 256 : 1,
41 : "The random number generator seed used for creating replicates "
42 : "while computing confidence level intervals.");
43 128 : return params;
44 0 : }
45 :
46 64 : MorrisReporter::MorrisReporter(const InputParameters & parameters)
47 : : GeneralReporter(parameters),
48 64 : _sampler(getSampler("sampler")),
49 128 : _ci_levels(getParam<std::vector<Real>>("ci_levels")),
50 128 : _ci_replicates(getParam<unsigned int>("ci_replicates")),
51 128 : _ci_seed(getParam<unsigned int>("ci_seed")),
52 64 : _initialized(false)
53 : {
54 64 : if (!dynamic_cast<MorrisSampler *>(&_sampler))
55 0 : paramError("sampler", "Computing Morris sensitivities requires the use of a Morris sampler.");
56 :
57 : // CI levels error checking
58 64 : if (!_ci_levels.empty())
59 : {
60 64 : if (*std::min_element(_ci_levels.begin(), _ci_levels.end()) <= 0)
61 0 : paramError("ci_levels", "The supplied levels must be greater than zero.");
62 64 : else if (*std::max_element(_ci_levels.begin(), _ci_levels.end()) >= 1)
63 0 : paramError("ci_levels", "The supplied levels must be less than 1.0");
64 : }
65 :
66 256 : if ((!isParamValid("reporters") && !isParamValid("vectorpostprocessors")) ||
67 160 : (getParam<std::vector<ReporterName>>("reporters").empty() &&
68 128 : getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors").empty()))
69 0 : mooseError(
70 : "The 'vectorpostprocessors' and/or 'reporters' parameters must be defined and non-empty.");
71 64 : }
72 :
73 : void
74 64 : MorrisReporter::initialize()
75 : {
76 64 : if (_initialized)
77 : return;
78 :
79 : // Stats for Reporters
80 128 : if (isParamValid("reporters"))
81 : {
82 : std::vector<std::string> unsupported_types;
83 224 : for (const auto & r_name : getParam<std::vector<ReporterName>>("reporters"))
84 : {
85 96 : if (hasReporterValueByName<std::vector<Real>>(r_name))
86 64 : declareValueHelper<Real>(r_name);
87 32 : else if (hasReporterValueByName<std::vector<std::vector<Real>>>(r_name))
88 32 : declareValueHelper<std::vector<Real>>(r_name);
89 : else
90 0 : unsupported_types.push_back(r_name.getCombinedName());
91 : }
92 :
93 64 : if (!unsupported_types.empty())
94 0 : paramError("reporters",
95 : "The following reporter value(s) do not have a type supported by the "
96 : "MorrisReporter:\n",
97 0 : MooseUtils::join(unsupported_types, ", "));
98 64 : }
99 :
100 : // Stats for VPP
101 128 : if (isParamValid("vectorpostprocessors"))
102 32 : for (const auto & vpp_name :
103 96 : getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors"))
104 32 : for (const auto & vec_name :
105 96 : _fe_problem.getVectorPostprocessorObjectByName(vpp_name).getVectorNames())
106 64 : declareValueHelper<Real>(ReporterName(vpp_name, vec_name));
107 64 : _initialized = true;
108 : }
109 :
110 : void
111 40 : MorrisReporter::store(nlohmann::json & json) const
112 : {
113 40 : Reporter::store(json);
114 40 : if (!_ci_levels.empty())
115 40 : json["confidence_intervals"] = {{"method", "percentile"},
116 : {"levels", _ci_levels},
117 : {"replicates", _ci_replicates},
118 520 : {"seed", _ci_seed}};
119 80 : json["num_params"] = _sampler.getNumberOfCols();
120 40 : json["num_trajectories"] = _sampler.getNumberOfRows() / (_sampler.getNumberOfCols() + 1);
121 80 : if (_sampler.parameters().isParamValid("levels"))
122 120 : json["levels"] = _sampler.parameters().get<unsigned int>("levels");
123 560 : }
124 :
125 : template <typename DataType>
126 : void
127 128 : MorrisReporter::declareValueHelper(const ReporterName & r_name)
128 : {
129 128 : const auto & data = getReporterValueByName<std::vector<DataType>>(r_name);
130 128 : const std::string s_name = r_name.getObjectName() + "_" + r_name.getValueName();
131 128 : if (!_ci_levels.empty())
132 384 : declareValueByName<MorrisState<DataType>, MorrisReporterContext<DataType>>(
133 : s_name,
134 : REPORTER_MODE_ROOT,
135 : _sampler,
136 : data,
137 384 : MooseEnum("percentile", "percentile"),
138 : _ci_levels,
139 : _ci_replicates,
140 : _ci_seed);
141 : else
142 0 : declareValueByName<MorrisState<DataType>, MorrisReporterContext<DataType>>(
143 : s_name, REPORTER_MODE_ROOT, _sampler, data);
144 128 : }
145 :
146 : template <typename DataType>
147 128 : MorrisReporterContext<DataType>::MorrisReporterContext(const libMesh::ParallelObject & other,
148 : const MooseObject & producer,
149 : ReporterState<MorrisState<DataType>> & state,
150 : Sampler & sampler,
151 : const std::vector<DataType> & data)
152 : : ReporterGeneralContext<MorrisState<DataType>>(other, producer, state),
153 128 : _sampler(sampler),
154 128 : _data(data)
155 : {
156 128 : MultiMooseEnum items("mean meanabs stddev", "mean meanabs stddev", true);
157 256 : _mu_calc = StochasticTools::makeCalculator<std::vector<DataType>, DataType>(items[0], other);
158 256 : _mustar_calc = StochasticTools::makeCalculator<std::vector<DataType>, DataType>(items[1], other);
159 256 : _sig_calc = StochasticTools::makeCalculator<std::vector<DataType>, DataType>(items[2], other);
160 :
161 : // Initialize state
162 128 : auto & mu = this->_state.value()["mu"].first;
163 128 : mu.resize(_sampler.getNumberOfCols());
164 128 : auto & mu_star = this->_state.value()["mu_star"].first;
165 128 : mu_star.resize(_sampler.getNumberOfCols());
166 128 : auto & sig = this->_state.value()["sigma"].first;
167 128 : sig.resize(_sampler.getNumberOfCols());
168 128 : }
169 :
170 : template <typename DataType>
171 128 : MorrisReporterContext<DataType>::MorrisReporterContext(const libMesh::ParallelObject & other,
172 : const MooseObject & producer,
173 : ReporterState<MorrisState<DataType>> & state,
174 : Sampler & sampler,
175 : const std::vector<DataType> & data,
176 : const MooseEnum & ci_method,
177 : const std::vector<Real> & ci_levels,
178 : unsigned int ci_replicates,
179 : unsigned int ci_seed)
180 128 : : MorrisReporterContext<DataType>(other, producer, state, sampler, data)
181 : {
182 128 : _ci_mu_calc = StochasticTools::makeBootstrapCalculator<std::vector<DataType>, DataType>(
183 : ci_method, other, ci_levels, ci_replicates, ci_seed, *_mu_calc);
184 128 : _ci_mustar_calc = StochasticTools::makeBootstrapCalculator<std::vector<DataType>, DataType>(
185 : ci_method, other, ci_levels, ci_replicates, ci_seed, *_mustar_calc);
186 128 : _ci_sig_calc = StochasticTools::makeBootstrapCalculator<std::vector<DataType>, DataType>(
187 : ci_method, other, ci_levels, ci_replicates, ci_seed, *_sig_calc);
188 :
189 : // Initialize state
190 128 : auto & mu = this->_state.value()["mu"].second;
191 128 : mu.resize(_sampler.getNumberOfCols());
192 128 : auto & mu_star = this->_state.value()["mu_star"].second;
193 128 : mu_star.resize(_sampler.getNumberOfCols());
194 128 : auto & sig = this->_state.value()["sigma"].second;
195 128 : sig.resize(_sampler.getNumberOfCols());
196 128 : }
197 :
198 : template <typename DataType>
199 : void
200 128 : MorrisReporterContext<DataType>::finalize()
201 : {
202 : dof_id_type offset;
203 128 : if (_data.size() == _sampler.getNumberOfRows())
204 56 : offset = _sampler.getLocalRowBegin();
205 72 : else if (_data.size() == _sampler.getNumberOfLocalRows())
206 : offset = 0;
207 : else
208 0 : mooseError("Data size inconsistency. Expected data vector to have ",
209 0 : _sampler.getNumberOfLocalRows(),
210 : " (local) or ",
211 0 : _sampler.getNumberOfRows(),
212 : " (global) elements, but actually has ",
213 0 : _data.size(),
214 : " elements. Are you using the same sampler?");
215 :
216 128 : const dof_id_type nc = _sampler.getNumberOfCols(); // convenience
217 128 : if (_sampler.getNumberOfLocalRows() % (nc + 1) > 0)
218 0 : mooseError(
219 : "Sampler does not have an appropriate number of rows. Are you using a Morris sampler?");
220 :
221 128 : std::vector<std::vector<DataType>> elem_effects(
222 160 : nc, std::vector<DataType>(_sampler.getNumberOfLocalRows() / (nc + 1)));
223 128 : RealEigenMatrix x(nc + 1, nc);
224 128 : std::vector<DataType> y(nc + 1);
225 :
226 147688 : for (dof_id_type p = 0; p < _sampler.getNumberOfLocalRows(); ++p)
227 : {
228 147560 : dof_id_type traj = p / (nc + 1);
229 147560 : dof_id_type tind = p % (nc + 1);
230 147560 : const std::vector<Real> row = _sampler.getNextLocalRow();
231 1032920 : for (unsigned int k = 0; k < nc; ++k)
232 885360 : x(tind, k) = row[k];
233 147560 : y[tind] = _data[p + offset];
234 :
235 147560 : if (tind == nc)
236 : {
237 21080 : const std::vector<DataType> ee = computeElementaryEffects(x, y);
238 147560 : for (unsigned int k = 0; k < nc; ++k)
239 126480 : elem_effects[k][traj] = ee[k];
240 200 : }
241 : }
242 :
243 128 : auto & mu = this->_state.value()["mu"];
244 128 : auto & mustar = this->_state.value()["mu_star"];
245 128 : auto & sig = this->_state.value()["sigma"];
246 896 : for (unsigned int k = 0; k < nc; ++k)
247 : {
248 960 : mu.first[k] = _mu_calc->compute(elem_effects[k], true);
249 960 : mustar.first[k] = _mustar_calc->compute(elem_effects[k], true);
250 960 : sig.first[k] = _sig_calc->compute(elem_effects[k], true);
251 :
252 768 : if (_ci_mu_calc)
253 1536 : mu.second[k] = _ci_mu_calc->compute(elem_effects[k], true);
254 768 : if (_ci_mustar_calc)
255 1536 : mustar.second[k] = _ci_mustar_calc->compute(elem_effects[k], true);
256 768 : if (_ci_sig_calc)
257 1536 : sig.second[k] = _ci_sig_calc->compute(elem_effects[k], true);
258 : }
259 128 : }
260 :
261 : template <>
262 : std::vector<Real>
263 20880 : MorrisReporterContext<Real>::computeElementaryEffects(const RealEigenMatrix & x,
264 : const std::vector<Real> & y) const
265 : {
266 20880 : const auto k = y.size() - 1;
267 20880 : const RealEigenMatrix dx = x.bottomRows(k) - x.topRows(k);
268 20880 : RealEigenVector dy(k);
269 146160 : for (unsigned int j = 0; j < k; ++j)
270 125280 : dy(j) = y[j + 1] - y[j];
271 :
272 41760 : const RealEigenVector u = dx.fullPivLu().solve(dy);
273 20880 : return std::vector<Real>(u.data(), u.data() + u.size());
274 : }
275 :
276 : template <>
277 : std::vector<std::vector<Real>>
278 200 : MorrisReporterContext<std::vector<Real>>::computeElementaryEffects(
279 : const RealEigenMatrix & x, const std::vector<std::vector<Real>> & y) const
280 : {
281 200 : const auto k = y.size() - 1;
282 200 : const RealEigenMatrix dx = x.bottomRows(k) - x.topRows(k);
283 : const auto solver = dx.fullPivLu();
284 200 : RealEigenVector dy(k);
285 200 : std::vector<std::vector<Real>> ee(k, std::vector<Real>(y[0].size()));
286 600 : for (unsigned int i = 0; i < y[0].size(); ++i)
287 : {
288 2800 : for (unsigned int j = 0; j < k; ++j)
289 2400 : dy(j) = y[j + 1][i] - y[j][i];
290 400 : const RealEigenVector u = solver.solve(dy);
291 2800 : for (unsigned int j = 0; j < k; ++j)
292 2400 : ee[j][i] = u(j);
293 : }
294 200 : return ee;
295 200 : }
296 :
297 : template void MorrisReporter::declareValueHelper<Real>(const ReporterName & r_name);
298 : template class MorrisReporterContext<Real>;
299 : template void MorrisReporter::declareValueHelper<std::vector<Real>>(const ReporterName & r_name);
300 : template class MorrisReporterContext<std::vector<Real>>;
|
