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
|
// 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 abs_eval.cpp}
abs_eval: Example and Test
##########################
Purpose
*******
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 |
is affine, except for its absolute value terms.
For this case, the abs_normal approximation should be equal
to the function itself.
Source
******
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end abs_eval.cpp}
-------------------------------------------------------------------------------
*/
// BEGIN C++
# include <cppad/cppad.hpp>
# include "abs_eval.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 abs_eval(void)
{ bool ok = true;
//
using CppAD::AD;
using CppAD::ADFun;
//
typedef CPPAD_TESTVECTOR(double) d_vector;
typedef CPPAD_TESTVECTOR( AD<double> ) ad_vector;
//
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
//
size_t n = 3; // size of x
size_t m = 1; // size of y
size_t s = 2; // number of absolute value terms
//
// record the function f(x)
ad_vector ad_x(n), ad_y(m);
for(size_t j = 0; j < n; j++)
ad_x[j] = double(j + 1);
Independent( ad_x );
// for this example, we ensure first absolute value is | x_0 + x_1 |
AD<double> ad_0 = abs( ad_x[0] + ad_x[1] );
// and second absolute value is | x_1 + x_2 |
AD<double> ad_1 = abs( ad_x[1] + ad_x[2] );
ad_y[0] = ad_0 + ad_1;
ADFun<double> f(ad_x, ad_y);
// 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;
// --------------------------------------------------------------------
// Choose a point x_hat
d_vector x_hat(n);
for(size_t j = 0; j < n; j++)
x_hat[j] = double(j - 1);
// value of a_hat = a(x_hat)
d_vector a_hat = a.Forward(0, x_hat);
// (x_hat, a_hat)
d_vector xu_hat = join(x_hat, a_hat);
// value of g[ x_hat, a_hat ]
d_vector g_hat = g.Forward(0, xu_hat);
// Jacobian of g[ x_hat, a_hat ]
d_vector g_jac = g.Jacobian(xu_hat);
// value of delta_x
d_vector delta_x(n);
delta_x[0] = 1.0;
delta_x[1] = -2.0;
delta_x[2] = +2.0;
// value of x
d_vector x(n);
for(size_t j = 0; j < n; j++)
x[j] = x_hat[j] + delta_x[j];
// value of f(x)
d_vector y = f.Forward(0, x);
// value of g_tilde
d_vector g_tilde = CppAD::abs_eval(n, m, s, g_hat, g_jac, delta_x);
// should be equal because f is affine, except for abs terms
ok &= CppAD::NearEqual(y[0], g_tilde[0], eps99, eps99);
return ok;
}
// END C++
|
