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 "IndependentMHDecision.h"
      11             : 
      12             : registerMooseObject("StochasticToolsApp", IndependentMHDecision);
      13             : 
      14             : InputParameters
      15          42 : IndependentMHDecision::validParams()
      16             : {
      17          42 :   InputParameters params = PMCMCDecision::validParams();
      18          42 :   params.addClassDescription("Perform decision making for independent Metropolis-Hastings MCMC.");
      19          84 :   params.addParam<ReporterValueName>(
      20             :       "seed_input", "seed_input", "The seed vector input for proposing new samples.");
      21          42 :   return params;
      22           0 : }
      23             : 
      24          20 : IndependentMHDecision::IndependentMHDecision(const InputParameters & parameters)
      25             :   : PMCMCDecision(parameters),
      26          20 :     _igmh(dynamic_cast<const IndependentGaussianMH *>(&_sampler)),
      27          40 :     _seed_input(declareValue<std::vector<Real>>("seed_input"))
      28             : {
      29             :   // Check whether the selected sampler is a M-H sampler or not
      30          20 :   if (!_igmh)
      31           0 :     paramError("sampler", "The selected sampler is not of type IndependentGaussianMH.");
      32             : 
      33          20 :   _seed_outputs.resize(_num_confg_values);
      34          20 :   _tpm_modified.assign(_props + 1, 1.0 / (_props + 1));
      35          20 : }
      36             : 
      37             : void
      38          60 : IndependentMHDecision::computeEvidence(std::vector<Real> & evidence,
      39             :                                        const DenseMatrix<Real> & input_matrix)
      40             : {
      41          60 :   std::vector<Real> out(_num_confg_values);
      42         360 :   for (unsigned int i = 0; i < evidence.size(); ++i)
      43             :   {
      44         300 :     evidence[i] = 0.0;
      45         900 :     for (unsigned int j = 0; j < _priors.size(); ++j)
      46         600 :       evidence[i] += (std::log(_priors[j]->pdf(input_matrix(i, j))) -
      47         600 :                       std::log(_priors[j]->pdf(_seed_input[j])));
      48         900 :     for (unsigned int j = 0; j < _num_confg_values; ++j)
      49         600 :       out[j] = (*_outputs_required)[j * _props + i];
      50         600 :     for (unsigned int j = 0; j < _likelihoods.size(); ++j)
      51         300 :       evidence[i] += (_likelihoods[j]->function(out) - _likelihoods[j]->function(_seed_outputs));
      52             :   }
      53          60 : }
      54             : 
      55             : void
      56          60 : IndependentMHDecision::computeTransitionVector(std::vector<Real> & tv,
      57             :                                                const std::vector<Real> & evidence)
      58             : {
      59         360 :   for (unsigned int i = 0; i < tv.size(); ++i)
      60         420 :     tv[i] = (1.0 / tv.size()) * std::exp(std::min(evidence[i], 0.0));
      61          60 :   _tpm_modified = tv;
      62          60 :   _tpm_modified.push_back((1.0 - std::accumulate(tv.begin(), tv.end(), 0.0)));
      63          60 :   _outputs_sto = (*_outputs_required);
      64          60 : }
      65             : 
      66             : void
      67         500 : IndependentMHDecision::nextSamples(std::vector<Real> & req_inputs,
      68             :                                    DenseMatrix<Real> & input_matrix,
      69             :                                    const std::vector<Real> & /*tv*/,
      70             :                                    const unsigned int & parallel_index)
      71             : {
      72         500 :   const bool value = (_tpm_modified[0] == 1.0 / (_props + 1));
      73         500 :   if (!value)
      74             :   {
      75             :     unsigned int index =
      76         300 :         AdaptiveMonteCarloUtils::weightedResample(_tpm_modified, _rnd_vec[parallel_index]);
      77         300 :     if (index < _props)
      78             :     {
      79         240 :       for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
      80         160 :         req_inputs[k] = input_matrix(index, k);
      81         240 :       for (unsigned int k = 0; k < _num_confg_values; ++k)
      82         160 :         (*_outputs_required)[k * _props + parallel_index] = _outputs_sto[k * _props + index];
      83             :     }
      84             :     else
      85             :     {
      86         220 :       req_inputs = _seed_input;
      87         660 :       for (unsigned int k = 0; k < _num_confg_values; ++k)
      88         440 :         (*_outputs_required)[k * _props + parallel_index] = _seed_outputs[k];
      89             :     }
      90             :   }
      91             :   else
      92             :   {
      93         600 :     for (unsigned int k = 0; k < _sampler.getNumberOfCols() - _num_confg_params; ++k)
      94         400 :       req_inputs[k] = input_matrix(parallel_index, k);
      95         200 :     _variance[parallel_index] = _new_var_samples[parallel_index];
      96             :   }
      97         500 : }
      98             : 
      99             : void
     100         100 : IndependentMHDecision::nextSeeds()
     101             : {
     102         100 :   _seed_input = _inputs[_props - 1];
     103         300 :   for (unsigned int k = 0; k < _num_confg_values; ++k)
     104         200 :     _seed_outputs[k] = (*_outputs_required)[(k + 1) * _props - 1];
     105         100 : }