1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin optimize_forward_active.cpp}
Optimize Forward Activity Analysis: Example and Test
####################################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end optimize_forward_active.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp>
namespace {
struct tape_size { size_t n_var; size_t n_op; };
template <class Vector> void fun(
const Vector& x, Vector& y, tape_size& before, tape_size& after
)
{ typedef typename Vector::value_type scalar;
// phantom variable with index 0 and independent variables
// begin operator, independent variable operators and end operator
before.n_var = 1 + x.size(); before.n_op = 2 + x.size();
after.n_var = 1 + x.size(); after.n_op = 2 + x.size();
// adding the constant zero does not take any operations
scalar zero = 0.0 + x[0];
before.n_var += 0; before.n_op += 0;
after.n_var += 0; after.n_op += 0;
// multiplication by the constant one does not take any operations
scalar one = 1.0 * x[1];
before.n_var += 0; before.n_op += 0;
after.n_var += 0; after.n_op += 0;
// multiplication by the constant zero does not take any operations
// and results in the constant zero.
scalar two = 0.0 * x[0];
// operations that only involve constants do not take any operations
scalar three = (1.0 + two) * 3.0;
before.n_var += 0; before.n_op += 0;
after.n_var += 0; after.n_op += 0;
// The optimizer will reconize that zero + one = one + zero
// for all values of x.
scalar four = zero + one;
scalar five = one + zero;
before.n_var += 2; before.n_op += 2;
after.n_var += 1; after.n_op += 1;
// The optimizer will reconize that sin(x[3]) = sin(x[3])
// for all values of x. Note that, for computation of derivatives,
// sin(x[3]) and cos(x[3]) are stored on the tape as a pair.
scalar six = sin(x[2]);
scalar seven = sin(x[2]);
before.n_var += 4; before.n_op += 2;
after.n_var += 2; after.n_op += 1;
// If we used addition here, five + seven = zero + one + seven
// which would get converted to a cumulative summation operator.
scalar eight = five * seven;
before.n_var += 1; before.n_op += 1;
after.n_var += 1; after.n_op += 1;
// Use two, three, four and six in order to avoid a compiler warning
// Note that addition of two and three does not take any operations.
// Also note that optimizer reconizes four * six == five * seven.
scalar nine = eight + four * six * (two + three);
before.n_var += 3; before.n_op += 3;
after.n_var += 2; after.n_op += 2;
// results for this operation sequence
y[0] = nine;
before.n_var += 0; before.n_op += 0;
after.n_var += 0; after.n_op += 0;
}
}
bool forward_active(void)
{ bool ok = true;
using CppAD::AD;
using CppAD::NearEqual;
double eps10 = 10.0 * std::numeric_limits<double>::epsilon();
// domain space vector
size_t n = 3;
CPPAD_TESTVECTOR(AD<double>) ax(n);
ax[0] = 0.5;
ax[1] = 1.5;
ax[2] = 2.0;
// declare independent variables and start tape recording
CppAD::Independent(ax);
// range space vector
size_t m = 1;
CPPAD_TESTVECTOR(AD<double>) ay(m);
tape_size before, after;
fun(ax, ay, before, after);
// create f: x -> y and stop tape recording
CppAD::ADFun<double> f(ax, ay);
ok &= f.size_order() == 1; // this constructor does 0 order forward
ok &= f.size_var() == before.n_var;
ok &= f.size_op() == before.n_op;
// Optimize the operation sequence
// Note that, for this case, all the optimization was done during
// the recording and there is no benifit to the optimization.
f.optimize();
ok &= f.size_order() == 0; // 0 order forward not present
ok &= f.size_var() == after.n_var;
ok &= f.size_op() == after.n_op;
// check zero order forward with different argument value
CPPAD_TESTVECTOR(double) x(n), y(m), check(m);
for(size_t i = 0; i < n; i++)
x[i] = double(i + 2);
y = f.Forward(0, x);
fun(x, check, before, after);
ok &= NearEqual(y[0], check[0], eps10, eps10);
return ok;
}
// END C++
|
