C++ ADFun类代码示例

CppAD::ADFun<T>* CondExp_vpvpFunc(const std::vector<CppAD::AD<T> >& X) {
    using namespace CppAD;
    using namespace std;

    assert(X.size() == 3);

    // parameter value
    AD<T> one = T(1.);
    AD<T> zero = T(0.);

    // dependent variable vector 
    std::vector< AD<T> > Y(5);

    // CondExp(variable, parameter, variable, variable)
    Y[0] = CondExpLt(X[0], one, X[1] * X[2], zero);
    Y[1] = CondExpLe(X[0], one, X[1] * X[2], zero);
    Y[2] = CondExpEq(X[0], one, X[1] * X[2], zero);
    Y[3] = CondExpGe(X[0], one, X[1] * X[2], zero);
    Y[4] = CondExpGt(X[0], one, X[1] * X[2], zero);

    // create f: X -> Y 
    ADFun<T>* fun = new ADFun<T> (X, Y);
    // The option 'no_conditional_skip' is essential in order to avoid an 
    // assertion failure in CppAD and so that branches are not skipped with NDEBUG defined
    fun->optimize();

    return fun;
}

    virtual void SetUp() {
        // use a special object for source code generation
        typedef double Base;
        typedef CG<Base> CGD;
        typedef AD<CGD> ADCG;

        x[0] = 0.5;
        x[1] = 1.5;

        // independent variables
        std::vector<ADCG> u(n);
        for (size_t j = 0; j < n; j++)
            u[j] = x[j];

        CppAD::Independent(u);

        // dependent variable vector 
        std::vector<ADCG> Z(m);

        /**
         * create the CppAD tape as usual
         */
        Z[0] = 1.5 * x[0] + 1;
        Z[1] = 1.0 * x[1] + 2;

        // create f: U -> Z and vectors used for derivative calculations
        _fun = new ADFun<CGD>(u, Z);

        /**
         * Create the dynamic library
         * (generate and compile source code)
         */
        ModelCSourceGen<double> compHelp(*_fun, _modelName);

        compHelp.setCreateForwardZero(true);
        compHelp.setCreateForwardOne(true);
        compHelp.setCreateReverseOne(true);
        compHelp.setCreateReverseTwo(true);
        compHelp.setCreateSparseJacobian(true);
        compHelp.setCreateSparseHessian(true);

        GccCompiler<double> compiler;

        ModelLibraryCSourceGen<double> compDynHelp(compHelp);

        DynamicModelLibraryProcessor<double> p(compDynHelp);

        _dynamicLib = p.createDynamicLibrary(compiler);
        _model = _dynamicLib->model(_modelName);

        // dimensions
        ASSERT_EQ(_model->Domain(), _fun->Domain());
        ASSERT_EQ(_model->Range(), _fun->Range());
    }

	bool FunCheck(
		ADFun<Base>  &f , 
		Fun          &g , 
		const Vector &x , 
		const Base   &r ,
		const Base   &a )
	{	bool ok = true;
	
		size_t m   = f.Range();
		Vector yf  = f.Forward(0, x); 
		Vector yg  = g(x);
	
		size_t i;
		for(i = 0; i < m; i++)
			ok  &= NearEqual(yf[i], yg[i], r, a);
		return ok;
	}

    void testModel(ADFun<CGD>& f,
                   size_t expectedTmp,
                   size_t expectedArraySize) {
        using CppAD::vector;

        size_t n = f.Domain();
        //size_t m = f.Range();

        CodeHandler<double> handler(10 + n * n);

        vector<CGD> indVars(n);
        handler.makeVariables(indVars);

        vector<CGD> dep = f.Forward(0, indVars);

        LanguageC<double> langC("double");
        LangCDefaultVariableNameGenerator<double> nameGen;

        handler.generateCode(std::cout, langC, dep, nameGen);

        ASSERT_EQ(handler.getTemporaryVariableCount(), expectedTmp);
        ASSERT_EQ(handler.getTemporaryArraySize(), expectedArraySize);
    }

TEST_F(CppADCGEvaluatorTest, Atomic) {
    using ADCG = AD<CGD>;

    std::vector<AD<double>> ax(2);
    std::vector<AD<double>> ay(3);
    for (size_t j = 0; j < ax.size(); j++) {
        ax[j] = j + 2;
    }

    checkpoint<double> atomicFun("func", testModel, ax, ay); // the normal atomic function
    CGAtomicFun<double> atomic(atomicFun, ax); // a wrapper used to tape with CG<Base>

    ModelType model = [&](const std::vector<CGD>& x) {
        // independent variables
        std::vector<ADCG> ax(x.size());
        for (size_t j = 0; j < ax.size(); j++) {
            ax[j] = x[j];
        }

        CppAD::Independent(ax);

        // dependent variable vector
        std::vector<ADCG> ay(3);

        atomic(ax, ay);

        ADFun<CGD> fun;
        fun.Dependent(ay);

        std::vector<CGD> y = fun.Forward(0, x);

        return y;
    };

    this->testCG(model, std::vector<double>{0.5, 1.5});
}

void JacobianRev(ADFun<Base> &f, const Vector &x, Vector &jac)
{	size_t i;
	size_t j;

	size_t n = f.Domain();
	size_t m = f.Range();

	CPPAD_ASSERT_UNKNOWN( size_t(x.size())   == f.Domain() );
	CPPAD_ASSERT_UNKNOWN( size_t(jac.size()) == f.Range() * f.Domain() );

	// argument and result for reverse mode calculations
	Vector u(n);
	Vector v(m);

	// initialize all the components
	for(i = 0; i < m; i++)
		v[i] = Base(0);

	// loop through the different coordinate directions
	for(i = 0; i < m; i++)
	{	if( f.Parameter(i) )
		{	// return zero for this component of f
			for(j = 0; j < n; j++)
				jac[ i * n + j ] = Base(0);
		}
		else
		{
			// set v to the i-th coordinate direction
			v[i] = Base(1);

			// compute the derivative of this component of f
			u = f.Reverse(1, v);

			// reset v to vector of all zeros
			v[i] = Base(0);

			// return the result
			for(j = 0; j < n; j++)
				jac[ i * n + j ] = u[j];
		}
	}
}

void JacobianFor(ADFun<Base> &f, const Vector &x, Vector &jac)
{	size_t i;
	size_t j;

	size_t n = f.Domain();
	size_t m = f.Range();

	// check Vector is Simple Vector class with Base type elements
	CheckSimpleVector<Base, Vector>();

	CPPAD_ASSERT_UNKNOWN( size_t(x.size())   == f.Domain() );
	CPPAD_ASSERT_UNKNOWN( size_t(jac.size()) == f.Range() * f.Domain() );

	// argument and result for forward mode calculations
	Vector u(n);
	Vector v(m);

	// initialize all the components
	for(j = 0; j < n; j++)
		u[j] = Base(0);

	// loop through the different coordinate directions
	for(j = 0; j < n; j++)
	{	// set u to the j-th coordinate direction
		u[j] = Base(1);

		// compute the partial of f w.r.t. this coordinate direction
		v = f.Forward(1, u);

		// reset u to vector of all zeros
		u[j] = Base(0);

		// return the result
		for(i = 0; i < m; i++)
			jac[ i * n + j ] = v[i];
	}
}

    void createAtomicLib() {
        typedef CG<double> CGD;
        typedef AD<CGD> ADCGD;

        /**
         * Tape model
         */
        std::vector<ADCGD> xa(na_);
        for (size_t j = 0; j < na_; j++)
            xa[j] = xa_[j];
        CppAD::Independent(xa);

        std::vector<ADCGD> ya(ns_);

        atomicFunction(xa, ya);

        ADFun<CGD> fun;
        fun.Dependent(ya);

        /**
         * Compile
         */
        std::string lName = getAtomicLibName()+(ignoreParameters_ ? "" : "All");
        ModelCSourceGen<double> compHelpL(fun, lName);
        compHelpL.setCreateForwardZero(true);
        compHelpL.setCreateForwardOne(true);
        compHelpL.setCreateReverseOne(true);
        compHelpL.setCreateReverseTwo(true);
        compHelpL.setTypicalIndependentValues(xa_);

        if (ignoreParameters_) {
            std::vector<std::set<size_t> > jacSparAll = extra::jacobianSparsitySet<std::vector<std::set<size_t> > >(fun);
            std::vector<std::set<size_t> > jacSpar(jacSparAll.size());
            for (size_t i = 0; i < jacSparAll.size(); i++) {
                // only differential information for states and controls
                std::set<size_t>::const_iterator itEnd = jacSparAll[i].upper_bound(ns_ + nm_ - 1);
                if (itEnd != jacSparAll[i].begin())
                    jacSpar[i].insert(jacSparAll[i].begin(), itEnd);
            }
            compHelpL.setCustomSparseJacobianElements(jacSpar);

            std::vector<std::set<size_t> > hessSparAll = extra::hessianSparsitySet<std::vector<std::set<size_t> > >(fun);
            std::vector<std::set<size_t> > hessSpar(hessSparAll.size());
            for (size_t i = 0; i < ns_ + nm_; i++) {
                std::set<size_t>::const_iterator it = hessSparAll[i].upper_bound(i); // only the lower left side
                if (it != hessSparAll[i].begin())
                    hessSpar[i].insert(hessSparAll[i].begin(), it);
            }
            compHelpL.setCustomSparseHessianElements(hessSpar);
        }

        ModelLibraryCSourceGen<double> compDynHelpL(compHelpL);
        compDynHelpL.setVerbose(verbose_);

        SaveFilesModelLibraryProcessor<double>::saveLibrarySourcesTo(compDynHelpL, "sources_" + lName);

        DynamicModelLibraryProcessor<double> p(compDynHelpL, lName);
        GccCompiler<double> compiler;
        if (!compilerFlags_.empty())
            compiler.setCompileFlags(compilerFlags_);
        atomicDynamicLib_.reset(p.createDynamicLibrary(compiler));

        /**
         * load the model
         */
        atomicModel_.reset(atomicDynamicLib_->model(lName));
    }

bool old_usead_1(void)
{	bool ok = true;
	using CppAD::NearEqual;
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();

	// --------------------------------------------------------------------
	// Create the ADFun<doulbe> r_
	create_r();

	// --------------------------------------------------------------------
	// Create the function f(x)
	//
	// domain space vector
	size_t n  = 1;
	double  x0 = 0.5;
	vector< AD<double> > ax(n);
	ax[0]     = x0;

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

	// range space vector
	size_t m = 1;
	vector< AD<double> > ay(m);

	// call user function and store reciprocal(x) in au[0]
	vector< AD<double> > au(m);
	size_t id = 0;           // not used
	reciprocal(id, ax, au);	// u = 1 / x

	// call user function and store reciprocal(u) in ay[0]
	reciprocal(id, au, ay);	// y = 1 / u = x

	// create f: x -> y and stop tape recording
	ADFun<double> f;
	f.Dependent(ax, ay);  // f(x) = x

	// --------------------------------------------------------------------
	// Check function value results
	//
	// check function value
	double check = x0;
	ok &= NearEqual( Value(ay[0]) , check,  eps, eps);

	// check zero order forward mode
	size_t q;
	vector<double> x_q(n), y_q(m);
	q      = 0;
	x_q[0] = x0;
	y_q    = f.Forward(q, x_q);
	ok &= NearEqual(y_q[0] , check,  eps, eps);

	// check first order forward mode
	q      = 1;
	x_q[0] = 1;
	y_q    = f.Forward(q, x_q);
	check  = 1.;
	ok &= NearEqual(y_q[0] , check,  eps, eps);

	// check second order forward mode
	q      = 2;
	x_q[0] = 0;
	y_q    = f.Forward(q, x_q);
	check  = 0.;
	ok &= NearEqual(y_q[0] , check,  eps, eps);

	// --------------------------------------------------------------------
	// Check reverse mode results
	//
	// third order reverse mode
	q     = 3;
	vector<double> w(m), dw(n * q);
	w[0]  = 1.;
	dw    = f.Reverse(q, w);
	check = 1.;
	ok &= NearEqual(dw[0] , check,  eps, eps);
	check = 0.;
	ok &= NearEqual(dw[1] , check,  eps, eps);
	ok &= NearEqual(dw[2] , check,  eps, eps);

	// --------------------------------------------------------------------
	// forward mode sparstiy pattern
	size_t p = n;
	CppAD::vectorBool r1(n * p), s1(m * p);
	r1[0] = true;          // compute sparsity pattern for x[0]
	s1    = f.ForSparseJac(p, r1);
	ok  &= s1[0] == true;  // f[0] depends on x[0]

	// --------------------------------------------------------------------
	// reverse mode sparstiy pattern
	q = m;
	CppAD::vectorBool s2(q * m), r2(q * n);
	s2[0] = true;          // compute sparsity pattern for f[0]
	r2    = f.RevSparseJac(q, s2);
	ok  &= r2[0] == true;  // f[0] depends on x[0]

	// --------------------------------------------------------------------
	// Hessian sparsity (using previous ForSparseJac call)
	CppAD::vectorBool s3(m), h(p * n);
	s3[0] = true;        // compute sparsity pattern for f[0]
//.........这里部分代码省略.........

    void testDynamicFull(std::vector<ADCG>& u,
                         const std::vector<double>& x,
                         const std::vector<double>& xNorm,
                         const std::vector<double>& eqNorm,
                         size_t maxAssignPerFunc = 100,
                         double epsilonR = 1e-14,
                         double epsilonA = 1e-14) {
        ASSERT_EQ(u.size(), x.size());
        ASSERT_EQ(x.size(), xNorm.size());

        using namespace std;

        // use a special object for source code generation
        CppAD::Independent(u);

        for (size_t i = 0; i < u.size(); i++)
            u[i] *= xNorm[i];

        // dependent variable vector 
        std::vector<ADCG> Z = model(u);

        if (eqNorm.size() > 0) {
            ASSERT_EQ(Z.size(), eqNorm.size());
            for (size_t i = 0; i < Z.size(); i++)
                Z[i] /= eqNorm[i];
        }

        /**
         * create the CppAD tape as usual
         */
        // create f: U -> Z and vectors used for derivative calculations
        ADFun<CGD> fun;
        fun.Dependent(Z);

        /**
         * Create the dynamic library
         * (generate and compile source code)
         */
        ModelCSourceGen<double> compHelp(fun, _name + "dynamic");

        compHelp.setCreateForwardZero(true);
        compHelp.setCreateJacobian(true);
        compHelp.setCreateHessian(true);
        compHelp.setCreateSparseJacobian(true);
        compHelp.setCreateSparseHessian(true);
        compHelp.setCreateForwardOne(true);
        compHelp.setCreateReverseOne(true);
        compHelp.setCreateReverseTwo(true);
        compHelp.setMaxAssignmentsPerFunc(maxAssignPerFunc);

        ModelLibraryCSourceGen<double> compDynHelp(compHelp);

        SaveFilesModelLibraryProcessor<double>::saveLibrarySourcesTo(compDynHelp, "sources_" + _name + "_1");

        DynamicModelLibraryProcessor<double> p(compDynHelp);
        GccCompiler<double> compiler;
        DynamicLib<double>* dynamicLib = p.createDynamicLibrary(compiler);

        /**
         * test the library
         */
        GenericModel<double>* model = dynamicLib->model(_name + "dynamic");
        ASSERT_TRUE(model != nullptr);

        testModelResults(*model, fun, x, epsilonR, epsilonA);

        delete model;
        delete dynamicLib;
    }

bool old_usead_2(void)
{	bool ok = true;
	using CppAD::NearEqual;
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();

	// --------------------------------------------------------------------
	// Create the ADFun<doulbe> r_
	create_r();

	// --------------------------------------------------------------------
	// domain and range space vectors
	size_t n = 3, m = 2;
	vector< AD<double> > au(n), ax(n), ay(m);
	au[0]         = 0.0;        // value of z_0 (t) = t, at t = 0
	ax[1]         = 0.0;        // value of z_1 (t) = t^2/2, at t = 0
	au[2]         = 1.0;        // final t
	CppAD::Independent(au);
	size_t M      = 2;          // number of r steps to take
	ax[0]         = au[0];      // value of z_0 (t) = t, at t = 0
	ax[1]         = au[1];      // value of z_1 (t) = t^2/2, at t = 0
	AD<double> dt = au[2] / double(M);  // size of each r step
	ax[2]         = dt;
	for(size_t i_step = 0; i_step < M; i_step++)
	{	size_t id = 0;               // not used
		solve_ode(id, ax, ay);
		ax[0] = ay[0];
		ax[1] = ay[1];
	}

	// create f: u -> y and stop tape recording
	// y_0(t) = u_0 + t                   = u_0 + u_2
	// y_1(t) = u_1 + u_0 * t + t^2 / 2   = u_1 + u_0 * u_2 + u_2^2 / 2
	// where t = u_2
	ADFun<double> f;
	f.Dependent(au, ay);

	// --------------------------------------------------------------------
	// Check forward mode results
	//
	// zero order forward
	vector<double> up(n), yp(m);
	size_t q  = 0;
	double u0 = 0.5;
	double u1 = 0.25;
	double u2 = 0.75;
	double check;
	up[0]     = u0;
	up[1]     = u1;
	up[2]     = u2;
	yp        = f.Forward(q, up);
	check     = u0 + u2;
	ok       &= NearEqual( yp[0], check,  eps, eps);
	check     = u1 + u0 * u2 + u2 * u2 / 2.0;
	ok       &= NearEqual( yp[1], check,  eps, eps);
	//
	// forward mode first derivative w.r.t t
	q         = 1;
	up[0]     = 0.0;
	up[1]     = 0.0;
	up[2]     = 1.0;
	yp        = f.Forward(q, up);
	check     = 1.0;
	ok       &= NearEqual( yp[0], check,  eps, eps);
	check     = u0 + u2;
	ok       &= NearEqual( yp[1], check,  eps, eps);
	//
	// forward mode second order Taylor coefficient w.r.t t
	q         = 2;
	up[0]     = 0.0;
	up[1]     = 0.0;
	up[2]     = 0.0;
	yp        = f.Forward(q, up);
	check     = 0.0;
	ok       &= NearEqual( yp[0], check,  eps, eps);
	check     = 1.0 / 2.0;
	ok       &= NearEqual( yp[1], check,  eps, eps);
	// --------------------------------------------------------------------
	// reverse mode derivatives of \partial_t y_1 (t)
	vector<double> w(m * q), dw(n * q);
	w[0 * q + 0]  = 0.0;
	w[1 * q + 0]  = 0.0;
	w[0 * q + 1]  = 0.0;
	w[1 * q + 1]  = 1.0;
	dw        = f.Reverse(q, w);
	// derivative of y_1(u) = u_1 + u_0 * u_2 + u_2^2 / 2,  w.r.t. u
	// is equal deritative of \partial_u2 y_1(u) w.r.t \partial_u2 u
	check     = u2;
	ok       &= NearEqual( dw[0 * q + 1], check,  eps, eps);
	check     = 1.0;
	ok       &= NearEqual( dw[1 * q + 1], check,  eps, eps);
	check     = u0 + u2;
	ok       &= NearEqual( dw[2 * q + 1], check,  eps, eps);
	// derivative of \partial_t y_1 w.r.t u = u_0 + t,  w.r.t u
	check     = 1.0;
	ok       &= NearEqual( dw[0 * q + 0], check,  eps, eps);
	check     = 0.0;
	ok       &= NearEqual( dw[1 * q + 0], check,  eps, eps);
	check     = 1.0;
	ok       &= NearEqual( dw[2 * q + 0], check,  eps, eps);
	// --------------------------------------------------------------------
//.........这里部分代码省略.........

void BenderQuad(
	const BAvector   &x     , 
	const BAvector   &y     , 
	Fun               fun   , 
	BAvector         &g     ,
	BAvector         &gx    ,
	BAvector         &gxx   )
{	// determine the base type
	typedef typename BAvector::value_type Base;

	// check that BAvector is a SimpleVector class
	CheckSimpleVector<Base, BAvector>();

	// declare the ADvector type
	typedef CPPAD_TESTVECTOR(AD<Base>) ADvector;

	// size of the x and y spaces
	size_t n = size_t(x.size());
	size_t m = size_t(y.size());

	// check the size of gx and gxx
	CPPAD_ASSERT_KNOWN(
		g.size() == 1,
		"BenderQuad: size of the vector g is not equal to 1"
	);
	CPPAD_ASSERT_KNOWN(
		size_t(gx.size()) == n,
		"BenderQuad: size of the vector gx is not equal to n"
	);
	CPPAD_ASSERT_KNOWN(
		size_t(gxx.size()) == n * n,
		"BenderQuad: size of the vector gxx is not equal to n * n"
	);

	// some temporary indices
	size_t i, j;

	// variable versions x
	ADvector vx(n);
	for(j = 0; j < n; j++)
		vx[j] = x[j];
	
	// declare the independent variables
	Independent(vx);

	// evaluate h = H(x, y) 
	ADvector h(m);
	h = fun.h(vx, y);

	// evaluate dy (x) = Newton step as a function of x through h only
	ADvector dy(m);
	dy = fun.dy(x, y, h);

	// variable version of y
	ADvector vy(m);
	for(j = 0; j < m; j++)
		vy[j] = y[j] + dy[j];

	// evaluate G~ (x) = F [ x , y + dy(x) ] 
	ADvector gtilde(1);
	gtilde = fun.f(vx, vy);

	// AD function object that corresponds to G~ (x)
	// We will make heavy use of this tape, so optimize it
	ADFun<Base> Gtilde;
	Gtilde.Dependent(vx, gtilde); 
	Gtilde.optimize();

	// value of G(x)
	g = Gtilde.Forward(0, x);

	// initial forward direction vector as zero
	BAvector dx(n);
	for(j = 0; j < n; j++)
		dx[j] = Base(0);

	// weight, first and second order derivative values
	BAvector dg(1), w(1), ddw(2 * n);
	w[0] = 1.;


	// Jacobian and Hessian of G(x) is equal Jacobian and Hessian of Gtilde
	for(j = 0; j < n; j++)
	{	// compute partials in x[j] direction
		dx[j] = Base(1);
		dg    = Gtilde.Forward(1, dx);
		gx[j] = dg[0];

		// restore the dx vector to zero
		dx[j] = Base(0);

		// compute second partials w.r.t x[j] and x[l]  for l = 1, n
		ddw = Gtilde.Reverse(2, w);
		for(i = 0; i < n; i++)
			gxx[ i * n + j ] = ddw[ i * 2 + 1 ];
	}

	return;
}

void ADFun<Base>::operator=(const ADFun<Base>& f)
{	size_t m = f.Range();
	size_t n = f.Domain();
	size_t i;

	// go through member variables in ad_fun.hpp order
	// 
	// size_t objects
	has_been_optimized_        = f.has_been_optimized_;
	check_for_nan_             = f.check_for_nan_;
	compare_change_count_      = f.compare_change_count_;
	compare_change_number_     = f.compare_change_number_;
	compare_change_op_index_   = f.compare_change_op_index_;
	num_order_taylor_          = f.num_order_taylor_;
	cap_order_taylor_          = f.cap_order_taylor_;
	num_direction_taylor_      = f.num_direction_taylor_;
	num_var_tape_              = f.num_var_tape_;
	//
	// CppAD::vector objects
	ind_taddr_.resize(n);
	ind_taddr_                 = f.ind_taddr_;
	dep_taddr_.resize(m);
	dep_taddr_                 = f.dep_taddr_;
	dep_parameter_.resize(m);
	dep_parameter_             = f.dep_parameter_;
	//
	// pod_vector objects
	taylor_                    = f.taylor_;
	cskip_op_                  = f.cskip_op_;
	load_op_                   = f.load_op_;
	//
	// player
	play_                      = f.play_;
	//
	// sparse_pack
	for_jac_sparse_pack_.resize(0, 0);
	size_t n_set = f.for_jac_sparse_pack_.n_set();
	size_t end   = f.for_jac_sparse_pack_.end();
	if( n_set > 0 )
	{	CPPAD_ASSERT_UNKNOWN( n_set == num_var_tape_  );
		CPPAD_ASSERT_UNKNOWN( f.for_jac_sparse_set_.n_set() == 0 );
		for_jac_sparse_pack_.resize(n_set, end);
		for(i = 0; i < num_var_tape_  ; i++)
		{	for_jac_sparse_pack_.assignment(
				i                       ,
				i                       ,
				f.for_jac_sparse_pack_
			);
		}
	}
	//
	// sparse_set
	for_jac_sparse_set_.resize(0, 0);
	n_set = f.for_jac_sparse_set_.n_set();
	end   = f.for_jac_sparse_set_.end();
	if( n_set > 0 )
	{	CPPAD_ASSERT_UNKNOWN( n_set == num_var_tape_  );
		CPPAD_ASSERT_UNKNOWN( f.for_jac_sparse_pack_.n_set() == 0 );
		for_jac_sparse_set_.resize(n_set, end);
		for(i = 0; i < num_var_tape_; i++)
		{	for_jac_sparse_set_.assignment(
				i                       ,
				i                       ,
				f.for_jac_sparse_set_
			);
		}
	}
}

void ADFun<Base>::operator=(const ADFun<Base>& f)
{   size_t m = f.Range();
    size_t n = f.Domain();

    // go through member variables in order
    // (see ad_fun.hpp for meaning of each variable)
    compare_change_            = 0;
    taylor_per_var_            = 0;
    taylor_col_dim_            = 0;
    total_num_var_             = f.total_num_var_;
    ind_taddr_.resize(n);
    ind_taddr_                 = f.ind_taddr_;
    dep_taddr_.resize(m);
    dep_taddr_                 = f.dep_taddr_;
    dep_parameter_.resize(m);
    dep_parameter_             = f.dep_parameter_;
    play_                      = f.play_;
    if( taylor_ != CPPAD_NULL )
        CPPAD_TRACK_DEL_VEC(taylor_);
    taylor_                    = CPPAD_NULL;
    for_jac_sparse_pack_.resize(0, 0);
    for_jac_sparse_set_.resize(0, 0);

    // allocate and copy the Taylor coefficients
    taylor_per_var_     = f.taylor_per_var_;
    taylor_col_dim_     = f.taylor_col_dim_;
    size_t length       = total_num_var_ * taylor_col_dim_;
    if( length > 0 ) taylor_   = CPPAD_TRACK_NEW_VEC(length, taylor_);
    size_t i, j;
    for(i = 0; i < total_num_var_; i++)
    {   for(j = 0; j < taylor_per_var_; j++)
        {   taylor_[ i * taylor_col_dim_ + j ] =
                f.taylor_[ i * taylor_col_dim_ + j ];
        }
    }

    // allocate and copy the forward sparsity information
    size_t n_set = f.for_jac_sparse_pack_.n_set();
    size_t end   = f.for_jac_sparse_pack_.end();
    if( n_set > 0 )
    {   CPPAD_ASSERT_UNKNOWN( n_set == total_num_var_ );
        CPPAD_ASSERT_UNKNOWN( f.for_jac_sparse_set_.n_set() == 0 );
        for_jac_sparse_pack_.resize(n_set, end);
        for(i = 0; i < total_num_var_ ; i++)
        {   for_jac_sparse_pack_.assignment(
                i                       ,
                i                       ,
                f.for_jac_sparse_pack_
            );
        }
    }
    n_set = f.for_jac_sparse_set_.n_set();
    end   = f.for_jac_sparse_set_.end();
    if( n_set > 0 )
    {   CPPAD_ASSERT_UNKNOWN( n_set == total_num_var_ );
        CPPAD_ASSERT_UNKNOWN( f.for_jac_sparse_pack_.n_set() == 0 );
        for_jac_sparse_set_.resize(n_set, end);
        for(i = 0; i < total_num_var_; i++)
        {   for_jac_sparse_set_.assignment(
                i                       ,
                i                       ,
                f.for_jac_sparse_set_
            );
        }
    }

}

	/*!
 	Link from user_atomic to forward sparse Jacobian 

	\copydetails atomic_base::rev_sparse_hes
 	*/
	virtual bool rev_sparse_hes(
		const vector<bool>&                     vx ,
		const vector<bool>&                     s  ,
		      vector<bool>&                     t  ,
		size_t                                  q  ,
		const vector< std::set<size_t> >&       r  ,
		const vector< std::set<size_t> >&       u  ,
		      vector< std::set<size_t> >&       v  )
	{	size_t n       = v.size();
		size_t m       = u.size();
		CPPAD_ASSERT_UNKNOWN( r.size() == v.size() );
		CPPAD_ASSERT_UNKNOWN( s.size() == m );
		CPPAD_ASSERT_UNKNOWN( t.size() == n );
		bool ok        = true;
		bool transpose = true;
		std::set<size_t>::const_iterator itr;

		// compute sparsity pattern for T(x) = S(x) * f'(x)
		t = f_.RevSparseJac(1, s);
# ifndef NDEBUG
		for(size_t j = 0; j < n; j++)
			CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
# endif

		// V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R 
		// U(x) = g''(y) * f'(x) * R
		// S(x) = g'(y)
		
		// compute sparsity pattern for A(x) = f'(x)^T * U(x)
		vector< std::set<size_t> > a(n);
		a = f_.RevSparseJac(q, u, transpose);

		// set version of s
		vector< std::set<size_t> > set_s(1);
		CPPAD_ASSERT_UNKNOWN( set_s[0].empty() );
		size_t i;
		for(i = 0; i < m; i++)
			if( s[i] )
				set_s[0].insert(i);

		// compute sparsity pattern for H(x) = (S(x) * F)''(x) * R
		// (store it in v)
		f_.ForSparseJac(q, r);
		v = f_.RevSparseHes(q, set_s, transpose);

		// compute sparsity pattern for V(x) = A(x) + H(x)
		for(i = 0; i < n; i++)
		{	for(itr = a[i].begin(); itr != a[i].end(); itr++)
			{	size_t j = *itr;
				CPPAD_ASSERT_UNKNOWN( j < q );
				v[i].insert(j);
			}
		}

		// no longer need the forward mode sparsity pattern
		// (have to reconstruct them every time)
		f_.size_forward_set(0);

		return ok;
	}