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 "ParallelSubsetSimulation.h"
      11             : #include "AdaptiveMonteCarloUtils.h"
      12             : #include "Normal.h"
      13             : #include "Uniform.h"
      14             : 
      15             : registerMooseObject("StochasticToolsApp", ParallelSubsetSimulation);
      16             : 
      17             : InputParameters
      18          38 : ParallelSubsetSimulation::validParams()
      19             : {
      20          38 :   InputParameters params = Sampler::validParams();
      21          38 :   params.addClassDescription("Parallel Subset Simulation sampler.");
      22          76 :   params.addRequiredParam<std::vector<DistributionName>>(
      23             :       "distributions",
      24             :       "The distribution names to be sampled, the number of distributions provided defines the "
      25             :       "number of columns per matrix.");
      26          76 :   params.addRequiredParam<ReporterName>("output_reporter",
      27             :                                         "Reporter with results of samples created by the SubApp.");
      28          76 :   params.addRequiredParam<ReporterName>("inputs_reporter", "Reporter with input parameters.");
      29         114 :   params.addRangeCheckedParam<Real>("subset_probability",
      30          76 :                                     0.1,
      31             :                                     "subset_probability>0 & subset_probability<=1",
      32             :                                     "Conditional probability of each subset");
      33          76 :   params.addRequiredParam<unsigned int>("num_samplessub", "Number of samples per subset");
      34          76 :   params.addRequiredParam<unsigned int>("num_subsets", "Number of desired subsets");
      35          76 :   params.addParam<unsigned int>("num_parallel_chains",
      36             :                                 "Number of Markov chains to run in parallel, default is based on "
      37             :                                 "the number of processors used.");
      38          76 :   params.addParam<bool>("use_absolute_value", false, "Use absolute value of the sub app output");
      39          76 :   params.addParam<unsigned int>(
      40             :       "num_random_seeds",
      41          76 :       100000,
      42             :       "Initialize a certain number of random seeds. Change from the default only if you have to.");
      43          38 :   return params;
      44           0 : }
      45             : 
      46          22 : ParallelSubsetSimulation::ParallelSubsetSimulation(const InputParameters & parameters)
      47             :   : Sampler(parameters),
      48          22 :     _num_samplessub(getParam<unsigned int>("num_samplessub")),
      49          44 :     _num_subsets(getParam<unsigned int>("num_subsets")),
      50          44 :     _use_absolute_value(getParam<bool>("use_absolute_value")),
      51          44 :     _subset_probability(getParam<Real>("subset_probability")),
      52          44 :     _num_random_seeds(getParam<unsigned int>("num_random_seeds")),
      53          22 :     _outputs(getReporterValue<std::vector<Real>>("output_reporter")),
      54          22 :     _inputs(getReporterValue<std::vector<std::vector<Real>>>("inputs_reporter")),
      55          44 :     _step(getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")->timeStep()),
      56          22 :     _count_max(std::floor(1 / _subset_probability)),
      57          22 :     _check_step(0),
      58          22 :     _subset(0),
      59          44 :     _is_sampling_completed(false)
      60             : {
      61             :   // Fixing the number of rows to the number of processors
      62          22 :   const dof_id_type nchains = isParamValid("num_parallel_chains")
      63          66 :                                   ? getParam<unsigned int>("num_parallel_chains")
      64           0 :                                   : n_processors() / _min_procs_per_row;
      65          22 :   setNumberOfRows(nchains);
      66          22 :   if ((_num_samplessub / nchains) % _count_max > 0)
      67           0 :     mooseError("Number of model evaluations per chain per subset (",
      68           0 :                _num_samplessub / nchains,
      69             :                ") should be a multiple of requested chain length (",
      70           0 :                _count_max,
      71             :                ").");
      72             : 
      73             :   // Filling the `distributions` vector with the user-provided distributions.
      74          88 :   for (const DistributionName & name : getParam<std::vector<DistributionName>>("distributions"))
      75          44 :     _distributions.push_back(&getDistributionByName(name));
      76             : 
      77             :   // Setting the number of columns in the sampler matrix (equal to the number of distributions).
      78          22 :   setNumberOfCols(_distributions.size());
      79             : 
      80             :   /* `inputs_sto` is a member variable that aids in deciding the next set of samples
      81             :   in the Subset Simulation algorithm by storing the input parameter values*/
      82          22 :   _inputs_sto.resize(_distributions.size(), std::vector<Real>(_num_samplessub, 0.0));
      83          22 :   _outputs_sto.resize(_num_samplessub, 0.0);
      84             : 
      85             :   /* `inputs_sorted` is a member variable which also aids in deciding the next set of samples
      86             :   in the Subset Simulation algorithm by storing the sorted input parameter values
      87             :   by their corresponding output values*/
      88          22 :   _inputs_sorted.resize(_distributions.size());
      89             : 
      90             :   /* `markov_seed` is a member variable to store the seed input values for proposing
      91             :   the next set of Markov chain samples.*/
      92          22 :   _markov_seed.resize(_distributions.size());
      93             : 
      94          22 :   setNumberOfRandomSeeds(_num_random_seeds);
      95          22 : }
      96             : 
      97             : const unsigned int &
      98         288 : ParallelSubsetSimulation::getNumSamplesSub() const
      99             : {
     100         288 :   return _num_samplessub;
     101             : }
     102             : 
     103             : const bool &
     104          96 : ParallelSubsetSimulation::getUseAbsoluteValue() const
     105             : {
     106          96 :   return _use_absolute_value;
     107             : }
     108             : 
     109             : const Real &
     110         160 : ParallelSubsetSimulation::getSubsetProbability() const
     111             : {
     112         160 :   return _subset_probability;
     113             : }
     114             : 
     115             : void
     116         192 : ParallelSubsetSimulation::sampleSetUp(const SampleMode mode)
     117             : {
     118         192 :   if (_step <= 1 || _check_step == _step)
     119         112 :     return;
     120          80 :   _check_step = _step;
     121             : 
     122          80 :   if (_is_sampling_completed)
     123           0 :     mooseError("Internal bug: the adaptive sampling is supposed to be completed but another sample "
     124             :                "has been requested.");
     125             : 
     126          80 :   _subset = ((_step - 1) * getNumberOfRows()) / _num_samplessub;
     127          80 :   const unsigned int sub_ind = (_step - 1) - (_num_samplessub / getNumberOfRows()) * _subset;
     128          80 :   const unsigned int offset = sub_ind * getNumberOfRows();
     129             : 
     130             :   // Get and store the accepted samples input across all the procs from the previous step
     131         240 :   for (dof_id_type j = 0; j < _distributions.size(); ++j)
     132         480 :     for (dof_id_type ss = 0; ss < getNumberOfRows(); ++ss)
     133         320 :       _inputs_sto[j][ss + offset] = Normal::quantile(_distributions[j]->cdf(_inputs[j][ss]), 0, 1);
     134             : 
     135             :   // Get the accepted sample outputs across all the procs from the previous step
     136             :   std::vector<Real> tmp =
     137          80 :       _use_absolute_value ? AdaptiveMonteCarloUtils::computeVectorABS(_outputs) : _outputs;
     138          80 :   if (mode == Sampler::SampleMode::GLOBAL)
     139           0 :     _communicator.allgather(tmp);
     140             :   else
     141          80 :     _local_comm.allgather(tmp);
     142         240 :   for (dof_id_type ss = 0; ss < getNumberOfRows(); ++ss)
     143         160 :     _outputs_sto[ss + offset] = tmp[ss];
     144             : 
     145             :   // These are the subsequent subsets which use Markov Chain Monte Carlo sampling scheme
     146          80 :   if (_subset > 0)
     147             :   {
     148             :     // Check whether the subset index has changed
     149          64 :     if (sub_ind == 0)
     150             :     {
     151             :       // _inputs_sorted contains the input values corresponding to the largest po percentile
     152             :       // output values
     153          32 :       _inputs_sorted = AdaptiveMonteCarloUtils::sortInput(
     154          32 :           _inputs_sto, _outputs_sto, _num_samplessub, _subset_probability);
     155             :     }
     156             : 
     157             :     // Reinitialize the starting inputs values for the next set of Markov chains
     158          64 :     if (sub_ind % _count_max == 0)
     159             :     {
     160          32 :       const unsigned int soffset = (sub_ind / _count_max) * getNumberOfRows();
     161          96 :       for (dof_id_type j = 0; j < _distributions.size(); ++j)
     162          64 :         _markov_seed[j].assign(_inputs_sorted[j].begin() + soffset + getLocalRowBegin(),
     163          64 :                                _inputs_sorted[j].begin() + soffset + getLocalRowEnd());
     164             :     }
     165             :     // Otherwise, use the previously accepted input values to propose the next set of input
     166             :     // values
     167             :     else
     168             :     {
     169          96 :       for (dof_id_type j = 0; j < _distributions.size(); ++j)
     170          64 :         _markov_seed[j].assign(_inputs_sto[j].begin() + offset + getLocalRowBegin(),
     171          64 :                                _inputs_sto[j].begin() + offset + getLocalRowEnd());
     172             :     }
     173             :   }
     174             : 
     175             :   // check if we have completed the last sample (sub_ind == _num_samplessub /getNumberOfRows() - 1)
     176             :   // of the last subset (_subset == _num_subsets - 1)
     177          80 :   if (_subset == _num_subsets - 1 && sub_ind == _num_samplessub / getNumberOfRows() - 1)
     178          16 :     _is_sampling_completed = true;
     179             : }
     180             : 
     181             : Real
     182         480 : ParallelSubsetSimulation::computeSample(dof_id_type row_index, dof_id_type col_index)
     183             : {
     184         480 :   unsigned int seed_value = _step > 0 ? (_step - 1) * 2 : 0;
     185             :   Real val;
     186             : 
     187         480 :   if (_subset == 0)
     188         160 :     val = getRand(seed_value);
     189             :   else
     190             :   {
     191         320 :     const dof_id_type loc_ind = row_index - getLocalRowBegin();
     192         320 :     const Real rv = Normal::quantile(getRand(seed_value), _markov_seed[col_index][loc_ind], 1.0);
     193         320 :     const Real acceptance_ratio = std::log(Normal::pdf(rv, 0, 1)) -
     194         320 :                                   std::log(Normal::pdf(_markov_seed[col_index][loc_ind], 0, 1));
     195         320 :     const Real new_sample = acceptance_ratio > std::log(getRand(seed_value + 1))
     196         320 :                                 ? rv
     197         320 :                                 : _markov_seed[col_index][loc_ind];
     198         320 :     val = Normal::cdf(new_sample, 0, 1);
     199             :   }
     200             : 
     201         480 :   return _distributions[col_index]->quantile(val);
     202             : }