Redesign of MC working

(refs idaholab#13906)

framework/include/samplers/Sampler.h

@@ -153,12 +153,12 @@ class Sampler : public MooseObject, public SetupInterface, public DistributionIn
    *
    * @return The list of samples for the Sampler.
    */
-  virtual std::vector<DenseMatrix<Real>> sample() = 0;
+  //virtual std::vector<DenseMatrix<Real>> sample() = 0;
 
 
   virtual std::vector<DenseMatrix<Real>> computeSampleMatrices();
   //virtual DenseMatrix<Real> createSample(const dof_id_type & matrix_index);
-  virtual Real computeSample(const dof_id_type & matrix_index, const dof_id_type & row_index, const dof_id_type & col_index) = 0;
+  virtual Real computeSample(dof_id_type matrix_index, dof_id_type row_index, dof_id_type col_index) = 0;
 
 
   /**
@@ -167,7 +167,9 @@ class Sampler : public MooseObject, public SetupInterface, public DistributionIn
    * of child objects.
    * @param number The required number of random seeds, by default this is called with 1.
    */
-  void setNumberOfRequiedRandomSeeds(const std::size_t & number);
+  void setNumberOfRandomSeeds(std::size_t number);
+  void setNumberOfMatrices(dof_id_type number);
+
 
   /**
    * Reinitialize the offsets and row counts.
@@ -185,6 +187,9 @@ class Sampler : public MooseObject, public SetupInterface, public DistributionIn
   std::vector<std::string> _sample_names;
 
 private:
+
+  void reinit();
+
   /// Random number generator, don't give users access we want to control it via the interface
   /// from this class.
   MooseRandom _generator;
@@ -207,13 +212,18 @@ class Sampler : public MooseObject, public SetupInterface, public DistributionIn
   /// Global row index for end of data for this processor
   dof_id_type _local_row_end;
 
-  const dof_id_type _num_matrices;
+protected:
+  dof_id_type _num_matrices;
   const dof_id_type _num_rows_per_matrix;
   const dof_id_type _num_cols_per_matrix;
 
   /// Total number of rows
-  const dof_id_type _total_rows;
-
+  dof_id_type _total_rows;
 
+  // Create protected versions that are const references to these instead of passing the items in???
+  //dof_id_type _private_matrix_index;
+  //dof_id_type _private_row_index;
+  //dof_id_type _private_col_index;
 
+  bool _rand_called = false;
 };

framework/src/samplers/Sampler.C

@@ -57,17 +57,20 @@ Sampler::Sampler(const InputParameters & parameters)
 {
   for (const DistributionName & name : _distribution_names)
     _distributions.push_back(&getDistributionByName(name));
-  setNumberOfRequiedRandomSeeds(1);
+  setNumberOfRandomSeeds(1);
+  reinit();
+}
 
+void
+Sampler::reinit()
+{
+  _total_rows = _num_matrices * _num_rows_per_matrix;
   MooseUtils::linearPartitionItems(
       _total_rows, n_processors(), processor_id(), _local_rows, _local_row_begin, _local_row_end);
 
-  if (_sample_names.empty())
-  {
-    _sample_names.resize(_num_matrices);
-    for (dof_id_type i = 0; i < _num_matrices; ++i)
-      _sample_names[i] = "sample_" + std::to_string(i);
-  }
+  _sample_names.resize(_num_matrices);
+  for (dof_id_type i = 0; i < _num_matrices; ++i)
+    _sample_names[i] = "sample_" + std::to_string(i);
 }
 
 void
@@ -123,7 +126,6 @@ Sampler::computeSampleMatrices()
   }
 
   sampleTearDown();
-
   return output;
 }
 
@@ -135,7 +137,7 @@ Sampler::rand(const unsigned int index)
 }
 
 void
-Sampler::setNumberOfRequiedRandomSeeds(const std::size_t & number)
+Sampler::setNumberOfRandomSeeds(std::size_t number)
 {
   if (number == 0)
     mooseError("The number of seeds must be larger than zero.");
@@ -150,6 +152,14 @@ Sampler::setNumberOfRequiedRandomSeeds(const std::size_t & number)
   _generator.saveState();
 }
 
+void
+Sampler::setNumberOfMatrices(dof_id_type number)
+{
+  _num_matrices;
+  reinit();
+}
+
+
 void
 Sampler::setSampleNames(const std::vector<std::string> & names)
 {

modules/stochastic_tools/include/samplers/MonteCarloSampler.h

@@ -24,11 +24,6 @@ class MonteCarloSampler : public Sampler
   MonteCarloSampler(const InputParameters & parameters);
 
 protected:
-  virtual std::vector<DenseMatrix<Real>> sample() override;
-
-  /// Number of matrices
-  const dof_id_type _num_matrices;
-
-  /// Number of monte carlo samples to create for each distribution
-  const dof_id_type _num_samples;
+  //virtual std::vector<DenseMatrix<Real>> sample() override;
+  virtual Real computeSample(dof_id_type matrix_index, dof_id_type row_index, dof_id_type col_index);
 };

modules/stochastic_tools/include/samplers/SobolSampler.h

@@ -24,17 +24,15 @@ class SobolSampler : public Sampler
   SobolSampler(const InputParameters & parameters);
 
 protected:
-  virtual std::vector<DenseMatrix<Real>> sample() override;
+  //virtual std::vector<DenseMatrix<Real>> sample() override;
+
+  virtual Real computeSample(dof_id_type matrix_index, dof_id_type row_index, dof_id_type col_index) override;
   virtual void sampleSetUp() override;
   virtual void sampleTearDown() override;
 
-  /// Number of Monte Carlo samples to create for each Sobol matrix
-  const std::size_t _num_samples;
-
   ///@{
   /// Sobol Monte Carlo matrices, these are sized and cleared to avoid keeping large matrices around
   DenseMatrix<Real> _a_matrix;
   DenseMatrix<Real> _b_matrix;
   ///@}
 };
-

modules/stochastic_tools/src/samplers/MonteCarloSampler.C

@@ -17,21 +17,23 @@ validParams<MonteCarloSampler>()
 {
   InputParameters params = validParams<Sampler>();
   params.addClassDescription("Monte Carlo Sampler.");
-  params.addRequiredParam<dof_id_type>(
-      "n_samples",
-      "Number of Monte Carlo samples to perform for each distribution within each matrix.");
-  params.addParam<dof_id_type>(
-      "n_matrices", 1, "Number of matrices to create, each matrix will contain 'n_samples'.");
   return params;
 }
 
 MonteCarloSampler::MonteCarloSampler(const InputParameters & parameters)
-  : Sampler(parameters),
-    _num_matrices(getParam<dof_id_type>("n_matrices")),
-    _num_samples(getParam<dof_id_type>("n_samples"))
+  : Sampler(parameters)
 {
 }
 
+Real
+MonteCarloSampler::computeSample(dof_id_type matrix_index, dof_id_type row_index, dof_id_type col_index)
+{
+  return _distributions[col_index]->quantile(rand());
+}
+
+
+
+/*
 std::vector<DenseMatrix<Real>>
 MonteCarloSampler::sample()
 {
@@ -45,3 +47,4 @@ MonteCarloSampler::sample()
   }
   return output;
 }
+*/

modules/stochastic_tools/src/samplers/SobolSampler.C

@@ -17,40 +17,54 @@ validParams<SobolSampler>()
 {
   InputParameters params = validParams<Sampler>();
   params.addClassDescription("Sobol variance-based sensitivity analysis Sampler.");
-  params.addRequiredParam<unsigned int>(
-      "n_samples", "Number of Monte Carlo samples to perform for each distribution.");
   return params;
 }
 
 SobolSampler::SobolSampler(const InputParameters & parameters)
   : Sampler(parameters),
-    _num_samples(getParam<unsigned int>("n_samples")),
     _a_matrix(0, 0),
     _b_matrix(0, 0)
 {
-  setNumberOfRequiedRandomSeeds(2);
+  setNumberOfRandomSeeds(2);
+  setNumberOfMatrices(_num_cols_per_matrix + 2);
 }
 
 void
 SobolSampler::sampleSetUp()
 {
-  _a_matrix.resize(_num_samples, _distributions.size());
-  _b_matrix.resize(_num_samples, _distributions.size());
-  for (std::size_t i = 0; i < _num_samples; ++i)
-    for (MooseIndex(_distributions) j = 0; j < _distributions.size(); ++j)
+  _a_matrix.resize(_num_rows_per_matrix, _num_cols_per_matrix);
+  _b_matrix.resize(_num_rows_per_matrix, _num_cols_per_matrix);
+  for (std::size_t i = 0; i < _num_cols_per_matrix; ++i)
+    for (dof_id_type j = 0; j < _num_cols_per_matrix; ++j)
     {
       _a_matrix(i, j) = _distributions[j]->quantile(this->rand(0));
       _b_matrix(i, j) = _distributions[j]->quantile(this->rand(1));
     }
 }
 
+Real
+SobolSampler::computeSample(dof_id_type matrix_index, dof_id_type row_index, dof_id_type col_index)
+{
+  if (matrix_index == 0)
+    return _a_matrix(row_index, col_index);
+  else if (matrix_index == 1)
+    return _b_matrix(row_index, col_index);
+  else if (col_index == matrix_index - 2)
+    return _b_matrix(row_index, col_index);
+  else
+    return _a_matrix(row_index, col_index);
+
+}
+
 void
 SobolSampler::sampleTearDown()
 {
   _a_matrix.resize(0, 0);
   _b_matrix.resize(0, 0);
 }
 
+
+/*
 std::vector<DenseMatrix<Real>>
 SobolSampler::sample()
 {
@@ -72,3 +86,4 @@ SobolSampler::sample()
 
   return output;
 }
+*/

modules/stochastic_tools/test/tests/samplers/monte_carlo/monte_carlo_uniform.i

@@ -21,7 +21,8 @@
 [Samplers]
   [./sample]
     type = MonteCarloSampler
-    n_samples = 10
+    #n_samples = 10: TODO: deprecate
+    num_rows = 10
     distributions = 'uniform'
     execute_on = 'initial timestep_end'
   [../]

modules/stochastic_tools/test/tests/samplers/monte_carlo/monte_carlo_weibull.i

@@ -21,7 +21,7 @@
 [Samplers]
   [./sample]
     type = MonteCarloSampler
-    n_samples = 10
+    num_rows = 10
     distributions = 'weibull'
     execute_on = 'initial timestep_end'
   [../]