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 "PMCMCDecision.h"
11 : #include "Sampler.h"
12 : #include "DenseMatrix.h"
13 :
14 : registerMooseObject("StochasticToolsApp", PMCMCDecision);
15 :
16 : InputParameters
17 250 : PMCMCDecision::validParams()
18 : {
19 250 : InputParameters params = GeneralReporter::validParams();
20 250 : params += LikelihoodInterface::validParams();
21 250 : params.addClassDescription("Generic reporter which decides whether or not to accept a proposed "
22 : "sample in parallel Markov chain Monte Carlo type of algorithms.");
23 500 : params.addParam<ReporterName>("output_value", "Value of the model output from the SubApp.");
24 500 : params.addParam<ReporterValueName>(
25 : "outputs_required",
26 : "outputs_required",
27 : "Modified value of the model output from this reporter class.");
28 500 : params.addParam<ReporterValueName>("inputs", "inputs", "Uncertain inputs to the model.");
29 500 : params.addParam<ReporterValueName>("tpm", "tpm", "The transition probability matrix.");
30 500 : params.addParam<ReporterValueName>("variance", "variance", "Model variance term.");
31 500 : params.addParam<ReporterValueName>(
32 : "noise", "noise", "Model noise term to pass to Likelihoods object.");
33 500 : params.addRequiredParam<SamplerName>("sampler", "The sampler object.");
34 500 : params.addRequiredParam<std::vector<UserObjectName>>("likelihoods", "Names of likelihoods.");
35 250 : return params;
36 0 : }
37 :
38 120 : PMCMCDecision::PMCMCDecision(const InputParameters & parameters)
39 : : GeneralReporter(parameters),
40 : LikelihoodInterface(parameters),
41 120 : _inputs(declareValue<std::vector<std::vector<Real>>>("inputs")),
42 120 : _tpm(declareValue<std::vector<Real>>("tpm")),
43 120 : _variance(declareValue<std::vector<Real>>("variance")),
44 120 : _noise(declareValue<Real>("noise")),
45 120 : _sampler(getSampler("sampler")),
46 120 : _pmcmc(dynamic_cast<const PMCMCBase *>(&_sampler)),
47 120 : _rnd_vec(_pmcmc->getRandomNumbers()),
48 120 : _new_var_samples(_pmcmc->getVarSamples()),
49 120 : _priors(_pmcmc->getPriors()),
50 120 : _var_prior(_pmcmc->getVarPrior()),
51 120 : _outputs_required(
52 120 : isParamValid("output_value")
53 480 : ? &declareValue<std::vector<Real>>("outputs_required", REPORTER_MODE_DISTRIBUTED)
54 : : nullptr),
55 480 : _output_value(isParamValid("output_value") ? &getReporterValue<std::vector<Real>>(
56 : "output_value", REPORTER_MODE_DISTRIBUTED)
57 : : nullptr),
58 120 : _local_comm(_sampler.getLocalComm()),
59 240 : _check_step(std::numeric_limits<int>::max())
60 : {
61 : // Filling the `likelihoods` vector with the user-provided distributions.
62 360 : for (const UserObjectName & name : getParam<std::vector<UserObjectName>>("likelihoods"))
63 120 : _likelihoods.push_back(getLikelihoodFunctionByName(name));
64 :
65 : // Check whether the selected sampler is an MCMC sampler or not
66 120 : if (!_pmcmc)
67 0 : paramError("sampler", "The selected sampler is not of type MCMC.");
68 :
69 : // Fetching the sampler characteristics
70 120 : _props = _pmcmc->getNumParallelProposals();
71 120 : _num_confg_values = _pmcmc->getNumberOfConfigValues();
72 120 : _num_confg_params = _pmcmc->getNumberOfConfigParams();
73 :
74 : // Resizing the data arrays to transmit to the output file
75 120 : _inputs.resize(_props);
76 600 : for (unsigned int i = 0; i < _props; ++i)
77 480 : _inputs[i].resize(_sampler.getNumberOfCols() - _num_confg_params);
78 120 : if (_outputs_required)
79 120 : _outputs_required->resize(_sampler.getNumberOfRows());
80 120 : _tpm.resize(_props);
81 120 : _variance.resize(_props);
82 120 : }
83 :
84 : void
85 600 : PMCMCDecision::initialize()
86 : {
87 1800 : if (!isParamValid("output_value") && !usingGP())
88 0 : paramError("output_value", "Value of the model output from the SubApp should be specified.");
89 600 : }
90 :
91 : void
92 340 : PMCMCDecision::computeEvidence(std::vector<Real> & evidence, const DenseMatrix<Real> & input_matrix)
93 : {
94 340 : std::vector<Real> out1(_num_confg_values);
95 340 : std::vector<Real> out2(_num_confg_values);
96 1560 : for (unsigned int i = 0; i < evidence.size(); ++i)
97 : {
98 1220 : evidence[i] = 0.0;
99 3660 : for (unsigned int j = 0; j < _priors.size(); ++j)
100 2440 : evidence[i] += (std::log(_priors[j]->pdf(input_matrix(i, j))) -
101 2440 : std::log(_priors[j]->pdf(_data_prev(i, j))));
102 3660 : for (unsigned int j = 0; j < _num_confg_values; ++j)
103 : {
104 2440 : out1[j] = (*_outputs_required)[j * _props + i];
105 2440 : out2[j] = _outputs_prev[j * _props + i];
106 : }
107 1220 : if (_var_prior)
108 : {
109 300 : evidence[i] += (std::log(_var_prior->pdf(_new_var_samples[i])) -
110 300 : std::log(_var_prior->pdf(_var_prev[i])));
111 300 : _noise = std::sqrt(_new_var_samples[i]);
112 600 : for (unsigned int j = 0; j < _likelihoods.size(); ++j)
113 300 : evidence[i] += _likelihoods[j]->function(out1);
114 300 : _noise = std::sqrt(_var_prev[i]);
115 600 : for (unsigned int j = 0; j < _likelihoods.size(); ++j)
116 300 : evidence[i] -= _likelihoods[j]->function(out2);
117 : }
118 : else
119 1840 : for (unsigned int j = 0; j < _likelihoods.size(); ++j)
120 920 : evidence[i] += (_likelihoods[j]->function(out1) - _likelihoods[j]->function(out2));
121 : }
122 340 : }
123 :
124 : void
125 160 : PMCMCDecision::computeTransitionVector(std::vector<Real> & tv,
126 : const std::vector<Real> & /*evidence*/)
127 : {
128 160 : tv.assign(_props, 1.0);
129 160 : }
130 :
131 : void
132 1900 : PMCMCDecision::nextSamples(std::vector<Real> & req_inputs,
133 : DenseMatrix<Real> & input_matrix,
134 : const std::vector<Real> & tv,
135 : const unsigned int & parallel_index)
136 : {
137 1900 : if (tv[parallel_index] >= _rnd_vec[parallel_index])
138 : {
139 4020 : for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
140 2680 : req_inputs[k] = input_matrix(parallel_index, k);
141 1340 : _variance[parallel_index] = _new_var_samples[parallel_index];
142 : }
143 : else
144 : {
145 1680 : for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
146 : {
147 1120 : req_inputs[k] = _data_prev(parallel_index, k);
148 1120 : input_matrix(parallel_index, k) = _data_prev(parallel_index, k);
149 : }
150 560 : if (_var_prior)
151 220 : _variance[parallel_index] = _var_prev[parallel_index];
152 1680 : for (unsigned int k = 0; k < _num_confg_values; ++k)
153 1120 : (*_outputs_required)[k * _props + parallel_index] =
154 1120 : _outputs_prev[k * _props + parallel_index];
155 : }
156 1900 : }
157 :
158 : void
159 600 : PMCMCDecision::execute()
160 : {
161 600 : if (_sampler.getNumberOfLocalRows() == 0 || _check_step == _t_step)
162 : {
163 0 : _check_step = _t_step;
164 0 : return;
165 : }
166 :
167 : // Gather inputs and outputs from the sampler and subApps
168 600 : DenseMatrix<Real> data_in(_sampler.getNumberOfRows(), _sampler.getNumberOfCols());
169 1800 : for (dof_id_type ss = _sampler.getLocalRowBegin(); ss < _sampler.getLocalRowEnd(); ++ss)
170 : {
171 1200 : const auto data = _sampler.getNextLocalRow();
172 4800 : for (unsigned int j = 0; j < _sampler.getNumberOfCols(); ++j)
173 3600 : data_in(ss, j) = data[j];
174 1200 : }
175 600 : _local_comm.sum(data_in.get_values());
176 600 : if (!usingGP())
177 : {
178 600 : (*_outputs_required) = *_output_value;
179 600 : _local_comm.allgather((*_outputs_required));
180 : }
181 :
182 : // Compute the evidence and transitimkon vectors
183 600 : std::vector<Real> evidence(_props);
184 600 : if (_t_step > _pmcmc->decisionStep())
185 : {
186 400 : computeEvidence(evidence, data_in);
187 400 : computeTransitionVector(_tpm, evidence);
188 : }
189 : else
190 200 : _tpm.assign(_props, 1.0);
191 :
192 : // Accept/reject the proposed samples and assign the correct outputs
193 600 : std::vector<Real> req_inputs(_sampler.getNumberOfCols() - _num_confg_params);
194 3000 : for (unsigned int i = 0; i < _props; ++i)
195 : {
196 2400 : nextSamples(req_inputs, data_in, _tpm, i);
197 2400 : _inputs[i] = req_inputs;
198 : }
199 :
200 : // Compute the next seeds to facilitate proposals (not always required)
201 600 : nextSeeds();
202 :
203 : // Store data from previous step
204 600 : _data_prev = data_in;
205 600 : _var_prev = _variance;
206 600 : if (!usingGP())
207 600 : _outputs_prev = (*_outputs_required);
208 :
209 : // Track the current step
210 600 : _check_step = _t_step;
211 600 : }
|
