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 "AdaptiveImportanceStats.h"
      11             : #include "Sampler.h"
      12             : #include "Distribution.h"
      13             : #include "Normal.h"
      14             : #include "libmesh/utility.h"
      15             : 
      16             : registerMooseObject("StochasticToolsApp", AdaptiveImportanceStats);
      17             : 
      18             : InputParameters
      19          68 : AdaptiveImportanceStats::validParams()
      20             : {
      21          68 :   InputParameters params = GeneralReporter::validParams();
      22          68 :   params.addClassDescription("Reporter to compute statistics corresponding to the "
      23             :                              "AdaptiveImportanceSampler.");
      24         136 :   params.addRequiredParam<ReporterName>("output_value",
      25             :                                         "Value of the model output from the SubApp.");
      26         136 :   params.addParam<ReporterValueName>("mu_imp", "mu_imp", "Means of the importance distributions.");
      27         136 :   params.addParam<ReporterValueName>(
      28             :       "std_imp", "std_imp", "Standard deviations of the importance distributions.");
      29         136 :   params.addParam<ReporterValueName>("pf", "pf", "Failure probability estimate.");
      30         136 :   params.addParam<ReporterValueName>(
      31             :       "cov_pf", "cov_pf", "Coefficient of variation of failure probability.");
      32         136 :   params.addParam<ReporterName>("flag_sample",
      33             :                                 "Flag samples if the surrogate prediction was inadequate.");
      34         136 :   params.addRequiredParam<SamplerName>("sampler", "The sampler object.");
      35          68 :   return params;
      36           0 : }
      37             : 
      38          34 : AdaptiveImportanceStats::AdaptiveImportanceStats(const InputParameters & parameters)
      39             :   : GeneralReporter(parameters),
      40          68 :     _output_value(isParamValid("flag_sample") ? getReporterValue<std::vector<Real>>("output_value")
      41          68 :                                               : getReporterValue<std::vector<Real>>(
      42             :                                                     "output_value", REPORTER_MODE_DISTRIBUTED)),
      43          34 :     _mu_imp(declareValue<std::vector<Real>>("mu_imp")),
      44          34 :     _std_imp(declareValue<std::vector<Real>>("std_imp")),
      45          34 :     _pf(declareValue<std::vector<Real>>("pf")),
      46          34 :     _cov_pf(declareValue<std::vector<Real>>("cov_pf")),
      47          68 :     _step(getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")->timeStep()),
      48          34 :     _ais(getSampler<AdaptiveImportanceSampler>("sampler")),
      49          85 :     _gp_flag(isParamValid("flag_sample") ? &getReporterValue<std::vector<bool>>("flag_sample")
      50             :                                          : nullptr),
      51          68 :     _check_step(std::numeric_limits<int>::max())
      52             : {
      53             :   // Initialize variables
      54          34 :   const auto rows = _ais.getNumberOfRows();
      55          34 :   _mu_imp.resize(rows);
      56          34 :   _std_imp.resize(rows);
      57          34 :   _pf.resize(1);
      58          34 :   _cov_pf.resize(1);
      59          34 :   _pf_sum = 0.0;
      60          34 :   _var_sum = 0.0;
      61          34 :   _distributions_store = _ais.getDistributionNames();
      62          34 :   _factor = _ais.getStdFactor();
      63          34 : }
      64             : 
      65             : void
      66        1377 : AdaptiveImportanceStats::execute()
      67             : {
      68        1377 :   if (_ais.getNumberOfLocalRows() == 0 || _check_step == _step)
      69             :   {
      70         486 :     _check_step = _step;
      71         486 :     return;
      72             :   }
      73             : 
      74         891 :   const bool gp_flag = _gp_flag ? (*_gp_flag)[0] : false;
      75             :   // Compute AdaptiveImportanceSampler statistics at each sample during the evaluation phase only.
      76         891 :   if (_step > _ais.getNumSamplesTrain() && !gp_flag)
      77             :   {
      78             :     // Get the statistics of the importance distributions in the standard Normal space.
      79         385 :     _mu_imp = _ais.getImportanceVectorMean();
      80         385 :     _std_imp = _ais.getImportanceVectorStd();
      81             : 
      82             :     // Get the failure probability estimate.
      83         385 :     const Real tmp = _ais.getUseAbsoluteValue() ? std::abs(_output_value[0]) : _output_value[0];
      84         385 :     const bool output_limit_reached = tmp >= _ais.getOutputLimit();
      85         385 :     Real prod1 = output_limit_reached ? 1.0 : 0.0;
      86         385 :     std::vector<Real> input1 = _ais.getNextLocalRow();
      87         385 :     Real input_tmp = 0.0;
      88        1155 :     for (dof_id_type ss = 0; ss < _ais.getNumberOfCols(); ++ss)
      89             :     {
      90         770 :       input_tmp = Normal::quantile(_distributions_store[ss]->cdf(input1[ss]), 0.0, 1.0);
      91         770 :       prod1 = prod1 * (Normal::pdf(input_tmp, 0.0, 1.0) /
      92         770 :                        Normal::pdf(input_tmp, _mu_imp[ss], _factor * _std_imp[ss]));
      93             :     }
      94         385 :     _pf_sum += prod1;
      95         385 :     _var_sum += Utility::pow<2>(prod1);
      96         385 :     _pf[0] = _pf_sum / (_step - _ais.getNumSamplesTrain());
      97             : 
      98             :     // Get coefficient of variation of failure probability.
      99             :     Real tmp_var =
     100         385 :         std::sqrt(1.0 / (_step - _ais.getNumSamplesTrain()) *
     101         385 :                   (_var_sum / (_step - _ais.getNumSamplesTrain()) - Utility::pow<2>(_pf[0])));
     102         385 :     _cov_pf[0] = tmp_var / _pf[0];
     103         385 :   }
     104             : }