$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
# 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 small, ai, yp, ym, dy_da; // We assume that the type float has at least 7 digits of // precision, so we choose small to be about pow(10., -7./2.). small = (float) 3e-4; // value of this component of a ai = a[i]; // evaluate F( a + small * ei ) a[i] = ai + small; sumGauss(x, a, &yp, tmp, na); // evaluate F( a - small * ei ) a[i] = ai - small; sumGauss(x, a, &ym, tmp, na); // evaluate central difference approximates for partial dy_da = (yp - ym) / (2 * small); // restore this component of a a[i] = ai; ok &= NearEqual(dyda[i], dy_da, small, small); } return ok; }