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 "PolynomialRegressionTrainer.h"
11 : #include "Sampler.h"
12 :
13 : registerMooseObject("StochasticToolsApp", PolynomialRegressionTrainer);
14 :
15 : InputParameters
16 508 : PolynomialRegressionTrainer::validParams()
17 : {
18 508 : InputParameters params = SurrogateTrainer::validParams();
19 :
20 508 : params.addClassDescription("Computes coefficients for polynomial regession model.");
21 :
22 1016 : MooseEnum rtype("ols=0 ridge=1");
23 1016 : params.addRequiredParam<MooseEnum>(
24 : "regression_type", rtype, "The type of regression to perform.");
25 1016 : params.addRequiredParam<unsigned int>("max_degree",
26 : "Maximum polynomial degree to use for the regression.");
27 1016 : params.addParam<Real>("penalty", 0.0, "Penalty for Ridge regularization.");
28 :
29 508 : return params;
30 508 : }
31 :
32 256 : PolynomialRegressionTrainer::PolynomialRegressionTrainer(const InputParameters & parameters)
33 : : SurrogateTrainer(parameters), // SurrogateTrainer(parameters),
34 252 : _predictor_row(getPredictorData()),
35 504 : _regression_type(getParam<MooseEnum>("regression_type")),
36 504 : _penalty(getParam<Real>("penalty")),
37 504 : _coeff(declareModelData<std::vector<std::vector<Real>>>("_coeff")),
38 252 : _max_degree(
39 756 : declareModelData<unsigned int>("_max_degree", getParam<unsigned int>("max_degree"))),
40 504 : _power_matrix(declareModelData<std::vector<std::vector<unsigned int>>>(
41 : "_power_matrix",
42 504 : StochasticTools::MultiDimPolynomialGenerator::generateTuple(_n_dims, _max_degree + 1))),
43 252 : _n_poly_terms(_power_matrix.size()),
44 252 : _matrix(_n_poly_terms, _n_poly_terms),
45 252 : _rhs(1, DenseVector<Real>(_n_poly_terms, 0.0)),
46 508 : _r_sum(1, 0.0)
47 : {
48 252 : _coeff.resize(_n_poly_terms);
49 :
50 : // Throwing a warning if the penalty parameter is set with OLS regression
51 252 : if (_regression_type == "ols" && _penalty != 0)
52 0 : mooseWarning("Penalty parameter is not used for OLS regression, found penalty=", _penalty);
53 :
54 : // Check if we have enough data points to solve the problem
55 252 : if (_sampler.getNumberOfRows() < _n_poly_terms)
56 0 : mooseError("Number of data points must be greater than the number of terms in the polynomial.");
57 :
58 : // Creating calculators needed for feature standardization.
59 252 : _calculators.resize(_n_poly_terms * 3);
60 1804 : for (const auto & term : make_range(_n_poly_terms))
61 : {
62 3104 : _calculators[3 * term] = StochasticTools::makeCalculator(MooseEnumItem("mean"), *this);
63 3104 : _calculators[3 * term + 1] = StochasticTools::makeCalculator(MooseEnumItem("stddev"), *this);
64 4656 : _calculators[3 * term + 2] = StochasticTools::makeCalculator(MooseEnumItem("sum"), *this);
65 : }
66 252 : }
67 :
68 : void
69 536 : PolynomialRegressionTrainer::preTrain()
70 : {
71 : _matrix.zero();
72 1936 : for (unsigned int r = 0; r < _rhs.size(); ++r)
73 : _rhs[r].zero();
74 :
75 : /// Init calculators.
76 8528 : for (auto & calc : _calculators)
77 : calc->initializeCalculator();
78 :
79 536 : _r_sum.assign(_rhs.size(), 0.0);
80 536 : }
81 :
82 : void
83 16300 : PolynomialRegressionTrainer::train()
84 : {
85 : // Caching the different powers of data to accelerate the assembly of the
86 : // system
87 16300 : DenseMatrix<Real> data_pow(_n_dims, _max_degree + 1);
88 57000 : for (unsigned int d = 0; d < _n_dims; ++d)
89 193220 : for (unsigned int i = 0; i <= _max_degree; ++i)
90 152520 : data_pow(d, i) = pow(_predictor_row[d], i);
91 :
92 : // Emplace new values if necessary
93 16300 : if (_rvecval)
94 914 : for (unsigned int r = _rhs.size(); r < _rvecval->size(); ++r)
95 : {
96 254 : _rhs.emplace_back(_n_poly_terms, 0.0);
97 254 : _r_sum.emplace_back(0.0);
98 : }
99 :
100 247880 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
101 : {
102 231580 : Real i_value(1.0);
103 899180 : for (unsigned int ii = 0; ii < _n_dims; ++ii)
104 667600 : i_value *= data_pow(ii, _power_matrix[i][ii]);
105 :
106 4446840 : for (unsigned int j = 0; j < _n_poly_terms; ++j)
107 : {
108 : Real j_value(1.0);
109 16833740 : for (unsigned int jj = 0; jj < _n_dims; ++jj)
110 12618480 : j_value *= data_pow(jj, _power_matrix[j][jj]);
111 :
112 4215260 : _matrix(i, j) += i_value * j_value;
113 : }
114 :
115 : // Update calculators.
116 926320 : for (unsigned int c = i * 3; c < (i + 1) * 3; ++c)
117 694740 : _calculators[c]->updateCalculator(i_value);
118 :
119 231580 : if (_rval)
120 229300 : _rhs[0](i) += i_value * (*_rval);
121 2280 : else if (_rvecval)
122 21240 : for (unsigned int r = 0; r < _rvecval->size(); ++r)
123 18960 : _rhs[r](i) += i_value * (*_rvecval)[r];
124 : }
125 :
126 16300 : if (_rval)
127 15640 : _r_sum[0] += (*_rval);
128 660 : else if (_rvecval)
129 5980 : for (unsigned int r = 0; r < _rvecval->size(); ++r)
130 5320 : _r_sum[r] += (*_rvecval)[r];
131 :
132 : // Adding penalty term for Ridge regularization
133 16300 : if (_regression_type == "ridge")
134 0 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
135 0 : _matrix(i, i) += _penalty;
136 16300 : }
137 :
138 : void
139 536 : PolynomialRegressionTrainer::postTrain()
140 : {
141 : // Make sure _rhs are all the same size
142 536 : unsigned int nrval = _rhs.size();
143 : gatherMax(nrval);
144 610 : for (unsigned int r = _rhs.size(); r < nrval; ++r)
145 : {
146 74 : _rhs.emplace_back(_n_poly_terms, 0.0);
147 74 : _r_sum.emplace_back(0.0);
148 : }
149 :
150 : // Gather regression data
151 : gatherSum(_matrix.get_values());
152 2264 : for (auto & it : _rhs)
153 : gatherSum(it.get_values());
154 :
155 : // Gather response sums.
156 536 : gatherSum(_r_sum);
157 :
158 : // Finalize calculators.
159 8528 : for (auto & calc : _calculators)
160 7992 : calc->finalizeCalculator(true);
161 :
162 536 : std::vector<Real> mu(_n_poly_terms);
163 536 : std::vector<Real> sig(_n_poly_terms);
164 536 : std::vector<Real> sum_pf(_n_poly_terms);
165 :
166 3200 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
167 : {
168 : // To handle intercept, use mu = 0, sig = 1.
169 2664 : mu[i] = i > 0 ? _calculators[3 * i]->getValue() : 0.0;
170 2664 : sig[i] = i > 0 ? _calculators[3 * i + 1]->getValue() : 1.0;
171 2664 : sum_pf[i] = _calculators[3 * i + 2]->getValue();
172 : }
173 :
174 : // Transform _matrix and _rhs to match standardized features.
175 536 : const Real n = sum_pf[0]; // Sum of intercept term = n.
176 3200 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
177 : {
178 21920 : for (unsigned int j = 0; j < _n_poly_terms; ++j)
179 : {
180 19256 : _matrix(i, j) -= (mu[j] * sum_pf[i] + mu[i] * sum_pf[j]);
181 19256 : _matrix(i, j) += n * mu[i] * mu[j];
182 19256 : _matrix(i, j) /= (sig[i] * sig[j]);
183 : }
184 9336 : for (unsigned int r = 0; r < _rhs.size(); ++r)
185 6672 : _rhs[r](i) = (_rhs[r](i) - mu[i] * _r_sum[r]) / sig[i];
186 : }
187 :
188 536 : DenseVector<Real> sol;
189 :
190 536 : _coeff.resize(_rhs.size());
191 2264 : for (unsigned int r = 0; r < _rhs.size(); ++r)
192 : {
193 1728 : _matrix.lu_solve(_rhs[r], sol);
194 1728 : _coeff[r] = sol.get_values();
195 :
196 : // Transform coefficients to match unstandardized features.
197 6672 : for (unsigned int i = 1; i < _n_poly_terms; ++i)
198 : {
199 4944 : _coeff[r][i] /= sig[i];
200 4944 : _coeff[r][0] -= _coeff[r][i] * mu[i];
201 : }
202 : }
203 536 : }
|
