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 244 : PolynomialRegressionTrainer::validParams()
17 : {
18 244 : InputParameters params = SurrogateTrainer::validParams();
19 :
20 244 : params.addClassDescription("Computes coefficients for polynomial regession model.");
21 :
22 488 : MooseEnum rtype("ols=0 ridge=1");
23 488 : params.addRequiredParam<MooseEnum>(
24 : "regression_type", rtype, "The type of regression to perform.");
25 488 : params.addRequiredParam<unsigned int>("max_degree",
26 : "Maximum polynomial degree to use for the regression.");
27 488 : params.addParam<Real>("penalty", 0.0, "Penalty for Ridge regularization.");
28 :
29 244 : return params;
30 244 : }
31 :
32 123 : PolynomialRegressionTrainer::PolynomialRegressionTrainer(const InputParameters & parameters)
33 : : SurrogateTrainer(parameters), // SurrogateTrainer(parameters),
34 121 : _predictor_row(getPredictorData()),
35 242 : _regression_type(getParam<MooseEnum>("regression_type")),
36 242 : _penalty(getParam<Real>("penalty")),
37 242 : _coeff(declareModelData<std::vector<std::vector<Real>>>("_coeff")),
38 121 : _max_degree(
39 363 : declareModelData<unsigned int>("_max_degree", getParam<unsigned int>("max_degree"))),
40 242 : _power_matrix(declareModelData<std::vector<std::vector<unsigned int>>>(
41 : "_power_matrix",
42 242 : StochasticTools::MultiDimPolynomialGenerator::generateTuple(_n_dims, _max_degree + 1))),
43 121 : _n_poly_terms(_power_matrix.size()),
44 121 : _matrix(_n_poly_terms, _n_poly_terms),
45 242 : _rhs(1, DenseVector<Real>(_n_poly_terms, 0.0)),
46 244 : _r_sum(1, 0.0)
47 : {
48 121 : _coeff.resize(_n_poly_terms);
49 :
50 : // Throwing a warning if the penalty parameter is set with OLS regression
51 121 : 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 121 : 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 121 : _calculators.resize(_n_poly_terms * 3);
60 850 : for (const auto & term : make_range(_n_poly_terms))
61 : {
62 1458 : _calculators[3 * term] = StochasticTools::makeCalculator(MooseEnumItem("mean"), *this);
63 1458 : _calculators[3 * term + 1] = StochasticTools::makeCalculator(MooseEnumItem("stddev"), *this);
64 2187 : _calculators[3 * term + 2] = StochasticTools::makeCalculator(MooseEnumItem("sum"), *this);
65 : }
66 121 : }
67 :
68 : void
69 293 : PolynomialRegressionTrainer::preTrain()
70 : {
71 : _matrix.zero();
72 964 : for (unsigned int r = 0; r < _rhs.size(); ++r)
73 : _rhs[r].zero();
74 :
75 : /// Init calculators.
76 4652 : for (auto & calc : _calculators)
77 : calc->initializeCalculator();
78 :
79 293 : _r_sum.assign(_rhs.size(), 0.0);
80 293 : }
81 :
82 : void
83 8350 : PolynomialRegressionTrainer::train()
84 : {
85 : // Caching the different powers of data to accelerate the assembly of the
86 : // system
87 8350 : DenseMatrix<Real> data_pow(_n_dims, _max_degree + 1);
88 29500 : for (unsigned int d = 0; d < _n_dims; ++d)
89 99010 : for (unsigned int i = 0; i <= _max_degree; ++i)
90 77860 : data_pow(d, i) = pow(_predictor_row[d], i);
91 :
92 : // Emplace new values if necessary
93 8350 : if (_rvecval)
94 448 : for (unsigned int r = _rhs.size(); r < _rvecval->size(); ++r)
95 : {
96 118 : _rhs.emplace_back(_n_poly_terms, 0.0);
97 118 : _r_sum.emplace_back(0.0);
98 : }
99 :
100 125140 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
101 : {
102 116790 : Real i_value(1.0);
103 454590 : for (unsigned int ii = 0; ii < _n_dims; ++ii)
104 337800 : i_value *= data_pow(ii, _power_matrix[i][ii]);
105 :
106 2229420 : for (unsigned int j = 0; j < _n_poly_terms; ++j)
107 : {
108 : Real j_value(1.0);
109 8441870 : for (unsigned int jj = 0; jj < _n_dims; ++jj)
110 6329240 : j_value *= data_pow(jj, _power_matrix[j][jj]);
111 :
112 2112630 : _matrix(i, j) += i_value * j_value;
113 : }
114 :
115 : // Update calculators.
116 467160 : for (unsigned int c = i * 3; c < (i + 1) * 3; ++c)
117 350370 : _calculators[c]->updateCalculator(i_value);
118 :
119 116790 : if (_rval)
120 115650 : _rhs[0](i) += i_value * (*_rval);
121 1140 : else if (_rvecval)
122 10620 : for (unsigned int r = 0; r < _rvecval->size(); ++r)
123 9480 : _rhs[r](i) += i_value * (*_rvecval)[r];
124 : }
125 :
126 8350 : if (_rval)
127 8020 : _r_sum[0] += (*_rval);
128 330 : else if (_rvecval)
129 2990 : for (unsigned int r = 0; r < _rvecval->size(); ++r)
130 2660 : _r_sum[r] += (*_rvecval)[r];
131 :
132 : // Adding penalty term for Ridge regularization
133 8350 : if (_regression_type == "ridge")
134 0 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
135 0 : _matrix(i, i) += _penalty;
136 8350 : }
137 :
138 : void
139 293 : PolynomialRegressionTrainer::postTrain()
140 : {
141 : // Make sure _rhs are all the same size
142 293 : unsigned int nrval = _rhs.size();
143 : gatherMax(nrval);
144 321 : for (unsigned int r = _rhs.size(); r < nrval; ++r)
145 : {
146 28 : _rhs.emplace_back(_n_poly_terms, 0.0);
147 28 : _r_sum.emplace_back(0.0);
148 : }
149 :
150 : // Gather regression data
151 : gatherSum(_matrix.get_values());
152 1110 : for (auto & it : _rhs)
153 : gatherSum(it.get_values());
154 :
155 : // Gather response sums.
156 293 : gatherSum(_r_sum);
157 :
158 : // Finalize calculators.
159 4652 : for (auto & calc : _calculators)
160 4359 : calc->finalizeCalculator(true);
161 :
162 293 : std::vector<Real> mu(_n_poly_terms);
163 293 : std::vector<Real> sig(_n_poly_terms);
164 293 : std::vector<Real> sum_pf(_n_poly_terms);
165 :
166 1746 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
167 : {
168 : // To handle intercept, use mu = 0, sig = 1.
169 1453 : mu[i] = i > 0 ? _calculators[3 * i]->getValue() : 0.0;
170 1453 : sig[i] = i > 0 ? _calculators[3 * i + 1]->getValue() : 1.0;
171 1453 : sum_pf[i] = _calculators[3 * i + 2]->getValue();
172 : }
173 :
174 : // Transform _matrix and _rhs to match standardized features.
175 293 : const Real n = sum_pf[0]; // Sum of intercept term = n.
176 1746 : for (unsigned int i = 0; i < _n_poly_terms; ++i)
177 : {
178 11300 : for (unsigned int j = 0; j < _n_poly_terms; ++j)
179 : {
180 9847 : _matrix(i, j) -= (mu[j] * sum_pf[i] + mu[i] * sum_pf[j]);
181 9847 : _matrix(i, j) += n * mu[i] * mu[j];
182 9847 : _matrix(i, j) /= (sig[i] * sig[j]);
183 : }
184 4667 : for (unsigned int r = 0; r < _rhs.size(); ++r)
185 3214 : _rhs[r](i) = (_rhs[r](i) - mu[i] * _r_sum[r]) / sig[i];
186 : }
187 :
188 293 : DenseVector<Real> sol;
189 :
190 293 : _coeff.resize(_rhs.size());
191 1110 : for (unsigned int r = 0; r < _rhs.size(); ++r)
192 : {
193 817 : _matrix.lu_solve(_rhs[r], sol);
194 817 : _coeff[r] = sol.get_values();
195 :
196 : // Transform coefficients to match unstandardized features.
197 3214 : for (unsigned int i = 1; i < _n_poly_terms; ++i)
198 : {
199 2397 : _coeff[r][i] /= sig[i];
200 2397 : _coeff[r][0] -= _coeff[r][i] * mu[i];
201 : }
202 : }
203 293 : }
|
