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 "PMCMCDecision.h"
      11             : #include "Sampler.h"
      12             : #include "DenseMatrix.h"
      13             : 
      14             : registerMooseObject("StochasticToolsApp", PMCMCDecision);
      15             : 
      16             : InputParameters
      17         548 : PMCMCDecision::validParams()
      18             : {
      19         548 :   InputParameters params = GeneralReporter::validParams();
      20         548 :   params += LikelihoodInterface::validParams();
      21         548 :   params.addClassDescription("Generic reporter which decides whether or not to accept a proposed "
      22             :                              "sample in parallel Markov chain Monte Carlo type of algorithms.");
      23        1096 :   params.addRequiredParam<ReporterName>("output_value",
      24             :                                         "Value of the model output from the SubApp.");
      25        1096 :   params.addParam<ReporterValueName>(
      26             :       "outputs_required",
      27             :       "outputs_required",
      28             :       "Modified value of the model output from this reporter class.");
      29        1096 :   params.addParam<ReporterValueName>("inputs", "inputs", "Uncertain inputs to the model.");
      30        1096 :   params.addParam<ReporterValueName>("tpm", "tpm", "The transition probability matrix.");
      31        1096 :   params.addParam<ReporterValueName>("variance", "variance", "Model variance term.");
      32        1096 :   params.addParam<ReporterValueName>(
      33             :       "noise", "noise", "Model noise term to pass to Likelihoods object.");
      34        1096 :   params.addRequiredParam<SamplerName>("sampler", "The sampler object.");
      35        1096 :   params.addRequiredParam<std::vector<UserObjectName>>("likelihoods", "Names of likelihoods.");
      36         548 :   return params;
      37           0 : }
      38             : 
      39         264 : PMCMCDecision::PMCMCDecision(const InputParameters & parameters)
      40             :   : GeneralReporter(parameters),
      41             :     LikelihoodInterface(parameters),
      42         264 :     _output_value(getReporterValue<std::vector<Real>>("output_value", REPORTER_MODE_DISTRIBUTED)),
      43         264 :     _outputs_required(declareValue<std::vector<Real>>("outputs_required")),
      44         264 :     _inputs(declareValue<std::vector<std::vector<Real>>>("inputs")),
      45         264 :     _tpm(declareValue<std::vector<Real>>("tpm")),
      46         264 :     _variance(declareValue<std::vector<Real>>("variance")),
      47         264 :     _noise(declareValue<Real>("noise")),
      48         264 :     _sampler(getSampler("sampler")),
      49         264 :     _pmcmc(dynamic_cast<const PMCMCBase *>(&_sampler)),
      50         264 :     _rnd_vec(_pmcmc->getRandomNumbers()),
      51         264 :     _new_var_samples(_pmcmc->getVarSamples()),
      52         264 :     _priors(_pmcmc->getPriors()),
      53         264 :     _var_prior(_pmcmc->getVarPrior()),
      54         264 :     _local_comm(_sampler.getLocalComm()),
      55         528 :     _check_step(std::numeric_limits<int>::max())
      56             : {
      57             :   // Filling the `likelihoods` vector with the user-provided distributions.
      58         792 :   for (const UserObjectName & name : getParam<std::vector<UserObjectName>>("likelihoods"))
      59         264 :     _likelihoods.push_back(getLikelihoodFunctionByName(name));
      60             : 
      61             :   // Check whether the selected sampler is an MCMC sampler or not
      62         264 :   if (!_pmcmc)
      63           0 :     paramError("sampler", "The selected sampler is not of type MCMC.");
      64             : 
      65             :   // Fetching the sampler characteristics
      66         264 :   _props = _pmcmc->getNumParallelProposals();
      67         264 :   _num_confg_values = _pmcmc->getNumberOfConfigValues();
      68         264 :   _num_confg_params = _pmcmc->getNumberOfConfigParams();
      69             : 
      70             :   // Resizing the data arrays to transmit to the output file
      71         264 :   _inputs.resize(_props);
      72        1320 :   for (unsigned int i = 0; i < _props; ++i)
      73        1056 :     _inputs[i].resize(_sampler.getNumberOfCols() - _num_confg_params);
      74         264 :   _outputs_required.resize(_sampler.getNumberOfRows());
      75         264 :   _tpm.resize(_props);
      76         264 :   _variance.resize(_props);
      77         264 : }
      78             : 
      79             : void
      80         748 : PMCMCDecision::computeEvidence(std::vector<Real> & evidence, const DenseMatrix<Real> & input_matrix)
      81             : {
      82         748 :   std::vector<Real> out1(_num_confg_values);
      83         748 :   std::vector<Real> out2(_num_confg_values);
      84        3432 :   for (unsigned int i = 0; i < evidence.size(); ++i)
      85             :   {
      86        2684 :     evidence[i] = 0.0;
      87        8052 :     for (unsigned int j = 0; j < _priors.size(); ++j)
      88        5368 :       evidence[i] += (std::log(_priors[j]->pdf(input_matrix(i, j))) -
      89        5368 :                       std::log(_priors[j]->pdf(_data_prev(i, j))));
      90        8052 :     for (unsigned int j = 0; j < _num_confg_values; ++j)
      91             :     {
      92        5368 :       out1[j] = _outputs_required[j * _props + i];
      93        5368 :       out2[j] = _outputs_prev[j * _props + i];
      94             :     }
      95        2684 :     if (_var_prior)
      96             :     {
      97         660 :       evidence[i] += (std::log(_var_prior->pdf(_new_var_samples[i])) -
      98         660 :                       std::log(_var_prior->pdf(_var_prev[i])));
      99         660 :       _noise = std::sqrt(_new_var_samples[i]);
     100        1320 :       for (unsigned int j = 0; j < _likelihoods.size(); ++j)
     101         660 :         evidence[i] += _likelihoods[j]->function(out1);
     102         660 :       _noise = std::sqrt(_var_prev[i]);
     103        1320 :       for (unsigned int j = 0; j < _likelihoods.size(); ++j)
     104         660 :         evidence[i] -= _likelihoods[j]->function(out2);
     105             :     }
     106             :     else
     107        4048 :       for (unsigned int j = 0; j < _likelihoods.size(); ++j)
     108        2024 :         evidence[i] += (_likelihoods[j]->function(out1) - _likelihoods[j]->function(out2));
     109             :   }
     110         748 : }
     111             : 
     112             : void
     113         352 : PMCMCDecision::computeTransitionVector(std::vector<Real> & tv,
     114             :                                        const std::vector<Real> & /*evidence*/)
     115             : {
     116         352 :   tv.assign(_props, 1.0);
     117         352 : }
     118             : 
     119             : void
     120        4180 : PMCMCDecision::nextSamples(std::vector<Real> & req_inputs,
     121             :                            DenseMatrix<Real> & input_matrix,
     122             :                            const std::vector<Real> & tv,
     123             :                            const unsigned int & parallel_index)
     124             : {
     125        4180 :   if (tv[parallel_index] >= _rnd_vec[parallel_index])
     126             :   {
     127        8844 :     for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
     128        5896 :       req_inputs[k] = input_matrix(parallel_index, k);
     129        2948 :     _variance[parallel_index] = _new_var_samples[parallel_index];
     130             :   }
     131             :   else
     132             :   {
     133        3696 :     for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
     134             :     {
     135        2464 :       req_inputs[k] = _data_prev(parallel_index, k);
     136        2464 :       input_matrix(parallel_index, k) = _data_prev(parallel_index, k);
     137             :     }
     138        1232 :     if (_var_prior)
     139         484 :       _variance[parallel_index] = _var_prev[parallel_index];
     140        3696 :     for (unsigned int k = 0; k < _num_confg_values; ++k)
     141        2464 :       _outputs_required[k * _props + parallel_index] = _outputs_prev[k * _props + parallel_index];
     142             :   }
     143        4180 : }
     144             : 
     145             : void
     146        1320 : PMCMCDecision::execute()
     147             : {
     148        1320 :   if (_sampler.getNumberOfLocalRows() == 0 || _check_step == _t_step)
     149             :   {
     150           0 :     _check_step = _t_step;
     151           0 :     return;
     152             :   }
     153             : 
     154             :   // Gather inputs and outputs from the sampler and subApps
     155        1320 :   DenseMatrix<Real> data_in(_sampler.getNumberOfRows(), _sampler.getNumberOfCols());
     156        3960 :   for (dof_id_type ss = _sampler.getLocalRowBegin(); ss < _sampler.getLocalRowEnd(); ++ss)
     157             :   {
     158        2640 :     const auto data = _sampler.getNextLocalRow();
     159       10560 :     for (unsigned int j = 0; j < _sampler.getNumberOfCols(); ++j)
     160        7920 :       data_in(ss, j) = data[j];
     161        2640 :   }
     162        1320 :   _local_comm.sum(data_in.get_values());
     163        1320 :   _outputs_required = _output_value;
     164        1320 :   _local_comm.allgather(_outputs_required);
     165             : 
     166             :   // Compute the evidence and transition vectors
     167        1320 :   std::vector<Real> evidence(_props);
     168        1320 :   if (_t_step > _pmcmc->decisionStep())
     169             :   {
     170         880 :     computeEvidence(evidence, data_in);
     171         880 :     computeTransitionVector(_tpm, evidence);
     172             :   }
     173             :   else
     174         440 :     _tpm.assign(_props, 1.0);
     175             : 
     176             :   // Accept/reject the proposed samples and assign the correct outputs
     177        1320 :   std::vector<Real> req_inputs(_sampler.getNumberOfCols() - _num_confg_params);
     178        6600 :   for (unsigned int i = 0; i < _props; ++i)
     179             :   {
     180        5280 :     nextSamples(req_inputs, data_in, _tpm, i);
     181        5280 :     _inputs[i] = req_inputs;
     182             :   }
     183             : 
     184             :   // Compute the next seeds to facilitate proposals (not always required)
     185        1320 :   nextSeeds();
     186             : 
     187             :   // Store data from previous step
     188        1320 :   _data_prev = data_in;
     189        1320 :   _outputs_prev = _outputs_required;
     190        1320 :   _var_prev = _variance;
     191             : 
     192             :   // Track the current step
     193        1320 :   _check_step = _t_step;
     194        1320 : }