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         500 : PMCMCDecision::validParams()
      18             : {
      19         500 :   InputParameters params = GeneralReporter::validParams();
      20         500 :   params += LikelihoodInterface::validParams();
      21         500 :   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        1000 :   params.addRequiredParam<ReporterName>("output_value",
      24             :                                         "Value of the model output from the SubApp.");
      25        1000 :   params.addParam<ReporterValueName>(
      26             :       "outputs_required",
      27             :       "outputs_required",
      28             :       "Modified value of the model output from this reporter class.");
      29        1000 :   params.addParam<ReporterValueName>("inputs", "inputs", "Uncertain inputs to the model.");
      30        1000 :   params.addParam<ReporterValueName>("tpm", "tpm", "The transition probability matrix.");
      31        1000 :   params.addParam<ReporterValueName>("variance", "variance", "Model variance term.");
      32        1000 :   params.addParam<ReporterValueName>(
      33             :       "noise", "noise", "Model noise term to pass to Likelihoods object.");
      34        1000 :   params.addRequiredParam<SamplerName>("sampler", "The sampler object.");
      35        1000 :   params.addRequiredParam<std::vector<UserObjectName>>("likelihoods", "Names of likelihoods.");
      36         500 :   return params;
      37           0 : }
      38             : 
      39         240 : PMCMCDecision::PMCMCDecision(const InputParameters & parameters)
      40             :   : GeneralReporter(parameters),
      41             :     LikelihoodInterface(parameters),
      42         240 :     _output_value(getReporterValue<std::vector<Real>>("output_value", REPORTER_MODE_DISTRIBUTED)),
      43         240 :     _outputs_required(declareValue<std::vector<Real>>("outputs_required")),
      44         240 :     _inputs(declareValue<std::vector<std::vector<Real>>>("inputs")),
      45         240 :     _tpm(declareValue<std::vector<Real>>("tpm")),
      46         240 :     _variance(declareValue<std::vector<Real>>("variance")),
      47         240 :     _noise(declareValue<Real>("noise")),
      48         240 :     _sampler(getSampler("sampler")),
      49         240 :     _pmcmc(dynamic_cast<const PMCMCBase *>(&_sampler)),
      50         240 :     _rnd_vec(_pmcmc->getRandomNumbers()),
      51         240 :     _new_var_samples(_pmcmc->getVarSamples()),
      52         240 :     _priors(_pmcmc->getPriors()),
      53         240 :     _var_prior(_pmcmc->getVarPrior()),
      54         240 :     _local_comm(_sampler.getLocalComm()),
      55         720 :     _check_step(std::numeric_limits<int>::max())
      56             : {
      57             :   // Filling the `likelihoods` vector with the user-provided distributions.
      58         720 :   for (const UserObjectName & name : getParam<std::vector<UserObjectName>>("likelihoods"))
      59         240 :     _likelihoods.push_back(getLikelihoodFunctionByName(name));
      60             : 
      61             :   // Check whether the selected sampler is an MCMC sampler or not
      62         240 :   if (!_pmcmc)
      63           0 :     paramError("sampler", "The selected sampler is not of type MCMC.");
      64             : 
      65             :   // Fetching the sampler characteristics
      66         240 :   _props = _pmcmc->getNumParallelProposals();
      67         240 :   _num_confg_values = _pmcmc->getNumberOfConfigValues();
      68         240 :   _num_confg_params = _pmcmc->getNumberOfConfigParams();
      69             : 
      70             :   // Resizing the data arrays to transmit to the output file
      71         240 :   _inputs.resize(_props);
      72        1200 :   for (unsigned int i = 0; i < _props; ++i)
      73         960 :     _inputs[i].resize(_sampler.getNumberOfCols() - _num_confg_params);
      74         240 :   _outputs_required.resize(_sampler.getNumberOfRows());
      75         240 :   _tpm.resize(_props);
      76         240 :   _variance.resize(_props);
      77         240 : }
      78             : 
      79             : void
      80         680 : PMCMCDecision::computeEvidence(std::vector<Real> & evidence, const DenseMatrix<Real> & input_matrix)
      81             : {
      82         680 :   std::vector<Real> out1(_num_confg_values);
      83         680 :   std::vector<Real> out2(_num_confg_values);
      84        3120 :   for (unsigned int i = 0; i < evidence.size(); ++i)
      85             :   {
      86        2440 :     evidence[i] = 0.0;
      87        7320 :     for (unsigned int j = 0; j < _priors.size(); ++j)
      88        4880 :       evidence[i] += (std::log(_priors[j]->pdf(input_matrix(i, j))) -
      89        4880 :                       std::log(_priors[j]->pdf(_data_prev(i, j))));
      90        7320 :     for (unsigned int j = 0; j < _num_confg_values; ++j)
      91             :     {
      92        4880 :       out1[j] = _outputs_required[j * _props + i];
      93        4880 :       out2[j] = _outputs_prev[j * _props + i];
      94             :     }
      95        2440 :     if (_var_prior)
      96             :     {
      97         600 :       evidence[i] += (std::log(_var_prior->pdf(_new_var_samples[i])) -
      98         600 :                       std::log(_var_prior->pdf(_var_prev[i])));
      99         600 :       _noise = std::sqrt(_new_var_samples[i]);
     100        1200 :       for (unsigned int j = 0; j < _likelihoods.size(); ++j)
     101         600 :         evidence[i] += _likelihoods[j]->function(out1);
     102         600 :       _noise = std::sqrt(_var_prev[i]);
     103        1200 :       for (unsigned int j = 0; j < _likelihoods.size(); ++j)
     104         600 :         evidence[i] -= _likelihoods[j]->function(out2);
     105             :     }
     106             :     else
     107        3680 :       for (unsigned int j = 0; j < _likelihoods.size(); ++j)
     108        1840 :         evidence[i] += (_likelihoods[j]->function(out1) - _likelihoods[j]->function(out2));
     109             :   }
     110         680 : }
     111             : 
     112             : void
     113         320 : PMCMCDecision::computeTransitionVector(std::vector<Real> & tv,
     114             :                                        const std::vector<Real> & /*evidence*/)
     115             : {
     116         320 :   tv.assign(_props, 1.0);
     117         320 : }
     118             : 
     119             : void
     120        3800 : 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        3800 :   if (tv[parallel_index] >= _rnd_vec[parallel_index])
     126             :   {
     127        8040 :     for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
     128        5360 :       req_inputs[k] = input_matrix(parallel_index, k);
     129        2680 :     _variance[parallel_index] = _new_var_samples[parallel_index];
     130             :   }
     131             :   else
     132             :   {
     133        3360 :     for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
     134             :     {
     135        2240 :       req_inputs[k] = _data_prev(parallel_index, k);
     136        2240 :       input_matrix(parallel_index, k) = _data_prev(parallel_index, k);
     137             :     }
     138        1120 :     if (_var_prior)
     139         440 :       _variance[parallel_index] = _var_prev[parallel_index];
     140        3360 :     for (unsigned int k = 0; k < _num_confg_values; ++k)
     141        2240 :       _outputs_required[k * _props + parallel_index] = _outputs_prev[k * _props + parallel_index];
     142             :   }
     143        3800 : }
     144             : 
     145             : void
     146        1200 : PMCMCDecision::execute()
     147             : {
     148        1200 :   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        1200 :   DenseMatrix<Real> data_in(_sampler.getNumberOfRows(), _sampler.getNumberOfCols());
     156        3600 :   for (dof_id_type ss = _sampler.getLocalRowBegin(); ss < _sampler.getLocalRowEnd(); ++ss)
     157             :   {
     158        2400 :     const auto data = _sampler.getNextLocalRow();
     159        9600 :     for (unsigned int j = 0; j < _sampler.getNumberOfCols(); ++j)
     160        7200 :       data_in(ss, j) = data[j];
     161             :   }
     162        1200 :   _local_comm.sum(data_in.get_values());
     163        1200 :   _outputs_required = _output_value;
     164        1200 :   _local_comm.allgather(_outputs_required);
     165             : 
     166             :   // Compute the evidence and transition vectors
     167        1200 :   std::vector<Real> evidence(_props);
     168        1200 :   if (_t_step > _pmcmc->decisionStep())
     169             :   {
     170         800 :     computeEvidence(evidence, data_in);
     171         800 :     computeTransitionVector(_tpm, evidence);
     172             :   }
     173             :   else
     174         400 :     _tpm.assign(_props, 1.0);
     175             : 
     176             :   // Accept/reject the proposed samples and assign the correct outputs
     177        1200 :   std::vector<Real> req_inputs(_sampler.getNumberOfCols() - _num_confg_params);
     178        6000 :   for (unsigned int i = 0; i < _props; ++i)
     179             :   {
     180        4800 :     nextSamples(req_inputs, data_in, _tpm, i);
     181        4800 :     _inputs[i] = req_inputs;
     182             :   }
     183             : 
     184             :   // Compute the next seeds to facilitate proposals (not always required)
     185        1200 :   nextSeeds();
     186             : 
     187             :   // Store data from previous step
     188        1200 :   _data_prev = data_in;
     189        1200 :   _outputs_prev = _outputs_required;
     190        1200 :   _var_prev = _variance;
     191             : 
     192             :   // Track the current step
     193        1200 :   _check_step = _t_step;
     194        1200 : }