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
134
135
136
137
138
139
140
141
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-24 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin abs_get_started.cpp}
{xrst_spell
rclrcl
}
abs_normal Getting Started: Example and Test
############################################
Purpose
*******
Creates an :ref:`abs_normal<abs_normal_fun-name>`
representation :math:`g` for the function
:math:`f : \B{R}^3 \rightarrow \B{R}` defined by
.. math::
f( x_0, x_1, x_2 ) = | x_0 + x_1 | + | x_1 + x_2 |
The corresponding
:ref:`abs_normal_fun@g` :math:`: \B{R}^5 \rightarrow \B{R}^3` is
given by
.. math::
\begin{array}{rclrcl}
g_0 ( x_0, x_1, x_2, u_0, u_1 ) & = & u_0 + u_1 & = & y_0 (x, u)
\\
g_1 ( x_0, x_1, x_2, u_0, u_1 ) & = & x_0 + x_1 & = & z_0 (x, u)
\\
g_1 ( x_0, x_1, x_2, u_0, u_1 ) & = & x_1 + x_2 & = & z_1 (x, u)
\end{array}
Source
******
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end abs_get_started.cpp}
-------------------------------------------------------------------------------
*/
// BEGIN C++
# include <cppad/cppad.hpp>
namespace {
CPPAD_TESTVECTOR(double) join(
const CPPAD_TESTVECTOR(double)& x ,
const CPPAD_TESTVECTOR(double)& u )
{ size_t n = x.size();
size_t s = u.size();
CPPAD_TESTVECTOR(double) xu(n + s);
for(size_t j = 0; j < n; j++)
xu[j] = x[j];
for(size_t j = 0; j < s; j++)
xu[n + j] = u[j];
return xu;
}
}
bool get_started(void)
{ bool ok = true;
//
using CppAD::AD;
using CppAD::ADFun;
//
size_t n = 3; // size of x
size_t m = 1; // size of y
size_t s = 2; // size of u and z
//
// record the function f(x)
CPPAD_TESTVECTOR( AD<double> ) ax(n), ay(m);
for(size_t j = 0; j < n; j++)
ax[j] = double(j + 1);
Independent( ax );
// for this example, we ensure first absolute value is | x_0 + x_1 |
AD<double> a0 = abs( ax[0] + ax[1] );
// and second absolute value is | x_1 + x_2 |
AD<double> a1 = abs( ax[1] + ax[2] );
ay[0] = a0 + a1;
ADFun<double> f(ax, ay);
// create its abs_normal representation in g, a
ADFun<double> g, a;
f.abs_normal_fun(g, a);
// check dimension of domain and range space for g
ok &= g.Domain() == n + s;
ok &= g.Range() == m + s;
// check dimension of domain and range space for a
ok &= a.Domain() == n;
ok &= a.Range() == s;
// --------------------------------------------------------------------
// a(x) has all the operations used to compute f(x), but the sum of the
// absolute values is not needed for a(x), so optimize it out.
size_t n_op = f.size_op();
ok &= a.size_op() == n_op;
a.optimize();
ok &= a.size_op() < n_op;
// --------------------------------------------------------------------
// zero order forward mode calculation using g(x, u)
CPPAD_TESTVECTOR(double) x(n), u(s), xu(n+s), yz(m+s);
for(size_t j = 0; j < n; j++)
x[j] = double(j + 2);
for(size_t j = 0; j < s; j++)
u[j] = double(j + n + 2);
xu = join(x, u);
yz = g.Forward(0, xu);
// check y_0(x, u)
double y0 = u[0] + u[1];
ok &= y0 == yz[0];
// check z_0 (x, u)
double z0 = x[0] + x[1];
ok &= z0 == yz[1];
// check z_1 (x, u)
double z1 = x[1] + x[2];
ok &= z1 == yz[2];
// --------------------------------------------------------------------
// check that y(x, a(x) ) == f(x)
CPPAD_TESTVECTOR(double) y(m);
y = f.Forward(0, x); // y = f(x)
u = a.Forward(0, x); // u = a(x)
xu = join(x, u); // xu = ( x, a(x) )
yz = g.Forward(0, xu); // yz = ( y(x, a(x)), z(x, a(x)) )
ok &= yz[0] == y[0];
return ok;
}
// END C++
|
