1// Copyright (C) 2016-2018 T. Zachary Laine
2//
3// Distributed under the Boost Software License, Version 1.0. (See
4// accompanying file LICENSE_1_0.txt or copy at
5// http://www.boost.org/LICENSE_1_0.txt)
6#include "autodiff.h"
7
8#include <iostream>
9
10#include <boost/yap/algorithm.hpp>
11#include <boost/polymorphic_cast.hpp>
12#include <boost/hana/for_each.hpp>
13
14#define BOOST_TEST_MODULE autodiff_test
15#include <boost/test/included/unit_test.hpp>
16
17
18double const Epsilon = 10.0e-6;
19#define CHECK_CLOSE(A,B) do { BOOST_CHECK_CLOSE(A,B,Epsilon); } while(0)
20
21using namespace AutoDiff;
22
23//[ autodiff_expr_template_decl
24template <boost::yap::expr_kind Kind, typename Tuple>
25struct autodiff_expr
26{
27 static boost::yap::expr_kind const kind = Kind;
28
29 Tuple elements;
30};
31
32BOOST_YAP_USER_UNARY_OPERATOR(negate, autodiff_expr, autodiff_expr)
33
34BOOST_YAP_USER_BINARY_OPERATOR(plus, autodiff_expr, autodiff_expr)
35BOOST_YAP_USER_BINARY_OPERATOR(minus, autodiff_expr, autodiff_expr)
36BOOST_YAP_USER_BINARY_OPERATOR(multiplies, autodiff_expr, autodiff_expr)
37BOOST_YAP_USER_BINARY_OPERATOR(divides, autodiff_expr, autodiff_expr)
38//]
39
40//[ autodiff_expr_literals_decl
41namespace autodiff_placeholders {
42
43 // This defines a placeholder literal operator that creates autodiff_expr
44 // placeholders.
45 BOOST_YAP_USER_LITERAL_PLACEHOLDER_OPERATOR(autodiff_expr)
46
47}
48//]
49
50//[ autodiff_function_terminals
51template <OPCODE Opcode>
52struct autodiff_fn_expr :
53 autodiff_expr<boost::yap::expr_kind::terminal, boost::hana::tuple<OPCODE>>
54{
55 autodiff_fn_expr () :
56 autodiff_expr {boost::hana::tuple<OPCODE>{Opcode}}
57 {}
58
59 BOOST_YAP_USER_CALL_OPERATOR_N(::autodiff_expr, 1);
60};
61
62// Someone included <math.h>, so we have to add trailing underscores.
63autodiff_fn_expr<OP_SIN> const sin_;
64autodiff_fn_expr<OP_COS> const cos_;
65autodiff_fn_expr<OP_SQRT> const sqrt_;
66//]
67
68//[ autodiff_xform
69struct xform
70{
71 // Create a var-node for each placeholder when we see it for the first
72 // time.
73 template <long long I>
74 Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
75 boost::yap::placeholder<I>)
76 {
77 if (list_.size() < I)
78 list_.resize(I);
79 auto & retval = list_[I - 1];
80 if (retval == nullptr)
81 retval = create_var_node();
82 return retval;
83 }
84
85 // Create a param-node for every numeric terminal in the expression.
86 Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>, double x)
87 { return create_param_node(x); }
88
89 // Create a "uary" node for each call expression, using its OPCODE.
90 template <typename Expr>
91 Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::call>,
92 OPCODE opcode, Expr const & expr)
93 {
94 return create_uary_op_node(
95 opcode,
96 boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr), *this)
97 );
98 }
99
100 template <typename Expr>
101 Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::negate>,
102 Expr const & expr)
103 {
104 return create_uary_op_node(
105 OP_NEG,
106 boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr), *this)
107 );
108 }
109
110 // Define a mapping from binary arithmetic expr_kind to OPCODE...
111 static OPCODE op_for_kind (boost::yap::expr_kind kind)
112 {
113 switch (kind) {
114 case boost::yap::expr_kind::plus: return OP_PLUS;
115 case boost::yap::expr_kind::minus: return OP_MINUS;
116 case boost::yap::expr_kind::multiplies: return OP_TIMES;
117 case boost::yap::expr_kind::divides: return OP_DIVID;
118 default: assert(!"This should never execute"); return OPCODE{};
119 }
120 assert(!"This should never execute");
121 return OPCODE{};
122 }
123
124 // ... and use it to handle all the binary arithmetic operators.
125 template <boost::yap::expr_kind Kind, typename Expr1, typename Expr2>
126 Node * operator() (boost::yap::expr_tag<Kind>, Expr1 const & expr1, Expr2 const & expr2)
127 {
128 return create_binary_op_node(
129 op_for_kind(Kind),
130 boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr1), *this),
131 boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr2), *this)
132 );
133 }
134
135 vector<Node *> & list_;
136};
137//]
138
139//[ autodiff_to_node
140template <typename Expr, typename ...T>
141Node * to_auto_diff_node (Expr const & expr, vector<Node *> & list, T ... args)
142{
143 Node * retval = nullptr;
144
145 // This fills in list as a side effect.
146 retval = boost::yap::transform(expr, xform{list});
147
148 assert(list.size() == sizeof...(args));
149
150 // Fill in the values of the value-nodes in list with the "args"
151 // parameter pack.
152 auto it = list.begin();
153 boost::hana::for_each(
154 boost::hana::make_tuple(args ...),
155 [&it](auto x) {
156 Node * n = *it;
157 VNode * v = boost::polymorphic_downcast<VNode *>(n);
158 v->val = x;
159 ++it;
160 }
161 );
162
163 return retval;
164}
165//]
166
167struct F{
168 F() { AutoDiff::autodiff_setup(); }
169 ~F(){ AutoDiff::autodiff_cleanup(); }
170};
171
172
173BOOST_FIXTURE_TEST_SUITE(all, F)
174
175//[ autodiff_original_node_builder
176Node* build_linear_fun1_manually(vector<Node*>& list)
177{
178 //f(x1,x2,x3) = -5*x1+sin(10)*x1+10*x2-x3/6
179 PNode* v5 = create_param_node(-5);
180 PNode* v10 = create_param_node(10);
181 PNode* v6 = create_param_node(6);
182 VNode* x1 = create_var_node();
183 VNode* x2 = create_var_node();
184 VNode* x3 = create_var_node();
185
186 OPNode* op1 = create_binary_op_node(OP_TIMES,v5,x1); //op1 = v5*x1
187 OPNode* op2 = create_uary_op_node(OP_SIN,v10); //op2 = sin(v10)
188 OPNode* op3 = create_binary_op_node(OP_TIMES,op2,x1); //op3 = op2*x1
189 OPNode* op4 = create_binary_op_node(OP_PLUS,op1,op3); //op4 = op1 + op3
190 OPNode* op5 = create_binary_op_node(OP_TIMES,v10,x2); //op5 = v10*x2
191 OPNode* op6 = create_binary_op_node(OP_PLUS,op4,op5); //op6 = op4+op5
192 OPNode* op7 = create_binary_op_node(OP_DIVID,x3,v6); //op7 = x3/v6
193 OPNode* op8 = create_binary_op_node(OP_MINUS,op6,op7); //op8 = op6 - op7
194 x1->val = -1.9;
195 x2->val = 2;
196 x3->val = 5./6.;
197 list.push_back(x1);
198 list.push_back(x2);
199 list.push_back(x3);
200 return op8;
201}
202//]
203
204//[ autodiff_yap_node_builder
205Node* build_linear_fun1(vector<Node*>& list)
206{
207 //f(x1,x2,x3) = -5*x1+sin(10)*x1+10*x2-x3/6
208 using namespace autodiff_placeholders;
209 return to_auto_diff_node(
210 -5 * 1_p + sin_(10) * 1_p + 10 * 2_p - 3_p / 6,
211 list,
212 -1.9,
213 2,
214 5./6.
215 );
216}
217//]
218
219Node* build_linear_function2_manually(vector<Node*>& list)
220{
221 //f(x1,x2,x3) = -5*x1+-10*x1+10*x2-x3/6
222 PNode* v5 = create_param_node(-5);
223 PNode* v10 = create_param_node(10);
224 PNode* v6 = create_param_node(6);
225 VNode* x1 = create_var_node();
226 VNode* x2 = create_var_node();
227 VNode* x3 = create_var_node();
228 list.push_back(x1);
229 list.push_back(x2);
230 list.push_back(x3);
231 OPNode* op1 = create_binary_op_node(OP_TIMES,v5,x1); //op1 = v5*x1
232 OPNode* op2 = create_uary_op_node(OP_NEG,v10); //op2 = -v10
233 OPNode* op3 = create_binary_op_node(OP_TIMES,op2,x1);//op3 = op2*x1
234 OPNode* op4 = create_binary_op_node(OP_PLUS,op1,op3);//op4 = op1 + op3
235 OPNode* op5 = create_binary_op_node(OP_TIMES,v10,x2);//op5 = v10*x2
236 OPNode* op6 = create_binary_op_node(OP_PLUS,op4,op5);//op6 = op4+op5
237 OPNode* op7 = create_binary_op_node(OP_DIVID,x3,v6); //op7 = x3/v6
238 OPNode* op8 = create_binary_op_node(OP_MINUS,op6,op7);//op8 = op6 - op7
239 x1->val = -1.9;
240 x2->val = 2;
241 x3->val = 5./6.;
242 return op8;
243}
244
245Node* build_linear_function2(vector<Node*>& list)
246{
247 //f(x1,x2,x3) = -5*x1+-10*x1+10*x2-x3/6
248 using namespace autodiff_placeholders;
249 auto ten = boost::yap::make_terminal<autodiff_expr>(t: 10);
250 return to_auto_diff_node(
251 -5 * 1_p + -ten * 1_p + 10 * 2_p - 3_p / 6,
252 list,
253 -1.9,
254 2,
255 5./6.
256 );
257}
258
259Node* build_nl_function1_manually(vector<Node*>& list)
260{
261// (x1*x2 * sin(x1))/x3 + x2*x4 - x1/x2
262 VNode* x1 = create_var_node();
263 VNode* x2 = create_var_node();
264 VNode* x3 = create_var_node();
265 VNode* x4 = create_var_node();
266 x1->val = -1.23;
267 x2->val = 7.1231;
268 x3->val = 2;
269 x4->val = -10;
270 list.push_back(x1);
271 list.push_back(x2);
272 list.push_back(x3);
273 list.push_back(x4);
274
275 OPNode* op1 = create_binary_op_node(OP_TIMES,x2,x1);
276 OPNode* op2 = create_uary_op_node(OP_SIN,x1);
277 OPNode* op3 = create_binary_op_node(OP_TIMES,op1,op2);
278 OPNode* op4 = create_binary_op_node(OP_DIVID,op3,x3);
279 OPNode* op5 = create_binary_op_node(OP_TIMES,x2,x4);
280 OPNode* op6 = create_binary_op_node(OP_PLUS,op4,op5);
281 OPNode* op7 = create_binary_op_node(OP_DIVID,x1,x2);
282 OPNode* op8 = create_binary_op_node(OP_MINUS,op6,op7);
283 return op8;
284}
285
286Node* build_nl_function1(vector<Node*>& list)
287{
288 // (x1*x2 * sin(x1))/x3 + x2*x4 - x1/x2
289 using namespace autodiff_placeholders;
290 return to_auto_diff_node(
291 (1_p * 2_p * sin_(1_p)) / 3_p + 2_p * 4_p - 1_p / 2_p,
292 list,
293 -1.23,
294 7.1231,
295 2,
296 -10
297 );
298}
299
300BOOST_AUTO_TEST_CASE( test_linear_fun1 )
301{
302 BOOST_TEST_MESSAGE("test_linear_fun1");
303 vector<Node*> list;
304 Node* root = build_linear_fun1(list);
305 vector<double> grad;
306 double val1 = grad_reverse(root,list,grad);
307 double val2 = eval_function(root);
308 double x1g[] = {-5.5440211108893697744548489936278,10.0,-0.16666666666666666666666666666667};
309
310 for(unsigned int i=0;i<3;i++){
311 CHECK_CLOSE(grad[i],x1g[i]);
312 }
313
314 double eval = 30.394751221800913;
315
316 CHECK_CLOSE(val1,eval);
317 CHECK_CLOSE(val2,eval);
318
319
320 EdgeSet s;
321 nonlinearEdges(root,s);
322 unsigned int n = nzHess(s);
323 BOOST_CHECK_EQUAL(n,0);
324}
325
326BOOST_AUTO_TEST_CASE( test_grad_sin )
327{
328 BOOST_TEST_MESSAGE("test_grad_sin");
329 VNode* x1 = create_var_node();
330 x1->val = 10;
331 OPNode* root = create_uary_op_node(OP_SIN,x1);
332 vector<Node*> nodes;
333 nodes.push_back(x1);
334 vector<double> grad;
335 grad_reverse(root,nodes,grad);
336 double x1g = -0.83907152907645244;
337 //the matlab give cos(10) = -0.839071529076452
338
339 CHECK_CLOSE(grad[0],x1g);
340 BOOST_CHECK_EQUAL(nodes.size(),1);
341
342 EdgeSet s;
343 nonlinearEdges(root,s);
344 unsigned int n = nzHess(s);
345 BOOST_CHECK_EQUAL(n,1);
346}
347
348BOOST_AUTO_TEST_CASE(test_grad_single_node)
349{
350 VNode* x1 = create_var_node();
351 x1->val = -2;
352 vector<Node*> nodes;
353 nodes.push_back(x1);
354 vector<double> grad;
355 double val = grad_reverse(x1,nodes,grad);
356 CHECK_CLOSE(grad[0],1);
357 CHECK_CLOSE(val,-2);
358
359 EdgeSet s;
360 unsigned int n = 0;
361 nonlinearEdges(x1,s);
362 n = nzHess(s);
363 BOOST_CHECK_EQUAL(n,0);
364
365 grad.clear();
366 nodes.clear();
367 PNode* p = create_param_node(-10);
368 //OPNode* op = create_binary_op_node(TIMES,p,create_param_node(2));
369 val = grad_reverse(p,nodes,grad);
370 BOOST_CHECK_EQUAL(grad.size(),0);
371 CHECK_CLOSE(val,-10);
372
373 s.clear();
374 nonlinearEdges(p,s);
375 n = nzHess(s);
376 BOOST_CHECK_EQUAL(n,0);
377}
378
379BOOST_AUTO_TEST_CASE(test_grad_neg)
380{
381 VNode* x1 = create_var_node();
382 x1->val = 10;
383 PNode* p2 = create_param_node(-1);
384 vector<Node*> nodes;
385 vector<double> grad;
386 nodes.push_back(x1);
387 Node* root = create_binary_op_node(OP_TIMES,x1,p2);
388 grad_reverse(root,nodes,grad);
389 CHECK_CLOSE(grad[0],-1);
390 BOOST_CHECK_EQUAL(nodes.size(),1);
391 nodes.clear();
392 grad.clear();
393 nodes.push_back(x1);
394 root = create_uary_op_node(OP_NEG,x1);
395 grad_reverse(root,nodes,grad);
396 CHECK_CLOSE(grad[0],-1);
397
398 EdgeSet s;
399 unsigned int n = 0;
400 nonlinearEdges(root,s);
401 n = nzHess(s);
402 BOOST_CHECK_EQUAL(n,0);
403}
404
405BOOST_AUTO_TEST_CASE( test_nl_function)
406{
407 vector<Node*> list;
408 Node* root = build_nl_function1(list);
409 double val = eval_function(root);
410 vector<double> grad;
411 grad_reverse(root,list,grad);
412 double eval =-66.929555552886214;
413 double gx[] = {-4.961306690356109,-9.444611307649055,-2.064383410399700,7.123100000000000};
414 CHECK_CLOSE(val,eval);
415
416 for(unsigned int i=0;i<4;i++)
417 {
418 CHECK_CLOSE(grad[i],gx[i]);
419 }
420 unsigned int nzgrad = nzGrad(root);
421 unsigned int tol = numTotalNodes(root);
422 BOOST_CHECK_EQUAL(nzgrad,4);
423 BOOST_CHECK_EQUAL(tol,16);
424
425 EdgeSet s;
426 nonlinearEdges(root,s);
427 unsigned int n = nzHess(s);
428 BOOST_CHECK_EQUAL(n,11);
429}
430
431BOOST_AUTO_TEST_CASE( test_hess_reverse_1)
432{
433 vector<Node*> nodes;
434 Node* root = build_linear_fun1(nodes);
435 vector<double> grad;
436 double val = grad_reverse(root,nodes,grad);
437 double eval = eval_function(root);
438// cout<<eval<<"\t"<<grad[0]<<"\t"<<grad[1]<<"\t"<<grad[2]<<"\t"<<endl;
439
440 CHECK_CLOSE(val,eval);
441
442 for(unsigned int i=0;i<nodes.size();i++)
443 {
444 static_cast<VNode*>(nodes[i])->u = 0;
445 }
446
447 static_cast<VNode*>(nodes[0])->u = 1;
448 double hval = 0;
449 vector<double> dhess;
450 hval = hess_reverse(root,nodes,dhess);
451 CHECK_CLOSE(hval,eval);
452 for(unsigned int i=0;i<dhess.size();i++)
453 {
454 CHECK_CLOSE(dhess[i],0);
455 }
456}
457
458BOOST_AUTO_TEST_CASE( test_hess_reverse_2)
459{
460 vector<Node*> nodes;
461 Node* root = build_linear_function2(nodes);
462 vector<double> grad;
463 double val = grad_reverse(root,nodes,grad);
464 double eval = eval_function(root);
465
466 CHECK_CLOSE(val,eval);
467 for(unsigned int i=0;i<nodes.size();i++)
468 {
469 static_cast<VNode*>(nodes[i])->u = 0;
470 }
471
472 static_cast<VNode*>(nodes[0])->u = 1;
473 double hval = 0;
474 vector<double> dhess;
475 hval = hess_reverse(root,nodes,dhess);
476 CHECK_CLOSE(hval,eval);
477
478 for(unsigned int i=0;i<dhess.size();i++)
479 {
480 CHECK_CLOSE(dhess[i],0);
481 }
482
483 EdgeSet s;
484 nonlinearEdges(root,s);
485 unsigned int n = nzHess(s);
486 BOOST_CHECK_EQUAL(n,0);
487}
488
489BOOST_AUTO_TEST_CASE( test_hess_reverse_4)
490{
491 vector<Node*> nodes;
492// Node* root = build_nl_function1(nodes);
493
494 VNode* x1 = create_var_node();
495 nodes.push_back(x1);
496 x1->val = 1;
497 x1->u =1;
498 Node* op = create_uary_op_node(OP_SIN,x1);
499 Node* root = create_uary_op_node(OP_SIN,op);
500 vector<double> grad;
501 double eval = eval_function(root);
502 vector<double> dhess;
503 double hval = hess_reverse(root,nodes,dhess);
504 CHECK_CLOSE(hval,eval);
505 BOOST_CHECK_EQUAL(dhess.size(),1);
506 CHECK_CLOSE(dhess[0], -0.778395788418109);
507
508 EdgeSet s;
509 nonlinearEdges(root,s);
510 unsigned int n = nzHess(s);
511 BOOST_CHECK_EQUAL(n,1);
512}
513
514BOOST_AUTO_TEST_CASE( test_hess_reverse_3)
515{
516 vector<Node*> nodes;
517 VNode* x1 = create_var_node();
518 VNode* x2 = create_var_node();
519 nodes.push_back(x1);
520 nodes.push_back(x2);
521 x1->val = 2.5;
522 x2->val = -9;
523 Node* op1 = create_binary_op_node(OP_TIMES,x1,x2);
524 Node* root = create_binary_op_node(OP_TIMES,x1,op1);
525 double eval = eval_function(root);
526 for(unsigned int i=0;i<nodes.size();i++)
527 {
528 static_cast<VNode*>(nodes[i])->u = 0;
529 }
530 static_cast<VNode*>(nodes[0])->u = 1;
531
532 vector<double> dhess;
533 double hval = hess_reverse(root,nodes,dhess);
534 BOOST_CHECK_EQUAL(dhess.size(),2);
535 CHECK_CLOSE(hval,eval);
536 double hx[]={-18,5};
537 for(unsigned int i=0;i<dhess.size();i++)
538 {
539 //Print("\t["<<i<<"]="<<dhess[i]);
540 CHECK_CLOSE(dhess[i],hx[i]);
541 }
542
543 EdgeSet s;
544 nonlinearEdges(root,s);
545 unsigned int n = nzHess(s);
546 BOOST_CHECK_EQUAL(n,3);
547}
548
549BOOST_AUTO_TEST_CASE( test_hess_reverse_5)
550{
551 vector<Node*> nodes;
552 VNode* x1 = create_var_node();
553 VNode* x2 = create_var_node();
554 nodes.push_back(x1);
555 nodes.push_back(x2);
556 x1->val = 2.5;
557 x2->val = -9;
558 Node* op1 = create_binary_op_node(OP_TIMES,x1,x1);
559 Node* op2 = create_binary_op_node(OP_TIMES,x2,x2);
560 Node* op3 = create_binary_op_node(OP_MINUS,op1,op2);
561 Node* op4 = create_binary_op_node(OP_PLUS,op1,op2);
562 Node* root = create_binary_op_node(OP_TIMES,op3,op4);
563
564 double eval = eval_function(root);
565
566 for(unsigned int i=0;i<nodes.size();i++)
567 {
568 static_cast<VNode*>(nodes[i])->u = 0;
569 }
570 static_cast<VNode*>(nodes[0])->u = 1;
571
572 vector<double> dhess;
573 double hval = hess_reverse(root,nodes,dhess);
574 CHECK_CLOSE(hval,eval);
575 double hx[] ={75,0};
576 for(unsigned int i=0;i<dhess.size();i++)
577 {
578 CHECK_CLOSE(dhess[i],hx[i]);
579 }
580
581 for(unsigned int i=0;i<nodes.size();i++)
582 {
583 static_cast<VNode*>(nodes[i])->u = 0;
584 }
585 static_cast<VNode*>(nodes[1])->u = 1;
586
587 double hx2[] = {0, -972};
588 hval = hess_reverse(root,nodes,dhess);
589 for(unsigned int i=0;i<dhess.size();i++)
590 {
591 CHECK_CLOSE(dhess[i],hx2[i]);
592 }
593
594 EdgeSet s;
595 nonlinearEdges(root,s);
596 unsigned int n = nzHess(s);
597 BOOST_CHECK_EQUAL(n,4);
598}
599BOOST_AUTO_TEST_CASE( test_hess_reverse_6)
600{
601 vector<Node*> nodes;
602// Node* root = build_nl_function1(nodes);
603
604 VNode* x1 = create_var_node();
605 VNode* x2 = create_var_node();
606 nodes.push_back(x1);
607 nodes.push_back(x2);
608 x1->val = 2.5;
609 x2->val = -9;
610 Node* root = create_binary_op_node(OP_POW,x1,x2);
611
612 double eval = eval_function(root);
613
614 static_cast<VNode*>(nodes[0])->u=1;static_cast<VNode*>(nodes[1])->u=0;
615 vector<double> dhess;
616 double hval = hess_reverse(root,nodes,dhess);
617 CHECK_CLOSE(hval,eval);
618 double hx1[] ={0.003774873600000 , -0.000759862823419};
619 double hx2[] ={-0.000759862823419, 0.000220093141567};
620 for(unsigned int i=0;i<dhess.size();i++)
621 {
622 CHECK_CLOSE(dhess[i],hx1[i]);
623 }
624 static_cast<VNode*>(nodes[0])->u=0;static_cast<VNode*>(nodes[1])->u=1;
625 hess_reverse(root,nodes,dhess);
626 for(unsigned int i=0;i<dhess.size();i++)
627 {
628 CHECK_CLOSE(dhess[i],hx2[i]);
629 }
630
631 EdgeSet s;
632 nonlinearEdges(root,s);
633 unsigned int n = nzHess(s);
634 BOOST_CHECK_EQUAL(n,4);
635}
636
637BOOST_AUTO_TEST_CASE( test_hess_reverse_7)
638{
639 vector<Node*> nodes;
640 Node* root = build_nl_function1(nodes);
641
642 double eval = eval_function(root);
643
644
645 vector<double> dhess;
646 double hx0[] ={-1.747958066718855,
647 -0.657091724418110,
648 2.410459188139686,
649 0};
650 double hx1[] ={ -0.657091724418110,
651 0.006806564792590,
652 -0.289815306593997,
653 1.000000000000000};
654 double hx2[] ={ 2.410459188139686,
655 -0.289815306593997,
656 2.064383410399700,
657 0};
658 double hx3[] ={0,1,0,0};
659 for(unsigned int i=0;i<nodes.size();i++)
660 {
661 static_cast<VNode*>(nodes[i])->u = 0;
662 }
663 static_cast<VNode*>(nodes[0])->u = 1;
664 double hval = hess_reverse(root,nodes,dhess);
665 CHECK_CLOSE(hval,eval);
666 for(unsigned int i=0;i<dhess.size();i++)
667 {
668 CHECK_CLOSE(dhess[i],hx0[i]);
669 }
670
671 for (unsigned int i = 0; i < nodes.size(); i++) {
672 static_cast<VNode*>(nodes[i])->u = 0;
673 }
674 static_cast<VNode*>(nodes[1])->u = 1;
675 hess_reverse(root, nodes, dhess);
676 for (unsigned int i = 0; i < dhess.size(); i++) {
677 CHECK_CLOSE(dhess[i], hx1[i]);
678 }
679
680 for (unsigned int i = 0; i < nodes.size(); i++) {
681 static_cast<VNode*>(nodes[i])->u = 0;
682 }
683 static_cast<VNode*>(nodes[2])->u = 1;
684 hess_reverse(root, nodes, dhess);
685 for (unsigned int i = 0; i < dhess.size(); i++) {
686 CHECK_CLOSE(dhess[i], hx2[i]);
687 }
688
689 for (unsigned int i = 0; i < nodes.size(); i++) {
690 static_cast<VNode*>(nodes[i])->u = 0;
691 }
692 static_cast<VNode*>(nodes[3])->u = 1;
693 hess_reverse(root, nodes, dhess);
694 for (unsigned i = 0; i < dhess.size(); i++) {
695 CHECK_CLOSE(dhess[i], hx3[i]);
696 }
697}
698
699#if FORWARD_ENABLED
700void test_hess_forward(Node* root, unsigned int& nvar)
701{
702 AutoDiff::num_var = nvar;
703 unsigned int len = (nvar+3)*nvar/2;
704 double* hess = new double[len];
705 hess_forward(root,nvar,&hess);
706 for(unsigned int i=0;i<len;i++){
707 cout<<"hess["<<i<<"]="<<hess[i]<<endl;
708 }
709 delete[] hess;
710}
711#endif
712
713BOOST_AUTO_TEST_CASE( test_hess_reverse_8)
714{
715 vector<Node*> list;
716 vector<double> dhess;
717
718 VNode* x1 = create_var_node();
719 list.push_back(x1);
720 static_cast<VNode*>(list[0])->val = -10.5;
721 static_cast<VNode*>(list[0])->u = 1;
722 double deval = hess_reverse(x1,list,dhess);
723 CHECK_CLOSE(deval,-10.5);
724 BOOST_CHECK_EQUAL(dhess.size(),1);
725 BOOST_CHECK(isnan(dhess[0]));
726
727 EdgeSet s;
728 nonlinearEdges(x1,s);
729 unsigned int n = nzHess(s);
730 BOOST_CHECK_EQUAL(n,0);
731
732 PNode* p1 = create_param_node(-1.5);
733 list.clear();
734 deval = hess_reverse(p1,list,dhess);
735 CHECK_CLOSE(deval,-1.5);
736 BOOST_CHECK_EQUAL(dhess.size(),0);
737
738 s.clear();
739 nonlinearEdges(p1,s);
740 n = nzHess(s);
741 BOOST_CHECK_EQUAL(n,0);
742}
743
744BOOST_AUTO_TEST_CASE( test_hess_revers9)
745{
746 vector<Node*> list;
747 vector<double> dhess;
748 VNode* x1 = create_var_node();
749 list.push_back(x1);
750 static_cast<VNode*>(list[0])->val = 2.5;
751 static_cast<VNode*>(list[0])->u =1;
752 Node* op1 = create_binary_op_node(OP_TIMES,x1,x1);
753 Node* root = create_binary_op_node(OP_TIMES,op1,op1);
754 double deval = hess_reverse(root,list,dhess);
755 double eval = eval_function(root);
756 CHECK_CLOSE(eval,deval);
757 BOOST_CHECK_EQUAL(dhess.size(),1);
758 CHECK_CLOSE(dhess[0],75);
759
760 EdgeSet s;
761 nonlinearEdges(root,s);
762 unsigned int n = nzHess(s);
763 BOOST_CHECK_EQUAL(n,1);
764
765}
766
767BOOST_AUTO_TEST_CASE( test_hess_revers10)
768{
769 vector<Node*> list;
770 vector<double> dhess;
771 VNode* x1 = create_var_node();
772 VNode* x2 = create_var_node();
773 list.push_back(x1);
774 list.push_back(x2);
775 Node* op1 = create_binary_op_node(OP_TIMES, x1,x2);
776 Node* op2 = create_uary_op_node(OP_SIN,op1);
777 Node* op3 = create_uary_op_node(OP_COS,op1);
778 Node* root = create_binary_op_node(OP_TIMES, op2, op3);
779 static_cast<VNode*>(list[0])->val = 2.1;
780 static_cast<VNode*>(list[1])->val = 1.8;
781 double eval = eval_function(root);
782
783 //second column
784 static_cast<VNode*>(list[0])->u = 0;
785 static_cast<VNode*>(list[1])->u = 1;
786 double deval = hess_reverse(root,list,dhess);
787 CHECK_CLOSE(eval,deval);
788 BOOST_CHECK_EQUAL(dhess.size(),2);
789 CHECK_CLOSE(dhess[0], -6.945893481707861);
790 CHECK_CLOSE(dhess[1], -8.441601940854081);
791
792 //first column
793 static_cast<VNode*>(list[0])->u = 1;
794 static_cast<VNode*>(list[1])->u = 0;
795 deval = hess_reverse(root,list,dhess);
796 CHECK_CLOSE(eval,deval);
797 BOOST_CHECK_EQUAL(dhess.size(),2);
798 CHECK_CLOSE(dhess[0], -6.201993262668304);
799 CHECK_CLOSE(dhess[1], -6.945893481707861);
800}
801
802BOOST_AUTO_TEST_CASE( test_grad_reverse11)
803{
804 vector<Node*> list;
805 VNode* x1 = create_var_node();
806 Node* p2 = create_param_node(2);
807 list.push_back(x1);
808 Node* op1 = create_binary_op_node(OP_POW,x1,p2);
809 static_cast<VNode*>(x1)->val = 0;
810 vector<double> grad;
811 grad_reverse(op1,list,grad);
812 BOOST_CHECK_EQUAL(grad.size(),1);
813 CHECK_CLOSE(grad[0],0);
814}
815
816BOOST_AUTO_TEST_CASE( test_hess_reverse12)
817{
818 vector<Node*> list;
819 VNode* x1 = create_var_node();
820 Node* p2 = create_param_node(2);
821 list.push_back(x1);
822 Node* op1 = create_binary_op_node(OP_POW,x1,p2);
823 x1->val = 0;
824 x1->u = 1;
825 vector<double> hess;
826 hess_reverse(op1,list,hess);
827 BOOST_CHECK_EQUAL(hess.size(),1);
828 CHECK_CLOSE(hess[0],2);
829}
830
831BOOST_AUTO_TEST_CASE( test_grad_reverse13)
832{
833 vector<Node*> list;
834 VNode* x1 = create_var_node();
835 PNode* p1 = create_param_node(0.090901);
836 VNode* x2 = create_var_node();
837 PNode* p2 = create_param_node(0.090901);
838 list.push_back(x1);
839 list.push_back(x2);
840 Node* op1 = create_binary_op_node(OP_TIMES,x1,p1);
841 Node* op2 = create_binary_op_node(OP_TIMES,x2,p2);
842 Node* root = create_binary_op_node(OP_PLUS,op1,op2);
843 x1->val = 1;
844 x2->val = 1;
845 vector<double> grad;
846 grad_reverse(root,list,grad);
847 BOOST_CHECK_EQUAL(grad.size(),2);
848 CHECK_CLOSE(grad[0],0.090901);
849 CHECK_CLOSE(grad[1],0.090901);
850}
851
852BOOST_AUTO_TEST_SUITE_END()
853

source code of boost/libs/yap/example/autodiff_example.cpp