LCOV - code coverage report
Current view: top level - src/reporters - AdaptiveMonteCarloDecision.C (source / functions) Hit Total Coverage
Test: idaholab/moose stochastic_tools: f45d79 Lines: 111 112 99.1 %
Date: 2025-07-25 05:00:46 Functions: 4 4 100.0 %
Legend: Lines: hit not hit

          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             : }

Generated by: LCOV version 1.14