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