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 "LMC.h"
11 : #include "MooseRandom.h"
12 : #include "MathUtils.h"
13 :
14 : registerMooseObject("StochasticToolsApp", LMC);
15 :
16 : InputParameters
17 102 : LMC::validParams()
18 : {
19 102 : InputParameters params = CovarianceFunctionBase::validParams();
20 102 : params.addClassDescription("Covariance function for multioutput Gaussian Processes based on the "
21 : "Linear Model of Coregionalization (LMC).");
22 204 : params.addParam<unsigned int>(
23 204 : "num_latent_funcs", 1., "The number of latent functions for the expansion of the outputs.");
24 102 : params.makeParamRequired<unsigned int>("num_outputs");
25 102 : params.makeParamRequired<std::vector<UserObjectName>>("covariance_functions");
26 102 : return params;
27 0 : }
28 :
29 51 : LMC::LMC(const InputParameters & parameters)
30 : : CovarianceFunctionBase(parameters),
31 102 : _num_expansion_terms(getParam<unsigned int>("num_latent_funcs"))
32 : {
33 : // We use a random number generator to obtain the initial guess for the
34 : // hyperparams
35 : MooseRandom generator_latent;
36 51 : generator_latent.seed(0, 1980);
37 :
38 : // First add and initialize the a A coefficients in the (aa^T+lambda*I) matrix
39 102 : for (const auto exp_i : make_range(_num_expansion_terms))
40 : {
41 51 : const std::string a_coeff_name = "acoeff_" + std::to_string(exp_i);
42 51 : const std::string a_coeff_name_with_prefix = _name + ":" + a_coeff_name;
43 51 : _a_coeffs.push_back(
44 51 : &addVectorRealHyperParameter(a_coeff_name, std::vector(_num_outputs, 1.0), true));
45 :
46 : auto & acoeff_vector = _hp_map_vector_real[a_coeff_name_with_prefix];
47 153 : for (const auto out_i : make_range(_num_outputs))
48 102 : acoeff_vector[out_i] = 3.0 * generator_latent.rand(0) + 1.0;
49 : }
50 :
51 : // Then add and initialize the lambda coefficients in the (aa^T+lambda*I) matrix
52 102 : for (const auto exp_i : make_range(_num_expansion_terms))
53 : {
54 51 : const std::string lambda_name = "lambda_" + std::to_string(exp_i);
55 51 : const std::string lambda_name_with_prefix = _name + ":" + lambda_name;
56 51 : _lambdas.push_back(
57 51 : &addVectorRealHyperParameter(lambda_name, std::vector(_num_outputs, 1.0), true));
58 :
59 : auto & lambda_vector = _hp_map_vector_real[lambda_name_with_prefix];
60 153 : for (const auto out_i : make_range(_num_outputs))
61 102 : lambda_vector[out_i] = 3.0 * generator_latent.rand(0) + 1.0;
62 : }
63 51 : }
64 :
65 : void
66 17364 : LMC::computeCovarianceMatrix(RealEigenMatrix & K,
67 : const RealEigenMatrix & x,
68 : const RealEigenMatrix & xp,
69 : const bool is_self_covariance) const
70 : {
71 : // Create temporary vectors for constructing the covariance matrix
72 17364 : RealEigenMatrix K_params = RealEigenMatrix::Zero(x.rows(), xp.rows());
73 34728 : RealEigenMatrix B = RealEigenMatrix::Zero(_num_outputs, _num_outputs);
74 17364 : K = RealEigenMatrix::Zero(x.rows() * _num_outputs, xp.rows() * _num_outputs);
75 : RealEigenMatrix K_working =
76 17364 : RealEigenMatrix::Zero(x.rows() * _num_outputs, xp.rows() * _num_outputs);
77 :
78 : // For every expansion term we add the contribution to the covariance matrix
79 34728 : for (const auto exp_i : make_range(_num_expansion_terms))
80 : {
81 17364 : _covariance_functions[exp_i]->computeCovarianceMatrix(K_params, x, xp, is_self_covariance);
82 17364 : computeBMatrix(B, exp_i);
83 17364 : MathUtils::kron(K_working, B, K_params);
84 : K += K_working;
85 : }
86 17364 : }
87 :
88 : bool
89 119000 : LMC::computedKdhyper(RealEigenMatrix & dKdhp,
90 : const RealEigenMatrix & x,
91 : const std::string & hyper_param_name,
92 : unsigned int ind) const
93 : {
94 : // Early return in the paramter name is longer than the expected [name] prefix.
95 : // We prefix the parameter names with the name of the covariance function.
96 119000 : if (name().length() + 1 > hyper_param_name.length())
97 : return false;
98 :
99 : // Strip the prefix from the given parameter name
100 119000 : const std::string name_without_prefix = hyper_param_name.substr(name().length() + 1);
101 :
102 : // Check if the parameter is tunable
103 119000 : if (_tunable_hp.find(hyper_param_name) != _tunable_hp.end())
104 : {
105 68000 : const std::string acoeff_prefix = "acoeff_";
106 68000 : const std::string lambda_prefix = "lambda_";
107 :
108 : // Allocate storage for the factors of the total gradient matrix
109 136000 : RealEigenMatrix dBdhp = RealEigenMatrix::Zero(_num_outputs, _num_outputs);
110 68000 : RealEigenMatrix K_params = RealEigenMatrix::Zero(x.rows(), x.rows());
111 :
112 68000 : if (name_without_prefix.find(acoeff_prefix) != std::string::npos)
113 : {
114 : // Automatically grab the expansion index
115 34000 : const int number = std::stoi(name_without_prefix.substr(acoeff_prefix.length()));
116 34000 : computeAGradient(dBdhp, number, ind);
117 34000 : _covariance_functions[number]->computeCovarianceMatrix(K_params, x, x, true);
118 : }
119 34000 : else if (name_without_prefix.find(lambda_prefix) != std::string::npos)
120 : {
121 : // Automatically grab the expansion index
122 34000 : const int number = std::stoi(name_without_prefix.substr(lambda_prefix.length()));
123 34000 : computeLambdaGradient(dBdhp, number, ind);
124 34000 : _covariance_functions[number]->computeCovarianceMatrix(K_params, x, x, true);
125 : }
126 68000 : MathUtils::kron(dKdhp, dBdhp, K_params);
127 : return true;
128 : }
129 : else
130 : {
131 : // Allocate storage for the matrix factors
132 51000 : RealEigenMatrix B_tmp = RealEigenMatrix::Zero(_num_outputs, _num_outputs);
133 102000 : RealEigenMatrix B = RealEigenMatrix::Zero(_num_outputs, _num_outputs);
134 51000 : RealEigenMatrix dKdhp_sub = RealEigenMatrix::Zero(x.rows(), x.rows());
135 :
136 : // First, check the dependent covariances
137 : bool found = false;
138 102000 : for (const auto dependent_covar : _covariance_functions)
139 51000 : if (!found)
140 51000 : found = dependent_covar->computedKdhyper(dKdhp_sub, x, hyper_param_name, ind);
141 :
142 51000 : if (!found)
143 0 : mooseError("Hyperparameter ", hyper_param_name, "not found!");
144 :
145 : // Then we compute the output covariance
146 102000 : for (const auto exp_i : make_range(_num_expansion_terms))
147 : {
148 51000 : computeBMatrix(B_tmp, exp_i);
149 : B += B_tmp;
150 : }
151 :
152 51000 : MathUtils::kron(dKdhp, B, dKdhp_sub);
153 :
154 : return true;
155 : }
156 :
157 : return false;
158 : }
159 :
160 : void
161 68364 : LMC::computeBMatrix(RealEigenMatrix & Bmat, const unsigned int exp_i) const
162 : {
163 68364 : const auto & a_coeffs = *_a_coeffs[exp_i];
164 68364 : const auto & lambda_coeffs = *_lambdas[exp_i];
165 :
166 205092 : for (const auto row_i : make_range(_num_outputs))
167 410184 : for (const auto col_i : make_range(_num_outputs))
168 : {
169 273456 : Bmat(row_i, col_i) = a_coeffs[row_i] * a_coeffs[col_i];
170 273456 : if (row_i == col_i)
171 136728 : Bmat(row_i, col_i) += lambda_coeffs[col_i];
172 : }
173 68364 : }
174 :
175 : void
176 34000 : LMC::computeAGradient(RealEigenMatrix & grad,
177 : const unsigned int exp_i,
178 : const unsigned int index) const
179 : {
180 34000 : const auto & a_coeffs = *_a_coeffs[exp_i];
181 : // Add asserts here
182 : grad.setZero();
183 102000 : for (const auto col_i : make_range(_num_outputs))
184 68000 : grad(index, col_i) = a_coeffs[col_i];
185 34000 : const RealEigenMatrix transpose = grad.transpose();
186 34000 : grad = grad + transpose;
187 34000 : }
188 :
189 : void
190 34000 : LMC::computeLambdaGradient(RealEigenMatrix & grad,
191 : const unsigned int /*exp_i*/,
192 : const unsigned int index) const
193 : {
194 : grad.setZero();
195 34000 : grad(index, index) = 1.0;
196 34000 : }
|
