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          96 : AffineInvariantStretchSampler::validParams()
      16             : {
      17          96 :   InputParameters params = PMCMCBase::validParams();
      18          96 :   params.addClassDescription("Perform Affine Invariant Ensemble MCMC with stretch sampler.");
      19         192 :   params.addRequiredParam<ReporterName>(
      20             :       "previous_state", "Reporter value with the previous state of all the walkers.");
      21         192 :   params.addRequiredParam<ReporterName>(
      22             :       "previous_state_var",
      23             :       "Reporter value with the previous state of all the walkers for variance.");
      24         192 :   params.addParam<Real>("step_size", 2.0, "Step size for each of the walkers.");
      25          96 :   return params;
      26           0 : }
      27             : 
      28          56 : AffineInvariantStretchSampler::AffineInvariantStretchSampler(const InputParameters & parameters)
      29             :   : PMCMCBase(parameters),
      30          56 :     _step_size(getParam<Real>("step_size")),
      31          56 :     _previous_state(getReporterValue<std::vector<std::vector<Real>>>("previous_state")),
      32         112 :     _previous_state_var(getReporterValue<std::vector<Real>>("previous_state_var"))
      33             : {
      34          56 :   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          56 :   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          56 :   _affine_step.resize(_num_parallel_proposals);
      45          56 : }
      46             : 
      47             : void
      48         200 : AffineInvariantStretchSampler::proposeSamples(const unsigned int seed_value)
      49             : {
      50         200 :   unsigned int j = 0;
      51             :   bool indicator;
      52         200 :   unsigned int index_req = 0;
      53        1400 :   while (j < _num_parallel_proposals)
      54             :   {
      55             :     indicator = 0;
      56        1000 :     _affine_step[j] = Utility::pow<2>((_step_size - 1.0) * getRand(seed_value) + 1.0) / _step_size;
      57        1000 :     randomIndex(_num_parallel_proposals, j, seed_value, index_req);
      58        3000 :     for (unsigned int i = 0; i < _priors.size(); ++i)
      59             :     {
      60        2000 :       _new_samples[j][i] =
      61        2000 :           (_t_step > decisionStep())
      62        2000 :               ? (_previous_state[index_req][i] +
      63        1200 :                  _affine_step[j] * (_previous_state[j][i] - _previous_state[index_req][i]))
      64         800 :               : _priors[i]->quantile(getRand(seed_value));
      65        2000 :       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        1000 :     if (_var_prior)
      72             :     {
      73           0 :       _new_var_samples[j] =
      74           0 :           (_t_step > 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(getRand(seed_value));
      78           0 :       if (_new_var_samples[j] < 0.0)
      79             :         indicator = 1;
      80             :     }
      81        1000 :     if (!indicator)
      82        1000 :       ++j;
      83             :   }
      84         200 : }
      85             : 
      86             : const std::vector<Real> &
      87          40 : AffineInvariantStretchSampler::getAffineStepSize() const
      88             : {
      89          40 :   return _affine_step;
      90             : }