1 /* _______________________________________________________________________
2
3 DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
4 Copyright 2014-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Dakota directory.
7 _______________________________________________________________________ */
8
9 #include "ExperimentDataUtils.hpp"
10 #include "dakota_global_defs.hpp"
11 #include "DataMethod.hpp"
12 #include "DakotaResponse.hpp"
13 #include "NonDBayesCalibration.hpp"
14 // Boost.Test
15 #define BOOST_TEST_MODULE dakota_field_covariance_utils
16 #include <boost/test/included/unit_test.hpp>
17
18 #include <cassert>
19 #include <iostream>
20
21
22 namespace Dakota {
23 namespace TestFieldCovariance {
24
test_multiple_scalar_covariance_matrix()25 void test_multiple_scalar_covariance_matrix()
26 {
27 std::vector<RealMatrix> matrices;
28 std::vector<RealVector> diagonals;
29 RealVector scalars;
30 IntVector matrix_map_indices, diagonal_map_indices, scalar_map_indices;
31
32 int num_scalars = 3;
33 Real scalar_array[] = {1.,2.,4.};
34 int scalar_map_index_array[] = {0, 1, 2};
35 scalars.sizeUninitialized( num_scalars );
36 scalar_map_indices.sizeUninitialized( num_scalars );
37 for ( int i=0; i<num_scalars; i++ ){
38 scalars[i] = scalar_array[i];
39 scalar_map_indices[i] = scalar_map_index_array[i];
40 }
41
42 ExperimentCovariance exper_cov;
43 exper_cov.set_covariance_matrices( matrices, diagonals, scalars,
44 matrix_map_indices,
45 diagonal_map_indices,
46 scalar_map_indices );
47
48 // Test determinant and log_determinant
49 BOOST_CHECK_CLOSE(exper_cov.determinant(), 8.0, 1.0e-12);
50 BOOST_CHECK_CLOSE(exper_cov.log_determinant(), std::log(8.0), 1.0e-12);
51
52 int num_residuals = 3;
53 Real residual_array[] = {1.,2.,4.};
54 RealVector residual( Teuchos::Copy, residual_array, num_residuals );
55
56 // Test application of the covariance inverse to residual vector
57 Real prod = exper_cov.apply_experiment_covariance( residual );
58 BOOST_CHECK( std::abs( prod - 7. ) <
59 10.*std::numeric_limits<double>::epsilon() );
60
61 // Test application of the sqrt of the covariance inverse to residual vector
62 RealVector result;
63 exper_cov.apply_experiment_covariance_inverse_sqrt( residual, result );
64 prod = result.dot( result );
65 BOOST_CHECK( std::abs( prod - 7. ) <
66 10.*std::numeric_limits<double>::epsilon() );
67
68 // Test application of the sqrt of the covariance inverse to matrix of
69 // gradient vectors
70 RealMatrix scaled_grads;
71 Real grads_array[] = {1.,2.,2.,4.,4.,8.};
72 RealMatrix grads( Teuchos::Copy, grads_array, 2, 2, 3 );
73 exper_cov.apply_experiment_covariance_inverse_sqrt_to_gradients( grads,
74 scaled_grads);
75
76 RealMatrix grammian( grads.numRows(), grads.numRows(), false );
77 grammian.multiply( Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, scaled_grads,
78 scaled_grads, 0. );
79 BOOST_CHECK( std::abs( grammian(0,0) - 7. ) <
80 10.*std::numeric_limits<double>::epsilon() );
81 BOOST_CHECK( std::abs( grammian(1,1) - 28. ) <
82 10.*std::numeric_limits<double>::epsilon() );
83
84 // Test application of the sqrt of the covariance inverse to matrix of
85 // Hessian matrices
86 Real hessian_0_array[] = {4., 2., 2., 2.};
87 RealSymMatrix hessian_0( Teuchos::Copy, false, hessian_0_array, 2, 2 );
88 RealSymMatrix hessian_1( hessian_0 );
89 RealSymMatrix hessian_2( hessian_0 );
90 hessian_1 *= 2.;
91 hessian_2 *= 4.;
92
93 RealSymMatrixArray hessians( 3 );
94 hessians[0] = hessian_0;
95 hessians[1] = hessian_1;
96 hessians[2] = hessian_2;
97 RealSymMatrixArray scaled_hessians;
98 exper_cov.apply_experiment_covariance_inverse_sqrt_to_hessians( hessians,
99 scaled_hessians );
100
101 Real exact_scaled_hessian_0_array[] = {4.,2.,2.,2.};
102 Real exact_scaled_hessian_1_array[] = {8./std::sqrt(2.),4./std::sqrt(2.),
103 4./std::sqrt(2.),4./std::sqrt(2.)};
104 Real exact_scaled_hessian_2_array[] = {8., 4., 4., 4.};
105 RealSymMatrix exact_scaled_hessian_0( Teuchos::View, false,
106 exact_scaled_hessian_0_array, 2, 2 );
107 RealSymMatrix exact_scaled_hessian_1( Teuchos::View, false,
108 exact_scaled_hessian_1_array, 2, 2 );
109 RealSymMatrix exact_scaled_hessian_2( Teuchos::View, false,
110 exact_scaled_hessian_2_array, 2, 2 );
111
112 scaled_hessians[0] -= exact_scaled_hessian_0;
113 scaled_hessians[1] -= exact_scaled_hessian_1;
114 scaled_hessians[2] -= exact_scaled_hessian_2;
115
116 BOOST_CHECK( scaled_hessians[0].normInf() <
117 10.*std::numeric_limits<double>::epsilon() );
118 BOOST_CHECK( scaled_hessians[1].normInf() <
119 10.*std::numeric_limits<double>::epsilon() );
120 BOOST_CHECK( scaled_hessians[2].normInf() <
121 10.*std::numeric_limits<double>::epsilon() );
122
123 // Test extraction of main diagonal
124 RealVector diagonal;
125 exper_cov.get_main_diagonal( diagonal );
126
127 Real exact_diagonal_array[] = {1.,2.,4.};
128 RealVector exact_diagonal(Teuchos::View, exact_diagonal_array, 3 );
129
130 exact_diagonal -= diagonal;
131 BOOST_CHECK( exact_diagonal.normInf() <
132 10.*std::numeric_limits<double>::epsilon() );
133
134 }
135
test_single_diagonal_block_covariance_matrix()136 void test_single_diagonal_block_covariance_matrix()
137 {
138 std::vector<RealMatrix> matrices;
139 std::vector<RealVector> diagonals;
140 RealVector scalars;
141 IntVector matrix_map_indices, diagonal_map_indices, scalar_map_indices;
142
143 int num_diags = 1;
144 int num_diag_entries = 3;
145 Real diagonal_array[] = {1.,2.,4.};
146 diagonal_map_indices.sizeUninitialized( num_diags );
147 diagonal_map_indices[0] = 0;
148 diagonals.resize( num_diags );
149 diagonals[0].sizeUninitialized( num_diag_entries );
150 for ( int i=0; i<num_diag_entries; i++ ){
151 diagonals[0][i] = diagonal_array[i];
152 }
153
154 ExperimentCovariance exper_cov;
155 exper_cov.set_covariance_matrices( matrices, diagonals, scalars,
156 matrix_map_indices,
157 diagonal_map_indices,
158 scalar_map_indices );
159
160
161 // Test determinant and log_determinant
162 BOOST_CHECK_CLOSE(exper_cov.determinant(), 8.0, 1.0e-12);
163 BOOST_CHECK_CLOSE(exper_cov.log_determinant(), std::log(8.0), 1.0e-12);
164
165 int num_residuals = 3;
166 Real residual_array[] = {1.,2.,4.};
167 RealVector residual( Teuchos::Copy, residual_array, num_residuals );
168
169 // Test application of the covariance inverse to residual vector
170 Real prod = exper_cov.apply_experiment_covariance( residual );
171 BOOST_CHECK( std::abs( prod - 7. ) <
172 10.*std::numeric_limits<double>::epsilon() );
173
174 // Test application of the sqrt of the covariance inverse to residual vector
175 RealVector result;
176 exper_cov.apply_experiment_covariance_inverse_sqrt( residual, result );
177 prod = result.dot( result );
178 BOOST_CHECK( std::abs( prod - 7. ) <
179 10.*std::numeric_limits<double>::epsilon() );
180
181 // Test application of the sqrt of the covariance inverse to matrix of
182 // gradient vectors
183 RealMatrix scaled_grads;
184 Real grads_array[] = {1.,2.,2.,4.,4.,8.};
185 RealMatrix grads( Teuchos::Copy, grads_array, 2, 2, 3 );
186 exper_cov.apply_experiment_covariance_inverse_sqrt_to_gradients( grads,
187 scaled_grads);
188
189 RealMatrix grammian( grads.numRows(), grads.numRows(), false );
190 grammian.multiply( Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, scaled_grads,
191 scaled_grads, 0. );
192 BOOST_CHECK( std::abs( grammian(0,0) - 7. ) <
193 10.*std::numeric_limits<double>::epsilon() );
194 BOOST_CHECK( std::abs( grammian(1,1) - 28. ) <
195 10.*std::numeric_limits<double>::epsilon() );
196
197 // Test application of the sqrt of the covariance inverse to matrix of
198 // Hessian matrices
199 Real hessian_0_array[] = {4., 2., 2., 2.};
200 RealSymMatrix hessian_0( Teuchos::Copy, false, hessian_0_array, 2, 2 );
201 RealSymMatrix hessian_1( hessian_0 );
202 RealSymMatrix hessian_2( hessian_0 );
203 hessian_1 *= 2.;
204 hessian_2 *= 4.;
205
206 RealSymMatrixArray hessians( 3 );
207 hessians[0] = hessian_0;
208 hessians[1] = hessian_1;
209 hessians[2] = hessian_2;
210 RealSymMatrixArray scaled_hessians;
211 exper_cov.apply_experiment_covariance_inverse_sqrt_to_hessians( hessians,
212 scaled_hessians );
213
214 Real exact_scaled_hessian_0_array[] = {4.,2.,2.,2.};
215 Real exact_scaled_hessian_1_array[] = {8./std::sqrt(2.),4./std::sqrt(2.),
216 4./std::sqrt(2.),4./std::sqrt(2.)};
217 Real exact_scaled_hessian_2_array[] = {8., 4., 4., 4.};
218 RealSymMatrix exact_scaled_hessian_0( Teuchos::View, false,
219 exact_scaled_hessian_0_array, 2, 2 );
220 RealSymMatrix exact_scaled_hessian_1( Teuchos::View, false,
221 exact_scaled_hessian_1_array, 2, 2 );
222 RealSymMatrix exact_scaled_hessian_2( Teuchos::View, false,
223 exact_scaled_hessian_2_array, 2, 2 );
224
225 scaled_hessians[0] -= exact_scaled_hessian_0;
226 scaled_hessians[1] -= exact_scaled_hessian_1;
227 scaled_hessians[2] -= exact_scaled_hessian_2;
228
229 BOOST_CHECK( scaled_hessians[0].normInf() <
230 10.*std::numeric_limits<double>::epsilon() );
231 BOOST_CHECK( scaled_hessians[1].normInf() <
232 10.*std::numeric_limits<double>::epsilon() );
233 BOOST_CHECK( scaled_hessians[2].normInf() <
234 10.*std::numeric_limits<double>::epsilon() );
235 }
236
test_single_full_block_covariance_matrix()237 void test_single_full_block_covariance_matrix()
238 {
239 std::vector<RealMatrix> matrices;
240 std::vector<RealVector> diagonals;
241 RealVector scalars;
242 IntVector matrix_map_indices, diagonal_map_indices, scalar_map_indices;
243
244 int num_matrices = 1;
245 int num_matrix_rows = 3;
246 Real matrix_array[] = {1.,0.5,0.25,0.5,2.,0.5,0.25,0.5,4.};
247 matrix_map_indices.sizeUninitialized( num_matrices );
248 matrix_map_indices[0] = 0;
249 matrices.resize( num_matrices );
250 matrices[0].shapeUninitialized( num_matrix_rows, num_matrix_rows );
251 for ( int j=0; j<num_matrix_rows; j++ ){
252 for ( int i=0; i<num_matrix_rows; i++ )
253 matrices[0](i,j) = matrix_array[j*num_matrix_rows+i];
254 }
255
256 ExperimentCovariance exper_cov;
257 exper_cov.set_covariance_matrices( matrices, diagonals, scalars,
258 matrix_map_indices,
259 diagonal_map_indices,
260 scalar_map_indices );
261
262 // Test determinant and log_determinant
263 BOOST_CHECK_CLOSE(exper_cov.determinant(), 6.75, 1.0e-12);
264 BOOST_CHECK_CLOSE(exper_cov.log_determinant(), std::log(6.75), 1.0e-12);
265
266 int num_residuals = 3;
267 Real residual_array[] = {1.,2.,4.};
268 RealVector residual( Teuchos::Copy, residual_array, num_residuals );
269
270 // Test application of the covariance inverse to residual vector
271 Real prod = exper_cov.apply_experiment_covariance( residual );
272 BOOST_CHECK( std::abs( prod - 16./3. ) <
273 10.*std::numeric_limits<double>::epsilon() );
274
275 // Test application of the sqrt of the covariance inverse to residual vector
276 RealVector result;
277 exper_cov.apply_experiment_covariance_inverse_sqrt( residual, result );
278 prod = result.dot( result );
279 BOOST_CHECK( std::abs( prod - 16./3. ) <
280 10.*std::numeric_limits<double>::epsilon() );
281
282 // Test application of the sqrt of the covariance inverse to matrix of
283 // gradient vectors
284 RealMatrix scaled_grads;
285 Real grads_array[] = {1.,2.,2.,4.,4.,8.};
286 RealMatrix grads( Teuchos::Copy, grads_array, 2, 2, 3 );
287 exper_cov.apply_experiment_covariance_inverse_sqrt_to_gradients( grads,
288 scaled_grads);
289
290 RealMatrix grammian( grads.numRows(), grads.numRows(), false );
291 grammian.multiply( Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, scaled_grads,
292 scaled_grads, 0. );
293 BOOST_CHECK( std::abs( grammian(0,0) - 16./3. ) <
294 10.*std::numeric_limits<double>::epsilon() );
295 BOOST_CHECK( std::abs( grammian(1,1) - 64./3. ) <
296 20.*std::numeric_limits<double>::epsilon() );
297
298 // Test application of the sqrt of the covariance inverse to matrix of
299 // Hessian matrices
300 Real hessian_0_array[] = {4., 2., 2., 2.};
301 RealSymMatrix hessian_0( Teuchos::Copy, false, hessian_0_array, 2, 2 );
302 RealSymMatrix hessian_1( hessian_0 );
303 RealSymMatrix hessian_2( hessian_0 );
304 hessian_1 *= 2.;
305 hessian_2 *= 4.;
306
307 RealSymMatrixArray hessians( 3 );
308 hessians[0] = hessian_0;
309 hessians[1] = hessian_1;
310 hessians[2] = hessian_2;
311 RealSymMatrixArray scaled_hessians;
312 exper_cov.apply_experiment_covariance_inverse_sqrt_to_hessians( hessians,
313 scaled_hessians );
314
315 Real exact_scaled_hessian_0_array[] = {4.,2.,2.,2.};
316 Real exact_scaled_hessian_1_array[] = {4.53557367611073,2.26778683805536,
317 2.26778683805536,2.26778683805536};
318 Real exact_scaled_hessian_2_array[] = {6.98297248755176,3.49148624377588,
319 3.49148624377588,3.49148624377588};
320 RealSymMatrix exact_scaled_hessian_0( Teuchos::View, false,
321 exact_scaled_hessian_0_array, 2, 2 );
322 RealSymMatrix exact_scaled_hessian_1( Teuchos::View, false,
323 exact_scaled_hessian_1_array, 2, 2 );
324 RealSymMatrix exact_scaled_hessian_2( Teuchos::View, false,
325 exact_scaled_hessian_2_array, 2, 2 );
326
327 scaled_hessians[0] -= exact_scaled_hessian_0;
328 scaled_hessians[1] -= exact_scaled_hessian_1;
329 scaled_hessians[2] -= exact_scaled_hessian_2;
330
331 BOOST_CHECK( scaled_hessians[0].normInf() <
332 100.*std::numeric_limits<double>::epsilon() );
333 BOOST_CHECK( scaled_hessians[1].normInf() <
334 100.*std::numeric_limits<double>::epsilon() );
335 BOOST_CHECK( scaled_hessians[2].normInf() <
336 100.*std::numeric_limits<double>::epsilon() );
337
338 // Test extraction of main diagonal
339 RealVector diagonal;
340 exper_cov.get_main_diagonal( diagonal );
341
342 Real exact_diagonal_array[] = {1.,2.,4.};
343 RealVector exact_diagonal(Teuchos::View, exact_diagonal_array, 3 );
344
345 exact_diagonal -= diagonal;
346 BOOST_CHECK( exact_diagonal.normInf() <
347 10.*std::numeric_limits<double>::epsilon() );
348
349 }
350
test_mixed_scalar_diagonal_full_block_covariance_matrix()351 void test_mixed_scalar_diagonal_full_block_covariance_matrix()
352 {
353 std::vector<RealMatrix> matrices;
354 std::vector<RealVector> diagonals;
355 RealVector scalars;
356 IntVector matrix_map_indices, diagonal_map_indices, scalar_map_indices;
357
358 // Experiment covariance matrix consists of the following blocks
359 // scalar_1, diagonal_1, matrix_1, scalar_2, scalar_3,
360
361 // MATLAB CODE:
362 // A = [1,0.5,0.25;0.5,2,0.5;0.25,0.5,4.]; B = eye(9); B(1,1)=1.;
363 // B(8,8) = 2.; B(9,9) = 4.; B(3,3) = 2.; B(4,4) = 4.; B(5:7,5:7)=A;
364 // d = [1., 1., 2., 4., 1., 2., 4., 2., 4.]';
365 // chol(inv(B))*d; r'*r
366
367 // Generate scalar matrix blocks
368 int num_scalars = 3;
369 Real scalar_array[] = {1.,2.,4.};
370 int scalar_map_index_array[] = {0, 3, 4};
371 scalars.sizeUninitialized( num_scalars );
372 scalar_map_indices.sizeUninitialized( num_scalars );
373 for ( int i=0; i<num_scalars; i++ ){
374 scalars[i] = scalar_array[i];
375 scalar_map_indices[i] = scalar_map_index_array[i];
376 }
377
378 // Generate diagonal covariance matrix blocks
379 int num_diags = 1;
380 int num_diag_entries = 3;
381 Real diagonal_array[] = {1.,2.,4.};
382 diagonal_map_indices.sizeUninitialized( num_diags );
383 diagonal_map_indices[0] = 1;
384 diagonals.resize( num_diags );
385 diagonals[0].sizeUninitialized( num_diag_entries );
386 for ( int i=0; i<num_diag_entries; i++ ){
387 diagonals[0][i] = diagonal_array[i];
388 }
389
390 // Generate full covariance matrix blocks
391 int num_matrices=1;
392 int num_matrix_rows = 3;
393 Real matrix_array[] = {1.,0.5,0.25,0.5,2.,0.5,0.25,0.5,4.};
394 matrix_map_indices.sizeUninitialized( num_matrices );
395 matrix_map_indices[0] = 2;
396 matrices.resize( num_matrices );
397 matrices[0].shapeUninitialized( num_matrix_rows, num_matrix_rows );
398 for ( int j=0; j<num_matrix_rows; j++ ){
399 for ( int i=0; i<num_matrix_rows; i++ )
400 matrices[0](i,j) = matrix_array[j*num_matrix_rows+i];
401 }
402
403 ExperimentCovariance exper_cov;
404 exper_cov.set_covariance_matrices( matrices, diagonals, scalars,
405 matrix_map_indices,
406 diagonal_map_indices,
407 scalar_map_indices );
408
409 // Test determinant and log_determinant
410 BOOST_CHECK_CLOSE(exper_cov.determinant(), 432.0, 1.0e-12);
411 BOOST_CHECK_CLOSE(exper_cov.log_determinant(), std::log(432.0), 1.0e-12);
412
413 int num_residuals = 9;
414 Real residual_array[] = {1., 1., 2., 4., 1., 2., 4., 2., 4.};
415 RealVector residual( Teuchos::Copy, residual_array, num_residuals );
416
417 // Test application of the covariance inverse to residual vector
418 Real prod = exper_cov.apply_experiment_covariance( residual );
419 BOOST_CHECK( std::abs( prod - 58./3. ) <
420 20.*std::numeric_limits<double>::epsilon() );
421
422 // Test application of the sqrt of the covariance inverse to residual vector
423 RealVector result;
424 exper_cov.apply_experiment_covariance_inverse_sqrt( residual, result );
425 prod = result.dot( result );
426 BOOST_CHECK( std::abs( prod - 58./3. ) <
427 10.*std::numeric_limits<double>::epsilon() );
428
429 // Test application of the sqrt of the covariance inverse to matrix of
430 // gradient vectors
431 RealMatrix scaled_grads;
432 Real grads_array[] = { 1., 2., 1., 2., 2., 4., 4., 8., 1., 2., 2., 4., 4., 8.,
433 2., 4., 4., 8. };
434 RealMatrix grads( Teuchos::Copy, grads_array, 2, 2, 9 );
435 exper_cov.apply_experiment_covariance_inverse_sqrt_to_gradients( grads,
436 scaled_grads);
437
438 RealMatrix grammian( grads.numRows(), grads.numRows(), false );
439 grammian.multiply( Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, scaled_grads,
440 scaled_grads, 0. );
441 BOOST_CHECK( std::abs( grammian(0,0) - 58./3. ) <
442 20.*std::numeric_limits<double>::epsilon() );
443 BOOST_CHECK( std::abs( grammian(1,1) - 232./3. ) <
444 20.*std::numeric_limits<double>::epsilon() );
445
446 // Test application of the sqrt of the covariance inverse to matrix of
447 // Hessian matrices
448 Real hessian_0_array[] = {4., 2., 2., 2.};
449 RealSymMatrix hessian_0( Teuchos::Copy, false, hessian_0_array, 2, 2 );
450 RealSymMatrix hessian_1( hessian_0 ), hessian_2( hessian_0 ),
451 hessian_3( hessian_0 ), hessian_4( hessian_0 ),
452 hessian_5( hessian_0 ), hessian_6( hessian_0 ),
453 hessian_7( hessian_0 ), hessian_8( hessian_0 );
454 hessian_1 *= 1.; hessian_2 *= 2.; hessian_3 *= 4.; hessian_4 *= 1.;
455 hessian_5 *= 2.; hessian_6 *= 4.; hessian_7 *= 2.; hessian_8 *= 4.;
456
457 RealSymMatrixArray hessians( 9 );
458 hessians[0] = hessian_0; hessians[1] = hessian_1;
459 hessians[2] = hessian_2; hessians[3] = hessian_3;
460 hessians[4] = hessian_4; hessians[5] = hessian_5;
461 hessians[6] = hessian_6; hessians[7] = hessian_7;
462 hessians[8] = hessian_8;
463
464 RealSymMatrixArray scaled_hessians;
465 exper_cov.apply_experiment_covariance_inverse_sqrt_to_hessians( hessians,
466 scaled_hessians );
467
468 Real exact_scaled_hessian_0_array[] = {4.,2.,2.,2.};
469 Real exact_scaled_hessian_1_array[] = {4.,2.,2.,2.};
470 Real exact_scaled_hessian_2_array[] = {8./std::sqrt(2.),4./std::sqrt(2.),
471 4./std::sqrt(2.),4./std::sqrt(2.)};
472 Real exact_scaled_hessian_3_array[] = {8., 4., 4., 4.};
473 Real exact_scaled_hessian_4_array[] = {4.,2.,2.,2.};
474 Real exact_scaled_hessian_5_array[] = {4.53557367611073,2.26778683805536,
475 2.26778683805536,2.26778683805536};
476 Real exact_scaled_hessian_6_array[] = {6.98297248755176,3.49148624377588,
477 3.49148624377588,3.49148624377588};
478 Real exact_scaled_hessian_7_array[] = {8./std::sqrt(2.),4./std::sqrt(2.),
479 4./std::sqrt(2.),4./std::sqrt(2.)};
480 Real exact_scaled_hessian_8_array[] = {8., 4., 4., 4.};
481
482 RealSymMatrix exact_scaled_hessian_0( Teuchos::View, false,
483 exact_scaled_hessian_0_array, 2, 2 );
484 RealSymMatrix exact_scaled_hessian_1( Teuchos::View, false,
485 exact_scaled_hessian_1_array, 2, 2 );
486 RealSymMatrix exact_scaled_hessian_2( Teuchos::View, false,
487 exact_scaled_hessian_2_array, 2, 2 );
488 RealSymMatrix exact_scaled_hessian_3( Teuchos::View, false,
489 exact_scaled_hessian_3_array, 2, 2 );
490 RealSymMatrix exact_scaled_hessian_4( Teuchos::View, false,
491 exact_scaled_hessian_4_array, 2, 2 );
492 RealSymMatrix exact_scaled_hessian_5( Teuchos::View, false,
493 exact_scaled_hessian_5_array, 2, 2 );
494 RealSymMatrix exact_scaled_hessian_6( Teuchos::View, false,
495 exact_scaled_hessian_6_array, 2, 2 );
496 RealSymMatrix exact_scaled_hessian_7( Teuchos::View, false,
497 exact_scaled_hessian_7_array, 2, 2 );
498 RealSymMatrix exact_scaled_hessian_8( Teuchos::View, false,
499 exact_scaled_hessian_8_array, 2, 2 );
500
501 scaled_hessians[0] -= exact_scaled_hessian_0;
502 scaled_hessians[1] -= exact_scaled_hessian_1;
503 scaled_hessians[2] -= exact_scaled_hessian_2;
504 scaled_hessians[3] -= exact_scaled_hessian_3;
505 scaled_hessians[4] -= exact_scaled_hessian_4;
506 scaled_hessians[5] -= exact_scaled_hessian_5;
507 scaled_hessians[6] -= exact_scaled_hessian_6;
508 scaled_hessians[7] -= exact_scaled_hessian_7;
509 scaled_hessians[8] -= exact_scaled_hessian_8;
510
511 BOOST_CHECK( scaled_hessians[0].normInf() <
512 100.*std::numeric_limits<double>::epsilon() );
513 BOOST_CHECK( scaled_hessians[1].normInf() <
514 100.*std::numeric_limits<double>::epsilon() );
515 BOOST_CHECK( scaled_hessians[2].normInf() <
516 100.*std::numeric_limits<double>::epsilon() );
517 BOOST_CHECK( scaled_hessians[3].normInf() <
518 100.*std::numeric_limits<double>::epsilon() );
519 BOOST_CHECK( scaled_hessians[4].normInf() <
520 100.*std::numeric_limits<double>::epsilon() );
521 BOOST_CHECK( scaled_hessians[5].normInf() <
522 100.*std::numeric_limits<double>::epsilon() );
523 BOOST_CHECK( scaled_hessians[6].normInf() <
524 100.*std::numeric_limits<double>::epsilon() );
525 BOOST_CHECK( scaled_hessians[7].normInf() <
526 100.*std::numeric_limits<double>::epsilon() );
527 BOOST_CHECK( scaled_hessians[8].normInf() <
528 100.*std::numeric_limits<double>::epsilon() );
529
530
531 // Test extraction of main diagonal
532 RealVector diagonal;
533 exper_cov.get_main_diagonal( diagonal );
534
535 Real exact_diagonal_array[] = {1.,1.,2.,4.,1.,2.,4.,2.,4.};
536 RealVector exact_diagonal(Teuchos::View, exact_diagonal_array, 9 );
537
538 exact_diagonal -= diagonal;
539 BOOST_CHECK( exact_diagonal.normInf() <
540 10.*std::numeric_limits<double>::epsilon() );
541
542 // Test conversion to correlation matrix
543 // Matlab for correlation matrix
544 // std_dev = sqrt(diag(B)); correl_mat = diag(1./std_dev)*B*diag(1./std_dev);
545 RealSymMatrix exact_correl(9);
546 for (int i=0; i<9; ++i)
547 exact_correl(i,i) = 1.0;
548 // update off-diagonal entries for the full matrix block
549 for (int i=0; i<3; ++i)
550 for (int j=0; j<i; ++j)
551 exact_correl(4+i,4+j) = matrices[0](i,j) /
552 std::sqrt(matrices[0](i,i)) / std::sqrt(matrices[0](j,j));
553
554 RealSymMatrix calc_correl;
555 exper_cov.as_correlation(calc_correl);
556
557 calc_correl -= exact_correl;
558 BOOST_CHECK( calc_correl.normInf() <
559 10.*std::numeric_limits<double>::epsilon() );
560 }
561
test_linear_interpolate_1d_no_extrapolation()562 void test_linear_interpolate_1d_no_extrapolation()
563 {
564 // Generate field coordinates on [0,1]
565 int num_field_pts = 6;
566 Real field_pts_array[] = {1./6.,2./5.,7./12.,2./3.,6./7.,1.};
567 RealMatrix field_pts( Teuchos::View, field_pts_array, 1, num_field_pts, 1 );
568 // Generate field data which is the 3rd degree Legendre polynomial 1/2(5x^3-3x)
569 RealVector field_vals( num_field_pts, false );
570 for ( int i=0; i<num_field_pts; i++)
571 field_vals[i] = 0.5*(5.*std::pow(field_pts(i,0),3)-3.*field_pts(i,0));
572
573 // Generate simulation coordinates on equally spaced on [0,1]
574 int num_sim_pts = 11;
575 RealMatrix sim_pts( num_sim_pts, 1, false );
576 Real dx = 1./(num_sim_pts-1);
577 for ( int i=0; i<num_sim_pts; i++)
578 sim_pts(i,0) = i*dx;
579 // Generate simulation data which is the 2nd degree Legendre polynomial
580 // 1/2(3x^2-1)
581 RealVector sim_vals( num_sim_pts, false );
582 for ( int i=0; i<num_sim_pts; i++)
583 sim_vals[i] = 0.5*(3.*std::pow(sim_pts(i,0),2)-1.);
584 // Generate gradients of simulation data which is [3x,4x]
585 int num_vars = 2;
586 RealMatrix sim_grads( num_vars, num_sim_pts, false );
587 for ( int i=0; i<num_sim_pts; i++){
588 sim_grads(0,i) = 3.*sim_pts(i,0);
589 sim_grads(1,i) = 4.*sim_pts(i,0);
590 }
591 // Generate gradients of simulation data which we assume is also
592 // [3x+0.2,3x+0.05;3x+0.05,3x+0.1]
593 // Note interp sets num_vars from gradients and so hessian must be
594 // consistent with grads
595 RealSymMatrixArray sim_hessians( num_sim_pts );
596 for ( int i=0; i<num_sim_pts; i++){
597 sim_hessians[i].shapeUninitialized(num_vars);
598 //symmetric so do not have to set all entries, just upper triangular ones
599 sim_hessians[i](0,0)=3.*sim_pts(i,0)+0.2;
600 sim_hessians[i](0,1)=3.*sim_pts(i,0)+0.05;
601 sim_hessians[i](1,1)=3.*sim_pts(i,0)+0.1;
602 }
603
604 // Interpolate the simulation data onto the coordinates of the field data
605 RealVector interp_vals;
606 RealMatrix interp_grads;
607 RealSymMatrixArray interp_hessians;
608 linear_interpolate_1d( sim_pts, sim_vals, sim_grads, sim_hessians,
609 field_pts, interp_vals, interp_grads, interp_hessians );
610 interp_vals -= field_vals;
611
612 Real diff_array[] = {-2.16574074074074e-1,1.8e-1,3.91261574074074e-1,
613 4.29259259259259e-1,3.17084548104956e-1,0.0};
614 RealVector diff( Teuchos::View, diff_array, num_field_pts );
615 diff -= interp_vals;
616 BOOST_CHECK( diff.normInf() < 10.*std::numeric_limits<double>::epsilon() );
617
618 RealMatrix true_grads(num_vars,num_field_pts,false);
619 for ( int i=0; i<num_field_pts; i++){
620 true_grads(0,i) = 3.*field_pts(i,0);
621 true_grads(1,i) = 4.*field_pts(i,0);
622 }
623 true_grads -= interp_grads;
624 BOOST_CHECK( true_grads.normInf()<10.*std::numeric_limits<double>::epsilon() );
625
626 RealSymMatrixArray true_hessians( num_sim_pts );
627 for ( int i=0; i<num_field_pts; i++){
628 true_hessians[i].shapeUninitialized(num_vars);
629 //symmetric so do not have to set all entries, just upper triangular ones
630 true_hessians[i](0,0)=3.*field_pts(i,0)+0.2;
631 true_hessians[i](0,1)=3.*field_pts(i,0)+0.05;
632 true_hessians[i](1,1)=3.*field_pts(i,0)+0.1;
633 true_hessians[i] -= interp_hessians[i];
634 BOOST_CHECK( true_hessians[i].normInf()<
635 10.*std::numeric_limits<double>::epsilon() );
636 }
637 }
638
test_linear_interpolate_1d_with_extrapolation()639 void test_linear_interpolate_1d_with_extrapolation()
640 {
641 // Linear interpolate uses constant extrapolation
642
643 // Generate field coordinates on [-1/6,11/10]
644 int num_field_pts = 6;
645 Real field_pts_array[] = {-1./6.,2./5.,7./12.,2./3.,6./7.,11./10.};
646 RealMatrix field_pts( Teuchos::View, field_pts_array, 1, num_field_pts, 1 );
647 // Generate field data which is the 3rd degree Legendre polynomial 1/2(5x^3-3x)
648 RealVector field_vals( num_field_pts, false );
649 for ( int i=0; i<num_field_pts; i++)
650 field_vals[i] = 0.5*(5.*std::pow(field_pts(i,0),3)-3.*field_pts(i,0));
651
652 // Generate simulation coordinates on equally spaced on [0,1]
653 int num_sim_pts = 11;
654 RealMatrix sim_pts( num_sim_pts, 1, false );
655 Real dx = 1./(num_sim_pts-1);
656 for ( int i=0; i<num_sim_pts; i++)
657 sim_pts(i,0) = i*dx;
658 // Generate simulation data which is the 2nd degree Legendre polynomial
659 // 1/2(3x^2-1)
660 RealVector sim_vals( num_sim_pts, false );
661 for ( int i=0; i<num_sim_pts; i++)
662 sim_vals[i] = 0.5*(3.*std::pow(sim_pts(i,0),2)-1.);
663
664 // Interpolate the simulation data onto the coordinates of the field data
665 RealVector interp_vals;
666 RealMatrix sim_grads, interp_grads;
667 RealSymMatrixArray sim_hessians, interp_hessians;
668 linear_interpolate_1d( sim_pts, sim_vals, sim_grads, sim_hessians,
669 field_pts, interp_vals, interp_grads, interp_hessians );
670 interp_vals -= field_vals;
671
672 Real diff_array[] = {-7.38425925925926e-1,1.8e-1,3.91261574074074e-1,
673 4.29259259259259e-1,3.17084548104956e-1,-6.775e-1};
674 RealVector diff( Teuchos::View, diff_array, num_field_pts );
675 diff -= interp_vals;
676 BOOST_CHECK( diff.normInf() < 10.*std::numeric_limits<double>::epsilon() );
677 }
678
679 /*void test_build_hessian_of_sum_square_residuals_from_function_hessians()
680 {
681 int num_residuals = 3;
682 RealSymMatrixArray func_hessians( num_residuals );
683 RealMatrix func_gradients( 2, num_residuals, false );
684 RealVector residuals( num_residuals );
685
686 Real pts_array[] = {-1,-1,-0.5,0.5,1./3.,2./3.};
687 RealMatrix pts( Teuchos::View, pts_array, 2, 2, 3 );
688
689 for ( int i=0; i<num_residuals; i++ ){
690 Real x = pts(0,i), y = pts(1,i);
691 Real x2 = x*x, y2 = y*y;
692 residuals[i] = (2.*x2-y2)*(2.*x2-y2)/10.-x2*y2;
693 // The following will not work. Build hessian assumes
694 // residual = (approx-data)
695 //residuals[i] = x2*y2-(2.*x2-y2)*(2.*x2-y2)/10.;
696 func_gradients(0,i) = 4./5.*x*(2.*x2-y2);
697 func_gradients(1,i) = 2./5.*y*(y2-2.*x2);
698 func_hessians[i].shape( 2 );
699 func_hessians[i](0,0) = -4./5.*(y2-6.*x2);
700 func_hessians[i](1,0) = -8.*x*y/5.;
701 func_hessians[i](1,1) = 2./5.*(3.*y2-2.*x2);
702 }
703
704 ActiveSet set(num_residuals, 2); set.request_values(7);
705 Response resp(SIMULATION_RESPONSE, set);
706 resp.function_values(residuals);
707 resp.function_gradients(func_gradients);
708 resp.function_hessians(func_hessians);
709
710 // -------------------------------------- //
711 // Build hessian without noise covariance
712 // -------------------------------------- //
713
714 // If no noise covariance specify exper_cov as empty
715 ExperimentCovariance exper_cov;
716
717 RealSymMatrix ssr_hessian;
718 build_hessian_of_sum_square_residuals_from_response(resp, exper_cov,
719 ssr_hessian);
720 // hessian computed for ssr= r'r/2
721 RealSymMatrix truth_ssr_hessian( 2 );
722 for ( int i=0; i<num_residuals; i++ ){
723 Real x = pts(0,i), y = pts(1,i);
724 Real x2 = x*x, x3 = x2*x, x4 = x2*x2, x5 = x3*x2, x6 = x4*x2,
725 y2 = y*y, y3 = y2*y, y4 = y2*y2, y5 = y3*y2, y6 = y4*y2;
726 truth_ssr_hessian(0,0) += 2./25.*( 10.*x2*y2*(y2-6.*x2)-
727 (y2-14.*x2)*(2.*x2-y2)*(2.*x2-y2) );
728 truth_ssr_hessian(1,0) += 4./25.*y*x*( 10.*x2*y2-3.*(2.*x2-y2)*(2.*x2-y2) );
729 truth_ssr_hessian(1,1) += 1./25.*( 10.*x2*y2*(2.*x2-3.*y2) +
730 (7.*y2-2.*x2)*(2.*x2-y2)*(2.*x2-y2) );
731 }
732 // hack until build hessian can use covariance for multiple experiments
733 truth_ssr_hessian *=2;
734
735 truth_ssr_hessian -= ssr_hessian;
736 BOOST_CHECK( truth_ssr_hessian.normInf() <
737 10.*std::numeric_limits<double>::epsilon() );
738
739 // -------------------------------------- //
740 // Build hessian with noise covariance
741 // -------------------------------------- //
742
743 // Fill exper_cov with noise covariance
744 std::vector<RealMatrix> matrices;
745 std::vector<RealVector> diagonals;
746 RealVector scalars;
747 IntVector matrix_map_indices, diagonal_map_indices, scalar_map_indices;
748
749 // Experiment covariance matrix consists of the following blocks
750 // scalar_1, matrix_1
751
752 // MATLAB CODE:
753 //
754 // A = [1,0.5;0.5,2.]; S = eye(3); S(1,1)=1.; S(2:3,2:3)=A;
755 // U = chol( inv(S) );
756
757 // r = [-0.9 -0.05625 -0.04444444444444444]';
758 // g = [-0.8, -0.1, -0.05925925925925925; 0.4, -0.05, 0.05925925925925925];
759 // h1 = [4,-1.6;-1.6,0.4]; h2=[1 0.4;0.4,0.1];
760 // h3 =[0.1777777777777778,-0.3555555555555555;
761 // -0.3555555555555555, 0.4444444444444445];
762
763 // gnewton_hess = g*inv(S)*g';
764 // rs = r'*inv(S);
765 // hess = gnewton_hess + rs(1)*h1+rs(2)*h2+rs(3)*h3
766
767 // Generate scalar matrix blocks
768 int num_scalars = 1;
769 Real scalar_array[] = {1.};
770 int scalar_map_index_array[] = {0};
771 scalars.sizeUninitialized( num_scalars );
772 scalar_map_indices.sizeUninitialized( num_scalars );
773 for ( int i=0; i<num_scalars; i++ ){
774 scalars[i] = scalar_array[i];
775 scalar_map_indices[i] = scalar_map_index_array[i];
776 }
777
778 // Generate full covariance matrix blocks
779 int num_matrices=1;
780 int num_matrix_rows = 2;
781 Real matrix_array[] = {1.,0.5,0.5,2.};
782 matrix_map_indices.sizeUninitialized( num_matrices );
783 matrix_map_indices[0] = 1;
784 matrices.resize( num_matrices );
785 matrices[0].shapeUninitialized( num_matrix_rows, num_matrix_rows );
786 for ( int j=0; j<num_matrix_rows; j++ ){
787 for ( int i=0; i<num_matrix_rows; i++ )
788 matrices[0](i,j) = matrix_array[j*num_matrix_rows+i];
789 }
790
791 exper_cov.set_covariance_matrices( matrices, diagonals, scalars,
792 matrix_map_indices,
793 diagonal_map_indices,
794 scalar_map_indices );
795
796 // must reset ssr_hessian because if it is the right size
797 // the build_hessians... function will assume we want to add to it
798 ssr_hessian.shape(0);
799 build_hessian_of_sum_square_residuals_from_response(resp, exper_cov,
800 ssr_hessian);
801
802 Real truth_noise_scaled_ssr_hessian_array[] = {-3.00319615912208e+00,
803 1.10723495982755e+00,
804 1.10723495982755e+00,
805 -2.02746423672350e-01};
806 RealSymMatrix truth_noise_scaled_ssr_hessian( Teuchos::View, true,
807 truth_noise_scaled_ssr_hessian_array,
808 2, 2 );
809
810 // hack until build hessian can use covariance for multiple experiments
811 truth_noise_scaled_ssr_hessian *=2;
812
813 truth_noise_scaled_ssr_hessian -= ssr_hessian;
814 BOOST_CHECK( truth_noise_scaled_ssr_hessian.normInf() <
815 100.*std::numeric_limits<double>::epsilon() );
816
817 // -------------------------------------- //
818 // Build Gauss-Newton hessian with noise covariance
819 // -------------------------------------- //
820 ActiveSet set1(num_residuals, 2); set1.request_values(3);
821 Response resp1(SIMULATION_RESPONSE, set1);
822 resp1.function_values(residuals);
823 resp1.function_gradients(func_gradients);
824 resp1.function_hessians(func_hessians);
825
826 // must reset ssr_hessian because if it is the right size
827 // the build_hessians... function will assume we want to add to it
828 ssr_hessian.shape(0);
829 build_hessian_of_sum_square_residuals_from_response(resp1, exper_cov,
830 ssr_hessian);
831
832 Real truth_noise_scaled_gn_ssr_hessian_array[] = {6.50048990789732e-01,
833 -3.15445816186557e-01,
834 -3.15445816186557e-01,
835 1.66556927297668e-01};
836 RealSymMatrix truth_noise_scaled_gn_ssr_hessian( Teuchos::View, true,
837 truth_noise_scaled_gn_ssr_hessian_array,
838 2, 2 );
839
840 // hack until build hessian can use covariance for multiple experiments
841 truth_noise_scaled_gn_ssr_hessian *=2;
842
843 truth_noise_scaled_gn_ssr_hessian -= ssr_hessian;
844 BOOST_CHECK( truth_noise_scaled_gn_ssr_hessian.normInf() <
845 100.*std::numeric_limits<double>::epsilon() );
846 }*/
847
test_symmetric_eigenvalue_decomposition()848 void test_symmetric_eigenvalue_decomposition()
849 {
850 Real matrix_array[] = { 1.64, 0.48, 0.48, 1.36 };
851 RealSymMatrix matrix( Teuchos::View, false, matrix_array, 2, 2 );
852
853 RealVector eigenvalues;
854 RealMatrix eigenvectors;
855 symmetric_eigenvalue_decomposition( matrix, eigenvalues, eigenvectors );
856
857 Real truth_eigenvalues_array[] = {1.,2.};
858 RealVector truth_eigenvalues( Teuchos::View, truth_eigenvalues_array, 2 );
859
860 truth_eigenvalues -= eigenvalues;
861 BOOST_CHECK( truth_eigenvalues.normInf() <
862 10.*std::numeric_limits<double>::epsilon() );
863
864
865 Real truth_eigenvectors_array[] ={ 0.6, -0.8, -0.8, -0.6 };
866 RealMatrix truth_eigenvectors( Teuchos::View, truth_eigenvectors_array, 2,2,2);
867
868 truth_eigenvectors -= eigenvectors;
869 BOOST_CHECK( truth_eigenvectors.normInf() <
870 10.*std::numeric_limits<double>::epsilon() );
871 }
872
test_get_positive_definite_covariance_from_hessian()873 void test_get_positive_definite_covariance_from_hessian()
874 {
875 // uncorrelated prior
876 Real prior_chol_fact1[] = { 0.2, 0., 0., 0.2 };
877 RealMatrix prior_L1(Teuchos::View, prior_chol_fact1, 2, 2, 2);
878
879 // non positive definite matrix
880 Real misfit_h_array1[] = { 0.92, 1.44, 1.44, 0.08 };
881 RealSymMatrix misfit_hessian1(Teuchos::View, false, misfit_h_array1, 2, 2 );
882
883 RealSymMatrix covariance1;
884 NonDBayesCalibration::
885 get_positive_definite_covariance_from_hessian(misfit_hessian1, prior_L1,
886 covariance1, NORMAL_OUTPUT);
887
888 // MATLAB result (no truncation of eigenvalues)
889 //Real truth_cov1_array[] ={ 0.038703703703704, -0.002222222222222,
890 // -0.002222222222222, 0.04 };
891 // Dakota result (truncation of 1 eigenvalue)
892 Real truth_cov1_array[] ={ 3.81037037037037e-02, -1.42222222222222e-03,
893 -1.42222222222222e-03, 3.89333333333333e-02 };
894 RealSymMatrix truth_covariance1(Teuchos::View, false, truth_cov1_array, 2, 2);
895
896 truth_covariance1 -= covariance1;
897 BOOST_CHECK( truth_covariance1.normInf() <
898 10.*std::numeric_limits<double>::epsilon() );
899
900 //////////////////////////////////////////////////////////////////////////////
901
902 // correlated prior
903 Real prior_chol_fact2[] = { 0.2, 0.05, 0., 0.2 };
904 RealMatrix prior_L2(Teuchos::View, prior_chol_fact2, 2, 2, 2);
905
906 // positive definite matrix
907 Real misfit_h_array2[] = { 1.64, 0.48, 0.48, 1.36 };
908 RealSymMatrix misfit_hessian2(Teuchos::View, false, misfit_h_array2, 2, 2);
909
910 RealSymMatrix covariance2;
911 NonDBayesCalibration::
912 get_positive_definite_covariance_from_hessian(misfit_hessian2, prior_L2,
913 covariance2, NORMAL_OUTPUT);
914
915 // MATLAB and Dakota result (no truncation of eigenvalues):
916 Real truth_cov2_array[] = { 0.037120225312445, 0.008125330047527,
917 0.008125330047527, 0.039714838936807 };
918 RealSymMatrix truth_covariance2(Teuchos::View, false, truth_cov2_array, 2, 2);
919
920 truth_covariance2 -= covariance2;
921 BOOST_CHECK( truth_covariance2.normInf() <
922 10.*std::numeric_limits<double>::epsilon() );
923
924 //////////////////////////////////////////////////////////////////////////////
925 /*
926 // zero misfit matrix: proposal covariance = prior covariance
927 // can't perform symmetric eigen-deomposition for L'HL = 0
928 Real misfit_h_array3[] = { 0., 0., 0., 0. };
929 RealSymMatrix misfit_hessian3(Teuchos::View, false, misfit_h_array3, 2, 2);
930
931 RealSymMatrix covariance3;
932 NonDBayesCalibration::
933 get_positive_definite_covariance_from_hessian(misfit_hessian3, prior_L2,
934 covariance3, NORMAL_OUTPUT);
935
936 Real truth_cov3_array[] = { 4., 0.1, 0.1, 4.0025 };
937 RealSymMatrix truth_covariance3(Teuchos::View, false, truth_cov3_array, 2, 2);
938
939 truth_covariance3 -= covariance3;
940 BOOST_CHECK( truth_covariance3.normInf() <
941 10.*std::numeric_limits<double>::epsilon() );
942 */
943 }
944
test_matrix_symmetry()945 void test_matrix_symmetry()
946 {
947 // Test non-square matrix
948 RealMatrix test_rect_matrix(4,3, true);
949 bool is_symm = is_matrix_symmetric(test_rect_matrix);
950 BOOST_CHECK( !is_symm );
951
952 Real matrix_array[] = { 1.0, 0.7, 0.7, 0.7,
953 0.7, 1.0, 0.7, 0.7,
954 0.7, 0.7, 1.0, 0.7,
955 0.7, 0.7, 0.7, 1.0 };
956
957 // Test symmetric matrix
958 RealMatrix test_symm_mat(Teuchos::Copy, matrix_array, 4, 4, 4);
959 is_symm = is_matrix_symmetric(test_symm_mat);
960 BOOST_CHECK( is_symm );
961
962 // Test non-symmetric square matrix
963 RealMatrix test_nonsymm_mat(test_symm_mat);
964 test_nonsymm_mat(1,0) = 0.5;
965 is_symm = is_matrix_symmetric(test_nonsymm_mat);
966 BOOST_CHECK( !is_symm );
967 }
968
969 } // end namespace TestFieldCovariance
970 } // end namespace Dakota
971
972 // NOTE: Boost.Test framework provides the main progran driver
973
974 //____________________________________________________________________________//
975
BOOST_AUTO_TEST_CASE(test_main)976 BOOST_AUTO_TEST_CASE( test_main )
977 //int test_main( int argc, char* argv[] ) // note the name!
978 {
979 using namespace Dakota::TestFieldCovariance;
980
981 // Test ExperimentData covariance matrix
982 test_multiple_scalar_covariance_matrix();
983 test_single_diagonal_block_covariance_matrix();
984 test_single_full_block_covariance_matrix();
985 test_mixed_scalar_diagonal_full_block_covariance_matrix();
986
987 // Test field interpolation functions
988 test_linear_interpolate_1d_no_extrapolation();
989 test_linear_interpolate_1d_with_extrapolation();
990
991 // Test hessian functions
992 // Turn following test off until I can create an ExperimentData object
993 //test_build_hessian_of_sum_square_residuals_from_function_hessians();
994 test_get_positive_definite_covariance_from_hessian();
995
996 // Test linear algebra routines
997 test_symmetric_eigenvalue_decomposition();
998 test_matrix_symmetry();
999
1000 int run_result = 0;
1001 BOOST_CHECK( run_result == 0 || run_result == boost::exit_success );
1002
1003 // return boost::exit_success;
1004 }
1005