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         382 : MorrisSampler::validParams()
      17             : {
      18         382 :   InputParameters params = Sampler::validParams();
      19         382 :   params.addClassDescription("Morris variance-based sensitivity analysis Sampler.");
      20         764 :   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         764 :   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        1146 :   params.addRangeCheckedParam<unsigned int>(
      31             :       "levels",
      32         764 :       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         382 :   return params;
      37           0 : }
      38             : 
      39         222 : MorrisSampler::MorrisSampler(const InputParameters & parameters)
      40             :   : Sampler(parameters),
      41         222 :     _num_trajectories(getParam<dof_id_type>("trajectories")),
      42         666 :     _num_levels(getParam<unsigned int>("levels"))
      43             : 
      44             : {
      45        1374 :   for (const auto & name : getParam<std::vector<DistributionName>>("distributions"))
      46         930 :     _distributions.push_back(&getDistributionByName(name));
      47             : 
      48         222 :   setNumberOfCols(_distributions.size());
      49         222 :   setNumberOfRows(_num_trajectories * (_distributions.size() + 1));
      50         222 : }
      51             : 
      52             : void
      53         364 : MorrisSampler::sampleSetUp(const Sampler::SampleMode /*mode*/)
      54             : {
      55         364 :   const dof_id_type nc = getNumberOfCols();
      56         364 :   _b = RealEigenMatrix::Ones(nc + 1, nc).triangularView<Eigen::StrictlyLower>();
      57         364 :   _pstar.resize(nc, nc);
      58         364 :   _j.setOnes(nc + 1, nc);
      59         364 :   _dstar.resize(nc, nc);
      60         364 :   _xstar.resize(nc + 1, nc);
      61         364 :   _bstar.resize(nc + 1, nc);
      62         364 : }
      63             : 
      64             : Real
      65     1756800 : MorrisSampler::computeSample(dof_id_type row_index, dof_id_type col_index)
      66             : {
      67     1756800 :   const dof_id_type traj_ind = row_index % (getNumberOfCols() + 1);
      68     1756800 :   if (traj_ind == 0 && col_index == 0)
      69       42000 :     updateBstar();
      70     1756800 :   return _distributions[col_index]->quantile(_bstar(traj_ind, col_index));
      71             : }
      72             : 
      73             : void
      74       42000 : MorrisSampler::updateBstar()
      75             : {
      76       42000 :   const dof_id_type nc = getNumberOfCols(); // convenience
      77       42000 :   _pstar.setZero();
      78             :   // Which parameter to perturb
      79       42000 :   std::vector<dof_id_type> pchoice(nc);
      80             :   std::iota(pchoice.begin(), pchoice.end(), 0);
      81      293280 :   for (dof_id_type c = 0; c < nc; ++c)
      82             :   {
      83      251280 :     const unsigned int ind = nc > 1 ? getRandl(0, 0, pchoice.size()) : 0;
      84      251280 :     _pstar(pchoice[ind], c) = 1.0;
      85             :     pchoice.erase(pchoice.begin() + ind);
      86             :   }
      87             : 
      88       42000 :   _dstar.setZero();
      89             :   // Direction of perturbation
      90      293280 :   for (dof_id_type c = 0; c < nc; ++c)
      91      379900 :     _dstar(c, c) = getRand() < 0.5 ? -1.0 : 1.0;
      92             : 
      93             :   // Initial value
      94      293280 :   for (dof_id_type c = 0; c < nc; ++c)
      95             :   {
      96      251280 :     const auto lind = getRandl(0, 0, _num_levels / 2);
      97      251280 :     _xstar.col(c).setConstant((Real)lind * 1.0 / ((Real)_num_levels - 1));
      98             :   }
      99             : 
     100             :   _bstar =
     101       42000 :       _xstar + _num_levels / 4.0 / (_num_levels - 1) * ((2.0 * _b * _pstar - _j) * _dstar + _j);
     102             : 
     103             :   // This matrix represent _n_cols * (_n_cols + 1) samples, but so far we have only
     104             :   // advanced the generator 3 * _n_cols times. For the generator state restore
     105             :   // to work properly, we need to finish advancing the generator
     106       42000 :   if (nc > 2)
     107       42000 :     advanceGenerator(0, nc * (nc - 2));
     108       42000 : }
     109             : 
     110             : LocalRankConfig
     111         444 : MorrisSampler::constructRankConfig(bool batch_mode) const
     112             : {
     113         444 :   std::vector<LocalRankConfig> all_rc(processor_id() + 1);
     114        1332 :   for (processor_id_type r = 0; r <= processor_id(); ++r)
     115         888 :     all_rc[r] = rankConfig(
     116         888 :         r, n_processors(), _num_trajectories, _min_procs_per_row, _max_procs_per_row, batch_mode);
     117             :   LocalRankConfig & rc = all_rc.back();
     118             : 
     119         444 :   rc.num_local_sims *= _distributions.size() + 1;
     120             :   bool found_first = false;
     121        1332 :   for (auto it = all_rc.rbegin(); it != all_rc.rend(); ++it)
     122         888 :     if (it->is_first_local_rank)
     123             :     {
     124         776 :       if (found_first)
     125         332 :         rc.first_local_sim_index += it->num_local_sims * _distributions.size();
     126             :       else
     127             :         found_first = true;
     128             :     }
     129             : 
     130         444 :   if (!batch_mode)
     131             :   {
     132         222 :     rc.num_local_apps = rc.num_local_sims;
     133         222 :     rc.first_local_app_index = rc.first_local_sim_index;
     134             :   }
     135             : 
     136         444 :   return rc;
     137             : }