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