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 "MaternHalfIntCovariance.h"
11 : #include <cmath>
12 :
13 : registerMooseObject("StochasticToolsApp", MaternHalfIntCovariance);
14 :
15 : InputParameters
16 96 : MaternHalfIntCovariance::validParams()
17 : {
18 96 : InputParameters params = CovarianceFunctionBase::validParams();
19 96 : params.addClassDescription("Matern half-integer covariance function.");
20 192 : params.addRequiredParam<std::vector<Real>>("length_factor",
21 : "Length factors to use for Covariance Kernel");
22 192 : params.addRequiredParam<Real>("signal_variance",
23 : "Signal Variance ($\\sigma_f^2$) to use for kernel calculation.");
24 192 : params.addParam<Real>(
25 192 : "noise_variance", 0.0, "Noise Variance ($\\sigma_n^2$) to use for kernel calculation.");
26 192 : params.addRequiredParam<unsigned int>(
27 : "p", "Integer p to use for Matern Half Integer Covariance Kernel");
28 96 : return params;
29 0 : }
30 :
31 48 : MaternHalfIntCovariance::MaternHalfIntCovariance(const InputParameters & parameters)
32 : : CovarianceFunctionBase(parameters),
33 48 : _length_factor(addVectorRealHyperParameter(
34 : "length_factor", getParam<std::vector<Real>>("length_factor"), true)),
35 48 : _sigma_f_squared(
36 144 : addRealHyperParameter("signal_variance", getParam<Real>("signal_variance"), true)),
37 48 : _sigma_n_squared(
38 144 : addRealHyperParameter("noise_variance", getParam<Real>("noise_variance"), true)),
39 144 : _p(addRealHyperParameter("p", getParam<unsigned int>("p"), false))
40 : {
41 48 : }
42 :
43 : void
44 39048 : MaternHalfIntCovariance::computeCovarianceMatrix(RealEigenMatrix & K,
45 : const RealEigenMatrix & x,
46 : const RealEigenMatrix & xp,
47 : const bool is_self_covariance) const
48 : {
49 39048 : if ((unsigned)x.cols() != _length_factor.size())
50 0 : mooseError("length_factor size does not match dimension of trainer input.");
51 :
52 39048 : maternHalfIntFunction(
53 39048 : K, x, xp, _length_factor, _sigma_f_squared, _sigma_n_squared, _p, is_self_covariance);
54 39048 : }
55 :
56 : void
57 71048 : MaternHalfIntCovariance::maternHalfIntFunction(RealEigenMatrix & K,
58 : const RealEigenMatrix & x,
59 : const RealEigenMatrix & xp,
60 : const std::vector<Real> & length_factor,
61 : const Real sigma_f_squared,
62 : const Real sigma_n_squared,
63 : const unsigned int p,
64 : const bool is_self_covariance)
65 : {
66 71048 : unsigned int num_samples_x = x.rows();
67 71048 : unsigned int num_samples_xp = xp.rows();
68 71048 : unsigned int num_params_x = x.cols();
69 :
70 : mooseAssert(num_params_x == xp.cols(),
71 : "Number of parameters do not match in covariance kernel calculation");
72 :
73 : // This factor is used over and over, don't calculate each time
74 71048 : Real factor = sqrt(2 * p + 1);
75 :
76 1094348 : for (unsigned int ii = 0; ii < num_samples_x; ++ii)
77 : {
78 17100200 : for (unsigned int jj = 0; jj < num_samples_xp; ++jj)
79 : {
80 : // Compute distance per parameter, scaled by length factor
81 : Real r_scaled = 0;
82 48230700 : for (unsigned int kk = 0; kk < num_params_x; ++kk)
83 32153800 : r_scaled += pow((x(ii, kk) - xp(jj, kk)) / length_factor[kk], 2);
84 16076900 : r_scaled = sqrt(r_scaled);
85 : // tgamma(x+1) == x! when x is a natural number, which should always be the case for
86 : // MaternHalfInt
87 : Real summation = 0;
88 64307600 : for (unsigned int tt = 0; tt < p + 1; ++tt)
89 48230700 : summation += (tgamma(p + tt + 1) / (tgamma(tt + 1) * tgamma(p - tt + 1))) *
90 48230700 : pow(2 * factor * r_scaled, p - tt);
91 16076900 : K(ii, jj) = sigma_f_squared * std::exp(-factor * r_scaled) *
92 16076900 : (tgamma(p + 1) / (tgamma(2 * p + 1))) * summation;
93 : }
94 1023300 : if (is_self_covariance)
95 484300 : K(ii, ii) += sigma_n_squared;
96 : }
97 71048 : }
98 :
99 : bool
100 96000 : MaternHalfIntCovariance::computedKdhyper(RealEigenMatrix & dKdhp,
101 : const RealEigenMatrix & x,
102 : const std::string & hyper_param_name,
103 : unsigned int ind) const
104 : {
105 96000 : if (name().length() + 1 > hyper_param_name.length())
106 : return false;
107 :
108 96000 : const std::string name_without_prefix = hyper_param_name.substr(name().length() + 1);
109 :
110 96000 : if (name_without_prefix == "noise_variance")
111 : {
112 0 : maternHalfIntFunction(dKdhp, x, x, _length_factor, 0, 1, _p, true);
113 : return true;
114 : }
115 :
116 96000 : if (name_without_prefix == "signal_variance")
117 : {
118 32000 : maternHalfIntFunction(dKdhp, x, x, _length_factor, 1, 0, _p, false);
119 : return true;
120 : }
121 :
122 64000 : if (name_without_prefix == "length_factor")
123 : {
124 64000 : computedKdlf(dKdhp, x, _length_factor, _sigma_f_squared, _p, ind);
125 : return true;
126 : }
127 :
128 : return false;
129 : }
130 :
131 : void
132 64000 : MaternHalfIntCovariance::computedKdlf(RealEigenMatrix & K,
133 : const RealEigenMatrix & x,
134 : const std::vector<Real> & length_factor,
135 : const Real sigma_f_squared,
136 : const unsigned int p,
137 : const int ind)
138 : {
139 64000 : unsigned int num_samples_x = x.rows();
140 64000 : unsigned int num_params_x = x.cols();
141 :
142 : mooseAssert(ind < x.cols(), "Incorrect length factor index");
143 :
144 : // This factor is used over and over, don't calculate each time
145 64000 : Real factor = sqrt(2 * p + 1);
146 :
147 1024000 : for (unsigned int ii = 0; ii < num_samples_x; ++ii)
148 : {
149 16960000 : for (unsigned int jj = 0; jj < num_samples_x; ++jj)
150 : {
151 : // Compute distance per parameter, scaled by length factor
152 : Real r_scaled = 0;
153 48000000 : for (unsigned int kk = 0; kk < num_params_x; ++kk)
154 32000000 : r_scaled += pow((x(ii, kk) - x(jj, kk)) / length_factor[kk], 2);
155 16000000 : r_scaled = sqrt(r_scaled);
156 16000000 : if (r_scaled != 0)
157 : {
158 : // product rule to compute dK/dr_scaled
159 : // u'v
160 : Real summation = 0;
161 60160000 : for (unsigned int tt = 0; tt < p + 1; ++tt)
162 45120000 : summation += (tgamma(p + tt + 1) / (tgamma(tt + 1) * tgamma(p - tt + 1))) *
163 45120000 : pow(2 * factor * r_scaled, p - tt);
164 15040000 : K(ii, jj) = -factor * std::exp(-factor * r_scaled) * summation;
165 : // uv'
166 : // dont need tt=p, (p-tt) factor ->0. Also avoids unsigned integer subtraction wraparound
167 : summation = 0;
168 45120000 : for (unsigned int tt = 0; tt < p; ++tt)
169 30080000 : summation += (tgamma(p + tt + 1) / (tgamma(tt + 1) * tgamma(p - tt + 1))) * 2 * factor *
170 30080000 : (p - tt) * pow(2 * factor * r_scaled, p - tt - 1);
171 15040000 : K(ii, jj) += std::exp(-factor * r_scaled) * summation;
172 : // Apply chain rule for dr_scaled/dl_i
173 15040000 : K(ii, jj) *= -std::pow(x(ii, ind) - x(jj, ind), 2) /
174 15040000 : (std::pow(length_factor[ind], 3) * r_scaled) * sigma_f_squared *
175 15040000 : (tgamma(p + 1) / (tgamma(2 * p + 1)));
176 : }
177 : else // avoid div by 0. 0/0=0 scenario.
178 960000 : K(ii, jj) = 0;
179 : }
180 : }
181 64000 : }
|
