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         196 : AffineInvariantDES::validParams()
      25             : {
      26         196 :   InputParameters params = PMCMCBase::validParams();
      27         196 :   params.addClassDescription("Perform Affine Invariant Ensemble MCMC with differential sampler.");
      28         392 :   params.addRequiredParam<ReporterName>(
      29             :       "previous_state", "Reporter value with the previous state of all the walkers.");
      30         392 :   params.addRequiredParam<ReporterName>(
      31             :       "previous_state_var",
      32             :       "Reporter value with the previous state of all the walkers for variance.");
      33         392 :   MooseEnum tuning_option("Braak2006_static", "Braak2006_static");
      34         392 :   params.addParam<MooseEnum>(
      35             :       "tuning_option", tuning_option, "The tuning option for internal parameters.");
      36         392 :   params.addParam<std::vector<Real>>("scales", "Scales for the parameters.");
      37         196 :   return params;
      38         196 : }
      39             : 
      40         116 : AffineInvariantDES::AffineInvariantDES(const InputParameters & parameters)
      41             :   : PMCMCBase(parameters),
      42         116 :     _previous_state(getReporterValue<std::vector<std::vector<Real>>>("previous_state")),
      43         116 :     _previous_state_var(getReporterValue<std::vector<Real>>("previous_state_var")),
      44         348 :     _tuning_option(getParam<MooseEnum>("tuning_option"))
      45             : {
      46         116 :   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         116 :   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         232 :   if (isParamValid("scales"))
      57             :   {
      58          12 :     _scales = getParam<std::vector<Real>>("scales");
      59           4 :     if (_scales.size() != _priors.size())
      60           4 :       paramError("scales",
      61             :                  "The number of scales provided should match the number of tunable params.");
      62             :   }
      63             :   else
      64         112 :     _scales.assign(_priors.size(), 1.0);
      65         112 : }
      66             : 
      67             : void
      68        5360 : 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        5360 :   tuneParams(gamma, b, scale);
      74        5360 :   diff = gamma * (state1 - state2) + Normal::quantile(rnd, 0.0, b);
      75        5360 : }
      76             : 
      77             : void
      78        5360 : AffineInvariantDES::tuneParams(Real & gamma, Real & b, const Real & scale)
      79             : {
      80        5360 :   if (_tuning_option == "Braak2006_static")
      81             :   {
      82        5360 :     gamma = 2.38 / std::sqrt(2 * _priors.size());
      83        5360 :     b = 1e-6 * scale;
      84             :   }
      85        5360 : }
      86             : 
      87             : void
      88         400 : AffineInvariantDES::proposeSamples(const unsigned int seed_value)
      89             : {
      90         400 :   unsigned int j = 0;
      91             :   bool indicator;
      92             :   unsigned int index_req1, index_req2;
      93             :   Real diff;
      94        2920 :   while (j < _num_parallel_proposals)
      95             :   {
      96             :     indicator = 0;
      97        2120 :     randomIndexPair(_num_parallel_proposals, j, seed_value, index_req1, index_req2);
      98        6360 :     for (unsigned int i = 0; i < _priors.size(); ++i)
      99             :     {
     100        4240 :       computeDifferential(_previous_state[index_req1][i],
     101        4240 :                           _previous_state[index_req2][i],
     102        4240 :                           getRand(seed_value),
     103             :                           _scales[i],
     104             :                           diff);
     105        5840 :       _new_samples[j][i] = (_t_step > decisionStep()) ? (_previous_state[j][i] + diff)
     106        1600 :                                                       : _priors[i]->quantile(getRand(seed_value));
     107        4240 :       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        2120 :     if (_var_prior)
     114             :     {
     115        1120 :       computeDifferential(_previous_state_var[index_req1],
     116        1120 :                           _previous_state_var[index_req2],
     117        1120 :                           getRand(seed_value),
     118        1120 :                           1.0,
     119             :                           diff);
     120        1520 :       _new_var_samples[j] = (_t_step > decisionStep()) ? (_previous_state_var[j] + diff)
     121         400 :                                                        : _var_prior->quantile(getRand(seed_value));
     122        1120 :       if (_new_var_samples[j] < 0.0)
     123             :         indicator = 1;
     124             :     }
     125        2120 :     if (!indicator)
     126        2000 :       ++j;
     127             :   }
     128         400 : }