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
|
// 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_min_linear.cpp}
abs_min_linear: Example and Test
################################
Purpose
*******
The function
:math:`f : \B{R}^3 \rightarrow \B{R}` defined by
.. math::
:nowrap:
\begin{eqnarray}
f( x_0, x_1 )
& = &
| d_0 - x_0 | + | d_1 - x_0 | + | d_2 - x_0 |
\\
& + &
| d_3 - x_1 | + | d_4 - x_1 | + | d_5 - x_1 |
\\
\end{eqnarray}
is affine, except for its absolute value terms.
For this case, the abs_normal approximation should be equal
to the function itself.
In addition, the function is convex and
:ref:`abs_min_linear-name` should find its global minimizer.
The minimizer of this function is
:math:`x_0 = \R{median}( d_0, d_1, d_2 )`
and
:math:`x_1 = \R{median}( d_3, d_4, d_5 )`
Source
******
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end abs_min_linear.cpp}
-------------------------------------------------------------------------------
*/
// BEGIN C++
# include <cppad/cppad.hpp>
# include "abs_min_linear.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_min_linear(void)
{ bool ok = true;
//
using CppAD::AD;
using CppAD::ADFun;
//
typedef CPPAD_TESTVECTOR(size_t) s_vector;
typedef CPPAD_TESTVECTOR(double) d_vector;
typedef CPPAD_TESTVECTOR( AD<double> ) ad_vector;
//
size_t dpx = 3; // number of data points per x variable
size_t level = 0; // level of tracing
size_t n = 2; // size of x
size_t m = 1; // size of y
size_t s = dpx * n; // number of data points and absolute values
// data points
d_vector data(s);
for(size_t i = 0; i < s; i++)
data[i] = double(s - i) + 5.0 - double(i % 2) / 2.0;
//
// 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 );
AD<double> sum = 0.0;
for(size_t j = 0; j < n; j++)
for(size_t k = 0; k < dpx; k++)
sum += abs( data[j * dpx + k] - ad_x[j] );
ad_y[0] = sum;
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(0.0);
// 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);
// trust region bound (make large enough to include solutuion)
d_vector bound(n);
for(size_t j = 0; j < n; j++)
bound[j] = 10.0;
// convergence criteria
d_vector epsilon(2);
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
epsilon[0] = eps99;
epsilon[1] = eps99;
// maximum number of iterations
s_vector maxitr(2);
maxitr[0] = 10; // maximum number of abs_min_linear iterations
maxitr[1] = 35; // maximum number of qp_interior iterations
// minimize the approxiamtion for f, which is equal to f because
// f is affine, except for absolute value terms
d_vector delta_x(n);
ok &= CppAD::abs_min_linear(
level, n, m, s, g_hat, g_jac, bound, epsilon, maxitr, delta_x
);
// number of data points per variable is odd
ok &= dpx % 2 == 1;
// check that the solution is the median of the corresponding data`
for(size_t j = 0; j < n; j++)
{ // data[j * dpx + 0] , ... , data[j * dpx + dpx - 1] corresponds to x[j]
// the median of this data has index j * dpx + dpx / 2
size_t j_median = j * dpx + (dpx / 2);
//
ok &= CppAD::NearEqual( delta_x[j], data[j_median], eps99, eps99 );
}
return ok;
}
// END C++
|
