gprof -b speed/profile/profile.exe

 
	gprof -b speed/profile/profile | sed -f speed/profile/gprof.sed

 
	gprof -b speed/profile/profile.exe | sed -f speed/profile/gprof.sed

 
	./configure --prefix=/usr/local

 
	std::vector

 
	./configure --prefix=/usr/local POSTFIX_DIR=coin

 
	boost::numeric::ublas::vector

 
	make install

 
	which svn

 
	which omhelp

 
	cppad-20100902\introduction\get_started\get_started.sln

 
	introduction\get_started\Debug\get_started

 
	cppad-20100902\introduction\exp_apx\exp_apx.sln

 
	introduction\exp_apx\Debug\exp_apx

 
	cppad-20100902\example\example.sln

 
	example\Debug\example

 
	cppad-20100902\test_more\test_more.sln

 
	test_more\Debug\test_more

 
	cppad-20100902\print_for\print_for.sln

 
	print_for\Debug\print_for

 
	cppad-20100902\speed\cppad\cppad.sln

 
	speed\cppad\Debug\cppad correct 123
	speed\cppad\Debug\cppad speed 123

 
	cppad-20100902\speed\double\double.sln

 
	speed\double\Debug\double correct 123
	speed\double\Debug\double speed 123

 
	speed\example\Debug\example

 
#include <iostream>      // standard input/output 
#include <vector>        // standard vector
#include <cppad/cppad.hpp> // the CppAD package http://www.coin-or.org/CppAD/

namespace { 
      // define y(x) = Poly(a, x) in the empty namespace
      template <class Type>
      Type Poly(const std::vector<double> &a, const Type &x)
      {     size_t k  = a.size();
            Type y   = 0.;  // initialize summation
            Type x_i = 1.;  // initialize x^i
            size_t i;
            for(i = 0; i < k; i++)
            {     y   += a[i] * x_i;  // y   = y + a_i * x^i
                  x_i *= x;           // x_i = x_i * x
            }
            return y;
      }
}
// main program
int main(void)
{     using CppAD::AD;           // use AD as abbreviation for CppAD::AD
      using std::vector;         // use vector as abbreviation for std::vector
      size_t i;                  // a temporary index

      // vector of polynomial coefficients
      size_t k = 5;              // number of polynomial coefficients
      vector<double> a(k);       // vector of polynomial coefficients
      for(i = 0; i < k; i++)
            a[i] = 1.;           // value of polynomial coefficients

      // domain space vector
      size_t n = 1;              // number of domain space variables
      vector< AD<double> > X(n); // vector of domain space variables
      X[0] = 3.;                 // value corresponding to operation sequence

      // declare independent variables and start recording operation sequence
      CppAD::Independent(X);

      // range space vector
      size_t m = 1;              // number of ranges space variables
      vector< AD<double> > Y(m); // vector of ranges space variables
      Y[0] = Poly(a, X[0]);      // value during recording of operations

      // store operation sequence in f: X -> Y and stop recording
      CppAD::ADFun<double> f(X, Y);

      // compute derivative using operation sequence stored in f
      vector<double> jac(m * n); // Jacobian of f (m by n matrix)
      vector<double> x(n);       // domain space vector
      x[0] = 3.;                 // argument value for derivative
      jac  = f.Jacobian(x);      // Jacobian for operation sequence

      // print the results
      std::cout << "f'(3) computed by CppAD = " << jac[0] << std::endl;

      // check if the derivative is correct
      int error_code;
      if( jac[0] == 142. )
            error_code = 0;      // return code for correct case
      else  error_code = 1;      // return code for incorrect case

      return error_code;
}

 
	f'(3) computed by CppAD = 142

 
template <class Type>
Type exp_2(const Type &x) 
{       Type v1  = x;                // v1 = x
        Type v2  = Type(1) + v1;     // v2 = 1 + x
        Type v3  = v1 * v1;          // v3 = x^2
        Type v4  = v3 / Type(2);     // v4 = x^2 / 2 
        Type v5  = v2 + v4;          // v5 = 1 + x + x^2 / 2
        return v5;                   // exp_2(x) = 1 + x + x^2 / 2
}

 
# include <cmath>           // define fabs function
# include "exp_2.hpp"       // definition of exp_2 algorithm
bool exp_2(void)
{	double x     = .5;
	double check = 1 + x + x * x / 2.; 
	bool   ok    = std::fabs( exp_2(x) - check ) <= 1e-10; 
	return ok;
}

 # include <cmath>            // for fabs function
bool exp_2_for0(double *v0)  // double v0[6]
{	bool  ok = true;
	double x = .5;

	v0[1] = x;                                  // v1 = x
	ok  &= std::fabs( v0[1] - 0.5) < 1e-10;

	v0[2] = 1. + v0[1];                         // v2 = 1 + v1
	ok  &= std::fabs( v0[2] - 1.5) < 1e-10;

	v0[3] = v0[1] * v0[1];                      // v3 = v1 * v1
	ok  &= std::fabs( v0[3] - 0.25) < 1e-10;

	v0[4] = v0[3] / 2.;                         // v4 = v3 / 2
	ok  &= std::fabs( v0[4] - 0.125) < 1e-10;

	v0[5] = v0[2] + v0[4];                      // v5  = v2 + v4
	ok  &= std::fabs( v0[5] - 1.625) < 1e-10;

	return ok;
}
bool exp_2_for0(void)
{	double v0[6];
	return exp_2_for0(v0);
}

 # include <cmath>                   // prototype for fabs
extern bool exp_2_for0(double *v0); // computes zero order forward sweep
bool exp_2_for1(double *v1)         // double v1[6]
{	bool ok = true;
	double v0[6];

	// set the value of v0[j] for j = 1 , ... , 5
	ok &= exp_2_for0(v0);

	v1[1] = 1.;                                     // v1 = x
	ok    &= std::fabs( v1[1] - 1. ) <= 1e-10;

	v1[2] = v1[1];                                  // v2 = 1 + v1
	ok    &= std::fabs( v1[2] - 1. ) <= 1e-10;

	v1[3] = v1[1] * v0[1] + v0[1] * v1[1];          // v3 = v1 * v1
	ok    &= std::fabs( v1[3] - 1. ) <= 1e-10;

	v1[4] = v1[3] / 2.;                             // v4 = v3 / 2
	ok    &= std::fabs( v1[4] - 0.5) <= 1e-10;

	v1[5] = v1[2] + v1[4];                          // v5 = v2 + v4
	ok    &= std::fabs( v1[5] - 1.5) <= 1e-10;

	return ok;
}
bool exp_2_for1(void)
{	double v1[6];
	return exp_2_for1(v1);
}

 # include <cstddef>                 // define size_t
# include <cmath>                   // prototype for fabs
extern bool exp_2_for0(double *v0); // computes zero order forward sweep
bool exp_2_rev1(void)
{	bool ok = true;

	// set the value of v0[j] for j = 1 , ... , 5
	double v0[6];
	ok &= exp_2_for0(v0);

	// initial all partial derivatives as zero
	double f_v[6];
	size_t j;
	for(j = 0; j < 6; j++)
		f_v[j] = 0.;

	// set partial derivative for f5
	f_v[5] = 1.;
	ok &= std::fabs( f_v[5] - 1. ) <= 1e-10; // f5_v5

	// f4 = f5( v1 , v2 , v3 , v4 , v2 + v4 )
	f_v[2] += f_v[5] * 1.;
	f_v[4] += f_v[5] * 1.;
	ok &= std::fabs( f_v[2] - 1. ) <= 1e-10; // f4_v2
	ok &= std::fabs( f_v[4] - 1. ) <= 1e-10; // f4_v4

	// f3 = f4( v1 , v2 , v3 , v3 / 2 )
	f_v[3] += f_v[4] / 2.;
	ok &= std::fabs( f_v[3] - 0.5) <= 1e-10; // f3_v3

	// f2 = f3( v1 , v2 , v1 * v1 )
	f_v[1] += f_v[3] * 2. * v0[1];
	ok &= std::fabs( f_v[1] - 0.5) <= 1e-10; // f2_v1

	// f1 = f2( v1 , 1 + v1 )
	f_v[1] += f_v[2] * 1.;
	ok &= std::fabs( f_v[1] - 1.5) <= 1e-10; // f1_v1

	return ok;
}

 # include <cmath>                   // prototype for fabs
extern bool exp_2_for0(double *v0); // computes zero order forward sweep
extern bool exp_2_for1(double *v1); // computes first order forward sweep
bool exp_2_for2(void)
{	bool ok = true;
	double v0[6], v1[6], v2[6];

	// set the value of v0[j], v1[j], for j = 1 , ... , 5
	ok &= exp_2_for0(v0);
	ok &= exp_2_for1(v1);

	v2[1] = 0.;                                     // v1 = x
	ok    &= std::fabs( v2[1] - 0. ) <= 1e-10;

	v2[2] = v2[1];                                  // v2 = 1 + v1
	ok    &= std::fabs( v2[2] - 0. ) <= 1e-10;

	v2[3] = 2.*(v0[1]*v2[1] + v1[1]*v1[1]);         // v3 = v1 * v1
	ok    &= std::fabs( v2[3] - 2. ) <= 1e-10;

	v2[4] = v2[3] / 2.;                             // v4 = v3 / 2
	ok    &= std::fabs( v2[4] - 1. ) <= 1e-10;

	v2[5] = v2[2] + v2[4];                          // v5 = v2 + v4
	ok    &= std::fabs( v2[5] - 1. ) <= 1e-10;

	return ok;
}

 # include <cstddef>                 // define size_t
# include <cmath>                   // prototype for fabs
extern bool exp_2_for0(double *v0); // computes zero order forward sweep
extern bool exp_2_for1(double *v1); // computes first order forward sweep
bool exp_2_rev2(void)
{	bool ok = true;

	// set the value of v0[j], v1[j] for j = 1 , ... , 5
	double v0[6], v1[6];
	ok &= exp_2_for0(v0);
	ok &= exp_2_for1(v1);

	// initial all partial derivatives as zero
	double f_v0[6], f_v1[6];
	size_t j;
	for(j = 0; j < 6; j++)
	{	f_v0[j] = 0.;
		f_v1[j] = 0.;
	}

	// set partial derivative for f_5
	f_v1[5] = 1.;
	ok &= std::fabs( f_v1[5] - 1. ) <= 1e-10; // partial f_5 w.r.t v_5^1

	// f_4 = f_5( v_1^0 , ... , v_4^1 , v_2^0 + v_4^0 , v_2^1 + v_4^1 )
	f_v0[2] += f_v0[5] * 1.;
	f_v0[4] += f_v0[5] * 1.;
	f_v1[2] += f_v1[5] * 1.;
	f_v1[4] += f_v1[5] * 1.;
	ok &= std::fabs( f_v0[2] - 0. ) <= 1e-10; // partial f_4 w.r.t. v_2^0
	ok &= std::fabs( f_v0[4] - 0. ) <= 1e-10; // partial f_4 w.r.t. v_4^0
	ok &= std::fabs( f_v1[2] - 1. ) <= 1e-10; // partial f_4 w.r.t. v_2^1
	ok &= std::fabs( f_v1[4] - 1. ) <= 1e-10; // partial f_4 w.r.t. v_4^1

	// f_3 = f_4( v_1^0 , ... , v_3^1, v_3^0 / 2 , v_3^1 / 2 )
	f_v0[3] += f_v0[4] / 2.;
	f_v1[3] += f_v1[4] / 2.;
	ok &= std::fabs( f_v0[3] - 0.  ) <= 1e-10; // partial f_3 w.r.t. v_3^0
	ok &= std::fabs( f_v1[3] - 0.5 ) <= 1e-10; // partial f_3 w.r.t. v_3^1

	// f_2 = f_3(  v_1^0 , ... , v_2^1, v_1^0 * v_1^0 , 2 * v_1^0 * v_1^1 )
	f_v0[1] += f_v0[3] * 2. * v0[1];
	f_v0[1] += f_v1[3] * 2. * v1[1];
	f_v1[1] += f_v1[3] * 2. * v0[1];
	ok &= std::fabs( f_v0[1] - 1.  ) <= 1e-10; // partial f_2 w.r.t. v_1^0
	ok &= std::fabs( f_v1[1] - 0.5 ) <= 1e-10; // partial f_2 w.r.t. v_1^1

	// f_1 = f_2( v_1^0 , v_1^1 , 1 + v_1^0 , v_1^1 )
	f_v0[1] += f_v0[2] * 1.;
	f_v1[1] += f_v1[2] * 1.;
	ok &= std::fabs( f_v0[1] - 1. ) <= 1e-10; // partial f_1 w.r.t. v_1^0
	ok &= std::fabs( f_v1[1] - 1.5) <= 1e-10; // partial f_1 w.r.t. v_1^1

	return ok;
}

 

# include <cppad/cppad.hpp>  // http://www.coin-or.org/CppAD/ 
# include "exp_2.hpp"        // second order exponential approximation
bool exp_2_cppad(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::vector;    // can use any simple vector template class
	using CppAD::NearEqual; // checks if values are nearly equal

	// domain space vector
	size_t n = 1; // dimension of the domain space
	vector< AD<double> > X(n);
	X[0] = .5;    // value of x for this operation sequence

	// declare independent variables and start recording operation sequence
	CppAD::Independent(X);

	// evaluate our exponential approximation
	AD<double> x   = X[0];
	AD<double> apx = exp_2(x);  

	// range space vector
	size_t m = 1;  // dimension of the range space
	vector< AD<double> > Y(m);
	Y[0] = apx;    // variable that represents only range space component

	// Create f: X -> Y corresponding to this operation sequence
	// and stop recording. This also executes a zero order forward 
	// sweep using values in X for x.
	CppAD::ADFun<double> f(X, Y);

	// first order forward sweep that computes
	// partial of exp_2(x) with respect to x
	vector<double> dx(n);  // differential in domain space
	vector<double> dy(m);  // differential in range space
	dx[0] = 1.;            // direction for partial derivative
	dy    = f.Forward(1, dx);
	double check = 1.5;
	ok   &= NearEqual(dy[0], check, 1e-10, 1e-10);

	// first order reverse sweep that computes the derivative
	vector<double>  w(m);   // weights for components of the range
	vector<double> dw(n);   // derivative of the weighted function
	w[0] = 1.;              // there is only one weight
	dw   = f.Reverse(1, w); // derivative of w[0] * exp_2(x)
	check = 1.5;            // partial of exp_2(x) with respect to x
	ok   &= NearEqual(dw[0], check, 1e-10, 1e-10);

	// second order forward sweep that computes
	// second partial of exp_2(x) with respect to x
	vector<double> x2(n);     // second order Taylor coefficients 
	vector<double> y2(m);  
	x2[0] = 0.;               // evaluate second partial .w.r.t. x
	y2    = f.Forward(2, x2);
	check = 0.5 * 1.;         // Taylor coef is 1/2 second derivative 
	ok   &= NearEqual(y2[0], check, 1e-10, 1e-10);

	// second order reverse sweep that computes
	// derivative of partial of exp_2(x) w.r.t. x
	dw.resize(2 * n);         // space for first and second derivatives
	dw    = f.Reverse(2, w);
	check = 1.;               // result should be second derivative
	ok   &= NearEqual(dw[0*2+1], check, 1e-10, 1e-10);

	return ok;
}

 
template <class Type>
Type exp_eps(const Type &x, const Type &epsilon)
{	// abs_x = |x|
	Type abs_x = x;
	if( Type(0) > x )
		abs_x = - x;
	// initialize
	int  k    = 0;          // initial order 
	Type term = 1.;         // term = |x|^k / k !
	Type sum  = term;       // initial sum
	while(term > epsilon)
	{	k         = k + 1;          // order for next term
		Type temp = term * abs_x;   // term = |x|^k / (k-1)!
		term      = temp / Type(k); // term = |x|^k / k !
		sum       = sum + term;     // sum  = 1 + ... + |x|^k / k !
	}
	// In the case where x is negative, use exp(x) = 1 / exp(-|x|)
	if( Type(0) > x ) 
		sum = Type(1) / sum;
	return sum;
}

 
# include <cmath>             // for fabs function
# include "exp_eps.hpp"       // definition of exp_eps algorithm
bool exp_eps(void)
{	double x       = .5;
	double epsilon = .2;
	double check   = 1 + .5 + .125; // include 1 term less than epsilon
	bool   ok      = std::fabs( exp_eps(x, epsilon) - check ) <= 1e-10; 
	return ok;
}

 
	if( Type(0) > x ) 
	while(term > epsilon)

 # include <cmath>                // for fabs function
bool exp_eps_for0(double *v0)    // double v0[8]
{	bool  ok = true;
	double x = .5;

	v0[1] = x;                                  // abs_x = x;
	ok  &= std::fabs( v0[1] - 0.5) < 1e-10;

	v0[2] = 1. * v0[1];                         // temp = term * abs_x;
	ok  &= std::fabs( v0[2] - 0.5) < 1e-10;

	v0[3] = v0[2] / 1.;                         // term = temp / Type(k);
	ok  &= std::fabs( v0[3] - 0.5) < 1e-10;

	v0[4] = 1. + v0[3];                         // sum = sum + term;
	ok  &= std::fabs( v0[4] - 1.5) < 1e-10;

	v0[5] = v0[3] * v0[1];                      // temp = term * abs_x;
	ok  &= std::fabs( v0[5] - 0.25) < 1e-10;

	v0[6] = v0[5] / 2.;                         // term = temp / Type(k);
	ok  &= std::fabs( v0[6] - 0.125) < 1e-10;

	v0[7] = v0[4] + v0[6];                      // sum = sum + term;
	ok  &= std::fabs( v0[7] - 1.625) < 1e-10;

	return ok;
}
bool exp_eps_for0(void)
{	double v0[8];
	return exp_eps_for0(v0);
}

 # include <cmath>                     // for fabs function
extern bool exp_eps_for0(double *v0); // computes zero order forward sweep
bool exp_eps_for1(double *v1)         // double v[8]
{	bool ok = true;
	double v0[8];

	// set the value of v0[j] for j = 1 , ... , 7
	ok &= exp_eps_for0(v0);

	v1[1] = 1.;                                      // v1 = x
	ok    &= std::fabs( v1[1] - 1. ) <= 1e-10;

	v1[2] = 1. * v1[1];                              // v2 = 1 * v1
	ok    &= std::fabs( v1[2] - 1. ) <= 1e-10;

	v1[3] = v1[2] / 1.;                              // v3 = v2 / 1
	ok    &= std::fabs( v1[3] - 1. ) <= 1e-10;

	v1[4] = v1[3];                                   // v4 = 1 + v3
	ok    &= std::fabs( v1[4] - 1. ) <= 1e-10;

	v1[5] = v1[3] * v0[1] + v0[3] * v1[1];           // v5 = v3 * v1
	ok    &= std::fabs( v1[5] - 1. ) <= 1e-10;

	v1[6] = v1[5] / 2.;                              // v6 = v5 / 2
	ok    &= std::fabs( v1[6] - 0.5 ) <= 1e-10;

	v1[7] = v1[4] + v1[6];                           // v7 = v4 + v6
	ok    &= std::fabs( v1[7] - 1.5 ) <= 1e-10;

	return ok;
}
bool exp_eps_for1(void)
{	double v1[8];
	return exp_eps_for1(v1);
}

 # include <cstddef>                     // define size_t
# include <cmath>                       // for fabs function
extern bool exp_eps_for0(double *v0);   // computes zero order forward sweep
bool exp_eps_rev1(void)
{	bool ok = true;

	// set the value of v0[j] for j = 1 , ... , 7
	double v0[8];
	ok &= exp_eps_for0(v0);

	// initial all partial derivatives as zero
	double f_v[8];
	size_t j;
	for(j = 0; j < 8; j++)
		f_v[j] = 0.;

	// set partial derivative for f7
	f_v[7] = 1.;
	ok    &= std::fabs( f_v[7] - 1. ) <= 1e-10;     // f7_v7

	// f6( v1 , v2 , v3 , v4 , v5 , v6 )
	f_v[4] += f_v[7] * 1.;
	f_v[6] += f_v[7] * 1.;
	ok     &= std::fabs( f_v[4] - 1.  ) <= 1e-10;   // f6_v4
	ok     &= std::fabs( f_v[6] - 1.  ) <= 1e-10;   // f6_v6

	// f5( v1 , v2 , v3 , v4 , v5 )
	f_v[5] += f_v[6] / 2.;
	ok     &= std::fabs( f_v[5] - 0.5 ) <= 1e-10;   // f5_v5

	// f4( v1 , v2 , v3 , v4 )
	f_v[1] += f_v[5] * v0[3];
	f_v[3] += f_v[5] * v0[1];
	ok     &= std::fabs( f_v[1] - 0.25) <= 1e-10;   // f4_v1
	ok     &= std::fabs( f_v[3] - 0.25) <= 1e-10;   // f4_v3

	// f3( v1 , v2 , v3 )
	f_v[3] += f_v[4] * 1.;
	ok     &= std::fabs( f_v[3] - 1.25) <= 1e-10;   // f3_v3

	// f2( v1 , v2 )
	f_v[2] += f_v[3] / 1.;
	ok     &= std::fabs( f_v[2] - 1.25) <= 1e-10;   // f2_v2

	// f1( v1 )
	f_v[1] += f_v[2] * 1.;
	ok     &= std::fabs( f_v[1] - 1.5 ) <= 1e-10;   // f1_v2

	return ok;
}

 # include <cmath>                     // for fabs function
extern bool exp_eps_for0(double *v0); // computes zero order forward sweep
extern bool exp_eps_for1(double *v1); // computes first order forward sweep
bool exp_eps_for2(void)
{	bool ok = true;
	double v0[8], v1[8], v2[8];

	// set the value of v0[j], v1[j] for j = 1 , ... , 7
	ok &= exp_eps_for0(v0);
	ok &= exp_eps_for1(v1);

	v2[1] = 0.;                                      // v1 = x
	ok    &= std::fabs( v2[1] - 0. ) <= 1e-10;

	v2[2] = 1. * v2[1];                              // v2 = 1 * v1
	ok    &= std::fabs( v2[2] - 0. ) <= 1e-10;

	v2[3] = v2[2] / 1.;                              // v3 = v2 / 1
	ok    &= std::fabs( v2[3] - 0. ) <= 1e-10;

	v2[4] = v2[3];                                   // v4 = 1 + v3
	ok    &= std::fabs( v2[4] - 0. ) <= 1e-10;

	v2[5] = v2[3] * v0[1] + 2. * v1[3] * v1[1]       // v5 = v3 * v1
	      + v0[3] * v2[1];           
	ok    &= std::fabs( v2[5] - 2. ) <= 1e-10;

	v2[6] = v2[5] / 2.;                              // v6 = v5 / 2
	ok    &= std::fabs( v2[6] - 1. ) <= 1e-10;

	v2[7] = v2[4] + v2[6];                           // v7 = v4 + v6
	ok    &= std::fabs( v2[7] - 1. ) <= 1e-10;

	return ok;
}

 # include <cstddef>                     // define size_t
# include <cmath>                       // for fabs function
extern bool exp_eps_for0(double *v0);   // computes zero order forward sweep
extern bool exp_eps_for1(double *v1);   // computes first order forward sweep
bool exp_eps_rev2(void)
{	bool ok = true;

	// set the value of v0[j], v1[j] for j = 1 , ... , 7
	double v0[8], v1[8];
	ok &= exp_eps_for0(v0);
	ok &= exp_eps_for1(v1);

	// initial all partial derivatives as zero
	double f_v0[8], f_v1[8];
	size_t j;
	for(j = 0; j < 8; j++)
	{	f_v0[j] = 0.;
		f_v1[j] = 0.;
	}

	// set partial derivative for f_7
	f_v1[7] = 1.;
	ok &= std::fabs( f_v1[7] - 1.  ) <= 1e-10; // partial f_7 w.r.t. v_7^1

	// f_6 = f_7( v_1^0 , ... , v_6^1 , v_4^0 + v_6^0, v_4^1 , v_6^1 )
	f_v0[4] += f_v0[7];
	f_v0[6] += f_v0[7];
	f_v1[4] += f_v1[7];
	f_v1[6] += f_v1[7];
	ok &= std::fabs( f_v0[4] - 0.  ) <= 1e-10; // partial f_6 w.r.t. v_4^0
	ok &= std::fabs( f_v0[6] - 0.  ) <= 1e-10; // partial f_6 w.r.t. v_6^0
	ok &= std::fabs( f_v1[4] - 1.  ) <= 1e-10; // partial f_6 w.r.t. v_4^1
	ok &= std::fabs( f_v1[6] - 1.  ) <= 1e-10; // partial f_6 w.r.t. v_6^1

	// f_5 = f_6( v_1^0 , ... , v_5^1 , v_5^0 / 2 , v_5^1 / 2 )
	f_v0[5] += f_v0[6] / 2.;
	f_v1[5] += f_v1[6] / 2.;
	ok &= std::fabs( f_v0[5] - 0.  ) <= 1e-10; // partial f_5 w.r.t. v_5^0
	ok &= std::fabs( f_v1[5] - 0.5 ) <= 1e-10; // partial f_5 w.r.t. v_5^1

	// f_4 = f_5( v_1^0 , ... , v_4^1 , v_3^0 * v_1^0 , 
	//            v_3^1 * v_1^0 + v_3^0 * v_1^1 )
	f_v0[1] += f_v0[5] * v0[3] + f_v1[5] * v1[3];
	f_v0[3] += f_v0[5] * v0[1] + f_v1[5] * v1[1];
	f_v1[1] += f_v1[5] * v0[3];
	f_v1[3] += f_v1[5] * v0[1];
	ok &= std::fabs( f_v0[1] - 0.5  ) <= 1e-10; // partial f_4 w.r.t. v_1^0
	ok &= std::fabs( f_v0[3] - 0.5  ) <= 1e-10; // partial f_4 w.r.t. v_3^0
	ok &= std::fabs( f_v1[1] - 0.25 ) <= 1e-10; // partial f_4 w.r.t. v_1^1
	ok &= std::fabs( f_v1[3] - 0.25 ) <= 1e-10; // partial f_4 w.r.t. v_3^1

	// f_3 = f_4(  v_1^0 , ... , v_3^1 , 1 + v_3^0 , v_3^1 )
	f_v0[3] += f_v0[4];
	f_v1[3] += f_v1[4];
	ok &= std::fabs( f_v0[3] - 0.5 ) <= 1e-10;  // partial f_3 w.r.t. v_3^0
	ok &= std::fabs( f_v1[3] - 1.25) <= 1e-10;  // partial f_3 w.r.t. v_3^1

	// f_2 = f_3( v_1^0 , ... , v_2^1 , v_2^0 , v_2^1 )
	f_v0[2] += f_v0[3];
	f_v1[2] += f_v1[3];
	ok &= std::fabs( f_v0[2] - 0.5 ) <= 1e-10;  // partial f_2 w.r.t. v_2^0
	ok &= std::fabs( f_v1[2] - 1.25) <= 1e-10;  // partial f_2 w.r.t. v_2^1

	// f_1 = f_2 ( v_1^0 , v_2^0 , v_1^0 , v_2^0 )
	f_v0[1] += f_v0[2];
	f_v1[1] += f_v1[2];
	ok &= std::fabs( f_v0[1] - 1.  ) <= 1e-10;  // partial f_1 w.r.t. v_1^0
	ok &= std::fabs( f_v1[1] - 1.5 ) <= 1e-10;  // partial f_1 w.r.t. v_1^1

	return ok;
}

 
# include <cppad/cppad.hpp>  // http://www.coin-or.org/CppAD/ 
# include "exp_eps.hpp"      // our example exponential function approximation
bool exp_eps_cppad(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::vector;    // can use any simple vector template class
	using CppAD::NearEqual; // checks if values are nearly equal

	// domain space vector
	size_t n = 2; // dimension of the domain space
	vector< AD<double> > U(n);
	U[0] = .5;    // value of x for this operation sequence
	U[1] = .2;    // value of e for this operation sequence

	// declare independent variables and start recording operation sequence
	CppAD::Independent(U);

	// evaluate our exponential approximation
	AD<double> x       = U[0];
	AD<double> epsilon = U[1];
	AD<double> apx = exp_eps(x, epsilon);  

	// range space vector
	size_t m = 1;  // dimension of the range space
	vector< AD<double> > Y(m);
	Y[0] = apx;    // variable that represents only range space component

	// Create f: U -> Y corresponding to this operation sequence
	// and stop recording. This also executes a zero order forward 
	// mode sweep using values in U for x and e.
	CppAD::ADFun<double> f(U, Y);

	// first order forward mode sweep that computes partial w.r.t x
	vector<double> du(n);      // differential in domain space
	vector<double> dy(m);      // differential in range space
	du[0] = 1.;                // x direction in domain space
	du[1] = 0.;
	dy    = f.Forward(1, du);  // partial w.r.t. x
	double check = 1.5;
	ok   &= NearEqual(dy[0], check, 1e-10, 1e-10);

	// first order reverse mode sweep that computes the derivative
	vector<double>  w(m);     // weights for components of the range
	vector<double> dw(n);     // derivative of the weighted function
	w[0] = 1.;                // there is only one weight 
	dw   = f.Reverse(1, w);   // derivative of w[0] * exp_eps(x, epsilon)
	check = 1.5;              // partial w.r.t. x
	ok   &= NearEqual(dw[0], check, 1e-10, 1e-10);
	check = 0.;               // partial w.r.t. epsilon
	ok   &= NearEqual(dw[1], check, 1e-10, 1e-10);

	// second order forward sweep that computes
	// second partial of exp_eps(x, epsilon) w.r.t. x
	vector<double> x2(n);     // second order Taylor coefficients
	vector<double> y2(m);
	x2[0] = 0.;               // evaluate partial w.r.t x
	x2[1] = 0.;
	y2    = f.Forward(2, x2);
	check = 0.5 * 1.;         // Taylor coef is 1/2 second derivative
	ok   &= NearEqual(y2[0], check, 1e-10, 1e-10);

	// second order reverse sweep that computes
	// derivative of partial of exp_eps(x, epsilon) w.r.t. x
	dw.resize(2 * n);         // space for first and second derivative
	dw    = f.Reverse(2, w);
	check = 1.;               // result should be second derivative
	ok   &= NearEqual(dw[0*2+1], check, 1e-10, 1e-10);

	return ok;
}

 

// system include files used for I/O
# include <iostream>

// external complied tests
extern bool exp_2(void);
extern bool exp_2_cppad(void);
extern bool exp_2_for1(void);
extern bool exp_2_for2(void);
extern bool exp_2_rev1(void);
extern bool exp_2_rev2(void);
extern bool exp_2_for0(void);
extern bool exp_eps(void);
extern bool exp_eps_cppad(void);
extern bool exp_eps_for1(void);
extern bool exp_eps_for2(void);
extern bool exp_eps_for0(void);
extern bool exp_eps_rev1(void);
extern bool exp_eps_rev2(void);

namespace {
	static size_t Run_ok_count    = 0;
	static size_t Run_error_count = 0;
	bool Run(bool TestOk(void), std::string name)
	{	bool ok               = true;
		std::streamsize width =  20;         
		std::cout.width( width );
		std::cout.setf( std::ios_base::left );
		std::cout << name.c_str();
		//
		ok &= name.size() < size_t(width);
		ok &= TestOk();
		if( ok )
		{	std::cout << "OK" << std::endl;
			Run_ok_count++;
		}
		else
		{	std::cout << "Error" << std::endl;
			Run_error_count++;
		}
		return ok;
	}
}

// main program that runs all the tests
int main(void)
{	bool ok = true;
	using namespace std;

	// This comment is used by OneTest 

	// external compiled tests
	ok &= Run( exp_2,           "exp_2"          );
	ok &= Run( exp_2_cppad,     "exp_2_cppad"    );
	ok &= Run( exp_2_for0,      "exp_2_for0"     );
	ok &= Run( exp_2_for1,      "exp_2_for1"     );
	ok &= Run( exp_2_for2,      "exp_2_for2"     );
	ok &= Run( exp_2_rev1,      "exp_2_rev1"     );
	ok &= Run( exp_2_rev2,      "exp_2_rev2"     );
	ok &= Run( exp_eps,         "exp_eps"        );
	ok &= Run( exp_eps_cppad,   "exp_eps_cppad"  );
	ok &= Run( exp_eps_for0,    "exp_eps_for0"   );
	ok &= Run( exp_eps_for1,    "exp_eps_for1"   );
	ok &= Run( exp_eps_for2,    "exp_eps_for2"   );
	ok &= Run( exp_eps_rev1,    "exp_eps_rev1"   );
	ok &= Run( exp_eps_rev2,    "exp_eps_rev2"   );
	if( ok )
		cout << "All " << int(Run_ok_count) << " tests passed." << endl;
	else	cout << int(Run_error_count) << " tests failed." << endl;

	return static_cast<int>( ! ok );
}

 

# include <cppad/cppad.hpp>

bool Default(void)
{	bool ok = true;
	using CppAD::AD;

	// default AD constructor
	AD<double> x, y;

	// check that they are parameters
	ok &= Parameter(x);
	ok &= Parameter(y);

	// assign them values
	x = 3.; 
	y = 4.;

	// just check a simple operation
	ok &= (x + y == 7.);

	return ok;
}

 

# include <cppad/cppad.hpp>

bool CopyAD(void)
{	bool ok = true;   // initialize test result flag
	using CppAD::AD;  // so can use AD in place of CppAD::AD

	// domain space vector
	size_t n = 1;
	CPPAD_TEST_VECTOR< AD<double> > x(n);
	x[0]     = 2.;

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// create an AD<double> that does not depend on x
	AD<double> b = 3.;   

	// use copy constructor 
	AD<double> u(x[0]);    
	AD<double> v = b;

	// check which are parameters
	ok &= Variable(u);
	ok &= Parameter(v);

	// range space vector
	size_t m = 2;
	CPPAD_TEST_VECTOR< AD<double> > y(m);
	y[0]  = u;
	y[1]  = v;

	// create f: x -> y and vectors used for derivative calculations
	CppAD::ADFun<double> f(x, y);
	CPPAD_TEST_VECTOR<double> dx(n);
	CPPAD_TEST_VECTOR<double> dy(m);
 
 	// check parameters flags
 	ok &= ! f.Parameter(0);
 	ok &=   f.Parameter(1);

	// check function values
	ok &= ( y[0] == 2. );
	ok &= ( y[1] == 3. );

	// forward computation of partials w.r.t. x[0]
	dx[0] = 1.;
	dy    = f.Forward(1, dx);
	ok   &= ( dy[0] == 1. );   // du / dx
	ok   &= ( dy[1] == 0. );   // dv / dx

	return ok;
}

 

# include <cppad/cppad.hpp>

bool CopyBase(void)
{	bool ok = true;    // initialize test result flag
	using CppAD::AD;   // so can use AD in place of CppAD::AD

	// construct directly from Base where Base is double 
	AD<double> x(1.); 

	// construct from a type that converts to Base where Base is double
	AD<double> y = 2;

	// construct from a type that converts to Base where Base = AD<double>
	AD< AD<double> > z(3); 

	// check that resulting objects are parameters
	ok &= Parameter(x);
	ok &= Parameter(y);
	ok &= Parameter(z);

	// check values of objects (compare AD<double> with double)
	ok &= ( x == 1.);
	ok &= ( y == 2.);
	ok &= ( Value(z) == 3.);

	// user constructor through the static_cast template function
	x   = static_cast < AD<double> >( 4 );
	z  = static_cast < AD< AD<double> > >( 5 );

	ok &= ( x == 4. );
	ok &= ( Value(z) == 5. );

	return ok;
}

 

# include <cppad/cppad.hpp>

bool Eq(void)
{	bool ok = true;
	using CppAD::AD;

	// domain space vector
	size_t n = 3;
	CPPAD_TEST_VECTOR< AD<double> > x(n);
	x[0]     = 2;      // AD<double> = int
	x[1]     = 3.;     // AD<double> = double
	x[2]     = x[1];   // AD<double> = AD<double>

	// declare independent variables and start tape recording
	CppAD::Independent(x);
	
	// range space vector 
	size_t m = 3;
	CPPAD_TEST_VECTOR< AD<double> > y(m);

	// assign an AD<Base> object equal to an independent variable
	// (choose the first independent variable to check a special case)
	// use the value returned by the assignment (for another assignment)
	y[0] = y[1] = x[0];  

	// assign an AD<Base> object equal to an expression 
	y[1] = x[1] + 1.;
	y[2] = x[2] + 2.;

	// check that all the resulting components of y depend on x
	ok &= Variable(y[0]);  // y[0] = x[0]
	ok &= Variable(y[1]);  // y[1] = x[1] + 1
	ok &= Variable(y[2]);  // y[2] = x[2] + 2

	// construct f : x -> y and stop the tape recording
	CppAD::ADFun<double> f(x, y);

	// check variable values
	ok &= ( y[0] == 2.);
	ok &= ( y[1] == 4.);
	ok &= ( y[2] == 5.);

	// compute partials w.r.t x[1]
	CPPAD_TEST_VECTOR<double> dx(n);
	CPPAD_TEST_VECTOR<double> dy(m);
	dx[0] = 0.;
	dx[1] = 1.;
	dx[2] = 0.;
	dy   = f.Forward(1, dx);
	ok  &= (dy[0] == 0.);  // dy[0] / dx[1]
	ok  &= (dy[1] == 1.);  // dy[1] / dx[1]
	ok  &= (dy[2] == 0.);  // dy[2] / dx[1]

	// compute the derivative y[2]
	CPPAD_TEST_VECTOR<double>  w(m);
	CPPAD_TEST_VECTOR<double> dw(n);
	w[0] = 0.;
	w[1] = 0.;
	w[2] = 1.;
	dw   = f.Reverse(1, w);
	ok  &= (dw[0] == 0.);  // dy[2] / dx[0]
	ok  &= (dw[1] == 0.);  // dy[2] / dx[1]
	ok  &= (dw[2] == 1.);  // dy[2] / dx[2]

	// assign a VecAD<Base>::reference
	CppAD::VecAD<double> v(1);
	AD<double> zero(0);
	v[zero] = 5.;
	ok     &= (v[0] == 5.);

	return ok;
}

 

# include <cppad/cppad.hpp>

bool Value(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::Value;

	// domain space vector
	size_t n = 2;
	CPPAD_TEST_VECTOR< AD<double> > x(n);
	x[0] = 3.;
	x[1] = 4.;

	// check value before recording
	ok &= (Value(x[0]) == 3.);
	ok &= (Value(x[1]) == 4.);

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// range space vector 
	size_t m = 1;
	CPPAD_TEST_VECTOR< AD<double> > y(m);
	y[0] = - x[1];

	// cannot call Value(x[j]) or Value(y[0]) here (currently variables)
	AD<double> p = 5.;        // p is a parameter (does not depend on x)
	ok &= (Value(p) == 5.);

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f(x, y);

	// can call Value(x[j]) or Value(y[0]) here (currently parameters)
	ok &= (Value(x[0]) ==  3.);
	ok &= (Value(x[1]) ==  4.);
	ok &= (Value(y[0]) == -4.);

	return ok;
}

 

# include <cppad/cppad.hpp>

bool Integer(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::Integer;

	// domain space vector
	size_t n = 2;
	CPPAD_TEST_VECTOR< AD<double> > x(n);
	x[0] = 3.5;
	x[1] = 4.5;

	// check integer before recording
	ok &= (Integer(x[0]) == 3);
	ok &= (Integer(x[1]) == 4);

	// start recording

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// check integer during recording
	ok &= (Integer(x[0]) == 3);
	ok &= (Integer(x[1]) == 4);

	// check integer for VecAD element
	CppAD::VecAD<double> v(1);
	AD<double> zero(0);
	v[zero] = 2;
	ok &= (Integer(v[zero]) == 2);

	// range space vector 
	size_t m = 1;
	CPPAD_TEST_VECTOR< AD<double> > y(m);
	y[0] = - x[1];

	// create f: x -> y and stop recording
	CppAD::ADFun<double> f(x, y);

	// check integer after recording
	ok &= (Integer(x[0]) ==  3.);
	ok &= (Integer(x[1]) ==  4.);
	ok &= (Integer(y[0]) == -4.);

	return ok;
}

 

# include <cppad/cppad.hpp>

# include <sstream>  // std::ostringstream
# include <string>   // std::string
# include <iomanip>  // std::setprecision, setw, setfill, right

namespace {
	template <class S>
	void set_ostream(S &os)
	{	os 
		<< std::setprecision(4) // 4 digits of precision
		<< std::setw(6)         // 6 characters per field
		<< std::setfill(' ')    // fill with spaces
		<< std::right;          // adjust value to the right
	}
}

bool Output(void)
{	bool ok = true;

	// This output stream is an ostringstream for testing purposes.
	// You can use << with other types of streams; i.e., std::cout.
	std::ostringstream stream;

	// ouput an AD<double> object
	CppAD::AD<double>  pi = 4. * atan(1.); // 3.1415926536
	set_ostream(stream);
	stream << pi;

	// ouput a VecAD<double>::reference object
	CppAD::VecAD<double> v(1);
	CppAD::AD<double> zero(0);
	v[zero]   = exp(1.);                  // 2.7182818285
	set_ostream(stream); 
	stream << v[zero];

	// convert output from stream to string
	std::string str = stream.str();

	// check the output
	ok      &= (str == " 3.142 2.718");

	return ok;
}

 
# include <cppad/cppad.hpp>

int main(void)
{
	using namespace CppAD;
	using std::cout;
	using std::endl;

	// independent variable vector
	size_t n = 1;
	CPPAD_TEST_VECTOR< AD<double> > X(n);
	X[0] = 1.;
	Independent(X);

	// print a VecAD<double>::reference object that is a parameter
	VecAD<double> V(1);
	AD<double> Zero(0);
	V[Zero] = 0.;
	PrintFor("v[0] = ", V[Zero]); 

	// dependent variable vector 
	size_t m = 1;
	CPPAD_TEST_VECTOR< AD<double> > Y(m);
	Y[0] = V[Zero] + X[0];

	// First print a newline to separate this from previous output,
	// then print an AD<double> object that is a variable.
	PrintFor(  "\nv[0] + x[0] = ", Y[0]); 

	// define f: x -> y and stop tape recording
	ADFun<double> f(X, Y); 

	// zero order forward with x[0] = 2 
	CPPAD_TEST_VECTOR<double> x(n);
	x[0] = 2.;


	cout << "v[0] = 0" << endl; 
	cout << "v[0] + x[0] = 2" << endl; 
	// Developer Note: ./makefile.am "Test passes" to begin next output line
	cout << "Test passes if two lines above repeat below:" << endl;
	f.Forward(0, x);	

	// print a new line after output
	std::cout << std::endl;

	return 0;
}

 
	v[0] = 0
	v[0] + x[0] = 2
	Test passes if two lines above repeat below:
	v[0] = 0
	v[0] + x[0] = 2

 

# include <cppad/cppad.hpp>


bool Var2Par(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::Value;
	using CppAD::Var2Par;

	// domain space vector
	size_t n = 2;
	CPPAD_TEST_VECTOR< AD<double> > x(n);
	x[0] = 3.;
	x[1] = 4.;

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// range space vector 
	size_t m = 1;
	CPPAD_TEST_VECTOR< AD<double> > y(m);
	y[0] = - x[1] * Var2Par(x[0]);    // same as y[0] = -x[1] * 3.;

	// cannot call Value(x[j]) or Value(y[0]) here (currently variables)
	ok &= ( Value( Var2Par(x[0]) ) == 3. );
	ok &= ( Value( Var2Par(x[1]) ) == 4. ); 
	ok &= ( Value( Var2Par(y[0]) ) == -12. ); 

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f(x, y);

	// can call Value(x[j]) or Value(y[0]) here (currently parameters)
	ok &= (Value(x[0]) ==  3.);
	ok &= (Value(x[1]) ==  4.);
	ok &= (Value(y[0]) == -12.);

	// evaluate derivative of y w.r.t x
	CPPAD_TEST_VECTOR<double> w(m);
	CPPAD_TEST_VECTOR<double> dw(n);
	w[0] = 1.;
	dw   = f.Reverse(1, w);
	ok  &= (dw[0] == 0.);  // derivative of y[0] w.r.t x[0] is zero
	ok  &= (dw[1] == -3.); // derivative of y[0] w.r.t x[1] is 3

	return ok;
}

 

# include <cppad/cppad.hpp>

bool UnaryPlus(void)
{	bool ok = true;
	using CppAD::AD;


	// domain space vector
	size_t n = 1;
	CPPAD_TEST_VECTOR< AD<double> > x(n);
	x[0]      = 3.;

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// range space vector 
	size_t m = 1;
	CPPAD_TEST_VECTOR< AD<double> > y(m);
	y[0] = + x[0];

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f(x, y);

	// check values
	ok &= ( y[0] == 3. );

	// forward computation of partials w.r.t. x[0]
	CPPAD_TEST_VECTOR<double> dx(n);
	CPPAD_TEST_VECTOR<double> dy(m);
	size_t p = 1;
	dx[0]    = 1.;
	dy       = f.Forward(p, dx);
	ok      &= ( dy[0] == 1. );   // dy[0] / dx[0]

	// reverse computation of dertivative of y[0]
	CPPAD_TEST_VECTOR<double>  w(m);
	CPPAD_TEST_VECTOR<double> dw(n);
	w[0] = 1.;
	dw   = f.Reverse(p, w);
	ok &= ( dw[0] == 1. );       // dy[0] / dx[0]

	// use a VecAD<Base>::reference object with unary plus
	CppAD::VecAD<double> v(1);
	AD<double> zero(0);
	v[zero] = x[0];
	AD<double> result = + v[zero];
	ok     &= (result == y[0]);
	 
	return ok;
}