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          46 : AdaptiveMonteCarloDecision::validParams()
      20             : {
      21          46 :   InputParameters params = GeneralReporter::validParams();
      22          46 :   params.addClassDescription("Generic reporter which decides whether or not to accept a proposed "
      23             :                              "sample in Adaptive Monte Carlo type of algorithms.");
      24          92 :   params.addRequiredParam<ReporterName>("output_value",
      25             :                                         "Value of the model output from the SubApp.");
      26          92 :   params.addParam<ReporterValueName>(
      27             :       "output_required",
      28             :       "output_required",
      29             :       "Modified value of the model output from this reporter class.");
      30          92 :   params.addParam<ReporterValueName>("inputs", "inputs", "Uncertain inputs to the model.");
      31          92 :   params.addRequiredParam<SamplerName>("sampler", "The sampler object.");
      32          92 :   params.addParam<UserObjectName>("gp_decision", "The Gaussian Process decision reporter.");
      33          46 :   return params;
      34           0 : }
      35             : 
      36          25 : AdaptiveMonteCarloDecision::AdaptiveMonteCarloDecision(const InputParameters & parameters)
      37             :   : GeneralReporter(parameters),
      38          50 :     _output_value(isParamValid("gp_decision") ? getReporterValue<std::vector<Real>>("output_value")
      39          61 :                                               : getReporterValue<std::vector<Real>>(
      40             :                                                     "output_value", REPORTER_MODE_DISTRIBUTED)),
      41          25 :     _output_required(declareValue<std::vector<Real>>("output_required")),
      42          25 :     _inputs(declareValue<std::vector<std::vector<Real>>>("inputs")),
      43          25 :     _sampler(getSampler("sampler")),
      44          25 :     _ais(dynamic_cast<const AdaptiveImportanceSampler *>(&_sampler)),
      45          25 :     _pss(dynamic_cast<const ParallelSubsetSimulation *>(&_sampler)),
      46          25 :     _check_step(std::numeric_limits<int>::max()),
      47          25 :     _local_comm(_sampler.getLocalComm()),
      48          50 :     _gp_used(isParamValid("gp_decision")),
      49          25 :     _gp_training_samples(
      50          32 :         _gp_used ? &getUserObject<ActiveLearningGPDecision>("gp_decision").getTrainingSamples()
      51          25 :                  : nullptr)
      52             : {
      53             : 
      54             :   // Check whether the selected sampler is an adaptive sampler or not
      55          25 :   if (!_ais && !_pss)
      56           4 :     paramError("sampler", "The selected sampler is not an adaptive sampler.");
      57             : 
      58          21 :   const auto rows = _sampler.getNumberOfRows();
      59          21 :   const auto cols = _sampler.getNumberOfCols();
      60             : 
      61             :   // Initialize the required variables depending upon the type of adaptive Monte Carlo algorithm
      62          21 :   _inputs.resize(cols, std::vector<Real>(rows));
      63          21 :   _prev_val.resize(cols, std::vector<Real>(rows));
      64          21 :   _output_required.resize(rows);
      65          21 :   _prev_val_out.resize(rows);
      66             : 
      67          21 :   if (_ais)
      68             :   {
      69          42 :     for (dof_id_type j = 0; j < _sampler.getNumberOfCols(); ++j)
      70          28 :       _prev_val[j][0] = _ais->getInitialValues()[j];
      71          14 :     _output_limit = _ais->getOutputLimit();
      72             :   }
      73           7 :   else if (_pss)
      74             :   {
      75           7 :     _inputs_sto.resize(cols, std::vector<Real>(_pss->getNumSamplesSub()));
      76           7 :     _inputs_sorted.resize(cols);
      77           7 :     _outputs_sto.resize(_pss->getNumSamplesSub());
      78           7 :     _output_limit = -std::numeric_limits<Real>::max();
      79             :   }
      80          21 : }
      81             : 
      82             : void
      83           5 : AdaptiveMonteCarloDecision::reinitChain()
      84             : {
      85           5 :   const std::vector<Real> & tmp1 = _ais->getInitialValues();
      86          15 :   for (dof_id_type j = 0; j < tmp1.size(); ++j)
      87          10 :     _inputs[j][0] = tmp1[j];
      88           5 :   _prev_val = _inputs;
      89           5 :   _prev_val_out[0] = 1.0;
      90           5 : }
      91             : 
      92             : void
      93         609 : AdaptiveMonteCarloDecision::execute()
      94             : {
      95         609 :   if (_sampler.getNumberOfLocalRows() == 0 || _check_step == _t_step)
      96             :   {
      97         162 :     _check_step = _t_step;
      98         162 :     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         447 :   if (_ais)
     104             :   {
     105         405 :     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         405 :     const bool restart_gp = _gp_used && _t_step == *_gp_training_samples;
     113         405 :     const bool output_limit_reached = _gp_used || tmp >= _output_limit;
     114         405 :     if (restart_gp)
     115           5 :       reinitChain();
     116             : 
     117         550 :     _output_required[0] = output_limit_reached ? 1.0 : 0.0;
     118             : 
     119         405 :     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         220 :       _inputs = output_limit_reached
     125         440 :                     ? StochasticTools::reshapeVector(_sampler.getNextLocalRow(), 1, true)
     126          75 :                     : _prev_val;
     127         220 :       if (output_limit_reached)
     128         145 :         _prev_val = _inputs;
     129         220 :       _prev_val_out = _output_required;
     130             :     }
     131         185 :     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         180 :       _inputs = StochasticTools::reshapeVector(_sampler.getNextLocalRow(), 1, true);
     137         180 :       _prev_val_out[0] = tmp;
     138             :     }
     139             :   }
     140          42 :   else if (_pss)
     141             :   {
     142             :     // Track the current subset
     143             :     const unsigned int subset =
     144          42 :         ((_t_step - 1) * _sampler.getNumberOfRows()) / _pss->getNumSamplesSub();
     145             :     const unsigned int sub_ind =
     146          42 :         (_t_step - 1) - (_pss->getNumSamplesSub() / _sampler.getNumberOfRows()) * subset;
     147          42 :     const unsigned int offset = sub_ind * _sampler.getNumberOfRows();
     148          42 :     const unsigned int count_max = 1 / _pss->getSubsetProbability();
     149             : 
     150          42 :     DenseMatrix<Real> data_in(_sampler.getNumberOfRows(), _sampler.getNumberOfCols());
     151         102 :     for (dof_id_type ss = _sampler.getLocalRowBegin(); ss < _sampler.getLocalRowEnd(); ++ss)
     152             :     {
     153          60 :       const auto data = _sampler.getNextLocalRow();
     154         180 :       for (unsigned int j = 0; j < _sampler.getNumberOfCols(); ++j)
     155         120 :         data_in(ss, j) = data[j];
     156          60 :     }
     157          42 :     _local_comm.sum(data_in.get_values());
     158             : 
     159             :     // Get the accepted samples outputs across all the procs from the previous step
     160          42 :     _output_required = (_pss->getUseAbsoluteValue())
     161          84 :                            ? AdaptiveMonteCarloUtils::computeVectorABS(_output_value)
     162          42 :                            : _output_value;
     163          42 :     _local_comm.allgather(_output_required);
     164             : 
     165             :     // These are the subsequent subsets which use Markov Chain Monte Carlo sampling scheme
     166          42 :     if (subset > 0)
     167             :     {
     168          28 :       if (sub_ind == 0)
     169             :       {
     170             :         // _output_sorted contains largest po percentile output values
     171          14 :         _output_sorted = AdaptiveMonteCarloUtils::sortOutput(
     172          14 :             _outputs_sto, _pss->getNumSamplesSub(), _pss->getSubsetProbability());
     173             :         // _inputs_sorted contains the input values corresponding to the largest po percentile
     174             :         // output values
     175          14 :         _inputs_sorted = AdaptiveMonteCarloUtils::sortInput(
     176          14 :             _inputs_sto, _outputs_sto, _pss->getNumSamplesSub(), _pss->getSubsetProbability());
     177             :         // Get the subset's intermediate failure threshold values
     178          14 :         _output_limit = AdaptiveMonteCarloUtils::computeMin(_output_sorted);
     179             :       }
     180             :       // Check whether the number of samples in a Markov chain exceeded the limit
     181          28 :       if (sub_ind % count_max == 0)
     182             :       {
     183          14 :         const unsigned int soffset = (sub_ind / count_max) * _sampler.getNumberOfRows();
     184             :         // Reinitialize the starting input values for the next set of Markov chains
     185          42 :         for (dof_id_type j = 0; j < _sampler.getNumberOfCols(); ++j)
     186          28 :           _prev_val[j].assign(_inputs_sorted[j].begin() + soffset,
     187          28 :                               _inputs_sorted[j].begin() + soffset + _sampler.getNumberOfRows());
     188          14 :         _prev_val_out.assign(_output_sorted.begin() + soffset,
     189          14 :                              _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          42 :         for (dof_id_type j = 0; j < _sampler.getNumberOfCols(); ++j)
     196          56 :           _prev_val[j].assign(_inputs_sto[j].begin() + offset - _sampler.getNumberOfRows(),
     197             :                               _inputs_sto[j].begin() + offset);
     198          28 :         _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         126 :     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          84 :       const bool output_limit_reached = _output_required[ss] >= _output_limit;
     210         252 :       for (dof_id_type i = 0; i < _sampler.getNumberOfCols(); ++i)
     211             :       {
     212         168 :         _inputs[i][ss] = output_limit_reached ? data_in(ss, i) : _prev_val[i][ss];
     213         168 :         _inputs_sto[i][ss + offset] = _inputs[i][ss];
     214             :       }
     215          84 :       if (!output_limit_reached)
     216          35 :         _output_required[ss] = _prev_val_out[ss];
     217          84 :       _outputs_sto[ss + offset] = _output_required[ss];
     218             :     }
     219          42 :   }
     220             :   // Track the current step
     221         447 :   _check_step = _t_step;
     222             : }