Surrogates System

Overview

The Stochastic Tools module contains a system for training and running surrogate models. These objects inherit from the SurrogateModel class and are added to a simulation using the [Surrogates] block in the input file. SurrogateModel objects are standalone utilities designed to be directly called by other objects in similar fashion to MOOSE Function objects, refer to Using a SurrogateModel for more information.

Creating a SurrogateModel

A model is created by inheriting from SurrogateModel and overriding the evaluate method. This same method is also what should be called by other objects requiring data. However, the model classes can implement an arbitrary interface as needed.

Using a SurrogateModel

Model objects are obtained by other objects using the getSurrogateModel and getSurrogateModelByName methods. The first expects an input parameter name from which the desired object name is determined. The later function simply accepts the actual name of the object.

For example, the PolynomialChaosReporter object requires a PolynomialChaos object. In the header a reference to the desired model is declare and in the source this reference is initialized with a get method.

std::vector<const PolynomialChaos *> _pc;

_pc.push_back(&getSurrogateModelByName<PolynomialChaos>(nm));

Getting Training Data

Training data computed by a trainer object is gathered using the getModelData method. It has a template argument that gives the data type and a single argument, the name of the training data. This name should match the name given in the associated training object(s). For example, the "order" mentioned in the PolynomialChaosTrainer is needed in the actual model. In the header a constant reference to the desired type is declared and this reference is initialized in the source.

const unsigned int & _order;

_order(getModelData<unsigned int>("_order")),

The training data can be supplied in one of two ways: using a file or directly from a trainer.

File based data initialization happens by using the "filename". The supplied file must have declared data with the same names and types as the get methods used within this class. Refer to SurrogateTrainerOutput for more information.

It is also possible to use data directly from a trainer object by using the "trainer" parameter. This method directly shares the same training data, so if the trainer updates the data the model has a reference to the same data.

Example Input File Syntax

The following input file snippet adds a PolynomialChaos surrogate model for evaluation. Please refer to the documentation on the individual models for more details.

[Surrogates]
  [poly_chaos]
    type = PolynomialChaos
    trainer = poly_chaos
  []
[]

Available Objects

Stochastic Tools App
GaussianProcessSurrogateComputes and evaluates Gaussian Process surrogate model.
NearestPointSurrogateSurrogate that evaluates the value from the nearest point from data in NearestPointTrainer
PODReducedBasisSurrogateEvaluates POD-RB surrogate model with reduced operators computed from PODReducedBasisTrainer.
PolynomialChaosComputes and evaluates polynomial chaos surrogate model.
PolynomialRegressionSurrogateEvaluates polynomial regression model with coefficients computed from PolynomialRegressionTrainer.

Available Actions

Stochastic Tools App
AddSurrogateActionAdds SurrogateTrainer and SurrogateModel objects contained within the [Trainers] and [Surrogates] input blocks.

(contrib/moose/modules/stochastic_tools/include/reporters/PolynomialChaosReporter.h)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "GeneralReporter.h"
#include "SurrogateModelInterface.h"

#include "PolynomialChaos.h"
#include "nlohmann/json.h"

class PolynomialChaosReporter : public GeneralReporter, public SurrogateModelInterface
{
public:
  static InputParameters validParams();
  PolynomialChaosReporter(const InputParameters & parameters);

  virtual void initialize() override {}
  virtual void execute() override;
  virtual void finalize() override {}

private:
  /// Helper function for computing local sensitivity from a polynomial chaos model
  static std::vector<Real> computeLocalSensitivity(const PolynomialChaos & pc,
                                                   const std::vector<Real> & data);

  /// Points for local sensitivity calculation
  const std::vector<std::vector<Real>> & _loc_point;
  /// Samplers for local sensitivity calculation
  std::vector<Sampler *> _loc_sampler;

  /// Polynomial chaos models
  std::vector<const PolynomialChaos *> _pc;

  /// Local sensitivity from specified points
  std::vector<std::vector<std::vector<Real>> *> _loc_point_sense;
  /// Local sensitivity from sampled points
  std::vector<std::vector<std::vector<Real>> *> _loc_sampler_sense;
};

void to_json(nlohmann::json & json, const PolynomialChaos * const & pc);

/**
 * PCStatisticsContext is almost identical to ReporterStatisticsContext with
 * InType == Outype. Unfortunately, we cannot derive from ReporterStatisticsContext
 * since that class relies on the construction of a Calculator object, something
 * that is unnecessary for calculating statistics with polynomial chaos.
 */
template <typename OutType>
class PCStatisticsContext : public ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>
{
public:
  PCStatisticsContext(const libMesh::ParallelObject & other,
                      const MooseObject & producer,
                      ReporterState<std::pair<OutType, std::vector<OutType>>> & state,
                      const PolynomialChaos & pc,
                      const MooseEnumItem & stat);

  virtual void finalize() override;
  virtual void storeInfo(nlohmann::json & json) const override;

private:
  /// Polynomial chaos surrogate object
  const PolynomialChaos & _pc;
  /// The stat to compute
  const MooseEnumItem _stat;
};

/**
 * PCSobolContext is almost identical to SobolReporterContext with
 * InType == Outype. Unfortunately, we cannot derive from SobolReporterContext
 * since that class relies on the construction of a Calculator object, something
 * that is unnecessary for calculating statistics with polynomial chaos.
 */
template <typename OutType>
class PCSobolContext : public ReporterGeneralContext<
                           std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>>>
{
public:
  PCSobolContext(
      const libMesh::ParallelObject & other,
      const MooseObject & producer,
      ReporterState<std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>>> & state,
      const PolynomialChaos & pc);

  virtual void finalize() override;
  virtual std::string type() const override
  {
    return "SobolIndices<" + MooseUtils::prettyCppType<OutType>() + ">";
  }

protected:
  virtual void store(nlohmann::json & json) const override;

private:
  /// Polynomial chaos surrogate object
  const PolynomialChaos & _pc;
};

(contrib/moose/modules/stochastic_tools/src/reporters/PolynomialChaosReporter.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "PolynomialChaosReporter.h"

#include "SobolReporter.h"
#include "StochasticReporter.h"

registerMooseObject("StochasticToolsApp", PolynomialChaosReporter);

InputParameters
PolynomialChaosReporter::validParams()
{
  InputParameters params = GeneralReporter::validParams();
  params.addClassDescription(
      "Tool for extracting data from PolynomialChaos surrogates and computing statistics.");
  params.addRequiredParam<std::vector<UserObjectName>>(
      "pc_name", "Name(s) of PolynomialChaos surrogate object(s).");
  params.addParam<bool>("include_data",
                        false,
                        "True to output information on the polynomial chaos model, including "
                        "polynomial types, orders, and coefficients.");
  MultiMooseEnum stats("mean=1 stddev=2 skewness=3 kurtosis=4", "");
  params.addParam<MultiMooseEnum>("statistics", stats, "Statistics to compute.");
  params.addParam<bool>("include_sobol", false, "True to compute Sobol indices.");
  params.addParam<std::vector<std::vector<Real>>>(
      "local_sensitivity_points",
      std::vector<std::vector<Real>>(0),
      "Points for each polynomial chaos surrogate specifying desired location of sensitivity "
      "measurement.");
  params.addParam<std::vector<SamplerName>>(
      "local_sensitivity_sampler",
      std::vector<SamplerName>(0),
      "Sampler for each polynomial chaos surrogate specifying desired location of sensitivity "
      "measurement.");
  return params;
}

PolynomialChaosReporter::PolynomialChaosReporter(const InputParameters & parameters)
  : GeneralReporter(parameters),
    SurrogateModelInterface(this),
    _loc_point(getParam<std::vector<std::vector<Real>>>("local_sensitivity_points"))
{
  for (const auto & sn : getParam<std::vector<SamplerName>>("local_sensitivity_sampler"))
    _loc_sampler.push_back(&getSamplerByName(sn));

  for (const auto & nm : getParam<std::vector<UserObjectName>>("pc_name"))
  {
    // Gather polynomial chaos surrogates
    _pc.push_back(&getSurrogateModelByName<PolynomialChaos>(nm));
    // PolynomialChaos data reporter
    if (getParam<bool>("include_data"))
    {
      auto & pc_ptr = declareValueByName<const PolynomialChaos *>(nm, REPORTER_MODE_ROOT);
      pc_ptr = _pc.back();
    }
    // Statistics reporter
    for (const auto & stat : getParam<MultiMooseEnum>("statistics"))
      declareValueByName<std::pair<Real, std::vector<Real>>, PCStatisticsContext<Real>>(
          nm + "_" + stat.name(), REPORTER_MODE_ROOT, *_pc.back(), stat);
    // Sobol reporter
    if (getParam<bool>("include_sobol"))
      declareValueByName<std::pair<std::vector<Real>, std::vector<std::vector<Real>>>,
                         PCSobolContext<Real>>(nm + "_SOBOL", REPORTER_MODE_ROOT, *_pc.back());
    // Local sensitivity points reporter
    if (_loc_point.size() > 0)
    {
      if (_loc_point.size() < _pc.size())
        paramError("local_sensitivity_points",
                   "There must be a set of points for each inputted polynomial chaos model.");
      _loc_point_sense.push_back(&declareValueByName<std::vector<std::vector<Real>>>(
          nm + "_POINT_SENSITIVITY", REPORTER_MODE_ROOT));
    }
    // Local sensitivity sampler reporter
    if (_loc_sampler.size() > 0)
    {
      if (_loc_sampler.size() < _pc.size())
        paramError("local_sensitivity_sampler",
                   "There must be a sampler for each inputted polynomial chaos model.");
      _loc_sampler_sense.push_back(
          &declareValueByName<std::vector<std::vector<Real>>,
                              StochasticReporterContext<std::vector<Real>>>(
              nm + "_SAMPLE_SENSITIVITY", REPORTER_MODE_ROOT, *_loc_sampler[_pc.size() - 1]));
    }
  }
}

void
PolynomialChaosReporter::execute()
{
  if ((_loc_point.size() + _loc_sampler.size()) == 0)
    return;

  for (const auto & i : index_range(_pc))
  {
    const auto & pc = *_pc[i];
    const auto nparam = pc.getNumberOfParameters();
    if (_loc_point.size() > 0)
    {
      if (_loc_point[i].size() % pc.getNumberOfParameters() != 0)
        paramError("local_sensitivity_points",
                   "Number of values must be divisible by number of parameters in "
                   "Polynomial Chaos model.");

      const auto npoint = _loc_point[i].size() / nparam;
      _loc_point_sense[i]->resize(npoint);
      for (const auto & p : make_range(npoint))
      {
        const std::vector<Real> data(_loc_point[i].begin() + p * nparam,
                                     _loc_point[i].begin() + (p + 1) * nparam);
        (*_loc_point_sense[i])[p] = computeLocalSensitivity(pc, data);
      }

      if (_loc_sampler.size() > 0)
      {
        if (_loc_sampler[i]->getNumberOfCols() != nparam)
          paramError(
              "local_sensitivity_sampler",
              "Sampler ",
              _loc_sampler[i]->name(),
              " does not have the same number of columns as the number of dimensions in model ",
              pc.name(),
              ".");
        for (const auto & p : make_range(_loc_sampler[i]->getNumberOfLocalRows()))
          (*_loc_sampler_sense[i])[p] =
              computeLocalSensitivity(pc, _loc_sampler[i]->getNextLocalRow());
      }
    }
  }
}

std::vector<Real>
PolynomialChaosReporter::computeLocalSensitivity(const PolynomialChaos & pc,
                                                 const std::vector<Real> & data)
{
  std::vector<Real> sense(data.size());
  const auto val = pc.evaluate(data);
  for (const auto & d : index_range(data))
    sense[d] = data[d] / val * pc.computeDerivative(d, data);
  return sense;
}

void
to_json(nlohmann::json & json, const PolynomialChaos * const & pc)
{
  pc->store(json);
}

template <typename OutType>
PCStatisticsContext<OutType>::PCStatisticsContext(
    const libMesh::ParallelObject & other,
    const MooseObject & producer,
    ReporterState<std::pair<OutType, std::vector<OutType>>> & state,
    const PolynomialChaos & pc,
    const MooseEnumItem & stat)
  : ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>(other, producer, state),
    _pc(pc),
    _stat(stat)
{
}

template <typename OutType>
void
PCStatisticsContext<OutType>::finalize()
{
  // Compute standard deviation
  OutType sig = OutType();
  if (_stat == "stddev" || _stat == "skewness" || _stat == "kurtosis")
    sig = _pc.computeStandardDeviation();

  OutType & val = this->_state.value().first;
  val = OutType();
  // Mean
  if (_stat == "mean" && this->processor_id() == 0)
    val = _pc.computeMean();
  // Standard Deviation
  else if (_stat == "stddev" && this->processor_id() == 0)
    val = sig;
  // Skewness
  else if (_stat == "skewness")
    val = _pc.powerExpectation(3) / (sig * sig * sig);
  // Kurtosis
  else if (_stat == "kurtosis")
    val = _pc.powerExpectation(4) / (sig * sig * sig * sig);
  this->_communicator.sum(val);

  ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>::finalize();
}

template <typename OutType>
void
PCStatisticsContext<OutType>::storeInfo(nlohmann::json & json) const
{
  ReporterGeneralContext<std::pair<OutType, std::vector<OutType>>>::storeInfo(json);
  json["stat"] = _stat.name();
}

template <typename OutType>
PCSobolContext<OutType>::PCSobolContext(
    const libMesh::ParallelObject & other,
    const MooseObject & producer,
    ReporterState<std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>>> & state,
    const PolynomialChaos & pc)
  : ReporterGeneralContext<std::pair<std::vector<OutType>, std::vector<std::vector<OutType>>>>(
        other, producer, state),
    _pc(pc)
{
}

template <typename OutType>
void
PCSobolContext<OutType>::finalize()
{
  const unsigned int nparam = _pc.getNumberOfParameters();
  std::vector<OutType> & val = this->_state.value().first;
  val.clear();

  // Compute variance
  auto var = _pc.computeStandardDeviation();
  var *= var;

  // First order
  for (const auto & i : make_range(nparam))
    val.push_back(_pc.computeSobolIndex({i}) / var);
  // Total
  for (const auto & i : make_range(nparam))
    val.push_back(_pc.computeSobolTotal(i) / var);
  // Second order
  for (const auto & i : make_range(nparam))
    for (const auto & j : make_range(i + 1, nparam))
      val.push_back(_pc.computeSobolIndex({i, j}) / var);
}

template <typename OutType>
void
PCSobolContext<OutType>::store(nlohmann::json & json) const
{
  SobolReporterContext<std::vector<OutType>, OutType>::storeSobol(
      json, this->_state.value(), _pc.getNumberOfParameters());
}

template class PCStatisticsContext<Real>;
template class PCSobolContext<Real>;

(contrib/moose/modules/stochastic_tools/include/surrogates/PolynomialChaos.h)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "SurrogateModel.h"
#include "PolynomialQuadrature.h"
#include "QuadratureSampler.h"

#include "Distribution.h"
#include "nlohmann/json.h"

class PolynomialChaos : public SurrogateModel
{
public:
  static InputParameters validParams();
  PolynomialChaos(const InputParameters & parameters);
  using SurrogateModel::evaluate;
  virtual Real evaluate(const std::vector<Real> & x) const override;

  /// Access number of dimensions/parameters
  std::size_t getNumberOfParameters() const { return _poly.size(); }

  /// Number of terms in expansion
  std::size_t getNumberofCoefficients() const { return _tuple.size(); }

  /// Access polynomial orders from tuple
  ////@{
  const std::vector<std::vector<unsigned int>> & getPolynomialOrders() const;
  unsigned int getPolynomialOrder(const unsigned int dim, const unsigned int i) const;
  ///@}

  /// Access computed expansion coefficients
  const std::vector<Real> & getCoefficients() const;

  /// Evaluate mean: \mu = E[u]
  virtual Real computeMean() const;

  /// Evaluate standard deviation: \sigma = sqrt(E[(u-\mu)^2])
  virtual Real computeStandardDeviation() const;

  /// Compute expectation of a certain power of the QoI: E[(u-\mu)^n]
  Real powerExpectation(const unsigned int n) const;

  /// Evaluates partial derivative of expansion: du(x)/dx_dim
  Real computeDerivative(const unsigned int dim, const std::vector<Real> & x) const;
  /**
   * Evaluates sum of partial derivative of expansion. Example:
   * computeGradient({0, 2, 3}, x) = du(x)/dx_0dx_2dx_3
   */
  Real computePartialDerivative(const std::vector<unsigned int> & dim,
                                const std::vector<Real> & x) const;

  /// Computes Sobol sensitivities S_{i_1,i_2,...,i_s}, where ind = i_1,i_2,...,i_s
  Real computeSobolIndex(const std::set<unsigned int> & ind) const;
  Real computeSobolTotal(const unsigned int dim) const;

  void store(nlohmann::json & json) const;

private:
  /// Variables calculation and for looping over the computed coefficients in parallel
  ///
  /// The various utility methods in this class require the coefficients be partitioned in parallel,
  /// but the data being partitioned is loaded from the trainer so it might not be available. Thus,
  /// the partitioning is done on demand, if needed.
  ///
  /// The methods are marked const because they do not modify the loaded data, to keep this interface
  /// the partitioning uses mutable variables.
  ///@{
  mutable dof_id_type _n_local_coeff = std::numeric_limits<dof_id_type>::max();
  mutable dof_id_type _local_coeff_begin = 0;
  mutable dof_id_type _local_coeff_end = 0;
  void linearPartitionCoefficients() const;
  ///@}

  // The following items are loaded from a SurrogateTrainer using getModelData

  /// Maximum polynomial order. The sum of 1D polynomial orders does not go above this value.
  const unsigned int & _order;

  /// Total number of parameters/dimensions
  const unsigned int & _ndim;

  /// Total number of coefficient (defined by size of _tuple)
  const std::size_t & _ncoeff;

  /// A _ndim-by-_ncoeff matrix containing the appropriate one-dimensional polynomial order
  const std::vector<std::vector<unsigned int>> & _tuple;

  /// These are the coefficients we are after in the PC expansion
  const std::vector<Real> & _coeff;

  /// The distributions used for sampling
  const std::vector<std::unique_ptr<const PolynomialQuadrature::Polynomial>> & _poly;

  friend void to_json(nlohmann::json & json, const PolynomialChaos * const & pc);
};

(contrib/moose/modules/stochastic_tools/src/surrogates/PolynomialChaos.C)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "PolynomialChaos.h"
#include "Sampler.h"
#include "CartesianProduct.h"

registerMooseObject("StochasticToolsApp", PolynomialChaos);

InputParameters
PolynomialChaos::validParams()
{
  InputParameters params = SurrogateModel::validParams();
  params.addClassDescription("Computes and evaluates polynomial chaos surrogate model.");
  return params;
}

PolynomialChaos::PolynomialChaos(const InputParameters & parameters)
  : SurrogateModel(parameters),
    _order(getModelData<unsigned int>("_order")),
    _ndim(getModelData<unsigned int>("_ndim")),
    _ncoeff(getModelData<std::size_t>("_ncoeff")),
    _tuple(getModelData<std::vector<std::vector<unsigned int>>>("_tuple")),
    _coeff(getModelData<std::vector<Real>>("_coeff")),
    _poly(
        getModelData<std::vector<std::unique_ptr<const PolynomialQuadrature::Polynomial>>>("_poly"))
{
}

Real
PolynomialChaos::evaluate(const std::vector<Real> & x) const
{
  mooseAssert(x.size() == _ndim, "Number of inputted parameters does not match PC model.");

  DenseMatrix<Real> poly_val(_ndim, _order);

  // Evaluate polynomials to avoid duplication
  for (unsigned int d = 0; d < _ndim; ++d)
    for (unsigned int i = 0; i < _order; ++i)
      poly_val(d, i) = _poly[d]->compute(i, x[d], /*normalize =*/false);

  Real val = 0;
  for (std::size_t i = 0; i < _ncoeff; ++i)
  {
    Real tmp = _coeff[i];
    for (unsigned int d = 0; d < _ndim; ++d)
      tmp *= poly_val(d, _tuple[i][d]);
    val += tmp;
  }

  return val;
}

const std::vector<std::vector<unsigned int>> &
PolynomialChaos::getPolynomialOrders() const
{
  return _tuple;
}

unsigned
PolynomialChaos::getPolynomialOrder(const unsigned int dim, const unsigned int i) const
{
  return _tuple[i][dim];
}

const std::vector<Real> &
PolynomialChaos::getCoefficients() const
{
  return _coeff;
}

Real
PolynomialChaos::computeMean() const
{
  mooseAssert(_coeff.size() > 0, "The coefficient matrix is empty.");
  return _coeff[0];
}

Real
PolynomialChaos::computeStandardDeviation() const
{
  Real var = 0;
  for (std::size_t i = 1; i < _ncoeff; ++i)
  {
    Real norm = 1.0;
    for (std::size_t d = 0; d < _ndim; ++d)
      norm *= _poly[d]->innerProduct(_tuple[i][d]);
    var += _coeff[i] * _coeff[i] * norm;
  }

  return std::sqrt(var);
}

Real
PolynomialChaos::powerExpectation(const unsigned int n) const
{
  std::vector<StochasticTools::WeightedCartesianProduct<unsigned int, Real>> order;
  order.reserve(_ndim);
  std::vector<std::vector<Real>> c_1d(n, std::vector<Real>(_coeff.begin() + 1, _coeff.end()));
  for (unsigned int d = 0; d < _ndim; ++d)
  {
    std::vector<std::vector<unsigned int>> order_1d(n, std::vector<unsigned int>(_ncoeff - 1));
    for (std::size_t i = 1; i < _ncoeff; ++i)
      for (unsigned int m = 0; m < n; ++m)
        order_1d[m][i - 1] = _tuple[i][d];
    order.push_back(StochasticTools::WeightedCartesianProduct<unsigned int, Real>(order_1d, c_1d));
  }

  dof_id_type n_local, st_local, end_local;
  MooseUtils::linearPartitionItems(
      order[0].numRows(), n_processors(), processor_id(), n_local, st_local, end_local);

  Real val = 0;
  for (dof_id_type i = st_local; i < end_local; ++i)
  {
    Real tmp = order[0].computeWeight(i);
    for (unsigned int d = 0; d < _ndim; ++d)
    {
      if (MooseUtils::absoluteFuzzyEqual(tmp, 0.0))
        break;
      std::vector<unsigned int> comb(n);
      for (unsigned int m = 0; m < n; ++m)
        comb[m] = order[d].computeValue(i, m);
      tmp *= _poly[d]->productIntegral(comb);
    }
    val += tmp;
  }

  return val;
}

Real
PolynomialChaos::computeDerivative(const unsigned int dim, const std::vector<Real> & x) const
{
  return computePartialDerivative({dim}, x);
}

Real
PolynomialChaos::computePartialDerivative(const std::vector<unsigned int> & dim,
                                          const std::vector<Real> & x) const
{
  mooseAssert(x.size() == _ndim, "Number of inputted parameters does not match PC model.");

  std::vector<unsigned int> grad(_ndim);
  for (const auto & d : dim)
  {
    mooseAssert(d < _ndim, "Specified dimension is greater than total number of parameters.");
    grad[d]++;
  }

  DenseMatrix<Real> poly_val(_ndim, _order);

  // Evaluate polynomials to avoid duplication
  for (unsigned int d = 0; d < _ndim; ++d)
    for (unsigned int i = 0; i < _order; ++i)
      poly_val(d, i) = _poly[d]->computeDerivative(i, x[d], grad[d]);

  Real val = 0;
  for (std::size_t i = 0; i < _ncoeff; ++i)
  {
    Real tmp = _coeff[i];
    for (unsigned int d = 0; d < _ndim; ++d)
      tmp *= poly_val(d, _tuple[i][d]);
    val += tmp;
  }

  return val;
}

Real
PolynomialChaos::computeSobolIndex(const std::set<unsigned int> & ind) const
{
  linearPartitionCoefficients();

  // If set is empty, compute mean
  if (ind.empty())
    return computeMean();

  // Do some sanity checks in debug
  mooseAssert(ind.size() <= _ndim, "Number of indices is greater than number of parameters.");
  mooseAssert(*ind.rbegin() < _ndim, "Maximum index provided exceeds number of parameters.");

  Real val = 0.0;
  for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
  {
    Real tmp = _coeff[i] * _coeff[i];
    for (unsigned int d = 0; d < _ndim; ++d)
    {
      if ((ind.find(d) != ind.end() && _tuple[i][d] > 0) ||
          (ind.find(d) == ind.end() && _tuple[i][d] == 0))
      {
        tmp *= _poly[d]->innerProduct(_tuple[i][d]);
      }
      else
      {
        tmp = 0.0;
        break;
      }
    }
    val += tmp;
  }

  return val;
}

Real
PolynomialChaos::computeSobolTotal(const unsigned int dim) const
{
  linearPartitionCoefficients();

  // Do some sanity checks in debug
  mooseAssert(dim < _ndim, "Requested dimension is greater than number of parameters.");

  Real val = 0.0;
  for (dof_id_type i = _local_coeff_begin; i < _local_coeff_end; ++i)
    if (_tuple[i][dim] > 0)
      val += _coeff[i] * _coeff[i] * _poly[dim]->innerProduct(_tuple[i][dim]);

  return val;
}

void
PolynomialChaos::linearPartitionCoefficients() const
{
  if (_n_local_coeff == std::numeric_limits<dof_id_type>::max())
    MooseUtils::linearPartitionItems(_ncoeff,
                                     n_processors(),
                                     processor_id(),
                                     _n_local_coeff,
                                     _local_coeff_begin,
                                     _local_coeff_end);
}

void
PolynomialChaos::store(nlohmann::json & json) const
{
  json["order"] = _order;
  json["ndim"] = getNumberOfParameters();
  json["ncoeff"] = getNumberofCoefficients();
  json["tuple"] = getPolynomialOrders();
  json["coeff"] = getCoefficients();
  for (const auto & p : _poly)
  {
    nlohmann::json jsonp;
    p->store(jsonp);
    json["poly"].push_back(jsonp);
  }
}

(contrib/moose/modules/stochastic_tools/test/tests/surrogates/poly_chaos/main_2d_mc.i)

[StochasticTools]
[]

[Distributions]
  [D_dist]
    type = Uniform
    lower_bound = 2.5
    upper_bound = 7.5
  []
  [S_dist]
    type = Uniform
    lower_bound = 2.5
    upper_bound = 7.5
  []
[]

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 100
    distributions = 'D_dist S_dist'
    execute_on = initial
  []
[]

[MultiApps]
  [quad_sub]
    type = SamplerFullSolveMultiApp
    input_files = sub.i
    sampler = sample
    mode = batch-restore
  []
[]

[Transfers]
  [quad]
    type = SamplerParameterTransfer
    to_multi_app = quad_sub
    sampler = sample
    parameters = 'Materials/diffusivity/prop_values Materials/xs/prop_values'
  []
  [data]
    type = SamplerReporterTransfer
    from_multi_app = quad_sub
    sampler = sample
    stochastic_reporter = storage
    from_reporter = avg/value
  []
[]

[Reporters]
  [storage]
    type = StochasticReporter
    outputs = none
  []
  [pc_samp]
    type = EvaluateSurrogate
    model = poly_chaos
    sampler = sample
    parallel_type = ROOT
    execute_on = final
  []
[]

[Surrogates]
  [poly_chaos]
    type = PolynomialChaos
    trainer = poly_chaos
  []
[]

[Trainers]
  [poly_chaos]
    type = PolynomialChaosTrainer
    execute_on = timestep_end
    order = 5
    distributions = 'D_dist S_dist'
    sampler = sample
    response = storage/data:avg:value
    regression_type = integration
  []
[]

[Outputs]
  [out]
    type = CSV
    execute_on = FINAL
  []
[]

Overview
Creating a SurrogateModel
Using a SurrogateModel
Getting Training Data
Example Input File Syntax
Available Objects
Available Actions