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
|
// 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 interp_retape.cpp}
{xrst_spell
retaping
}
Interpolation With Retaping: Example and Test
#############################################
See Also
********
:ref:`interp_onetape.cpp-name`
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end interp_retape.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp>
# include <cassert>
# include <cmath>
namespace {
double ArgumentValue[] = {
.0 ,
.2 ,
.4 ,
.8 ,
1.
};
double FunctionValue[] = {
std::sin( ArgumentValue[0] ) ,
std::sin( ArgumentValue[1] ) ,
std::sin( ArgumentValue[2] ) ,
std::sin( ArgumentValue[3] ) ,
std::sin( ArgumentValue[4] )
};
size_t TableLength = 5;
size_t Index(const CppAD::AD<double> &x)
{ // determine the index j such that x is between
// ArgumentValue[j] and ArgumentValue[j+1]
static size_t j = 0;
while ( x < ArgumentValue[j] && j > 0 )
j--;
while ( x > ArgumentValue[j+1] && j < TableLength - 2)
j++;
// assert conditions that must be true given logic above
assert( j >= 0 && j < TableLength - 1 );
return j;
}
double Argument(const CppAD::AD<double> &x)
{ size_t j = Index(x);
return ArgumentValue[j];
}
double Function(const CppAD::AD<double> &x)
{ size_t j = Index(x);
return FunctionValue[j];
}
double Slope(const CppAD::AD<double> &x)
{ size_t j = Index(x);
double dx = ArgumentValue[j+1] - ArgumentValue[j];
double dy = FunctionValue[j+1] - FunctionValue[j];
return dy / dx;
}
}
bool interp_retape(void)
{ bool ok = true;
using CppAD::AD;
using CppAD::NearEqual;
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
// domain space vector
size_t n = 1;
CPPAD_TESTVECTOR(AD<double>) X(n);
// loop over argument values
size_t k;
for(k = 0; k < TableLength - 1; k++)
{
X[0] = .4 * ArgumentValue[k] + .6 * ArgumentValue[k+1];
// declare independent variables and start tape recording
// (use a different tape for each argument value)
CppAD::Independent(X);
// evaluate piecewise linear interpolant at X[0]
AD<double> A = Argument(X[0]);
AD<double> F = Function(X[0]);
AD<double> S = Slope(X[0]);
AD<double> I = F + (X[0] - A) * S;
// range space vector
size_t m = 1;
CPPAD_TESTVECTOR(AD<double>) Y(m);
Y[0] = I;
// create f: X -> Y and stop tape recording
CppAD::ADFun<double> f(X, Y);
// vectors for arguments to the function object f
CPPAD_TESTVECTOR(double) x(n); // argument values
CPPAD_TESTVECTOR(double) y(m); // function values
CPPAD_TESTVECTOR(double) dx(n); // differentials in x space
CPPAD_TESTVECTOR(double) dy(m); // differentials in y space
// to check function value we use the fact that X[0] is between
// ArgumentValue[k] and ArgumentValue[k+1]
double delta, check;
x[0] = Value(X[0]);
delta = ArgumentValue[k+1] - ArgumentValue[k];
check = FunctionValue[k+1] * (x[0]-ArgumentValue[k]) / delta
+ FunctionValue[k] * (ArgumentValue[k+1]-x[0]) / delta;
ok &= NearEqual(Y[0], check, eps99, eps99);
// evaluate partials w.r.t. x[0]
dx[0] = 1.;
dy = f.Forward(1, dx);
// check that the derivative is the slope
check = (FunctionValue[k+1] - FunctionValue[k])
/ (ArgumentValue[k+1] - ArgumentValue[k]);
ok &= NearEqual(dy[0], check, eps99, eps99);
}
return ok;
}
// END C++
|
