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 136 : MorrisReporter::validParams()
18 : {
19 136 : InputParameters params = GeneralReporter::validParams();
20 136 : params.addClassDescription("Compute global sensitivities using the Morris method.");
21 :
22 272 : params.addRequiredParam<SamplerName>("sampler", "Morris sampler used to generate samples.");
23 272 : params.addParam<std::vector<VectorPostprocessorName>>(
24 : "vectorpostprocessors",
25 : "List of VectorPostprocessor(s) to utilized for statistic computations.");
26 272 : params.addParam<std::vector<ReporterName>>(
27 : "reporters", {}, "List of Reporter values to utilized for statistic computations.");
28 :
29 : // Confidence Levels
30 136 : params.addParam<std::vector<Real>>(
31 : "ci_levels",
32 136 : std::vector<Real>(),
33 : "A vector of confidence levels to consider, values must be in (0, 1).");
34 272 : params.addParam<unsigned int>("ci_replicates",
35 272 : 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 272 : params.addParam<unsigned int>("ci_seed",
40 272 : 1,
41 : "The random number generator seed used for creating replicates "
42 : "while computing confidence level intervals.");
43 136 : return params;
44 0 : }
45 :
46 68 : MorrisReporter::MorrisReporter(const InputParameters & parameters)
47 : : GeneralReporter(parameters),
48 68 : _sampler(getSampler("sampler")),
49 136 : _ci_levels(getParam<std::vector<Real>>("ci_levels")),
50 136 : _ci_replicates(getParam<unsigned int>("ci_replicates")),
51 136 : _ci_seed(getParam<unsigned int>("ci_seed")),
52 68 : _initialized(false)
53 : {
54 68 : 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 68 : if (!_ci_levels.empty())
59 : {
60 68 : 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 68 : 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 272 : if ((!isParamValid("reporters") && !isParamValid("vectorpostprocessors")) ||
67 170 : (getParam<std::vector<ReporterName>>("reporters").empty() &&
68 136 : getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors").empty()))
69 0 : mooseError(
70 : "The 'vectorpostprocessors' and/or 'reporters' parameters must be defined and non-empty.");
71 68 : }
72 :
73 : void
74 68 : MorrisReporter::initialize()
75 : {
76 68 : if (_initialized)
77 : return;
78 :
79 : // Stats for Reporters
80 136 : if (isParamValid("reporters"))
81 : {
82 : std::vector<std::string> unsupported_types;
83 238 : for (const auto & r_name : getParam<std::vector<ReporterName>>("reporters"))
84 : {
85 102 : if (hasReporterValueByName<std::vector<Real>>(r_name))
86 68 : declareValueHelper<Real>(r_name);
87 34 : else if (hasReporterValueByName<std::vector<std::vector<Real>>>(r_name))
88 34 : declareValueHelper<std::vector<Real>>(r_name);
89 : else
90 0 : unsupported_types.push_back(r_name.getCombinedName());
91 : }
92 :
93 68 : 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 68 : }
99 :
100 : // Stats for VPP
101 136 : if (isParamValid("vectorpostprocessors"))
102 34 : for (const auto & vpp_name :
103 102 : getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors"))
104 34 : for (const auto & vec_name :
105 102 : _fe_problem.getVectorPostprocessorObjectByName(vpp_name).getVectorNames())
106 68 : declareValueHelper<Real>(ReporterName(vpp_name, vec_name));
107 68 : _initialized = true;
108 : }
109 :
110 : void
111 44 : MorrisReporter::store(nlohmann::json & json) const
112 : {
113 44 : Reporter::store(json);
114 44 : if (!_ci_levels.empty())
115 44 : json["confidence_intervals"] = {{"method", "percentile"},
116 : {"levels", _ci_levels},
117 : {"replicates", _ci_replicates},
118 572 : {"seed", _ci_seed}};
119 88 : json["num_params"] = _sampler.getNumberOfCols();
120 44 : json["num_trajectories"] = _sampler.getNumberOfRows() / (_sampler.getNumberOfCols() + 1);
121 88 : if (_sampler.parameters().isParamValid("levels"))
122 132 : json["levels"] = _sampler.parameters().get<unsigned int>("levels");
123 616 : }
124 :
125 : template <typename DataType>
126 : void
127 136 : MorrisReporter::declareValueHelper(const ReporterName & r_name)
128 : {
129 136 : const auto & data = getReporterValueByName<std::vector<DataType>>(r_name);
130 136 : const std::string s_name = r_name.getObjectName() + "_" + r_name.getValueName();
131 136 : if (!_ci_levels.empty())
132 408 : declareValueByName<MorrisState<DataType>, MorrisReporterContext<DataType>>(
133 : s_name,
134 : REPORTER_MODE_ROOT,
135 : _sampler,
136 : data,
137 408 : 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 136 : }
145 :
146 : template <typename DataType>
147 136 : 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 136 : _sampler(sampler),
154 136 : _data(data)
155 : {
156 136 : MultiMooseEnum items("mean meanabs stddev", "mean meanabs stddev", true);
157 272 : _mu_calc = StochasticTools::makeCalculator<std::vector<DataType>, DataType>(items[0], other);
158 272 : _mustar_calc = StochasticTools::makeCalculator<std::vector<DataType>, DataType>(items[1], other);
159 272 : _sig_calc = StochasticTools::makeCalculator<std::vector<DataType>, DataType>(items[2], other);
160 :
161 : // Initialize state
162 136 : auto & mu = this->_state.value()["mu"].first;
163 136 : mu.resize(_sampler.getNumberOfCols());
164 136 : auto & mu_star = this->_state.value()["mu_star"].first;
165 136 : mu_star.resize(_sampler.getNumberOfCols());
166 136 : auto & sig = this->_state.value()["sigma"].first;
167 136 : sig.resize(_sampler.getNumberOfCols());
168 136 : }
169 :
170 : template <typename DataType>
171 136 : 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 136 : : MorrisReporterContext<DataType>(other, producer, state, sampler, data)
181 : {
182 136 : _ci_mu_calc = StochasticTools::makeBootstrapCalculator<std::vector<DataType>, DataType>(
183 : ci_method, other, ci_levels, ci_replicates, ci_seed, *_mu_calc);
184 136 : _ci_mustar_calc = StochasticTools::makeBootstrapCalculator<std::vector<DataType>, DataType>(
185 : ci_method, other, ci_levels, ci_replicates, ci_seed, *_mustar_calc);
186 136 : _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 136 : auto & mu = this->_state.value()["mu"].second;
191 136 : mu.resize(_sampler.getNumberOfCols());
192 136 : auto & mu_star = this->_state.value()["mu_star"].second;
193 136 : mu_star.resize(_sampler.getNumberOfCols());
194 136 : auto & sig = this->_state.value()["sigma"].second;
195 136 : sig.resize(_sampler.getNumberOfCols());
196 136 : }
197 :
198 : template <typename DataType>
199 : void
200 136 : MorrisReporterContext<DataType>::finalize()
201 : {
202 : dof_id_type offset;
203 136 : if (_data.size() == _sampler.getNumberOfRows())
204 64 : 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 136 : const dof_id_type nc = _sampler.getNumberOfCols(); // convenience
217 136 : 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 136 : std::vector<std::vector<DataType>> elem_effects(
222 272 : nc, std::vector<DataType>(_sampler.getNumberOfLocalRows() / (nc + 1)));
223 136 : RealEigenMatrix x(nc + 1, nc);
224 136 : std::vector<DataType> y(nc + 1);
225 :
226 162452 : for (dof_id_type p = 0; p < _sampler.getNumberOfLocalRows(); ++p)
227 : {
228 162316 : dof_id_type traj = p / (nc + 1);
229 162316 : dof_id_type tind = p % (nc + 1);
230 162316 : const std::vector<Real> row = _sampler.getNextLocalRow();
231 1136212 : for (unsigned int k = 0; k < nc; ++k)
232 973896 : x(tind, k) = row[k];
233 162316 : y[tind] = _data[p + offset];
234 :
235 162316 : if (tind == nc)
236 : {
237 23188 : const std::vector<DataType> ee = computeElementaryEffects(x, y);
238 162316 : for (unsigned int k = 0; k < nc; ++k)
239 139128 : elem_effects[k][traj] = ee[k];
240 23188 : }
241 : }
242 :
243 136 : auto & mu = this->_state.value()["mu"];
244 136 : auto & mustar = this->_state.value()["mu_star"];
245 136 : auto & sig = this->_state.value()["sigma"];
246 952 : for (unsigned int k = 0; k < nc; ++k)
247 : {
248 1020 : mu.first[k] = _mu_calc->compute(elem_effects[k], true);
249 1020 : mustar.first[k] = _mustar_calc->compute(elem_effects[k], true);
250 1020 : sig.first[k] = _sig_calc->compute(elem_effects[k], true);
251 :
252 816 : if (_ci_mu_calc)
253 1632 : mu.second[k] = _ci_mu_calc->compute(elem_effects[k], true);
254 816 : if (_ci_mustar_calc)
255 1632 : mustar.second[k] = _ci_mustar_calc->compute(elem_effects[k], true);
256 816 : if (_ci_sig_calc)
257 1632 : sig.second[k] = _ci_sig_calc->compute(elem_effects[k], true);
258 : }
259 136 : }
260 :
261 : template <>
262 : std::vector<Real>
263 22968 : MorrisReporterContext<Real>::computeElementaryEffects(const RealEigenMatrix & x,
264 : const std::vector<Real> & y) const
265 : {
266 22968 : const auto k = y.size() - 1;
267 22968 : const RealEigenMatrix dx = x.bottomRows(k) - x.topRows(k);
268 22968 : RealEigenVector dy(k);
269 160776 : for (unsigned int j = 0; j < k; ++j)
270 137808 : dy(j) = y[j + 1] - y[j];
271 :
272 45936 : const RealEigenVector u = dx.fullPivLu().solve(dy);
273 22968 : return std::vector<Real>(u.data(), u.data() + u.size());
274 : }
275 :
276 : template <>
277 : std::vector<std::vector<Real>>
278 220 : MorrisReporterContext<std::vector<Real>>::computeElementaryEffects(
279 : const RealEigenMatrix & x, const std::vector<std::vector<Real>> & y) const
280 : {
281 220 : const auto k = y.size() - 1;
282 220 : const RealEigenMatrix dx = x.bottomRows(k) - x.topRows(k);
283 : const auto solver = dx.fullPivLu();
284 220 : RealEigenVector dy(k);
285 220 : std::vector<std::vector<Real>> ee(k, std::vector<Real>(y[0].size()));
286 660 : for (unsigned int i = 0; i < y[0].size(); ++i)
287 : {
288 3080 : for (unsigned int j = 0; j < k; ++j)
289 2640 : dy(j) = y[j + 1][i] - y[j][i];
290 440 : const RealEigenVector u = solver.solve(dy);
291 3080 : for (unsigned int j = 0; j < k; ++j)
292 2640 : ee[j][i] = u(j);
293 : }
294 220 : return ee;
295 220 : }
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>>;
|
