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