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 "ParallelSubsetSimulation.h"
11 : #include "AdaptiveMonteCarloUtils.h"
12 : #include "Normal.h"
13 : #include "Uniform.h"
14 :
15 : registerMooseObject("StochasticToolsApp", ParallelSubsetSimulation);
16 :
17 : InputParameters
18 30 : ParallelSubsetSimulation::validParams()
19 : {
20 30 : InputParameters params = Sampler::validParams();
21 30 : params.addClassDescription("Parallel Subset Simulation sampler.");
22 60 : params.addRequiredParam<std::vector<DistributionName>>(
23 : "distributions",
24 : "The distribution names to be sampled, the number of distributions provided defines the "
25 : "number of columns per matrix.");
26 60 : params.addRequiredParam<ReporterName>("output_reporter",
27 : "Reporter with results of samples created by the SubApp.");
28 60 : params.addRequiredParam<ReporterName>("inputs_reporter", "Reporter with input parameters.");
29 90 : params.addRangeCheckedParam<Real>("subset_probability",
30 60 : 0.1,
31 : "subset_probability>0 & subset_probability<=1",
32 : "Conditional probability of each subset");
33 60 : params.addRequiredParam<unsigned int>("num_samplessub", "Number of samples per subset");
34 60 : params.addRequiredParam<unsigned int>("num_subsets", "Number of desired subsets");
35 60 : params.addParam<unsigned int>("num_parallel_chains",
36 : "Number of Markov chains to run in parallel, default is based on "
37 : "the number of processors used.");
38 60 : params.addParam<bool>("use_absolute_value", false, "Use absolute value of the sub app output");
39 60 : params.addParam<unsigned int>(
40 : "num_random_seeds",
41 60 : 100000,
42 : "Initialize a certain number of random seeds. Change from the default only if you have to.");
43 30 : return params;
44 0 : }
45 :
46 16 : ParallelSubsetSimulation::ParallelSubsetSimulation(const InputParameters & parameters)
47 : : Sampler(parameters),
48 16 : _num_samplessub(getParam<unsigned int>("num_samplessub")),
49 32 : _num_subsets(getParam<unsigned int>("num_subsets")),
50 32 : _use_absolute_value(getParam<bool>("use_absolute_value")),
51 32 : _subset_probability(getParam<Real>("subset_probability")),
52 32 : _num_random_seeds(getParam<unsigned int>("num_random_seeds")),
53 16 : _outputs(getReporterValue<std::vector<Real>>("output_reporter")),
54 16 : _inputs(getReporterValue<std::vector<std::vector<Real>>>("inputs_reporter")),
55 32 : _step(getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")->timeStep()),
56 16 : _count_max(std::floor(1 / _subset_probability)),
57 16 : _subset(0),
58 16 : _is_sampling_completed(false)
59 : {
60 : // Fixing the number of rows to the number of processors
61 16 : const dof_id_type nchains = isParamValid("num_parallel_chains")
62 48 : ? getParam<unsigned int>("num_parallel_chains")
63 0 : : n_processors() / _min_procs_per_row;
64 16 : setNumberOfRows(nchains);
65 16 : if ((_num_samplessub / nchains) % _count_max > 0)
66 0 : mooseError("Number of model evaluations per chain per subset (",
67 0 : _num_samplessub / nchains,
68 : ") should be a multiple of requested chain length (",
69 0 : _count_max,
70 : ").");
71 :
72 : // Filling the `distributions` vector with the user-provided distributions.
73 64 : for (const DistributionName & name : getParam<std::vector<DistributionName>>("distributions"))
74 32 : _distributions.push_back(&getDistributionByName(name));
75 :
76 : // Setting the number of columns in the sampler matrix (equal to the number of distributions).
77 16 : setNumberOfCols(_distributions.size());
78 :
79 : /* `inputs_sto` is a member variable that aids in deciding the next set of samples
80 : in the Subset Simulation algorithm by storing the input parameter values*/
81 16 : _inputs_sto.resize(_distributions.size(), std::vector<Real>(_num_samplessub, 0.0));
82 16 : _outputs_sto.resize(_num_samplessub, 0.0);
83 :
84 : /* `inputs_sorted` is a member variable which also aids in deciding the next set of samples
85 : in the Subset Simulation algorithm by storing the sorted input parameter values
86 : by their corresponding output values*/
87 16 : _inputs_sorted.resize(_distributions.size());
88 :
89 : /* `markov_seed` is a member variable to store the seed input values for proposing
90 : the next set of Markov chain samples.*/
91 16 : _markov_seed.resize(_distributions.size());
92 :
93 16 : setNumberOfRandomSeeds(_num_random_seeds);
94 16 : setAutoAdvanceGenerators(false);
95 16 : }
96 :
97 : const unsigned int &
98 252 : ParallelSubsetSimulation::getNumSamplesSub() const
99 : {
100 252 : return _num_samplessub;
101 : }
102 :
103 : const bool &
104 84 : ParallelSubsetSimulation::getUseAbsoluteValue() const
105 : {
106 84 : return _use_absolute_value;
107 : }
108 :
109 : const Real &
110 140 : ParallelSubsetSimulation::getSubsetProbability() const
111 : {
112 140 : return _subset_probability;
113 : }
114 :
115 : void
116 112 : ParallelSubsetSimulation::executeSetUp()
117 : {
118 112 : if (_is_sampling_completed)
119 0 : mooseError("Internal bug: the adaptive sampling is supposed to be completed but another sample "
120 : "has been requested.");
121 :
122 112 : _subset = (_step * getNumberOfRows()) / _num_samplessub;
123 112 : const unsigned int sub_ind = _step - (_num_samplessub / getNumberOfRows()) * _subset;
124 112 : const unsigned int offset = sub_ind * getNumberOfRows();
125 :
126 : // check if we have completed the last sample
127 112 : if (_subset >= _num_subsets)
128 : {
129 16 : _is_sampling_completed = true;
130 16 : return;
131 : }
132 :
133 : // Get and store the accepted samples input across all the procs from the previous step
134 288 : for (dof_id_type j = 0; j < _distributions.size(); ++j)
135 576 : for (dof_id_type ss = 0; ss < getNumberOfRows(); ++ss)
136 384 : _inputs_sto[j][ss + offset] = Normal::quantile(_distributions[j]->cdf(_inputs[j][ss]), 0, 1);
137 :
138 : // Get the accepted sample outputs across all the procs from the previous step
139 : std::vector<Real> tmp =
140 96 : _use_absolute_value ? AdaptiveMonteCarloUtils::computeVectorABS(_outputs) : _outputs;
141 96 : _communicator.allgather(tmp);
142 96 : if (tmp.empty())
143 8 : tmp.resize(getNumberOfRows(), 0.0);
144 288 : for (dof_id_type ss = 0; ss < getNumberOfRows(); ++ss)
145 192 : _outputs_sto[ss + offset] = tmp[ss];
146 :
147 : // These are the subsequent subsets which use Markov Chain Monte Carlo sampling scheme
148 96 : if (_subset > 0)
149 : {
150 : // Check whether the subset index has changed
151 64 : if (sub_ind == 0)
152 : {
153 : // _inputs_sorted contains the input values corresponding to the largest po percentile
154 : // output values
155 32 : _inputs_sorted = AdaptiveMonteCarloUtils::sortInput(
156 32 : _inputs_sto, _outputs_sto, _num_samplessub, _subset_probability);
157 : }
158 :
159 : // Reinitialize the starting inputs values for the next set of Markov chains
160 64 : if (sub_ind % _count_max == 0)
161 : {
162 32 : const unsigned int soffset = (sub_ind / _count_max) * getNumberOfRows();
163 96 : for (dof_id_type j = 0; j < _distributions.size(); ++j)
164 64 : _markov_seed[j].assign(_inputs_sorted[j].begin() + soffset,
165 64 : _inputs_sorted[j].begin() + soffset + getNumberOfRows());
166 : }
167 : // Otherwise, use the previously accepted input values to propose the next set of input
168 : // values
169 : else
170 : {
171 96 : for (dof_id_type j = 0; j < _distributions.size(); ++j)
172 64 : _markov_seed[j].assign(_inputs_sto[j].begin() + offset,
173 64 : _inputs_sto[j].begin() + offset + getNumberOfRows());
174 : }
175 : }
176 96 : }
177 :
178 : Real
179 480 : ParallelSubsetSimulation::computeSample(dof_id_type row_index, dof_id_type col_index)
180 : {
181 480 : unsigned int seed_value = _step > 0 ? (_step - 1) * 2 : 0;
182 480 : const dof_id_type n = row_index * getNumberOfCols() + col_index;
183 : Real val;
184 :
185 480 : if (_subset == 0)
186 160 : val = getRand(n, seed_value);
187 : else
188 : {
189 : const Real rv =
190 320 : Normal::quantile(getRand(n, seed_value), _markov_seed[col_index][row_index], 1.0);
191 320 : const Real acceptance_ratio = std::log(Normal::pdf(rv, 0, 1)) -
192 320 : std::log(Normal::pdf(_markov_seed[col_index][row_index], 0, 1));
193 320 : const Real new_sample = acceptance_ratio > std::log(getRand(n, seed_value + 1))
194 320 : ? rv
195 320 : : _markov_seed[col_index][row_index];
196 320 : val = Normal::cdf(new_sample, 0, 1);
197 : }
198 :
199 480 : return _distributions[col_index]->quantile(val);
200 : }
|
