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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
|
# ifndef CPPAD_EXAMPLE_ABS_NORMAL_ABS_EVAL_HPP
# define CPPAD_EXAMPLE_ABS_NORMAL_ABS_EVAL_HPP
// 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_eval}
abs_normal: Evaluate First Order Approximation
##############################################
Syntax
******
| *g_tilde* = ``abs_eval`` ( *n* , *m* , *s* , *g_hat* , *g_jac* , *delta_x* )
Prototype
*********
{xrst_literal
// BEGIN PROTOTYPE
// END PROTOTYPE
}
Source
******
This following is a link to the source code for this example:
:ref:`abs_eval.hpp-name` .
Purpose
*******
Given a current that abs-normal representation at a point
:math:`\hat{x} \in \B{R}^n`,
and a :math:`\Delta x \in \B{R}^n`,
this routine evaluates the abs-normal
:ref:`approximation for f(x)<abs_normal_fun@Abs-normal Approximation@Approximating f(x)>`
where :math:`x = \hat{x} + \Delta x`.
Vector
******
The type *Vector* is a
simple vector with elements of type ``double`` .
f
*
We use the notation *f* for the original function; see
:ref:`abs_normal_fun@f` .
n
*
This is the dimension of the domain space for *f* ; see
:ref:`abs_normal_fun@f@n` .
m
*
This is the dimension of the range space for *f* ; see
:ref:`abs_normal_fun@f@m` .
s
*
This is the number of absolute value terms in *f* ; see
g
*
We use the notation *g* for the abs-normal representation of *f* ;
see :ref:`abs_normal_fun@g` .
g_hat
*****
This vector has size *m* + *s* and is the value of
*g* ( *x* , *u* ) at :math:`x = \hat{x}` and :math:`u = a( \hat{x} )`.
g_jac
*****
This vector has size ( *m* + *s* ) * ( *n* + *s* ) and is the Jacobian of
:math:`g(x, u)` at :math:`x = \hat{x}` and :math:`u = a( \hat{x} )`.
delta_x
*******
This vector has size *n* and is the difference
:math:`\Delta x = x - \hat{x}`,
where :math:`x` is the point that we are approximating :math:`f(x)`.
g_tilde
*******
This vector has size *m* + *s* and is a the
first order approximation for
:ref:`abs_normal_fun@g`
that corresponds to the point
:math:`x = \hat{x} + \Delta x` and :math:`u = a(x)`.
{xrst_toc_hidden
example/abs_normal/abs_eval.cpp
example/abs_normal/abs_eval.xrst
}
Example
*******
The file :ref:`abs_eval.cpp-name` contains an example and test of
``abs_eval`` .
{xrst_end abs_eval}
-----------------------------------------------------------------------------
*/
// BEGIN C++
namespace CppAD { // BEGIN_CPPAD_NAMESPACE
// BEGIN PROTOTYPE
template <class Vector>
Vector abs_eval(
size_t n ,
size_t m ,
size_t s ,
const Vector& g_hat ,
const Vector& g_jac ,
const Vector& delta_x )
// END PROTOTYPE
{ using std::fabs;
//
CPPAD_ASSERT_KNOWN(
size_t(delta_x.size()) == n,
"abs_eval: size of delta_x not equal to n"
);
CPPAD_ASSERT_KNOWN(
size_t(g_hat.size()) == m + s,
"abs_eval: size of g_hat not equal to m + s"
);
CPPAD_ASSERT_KNOWN(
size_t(g_jac.size()) == (m + s) * (n + s),
"abs_eval: size of g_jac not equal to (m + s)*(n + s)"
);
# ifndef NDEBUG
// Check that partial_u z(x, u) is strictly lower triangular
for(size_t i = 0; i < s; i++)
{ for(size_t j = i; j < s; j++)
{ // index in g_jac of partial of z_i w.r.t u_j
// (note that g_jac has n + s elements in each row)
size_t index = (m + i) * (n + s) + (n + j);
CPPAD_ASSERT_KNOWN(
g_jac[index] == 0.0,
"abs_eval: partial z_i w.r.t u_j non-zero for i <= j"
);
}
}
# endif
// return value
Vector g_tilde(m + s);
//
// compute z_tilde, the last s components of g_tilde
for(size_t i = 0; i < s; i++)
{ // start at z_hat_i
g_tilde[m + i] = g_hat[m + i];
// contribution for change x
for(size_t j = 0; j < n; j++)
{ // index in g_jac of partial of z_i w.r.t x_j
size_t index = (m + i) * (n + s) + j;
// add contribution for delta_x_j to z_tilde_i
g_tilde[m + i] += g_jac[index] * delta_x[j];
}
// contribution for change in u_j for j < i
for(size_t j = 0; j < i; j++)
{ // approixmation for change in absolute value
double delta_a_j = fabs(g_tilde[m + j]) - fabs(g_hat[m + j]);
// index in g_jac of partial of z_i w.r.t u_j
size_t index = (m + i) * (n + s) + n + j;
// add contribution for delta_a_j to s_tilde_i
g_tilde[m + i] += g_jac[index] * delta_a_j;
}
}
//
// compute y_tilde, the first m components of g_tilde
for(size_t i = 0; i < m; i++)
{ // start at y_hat_i
g_tilde[i] = g_hat[i];
// contribution for change x
for(size_t j = 0; j < n; j++)
{ // index in g_jac of partial of y_i w.r.t x_j
size_t index = i * (n + s) + j;
// add contribution for delta_x_j to y_tilde_i
g_tilde[i] += g_jac[index] * delta_x[j];
}
// contribution for change in u_j
for(size_t j = 0; j < s; j++)
{ // approximation for change in absolute value
double delta_a_j = fabs(g_tilde[m + j]) - fabs(g_hat[m + j]);
// index in g_jac of partial of y_i w.r.t u_j
size_t index = i * (n + s) + n + j;
// add contribution for delta_a_j to s_tilde_i
g_tilde[i] += g_jac[index] * delta_a_j;
}
}
return g_tilde;
}
} // END_CPPAD_NAMESPACE
// END C++
# endif
|
