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 "PolynomialChaosTrainer.h"
11 : #include "Sampler.h"
12 : #include "CartesianProduct.h"
13 :
14 : registerMooseObject("StochasticToolsApp", PolynomialChaosTrainer);
15 :
16 : InputParameters
17 1020 : PolynomialChaosTrainer::validParams()
18 : {
19 1020 : InputParameters params = SurrogateTrainer::validParams();
20 1020 : params.addClassDescription("Computes and evaluates polynomial chaos surrogate model.");
21 2040 : params.addRequiredParam<unsigned int>("order", "Maximum polynomial order.");
22 2040 : params.addRequiredParam<std::vector<DistributionName>>(
23 : "distributions", "Names of the distributions samples were taken from.");
24 2040 : MooseEnum rtype("integration ols auto", "auto");
25 2040 : params.addParam<MooseEnum>(
26 : "regression_type",
27 : rtype,
28 : "The type of regression to perform for finding polynomial coefficents.");
29 2040 : params.addParam<Real>("penalty", 0.0, "Ridge regularization penalty factor for OLS regression.");
30 :
31 1020 : params.suppressParameter<MooseEnum>("response_type");
32 1020 : return params;
33 1020 : }
34 :
35 510 : PolynomialChaosTrainer::PolynomialChaosTrainer(const InputParameters & parameters)
36 : : SurrogateTrainer(parameters),
37 510 : _predictor_row(getPredictorData()),
38 1530 : _order(declareModelData<unsigned int>("_order", getParam<unsigned int>("order"))),
39 1020 : _ndim(declareModelData<unsigned int>("_ndim", _sampler.getNumberOfCols())),
40 1020 : _tuple(declareModelData<std::vector<std::vector<unsigned int>>>(
41 1020 : "_tuple", StochasticTools::MultiDimPolynomialGenerator::generateTuple(_ndim, _order))),
42 1020 : _ncoeff(declareModelData<std::size_t>("_ncoeff", _tuple.size())),
43 1020 : _coeff(declareModelData<std::vector<Real>>("_coeff")),
44 1020 : _poly(declareModelData<std::vector<std::unique_ptr<const PolynomialQuadrature::Polynomial>>>(
45 : "_poly")),
46 1020 : _ridge_penalty(getParam<Real>("penalty")),
47 1020 : _quad_sampler(dynamic_cast<QuadratureSampler *>(&_sampler))
48 : {
49 : // Check if number of distributions is correct
50 510 : if (_ndim != _sampler.getNumberOfCols())
51 0 : paramError("distributions",
52 : "Sampler number of columns does not match number of inputted distributions.");
53 :
54 1020 : auto rtype_enum = getParam<MooseEnum>("regression_type");
55 510 : if (rtype_enum == "auto")
56 430 : _rtype = _quad_sampler ? 0 : 1;
57 : else
58 80 : _rtype = rtype_enum == "integration" ? 0 : 1;
59 :
60 510 : if (_rtype == 0 && _quad_sampler &&
61 350 : (!_pvals.empty() || _pcols.size() != _sampler.getNumberOfCols()))
62 0 : paramError("sampler",
63 : "QuadratureSampler must use all Sampler columns for training, and cannot be"
64 : " used with other Reporters - otherwise, quadrature integration does not work.");
65 510 : if (_rtype == 0 && _ridge_penalty != 0.0)
66 0 : paramError("penalty",
67 : "Ridge regularization penalty is only relevant if 'regression_type = ols'.");
68 :
69 : // Make polynomials
70 2744 : for (const auto & nm : getParam<std::vector<DistributionName>>("distributions"))
71 3448 : _poly.push_back(PolynomialQuadrature::makePolynomial(&getDistributionByName(nm)));
72 :
73 : // Create calculators for standardization
74 510 : if (_rtype == 1)
75 : {
76 80 : _calculators.resize(_ncoeff * 3);
77 1056 : for (const auto & term : make_range(_ncoeff))
78 : {
79 1952 : _calculators[3 * term] = StochasticTools::makeCalculator(MooseEnumItem("mean"), *this);
80 1952 : _calculators[3 * term + 1] = StochasticTools::makeCalculator(MooseEnumItem("stddev"), *this);
81 2928 : _calculators[3 * term + 2] = StochasticTools::makeCalculator(MooseEnumItem("sum"), *this);
82 : }
83 : }
84 510 : }
85 :
86 : void
87 606 : PolynomialChaosTrainer::preTrain()
88 : {
89 606 : _coeff.assign(_ncoeff, 0.0);
90 :
91 606 : if (_rtype == 1)
92 : {
93 176 : if (getCurrentSampleSize() < _ncoeff)
94 0 : paramError("order",
95 : "Number of data points (",
96 : getCurrentSampleSize(),
97 : ") must be greater than the number of terms in the polynomial (",
98 : _ncoeff,
99 : ").");
100 176 : _matrix.resize(_ncoeff, _ncoeff);
101 176 : _rhs.resize(_ncoeff);
102 3392 : for (auto & calc : _calculators)
103 : calc->initializeCalculator();
104 176 : _r_sum = 0.0;
105 : }
106 606 : }
107 :
108 : void
109 180980 : PolynomialChaosTrainer::train()
110 : {
111 : // Evaluate polynomials to avoid duplication
112 180980 : DenseMatrix<Real> poly_val(_ndim, _order);
113 1231260 : for (unsigned int d = 0; d < _ndim; ++d)
114 5321880 : for (unsigned int i = 0; i < _order; ++i)
115 4271600 : poly_val(d, i) = _poly[d]->compute(i, _predictor_row[d], /*normalize=*/_rtype == 0);
116 :
117 : // Evaluate multi-dimensional polynomials
118 180980 : std::vector<Real> basis(_ncoeff, 1.0);
119 17013700 : for (const auto i : make_range(_ncoeff))
120 112225920 : for (const auto d : make_range(_ndim))
121 95393200 : basis[i] *= poly_val(d, _tuple[i][d]);
122 :
123 : // For integration
124 180980 : if (_rtype == 0)
125 : {
126 176580 : const Real fact = (*_rval) * (_quad_sampler ? _quad_sampler->getQuadratureWeight(_row) : 1.0);
127 16948900 : for (const auto i : make_range(_ncoeff))
128 16772320 : _coeff[i] += fact * basis[i];
129 : }
130 : // For least-squares
131 : else
132 : {
133 : // Loop over coefficients
134 64800 : for (const auto i : make_range(_ncoeff))
135 : {
136 : // Matrix is symmetric, so we'll add the upper diagonal later
137 540800 : for (const auto j : make_range(i + 1))
138 480400 : _matrix(i, j) += basis[i] * basis[j];
139 60400 : _rhs(i) += basis[i] * (*_rval);
140 :
141 241600 : for (unsigned int c = i * 3; c < (i + 1) * 3; ++c)
142 : _calculators[c]->updateCalculator(basis[i]);
143 : }
144 4400 : _r_sum += (*_rval);
145 :
146 4400 : if (_ridge_penalty != 0.0)
147 32000 : for (const auto i : make_range(_ncoeff))
148 30000 : _matrix(i, i) += _ridge_penalty;
149 : }
150 180980 : }
151 :
152 : void
153 606 : PolynomialChaosTrainer::postTrain()
154 : {
155 606 : if (_rtype == 0)
156 : {
157 430 : gatherSum(_coeff);
158 430 : if (!_quad_sampler)
159 46080 : for (std::size_t i = 0; i < _ncoeff; ++i)
160 46000 : _coeff[i] /= getCurrentSampleSize();
161 : }
162 : else
163 : {
164 : gatherSum(_matrix.get_values());
165 : gatherSum(_rhs.get_values());
166 3392 : for (auto & calc : _calculators)
167 3216 : calc->finalizeCalculator(true);
168 176 : gatherSum(_r_sum);
169 :
170 176 : std::vector<Real> mu(_ncoeff);
171 176 : std::vector<Real> sig(_ncoeff);
172 176 : std::vector<Real> sum_pf(_ncoeff);
173 1248 : for (const auto i : make_range(_ncoeff))
174 : {
175 1072 : mu[i] = i > 0 ? _calculators[3 * i]->getValue() : 0.0;
176 1072 : sig[i] = i > 0 ? _calculators[3 * i + 1]->getValue() : 1.0;
177 1072 : sum_pf[i] = _calculators[3 * i + 2]->getValue();
178 : }
179 :
180 176 : const Real n = getCurrentSampleSize();
181 1248 : for (const auto i : make_range(_ncoeff))
182 : {
183 8864 : for (const auto j : make_range(i + 1))
184 : {
185 7792 : _matrix(i, j) -= (mu[j] * sum_pf[i] + mu[i] * sum_pf[j]);
186 7792 : _matrix(i, j) += n * mu[i] * mu[j];
187 7792 : _matrix(i, j) /= (sig[i] * sig[j]);
188 7792 : _matrix(j, i) = _matrix(i, j);
189 : }
190 1072 : _rhs(i) = (_rhs(i) - mu[i] * _r_sum) / sig[i];
191 : }
192 :
193 176 : DenseVector<Real> sol;
194 176 : _matrix.lu_solve(_rhs, sol);
195 176 : _coeff = sol.get_values();
196 :
197 1072 : for (unsigned int i = 1; i < _ncoeff; ++i)
198 : {
199 896 : _coeff[i] /= sig[i];
200 896 : _coeff[0] -= _coeff[i] * mu[i];
201 : }
202 : }
203 606 : }
|
