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