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 "AffineInvariantDES.h"
      11             : #include "Normal.h"
      12             : #include "Uniform.h"
      13             : 
      14             : registerMooseObject("StochasticToolsApp", AffineInvariantDES);
      15             : 
      16             : /*
      17             :  Tuning options for the internal parameters
      18             :   1. Braak2006_static:
      19             :   - the gamma param is set to 2.38 / sqrt(2 * dim)
      20             :   - the b param is set to (scale) * 1e-6
      21             : */
      22             : 
      23             : InputParameters
      24          90 : AffineInvariantDES::validParams()
      25             : {
      26          90 :   InputParameters params = PMCMCBase::validParams();
      27          90 :   params.addClassDescription("Perform Affine Invariant Ensemble MCMC with differential sampler.");
      28         180 :   params.addRequiredParam<ReporterName>(
      29             :       "previous_state", "Reporter value with the previous state of all the walkers.");
      30         180 :   params.addRequiredParam<ReporterName>(
      31             :       "previous_state_var",
      32             :       "Reporter value with the previous state of all the walkers for variance.");
      33         180 :   MooseEnum tuning_option("Braak2006_static", "Braak2006_static");
      34         180 :   params.addParam<MooseEnum>(
      35             :       "tuning_option", tuning_option, "The tuning option for internal parameters.");
      36         180 :   params.addParam<std::vector<Real>>("scales", "Scales for the parameters.");
      37          90 :   return params;
      38          90 : }
      39             : 
      40          50 : AffineInvariantDES::AffineInvariantDES(const InputParameters & parameters)
      41             :   : PMCMCBase(parameters),
      42          50 :     _previous_state(getReporterValue<std::vector<std::vector<Real>>>("previous_state")),
      43          50 :     _previous_state_var(getReporterValue<std::vector<Real>>("previous_state_var")),
      44         150 :     _tuning_option(getParam<MooseEnum>("tuning_option"))
      45             : {
      46          50 :   if (_num_parallel_proposals < 5)
      47           0 :     paramError(
      48             :         "num_parallel_proposals",
      49             :         "At least five parallel proposals should be used for the Differential Evolution Sampler.");
      50             : 
      51          50 :   if (_num_parallel_proposals < _priors.size())
      52           0 :     mooseWarning(
      53             :         "It is recommended that the parallel proposals be greater than or equal to the "
      54             :         "inferred parameters. This will allow the sampler to not get stuck on a hyper-plane.");
      55             : 
      56         100 :   if (isParamValid("scales"))
      57             :   {
      58           6 :     _scales = getParam<std::vector<Real>>("scales");
      59           2 :     if (_scales.size() != _priors.size())
      60           2 :       paramError("scales",
      61             :                  "The number of scales provided should match the number of tunable params.");
      62             :   }
      63             :   else
      64          48 :     _scales.assign(_priors.size(), 1.0);
      65          48 : }
      66             : 
      67             : void
      68        3888 : AffineInvariantDES::computeDifferential(
      69             :     const Real & state1, const Real & state2, const Real & rnd, const Real & scale, Real & diff)
      70             : {
      71             :   Real gamma;
      72             :   Real b;
      73        3888 :   tuneParams(gamma, b, scale);
      74        3888 :   diff = gamma * (state1 - state2) + Normal::quantile(rnd, 0.0, b);
      75        3888 : }
      76             : 
      77             : void
      78        3888 : AffineInvariantDES::tuneParams(Real & gamma, Real & b, const Real & scale)
      79             : {
      80        3888 :   if (_tuning_option == "Braak2006_static")
      81             :   {
      82        3888 :     gamma = 2.38 / std::sqrt(2 * _priors.size());
      83        3888 :     b = 1e-6 * scale;
      84             :   }
      85        3888 : }
      86             : 
      87             : void
      88         288 : AffineInvariantDES::proposeSamples()
      89             : {
      90         288 :   unsigned int j = 0;
      91             :   bool indicator;
      92             :   std::pair<unsigned int, unsigned int> index_req;
      93             :   Real diff;
      94        2112 :   while (j < _num_parallel_proposals)
      95             :   {
      96             :     indicator = 0;
      97        1536 :     index_req = randomIndexPair(_num_parallel_proposals, j);
      98        4608 :     for (unsigned int i = 0; i < _priors.size(); ++i)
      99             :     {
     100        3072 :       computeDifferential(_previous_state[index_req.first][i],
     101        3072 :                           _previous_state[index_req.second][i],
     102        3072 :                           random(),
     103             :                           _scales[i],
     104             :                           diff);
     105        4032 :       _new_samples[j][i] = (_t_step + 1 > decisionStep()) ? (_previous_state[j][i] + diff)
     106         960 :                                                           : _priors[i]->quantile(random());
     107        3072 :       if (_lower_bound)
     108           0 :         indicator =
     109           0 :             (_new_samples[j][i] < (*_lower_bound)[i] || _new_samples[j][i] > (*_upper_bound)[i])
     110           0 :                 ? 1
     111             :                 : indicator;
     112             :     }
     113        1536 :     if (_var_prior)
     114             :     {
     115         816 :       computeDifferential(_previous_state_var[index_req.first],
     116         816 :                           _previous_state_var[index_req.second],
     117         816 :                           random(),
     118         816 :                           1.0,
     119             :                           diff);
     120        1056 :       _new_var_samples[j] = (_t_step + 1 > decisionStep()) ? (_previous_state_var[j] + diff)
     121         240 :                                                            : _var_prior->quantile(random());
     122         816 :       if (_new_var_samples[j] < 0.0 || _new_var_samples[j] > _variance_bound)
     123             :         indicator = 1;
     124             :     }
     125        1536 :     if (!indicator)
     126        1440 :       ++j;
     127             :   }
     128         288 : }