1 //
2 // Created by friedrich on 16.08.16.
3 //
4
5 #include <AugmentedSubsetOfRegressors.hpp>
6 #include <assert.h>
7 #include <random>
8 #include <algorithm>
9 #include <iostream>
10
11 //--------------------------------------------------------------------------------
AugmentedSubsetOfRegressors(int n,double & delta_input)12 AugmentedSubsetOfRegressors::AugmentedSubsetOfRegressors(int n,
13 double &delta_input) :
14 SubsetOfRegressors(n, delta_input) {
15 }
16 //--------------------------------------------------------------------------------
17
build(std::vector<std::vector<double>> const & nodes,std::vector<double> const & values,std::vector<double> const & noise)18 void AugmentedSubsetOfRegressors::build(
19 std::vector<std::vector<double> > const &nodes,
20 std::vector<double> const &values, std::vector<double> const &noise) {
21
22 int nb_u_nodes;
23 if(nodes.size() < evaluations.active_index.size()){
24 nb_u_nodes = nodes.size();
25 }else if(nodes.size()*u_ratio < evaluations.active_index.size()){
26 nb_u_nodes = evaluations.active_index.size();
27 }else{
28 nb_u_nodes = nodes.size()*u_ratio;
29 }
30 std::cout << "In ApproxGaussianProcessAugmentedSoR_build ["
31 << nodes.size() << "," << nb_u_nodes << "]" << std::endl;
32
33 if (nb_u_nodes >= min_nb_u_nodes) {
34 this->sample_u(nb_u_nodes);
35 //std::cout << "U created: ";
36 //VectorOperations::print_matrix(u);
37
38 this->build_with_u(nodes, values, noise);
39 } else {
40 GaussianProcess::build(nodes, values, noise);
41 }
42 std::cout << "In ApproxGaussianProcessAugmentedSoR_build DONE" << std::endl;
43 return;
44 }
45 //--------------------------------------------------------------------------------
build_with_u(std::vector<std::vector<double>> const & nodes,std::vector<double> const & values,std::vector<double> const & noise)46 void AugmentedSubsetOfRegressors::build_with_u(
47 std::vector<std::vector<double> > const &nodes,
48 std::vector<double> const &values, std::vector<double> const &noise) {
49
50 int nb_u_nodes = u.size();
51 //std::cout << "In ApproxGaussianProcessAugmentedSoR_buil_with_u [" << nodes.size() << "," << nb_u_nodes <<"]" << std::endl;
52 if (nb_u_nodes >= min_nb_u_nodes) {
53
54 nb_gp_nodes = nodes.size();
55 gp_nodes.clear();
56 gp_noise.clear();
57 for (int i = 0; i < nb_gp_nodes; ++i) {
58 gp_nodes.push_back(nodes.at(i));
59 gp_noise.push_back(noise.at(i));
60 }
61
62 //std::cout << "U created: ";
63 //VectorOperations::print_matrix(u);
64 //std::cout << "Creating u done!" << std::endl;
65
66 //std::cout << "Noise: " << std::endl;
67 //VectorOperations::print_vector(gp_noise);
68
69 //Set up matrix K_u_f and K_f_u
70 K_u_f.clear();
71 K_u_f.resize(nb_u_nodes);
72 for (int i = 0; i < nb_u_nodes; ++i) {
73 for (int j = 0; j < nb_gp_nodes; ++j) {
74 K_u_f.at(i).push_back(evaluate_kernel(u.at(i), gp_nodes.at(j)));
75 }
76 }
77 K_f_u.clear();
78 K_f_u.resize(nb_gp_nodes);
79 for (int i = 0; i < nb_gp_nodes; ++i) {
80 K_f_u.at(i).resize(nb_u_nodes);
81 }
82 VectorOperations::mat_transpose(K_u_f, K_f_u);
83 //std::cout << "K_u_f and K_f_u done!" << std::endl;
84 //VectorOperations::print_matrix(K_u_f);
85
86 //Set up matrix K_u_u
87 L.clear();
88 L.resize(nb_u_nodes);
89 for (int i = 0; i < nb_u_nodes; ++i) {
90 for (int j = 0; j <= i; ++j) {
91 L.at(i).push_back(evaluate_kernel(u.at(i), u.at(j)));
92 }
93 }
94 //std::cout << "K_u_u done!" << std::endl;
95 //VectorOperations::print_matrix(L);
96
97 //Compute 1/noise^2*eye(length(f))*K_f_u (we save it directly in K_f_u since we do not change it later)
98 for (int i = 0; i < nb_gp_nodes; ++i) {
99 for (int j = 0; j < nb_u_nodes; ++j) {
100 K_f_u[i][j] *= 1.0
101 / pow(gp_noise.at(i) / 2e0 + noise_regularization, 2e0);
102 }
103 }
104 //std::cout << "eye(noise)*K_f_u" << std::endl;
105 //VectorOperations::print_matrix(K_f_u);
106
107 //Compute K_u_f*K_f_u = K_u_f_f_u
108 std::vector<std::vector<double> > K_u_f_f_u;
109 K_u_f_f_u.resize(nb_u_nodes);
110 for (int i = 0; i < nb_u_nodes; ++i) {
111 K_u_f_f_u[i].resize(nb_u_nodes);
112 }
113 VectorOperations::mat_product(K_u_f, K_f_u, K_u_f_f_u);
114
115 //Add K_u_f_f_u + K_u_u = Sigma_not_inv
116 for (int i = 0; i < nb_u_nodes; ++i) {
117 for (int j = 0; j <= i; ++j) {
118 L[i][j] += K_u_f_f_u[i][j];
119 //Nugget
120 if (i == j) {
121 L[i][j] += 0.0001;
122 }
123 }
124 }
125
126 //Compute Cholesky of K_u_u
127 CholeskyFactorization::compute(L, pos, rho, nb_u_nodes);
128 assert(pos == 0);
129
130 //std::cout << "f:" << std::endl;
131 //VectorOperations::print_vector(values);
132
133 //Set f and compute 1/noise^2*eye(length(f)) * f
134 scaled_function_values.clear();
135 scaled_function_values.resize(nb_gp_nodes);
136 for (int i = 0; i < nb_gp_nodes; i++) {
137 scaled_function_values.at(i) = values.at(i);
138 // scaled_function_values.at(i) = values.at(i) - min_function_value;
139 // scaled_function_values.at(i) /= 5e-1*( max_function_value-min_function_value );
140 // scaled_function_values.at(i) -= 1e0;
141 }
142 std::vector<double> noisy_values(scaled_function_values.size());
143 for (int i = 0; i < nb_gp_nodes; ++i) {
144 noisy_values[i] = scaled_function_values[i] * 1.0
145 / pow(gp_noise.at(i) / 2e0 + noise_regularization, 2e0);
146 }
147 //std::cout << "eye(noise) * f:" << std::endl;
148 //VectorOperations::print_vector(noisy_values);
149
150 //Compute K_u_f * f = alpha
151 alpha.clear();
152 alpha.resize(nb_u_nodes);
153 VectorOperations::mat_vec_product(K_u_f, noisy_values, alpha);
154
155 //std::cout << "Alpha:" << std::endl;
156 //VectorOperations::print_vector(alpha);
157
158 //Solve Sigma_not_inv^(-1)*alpha
159 forward_substitution(L, alpha);
160 backward_substitution(L, alpha);
161
162 //std::cout << "Sigma*K_u_f*noise*y:" << std::endl;
163 //VectorOperations::print_vector(alpha);
164
165 } else {
166 GaussianProcess::build(nodes, values, noise);
167 }
168 //std::cout << "In ApproxGaussianProcessAugmentedSoR_buil_with_u DONE" << std::endl;;
169 return;
170 }
171 //--------------------------------------------------------------------------------
172
173 //--------------------------------------------------------------------------------
update(std::vector<double> const & x,double & value,double & noise)174 void AugmentedSubsetOfRegressors::update(std::vector<double> const &x,
175 double &value, double &noise) {
176
177 int nb_u_nodes;
178 if(gp_nodes.size() < evaluations.active_index.size()){
179 nb_u_nodes = (gp_nodes.size()+1);
180 }else if((gp_nodes.size()+1)*u_ratio < evaluations.active_index.size()){
181 nb_u_nodes = evaluations.active_index.size();
182 }else{
183 nb_u_nodes = (gp_nodes.size()+1)*u_ratio;
184 }
185 //int nb_u_nodes = (gp_nodes.size()+1)*u_ratio > evaluations.active_index.size() ?
186 // (gp_nodes.size()+1)*u_ratio : evaluations.active_index.size();
187
188 augmented_u.clear();
189 std::cout << "In update Gaussian with [" << gp_nodes.size()+1 << ","
190 << nb_u_nodes << "]" << std::endl;
191 if (nb_u_nodes >= min_nb_u_nodes) {
192 //std::cout << "#Update" << std::endl;
193 std::vector<std::vector<double> > temp_nodes;
194 temp_nodes.resize(gp_nodes.size());
195 std::vector<double> temp_values;
196 std::vector<double> temp_noise;
197 for (int i = 0; i < gp_nodes.size(); ++i) {
198 for (int j = 0; j < gp_nodes.at(i).size(); ++j) {
199 temp_nodes.at(i).push_back(gp_nodes.at(i).at(j));
200 }
201 temp_values.push_back(scaled_function_values.at(i));
202 temp_noise.push_back(gp_noise.at(i));
203 }
204 temp_nodes.push_back(x);
205 temp_noise.push_back(noise);
206 temp_values.push_back(value);
207 // scaled_function_values.push_back ( ( value - min_function_value ) /
208 // ( 5e-1*( max_function_value-min_function_value ) ) - 1e0 );
209
210 this->build(temp_nodes, temp_values, temp_noise);
211 } else {
212 GaussianProcess::update(x, value, noise);
213 }
214 std::cout << "In update Gaussian with DONE" << std::endl;
215 return;
216 }
217 //--------------------------------------------------------------------------------
218
evaluate(std::vector<double> const & x,double & mean,double & variance)219 void AugmentedSubsetOfRegressors::evaluate(std::vector<double> const &x,
220 double &mean, double &variance) {
221 if (u.size() > 0) {
222 //std::cout << "ApproxAugmented::evaluate" <<"[" << gp_nodes.size() << "," << u.size() << "]" << std::endl;
223 augmented_u.clear();
224 for (int i = 0; i < x.size(); ++i) {
225 augmented_u.push_back(x[i]);
226 }
227
228 /*std::cout << "Augmented_u:";
229 VectorOperations::print_vector(augmented_u);
230 std::cout << "U:";
231 VectorOperations::print_matrix(u);*/
232 u.push_back(augmented_u);
233
234 std::vector<std::vector<double> > temp_nodes;
235 temp_nodes.resize(gp_nodes.size());
236 std::vector<double> temp_values;
237 std::vector<double> temp_noise;
238 for (int i = 0; i < gp_nodes.size(); ++i) {
239 for (int j = 0; j < gp_nodes.at(i).size(); ++j) {
240 temp_nodes.at(i).push_back(gp_nodes.at(i).at(j));
241 }
242 temp_values.push_back(scaled_function_values.at(i));
243 temp_noise.push_back(gp_noise.at(i));
244 }
245 this->build_with_u(temp_nodes, temp_values, temp_noise);
246
247 int nb_u_nodes = u.size();
248
249 K0.clear();
250 K0.resize(nb_u_nodes);
251 for (int i = 0; i < nb_u_nodes; i++) {
252 K0.at(i) = evaluate_kernel(x, u[i]);
253 }
254 /*
255 std::cout << "##############Evaluate################" << std::endl;
256 std::cout << "At x:"<< std::endl;
257 VectorOperations::print_vector(x);
258 std::cout << "K_star_u:" << std::endl;
259 VectorOperations::print_vector(K0);
260 */
261 mean = VectorOperations::dot_product(K0, alpha);
262 //std::cout << "Mean: " << mean << std::endl;
263
264 TriangularMatrixOperations::forward_substitution(L, K0);
265 variance = VectorOperations::dot_product(K0, K0);
266 /*
267 std::cout << "Variance: " << variance << std::endl;
268 std::cout << "######################################" << std::endl;
269 */
270 /*
271 evaluate_counter++;
272 assert(evaluate_counter<20*3);
273 */
274 u.pop_back();
275 } else {
276 GaussianProcess::evaluate(x, mean, variance);
277 }
278
279 //std::cout << "ApproxAugmented::evaluate DONE" << std::endl;
280 return;
281
282 }
283
get_induced_nodes(std::vector<std::vector<double>> & induced_nodes) const284 void AugmentedSubsetOfRegressors::get_induced_nodes(
285 std::vector<std::vector<double> > &induced_nodes) const {
286 //assert(induced_nodes.size()==u_idx.size());
287 for (int i = 0; i < u.size(); ++i) {
288 induced_nodes.push_back(u[i]);
289 }
290 if (augmented_u.size() > 0) {
291 induced_nodes.push_back(augmented_u);
292 }
293 return;
294 }
295