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 "LatinHypercubeSampler.h"
      11             : #include "Distribution.h"
      12             : #include "MooseRandomPerturbation.h"
      13             : 
      14             : registerMooseObjectAliased("StochasticToolsApp", LatinHypercubeSampler, "LatinHypercube");
      15             : 
      16             : InputParameters
      17         741 : LatinHypercubeSampler::validParams()
      18             : {
      19         741 :   InputParameters params = Sampler::validParams();
      20         741 :   params.addClassDescription("Latin Hypercube Sampler.");
      21        1482 :   params.addRequiredParam<dof_id_type>("num_rows", "The size of the square matrix to generate.");
      22        1482 :   params.addRequiredParam<std::vector<DistributionName>>(
      23             :       "distributions",
      24             :       "The distribution names to be sampled, the number of distributions provided defines the "
      25             :       "number of columns per matrix.");
      26         741 :   return params;
      27           0 : }
      28             : 
      29         401 : LatinHypercubeSampler::LatinHypercubeSampler(const InputParameters & parameters)
      30         401 :   : Sampler(parameters)
      31             : {
      32         802 :   const auto & distribution_names = getParam<std::vector<DistributionName>>("distributions");
      33        1663 :   for (const DistributionName & name : distribution_names)
      34        1262 :     _distributions.push_back(&getDistributionByName(name));
      35             : 
      36        1203 :   setNumberOfRows(getParam<dof_id_type>("num_rows"));
      37         401 :   setNumberOfCols(distribution_names.size());
      38             :   // Generator 0: within-bin uniform draws. Generator 1: column shuffler seeds.
      39         401 :   setNumberOfRandomSeeds(2);
      40         401 : }
      41             : 
      42             : void
      43         473 : LatinHypercubeSampler::executeTearDown()
      44             : {
      45         473 :   _shufflers.clear();
      46        1855 :   for (const auto col : make_range(getNumberOfCols()))
      47             :   {
      48        1382 :     const auto seed = getRandl(col, 0, std::numeric_limits<uint32_t>::max(), 1);
      49        2764 :     _shufflers.push_back(std::make_unique<MooseRandomPerturbation>(seed, getNumberOfRows()));
      50             :   }
      51         473 : }
      52             : 
      53             : Real
      54       31868 : LatinHypercubeSampler::computeSample(dof_id_type row_index, dof_id_type col_index)
      55             : {
      56             :   mooseAssert(_shufflers.size() > 0, "Shufflers have not been initialized.");
      57             : 
      58             :   // Divide [0,1] into N equal bins of width 1/N.
      59       31868 :   const Real bin_size = 1. / getNumberOfRows();
      60             : 
      61             :   // Map row_index to a shuffled bin via the column's bijective permutation.
      62             :   // Because permute() is a bijection on [0, N), each row gets a distinct bin,
      63             :   // which is the core LHS stratification guarantee.
      64       31868 :   const auto bin = _shufflers[col_index]->permute(row_index);
      65             : 
      66             :   // Draw a uniform random point within the selected bin.
      67       31868 :   const auto lower = bin * bin_size;
      68       31868 :   const auto upper = (bin + 1) * bin_size;
      69             :   const Real probability =
      70       31868 :       getRand(row_index * getNumberOfCols() + col_index) * (upper - lower) + lower;
      71             : 
      72             :   // Transform the probability through the inverse CDF to obtain the sample value.
      73       31868 :   return _distributions[col_index]->quantile(probability);
      74             : }