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 "AdaptiveMonteCarloDecision.h" 11 : #include "Sampler.h" 12 : #include "DenseMatrix.h" 13 : #include "AdaptiveMonteCarloUtils.h" 14 : #include "StochasticToolsUtils.h" 15 : 16 : registerMooseObject("StochasticToolsApp", AdaptiveMonteCarloDecision); 17 : 18 : InputParameters 19 104 : AdaptiveMonteCarloDecision::validParams() 20 : { 21 104 : InputParameters params = GeneralReporter::validParams(); 22 104 : params.addClassDescription("Generic reporter which decides whether or not to accept a proposed " 23 : "sample in Adaptive Monte Carlo type of algorithms."); 24 208 : params.addRequiredParam<ReporterName>("output_value", 25 : "Value of the model output from the SubApp."); 26 208 : params.addParam<ReporterValueName>( 27 : "output_required", 28 : "output_required", 29 : "Modified value of the model output from this reporter class."); 30 208 : params.addParam<ReporterValueName>("inputs", "inputs", "Uncertain inputs to the model."); 31 208 : params.addRequiredParam<SamplerName>("sampler", "The sampler object."); 32 208 : params.addParam<UserObjectName>("gp_decision", "The Gaussian Process decision reporter."); 33 104 : return params; 34 0 : } 35 : 36 56 : AdaptiveMonteCarloDecision::AdaptiveMonteCarloDecision(const InputParameters & parameters) 37 : : GeneralReporter(parameters), 38 112 : _output_value(isParamValid("gp_decision") ? getReporterValue<std::vector<Real>>("output_value") 39 136 : : getReporterValue<std::vector<Real>>( 40 : "output_value", REPORTER_MODE_DISTRIBUTED)), 41 56 : _output_required(declareValue<std::vector<Real>>("output_required")), 42 56 : _inputs(declareValue<std::vector<std::vector<Real>>>("inputs")), 43 56 : _sampler(getSampler("sampler")), 44 56 : _ais(dynamic_cast<const AdaptiveImportanceSampler *>(&_sampler)), 45 56 : _pss(dynamic_cast<const ParallelSubsetSimulation *>(&_sampler)), 46 56 : _check_step(std::numeric_limits<int>::max()), 47 56 : _local_comm(_sampler.getLocalComm()), 48 112 : _gp_used(isParamValid("gp_decision")), 49 56 : _gp_training_samples( 50 72 : _gp_used ? &getUserObject<ActiveLearningGPDecision>("gp_decision").getTrainingSamples() 51 112 : : nullptr) 52 : { 53 : 54 : // Check whether the selected sampler is an adaptive sampler or not 55 56 : if (!_ais && !_pss) 56 8 : paramError("sampler", "The selected sampler is not an adaptive sampler."); 57 : 58 48 : const auto rows = _sampler.getNumberOfRows(); 59 48 : const auto cols = _sampler.getNumberOfCols(); 60 : 61 : // Initialize the required variables depending upon the type of adaptive Monte Carlo algorithm 62 48 : _inputs.resize(cols, std::vector<Real>(rows)); 63 48 : _prev_val.resize(cols, std::vector<Real>(rows)); 64 48 : _output_required.resize(rows); 65 48 : _prev_val_out.resize(rows); 66 : 67 48 : if (_ais) 68 : { 69 96 : for (dof_id_type j = 0; j < _sampler.getNumberOfCols(); ++j) 70 64 : _prev_val[j][0] = _ais->getInitialValues()[j]; 71 32 : _output_limit = _ais->getOutputLimit(); 72 : } 73 16 : else if (_pss) 74 : { 75 16 : _inputs_sto.resize(cols, std::vector<Real>(_pss->getNumSamplesSub())); 76 16 : _inputs_sorted.resize(cols); 77 16 : _outputs_sto.resize(_pss->getNumSamplesSub()); 78 16 : _output_limit = -std::numeric_limits<Real>::max(); 79 : } 80 48 : } 81 : 82 : void 83 10 : AdaptiveMonteCarloDecision::reinitChain() 84 : { 85 10 : const std::vector<Real> & tmp1 = _ais->getInitialValues(); 86 30 : for (dof_id_type j = 0; j < tmp1.size(); ++j) 87 20 : _inputs[j][0] = tmp1[j]; 88 10 : _prev_val = _inputs; 89 10 : _prev_val_out[0] = 1.0; 90 10 : } 91 : 92 : void 93 1392 : AdaptiveMonteCarloDecision::execute() 94 : { 95 1392 : if (_sampler.getNumberOfLocalRows() == 0 || _check_step == _t_step) 96 : { 97 486 : _check_step = _t_step; 98 486 : return; 99 : } 100 : 101 : /* Decision step to whether or not to accept the proposed sample by the sampler. 102 : This decision step changes with the type of adaptive Monte Carlo sampling algorithm. */ 103 906 : if (_ais) 104 : { 105 810 : const Real tmp = _ais->getUseAbsoluteValue() ? std::abs(_output_value[0]) : _output_value[0]; 106 : 107 : /* Checking whether a GP surrogate is used. If it is used, importance sampling is not performed 108 : during the training phase of the GP and all proposed samples are accepted until the training 109 : phase is completed. Once the training is completed, the importance sampling starts. 110 : If a GP surrogate is not used, the standard proposal and acceptance/rejection is performed as 111 : part of the importance sampling. */ 112 810 : const bool restart_gp = _gp_used && _t_step == *_gp_training_samples; 113 810 : const bool output_limit_reached = _gp_used || tmp >= _output_limit; 114 810 : if (restart_gp) 115 10 : reinitChain(); 116 : 117 1100 : _output_required[0] = output_limit_reached ? 1.0 : 0.0; 118 : 119 810 : if (_t_step <= _ais->getNumSamplesTrain() && !restart_gp) 120 : { 121 : /* This is the training phase of the Adaptive Importance Sampling algorithm. 122 : Here, it is decided whether or not to accept a proposed sample by the 123 : AdaptiveImportanceSampler.C sampler depending upon the model output_value. */ 124 440 : _inputs = output_limit_reached 125 880 : ? StochasticTools::reshapeVector(_sampler.getNextLocalRow(), 1, true) 126 150 : : _prev_val; 127 440 : if (output_limit_reached) 128 290 : _prev_val = _inputs; 129 440 : _prev_val_out = _output_required; 130 : } 131 370 : else if (_t_step > _ais->getNumSamplesTrain() && !restart_gp) 132 : { 133 : /* This is the sampling phase of the Adaptive Importance Sampling algorithm. 134 : Here, all proposed samples by the AdaptiveImportanceSampler.C sampler are accepted since 135 : the importance distribution traning phase is finished. */ 136 720 : _inputs = StochasticTools::reshapeVector(_sampler.getNextLocalRow(), 1, true); 137 360 : _prev_val_out[0] = tmp; 138 : } 139 : } 140 96 : else if (_pss) 141 : { 142 : // Track the current subset 143 : const unsigned int subset = 144 96 : ((_t_step - 1) * _sampler.getNumberOfRows()) / _pss->getNumSamplesSub(); 145 : const unsigned int sub_ind = 146 96 : (_t_step - 1) - (_pss->getNumSamplesSub() / _sampler.getNumberOfRows()) * subset; 147 96 : const unsigned int offset = sub_ind * _sampler.getNumberOfRows(); 148 96 : const unsigned int count_max = 1 / _pss->getSubsetProbability(); 149 : 150 96 : DenseMatrix<Real> data_in(_sampler.getNumberOfRows(), _sampler.getNumberOfCols()); 151 216 : for (dof_id_type ss = _sampler.getLocalRowBegin(); ss < _sampler.getLocalRowEnd(); ++ss) 152 : { 153 120 : const auto data = _sampler.getNextLocalRow(); 154 360 : for (unsigned int j = 0; j < _sampler.getNumberOfCols(); ++j) 155 240 : data_in(ss, j) = data[j]; 156 : } 157 96 : _local_comm.sum(data_in.get_values()); 158 : 159 : // Get the accepted samples outputs across all the procs from the previous step 160 96 : _output_required = (_pss->getUseAbsoluteValue()) 161 96 : ? AdaptiveMonteCarloUtils::computeVectorABS(_output_value) 162 96 : : _output_value; 163 96 : _local_comm.allgather(_output_required); 164 : 165 : // These are the subsequent subsets which use Markov Chain Monte Carlo sampling scheme 166 96 : if (subset > 0) 167 : { 168 64 : if (sub_ind == 0) 169 : { 170 : // _output_sorted contains largest po percentile output values 171 96 : _output_sorted = AdaptiveMonteCarloUtils::sortOutput( 172 32 : _outputs_sto, _pss->getNumSamplesSub(), _pss->getSubsetProbability()); 173 : // _inputs_sorted contains the input values corresponding to the largest po percentile 174 : // output values 175 32 : _inputs_sorted = AdaptiveMonteCarloUtils::sortInput( 176 32 : _inputs_sto, _outputs_sto, _pss->getNumSamplesSub(), _pss->getSubsetProbability()); 177 : // Get the subset's intermediate failure threshold values 178 32 : _output_limit = AdaptiveMonteCarloUtils::computeMin(_output_sorted); 179 : } 180 : // Check whether the number of samples in a Markov chain exceeded the limit 181 64 : if (sub_ind % count_max == 0) 182 : { 183 32 : const unsigned int soffset = (sub_ind / count_max) * _sampler.getNumberOfRows(); 184 : // Reinitialize the starting input values for the next set of Markov chains 185 96 : for (dof_id_type j = 0; j < _sampler.getNumberOfCols(); ++j) 186 64 : _prev_val[j].assign(_inputs_sorted[j].begin() + soffset, 187 64 : _inputs_sorted[j].begin() + soffset + _sampler.getNumberOfRows()); 188 32 : _prev_val_out.assign(_output_sorted.begin() + soffset, 189 32 : _output_sorted.begin() + soffset + _sampler.getNumberOfRows()); 190 : } 191 : else 192 : { 193 : // Otherwise, use the previously accepted input values to propose the next set of input 194 : // values 195 96 : for (dof_id_type j = 0; j < _sampler.getNumberOfCols(); ++j) 196 64 : _prev_val[j].assign(_inputs_sto[j].begin() + offset - _sampler.getNumberOfRows(), 197 : _inputs_sto[j].begin() + offset); 198 32 : _prev_val_out.assign(_outputs_sto.begin() + offset - _sampler.getNumberOfRows(), 199 : _outputs_sto.begin() + offset); 200 : } 201 : } 202 : 203 : // Check whether the outputs exceed the subset's intermediate failure threshold value 204 288 : for (dof_id_type ss = 0; ss < _sampler.getNumberOfRows(); ++ss) 205 : { 206 : // Check whether the outputs exceed the subset's intermediate failure threshold value 207 : // If so, accept the proposed input values by the Sampler object 208 : // Otherwise, use the previously accepted input values 209 192 : const bool output_limit_reached = _output_required[ss] >= _output_limit; 210 576 : for (dof_id_type i = 0; i < _sampler.getNumberOfCols(); ++i) 211 : { 212 384 : _inputs[i][ss] = output_limit_reached ? data_in(ss, i) : _prev_val[i][ss]; 213 384 : _inputs_sto[i][ss + offset] = _inputs[i][ss]; 214 : } 215 192 : if (!output_limit_reached) 216 80 : _output_required[ss] = _prev_val_out[ss]; 217 192 : _outputs_sto[ss + offset] = _output_required[ss]; 218 : } 219 96 : } 220 : // Track the current step 221 906 : _check_step = _t_step; 222 : }