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