1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2010 Sandia National Laboratories. Developed at the
5 University of Wisconsin--Madison under SNL contract number
6 624796. The U.S. Government and the University of Wisconsin
7 retain certain rights to this software.
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 (lgpl.txt) along with this library; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22
23 (2010) kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26
27
28 /** \file ObjectiveFunctionTests.cpp
29 * \brief
30 * \author Jason Kraftcheck
31 */
32
33 #include "ObjectiveFunctionTests.hpp"
34
35
36 /** get a patchdata to use for testing */
37 static bool init_pd( PatchData& pd );
patch()38 PatchData& ObjectiveFunctionTests::patch() {
39 static PatchData the_pd;
40 static bool did_init = init_pd( the_pd );
41 CPPUNIT_ASSERT(did_init);
42 return the_pd;
43 }
init_pd(PatchData & pd)44 static bool init_pd( PatchData& pd ) {
45 MsqError err;
46 const double coords[] = { 0,0,0, 1,1,1, 2,2,2 };
47 const bool fixed[] = { false, true, true };
48 const size_t indices[] = { 0, 1, 2 };
49 pd.fill( 3, coords, 1, TRIANGLE, indices, fixed, err );
50 return !err;
51 }
52
53 /** Internal helper function for test_eval_type */
evaluate_internal(ObjectiveFunction::EvalType type,OFTestMode test_mode,ObjectiveFunction * of)54 double ObjectiveFunctionTests::evaluate_internal(
55 ObjectiveFunction::EvalType type,
56 OFTestMode test_mode,
57 ObjectiveFunction* of )
58 {
59 MsqPrintError err(cout);
60 vector<Vector3D> grad;
61 vector<SymMatrix3D> diag;
62 MsqHessian hess;
63 bool valid = false;
64 double result;
65
66 switch (test_mode) {
67 case EVAL:
68 valid = of->evaluate( type, patch(), result, OF_FREE_EVALS_ONLY, err );
69 break;
70 case GRAD:
71 valid = of->evaluate_with_gradient( type, patch(), result, grad, err );
72 break;
73 case DIAG:
74 valid = of->evaluate_with_Hessian_diagonal( type, patch(), result, grad, diag, err );
75 break;
76 case HESS:
77 hess.initialize( patch(), err );
78 ASSERT_NO_ERROR( err );
79 valid = of->evaluate_with_Hessian( type, patch(), result, grad, hess, err );
80 break;
81 default:
82 CPPUNIT_ASSERT(false);
83 }
84
85 ASSERT_NO_ERROR( err );
86 CPPUNIT_ASSERT(valid);
87 return result;
88 }
89
test_eval_type(ObjectiveFunction::EvalType type,OFTestMode test_mode,ObjectiveFunctionTemplate * of)90 void ObjectiveFunctionTests::test_eval_type(
91 ObjectiveFunction::EvalType type,
92 OFTestMode test_mode,
93 ObjectiveFunctionTemplate* of )
94 {
95 // define two sets of quality metric values
96 const double vals1[] = { 1.0, 3.0, 4.0, 8.0 };
97 const size_t vals1_len = sizeof(vals1)/sizeof(vals1[0]);
98 const double vals2[] = { 21.5, 11.1, 30.0, 0.5 };
99 const size_t vals2_len = sizeof(vals2)/sizeof(vals2[0]);
100
101 // create a quality metric to use
102 OFTestQM metric;
103 of->set_quality_metric( &metric );
104
105 // get some initial values to compare to
106 of->clear();
107 metric.set_values( vals1, vals1_len );
108 const double init1 = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
109 of->clear();
110 metric.set_values( vals2, vals2_len );
111 const double init2 = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
112 of->clear();
113 metric.append_values( vals1, vals1_len );
114 const double inits = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
115 of->clear();
116
117 double val, expected;
118 switch (type) {
119 case ObjectiveFunction::CALCULATE:
120
121 // first make sure we get back the same values
122 metric.set_values( vals1, vals1_len );
123 val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
124 CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
125 metric.set_values( vals2, vals2_len );
126 val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
127 CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
128
129 // now do something that should modify the accumulated value of the OF
130 evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
131
132 // check that the values are unchanged
133 metric.set_values( vals1, vals1_len );
134 val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
135 CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
136 metric.set_values( vals2, vals2_len );
137 val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
138 CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
139
140 break;
141
142 case ObjectiveFunction::ACCUMULATE:
143
144 // begin with first set
145 metric.set_values( vals1, vals1_len );
146 val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
147 CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
148
149 // add in second set
150 metric.set_values( vals2, vals2_len );
151 val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
152 CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
153
154 // clear
155 of->clear();
156
157 // begin with second set
158 metric.set_values( vals2, vals2_len );
159 val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
160 CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
161
162 // add in first set
163 metric.set_values( vals1, vals1_len );
164 val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
165 CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
166
167 break;
168
169 case ObjectiveFunction::SAVE:
170
171 // calculate value for first set twice
172 metric.set_values( vals1, vals1_len );
173 metric.append_values( vals1, vals1_len );
174 expected = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
175
176 // begin with the first set
177 metric.set_values( vals1, vals1_len );
178 val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
179 CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
180
181 // add the second set
182 metric.set_values( vals2, vals2_len );
183 val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
184 CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
185
186 // now save the second set of values - OF value should not change
187 metric.set_values( vals2, vals2_len );
188 val = evaluate_internal( ObjectiveFunction::SAVE, test_mode, of );
189 CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
190
191 // now replace the second set with the first
192 metric.set_values( vals1, vals1_len );
193 val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
194 CPPUNIT_ASSERT_DOUBLES_EQUAL( expected, val, 1e-6 );
195
196 // check that saved values are cleared
197 of->clear();
198 metric.set_values( vals2, vals2_len );
199 val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
200 CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
201
202 break;
203
204 case ObjectiveFunction::UPDATE:
205
206 // With no saved data, an update should produce the same
207 // result as CALCULATE
208 metric.set_values( vals1, vals1_len );
209 val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
210 CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
211
212 // Doing an update with a second patch should change
213 // the global value to that of the second set
214 metric.set_values( vals2, vals2_len );
215 val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
216 CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
217
218 // Now add in the second set again
219 metric.set_values( vals2, vals2_len );
220 evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
221
222 // Now replace one accumulation of the second set with the first
223 // Should result in the first set + the second set
224 metric.set_values( vals1, vals1_len );
225 val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
226 CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
227
228 break;
229
230 case ObjectiveFunction::TEMPORARY:
231
232 // With no saved data, an TEMPORARY should produce the same
233 // result as CALCULATE
234 metric.set_values( vals1, vals1_len );
235 val = evaluate_internal( ObjectiveFunction::TEMPORARY, test_mode, of );
236 CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
237
238 // Begin with two instances of the second set in the
239 // accumulated value, with one instance saved so it
240 // can be removed later
241 metric.set_values( vals2, vals2_len );
242 evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
243 evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
244
245 // Now do a temporary eval, replacing one instance of the
246 // second set wiht the first set
247 metric.set_values( vals1, vals1_len );
248 val = evaluate_internal( ObjectiveFunction::TEMPORARY, test_mode, of );
249
250 // TEMPORARY should produce the same value as UPDATE, but without
251 // modifying any internal state.
252 expected = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
253 CPPUNIT_ASSERT_DOUBLES_EQUAL( expected, val, 1e-6 );
254
255 break;
256
257 default:
258 CPPUNIT_ASSERT_MESSAGE("No test for specified evaluation type", false);
259 break;
260 }
261 }
262
test_value(const double * input_values,unsigned num_input_values,double expected_value,OFTestMode test_mode,ObjectiveFunctionTemplate * of)263 void ObjectiveFunctionTests::test_value(
264 const double* input_values,
265 unsigned num_input_values,
266 double expected_value,
267 OFTestMode test_mode,
268 ObjectiveFunctionTemplate* of )
269 {
270 OFTestQM metric( input_values, num_input_values );
271 of->set_quality_metric( &metric );
272
273 double val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
274 CPPUNIT_ASSERT_DOUBLES_EQUAL( expected_value, val, 1e-6 );
275 }
276
test_clone(ObjectiveFunctionTemplate * of)277 void ObjectiveFunctionTests::test_clone( ObjectiveFunctionTemplate* of )
278 {
279 const double some_vals[] = { 1, 2, 3, 4, 5, 6, 0.5 };
280 const unsigned num_vals = sizeof(some_vals)/sizeof(some_vals[0]);
281 OFTestQM metric( some_vals, num_vals );
282 of->set_quality_metric(&metric);
283
284 // test that we get the same value from both
285 auto_ptr<ObjectiveFunction> of2( of->clone() );
286 double exp_val, val;
287 exp_val = evaluate_internal( ObjectiveFunction::CALCULATE, EVAL, of );
288 val = evaluate_internal( ObjectiveFunction::CALCULATE, EVAL, of2.get() );
289 CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_val, val, 1e-12 );
290
291 // check if OF supports BCD -- if not then done
292 MsqError err;
293 of->evaluate( ObjectiveFunction::UPDATE, patch(), val, false, err );
294 if (err) {
295 err.clear();
296 return;
297 }
298
299 // build up some saved state in the objective function
300 of->clear();
301 const double vals1[] = { 1.0, 3.0, 4.0, 8.0 };
302 const size_t vals1_len = sizeof(vals1)/sizeof(vals1[0]);
303 const double vals2[] = { 21.5, 11.1, 30.0, 0.5 };
304 const size_t vals2_len = sizeof(vals2)/sizeof(vals2[0]);
305 metric.set_values( vals1, vals1_len );
306 evaluate_internal( ObjectiveFunction::SAVE, EVAL, of );
307 metric.append_values( vals2, vals2_len );
308 evaluate_internal( ObjectiveFunction::ACCUMULATE, EVAL, of );
309
310 // check that clone has same accumulated data
311 of2 = auto_ptr<ObjectiveFunction>(of->clone());
312 metric.set_values( some_vals, num_vals );
313 exp_val = evaluate_internal( ObjectiveFunction::UPDATE, EVAL, of );
314 val = evaluate_internal( ObjectiveFunction::UPDATE, EVAL, of2.get() );
315 CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_val, val, 1e-12 );
316 }
317
318 /** Test correct handling of QM negate flag */
test_negate_flag(OFTestMode test_mode,ObjectiveFunctionTemplate * of)319 void ObjectiveFunctionTests::test_negate_flag(
320 OFTestMode test_mode,
321 ObjectiveFunctionTemplate* of )
322 {
323 const double some_vals[] = { 1, 2, 3, 4, 5, 6, 0.5 };
324 const unsigned num_vals = sizeof(some_vals)/sizeof(some_vals[0]);
325 OFTestQM metric( some_vals, num_vals );
326 of->set_quality_metric(&metric);
327
328 MsqPrintError err(cout);
329 bool rval = false;
330 double value[2];
331 vector<Vector3D> grad[2];
332 vector<SymMatrix3D> diag[2];
333 MsqHessian hess[2];
334
335 // Do twice, once w/out negate flag set and then once
336 // with negate flag == -1.
337 ObjectiveFunction::EvalType type = ObjectiveFunction::CALCULATE;
338 for (unsigned i = 0; i < 2; ++i ) {
339 switch (test_mode) {
340 case EVAL:
341 rval = of->evaluate( type, patch(), value[i], false, err );
342 break;
343 case GRAD:
344 rval = of->evaluate_with_gradient( type, patch(), value[i], grad[i], err );
345 break;
346 case DIAG:
347 rval = of->evaluate_with_Hessian_diagonal( type, patch(), value[i], grad[i], diag[i], err );
348 break;
349 case HESS:
350 hess[i].initialize(patch(),err);
351 ASSERT_NO_ERROR( err );
352 rval = of->evaluate_with_Hessian( type, patch(), value[i], grad[i], hess[i], err );
353 break;
354 default:
355 CPPUNIT_ASSERT_MESSAGE("Invalid enum value in test code",false);
356 break;
357 }
358 ASSERT_NO_ERROR( err );
359 CPPUNIT_ASSERT(rval);
360 metric.set_negate_flag(-1);
361 }
362
363 switch (test_mode) {
364 case HESS:
365 CPPUNIT_ASSERT_EQUAL( hess[0].size(), hess[1].size() );
366 for (size_t r = 0; r < hess[0].size(); ++r)
367 for (size_t c = r; c < hess[0].size(); ++c)
368 if (hess[0].get_block(r,c))
369 CPPUNIT_ASSERT_MATRICES_EQUAL( -*hess[0].get_block(r,c),
370 *hess[1].get_block(r,c),
371 1e-6 );
372 case DIAG:
373 // NOTE: When case HESS: falls through to here, diag[0] and diag[1]
374 // will be empty, making this a no-op.
375 CPPUNIT_ASSERT_EQUAL( diag[0].size(), diag[1].size() );
376 for (size_t j = 0; j < diag[0].size(); ++j)
377 CPPUNIT_ASSERT_MATRICES_EQUAL( -diag[0][j], diag[1][j], 1e-6 );
378 case GRAD:
379 CPPUNIT_ASSERT_EQUAL( grad[0].size(), grad[1].size() );
380 for (size_t j = 0; j < grad[0].size(); ++j)
381 CPPUNIT_ASSERT_VECTORS_EQUAL( -grad[0][j], grad[1][j], 1e-6 );
382 default:
383 CPPUNIT_ASSERT_DOUBLES_EQUAL( -value[0], value[1], 1e-6 );
384 }
385 }
386
compare_numerical_gradient(ObjectiveFunction * of)387 void ObjectiveFunctionTests::compare_numerical_gradient( ObjectiveFunction* of )
388 {
389 MsqPrintError err(std::cout);
390 PatchData pd;
391 create_twelve_hex_patch( pd, err );
392 ASSERT_NO_ERROR( err );
393
394 std::vector<Vector3D> num_grad, ana_grad;
395 double num_val, ana_val;
396 bool valid;
397
398 valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, ana_val, ana_grad, err );
399 ASSERT_NO_ERROR( err );
400 CPPUNIT_ASSERT(valid);
401 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), ana_grad.size() );
402
403 valid = of->ObjectiveFunction::evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, num_val, num_grad, err );
404 ASSERT_NO_ERROR( err );
405 CPPUNIT_ASSERT(valid);
406 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), num_grad.size() );
407
408 CPPUNIT_ASSERT_DOUBLES_EQUAL( ana_val, num_val, 1e-6 );
409 for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
410 CPPUNIT_ASSERT_VECTORS_EQUAL( num_grad[i], ana_grad[i], 1e-3 );
411 }
412 }
413
compare_hessian_gradient(ObjectiveFunction * of)414 void ObjectiveFunctionTests::compare_hessian_gradient( ObjectiveFunction* of )
415 {
416 MsqPrintError err(std::cout);
417 PatchData pd;
418 create_twelve_hex_patch( pd, err );
419 ASSERT_NO_ERROR( err );
420
421 std::vector<Vector3D> grad, hess_grad;
422 MsqHessian hess;
423 double grad_val, hess_val;
424 bool valid;
425
426 valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, grad_val, grad, err );
427 ASSERT_NO_ERROR( err );
428 CPPUNIT_ASSERT(valid);
429 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), grad.size() );
430
431 hess.initialize( pd, err );
432 ASSERT_NO_ERROR( err );
433 valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, hess_val, hess_grad, hess, err );
434 ASSERT_NO_ERROR( err );
435 CPPUNIT_ASSERT(valid);
436 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess_grad.size() );
437
438 CPPUNIT_ASSERT_DOUBLES_EQUAL( grad_val, hess_val, 1e-6 );
439 for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
440 CPPUNIT_ASSERT_VECTORS_EQUAL( grad[i], hess_grad[i], 1e-6 );
441 }
442 }
443
compare_diagonal_gradient(ObjectiveFunction * of)444 void ObjectiveFunctionTests::compare_diagonal_gradient( ObjectiveFunction* of )
445 {
446 MsqPrintError err(std::cout);
447 PatchData pd;
448 create_twelve_hex_patch( pd, err );
449 ASSERT_NO_ERROR( err );
450
451 std::vector<Vector3D> grad, hess_grad;
452 std::vector<SymMatrix3D> hess;
453 double grad_val, hess_val;
454 bool valid;
455
456 valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, grad_val, grad, err );
457 ASSERT_NO_ERROR( err );
458 CPPUNIT_ASSERT(valid);
459 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), grad.size() );
460
461 valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, hess_val, hess_grad, hess, err );
462 ASSERT_NO_ERROR( err );
463 CPPUNIT_ASSERT(valid);
464 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess_grad.size() );
465
466 CPPUNIT_ASSERT_DOUBLES_EQUAL( grad_val, hess_val, 1e-6 );
467 for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
468 CPPUNIT_ASSERT_VECTORS_EQUAL( grad[i], hess_grad[i], 1e-6 );
469 }
470 }
471
compare_hessian_diagonal(ObjectiveFunction * of)472 void ObjectiveFunctionTests::compare_hessian_diagonal( ObjectiveFunction* of )
473 {
474 MsqPrintError err(std::cout);
475 PatchData pd;
476 create_twelve_hex_patch( pd, err );
477 ASSERT_NO_ERROR( err );
478
479 std::vector<Vector3D> diag_grad, hess_grad;
480 std::vector<SymMatrix3D> diag;
481 MsqHessian hess;
482 double diag_val, hess_val;
483 bool valid;
484
485 valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, diag_val, diag_grad, diag, err );
486 ASSERT_NO_ERROR( err );
487 CPPUNIT_ASSERT(valid);
488 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), diag_grad.size() );
489 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), diag.size() );
490
491 hess.initialize( pd, err );
492 ASSERT_NO_ERROR( err );
493 valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, hess_val, hess_grad, hess, err );
494 ASSERT_NO_ERROR( err );
495 CPPUNIT_ASSERT(valid);
496 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess_grad.size() );
497 CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess.size() );
498
499 CPPUNIT_ASSERT_DOUBLES_EQUAL( hess_val, diag_val, 1e-6 );
500 for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
501 CPPUNIT_ASSERT_VECTORS_EQUAL( hess_grad[i], diag_grad[i], 1e-6 );
502 CPPUNIT_ASSERT_MATRICES_EQUAL( *hess.get_block(i,i), diag[i], 1e-6 );
503 }
504 }
505
MAT_EPS(const Matrix3D & A)506 double MAT_EPS(const Matrix3D& A)
507 { return std::max( 0.001 * sqrt(Frobenius_2(A)), 0.001 ); }
508
509 #define CHECK_EQUAL_MATRICES(A,B) \
510 CPPUNIT_ASSERT_MATRICES_EQUAL( (A), (B), MAT_EPS(A) )
511
compare_numerical_hessian(ObjectiveFunction * of,bool diagonal_only)512 void ObjectiveFunctionTests::compare_numerical_hessian( ObjectiveFunction* of,
513 bool diagonal_only )
514 {
515 const double delta = 0.0001;
516
517 MsqPrintError err(std::cout);
518 PatchData pd;
519 create_qm_two_tet_patch( pd, err );
520 ASSERT_NO_ERROR( err );
521 CPPUNIT_ASSERT( pd.num_free_vertices() != 0 );
522
523 // get analytical Hessian from objective function
524 std::vector<Vector3D> grad;
525 std::vector<SymMatrix3D> diag;
526 MsqHessian hess;
527 hess.initialize( pd, err );
528 ASSERT_NO_ERROR( err );
529 double value;
530 bool valid;
531 if (diagonal_only)
532 valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, value, grad, diag, err );
533 else
534 valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, hess, err );
535 ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
536
537
538 // do numerical approximation of each block and compare to analytical value
539 for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
540 const size_t j_end = diagonal_only ? i+1 : pd.num_free_vertices();
541 for (size_t j = i; j < j_end; ++j) {
542 // do numerical approximation for block corresponding to
543 // coorindates for ith and jth vertices.
544 Matrix3D block;
545 for (int k = 0; k < 3; ++k) {
546 for (int m = 0; m < 3; ++m) {
547 double dk, dm, dkm;
548 Vector3D ik = pd.vertex_by_index(i);
549 Vector3D im = pd.vertex_by_index(j);
550
551 Vector3D delta_k(0.0); delta_k[k] = delta;
552 pd.move_vertex( delta_k, i, err ); ASSERT_NO_ERROR(err);
553 valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dk, true, err );
554 ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
555
556 Vector3D delta_m(0.0); delta_m[m] = delta;
557 pd.move_vertex( delta_m, j, err ); ASSERT_NO_ERROR(err);
558 valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dkm, true, err );
559 ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
560
561 // be careful here that we do the right thing if i==j
562 pd.set_vertex_coordinates( ik, i, err ); ASSERT_NO_ERROR(err);
563 pd.set_vertex_coordinates( im, j, err ); ASSERT_NO_ERROR(err);
564 pd.move_vertex( delta_m, j, err ); ASSERT_NO_ERROR(err);
565 valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dm, true, err );
566 ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
567
568 pd.set_vertex_coordinates( ik, i, err ); ASSERT_NO_ERROR(err);
569 pd.set_vertex_coordinates( im, j, err ); ASSERT_NO_ERROR(err);
570
571 block[k][m] = (dkm - dk - dm + value)/(delta*delta);
572 }
573 }
574 // compare to analytical value
575 if (diagonal_only) {
576 CPPUNIT_ASSERT(i == j); // see j_end above
577 CPPUNIT_ASSERT(i < diag.size());
578 CHECK_EQUAL_MATRICES( block, Matrix3D(diag[i]) );
579 }
580 else {
581 Matrix3D* m = hess.get_block( i, j );
582 Matrix3D* mt = hess.get_block( j, i );
583 if (NULL != m) {
584 CHECK_EQUAL_MATRICES( block, *m );
585 }
586 if (NULL != mt) {
587 CHECK_EQUAL_MATRICES( transpose(block), *m );
588 }
589 if (NULL == mt && NULL == m) {
590 CHECK_EQUAL_MATRICES( Matrix3D(0.0), block );
591 }
592 }
593 }
594 }
595 }
596
597
598
599 const size_t HANDLE_VAL = 0xDEADBEEF;
600 class OFTestBadQM : public QualityMetric
601 {
602 public:
603
OFTestBadQM(bool return_error)604 OFTestBadQM( bool return_error ) : mError(return_error) {}
605
get_metric_type() const606 virtual MetricType get_metric_type() const
607 { return ELEMENT_BASED; }
608
get_name() const609 virtual string get_name() const
610 { return "ObjectiveFunctionTests"; }
611
get_negate_flag() const612 virtual int get_negate_flag() const
613 { return 1; }
614
get_evaluations(PatchData &,vector<size_t> & h,bool,MsqError &)615 virtual void get_evaluations( PatchData&, vector<size_t>& h, bool, MsqError& )
616 {
617 h.clear();
618 h.push_back(HANDLE_VAL);
619 }
620
evaluate(PatchData &,size_t h,double &,MsqError & err)621 virtual bool evaluate( PatchData&, size_t h, double&, MsqError& err )
622 {
623 CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
624 if (mError)
625 MSQ_SETERR(err)(MsqError::INVALID_STATE);
626 return false;
627 }
628
evaluate_with_indices(PatchData &,size_t h,double &,vector<size_t> &,MsqError & err)629 virtual bool evaluate_with_indices( PatchData&, size_t h, double&, vector<size_t>&, MsqError& err )
630 {
631 CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
632 if (mError)
633 MSQ_SETERR(err)(MsqError::INVALID_STATE);
634 return false;
635 }
636
evaluate_with_gradient(PatchData &,size_t h,double &,vector<size_t> &,vector<Vector3D> &,MsqError & err)637 virtual bool evaluate_with_gradient( PatchData&, size_t h, double&, vector<size_t>&, vector<Vector3D>&, MsqError& err )
638 {
639 CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
640 if (mError)
641 MSQ_SETERR(err)(MsqError::INVALID_STATE);
642 return false;
643 }
644
evaluate_with_Hessian(PatchData &,size_t h,double &,vector<size_t> &,vector<Vector3D> &,vector<Matrix3D> &,MsqError & err)645 virtual bool evaluate_with_Hessian( PatchData&, size_t h, double&, vector<size_t>&, vector<Vector3D>&, vector<Matrix3D>&, MsqError& err )
646 {
647 CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
648 if (mError)
649 MSQ_SETERR(err)(MsqError::INVALID_STATE);
650 return false;
651 }
652
653 private:
654
655 bool mError;
656 };
657
test_handles_invalid_qm(OFTestMode test_mode,ObjectiveFunctionTemplate * of)658 void ObjectiveFunctionTests::test_handles_invalid_qm(
659 OFTestMode test_mode,
660 ObjectiveFunctionTemplate* of )
661 {
662 OFTestBadQM metric(false);
663 of->set_quality_metric( &metric );
664
665 MsqPrintError err(cout);
666 vector<Vector3D> grad;
667 vector<SymMatrix3D> diag;
668 MsqHessian hess;
669 double result;
670 bool valid = false;
671
672 switch (test_mode) {
673 case EVAL:
674 valid = of->evaluate( ObjectiveFunction::CALCULATE, patch(), result, OF_FREE_EVALS_ONLY, err );
675 break;
676 case GRAD:
677 valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, patch(), result, grad, err );
678 break;
679 case DIAG:
680 valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, patch(), result, grad, diag, err );
681 break;
682 case HESS:
683 hess.initialize( patch(), err );
684 ASSERT_NO_ERROR( err );
685 valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, patch(), result, grad, hess, err );
686 break;
687 default:
688 CPPUNIT_ASSERT(false);
689 }
690
691 ASSERT_NO_ERROR( err );
692 CPPUNIT_ASSERT(!valid);
693 }
694
695
test_handles_qm_error(OFTestMode test_mode,ObjectiveFunctionTemplate * of)696 void ObjectiveFunctionTests::test_handles_qm_error(
697 OFTestMode test_mode,
698 ObjectiveFunctionTemplate* of )
699
700 {
701 OFTestBadQM metric(true);
702 of->set_quality_metric( &metric );
703
704 MsqError err;
705 vector<Vector3D> grad;
706 vector<SymMatrix3D> diag;
707 MsqHessian hess;
708 double result;
709 bool valid;
710
711 switch (test_mode) {
712 case EVAL:
713 valid = of->evaluate( ObjectiveFunction::CALCULATE, patch(), result, OF_FREE_EVALS_ONLY, err );
714 break;
715 case GRAD:
716 valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, patch(), result, grad, err );
717 break;
718 case DIAG:
719 valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, patch(), result, grad, diag, err );
720 break;
721 case HESS:
722 hess.initialize( patch(), err );
723 ASSERT_NO_ERROR( err );
724 valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, patch(), result, grad, hess, err );
725 break;
726 default:
727 CPPUNIT_ASSERT(false);
728 }
729
730 CPPUNIT_ASSERT(err);
731 }
732