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 "GaussianProcessSurrogate.h"
      11             : #include "Sampler.h"
      12             : 
      13             : #include "CovarianceFunctionBase.h"
      14             : 
      15             : registerMooseObject("StochasticToolsApp", GaussianProcessSurrogate);
      16             : 
      17             : InputParameters
      18         744 : GaussianProcessSurrogate::validParams()
      19             : {
      20         744 :   InputParameters params = SurrogateModel::validParams();
      21         744 :   params.addClassDescription("Computes and evaluates Gaussian Process surrogate model.");
      22         744 :   return params;
      23           0 : }
      24             : 
      25         370 : GaussianProcessSurrogate::GaussianProcessSurrogate(const InputParameters & parameters)
      26             :   : SurrogateModel(parameters),
      27             :     CovarianceInterface(parameters),
      28         370 :     _gp(declareModelData<StochasticTools::GaussianProcess>("_gp")),
      29        1110 :     _training_params(getModelData<RealEigenMatrix>("_training_params"))
      30             : {
      31         370 : }
      32             : 
      33             : void
      34          32 : GaussianProcessSurrogate::setupCovariance(UserObjectName covar_name)
      35             : {
      36          32 :   if (_gp.getCovarFunctionPtr() != nullptr)
      37           0 :     ::mooseError("Attempting to redefine covariance function using setupCovariance.");
      38          32 :   _gp.linkCovarianceFunction(getCovarianceFunctionByName(covar_name));
      39          32 : }
      40             : 
      41             : Real
      42         300 : GaussianProcessSurrogate::evaluate(const std::vector<Real> & x) const
      43             : {
      44             :   // Overlaod for evaluate to maintain general compatibility. Only returns mean
      45         300 :   Real dummy = 0;
      46         300 :   return this->evaluate(x, dummy);
      47             : }
      48             : 
      49             : Real
      50      115814 : GaussianProcessSurrogate::evaluate(const std::vector<Real> & x, Real & std_dev) const
      51             : {
      52             :   std::vector<Real> y;
      53             :   std::vector<Real> std;
      54      115814 :   this->evaluate(x, y, std);
      55      115814 :   std_dev = std[0];
      56      231628 :   return y[0];
      57             : }
      58             : 
      59             : void
      60           0 : GaussianProcessSurrogate::evaluate(const std::vector<Real> & x, std::vector<Real> & y) const
      61             : {
      62             :   // Overlaod for evaluate to maintain general compatibility. Only returns mean
      63             :   std::vector<Real> std_dummy;
      64           0 :   this->evaluate(x, y, std_dummy);
      65           0 : }
      66             : 
      67             : void
      68      115964 : GaussianProcessSurrogate::evaluate(const std::vector<Real> & x,
      69             :                                    std::vector<Real> & y,
      70             :                                    std::vector<Real> & std) const
      71             : {
      72      115964 :   const unsigned int n_dims = _training_params.cols();
      73             : 
      74             :   mooseAssert(x.size() == n_dims,
      75             :               "Number of parameters provided for evaluation does not match number of parameters "
      76             :               "used for training.");
      77      115964 :   const unsigned int n_outputs = _gp.getCovarFunction().numOutputs();
      78             : 
      79      115964 :   y = std::vector<Real>(n_outputs, 0.0);
      80      115964 :   std = std::vector<Real>(n_outputs, 0.0);
      81             : 
      82      115964 :   RealEigenMatrix test_points(1, n_dims);
      83      348410 :   for (unsigned int ii = 0; ii < n_dims; ++ii)
      84      232446 :     test_points(0, ii) = x[ii];
      85             : 
      86      115964 :   _gp.getParamStandardizer().getStandardized(test_points);
      87             : 
      88      115964 :   RealEigenMatrix K_train_test(_training_params.rows() * n_outputs, n_outputs);
      89             : 
      90      115964 :   _gp.getCovarFunction().computeCovarianceMatrix(
      91             :       K_train_test, _training_params, test_points, false);
      92      115964 :   RealEigenMatrix K_test(n_outputs, n_outputs);
      93      115964 :   _gp.getCovarFunction().computeCovarianceMatrix(K_test, test_points, test_points, true);
      94             : 
      95             :   // Compute the predicted mean value (centered)
      96      115964 :   RealEigenMatrix pred_value = (K_train_test.transpose() * _gp.getKResultsSolve()).transpose();
      97             :   // De-center/scale the value and store for return
      98      115964 :   _gp.getDataStandardizer().getDestandardized(pred_value);
      99             : 
     100             :   RealEigenMatrix pred_var =
     101      231928 :       K_test - (K_train_test.transpose() * _gp.getKCholeskyDecomp().solve(K_train_test));
     102             : 
     103             :   // Vairance computed, take sqrt for standard deviation, scale up by training data std and store
     104           0 :   RealEigenMatrix std_dev_mat = pred_var.array().sqrt();
     105      115964 :   _gp.getDataStandardizer().getDescaled(std_dev_mat);
     106             : 
     107      232078 :   for (const auto output_i : make_range(n_outputs))
     108             :   {
     109      116114 :     y[output_i] = pred_value(0, output_i);
     110      116114 :     std[output_i] = std_dev_mat(output_i, output_i);
     111             :   }
     112      115964 : }