自学教程：C++ CPPAD_ASSERT_KNOWN函数代码示例

void ode_evaluate(	CppAD::vector<Float> &x  , 	size_t m                 , 	CppAD::vector<Float> &fm ){	typedef CppAD::vector<Float> Vector;	size_t n = x.size();	size_t ell;	CPPAD_ASSERT_KNOWN( m == 0 || m == 1,		"ode_evaluate: m is not zero or one"	);	CPPAD_ASSERT_KNOWN( 		((m==0) & (fm.size()==n)) || ((m==1) & (fm.size()==n*n)),		"ode_evaluate: the size of fm is not correct"	);	if( m == 0 )		ell = n;	else	ell = n + n * n;	// set up the case we are integrating	Float  ti   = 0.;	Float  tf   = 1.;	Float  smin = 1e-5;	Float smax  = 1.;	Float scur  = 1.;	Float erel  = 0.;	vector<Float> yi(ell), eabs(ell);	size_t i, j;	for(i = 0; i < ell; i++)	{	eabs[i] = 1e-10;		if( i < n )			yi[i] = 1.;		else	yi[i]  = 0.;	}	// return values	Vector yf(ell), ef(ell), maxabs(ell);	size_t nstep;	// construct ode method for taking one step	ode_evaluate_method<Float> method(m, x);	// solve differential equation	yf = OdeErrControl(method, 		ti, tf, yi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);	if( m == 0 )	{	for(i = 0; i < n; i++)			fm[i] = yf[i];	}	else	{	for(i = 0; i < n; i++)			for(j = 0; j < n; j++)				fm[i * n + j] = yf[n + i * n + j];	}	return;}

void put_check_for_nan(const CppAD::vector<Base>& vec, std::string& file_name){	size_t char_size       = sizeof(Base) * vec.size();	const char* char_ptr   = reinterpret_cast<const char*>( vec.data() );# if CPPAD_HAS_MKSTEMP	char pattern[] = "/tmp/fileXXXXXX";	int fd = mkstemp(pattern);	file_name = pattern;	write(fd, char_ptr, char_size);	close(fd);# else# if CPPAD_HAS_TMPNAM_S		std::vector<char> name(L_tmpnam_s);		if( tmpnam_s( name.data(), L_tmpnam_s ) != 0 )		{	CPPAD_ASSERT_KNOWN(				false,				"Cannot create a temporary file name"			);		}		file_name = name.data();# else		file_name = tmpnam( CPPAD_NULL );# endif	std::fstream file_out(file_name.c_str(), std::ios::out|std::ios::binary );	file_out.write(char_ptr, char_size);	file_out.close();# endif	return;}

	void sparse_jac_fun(		size_t                       m    ,		size_t                       n    ,		const FloatVector&           x    ,		const CppAD::vector<size_t>& row  , 		const CppAD::vector<size_t>& col  , 		size_t                       p    ,		FloatVector&                 fp   )	{		// check numeric type specifications		CheckNumericType<Float>();		// check value of p		CPPAD_ASSERT_KNOWN(			p == 0 || p == 1,			"sparse_jac_fun: p != 0 and p != 1"		);		size_t K = row.size();		CPPAD_ASSERT_KNOWN(			K >= m,			"sparse_jac_fun: row.size() < m"		);		size_t i, j, k;		if( p == 0 )			for(i = 0; i < m; i++)				fp[i] = Float(0);		Float t;		for(k = 0; k < K; k++)		{	i    = row[k];			j    = col[k];			t    = exp( x[j] * x[j] / 2.0 );				switch(p)			{				case 0:				fp[i] += t;				break;				case 1:				fp[k] = t * x[j];				break;			}		}	}

	void sparse_hes_fun(		size_t                       n    ,		const FloatVector&           x    ,		const CppAD::vector<size_t>& row  , 		const CppAD::vector<size_t>& col  , 		size_t                       p    ,		FloatVector&                fp    )	{		// check numeric type specifications		CheckNumericType<Float>();		// check value of p		CPPAD_ASSERT_KNOWN(			p < 3,			"sparse_hes_fun: p > 2"		);		size_t i, j, k;		size_t size = 1;		for(k = 0; k < p; k++)			size *= n;		for(k = 0; k < size; k++)			fp[k] = Float(0);		size_t K = row.size();		Float t;		Float dt_i;		Float dt_j;		for(k = 0; k < K; k++)		{	i    = row[k];			j    = col[k];			t    = exp( x[i] * x[j] );				dt_i = t * x[j];			dt_j = t * x[i];			switch(p)			{				case 0:				fp[0] += t;				break;				case 1:				fp[i] += dt_i;				fp[j] += dt_j;				break;				case 2:				fp[i * n + i] += dt_i * x[j];				fp[i * n + j] += t + dt_j * x[j];				//				fp[j * n + i] += t + dt_i * x[i];				fp[j * n + j] += dt_j * x[i];				break;			}		}				}

Vector ADFun<Base>::ForOne(const Vector &x, size_t j){	size_t j1;	size_t n = Domain();	size_t m = Range();	// check Vector is Simple Vector class with Base type elements	CheckSimpleVector<Base, Vector>();	CPPAD_ASSERT_KNOWN(		x.size() == n,		"ForOne: Length of x not equal domain dimension for f"	);	CPPAD_ASSERT_KNOWN(		j < n,		"ForOne: the index j is not less than domain dimension for f"	);	// point at which we are evaluating the second partials	Forward(0, x);	// direction in which are are taking the derivative	Vector dx(n);	for(j1 = 0; j1 < n; j1++)		dx[j1] = Base(0);	dx[j] = Base(1);	// dimension the return value	Vector dy(m);	// compute the return value	dy = Forward(1, dx);	return dy;}

void parallel_ad(void){   CPPAD_ASSERT_KNOWN(        ! thread_alloc::in_parallel() ,        "parallel_ad must be called before entering parallel execution mode."    );    CPPAD_ASSERT_KNOWN(        AD<Base>::tape_ptr() == CPPAD_NULL ,        "parallel_ad cannot be called while a tape recording is in progress"    );    // ensure statics in following functions are initialized    elapsed_seconds();    ErrorHandler::Current();    local::NumArg(local::BeginOp);    local::NumRes(local::BeginOp);    local::one_element_std_set<size_t>();    local::two_element_std_set<size_t>();    // the sparse_pack class has member functions with static data    local::sparse_pack sp;    sp.resize(1, 1);       // so can call add_element    sp.add_element(0, 0);  // has static data    sp.clear(0);           // has static data    sp.is_element(0, 0);   // has static data    local::sparse_pack::const_iterator itr(sp, 0); // has static data    ++itr;                                  // has static data    // statics that depend on the value of Base    AD<Base>::tape_id_ptr(0);    AD<Base>::tape_handle(0);    discrete<Base>::List();    CheckSimpleVector< Base, CppAD::vector<Base> >();    CheckSimpleVector< AD<Base>, CppAD::vector< AD<Base> > >();}

void  AD<Base>::tape_delete(size_t id_old){	size_t thread = id_old % CPPAD_MAX_NUM_THREADS;	CPPAD_ASSERT_KNOWN(		thread == thread_alloc::thread_num(),		"AD tape recording must stop in same thread as it started in."	);	size_t        *id   = id_handle(thread);	ADTape<Base> **tape = tape_handle(thread);	CPPAD_ASSERT_UNKNOWN( *id   == id_old     );	CPPAD_ASSERT_UNKNOWN( *tape != CPPAD_NULL );	// increase the id for this thread in a way such that 	// thread = id % CPPAD_MAX_NUM_THREADS	CPPAD_ASSERT_KNOWN(	size_t(*id) + CPPAD_MAX_NUM_THREADS <= 		std::numeric_limits<CPPAD_TAPE_ID_TYPE>::max(),	"To many different tapes given the type used for CPPAD_TAPE_ID_TYPE"	);	*id  += CPPAD_MAX_NUM_THREADS;	// delete the old tape for this thread	delete ( *tape );	*tape = CPPAD_NULL;	return;}

Vector ADFun<Base>::RevOne(const Vector  &x, size_t i){	size_t i1;	size_t n = Domain();	size_t m = Range();	// check Vector is Simple Vector class with Base type elements	CheckSimpleVector<Base, Vector>();	CPPAD_ASSERT_KNOWN(		x.size() == n,		"RevOne: Length of x not equal domain dimension for f"	);	CPPAD_ASSERT_KNOWN(		i < m,		"RevOne: the index i is not less than range dimension for f"	);	// point at which we are evaluating the derivative	Forward(0, x);	// component which are are taking the derivative of	Vector w(m);	for(i1 = 0; i1 < m; i1++)		w[i1] = 0.;	w[i] = Base(1);	// dimension the return value	Vector dw(n);	// compute the return value	dw = Reverse(1, w);	return dw;}

void color_general_colpack(	const VectorSet&        pattern ,	const VectorSize&       row     ,	const VectorSize&       col     ,	CppAD::vector<size_t>&  color   ){	size_t i, j, k;	size_t m = pattern.n_set();	size_t n = pattern.end();	// Determine number of non-zero entries in each row	CppAD::vector<size_t> n_nonzero(m);	size_t n_nonzero_total = 0;	for(i = 0; i < m; i++)	{	n_nonzero[i] = 0;		typename VectorSet::const_iterator pattern_itr(pattern, i);		j = *pattern_itr;		while( j != pattern.end() )		{	n_nonzero[i]++;			j = *(++pattern_itr);		}		n_nonzero_total += n_nonzero[i];	}	// Allocate memory and fill in Adolc sparsity pattern	CppAD::vector<unsigned int*> adolc_pattern(m);	CppAD::vector<unsigned int>  adolc_memory(m + n_nonzero_total);	size_t i_memory = 0;	for(i = 0; i < m; i++)	{	adolc_pattern[i]    = adolc_memory.data() + i_memory;		CPPAD_ASSERT_KNOWN(			std::numeric_limits<unsigned int>::max() >= n_nonzero[i],			"Matrix is too large for colpack"		);		adolc_pattern[i][0] = static_cast<unsigned int>( n_nonzero[i] );		typename VectorSet::const_iterator pattern_itr(pattern, i);		j = *pattern_itr;		k = 1;		while(j != pattern.end() )		{			CPPAD_ASSERT_KNOWN(				std::numeric_limits<unsigned int>::max() >= j,				"Matrix is too large for colpack"			);			adolc_pattern[i][k++] = static_cast<unsigned int>( j );			j = *(++pattern_itr);		}		CPPAD_ASSERT_UNKNOWN( k == 1 + n_nonzero[i] );		i_memory += k;	}	CPPAD_ASSERT_UNKNOWN( i_memory == m + n_nonzero_total );	// Must use an external routine for this part of the calculation because	// ColPack/ColPackHeaders.h has as 'using namespace std' at global level.	cppad_colpack_general(color, m, n, adolc_pattern);	return;}

Vector ADFun<Base>::Hessian(const Vector &x, const Vector &w){	size_t j;	size_t k;	size_t n = Domain();	// check Vector is Simple Vector class with Base type elements	CheckSimpleVector<Base, Vector>();	CPPAD_ASSERT_KNOWN(		size_t(x.size()) == n,		"Hessian: length of x not equal domain dimension for f"	);	CPPAD_ASSERT_KNOWN(		size_t(w.size()) == Range(),		"Hessian: length of w not equal range dimension for f"	);	// point at which we are evaluating the Hessian	Forward(0, x);	// define the return value	Vector hes(n * n);	// direction vector for calls to forward	Vector u(n);	for(j = 0; j < n; j++)		u[j] = Base(0);	// location for return values from Reverse	Vector ddw(n * 2);	// loop over forward directions	for(j = 0; j < n; j++)	{	// evaluate partials of entire function w.r.t. j-th coordinate		u[j] = Base(1);		Forward(1, u);		u[j] = Base(0);		// evaluate derivative of partial corresponding to F_i		ddw = Reverse(2, w);		// return desired components		for(k = 0; k < n; k++)			hes[k * n + j] = ddw[k * 2 + 1];	}	return hes;}

void ADFun<Base>::RevSparseHesCase(	bool              set_type         ,	bool              transpose        ,  	size_t            q                ,  	const VectorSet&  s                ,	VectorSet&        h                ){	size_t n = Domain(); 		h.resize(q * n );	CPPAD_ASSERT_KNOWN( 		for_jac_sparse_pack_.n_set() > 0,		"RevSparseHes: previous stored call to ForSparseJac did not "		"use bool for the elements of r."	);	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 );	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == total_num_var_ );		// use sparse_pack for the calculation	CppAD::RevSparseHesBool( 		transpose                ,		q                        ,		s                        ,		h                        ,		total_num_var_           ,		dep_taddr_               ,		ind_taddr_               ,		play_                    ,		for_jac_sparse_pack_ 	);}

void fun_record(	cppad_ipopt_fg_info*                              fg_info ,	size_t                                            k       ,	const SizeVector&                                 p       ,	const SizeVector&                                 q       ,	size_t                                            n       ,	const NumVector&                                  x       ,	const SizeVector&                                 J       ,	CppAD::vector< CppAD::ADFun<Ipopt::Number> >&     r_fun   ){	size_t j;	// extract u from x	ADVector u(q[k]);	for(j = 0; j < q[k]; j++)	{	// when NDEBUG is not defined, this error should be caught		// during the cppad_ipopt_nlp constructor.		CPPAD_ASSERT_UNKNOWN( J[j] < n );		u[j] = x[ J[j] ];	}	// start the recording	CppAD::Independent(u);	// record the evaulation of r_k (u)	ADVector r_k = fg_info->eval_r(k, u);	CPPAD_ASSERT_KNOWN( r_k.size() == p[k] ,	"cppad_ipopt_nlp: eval_r return value size not equal to p[k]."	);	// stop the recording and store operation sequence in	r_fun[k].Dependent(u, r_k);}

 /// not a number static Float quiet_NaN(void) {   CPPAD_ASSERT_KNOWN(         false,         "numeric_limits<Float>::quiet_NaN() is not specialized for this Float"     );     return Float(0); }

	/*! 	Link from forward mode sweep to users routine.	/param index	index for this function in the list of all user_atomic objects	/param id	extra information vector that is just passed through by CppAD,	and possibly used by user's routines.	/param k	order for this forward mode calculation.	/param n	domain space size for this calcualtion.	/param m	range space size for this calculation.	/param tx	Taylor coefficients corresponding to /c x for this calculation.	/param ty	Taylor coefficient corresponding to /c y for this calculation	See the forward mode in user's documentation for user_atomic  	*/	static void forward(		size_t                index , 		size_t                   id ,		size_t                    k ,		size_t                    n , 		size_t                    m , 		const vector<Base>&      tx ,		vector<Base>&            ty )	{	# ifdef _OPENMP		vector<bool> empty(0);# else		static vector<bool> empty(0);# endif				CPPAD_ASSERT_UNKNOWN( tx.size() >= n * k );		CPPAD_ASSERT_UNKNOWN( ty.size() >= m * k );		CPPAD_ASSERT_UNKNOWN(index < List().size() );		user_atomic* op = List()[index];		bool ok = op->f_(id, k, n, m, empty, empty, tx, ty);		if( ! ok )		{	std::stringstream ss;			ss << k;			std::string msg = "user_atomic: ";			msg = msg + op->name_ + ": ok returned false from " + ss.str()			    + " order forward mode calculation";			CPPAD_ASSERT_KNOWN(false, msg.c_str());		}	}

inline void forward_sqrt_op_dir(	size_t q           ,	size_t r           ,	size_t i_z         ,	size_t i_x         ,	size_t cap_order   ,	Base*  taylor      ){	// check assumptions	CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );	CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );	CPPAD_ASSERT_UNKNOWN( 0 < q );	CPPAD_ASSERT_UNKNOWN( q < cap_order );	// Taylor coefficients corresponding to argument and result	size_t num_taylor_per_var = (cap_order-1) * r + 1;	Base* z = taylor + i_z * num_taylor_per_var;	Base* x = taylor + i_x * num_taylor_per_var;	CPPAD_ASSERT_KNOWN(		x[0] != Base(0),		"Forward: attempt to take derivatve of square root of zero"	)	size_t m = (q-1) * r + 1;	for(size_t ell = 0; ell < r; ell++)	{	z[m+ell] = Base(0);		for(size_t k = 1; k < q; k++)			z[m+ell] -= Base(k) * z[(k-1)*r+1+ell] * z[(q-k-1)*r+1+ell];		z[m+ell] /= Base(q);		z[m+ell] += x[m+ell] / Base(2);		z[m+ell] /= z[0];	}}

	void operator()(const ADVector& ax, ADVector& ay, size_t id = 0)	{	CPPAD_ASSERT_KNOWN(			id == 0,			"checkpoint: id is non-zero in afun(ax, ay, id)"		);		this->atomic_base<Base>::operator()(ax, ay, id);	}

	/*! 	Link from reverse mode sweep to users routine.	/param index	index in the list of all user_atomic objects	corresponding to this function.	/param id	extra information vector that is just passed through by CppAD,	and possibly used by user's routines.	/param k	order for this forward mode calculation.	/param n	domain space size for this calcualtion.	/param m	range space size for this calculation.	/param tx	Taylor coefficients corresponding to /c x for this calculation.	/param ty	Taylor coefficient corresponding to /c y for this calculation	/param px	Partials w.r.t. the /c x Taylor coefficients.	/param py	Partials w.r.t. the /c y Taylor coefficients.	See reverse mode documentation for user_atomic  	*/	static void reverse(		size_t               index , 		size_t                  id ,		size_t                   k ,		size_t                   n , 		size_t                   m , 		const vector<Base>&     tx ,		const vector<Base>&     ty ,		vector<Base>&           px ,		const vector<Base>&     py )	{		CPPAD_ASSERT_UNKNOWN(index < List().size() );		CPPAD_ASSERT_UNKNOWN( tx.size() >= n * k );		CPPAD_ASSERT_UNKNOWN( px.size() >= n * k );		CPPAD_ASSERT_UNKNOWN( ty.size() >= m * k );		CPPAD_ASSERT_UNKNOWN( py.size() >= m * k );		user_atomic* op = List()[index];		bool ok = op->r_(id, k, n, m, tx, ty, px, py);		if( ! ok )		{	std::stringstream ss;			ss << k;			std::string msg = "user_atomic: ";			msg = op->name_ + ": ok returned false from " + ss.str() 			    + " order reverse mode calculation";			CPPAD_ASSERT_KNOWN(false, msg.c_str());		}	}

 /// machine epsilon static Float epsilon(void) {   CPPAD_ASSERT_KNOWN(         false,         "numeric_limits<Float>::epsilon() is not specialized for this Float"     );     return Float(0); }

	/*! 	Link from reverse Hessian sparsity sweep to users routine.	/param index	index in the list of all user_atomic objects	corresponding to this function.	/param id	extra information vector that is just passed through by CppAD,	and possibly used by user's routines.	/param n	domain space size for this calcualtion.	/param m	range space size for this calculation.	/param q	is the column dimension for the sparsity partterns.	/param r	is the forward Jacobian sparsity pattern w.r.t the argument vector x	/param s	is the reverse Jacobian sparsity pattern w.r.t the result vector y.	/param t	is the reverse Jacobian sparsity pattern w.r.t the argument vector x.	/param u	is the Hessian sparsity pattern w.r.t the result vector y.	/param v	is the Hessian sparsity pattern w.r.t the argument vector x.	*/	static void rev_hes_sparse(		size_t                            index ,		size_t                               id ,		size_t                                n , 		size_t                                m , 		size_t                                q ,		vector< std::set<size_t> >&           r ,		const vector<bool>&                   s ,		vector<bool>&                         t ,		const vector< std::set<size_t> >&     u ,		vector< std::set<size_t> >&           v )	{		CPPAD_ASSERT_UNKNOWN(index < List().size() );		CPPAD_ASSERT_UNKNOWN( r.size() >= n );		CPPAD_ASSERT_UNKNOWN( s.size() >= m );		CPPAD_ASSERT_UNKNOWN( t.size() >= n );		CPPAD_ASSERT_UNKNOWN( u.size() >= m );		CPPAD_ASSERT_UNKNOWN( v.size() >= n );		user_atomic* op = List()[index];		bool ok = op->rhs_(id, n, m, q, r, s, t, u, v);		if( ! ok )		{	std::string msg = "user_atomic: ";			msg = msg + op->name_ 			    + ": ok returned false from rev_jac_sparse calculation";			CPPAD_ASSERT_KNOWN(false, msg.c_str());		}	}

void ADFun<Base>::RevSparseHesCase(	const std::set<size_t>&   set_type         ,	bool                      transpose        ,  	size_t                    q                ,  	const VectorSet&          s                ,	VectorSet&                h                ){	size_t n = Domain();	if( transpose )		h.resize(n);	else	h.resize(q);	CPPAD_ASSERT_KNOWN( 		for_jac_sparse_set_.n_set() > 0,		"RevSparseHes: previous stored call to ForSparseJac did not "		"use std::set<size_t> for the elements of r."	);	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 );	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == num_var_tape_  );		// use sparse_pack for the calculation	CppAD::RevSparseHesSet( 		transpose                ,		q                        ,		s                        ,		h                        ,		num_var_tape_            ,		dep_taddr_               ,		ind_taddr_               ,		play_                    ,		for_jac_sparse_set_ 	);}

inline void forward_sqrt_op(	size_t j           ,	size_t i_z         ,	size_t i_x         ,	size_t nc_taylor   , 	Base*  taylor      ){		// check assumptions	CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );	CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );	CPPAD_ASSERT_UNKNOWN( i_x < i_z );	CPPAD_ASSERT_UNKNOWN( j < nc_taylor );	// Taylor coefficients corresponding to argument and result	Base* x = taylor + i_x * nc_taylor;	Base* z = taylor + i_z * nc_taylor;	size_t k;	if( j == 0 )		z[j] = sqrt( x[0] );	else	{		CPPAD_ASSERT_KNOWN(			x[0] != Base(0),			"Forward: attempt to take derivatve of square root of zero"		)		z[j] = Base(0);		for(k = 1; k < j; k++)			z[j] -= Base(k) * z[k] * z[j-k];		z[j] /= Base(j);		z[j] += x[j] / Base(2);		z[j] /= z[0];	}}

	/// vector assignment operator	/// If the resulting length of the vector would be more than	/// /c max_length_, and /c NDEBUG is not defined,	/// a CPPAD_ASSERT is generated.	void operator=(		/// right hand size of the assingment operation		const pod_vector& x	)	{	size_t i;		if( x.length_ <= capacity_ )		{	// use existing allocation for this vector			length_ = x.length_;			CPPAD_ASSERT_KNOWN(				length_ <= max_length_ ,				"pod_vector.hpp: attempt to create to large a vector./n"				"If Type is CPPAD_TYPE_ADDR_TYPE, tape long for Type."			);		}		else		{	// free old memory and get new memory of sufficient length			if( capacity_ > 0 )			{	void* v_ptr = reinterpret_cast<void*>( data_ );				if( ! is_pod<Type>() )				{	// call destructor for each element					for(i = 0; i < capacity_; i++)						(data_ + i)->~Type();				}				thread_alloc::return_memory(v_ptr);			}			length_ = capacity_ = 0;			extend( x.length_ );		}		CPPAD_ASSERT_UNKNOWN( length_   == x.length_ );		for(i = 0; i < length_; i++)		{	data_[i] = x.data_[i]; }	}

void index_sort(const VectorKey& keys, VectorSize& ind){	typedef typename VectorKey::value_type Compare;	CheckSimpleVector<size_t, VectorSize>();	typedef index_sort_element<Compare> element;	CPPAD_ASSERT_KNOWN(		size_t(keys.size()) == size_t(ind.size()),		"index_sort: vector sizes do not match"	);	size_t size_work = size_t(keys.size());	size_t size_out;	element* work =		thread_alloc::create_array<element>(size_work, size_out);	// copy initial order into work	size_t i;	for(i = 0; i < size_work; i++)	{	work[i].set_key( keys[i] );		work[i].set_index( i );	}	// sort the work array	std::sort(work, work+size_work);	// copy the indices to the output vector	for(i = 0; i < size_work; i++)		ind[i] = work[i].get_index();	// we are done with this work array	thread_alloc::delete_array(work);	return;}

inline void forward_store_vp_op_0(	size_t         i_z         ,	const addr_t*  arg         , 	size_t         num_par     ,	size_t         nc_taylor   ,	Base*          taylor      ,	size_t         nc_combined ,	bool*          variable    ,	size_t*        combined    ){		CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) <= i_z );	size_t i_vec = Integer( taylor[ arg[1] * nc_taylor + 0 ] );	CPPAD_ASSERT_KNOWN( 		i_vec < combined[ arg[0] - 1 ] ,		"VecAD: index during zero order forward sweep is out of range"	);	CPPAD_ASSERT_UNKNOWN( variable != CPPAD_NULL );	CPPAD_ASSERT_UNKNOWN( combined != CPPAD_NULL );	CPPAD_ASSERT_UNKNOWN( NumArg(StvpOp) == 3 );	CPPAD_ASSERT_UNKNOWN( NumRes(StvpOp) == 0 );	CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );	CPPAD_ASSERT_UNKNOWN( arg[0] + i_vec < nc_combined );	CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );	variable[ arg[0] + i_vec ] = false;	combined[ arg[0] + i_vec ] = arg[2];}

void ADTape<Base>::Independent(VectorAD &x){	// check VectorAD is Simple Vector class with AD<Base> elements	CheckSimpleVector< AD<Base>, VectorAD>();	// dimension of the domain space	size_t n = x.size();	CPPAD_ASSERT_KNOWN(		n > 0,		"Indepdendent: the argument vector x has zero size"	);	CPPAD_ASSERT_UNKNOWN( Rec_.num_rec_var() == 0 );	// mark the beginning of the tape and skip the first variable index 	// (zero) because parameters use taddr zero	CPPAD_ASSERT_NARG_NRES(BeginOp, 0, 1);	Rec_.PutOp(BeginOp);	// place each of the independent variables in the tape	CPPAD_ASSERT_NARG_NRES(InvOp, 0, 1);	size_t j;	for(j = 0; j < n; j++)	{	// tape address for this independent variable		x[j].taddr_ = Rec_.PutOp(InvOp);		x[j].tape_id_    = id_;		CPPAD_ASSERT_UNKNOWN( size_t(x[j].taddr_) == j+1 );		CPPAD_ASSERT_UNKNOWN( Variable(x[j] ) );	}	// done specifying all of the independent variables	size_independent_ = n;}

Vector ADFun<Base>::Jacobian(const Vector &x){	size_t i;	size_t n = Domain();	size_t m = Range();	CPPAD_ASSERT_KNOWN(		size_t(x.size()) == n,		"Jacobian: length of x not equal domain dimension for F"	);	// point at which we are evaluating the Jacobian	Forward(0, x);	// work factor for forward mode	size_t workForward = n;	// work factor for reverse mode	size_t workReverse = 0;	for(i = 0; i < m; i++)	{	if( ! Parameter(i) )			++workReverse;	}	// choose the method with the least work	Vector jac( n * m );	if( workForward <= workReverse )		JacobianFor(*this, x, jac);	else	JacobianRev(*this, x, jac);	return jac;}

Float RombergOne(	Fun           &F , 	const Float   &a , 	const Float   &b , 	size_t         n , 	size_t         p ,	Float         &e ){	size_t ipow2 = 1;	size_t k, i;	Float pow2, sum, x; 	Float  zero  = Float(0);	Float  two   = Float(2);	// check specifications for a NumericType	CheckNumericType<Float>();	CPPAD_ASSERT_KNOWN(		n >= 2,		"RombergOne: n must be greater than or equal 2"	);	CppAD::vector<Float> r(n);	//  set r[i] = trapazoidal rule with 2^i intervals in [a, b]	r[0]  = ( F(a) + F(b) ) * (b - a) / two; 	for(i = 1; i < n; i++)	{	ipow2 *= 2;		// there must be a conversion from int to any numeric type		pow2   = Float(int(ipow2)); 		sum    = zero;		for(k = 1; k < ipow2; k += 2)		{	// start = a + (b-a)/pow2, increment = 2*(b-a)/pow2			x    = ( (pow2 - Float(int(k))) * a + k * b ) / pow2;			sum  = sum + F(x);		}		// combine function evaluations in sum with those in T[i-1]		r[i] = r[i-1] / two + sum * (b - a) / pow2;	}	// now compute the higher order estimates	size_t ipow4    = 1;   // order of accuract for previous estimate	Float pow4, pow4minus;	for(i = 0; i < p; i++)	{	// compute estimate accurate to O[ step^(2*(i+1)) ]		// put resutls in r[n-1], r[n-2], ... , r[n-i+1]		ipow4    *= 4;		pow4      = Float(int(ipow4));		pow4minus = Float(ipow4-1);		for(k = n-1; k > i; k--)			r[k] = ( pow4 * r[k] - r[k-1] ) / pow4minus;	}	// error estimate for r[n]	e = r[n-1] - r[n-2];	if( e < zero )		e = - e;	return r[n-1];}

	/// set row and column for a possibly non-zero element	void set(size_t k, size_t r, size_t c)	{	CPPAD_ASSERT_KNOWN(			k < nnz_,			"The index k is not less than nnz in sparse_rc::set"		);		CPPAD_ASSERT_KNOWN(			r < nr_,			"The index r is not less than nr in sparse_rc::set"		);		CPPAD_ASSERT_KNOWN(			c < nc_,			"The index c is to not less than nc in sparse_rc::set"		);		row_[k] = r;		col_[k] = c;		//	}

inline void forward_powpv_op(	size_t        p           ,	size_t        q           ,	size_t        i_z         ,	const addr_t* arg         ,	const Base*   parameter   ,	size_t        cap_order   ,	Base*         taylor      ){	// convert from final result to first result	i_z -= 2; // 2 = NumRes(PowpvOp) - 1;	// check assumptions	CPPAD_ASSERT_UNKNOWN( NumArg(PowpvOp) == 2 );	CPPAD_ASSERT_UNKNOWN( NumRes(PowpvOp) == 3 );	CPPAD_ASSERT_UNKNOWN( q < cap_order );	CPPAD_ASSERT_UNKNOWN( p <= q );	// Taylor coefficients corresponding to arguments and result	Base* z_0 = taylor + i_z    * cap_order;	// z_0 = log(x)	Base x    = parameter[ arg[0] ];	size_t d;	for(d = p; d <= q; d++)	{	if( d == 0 )			z_0[d] = log(x);		else	z_0[d] = Base(0.0);	}	// 2DO: remove requirement that i_z * cap_order <= max addr_t value	CPPAD_ASSERT_KNOWN(		std::numeric_limits<addr_t>::max() >= i_z * cap_order,		"cppad_tape_addr_type maximum value has been exceeded/n"		"This is due to a kludge in the pow operation and should be fixed."	);	// z_1 = z_0 * y	addr_t adr[2];	// offset of z_i in taylor (as if it were a parameter); i.e., log(x)	adr[0] = addr_t( i_z * cap_order );	// offset of y in taylor (as a variable)	adr[1] = arg[1];	// Trick: use taylor both for the parameter vector and variable values	forward_mulpv_op(p, q, i_z+1, adr, taylor, cap_order, taylor);	// z_2 = exp(z_1)	// zero order case exactly same as Base type operation	if( p == 0 )	{	Base* y   = taylor + arg[1]  * cap_order;		Base* z_2 = taylor + (i_z+2) * cap_order;		z_2[0] = pow(x, y[0]);		p++;	}	if( p <= q )		forward_exp_op(p, q, i_z+2, i_z+1, cap_order, taylor);}

 // -------------------------------------------------------------------- /// constant element access; i.e., we cannot change this element value const Type& operator[](     /// element index, must be less than length and convertable to size_t     size_t i ) const {   CPPAD_ASSERT_KNOWN(         i < length_,         "vector: index greater than or equal vector size"     );     return data_[i]; }