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
|
// 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 interface2c.cpp}
Interfacing to C: Example and Test
##################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end interface2c.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp> // CppAD utilities
# include <cassert> // assert macro
namespace { // Begin empty namespace
/*
Compute the value of a sum of Gaussians defined by a and evaluated at x
y = sum_{i=1}^n a[3*i] exp( (x - a[3*i+1])^2 / a[3*i+2])^2 )
where the floating point type is a template parameter
*/
template <class Float>
Float sumGauss(const Float &x, const CppAD::vector<Float> &a)
{
// number of components in a
size_t na = a.size();
// number of Gaussians
size_t n = na / 3;
// check the restricitons on na
assert( na == n * 3 );
// declare temporaries used inside of loop
Float ex, arg;
// initialize sum
Float y = 0.;
// loop with respect to Gaussians
size_t i;
for(i = 0; i < n; i++)
{
arg = (x - a[3*i+1]) / a[3*i+2];
ex = exp(-arg * arg);
y += a[3*i] * ex;
}
return y;
}
/*
Create a C function interface that computes both
y = sum_{i=1}^n a[3*i] exp( (x - a[3*i+1])^2 / a[3*i+2])^2 )
and its derivative with respect to the parameter vector a.
*/
extern "C"
void sumGauss(float x, float a[], float *y, float dyda[], size_t na)
{ // Note that any simple vector could replace CppAD::vector;
// for example, std::vector, std::valarray
// check the restrictions on na
assert( na % 3 == 0 ); // mod(na, 3) = 0
// use the shorthand ADfloat for the type CppAD::AD<float>
typedef CppAD::AD<float> ADfloat;
// vector for indpendent variables
CppAD::vector<ADfloat> A(na); // used with template function above
CppAD::vector<float> acopy(na); // used for derivative calculations
// vector for the dependent variables (there is only one)
CppAD::vector<ADfloat> Y(1);
// copy the independent variables from C vector to CppAD vectors
size_t i;
for(i = 0; i < na; i++)
A[i] = acopy[i] = a[i];
// declare that A is the independent variable vector
CppAD::Independent(A);
// value of x as an ADfloat object
ADfloat X = x;
// Evaluate template version of sumGauss with ADfloat as the template
// parameter. Set the independent variable to the resulting value
Y[0] = sumGauss(X, A);
// create the AD function object F : A -> Y
CppAD::ADFun<float> F(A, Y);
// use Value to convert Y[0] to float and return y = F(a)
*y = CppAD::Value(Y[0]);
// evaluate the derivative F'(a)
CppAD::vector<float> J(na);
J = F.Jacobian(acopy);
// return the value of dyda = F'(a) as a C vector
for(i = 0; i < na; i++)
dyda[i] = J[i];
return;
}
/*
Link CppAD::NearEqual so do not have to use namespace notation in Interface2C
*/
bool NearEqual(float x, float y, float r, float a)
{ return CppAD::NearEqual(x, y, r, a);
}
} // End empty namespace
bool Interface2C(void)
{ // This routine is intentionally coded as if it were a C routine
// except for the fact that it uses the predefined type bool.
bool ok = true;
// declare variables
float x, a[6], y, dyda[6], tmp[6];
size_t na, i;
// number of parameters (3 for each Gaussian)
na = 6;
// number of Gaussians: n = na / 3;
// value of x
x = 1.;
// value of the parameter vector a
for(i = 0; i < na; i++)
a[i] = (float) (i+1);
// evaulate function and derivative
sumGauss(x, a, &y, dyda, na);
// compare dyda to central difference approximation for deriative
for(i = 0; i < na; i++)
{ // local variables
float eps, ai, yp, ym, dy_da;
// We assume that the type float has at least 7 digits of
// precision, so we choose eps to be about pow(10., -7./2.).
eps = (float) 3e-4;
// value of this component of a
ai = a[i];
// evaluate F( a + eps * ei )
a[i] = ai + eps;
sumGauss(x, a, &yp, tmp, na);
// evaluate F( a - eps * ei )
a[i] = ai - eps;
sumGauss(x, a, &ym, tmp, na);
// evaluate central difference approximates for partial
dy_da = (yp - ym) / (2 * eps);
// restore this component of a
a[i] = ai;
ok &= NearEqual(dyda[i], dy_da, eps, eps);
}
return ok;
}
// END C++
|
