// Copyright (C) 2016-2018 T. Zachary Laine
//
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#include "autodiff.h"

#include <iostream>

#include <boost/yap/algorithm.hpp>
#include <boost/polymorphic_cast.hpp>
#include <boost/hana/for_each.hpp>

#define BOOST_TEST_MODULE autodiff_test
#include <boost/test/included/unit_test.hpp>


double const Epsilon = 10.0e-6;
#define CHECK_CLOSE(A,B) do { BOOST_CHECK_CLOSE(A,B,Epsilon); } while(0)

using namespace AutoDiff;

//[ autodiff_expr_template_decl
template <boost::yap::expr_kind Kind, typename Tuple>
struct autodiff_expr
{
    static boost::yap::expr_kind const kind = Kind;

    Tuple elements;
};

BOOST_YAP_USER_UNARY_OPERATOR(negate, autodiff_expr, autodiff_expr)

BOOST_YAP_USER_BINARY_OPERATOR(plus, autodiff_expr, autodiff_expr)
BOOST_YAP_USER_BINARY_OPERATOR(minus, autodiff_expr, autodiff_expr)
BOOST_YAP_USER_BINARY_OPERATOR(multiplies, autodiff_expr, autodiff_expr)
BOOST_YAP_USER_BINARY_OPERATOR(divides, autodiff_expr, autodiff_expr)
//]

//[ autodiff_expr_literals_decl
namespace autodiff_placeholders {

    // This defines a placeholder literal operator that creates autodiff_expr
    // placeholders.
    BOOST_YAP_USER_LITERAL_PLACEHOLDER_OPERATOR(autodiff_expr)

}
//]

//[ autodiff_function_terminals
template <OPCODE Opcode>
struct autodiff_fn_expr :
    autodiff_expr<boost::yap::expr_kind::terminal, boost::hana::tuple<OPCODE>>
{
    autodiff_fn_expr () :
        autodiff_expr {boost::hana::tuple<OPCODE>{Opcode}}
    {}

    BOOST_YAP_USER_CALL_OPERATOR_N(::autodiff_expr, 1);
};

// Someone included <math.h>, so we have to add trailing underscores.
autodiff_fn_expr<OP_SIN> const sin_;
autodiff_fn_expr<OP_COS> const cos_;
autodiff_fn_expr<OP_SQRT> const sqrt_;
//]

//[ autodiff_xform
struct xform
{
    // Create a var-node for each placeholder when we see it for the first
    // time.
    template <long long I>
    Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
                       boost::yap::placeholder<I>)
    {
        if (list_.size() < I)
            list_.resize(I);
        auto & retval = list_[I - 1];
        if (retval == nullptr)
            retval = create_var_node();
        return retval;
    }

    // Create a param-node for every numeric terminal in the expression.
    Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>, double x)
    { return create_param_node(x); }

    // Create a "uary" node for each call expression, using its OPCODE.
    template <typename Expr>
    Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::call>,
                       OPCODE opcode, Expr const & expr)
    {
        return create_uary_op_node(
            opcode,
            boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr), *this)
        );
    }

    template <typename Expr>
    Node * operator() (boost::yap::expr_tag<boost::yap::expr_kind::negate>,
                       Expr const & expr)
    {
        return create_uary_op_node(
            OP_NEG,
            boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr), *this)
        );
    }

    // Define a mapping from binary arithmetic expr_kind to OPCODE...
    static OPCODE op_for_kind (boost::yap::expr_kind kind)
    {
        switch (kind) {
        case boost::yap::expr_kind::plus: return OP_PLUS;
        case boost::yap::expr_kind::minus: return OP_MINUS;
        case boost::yap::expr_kind::multiplies: return OP_TIMES;
        case boost::yap::expr_kind::divides: return OP_DIVID;
        default: assert(!"This should never execute"); return OPCODE{};
        }
        assert(!"This should never execute");
        return OPCODE{};
    }

    // ... and use it to handle all the binary arithmetic operators.
    template <boost::yap::expr_kind Kind, typename Expr1, typename Expr2>
    Node * operator() (boost::yap::expr_tag<Kind>, Expr1 const & expr1, Expr2 const & expr2)
    {
        return create_binary_op_node(
            op_for_kind(Kind),
            boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr1), *this),
            boost::yap::transform(boost::yap::as_expr<autodiff_expr>(expr2), *this)
        );
    }

    vector<Node *> & list_;
};
//]

//[ autodiff_to_node
template <typename Expr, typename ...T>
Node * to_auto_diff_node (Expr const & expr, vector<Node *> & list, T ... args)
{
    Node * retval = nullptr;

    // This fills in list as a side effect.
    retval = boost::yap::transform(expr, xform{list});

    assert(list.size() == sizeof...(args));

    // Fill in the values of the value-nodes in list with the "args"
    // parameter pack.
    auto it = list.begin();
    boost::hana::for_each(
        boost::hana::make_tuple(args ...),
        [&it](auto x) {
            Node * n = *it;
            VNode * v = boost::polymorphic_downcast<VNode *>(n);
            v->val = x;
            ++it;
        }
    );

    return retval;
}
//]

struct F{
	F() {	AutoDiff::autodiff_setup();	}
	~F(){	AutoDiff::autodiff_cleanup();	}
};


BOOST_FIXTURE_TEST_SUITE(all, F)

//[ autodiff_original_node_builder
Node* build_linear_fun1_manually(vector<Node*>& list)
{
	//f(x1,x2,x3) = -5*x1+sin(10)*x1+10*x2-x3/6
	PNode* v5 = create_param_node(-5);
	PNode* v10 = create_param_node(10);
	PNode* v6 = create_param_node(6);
	VNode*	x1 = create_var_node();
	VNode*	x2 = create_var_node();
	VNode*	x3 = create_var_node();

	OPNode* op1 = create_binary_op_node(OP_TIMES,v5,x1);   //op1 = v5*x1
	OPNode* op2 = create_uary_op_node(OP_SIN,v10);         //op2 = sin(v10)
	OPNode* op3 = create_binary_op_node(OP_TIMES,op2,x1);  //op3 = op2*x1
	OPNode* op4 = create_binary_op_node(OP_PLUS,op1,op3);  //op4 = op1 + op3
	OPNode*	op5 = create_binary_op_node(OP_TIMES,v10,x2);  //op5 = v10*x2
	OPNode* op6 = create_binary_op_node(OP_PLUS,op4,op5);  //op6 = op4+op5
	OPNode* op7 = create_binary_op_node(OP_DIVID,x3,v6);   //op7 = x3/v6
	OPNode* op8 = create_binary_op_node(OP_MINUS,op6,op7); //op8 = op6 - op7
	x1->val = -1.9;
	x2->val = 2;
	x3->val = 5./6.;
	list.push_back(x1);
	list.push_back(x2);
	list.push_back(x3);
	return op8;
}
//]

//[ autodiff_yap_node_builder
Node* build_linear_fun1(vector<Node*>& list)
{
    //f(x1,x2,x3) = -5*x1+sin(10)*x1+10*x2-x3/6
    using namespace autodiff_placeholders;
    return to_auto_diff_node(
        -5 * 1_p + sin_(10) * 1_p + 10 * 2_p - 3_p / 6,
        list,
        -1.9,
        2,
        5./6.
    );
}
//]

Node* build_linear_function2_manually(vector<Node*>& list)
{
	//f(x1,x2,x3) = -5*x1+-10*x1+10*x2-x3/6
	PNode* v5 = create_param_node(-5);
	PNode* v10 = create_param_node(10);
	PNode* v6 = create_param_node(6);
	VNode*	x1 = create_var_node();
	VNode*	x2 = create_var_node();
	VNode*	x3 = create_var_node();
	list.push_back(x1);
	list.push_back(x2);
	list.push_back(x3);
	OPNode* op1 = create_binary_op_node(OP_TIMES,v5,x1); //op1 = v5*x1
	OPNode* op2 = create_uary_op_node(OP_NEG,v10);       //op2 = -v10
	OPNode* op3 = create_binary_op_node(OP_TIMES,op2,x1);//op3 = op2*x1
	OPNode* op4 = create_binary_op_node(OP_PLUS,op1,op3);//op4 = op1 + op3
	OPNode*	op5 = create_binary_op_node(OP_TIMES,v10,x2);//op5 = v10*x2
	OPNode* op6 = create_binary_op_node(OP_PLUS,op4,op5);//op6 = op4+op5
	OPNode* op7 = create_binary_op_node(OP_DIVID,x3,v6); //op7 = x3/v6
	OPNode* op8 = create_binary_op_node(OP_MINUS,op6,op7);//op8 = op6 - op7
	x1->val = -1.9;
	x2->val = 2;
	x3->val = 5./6.;
	return op8;
}

Node* build_linear_function2(vector<Node*>& list)
{
    //f(x1,x2,x3) = -5*x1+-10*x1+10*x2-x3/6
    using namespace autodiff_placeholders;
    auto ten = boost::yap::make_terminal<autodiff_expr>(10);
    return to_auto_diff_node(
        -5 * 1_p + -ten * 1_p + 10 * 2_p - 3_p / 6,
        list,
        -1.9,
        2,
        5./6.
    );
}

Node* build_nl_function1_manually(vector<Node*>& list)
{
//	(x1*x2 * sin(x1))/x3 + x2*x4 - x1/x2
	VNode* x1 = create_var_node();
	VNode* x2 = create_var_node();
	VNode* x3 = create_var_node();
	VNode* x4 = create_var_node();
	x1->val = -1.23;
	x2->val = 7.1231;
	x3->val = 2;
	x4->val = -10;
	list.push_back(x1);
	list.push_back(x2);
	list.push_back(x3);
	list.push_back(x4);

	OPNode* op1 = create_binary_op_node(OP_TIMES,x2,x1);
	OPNode* op2 = create_uary_op_node(OP_SIN,x1);
	OPNode* op3 = create_binary_op_node(OP_TIMES,op1,op2);
	OPNode* op4 = create_binary_op_node(OP_DIVID,op3,x3);
	OPNode* op5 = create_binary_op_node(OP_TIMES,x2,x4);
	OPNode* op6 = create_binary_op_node(OP_PLUS,op4,op5);
	OPNode* op7 = create_binary_op_node(OP_DIVID,x1,x2);
	OPNode* op8 = create_binary_op_node(OP_MINUS,op6,op7);
	return op8;
}

Node* build_nl_function1(vector<Node*>& list)
{
    // (x1*x2 * sin(x1))/x3 + x2*x4 - x1/x2
    using namespace autodiff_placeholders;
    return to_auto_diff_node(
        (1_p * 2_p * sin_(1_p)) / 3_p + 2_p * 4_p - 1_p / 2_p,
        list,
	-1.23,
	7.1231,
	2,
	-10
    );
}

BOOST_AUTO_TEST_CASE( test_linear_fun1 )
{
	BOOST_TEST_MESSAGE("test_linear_fun1");
	vector<Node*> list;
	Node* root = build_linear_fun1(list);
	vector<double> grad;
	double val1 = grad_reverse(root,list,grad);
	double val2 = eval_function(root);
	double x1g[] = {-5.5440211108893697744548489936278,10.0,-0.16666666666666666666666666666667};

	for(unsigned int i=0;i<3;i++){
		CHECK_CLOSE(grad[i],x1g[i]);
	}

	double eval = 30.394751221800913;

	CHECK_CLOSE(val1,eval);
	CHECK_CLOSE(val2,eval);


	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,0);
}

BOOST_AUTO_TEST_CASE( test_grad_sin )
{
	BOOST_TEST_MESSAGE("test_grad_sin");
	VNode* x1 = create_var_node();
	x1->val = 10;
	OPNode* root = create_uary_op_node(OP_SIN,x1);
	vector<Node*> nodes;
	nodes.push_back(x1);
	vector<double> grad;
	grad_reverse(root,nodes,grad);
	double x1g =  -0.83907152907645244;
	//the matlab give cos(10) = -0.839071529076452

	CHECK_CLOSE(grad[0],x1g);
	BOOST_CHECK_EQUAL(nodes.size(),1);

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,1);
}

BOOST_AUTO_TEST_CASE(test_grad_single_node)
{
	VNode* x1 = create_var_node();
	x1->val = -2;
	vector<Node*> nodes;
	nodes.push_back(x1);
	vector<double> grad;
	double val = grad_reverse(x1,nodes,grad);
	CHECK_CLOSE(grad[0],1);
	CHECK_CLOSE(val,-2);

	EdgeSet s;
	unsigned int n = 0;
	nonlinearEdges(x1,s);
	n = nzHess(s);
	BOOST_CHECK_EQUAL(n,0);

	grad.clear();
	nodes.clear();
	PNode* p = create_param_node(-10);
	//OPNode*	op = create_binary_op_node(TIMES,p,create_param_node(2));
	val = grad_reverse(p,nodes,grad);
	BOOST_CHECK_EQUAL(grad.size(),0);
	CHECK_CLOSE(val,-10);

	s.clear();
	nonlinearEdges(p,s);
	n = nzHess(s);
	BOOST_CHECK_EQUAL(n,0);
}

BOOST_AUTO_TEST_CASE(test_grad_neg)
{
	VNode* x1 = create_var_node();
	x1->val = 10;
	PNode* p2 = create_param_node(-1);
	vector<Node*> nodes;
	vector<double> grad;
	nodes.push_back(x1);
	Node* root = create_binary_op_node(OP_TIMES,x1,p2);
	grad_reverse(root,nodes,grad);
	CHECK_CLOSE(grad[0],-1);
	BOOST_CHECK_EQUAL(nodes.size(),1);
	nodes.clear();
	grad.clear();
	nodes.push_back(x1);
	root = create_uary_op_node(OP_NEG,x1);
	grad_reverse(root,nodes,grad);
	CHECK_CLOSE(grad[0],-1);

	EdgeSet s;
	unsigned int n = 0;
	nonlinearEdges(root,s);
	n = nzHess(s);
	BOOST_CHECK_EQUAL(n,0);
}

BOOST_AUTO_TEST_CASE( test_nl_function)
{
	vector<Node*> list;
	Node* root = build_nl_function1(list);
	double val = eval_function(root);
	vector<double> grad;
	grad_reverse(root,list,grad);
	double eval =-66.929555552886214;
	double gx[] = {-4.961306690356109,-9.444611307649055,-2.064383410399700,7.123100000000000};
	CHECK_CLOSE(val,eval);

	for(unsigned int i=0;i<4;i++)
	{
		CHECK_CLOSE(grad[i],gx[i]);
	}
	unsigned int nzgrad = nzGrad(root);
	unsigned int tol = numTotalNodes(root);
	BOOST_CHECK_EQUAL(nzgrad,4);
	BOOST_CHECK_EQUAL(tol,16);

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,11);
}

BOOST_AUTO_TEST_CASE( test_hess_reverse_1)
{
	vector<Node*> nodes;
	Node* root = build_linear_fun1(nodes);
	vector<double> grad;
	double val = grad_reverse(root,nodes,grad);
	double eval = eval_function(root);
//	cout<<eval<<"\t"<<grad[0]<<"\t"<<grad[1]<<"\t"<<grad[2]<<"\t"<<endl;

	CHECK_CLOSE(val,eval);

	for(unsigned int i=0;i<nodes.size();i++)
	{
		static_cast<VNode*>(nodes[i])->u = 0;
	}

	static_cast<VNode*>(nodes[0])->u = 1;
	double hval = 0;
	vector<double> dhess;
	hval = hess_reverse(root,nodes,dhess);
	CHECK_CLOSE(hval,eval);
	for(unsigned int i=0;i<dhess.size();i++)
	{
		CHECK_CLOSE(dhess[i],0);
	}
}

BOOST_AUTO_TEST_CASE( test_hess_reverse_2)
{
	vector<Node*> nodes;
	Node* root = build_linear_function2(nodes);
	vector<double> grad;
	double val = grad_reverse(root,nodes,grad);
	double eval = eval_function(root);

	CHECK_CLOSE(val,eval);
	for(unsigned int i=0;i<nodes.size();i++)
	{
		static_cast<VNode*>(nodes[i])->u = 0;
	}

	static_cast<VNode*>(nodes[0])->u = 1;
	double hval = 0;
	vector<double> dhess;
	hval = hess_reverse(root,nodes,dhess);
	CHECK_CLOSE(hval,eval);

	for(unsigned int i=0;i<dhess.size();i++)
	{
		CHECK_CLOSE(dhess[i],0);
	}

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,0);
}

BOOST_AUTO_TEST_CASE( test_hess_reverse_4)
{
	vector<Node*> nodes;
//	Node* root = build_nl_function1(nodes);

	VNode* x1 = create_var_node();
	nodes.push_back(x1);
	x1->val = 1;
	x1->u =1;
	Node* op = create_uary_op_node(OP_SIN,x1);
	Node* root = create_uary_op_node(OP_SIN,op);
	vector<double> grad;
	double eval = eval_function(root);
	vector<double> dhess;
	double hval = hess_reverse(root,nodes,dhess);
	CHECK_CLOSE(hval,eval);
	BOOST_CHECK_EQUAL(dhess.size(),1);
	CHECK_CLOSE(dhess[0], -0.778395788418109);

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,1);
}

BOOST_AUTO_TEST_CASE( test_hess_reverse_3)
{
	vector<Node*> nodes;
	VNode* x1 = create_var_node();
	VNode* x2 = create_var_node();
	nodes.push_back(x1);
	nodes.push_back(x2);
	x1->val = 2.5;
	x2->val = -9;
	Node* op1 = create_binary_op_node(OP_TIMES,x1,x2);
	Node* root = create_binary_op_node(OP_TIMES,x1,op1);
	double eval = eval_function(root);
	for(unsigned int i=0;i<nodes.size();i++)
	{
		static_cast<VNode*>(nodes[i])->u = 0;
	}
	static_cast<VNode*>(nodes[0])->u = 1;

	vector<double> dhess;
	double hval = hess_reverse(root,nodes,dhess);
	BOOST_CHECK_EQUAL(dhess.size(),2);
	CHECK_CLOSE(hval,eval);
	double hx[]={-18,5};
	for(unsigned int i=0;i<dhess.size();i++)
	{
		//Print("\t["<<i<<"]="<<dhess[i]);
		CHECK_CLOSE(dhess[i],hx[i]);
	}

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,3);
}

BOOST_AUTO_TEST_CASE( test_hess_reverse_5)
{
	vector<Node*> nodes;
	VNode* x1 = create_var_node();
	VNode* x2 = create_var_node();
	nodes.push_back(x1);
	nodes.push_back(x2);
	x1->val = 2.5;
	x2->val = -9;
	Node* op1 = create_binary_op_node(OP_TIMES,x1,x1);
	Node* op2 = create_binary_op_node(OP_TIMES,x2,x2);
	Node* op3 = create_binary_op_node(OP_MINUS,op1,op2);
	Node* op4 = create_binary_op_node(OP_PLUS,op1,op2);
	Node* root = create_binary_op_node(OP_TIMES,op3,op4);

	double eval = eval_function(root);

	for(unsigned int i=0;i<nodes.size();i++)
	{
		static_cast<VNode*>(nodes[i])->u = 0;
	}
	static_cast<VNode*>(nodes[0])->u = 1;

	vector<double> dhess;
	double hval = hess_reverse(root,nodes,dhess);
	CHECK_CLOSE(hval,eval);
	double hx[] ={75,0};
	for(unsigned int i=0;i<dhess.size();i++)
	{
		CHECK_CLOSE(dhess[i],hx[i]);
	}

	for(unsigned int i=0;i<nodes.size();i++)
	{
		static_cast<VNode*>(nodes[i])->u = 0;
	}
	static_cast<VNode*>(nodes[1])->u = 1;

	double hx2[] = {0, -972};
	hval = hess_reverse(root,nodes,dhess);
	for(unsigned int i=0;i<dhess.size();i++)
	{
		CHECK_CLOSE(dhess[i],hx2[i]);
	}

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,4);
}
BOOST_AUTO_TEST_CASE( test_hess_reverse_6)
{
	vector<Node*> nodes;
//	Node* root = build_nl_function1(nodes);

	VNode* x1 = create_var_node();
	VNode* x2 = create_var_node();
	nodes.push_back(x1);
	nodes.push_back(x2);
	x1->val = 2.5;
	x2->val = -9;
	Node* root = create_binary_op_node(OP_POW,x1,x2);

	double eval = eval_function(root);

	static_cast<VNode*>(nodes[0])->u=1;static_cast<VNode*>(nodes[1])->u=0;
	vector<double> dhess;
	double hval = hess_reverse(root,nodes,dhess);
	CHECK_CLOSE(hval,eval);
	double hx1[] ={0.003774873600000 , -0.000759862823419};
	double hx2[] ={-0.000759862823419, 0.000220093141567};
	for(unsigned int i=0;i<dhess.size();i++)
	{
		CHECK_CLOSE(dhess[i],hx1[i]);
	}
	static_cast<VNode*>(nodes[0])->u=0;static_cast<VNode*>(nodes[1])->u=1;
	hess_reverse(root,nodes,dhess);
	for(unsigned int i=0;i<dhess.size();i++)
	{
		CHECK_CLOSE(dhess[i],hx2[i]);
	}

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,4);
}

BOOST_AUTO_TEST_CASE( test_hess_reverse_7)
{
	vector<Node*> nodes;
	Node* root = build_nl_function1(nodes);

	double eval = eval_function(root);


	vector<double> dhess;
	double hx0[] ={-1.747958066718855,
			  -0.657091724418110,
			   2.410459188139686,
			                   0};
	double hx1[] ={ -0.657091724418110,
			   0.006806564792590,
			  -0.289815306593997,
			   1.000000000000000};
	double hx2[] ={  2.410459188139686,
			  -0.289815306593997,
			   2.064383410399700,
			                   0};
	double hx3[] ={0,1,0,0};
	for(unsigned int i=0;i<nodes.size();i++)
	{
		static_cast<VNode*>(nodes[i])->u = 0;
	}
	static_cast<VNode*>(nodes[0])->u = 1;
	double hval = hess_reverse(root,nodes,dhess);
	CHECK_CLOSE(hval,eval);
	for(unsigned int i=0;i<dhess.size();i++)
	{
		CHECK_CLOSE(dhess[i],hx0[i]);
	}

	for (unsigned int i = 0; i < nodes.size(); i++) {
		static_cast<VNode*>(nodes[i])->u = 0;
	}
	static_cast<VNode*>(nodes[1])->u = 1;
	hess_reverse(root, nodes, dhess);
	for (unsigned int i = 0; i < dhess.size(); i++) {
		CHECK_CLOSE(dhess[i], hx1[i]);
	}

	for (unsigned int i = 0; i < nodes.size(); i++) {
		static_cast<VNode*>(nodes[i])->u = 0;
	}
	static_cast<VNode*>(nodes[2])->u = 1;
	hess_reverse(root, nodes, dhess);
	for (unsigned int i = 0; i < dhess.size(); i++) {
		CHECK_CLOSE(dhess[i], hx2[i]);
	}

	for (unsigned int i = 0; i < nodes.size(); i++) {
		static_cast<VNode*>(nodes[i])->u = 0;
	}
	static_cast<VNode*>(nodes[3])->u = 1;
	hess_reverse(root, nodes, dhess);
	for (unsigned i = 0; i < dhess.size(); i++) {
		CHECK_CLOSE(dhess[i], hx3[i]);
	}
}

#if FORWARD_ENABLED
void test_hess_forward(Node* root, unsigned int& nvar)
{
	AutoDiff::num_var = nvar;
	unsigned int len = (nvar+3)*nvar/2;
	double* hess = new double[len];
	hess_forward(root,nvar,&hess);
	for(unsigned int i=0;i<len;i++){
		cout<<"hess["<<i<<"]="<<hess[i]<<endl;
	}
	delete[] hess;
}
#endif

BOOST_AUTO_TEST_CASE( test_hess_reverse_8)
{
	vector<Node*> list;
	vector<double> dhess;

	VNode* x1 = create_var_node();
	list.push_back(x1);
	static_cast<VNode*>(list[0])->val = -10.5;
	static_cast<VNode*>(list[0])->u = 1;
	double deval = hess_reverse(x1,list,dhess);
	CHECK_CLOSE(deval,-10.5);
	BOOST_CHECK_EQUAL(dhess.size(),1);
	BOOST_CHECK(isnan(dhess[0]));

	EdgeSet s;
	nonlinearEdges(x1,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,0);

	PNode* p1 = create_param_node(-1.5);
	list.clear();
	deval = hess_reverse(p1,list,dhess);
	CHECK_CLOSE(deval,-1.5);
	BOOST_CHECK_EQUAL(dhess.size(),0);

	s.clear();
	nonlinearEdges(p1,s);
	n = nzHess(s);
	BOOST_CHECK_EQUAL(n,0);
}

BOOST_AUTO_TEST_CASE( test_hess_revers9)
{
	vector<Node*> list;
	vector<double> dhess;
	VNode* x1 = create_var_node();
	list.push_back(x1);
	static_cast<VNode*>(list[0])->val = 2.5;
	static_cast<VNode*>(list[0])->u =1;
	Node* op1 = create_binary_op_node(OP_TIMES,x1,x1);
	Node* root = create_binary_op_node(OP_TIMES,op1,op1);
	double deval = hess_reverse(root,list,dhess);
	double eval = eval_function(root);
	CHECK_CLOSE(eval,deval);
	BOOST_CHECK_EQUAL(dhess.size(),1);
	CHECK_CLOSE(dhess[0],75);

	EdgeSet s;
	nonlinearEdges(root,s);
	unsigned int n = nzHess(s);
	BOOST_CHECK_EQUAL(n,1);

}

BOOST_AUTO_TEST_CASE( test_hess_revers10)
{
	vector<Node*> list;
	vector<double> dhess;
	VNode* x1 = create_var_node();
	VNode* x2 = create_var_node();
	list.push_back(x1);
	list.push_back(x2);
	Node* op1 = create_binary_op_node(OP_TIMES, x1,x2);
	Node* op2 = create_uary_op_node(OP_SIN,op1);
	Node* op3 = create_uary_op_node(OP_COS,op1);
	Node* root = create_binary_op_node(OP_TIMES, op2, op3);
	static_cast<VNode*>(list[0])->val = 2.1;
	static_cast<VNode*>(list[1])->val = 1.8;
	double eval = eval_function(root);

	//second column
	static_cast<VNode*>(list[0])->u = 0;
	static_cast<VNode*>(list[1])->u = 1;
	double deval = hess_reverse(root,list,dhess);
	CHECK_CLOSE(eval,deval);
	BOOST_CHECK_EQUAL(dhess.size(),2);
	CHECK_CLOSE(dhess[0], -6.945893481707861);
	CHECK_CLOSE(dhess[1],  -8.441601940854081);

	//first column
	static_cast<VNode*>(list[0])->u = 1;
	static_cast<VNode*>(list[1])->u = 0;
	deval = hess_reverse(root,list,dhess);
	CHECK_CLOSE(eval,deval);
	BOOST_CHECK_EQUAL(dhess.size(),2);
	CHECK_CLOSE(dhess[0], -6.201993262668304);
	CHECK_CLOSE(dhess[1], -6.945893481707861);
}

BOOST_AUTO_TEST_CASE( test_grad_reverse11)
{
	vector<Node*> list;
	VNode* x1 = create_var_node();
	Node* p2 = create_param_node(2);
	list.push_back(x1);
	Node* op1 = create_binary_op_node(OP_POW,x1,p2);
	static_cast<VNode*>(x1)->val = 0;
	vector<double> grad;
	grad_reverse(op1,list,grad);
	BOOST_CHECK_EQUAL(grad.size(),1);
	CHECK_CLOSE(grad[0],0);
}

BOOST_AUTO_TEST_CASE( test_hess_reverse12)
{
	vector<Node*> list;
	VNode* x1 = create_var_node();
	Node* p2 = create_param_node(2);
	list.push_back(x1);
	Node* op1 = create_binary_op_node(OP_POW,x1,p2);
	x1->val = 0;
	x1->u = 1;
	vector<double> hess;
	hess_reverse(op1,list,hess);
	BOOST_CHECK_EQUAL(hess.size(),1);
	CHECK_CLOSE(hess[0],2);
}

BOOST_AUTO_TEST_CASE( test_grad_reverse13)
{
	vector<Node*> list;
	VNode* x1 = create_var_node();
	PNode* p1 = create_param_node(0.090901);
	VNode* x2 = create_var_node();
	PNode* p2 = create_param_node(0.090901);
	list.push_back(x1);
	list.push_back(x2);
	Node* op1 = create_binary_op_node(OP_TIMES,x1,p1);
	Node* op2 = create_binary_op_node(OP_TIMES,x2,p2);
	Node* root = create_binary_op_node(OP_PLUS,op1,op2);
	x1->val = 1;
	x2->val = 1;
	vector<double> grad;
	grad_reverse(root,list,grad);
	BOOST_CHECK_EQUAL(grad.size(),2);
	CHECK_CLOSE(grad[0],0.090901);
	CHECK_CLOSE(grad[1],0.090901);
}

BOOST_AUTO_TEST_SUITE_END()