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 "AffineInvariantStretchSampler.h"
      11             : 
      12             : registerMooseObject("StochasticToolsApp", AffineInvariantStretchSampler);
      13             : 
      14             : InputParameters
      15          44 : AffineInvariantStretchSampler::validParams()
      16             : {
      17          44 :   InputParameters params = PMCMCBase::validParams();
      18          44 :   params.addClassDescription("Perform Affine Invariant Ensemble MCMC with stretch sampler.");
      19          88 :   params.addRequiredParam<ReporterName>(
      20             :       "previous_state", "Reporter value with the previous state of all the walkers.");
      21          88 :   params.addRequiredParam<ReporterName>(
      22             :       "previous_state_var",
      23             :       "Reporter value with the previous state of all the walkers for variance.");
      24          88 :   params.addParam<Real>("step_size", 2.0, "Step size for each of the walkers.");
      25          44 :   return params;
      26           0 : }
      27             : 
      28          24 : AffineInvariantStretchSampler::AffineInvariantStretchSampler(const InputParameters & parameters)
      29             :   : PMCMCBase(parameters),
      30          24 :     _step_size(getParam<Real>("step_size")),
      31          24 :     _previous_state(getReporterValue<std::vector<std::vector<Real>>>("previous_state")),
      32          48 :     _previous_state_var(getReporterValue<std::vector<Real>>("previous_state_var"))
      33             : {
      34          24 :   if (_num_parallel_proposals < 3)
      35           0 :     paramError("num_parallel_proposals",
      36             :                "At least three parallel proposals should be used for the Stretch Sampler.");
      37             : 
      38          24 :   if (_num_parallel_proposals < _priors.size())
      39           0 :     mooseWarning(
      40             :         "It is recommended that the parallel proposals be greater than or equal to the "
      41             :         "inferred parameters. This will allow the sampler to not get stuck on a hyper-plane.");
      42             : 
      43             :   // Assign the correct size to the step size vector
      44          24 :   _affine_step.resize(_num_parallel_proposals);
      45          24 : }
      46             : 
      47             : void
      48         144 : AffineInvariantStretchSampler::proposeSamples()
      49             : {
      50         144 :   unsigned int j = 0;
      51             :   bool indicator;
      52             :   unsigned int index_req = 0;
      53        1008 :   while (j < _num_parallel_proposals)
      54             :   {
      55             :     indicator = 0;
      56         720 :     _affine_step[j] = Utility::pow<2>((_step_size - 1.0) * random() + 1.0) / _step_size;
      57         720 :     index_req = randomIndex(_num_parallel_proposals, j);
      58        2160 :     for (unsigned int i = 0; i < _priors.size(); ++i)
      59             :     {
      60        1440 :       _new_samples[j][i] =
      61        1440 :           (_t_step + 1 > decisionStep())
      62        1440 :               ? (_previous_state[index_req][i] +
      63         960 :                  _affine_step[j] * (_previous_state[j][i] - _previous_state[index_req][i]))
      64         480 :               : _priors[i]->quantile(random());
      65        1440 :       if (_lower_bound)
      66           0 :         indicator =
      67           0 :             (_new_samples[j][i] < (*_lower_bound)[i] || _new_samples[j][i] > (*_upper_bound)[i])
      68           0 :                 ? 1
      69             :                 : indicator;
      70             :     }
      71         720 :     if (_var_prior)
      72             :     {
      73           0 :       _new_var_samples[j] =
      74           0 :           (_t_step + 1 > decisionStep())
      75           0 :               ? (_previous_state_var[index_req] +
      76           0 :                  _affine_step[j] * (_previous_state_var[j] - _previous_state_var[index_req]))
      77           0 :               : _var_prior->quantile(random());
      78           0 :       if (_new_var_samples[j] < 0.0 || _new_var_samples[j] > _variance_bound)
      79             :         indicator = 1;
      80             :     }
      81         720 :     if (!indicator)
      82         720 :       ++j;
      83             :   }
      84         144 : }
      85             : 
      86             : const std::vector<Real> &
      87          20 : AffineInvariantStretchSampler::getAffineStepSize() const
      88             : {
      89          20 :   return _affine_step;
      90             : }