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 "PMCMCBase.h"
      11             : #include "AdaptiveMonteCarloUtils.h"
      12             : #include "Uniform.h"
      13             : #include "DelimitedFileReader.h"
      14             : 
      15             : registerMooseObject("StochasticToolsApp", PMCMCBase);
      16             : 
      17             : InputParameters
      18         283 : PMCMCBase::validParams()
      19             : {
      20         283 :   InputParameters params = Sampler::validParams();
      21         283 :   params.addClassDescription("Parallel Markov chain Monte Carlo base.");
      22         566 :   params.addRequiredParam<std::vector<DistributionName>>(
      23             :       "prior_distributions", "The prior distributions of the parameters to be calibrated.");
      24         566 :   params.addParam<DistributionName>(
      25             :       "prior_variance", "The prior distribution of the variance parameter to be calibrated.");
      26         566 :   params.addRequiredParam<unsigned int>(
      27             :       "num_parallel_proposals",
      28             :       "Number of proposals to make and corresponding subApps executed in "
      29             :       "parallel.");
      30         566 :   params.addRequiredParam<FileName>("file_name", "Name of the CSV file with configuration values.");
      31         566 :   params.addParam<std::string>(
      32             :       "file_column_name", "Name of column in CSV file to use, by default first column is used.");
      33         566 :   params.addParam<unsigned int>(
      34             :       "num_columns", "Number of columns to be used in the CSV file with the configuration values.");
      35         566 :   params.addParam<std::vector<Real>>("lower_bound", "Lower bounds for making the next proposal.");
      36         566 :   params.addParam<std::vector<Real>>("upper_bound", "Upper bounds for making the next proposal.");
      37         566 :   params.addParam<Real>("variance_bound",
      38         566 :                         std::numeric_limits<Real>::max(),
      39             :                         "Upper bound for variance for making the next proposal.");
      40         566 :   params.addRequiredParam<std::vector<Real>>("initial_values",
      41             :                                              "The starting values of the inputs to be calibrated.");
      42         566 :   params.addParam<unsigned int>(
      43             :       "num_random_seeds",
      44         566 :       100000,
      45             :       "Initialize a certain number of random seeds. Change from the default only if you have to.");
      46         283 :   return params;
      47           0 : }
      48             : 
      49         159 : PMCMCBase::PMCMCBase(const InputParameters & parameters)
      50             :   : Sampler(parameters),
      51             :     TransientInterface(this),
      52         159 :     _num_parallel_proposals(getParam<unsigned int>("num_parallel_proposals")),
      53         374 :     _lower_bound(isParamValid("lower_bound") ? &getParam<std::vector<Real>>("lower_bound")
      54             :                                              : nullptr),
      55         370 :     _upper_bound(isParamValid("upper_bound") ? &getParam<std::vector<Real>>("upper_bound")
      56             :                                              : nullptr),
      57         318 :     _variance_bound(getParam<Real>("variance_bound")),
      58         318 :     _initial_values(getParam<std::vector<Real>>("initial_values")),
      59         318 :     _num_random_seeds(getParam<unsigned int>("num_random_seeds")),
      60         159 :     _seed_index(0),
      61         318 :     _rand_index(0)
      62             : {
      63             :   // Filling the `priors` vector with the user-provided distributions.
      64         159 :   for (const DistributionName & name :
      65         636 :        getParam<std::vector<DistributionName>>("prior_distributions"))
      66         318 :     _priors.push_back(&getDistributionByName(name));
      67             : 
      68             :   // Filling the `var_prior` object with the user-provided distribution for the variance.
      69         318 :   if (isParamValid("prior_variance"))
      70          48 :     _var_prior = &getDistributionByName(getParam<DistributionName>("prior_variance"));
      71             :   else
      72         135 :     _var_prior = nullptr;
      73             : 
      74             :   // Read the experimental configurations from a csv file
      75         318 :   MooseUtils::DelimitedFileReader reader(getParam<FileName>("file_name"));
      76         159 :   reader.read();
      77         159 :   _confg_values.resize(1);
      78         318 :   if (isParamValid("file_column_name"))
      79           0 :     _confg_values[0] = reader.getData(getParam<std::string>("file_column_name"));
      80         318 :   else if (isParamValid("num_columns"))
      81             :   {
      82           0 :     _confg_values.resize(getParam<unsigned int>("num_columns"));
      83           0 :     for (unsigned int i = 0; i < _confg_values.size(); ++i)
      84           0 :       _confg_values[i] = reader.getData(i);
      85             :   }
      86             :   else
      87         159 :     _confg_values[0] = reader.getData(0);
      88             : 
      89             :   // Setting the number of sampler rows to be equal to the number of parallel proposals
      90         159 :   setNumberOfRows(_num_parallel_proposals * _confg_values[0].size());
      91             : 
      92             :   // Setting the number of columns in the sampler matrix (equal to the number of distributions).
      93         159 :   setNumberOfCols(_priors.size() + _confg_values.size());
      94             : 
      95             :   // Resizing the vectors and vector of vectors
      96         159 :   _new_samples.resize(_num_parallel_proposals, _initial_values);
      97         159 :   _new_samples_confg.resize(_num_parallel_proposals * _confg_values[0].size(),
      98         318 :                             std::vector<Real>(_priors.size() + _confg_values.size(), 0.0));
      99         159 :   _rnd_vec.resize(_num_parallel_proposals);
     100         159 :   _new_var_samples.assign(_num_parallel_proposals, 0.0);
     101             : 
     102         159 :   setNumberOfRandomSeeds(_num_random_seeds);
     103         159 :   setAutoAdvanceGenerators(false);
     104             : 
     105             :   // Check whether both the lower and the upper bounds are specified and of same size
     106         159 :   bool bound_check1 = _lower_bound && !_upper_bound;
     107         159 :   bool bound_check2 = !_lower_bound && _upper_bound;
     108         159 :   if (bound_check1 || bound_check2)
     109           2 :     mooseError("Both lower and upper bounds should be specified.");
     110         157 :   bool size_check = _lower_bound ? ((*_lower_bound).size() != (*_upper_bound).size()) : 0;
     111             :   if (size_check)
     112           2 :     mooseError("Lower and upper bounds should be of the same size.");
     113             : 
     114             :   // Check whether the priors, bounds, and initial values are all of the same size
     115         155 :   if (_priors.size() != _initial_values.size())
     116           2 :     mooseError("The priors and initial values should be of the same size.");
     117         153 : }
     118             : 
     119             : void
     120         288 : PMCMCBase::proposeSamples()
     121             : {
     122         864 :   for (unsigned int j = 0; j < _num_parallel_proposals; ++j)
     123        1728 :     for (unsigned int i = 0; i < _priors.size(); ++i)
     124        1152 :       _new_samples[j][i] = _priors[i]->quantile(random());
     125         288 : }
     126             : 
     127             : void
     128         894 : PMCMCBase::executeSetUp()
     129             : {
     130         894 :   _seed_index = _t_step;
     131         894 :   _rand_index = 0;
     132             : 
     133             :   // Filling the new_samples vector of vectors with new proposal samples
     134         894 :   proposeSamples();
     135             : 
     136             :   // Draw random numbers to facilitate decision making later on
     137        4500 :   for (unsigned int j = 0; j < _num_parallel_proposals; ++j)
     138        3606 :     _rnd_vec[j] = random();
     139         894 : }
     140             : 
     141             : Real
     142       18536 : PMCMCBase::random()
     143             : {
     144       18536 :   return getRand(_rand_index++, _seed_index);
     145             : }
     146             : 
     147             : unsigned int
     148        4320 : PMCMCBase::randomIndex(const unsigned int & upper_bound, const unsigned int & exclude)
     149             : {
     150        4320 :   auto req_index = exclude;
     151        9936 :   while (req_index == exclude)
     152        5616 :     req_index = getRandl(_rand_index++, 0, upper_bound, _seed_index);
     153        4320 :   return req_index;
     154             : }
     155             : 
     156             : std::pair<unsigned int, unsigned int>
     157        1536 : PMCMCBase::randomIndexPair(const unsigned int & upper_bound, const unsigned int & exclude)
     158             : {
     159        1536 :   auto req_index1 = randomIndex(upper_bound, exclude);
     160             :   auto req_index2 = req_index1;
     161        3600 :   while (req_index1 == req_index2)
     162        2064 :     req_index2 = randomIndex(upper_bound, exclude);
     163        1536 :   return {req_index1, req_index2};
     164             : }
     165             : 
     166             : void
     167        8400 : PMCMCBase::combineWithExperimentalConfig()
     168             : {
     169             :   unsigned int index1;
     170             :   int index2 = -1;
     171             :   std::vector<Real> tmp;
     172       85200 :   for (unsigned int i = 0; i < _num_parallel_proposals * _confg_values[0].size(); ++i)
     173             :   {
     174       76800 :     index1 = i % _num_parallel_proposals;
     175       76800 :     if (index1 == 0)
     176       16800 :       ++index2;
     177       76800 :     tmp = _new_samples[index1];
     178      153600 :     for (unsigned int j = 0; j < _confg_values.size(); ++j)
     179       76800 :       tmp.push_back(_confg_values[j][index2]);
     180       76800 :     _new_samples_confg[i] = tmp;
     181             :   }
     182        8400 : }
     183             : 
     184             : const std::vector<Real> &
     185         120 : PMCMCBase::getRandomNumbers() const
     186             : {
     187         120 :   return _rnd_vec;
     188             : }
     189             : 
     190             : const std::vector<Real> &
     191         124 : PMCMCBase::getVarSamples() const
     192             : {
     193         124 :   return _new_var_samples;
     194             : }
     195             : 
     196             : const std::vector<std::vector<Real>> &
     197           0 : PMCMCBase::getSamples() const
     198             : {
     199           0 :   return _new_samples;
     200             : }
     201             : 
     202             : const std::vector<const Distribution *>
     203         120 : PMCMCBase::getPriors() const
     204             : {
     205         120 :   return _priors;
     206             : }
     207             : 
     208             : const Distribution *
     209         124 : PMCMCBase::getVarPrior() const
     210             : {
     211         124 :   return _var_prior;
     212             : }
     213             : 
     214             : Real
     215        8400 : PMCMCBase::computeSample(dof_id_type row_index, dof_id_type col_index)
     216             : {
     217        8400 :   if (_t_step < 1)
     218        4680 :     for (unsigned int i = 0; i < _num_parallel_proposals; ++i)
     219        3840 :       _new_samples[i] = _initial_values;
     220             : 
     221             :   // Combine the proposed samples with experimental configurations
     222        8400 :   combineWithExperimentalConfig();
     223             : 
     224        8400 :   return _new_samples_confg[row_index][col_index];
     225             : }