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 "IndependentGaussianMH.h"
      11             : #include "Normal.h"
      12             : #include "TruncatedNormal.h"
      13             : 
      14             : registerMooseObject("StochasticToolsApp", IndependentGaussianMH);
      15             : 
      16             : InputParameters
      17         100 : IndependentGaussianMH::validParams()
      18             : {
      19         100 :   InputParameters params = PMCMCBase::validParams();
      20         100 :   params.addClassDescription("Perform M-H MCMC sampling with independent Gaussian propoposals.");
      21         200 :   params.addRequiredParam<ReporterName>("seed_inputs",
      22             :                                         "Reporter with seed inputs values for the next proposals.");
      23         200 :   params.addRequiredParam<std::vector<Real>>("std_prop",
      24             :                                              "Standard deviations for making the next proposal.");
      25         100 :   return params;
      26           0 : }
      27             : 
      28          60 : IndependentGaussianMH::IndependentGaussianMH(const InputParameters & parameters)
      29             :   : PMCMCBase(parameters),
      30          60 :     _seed_inputs(getReporterValue<std::vector<Real>>("seed_inputs")),
      31         180 :     _std_prop(getParam<std::vector<Real>>("std_prop"))
      32             : {
      33             :   // Error check for sizes of proposal stds
      34          60 :   if (_std_prop.size() != _priors.size())
      35           4 :     paramError("std_prop",
      36             :                "The number of proposal stds, initial values, and priors should be the same.");
      37          56 : }
      38             : 
      39             : void
      40         200 : IndependentGaussianMH::proposeSamples(const unsigned int seed_value)
      41             : {
      42         200 :   std::vector<Real> old_sample = (_t_step > decisionStep()) ? _seed_inputs : _initial_values;
      43        1200 :   for (unsigned int j = 0; j < _num_parallel_proposals; ++j)
      44        3000 :     for (unsigned int i = 0; i < _priors.size(); ++i)
      45             :     {
      46        2000 :       if (_lower_bound)
      47           0 :         _new_samples[j][i] = TruncatedNormal::quantile(getRand(seed_value),
      48             :                                                        old_sample[i],
      49           0 :                                                        _std_prop[i],
      50             :                                                        (*_lower_bound)[i],
      51           0 :                                                        (*_upper_bound)[i]);
      52             :       else
      53        2000 :         _new_samples[j][i] = Normal::quantile(getRand(seed_value), old_sample[i], _std_prop[i]);
      54             :     }
      55         200 : }