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          76 : AdaptiveImportanceSampler::validParams()
      20             : {
      21          76 :   InputParameters params = Sampler::validParams();
      22          76 :   params.addClassDescription("Adaptive Importance Sampler.");
      23         152 :   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         152 :   params.addRequiredParam<ReporterName>("inputs_reporter", "Reporter with input parameters.");
      28         152 :   params.addRequiredParam<std::vector<Real>>("proposal_std",
      29             :                                              "Standard deviations of the proposal distributions");
      30         152 :   params.addRequiredParam<Real>("output_limit", "Limiting values of the VPPs");
      31         152 :   params.addRequiredParam<std::vector<Real>>(
      32             :       "initial_values", "Initial input values to get the importance sampler started");
      33         152 :   params.addRequiredRangeCheckedParam<int>(
      34             :       "num_samples_train",
      35             :       "num_samples_train>0",
      36             :       "Number of samples to learn the importance distribution");
      37         152 :   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         152 :   params.addRequiredParam<Real>(
      42             :       "std_factor", "Factor to be multiplied to the standard deviation of the importance samples");
      43         152 :   params.addParam<bool>("use_absolute_value", false, "Use absolute value of the sub app output");
      44         152 :   params.addParam<unsigned int>(
      45             :       "num_random_seeds",
      46         152 :       100000,
      47             :       "Initialize a certain number of random seeds. Change from the default only if you have to.");
      48         152 :   params.addParam<ReporterName>("flag_sample",
      49             :                                 "Flag samples if the surrogate prediction was inadequate.");
      50          76 :   return params;
      51           0 : }
      52             : 
      53          44 : AdaptiveImportanceSampler::AdaptiveImportanceSampler(const InputParameters & parameters)
      54             :   : Sampler(parameters),
      55             :     TransientInterface(this),
      56          88 :     _proposal_std(getParam<std::vector<Real>>("proposal_std")),
      57          88 :     _initial_values(getParam<std::vector<Real>>("initial_values")),
      58          88 :     _output_limit(getParam<Real>("output_limit")),
      59          88 :     _num_samples_train(getParam<int>("num_samples_train")),
      60          88 :     _num_importance_sampling_steps(getParam<int>("num_importance_sampling_steps")),
      61          88 :     _std_factor(getParam<Real>("std_factor")),
      62          88 :     _use_absolute_value(getParam<bool>("use_absolute_value")),
      63          88 :     _num_random_seeds(getParam<unsigned int>("num_random_seeds")),
      64          44 :     _is_sampling_completed(false),
      65          44 :     _inputs(getReporterValue<std::vector<std::vector<Real>>>("inputs_reporter")),
      66          44 :     _retraining_steps(0),
      67         110 :     _gp_flag(isParamValid("flag_sample") ? &getReporterValue<std::vector<bool>>("flag_sample")
      68          88 :                                          : nullptr)
      69             : {
      70             :   // Filling the `distributions` vector with the user-provided distributions.
      71         176 :   for (const DistributionName & name : getParam<std::vector<DistributionName>>("distributions"))
      72          88 :     _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          44 :   setNumberOfRows(1);
      80             : 
      81             :   // Setting the number of columns in the sampler matrix (equal to the number of distributions).
      82          44 :   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          44 :   _inputs_sto.resize(_distributions.size());
      88             : 
      89             :   // Mapping all the input distributions to a standard normal space
      90         132 :   for (unsigned int i = 0; i < _distributions.size(); ++i)
      91          88 :     _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          44 :   _prev_value.resize(_distributions.size());
      96             : 
      97             :   // `check_step` is a member variable for ensuring that the MCMC algorithm proceeds in a sequential
      98             :   // fashion.
      99          44 :   _check_step = 0;
     100             : 
     101             :   // Storage for means of input values for proposing the next sample
     102          44 :   _mean_sto.resize(_distributions.size());
     103             : 
     104             :   // Storage for standard deviations of input values for proposing the next sample
     105          44 :   _std_sto.resize(_distributions.size());
     106             : 
     107          44 :   setNumberOfRandomSeeds(_num_random_seeds);
     108          44 : }
     109             : 
     110             : Real
     111        4500 : AdaptiveImportanceSampler::computeSample(dof_id_type /*row_index*/, dof_id_type col_index)
     112             : {
     113        4500 :   const bool sample = _t_step > 1 && col_index == 0 && _check_step != _t_step;
     114        4500 :   const bool gp_flag = _gp_flag ? (*_gp_flag)[0] : false;
     115             : 
     116        4500 :   if (sample && _is_sampling_completed)
     117           0 :     mooseError("Internal bug: the adaptive sampling is supposed to be completed but another sample "
     118             :                "has been requested.");
     119             : 
     120        4500 :   if (_t_step <= _num_samples_train)
     121             :   {
     122             :     /* This is the importance distribution training step. Markov Chains are set up
     123             :       to sample from the importance region or the failure region using the Metropolis
     124             :       algorithm. Given that the previous sample resulted in a model failure, the next
     125             :       sample is proposed such that it is very likely to result in a model failure as well.
     126             :       The `initial_values` and `proposal_std` parameters provided by the user affects the
     127             :       formation of the importance distribution. */
     128        2120 :     if (sample && !gp_flag)
     129             :     {
     130        1290 :       for (dof_id_type j = 0; j < _distributions.size(); ++j)
     131         860 :         _prev_value[j] = Normal::quantile(_distributions[j]->cdf(_inputs[j][0]), 0.0, 1.0);
     132             :       Real acceptance_ratio = 0.0;
     133        1290 :       for (dof_id_type i = 0; i < _distributions.size(); ++i)
     134         860 :         acceptance_ratio += std::log(Normal::pdf(_prev_value[i], 0.0, 1.0)) -
     135         860 :                             std::log(Normal::pdf(_inputs_sto[i].back(), 0.0, 1.0));
     136         430 :       if (acceptance_ratio > std::log(getRand(_t_step)))
     137             :       {
     138         750 :         for (dof_id_type i = 0; i < _distributions.size(); ++i)
     139         500 :           _inputs_sto[i].push_back(_prev_value[i]);
     140             :       }
     141             :       else
     142             :       {
     143         540 :         for (dof_id_type i = 0; i < _distributions.size(); ++i)
     144         360 :           _inputs_sto[i].push_back(_inputs_sto[i].back());
     145             :       }
     146        1290 :       for (dof_id_type i = 0; i < _distributions.size(); ++i)
     147         860 :         _prev_value[i] =
     148         860 :             Normal::quantile(getRand(_t_step), _inputs_sto[i].back(), _proposal_std[i]);
     149             :     }
     150             :   }
     151        2380 :   else if (sample && !gp_flag)
     152             :   {
     153             :     /* This is the importance sampling step using the importance distribution created
     154             :       in the previous step. Once the importance distribution is known, sampling from
     155             :       it is similar to a regular Monte Carlo sampling. */
     156        1050 :     for (dof_id_type i = 0; i < _distributions.size(); ++i)
     157             :     {
     158         700 :       if (_t_step == _num_samples_train + 1)
     159             :       {
     160          40 :         _mean_sto[i] = AdaptiveMonteCarloUtils::computeMean(_inputs_sto[i], 1);
     161          40 :         _std_sto[i] = AdaptiveMonteCarloUtils::computeSTD(_inputs_sto[i], 1);
     162             :       }
     163         700 :       _prev_value[i] =
     164         700 :           (Normal::quantile(getRand(_t_step), _mean_sto[i], _std_factor * _std_sto[i]));
     165             :     }
     166             : 
     167             :     // check if we have performed all the importance sampling steps
     168         350 :     if (_t_step >= _num_samples_train + _num_importance_sampling_steps + _retraining_steps)
     169          20 :       _is_sampling_completed = true;
     170             :   }
     171             : 
     172             :   // When the GP fails, the current time step is 'wasted' and the retraining step doesn't
     173             :   // happen until the next time step. Therefore, keep track of the number of retraining steps
     174             :   // to increase the total number of steps taken.
     175        4500 :   if (sample && gp_flag && _t_step > _num_samples_train)
     176          10 :     ++_retraining_steps;
     177             : 
     178        4500 :   _check_step = _t_step;
     179        4500 :   return _distributions[col_index]->quantile(Normal::cdf(_prev_value[col_index], 0.0, 1.0));
     180             : }