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 15 : ParallelSubsetSimulation::validParams()
19 : {
20 15 : InputParameters params = Sampler::validParams();
21 15 : params.addClassDescription("Parallel Subset Simulation sampler.");
22 30 : 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 30 : params.addRequiredParam<ReporterName>("output_reporter",
27 : "Reporter with results of samples created by the SubApp.");
28 30 : params.addRequiredParam<ReporterName>("inputs_reporter", "Reporter with input parameters.");
29 45 : params.addRangeCheckedParam<Real>("subset_probability",
30 30 : 0.1,
31 : "subset_probability>0 & subset_probability<=1",
32 : "Conditional probability of each subset");
33 30 : params.addRequiredParam<unsigned int>("num_samplessub", "Number of samples per subset");
34 30 : params.addRequiredParam<unsigned int>("num_subsets", "Number of desired subsets");
35 30 : 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 30 : params.addParam<bool>("use_absolute_value", false, "Use absolute value of the sub app output");
39 30 : params.addParam<unsigned int>(
40 : "num_random_seeds",
41 30 : 100000,
42 : "Initialize a certain number of random seeds. Change from the default only if you have to.");
43 15 : return params;
44 0 : }
45 :
46 8 : ParallelSubsetSimulation::ParallelSubsetSimulation(const InputParameters & parameters)
47 : : Sampler(parameters),
48 8 : _num_samplessub(getParam<unsigned int>("num_samplessub")),
49 16 : _num_subsets(getParam<unsigned int>("num_subsets")),
50 16 : _use_absolute_value(getParam<bool>("use_absolute_value")),
51 16 : _subset_probability(getParam<Real>("subset_probability")),
52 16 : _num_random_seeds(getParam<unsigned int>("num_random_seeds")),
53 8 : _outputs(getReporterValue<std::vector<Real>>("output_reporter")),
54 8 : _inputs(getReporterValue<std::vector<std::vector<Real>>>("inputs_reporter")),
55 16 : _step(getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")->timeStep()),
56 8 : _count_max(std::floor(1 / _subset_probability)),
57 8 : _subset(0),
58 8 : _is_sampling_completed(false)
59 : {
60 : // Fixing the number of rows to the number of processors
61 8 : const dof_id_type nchains = isParamValid("num_parallel_chains")
62 24 : ? getParam<unsigned int>("num_parallel_chains")
63 0 : : n_processors() / _min_procs_per_row;
64 8 : setNumberOfRows(nchains);
65 8 : 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 32 : for (const DistributionName & name : getParam<std::vector<DistributionName>>("distributions"))
74 16 : _distributions.push_back(&getDistributionByName(name));
75 :
76 : // Setting the number of columns in the sampler matrix (equal to the number of distributions).
77 8 : 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 8 : _inputs_sto.resize(_distributions.size(), std::vector<Real>(_num_samplessub, 0.0));
82 8 : _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 8 : _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 8 : _markov_seed.resize(_distributions.size());
92 :
93 8 : setNumberOfRandomSeeds(_num_random_seeds);
94 8 : setAutoAdvanceGenerators(false);
95 8 : }
96 :
97 : const unsigned int &
98 126 : ParallelSubsetSimulation::getNumSamplesSub() const
99 : {
100 126 : return _num_samplessub;
101 : }
102 :
103 : const bool &
104 42 : ParallelSubsetSimulation::getUseAbsoluteValue() const
105 : {
106 42 : return _use_absolute_value;
107 : }
108 :
109 : const Real &
110 70 : ParallelSubsetSimulation::getSubsetProbability() const
111 : {
112 70 : return _subset_probability;
113 : }
114 :
115 : void
116 56 : ParallelSubsetSimulation::executeSetUp()
117 : {
118 56 : 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 56 : _subset = (_step * getNumberOfRows()) / _num_samplessub;
123 56 : const unsigned int sub_ind = _step - (_num_samplessub / getNumberOfRows()) * _subset;
124 56 : const unsigned int offset = sub_ind * getNumberOfRows();
125 :
126 : // check if we have completed the last sample
127 56 : if (_subset >= _num_subsets)
128 : {
129 8 : _is_sampling_completed = true;
130 8 : return;
131 : }
132 :
133 : // Get and store the accepted samples input across all the procs from the previous step
134 144 : for (dof_id_type j = 0; j < _distributions.size(); ++j)
135 288 : for (dof_id_type ss = 0; ss < getNumberOfRows(); ++ss)
136 192 : _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 48 : _use_absolute_value ? AdaptiveMonteCarloUtils::computeVectorABS(_outputs) : _outputs;
141 48 : _communicator.allgather(tmp);
142 144 : for (dof_id_type ss = 0; ss < getNumberOfRows(); ++ss)
143 96 : _outputs_sto[ss + offset] = tmp[ss];
144 :
145 : // These are the subsequent subsets which use Markov Chain Monte Carlo sampling scheme
146 48 : if (_subset > 0)
147 : {
148 : // Check whether the subset index has changed
149 32 : if (sub_ind == 0)
150 : {
151 : // _inputs_sorted contains the input values corresponding to the largest po percentile
152 : // output values
153 16 : _inputs_sorted = AdaptiveMonteCarloUtils::sortInput(
154 16 : _inputs_sto, _outputs_sto, _num_samplessub, _subset_probability);
155 : }
156 :
157 : // Reinitialize the starting inputs values for the next set of Markov chains
158 32 : if (sub_ind % _count_max == 0)
159 : {
160 16 : const unsigned int soffset = (sub_ind / _count_max) * getNumberOfRows();
161 48 : for (dof_id_type j = 0; j < _distributions.size(); ++j)
162 32 : _markov_seed[j].assign(_inputs_sorted[j].begin() + soffset,
163 32 : _inputs_sorted[j].begin() + soffset + getNumberOfRows());
164 : }
165 : // Otherwise, use the previously accepted input values to propose the next set of input
166 : // values
167 : else
168 : {
169 48 : for (dof_id_type j = 0; j < _distributions.size(); ++j)
170 32 : _markov_seed[j].assign(_inputs_sto[j].begin() + offset,
171 32 : _inputs_sto[j].begin() + offset + getNumberOfRows());
172 : }
173 : }
174 48 : }
175 :
176 : Real
177 240 : ParallelSubsetSimulation::computeSample(dof_id_type row_index, dof_id_type col_index)
178 : {
179 240 : unsigned int seed_value = _step > 0 ? (_step - 1) * 2 : 0;
180 240 : const dof_id_type n = row_index * getNumberOfCols() + col_index;
181 : Real val;
182 :
183 240 : if (_subset == 0)
184 80 : val = getRand(n, seed_value);
185 : else
186 : {
187 : const Real rv =
188 160 : Normal::quantile(getRand(n, seed_value), _markov_seed[col_index][row_index], 1.0);
189 160 : const Real acceptance_ratio = std::log(Normal::pdf(rv, 0, 1)) -
190 160 : std::log(Normal::pdf(_markov_seed[col_index][row_index], 0, 1));
191 160 : const Real new_sample = acceptance_ratio > std::log(getRand(n, seed_value + 1))
192 160 : ? rv
193 160 : : _markov_seed[col_index][row_index];
194 160 : val = Normal::cdf(new_sample, 0, 1);
195 : }
196 :
197 240 : return _distributions[col_index]->quantile(val);
198 : }
|
