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 "AdaptiveImportanceSampler.h"
      11             : #include "AdaptiveMonteCarloUtils.h"
      12             : #include "Distribution.h"
      13             : #include "Normal.h"
      14             : #include "Uniform.h"
      15             : 
      16             : registerMooseObjectAliased("StochasticToolsApp", AdaptiveImportanceSampler, "AdaptiveImportance");
      17             : 
      18             : InputParameters
      19          30 : AdaptiveImportanceSampler::validParams()
      20             : {
      21          30 :   InputParameters params = Sampler::validParams();
      22          30 :   params.addClassDescription("Adaptive Importance Sampler.");
      23          60 :   params.addRequiredParam<std::vector<DistributionName>>(
      24             :       "distributions",
      25             :       "The distribution names to be sampled, the number of distributions provided defines the "
      26             :       "number of columns per matrix.");
      27          60 :   params.addRequiredParam<ReporterName>("inputs_reporter", "Reporter with input parameters.");
      28          60 :   params.addRequiredParam<std::vector<Real>>("proposal_std",
      29             :                                              "Standard deviations of the proposal distributions");
      30          60 :   params.addRequiredParam<Real>("output_limit", "Limiting values of the VPPs");
      31          60 :   params.addRequiredParam<std::vector<Real>>(
      32             :       "initial_values", "Initial input values to get the importance sampler started");
      33          60 :   params.addRequiredRangeCheckedParam<int>(
      34             :       "num_samples_train",
      35             :       "num_samples_train>0",
      36             :       "Number of samples to learn the importance distribution");
      37          60 :   params.addRequiredRangeCheckedParam<int>(
      38             :       "num_importance_sampling_steps",
      39             :       "num_importance_sampling_steps>0",
      40             :       "Number of importance sampling steps (after the importance distribution has been trained)");
      41          60 :   params.addRequiredParam<Real>(
      42             :       "std_factor", "Factor to be multiplied to the standard deviation of the importance samples");
      43          60 :   params.addParam<bool>("use_absolute_value", false, "Use absolute value of the sub app output");
      44          60 :   params.addParam<unsigned int>(
      45             :       "num_random_seeds",
      46          60 :       100000,
      47             :       "Initialize a certain number of random seeds. Change from the default only if you have to.");
      48          60 :   params.addParam<ReporterName>("flag_sample",
      49             :                                 "Flag samples if the surrogate prediction was inadequate.");
      50          30 :   return params;
      51           0 : }
      52             : 
      53          16 : AdaptiveImportanceSampler::AdaptiveImportanceSampler(const InputParameters & parameters)
      54             :   : Sampler(parameters),
      55             :     TransientInterface(this),
      56          32 :     _proposal_std(getParam<std::vector<Real>>("proposal_std")),
      57          32 :     _initial_values(getParam<std::vector<Real>>("initial_values")),
      58          32 :     _output_limit(getParam<Real>("output_limit")),
      59          32 :     _num_samples_train(getParam<int>("num_samples_train")),
      60          32 :     _num_importance_sampling_steps(getParam<int>("num_importance_sampling_steps")),
      61          32 :     _std_factor(getParam<Real>("std_factor")),
      62          32 :     _use_absolute_value(getParam<bool>("use_absolute_value")),
      63          32 :     _num_random_seeds(getParam<unsigned int>("num_random_seeds")),
      64          16 :     _is_sampling_completed(false),
      65          16 :     _inputs(getReporterValue<std::vector<std::vector<Real>>>("inputs_reporter")),
      66          16 :     _retraining_steps(0),
      67          40 :     _gp_flag(isParamValid("flag_sample") ? &getReporterValue<std::vector<bool>>("flag_sample")
      68          16 :                                          : nullptr)
      69             : {
      70             :   // Filling the `distributions` vector with the user-provided distributions.
      71          64 :   for (const DistributionName & name : getParam<std::vector<DistributionName>>("distributions"))
      72          32 :     _distributions.push_back(&getDistributionByName(name));
      73             : 
      74             :   /* Adaptive Importance Sampling (AdaptiveImportanceSampler) relies on a Markov Chain Monte Carlo
      75             :      (MCMC) algorithm. As such, in MOOSE, any use of MCMC algorithms requires that the `num_steps`
      76             :      parameter in the main App's executioner would control the total number of samples. Therefore,
      77             :      the `num_rows` parameter typically used by exisiting non-MCMC samplers to set the total number
      78             :      of samples has no use here and is fixed to 1.*/
      79          16 :   setNumberOfRows(1);
      80             : 
      81             :   // Setting the number of columns in the sampler matrix (equal to the number of distributions).
      82          16 :   setNumberOfCols(_distributions.size());
      83             : 
      84             :   /* `inputs_sto` is a member variable that aids in forming the importance distribution.
      85             :      One dimension of this variable is equal to the number of distributions. The other dimension
      86             :      of the variable, at the last step, is equal to the number of samples the user desires.*/
      87          16 :   _inputs_sto.resize(_distributions.size());
      88             : 
      89             :   // Mapping all the input distributions to a standard normal space
      90          48 :   for (unsigned int i = 0; i < _distributions.size(); ++i)
      91          32 :     _inputs_sto[i].push_back(Normal::quantile(_distributions[i]->cdf(_initial_values[i]), 0, 1));
      92             : 
      93             :   /* `prev_value` is a member variable for tracking the previously accepted samples in the
      94             :      MCMC algorithm and proposing the next sample.*/
      95          16 :   _prev_value.resize(_distributions.size());
      96             : 
      97             :   // Storage for means of input values for proposing the next sample
      98          16 :   _mean_sto.resize(_distributions.size());
      99             : 
     100             :   // Storage for standard deviations of input values for proposing the next sample
     101          16 :   _std_sto.resize(_distributions.size());
     102             : 
     103          16 :   setNumberOfRandomSeeds(_num_random_seeds);
     104          16 :   setAutoAdvanceGenerators(false);
     105          16 : }
     106             : 
     107             : void
     108         664 : AdaptiveImportanceSampler::executeSetUp()
     109             : {
     110         664 :   const bool sample = _t_step > 1;
     111         664 :   const bool gp_flag = _gp_flag ? (*_gp_flag)[0] : false;
     112             : 
     113         664 :   if (sample && _is_sampling_completed)
     114           0 :     mooseError("Internal bug: the adaptive sampling is supposed to be completed but another sample "
     115             :                "has been requested.");
     116             : 
     117         664 :   if (_t_step <= _num_samples_train)
     118             :   {
     119             :     /* This is the importance distribution training step. Markov Chains are set up
     120             :       to sample from the importance region or the failure region using the Metropolis
     121             :       algorithm. Given that the previous sample resulted in a model failure, the next
     122             :       sample is proposed such that it is very likely to result in a model failure as well.
     123             :       The `initial_values` and `proposal_std` parameters provided by the user affects the
     124             :       formation of the importance distribution. */
     125         376 :     if (sample && !gp_flag)
     126             :     {
     127        1032 :       for (dof_id_type j = 0; j < _distributions.size(); ++j)
     128         688 :         _prev_value[j] = Normal::quantile(_distributions[j]->cdf(_inputs[j][0]), 0.0, 1.0);
     129             :       Real acceptance_ratio = 0.0;
     130        1032 :       for (dof_id_type i = 0; i < _distributions.size(); ++i)
     131         688 :         acceptance_ratio += std::log(Normal::pdf(_prev_value[i], 0.0, 1.0)) -
     132         688 :                             std::log(Normal::pdf(_inputs_sto[i].back(), 0.0, 1.0));
     133         344 :       if (acceptance_ratio > std::log(getRand(0, _t_step)))
     134             :       {
     135         678 :         for (dof_id_type i = 0; i < _distributions.size(); ++i)
     136         452 :           _inputs_sto[i].push_back(_prev_value[i]);
     137             :       }
     138             :       else
     139             :       {
     140         354 :         for (dof_id_type i = 0; i < _distributions.size(); ++i)
     141         236 :           _inputs_sto[i].push_back(_inputs_sto[i].back());
     142             :       }
     143        1032 :       for (dof_id_type i = 0; i < _distributions.size(); ++i)
     144         688 :         _prev_value[i] =
     145         688 :             Normal::quantile(getRand(i + 1, _t_step), _inputs_sto[i].back(), _proposal_std[i]);
     146             :     }
     147             :   }
     148         288 :   else if (sample && !gp_flag)
     149             :   {
     150             :     /* This is the importance sampling step using the importance distribution created
     151             :       in the previous step. Once the importance distribution is known, sampling from
     152             :       it is similar to a regular Monte Carlo sampling. */
     153         840 :     for (dof_id_type i = 0; i < _distributions.size(); ++i)
     154             :     {
     155         560 :       if (_t_step == _num_samples_train + 1)
     156             :       {
     157          32 :         _mean_sto[i] = AdaptiveMonteCarloUtils::computeMean(_inputs_sto[i], 1);
     158          32 :         _std_sto[i] = AdaptiveMonteCarloUtils::computeSTD(_inputs_sto[i], 1);
     159             :       }
     160         560 :       _prev_value[i] =
     161         560 :           (Normal::quantile(getRand(i, _t_step), _mean_sto[i], _std_factor * _std_sto[i]));
     162             :     }
     163             : 
     164             :     // check if we have performed all the importance sampling steps
     165         280 :     if (_t_step >= _num_samples_train + _num_importance_sampling_steps + _retraining_steps)
     166          16 :       _is_sampling_completed = true;
     167             :   }
     168             : 
     169             :   // When the GP fails, the current time step is 'wasted' and the retraining step doesn't
     170             :   // happen until the next time step. Therefore, keep track of the number of retraining steps
     171             :   // to increase the total number of steps taken.
     172         664 :   if (sample && gp_flag && _t_step > _num_samples_train)
     173           8 :     ++_retraining_steps;
     174         664 : }
     175             : 
     176             : Real
     177        2250 : AdaptiveImportanceSampler::computeSample(dof_id_type /*row_index*/, dof_id_type col_index)
     178             : {
     179        2250 :   return _distributions[col_index]->quantile(Normal::cdf(_prev_value[col_index], 0.0, 1.0));
     180             : }