# 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;
}