Creating StatisticsReporter to replace Statistics VPP

(refs idaholab#16584)

framework/include/reporters/ReporterContext.h

@@ -78,7 +78,7 @@ class ReporterContextBase : public libMesh::ParallelObject
   /**
    * Helper for enabling generic transfer of Reporter values
    * @param r_data The ReporterData on the app that this data is being transferred to
-   * @param r_name The name of the Report being transfered to
+   * @param r_name The name of the Report being transferred to
    *
    * @see MultiAppReporterTransfer
    */
@@ -181,9 +181,11 @@ class ReporterContext : public ReporterContextBase
   /// The state on which this context object operates
   ReporterState<T> & _state;
 
+  /// Output data to JSON, see JSONOutput
+  virtual void store(nlohmann::json & json) const override;
+
   // The following are called by the ReporterData and are not indented for external use
   virtual void copyValuesBack() override;
-  virtual void store(nlohmann::json & json) const override;
   friend class ReporterData;
 };
 
@@ -286,7 +288,10 @@ template <typename T>
 void
 ReporterContext<T>::store(nlohmann::json & json) const
 {
-  storeHelper(json, this->_state.value());
+  json["object_name"] = this->name().getObjectName();
+  json["value_name"] = this->name().getValueName();
+  json["type"] = this->type();
+  storeHelper(json["value"], this->_state.value());
 }
 
 template <typename T>

framework/include/reporters/ReporterData.h

@@ -193,7 +193,7 @@ class ReporterData
   void check() const;
 
   /**
-   * Return true if the supplied mode exists in the produced Repoter values
+   * Return true if the supplied mode exists in the produced Reporter values
    *
    * @see JSONOutput.C/h
    */
@@ -237,6 +237,13 @@ class ReporterData
    */
   const ReporterContextBase * getReporterContextBase(const ReporterName & reporter_name) const;
 
+  /**
+   * Return the ReporterProducerEnum for an existing ReporterValue
+   * @param reporter_name The name of the reporter value, which includes the object name and the
+   *                      data name.
+   */
+  const ReporterProducerEnum & getReporterMode(const ReporterName & reporter_name) const;
+
 private:
   /// For accessing the restart/recover system, which is where Reporter values are stored
   MooseApp & _app;

framework/src/outputs/JSONOutput.C

@@ -80,10 +80,7 @@ JSONOutput::outputReporters()
       ReporterName r_name(combined_name);
       const auto * context_ptr = _reporter_data.getReporterContextBase(r_name);
       auto & node = current_node["reporters"].emplace_back(); // non-accidental insertion
-      node["object_name"] = context_ptr->name().getObjectName();
-      node["value_name"] = context_ptr->name().getValueName();
-      node["type"] = context_ptr->type();
-      context_ptr->store(node["value"]);
+      context_ptr->store(node);
     }
   }
 }

framework/src/reporters/ReporterData.C

@@ -140,3 +140,12 @@ ReporterData::addConsumerMode(ReporterMode mode, const std::string & object_name
     mooseError("Unable to locate Reporter with name:", object_name);
   const_cast<ReporterContextBase *>(ptr)->addConsumerMode(mode, object_name);
 }
+
+const ReporterProducerEnum &
+ReporterData::getReporterMode(const ReporterName & reporter_name) const
+{
+  const ReporterContextBase * ptr = getReporterContextBase(reporter_name);
+  if (ptr == nullptr)
+    mooseError("Unable to locate Reporter with name:", reporter_name);
+  return ptr->getProducerModeEnum();
+}

modules/stochastic_tools/include/reporters/StatisticsReporter.h

@@ -0,0 +1,122 @@
+//* 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 "Calculators.h"
+#include "BootstrapCalculators.h"
+
+#include "nlohmann/json.h"
+
+/**
+ * ReporterContext that utilizes a Calculator object to compute its value and confidence levels
+ */
+template <typename T>
+class ReporterStatisticsContext : public ReporterContext<T>
+{
+public:
+  ReporterStatisticsContext(const libMesh::ParallelObject & other,
+                            ReporterState<T> & state,
+                            const std::vector<T> & data,
+                            const ReporterProducerEnum & mode,
+                            const MooseEnumItem & stat,
+                            const StochasticTools::BootstrapCalculator * ci_calc);
+  virtual void finalize() override;
+  virtual void store(nlohmann::json & json) const override;
+
+private:
+  /// Data used for the statistic calculation
+  const std::vector<T> & _data;
+
+  /// Mode in which the above data was produced
+  const ReporterProducerEnum & _data_mode;
+
+  /// Storage for the Calculator object for the desired stat, this is created in constructor
+  std::unique_ptr<const StochasticTools::Calculator> _calc_ptr;
+
+  /// Storage for the BootstrapCalculator for the desired confidence interval calculations (optional)
+  const StochasticTools::BootstrapCalculator * _ci_calc_ptr;
+
+  /// The results
+  std::vector<T> _ci_results;
+};
+
+template <typename T>
+ReporterStatisticsContext<T>::ReporterStatisticsContext(
+    const libMesh::ParallelObject & other,
+    ReporterState<T> & state,
+    const std::vector<T> & data,
+    const ReporterProducerEnum & mode,
+    const MooseEnumItem & stat,
+    const StochasticTools::BootstrapCalculator * ci_calc)
+  : ReporterContext<T>(other, state),
+    _data(data),
+    _data_mode(mode),
+    _calc_ptr(StochasticTools::makeCalculator(stat, other)),
+    _ci_calc_ptr(ci_calc)
+{
+}
+
+template <typename T>
+void
+ReporterStatisticsContext<T>::finalize()
+{
+  this->_state.value() = _calc_ptr->compute(_data, _data_mode == REPORTER_MODE_DISTRIBUTED);
+  ReporterContext<T>::finalize();
+
+  if (_ci_calc_ptr)
+    _ci_results = _ci_calc_ptr->compute(_data, *_calc_ptr, _data_mode == REPORTER_MODE_DISTRIBUTED);
+}
+
+template <typename T>
+void
+ReporterStatisticsContext<T>::store(nlohmann::json & json) const
+{
+  ReporterContext<T>::store(json);
+  json["stat"] = _calc_ptr->name();
+  if (_ci_calc_ptr)
+    json["confidence_intervals"] = {{"method", _ci_calc_ptr->name()},
+                                    {"values", _ci_results},
+                                    {"levels", _ci_calc_ptr->levels()},
+                                    {"replicates", _ci_calc_ptr->replicates()},
+                                    {"seed", _ci_calc_ptr->seed()}};
+}
+
+/**
+ * Compute several metrics for supplied data.
+ *
+ * This class uses Calculator objects defined in StatisticsReporter.h and is setup such that if a
+ * new calculation is needed it can be added in StatisticsReporter.h without modification of this
+ * object.
+ */
+class StatisticsReporter : public GeneralReporter
+{
+public:
+  static InputParameters validParams();
+  StatisticsReporter(const InputParameters & parameters);
+
+  /// Not used; all operations are  wrapped in the ReporterStatisticsContext
+  virtual void execute() final{};
+  virtual void initialize() final{};
+  virtual void finalize() final{};
+
+protected:
+  /// Confidence level calculator, this is shared by all reporters that are declared.
+  std::unique_ptr<const StochasticTools::BootstrapCalculator> _ci_calculator = nullptr;
+
+private:
+  /**
+   * Helper function for converting confidence levels given in (0, 0.5] into levels in (0, 1).
+   * For example, levels_in = {0.05, 0.1, 0.5} converts to {0.05 0.1, 0.5, 0.9, 0.95}.
+   *
+   * This also performs error checking on the supplied "ci_levels".
+   */
+  std::vector<Real> computeLevels(const std::vector<Real> & levels_in) const;
+};

modules/stochastic_tools/include/utils/BootstrapCalculators.h

@@ -49,6 +49,7 @@ class BootstrapCalculator : public libMesh::ParallelObject
 {
 public:
   BootstrapCalculator(const libMesh::ParallelObject & other,
+                      const std::string & method,
                       const std::vector<Real> & levels,
                       unsigned int replicates,
                       unsigned int seed);
@@ -64,6 +65,14 @@ class BootstrapCalculator : public libMesh::ParallelObject
                                     const Calculator & calc,
                                     const bool is_distributed) const = 0;
 
+  ///@{
+  /// Return the input items (see ReporterStatisticsContext)
+  const std::string & name() const;
+  const std::vector<Real> & levels() const;
+  unsigned int replicates() const;
+  unsigned int seed() const;
+  ///@}
+
 protected:
   // Compute Bootstrap estimates of a statistic
   std::vector<Real>
@@ -72,6 +81,9 @@ class BootstrapCalculator : public libMesh::ParallelObject
   // Randomly shuffle a vector of data
   std::vector<Real> shuffle(const std::vector<Real> &, MooseRandom &, const bool) const;
 
+  // Calculation method
+  const std::string _method;
+
   // Confidence levels to compute in range (0, 1)
   const std::vector<Real> _levels;
 

modules/stochastic_tools/include/utils/Calculators.h

@@ -46,9 +46,13 @@ MultiMooseEnum makeCalculatorEnum();
 class Calculator : public libMesh::ParallelObject
 {
 public:
-  Calculator(const libMesh::ParallelObject &);
+  Calculator(const libMesh::ParallelObject &, const std::string & stat);
   virtual ~Calculator() = default;
+  const std::string & name() const;
   virtual Real compute(const std::vector<Real> &, bool) const = 0;
+
+private:
+  const std::string _stat;
 };
 
 class Mean : public Calculator
@@ -83,6 +87,7 @@ class StdDev : public Calculator
 {
 public:
   StdDev(const libMesh::ParallelObject &);
+  StdDev(const libMesh::ParallelObject &, const std::string & stat);
   virtual Real compute(const std::vector<Real> &, bool) const override;
 };
 

modules/stochastic_tools/src/reporters/StatisticsReporter.C

@@ -0,0 +1,125 @@
+//* 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 "StatisticsReporter.h"
+
+// MOOSE includes
+#include "MooseVariable.h"
+#include "ThreadedElementLoopBase.h"
+#include "ThreadedNodeLoop.h"
+
+#include "libmesh/quadrature.h"
+
+#include <numeric>
+
+registerMooseObject("StochasticToolsApp", StatisticsReporter);
+
+InputParameters
+StatisticsReporter::validParams()
+{
+  InputParameters params = GeneralReporter::validParams();
+  params.addClassDescription(
+      "Compute statistical values of a given VectorPostprocessor objects and vectors.");
+
+  params.addRequiredParam<std::vector<VectorPostprocessorName>>(
+      "vectorpostprocessors",
+      "List of VectorPostprocessor(s) to utilized for statistic computations.");
+
+  MultiMooseEnum stats = StochasticTools::makeCalculatorEnum();
+  params.addRequiredParam<MultiMooseEnum>(
+      "compute",
+      stats,
+      "The statistic(s) to compute for each of the supplied vector postprocessors.");
+
+  // Confidence Levels
+  MooseEnum ci = StochasticTools::makeBootstrapCalculatorEnum();
+  params.addParam<MooseEnum>(
+      "ci_method", ci, "The method to use for computing confidence level intervals.");
+
+  params.addParam<std::vector<Real>>(
+      "ci_levels", "A vector of confidence levels to consider, values must be in (0, 0.5].");
+  params.addParam<unsigned int>(
+      "ci_replicates",
+      10000,
+      "The number of replicates to use when computing confidence level intervals.");
+  params.addParam<unsigned int>("ci_seed",
+                                1,
+                                "The random number generator seed used for creating replicates "
+                                "while computing confidence level intervals.");
+  return params;
+}
+
+StatisticsReporter::StatisticsReporter(const InputParameters & parameters)
+  : GeneralReporter(parameters)
+{
+  // Statistics to be computed
+  const auto & compute_stats = getParam<MultiMooseEnum>("compute");
+
+  // Bootstrap CI
+  std::unique_ptr<const StochasticTools::BootstrapCalculator> ci_calculator = nullptr;
+  const MooseEnum & ci_method = getParam<MooseEnum>("ci_method");
+  if (ci_method.isValid())
+  {
+    std::vector<Real> ci_levels = computeLevels(getParam<std::vector<Real>>("ci_levels"));
+    unsigned int replicates = getParam<unsigned int>("ci_replicates");
+    unsigned int seed = getParam<unsigned int>("ci_seed");
+    _ci_calculator =
+        StochasticTools::makeBootstrapCalculator(ci_method, *this, ci_levels, replicates, seed);
+  }
+
+  // Stats for VPP
+  const auto & vpp_names = getParam<std::vector<VectorPostprocessorName>>("vectorpostprocessors");
+  for (const auto & vpp_name : vpp_names)
+  {
+    const VectorPostprocessor & vpp_object =
+        _fe_problem.getVectorPostprocessorObjectByName(vpp_name);
+    const std::set<std::string> & vpp_vectors = vpp_object.getVectorNames();
+    for (const auto & vec_name : vpp_vectors)
+    {
+      ReporterName r_name(vpp_name, vec_name);
+      const auto & data = getReporterValueByName<std::vector<Real>>(r_name);
+      const auto & mode = _fe_problem.getReporterData().getReporterMode(r_name);
+      for (const auto & item : compute_stats)
+      {
+        const std::string s_name = vpp_name + "::" + vec_name + "_" + item.name();
+        declareValueByName<Real, ReporterStatisticsContext>(
+            s_name, REPORTER_MODE_ROOT, data, mode, item, _ci_calculator.get());
+      }
+    }
+  }
+}
+
+std::vector<Real>
+StatisticsReporter::computeLevels(const std::vector<Real> & levels_in) const
+{
+  if (levels_in.empty())
+    paramError("ci_levels",
+               "If the 'ci_method' parameter is supplied then the 'ci_levels' must also be "
+               "supplied with values in (0, 0.5].");
+
+  else if (*std::min_element(levels_in.begin(), levels_in.end()) <= 0)
+    paramError("ci_levels", "The supplied levels must be greater than zero.");
+
+  else if (*std::max_element(levels_in.begin(), levels_in.end()) > 0.5)
+    paramError("ci_levels", "The supplied levels must be less than or equal to 0.5");
+
+  std::list<Real> levels_out;
+  for (auto it = levels_in.rbegin(); it != levels_in.rend(); ++it)
+  {
+    if (*it == 0.5)
+      levels_out.push_back(*it);
+
+    else
+    {
+      levels_out.push_front(*it);
+      levels_out.push_back(1 - *it);
+    }
+  }
+  return std::vector<Real>(levels_out.begin(), levels_out.end());
+}