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