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 "MorrisSampler.h"
      11             : #include "Distribution.h"
      12             : 
      13             : registerMooseObject("StochasticToolsApp", MorrisSampler);
      14             : 
      15             : InputParameters
      16         163 : MorrisSampler::validParams()
      17             : {
      18         163 :   InputParameters params = Sampler::validParams();
      19         163 :   params.addClassDescription("Morris variance-based sensitivity analysis Sampler.");
      20         326 :   params.addRequiredParam<std::vector<DistributionName>>(
      21             :       "distributions",
      22             :       "The distribution names to be sampled, the number of distributions provided defines the "
      23             :       "number of columns per matrix.");
      24         326 :   params.addRequiredRangeCheckedParam<dof_id_type>(
      25             :       "trajectories",
      26             :       "trajectories > 0",
      27             :       "Number of unique trajectories to perform. The higher number of these usually means a more "
      28             :       "accurate sensitivity evaluation, but it is proportional to the number of required model "
      29             :       "evaluations: 'trajectoris' x (number of 'distributions' + 1).");
      30         489 :   params.addRangeCheckedParam<unsigned int>(
      31             :       "levels",
      32         326 :       4,
      33             :       "levels % 2 = 0 & levels > 0",
      34             :       "The number of levels in the sampling. This determines the discretization of the input "
      35             :       "space, more levels means finer discretization and more possible model perturbations.");
      36         163 :   return params;
      37           0 : }
      38             : 
      39          88 : MorrisSampler::MorrisSampler(const InputParameters & parameters)
      40             :   : Sampler(parameters),
      41          88 :     _num_trajectories(getParam<dof_id_type>("trajectories")),
      42         264 :     _num_levels(getParam<unsigned int>("levels"))
      43             : 
      44             : {
      45         536 :   for (const auto & name : getParam<std::vector<DistributionName>>("distributions"))
      46         360 :     _distributions.push_back(&getDistributionByName(name));
      47             :   const dof_id_type nc = _distributions.size();
      48             : 
      49          88 :   setNumberOfCols(nc);
      50          88 :   setNumberOfRows(_num_trajectories * (nc + 1));
      51             : 
      52          88 :   _b = RealEigenMatrix::Ones(nc + 1, nc).triangularView<Eigen::StrictlyLower>();
      53          88 :   _pstar.resize(nc, nc);
      54          88 :   _j.setOnes(nc + 1, nc);
      55          88 :   _dstar.resize(nc, nc);
      56          88 :   _xstar.resize(nc + 1, nc);
      57          88 :   _bstar.resize(nc + 1, nc);
      58          88 : }
      59             : 
      60             : Real
      61      878400 : MorrisSampler::computeSample(dof_id_type row_index, dof_id_type col_index)
      62             : {
      63      878400 :   const dof_id_type traj = row_index / (getNumberOfCols() + 1);
      64      878400 :   const dof_id_type traj_ind = row_index % (getNumberOfCols() + 1);
      65      878400 :   if (traj != _curr_trajectory)
      66             :   {
      67       20970 :     _curr_trajectory = traj;
      68       20970 :     updateBstar();
      69             :   }
      70      878400 :   return _distributions[col_index]->quantile(_bstar(traj_ind, col_index));
      71             : }
      72             : 
      73             : void
      74       20970 : MorrisSampler::updateBstar()
      75             : {
      76             :   mooseAssert(_curr_trajectory < _num_trajectories,
      77             :               "Current trajectory index is greater than the prescribed number of trajectories.");
      78             : 
      79       20970 :   const dof_id_type nc = getNumberOfCols(); // convenience
      80       20970 :   dof_id_type rn_ind = _curr_trajectory * nc * (nc + 1);
      81             : 
      82       20970 :   _pstar.setZero();
      83             :   // Which parameter to perturb
      84       20970 :   std::vector<dof_id_type> pchoice(nc);
      85             :   std::iota(pchoice.begin(), pchoice.end(), 0);
      86      146520 :   for (dof_id_type c = 0; c < nc; ++c)
      87             :   {
      88      125550 :     const unsigned int ind = nc > 1 ? getRandl(rn_ind++, 0, pchoice.size()) : 0;
      89      125550 :     _pstar(pchoice[ind], c) = 1.0;
      90             :     pchoice.erase(pchoice.begin() + ind);
      91             :   }
      92             : 
      93       20970 :   _dstar.setZero();
      94             :   // Direction of perturbation
      95      146520 :   for (dof_id_type c = 0; c < nc; ++c)
      96      189800 :     _dstar(c, c) = getRand(rn_ind++) < 0.5 ? -1.0 : 1.0;
      97             : 
      98             :   // Initial value
      99      146520 :   for (dof_id_type c = 0; c < nc; ++c)
     100             :   {
     101      125550 :     const auto lind = getRandl(rn_ind++, 0, _num_levels / 2);
     102      125550 :     _xstar.col(c).setConstant((Real)lind * 1.0 / ((Real)_num_levels - 1));
     103             :   }
     104             : 
     105             :   _bstar =
     106       20970 :       _xstar + _num_levels / 4.0 / (_num_levels - 1) * ((2.0 * _b * _pstar - _j) * _dstar + _j);
     107       20970 : }
     108             : 
     109             : LocalRankConfig
     110         176 : MorrisSampler::constructRankConfig(bool batch_mode) const
     111             : {
     112         176 :   std::vector<LocalRankConfig> all_rc(processor_id() + 1);
     113         528 :   for (processor_id_type r = 0; r <= processor_id(); ++r)
     114         352 :     all_rc[r] = rankConfig(
     115         352 :         r, n_processors(), _num_trajectories, _min_procs_per_row, _max_procs_per_row, batch_mode);
     116             :   LocalRankConfig & rc = all_rc.back();
     117             : 
     118         176 :   rc.num_local_sims *= _distributions.size() + 1;
     119             :   bool found_first = false;
     120         528 :   for (auto it = all_rc.rbegin(); it != all_rc.rend(); ++it)
     121         352 :     if (it->is_first_local_rank)
     122             :     {
     123         304 :       if (found_first)
     124         128 :         rc.first_local_sim_index += it->num_local_sims * _distributions.size();
     125             :       else
     126             :         found_first = true;
     127             :     }
     128             : 
     129         176 :   if (!batch_mode)
     130             :   {
     131          88 :     rc.num_local_apps = rc.num_local_sims;
     132          88 :     rc.first_local_app_index = rc.first_local_sim_index;
     133             :   }
     134             : 
     135         176 :   return rc;
     136         176 : }