Adding upstream version 14.2.21.upstream/14.2.21upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

27 files changed, 2676 insertions, 0 deletions

diff --git a/src/boost/libs/yap/example/autodiff_library/ActNode.cpp b/src/boost/libs/yap/example/autodiff_library/ActNode.cpp
new file mode 100644
index 00000000..c21594f8
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/ActNode.cpp
@@ -0,0 +1,34 @@
+/*
+ * ActNode.cpp
+ *
+ *  Created on: 13 Apr 2013
+ *      Author: s0965328
+ */
+
+#include "ActNode.h"
+
+namespace AutoDiff {
+
+ActNode::ActNode() : AutoDiff::Node(),adj(NaN_Double){
+
+}
+
+ActNode::~ActNode() {
+}
+
+
+void ActNode::update_adj(double& v)
+{
+	assert(!isnan(adj));
+	assert(!isnan(v));
+	adj+=v;
+}
+
+void ActNode::grad_reverse_1_init_adj()
+{
+	adj = 1;
+}
+
+}
+
+
diff --git a/src/boost/libs/yap/example/autodiff_library/ActNode.h b/src/boost/libs/yap/example/autodiff_library/ActNode.h
new file mode 100644
index 00000000..58e77cc3
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/ActNode.h
@@ -0,0 +1,29 @@
+/*
+ * ActNode.h
+ *
+ *  Created on: 13 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef ACTNODE_H_
+#define ACTNODE_H_
+
+#include "Node.h"
+
+
+namespace AutoDiff {
+
+class ActNode : public Node{
+public:
+	ActNode();
+	virtual ~ActNode();
+
+	void update_adj(double& v);
+	void grad_reverse_1_init_adj();
+
+	double adj;
+};
+
+} // end namespace foo
+
+#endif /* ACTNODE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/BinaryOPNode.cpp b/src/boost/libs/yap/example/autodiff_library/BinaryOPNode.cpp
new file mode 100644
index 00000000..992d692d
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/BinaryOPNode.cpp
@@ -0,0 +1,651 @@
+/*
+ * BinaryOPNode.cpp
+ *
+ *  Created on: 6 Nov 2013
+ *      Author: s0965328
+ */
+
+#include "auto_diff_types.h"
+#include "BinaryOPNode.h"
+#include "PNode.h"
+#include "Stack.h"
+#include "Tape.h"
+#include "EdgeSet.h"
+#include "Node.h"
+#include "VNode.h"
+#include "OPNode.h"
+#include "ActNode.h"
+#include "EdgeSet.h"
+
+namespace AutoDiff {
+
+BinaryOPNode::BinaryOPNode(OPCODE op_, Node* left_, Node* right_):OPNode(op_,left_),right(right_)
+{
+}
+
+OPNode* BinaryOPNode::createBinaryOpNode(OPCODE op, Node* left, Node* right)
+{
+	assert(left!=NULL && right!=NULL);
+	OPNode* node = NULL;
+	node = new BinaryOPNode(op,left,right);
+	return node;
+}
+
+BinaryOPNode::~BinaryOPNode() {
+	if(right->getType()!=VNode_Type)
+	{
+		delete right;
+		right = NULL;
+	}
+}
+
+void BinaryOPNode::inorder_visit(int level,ostream& oss){
+	if(left!=NULL){
+		left->inorder_visit(level+1,oss);
+	}
+	oss<<this->toString(level)<<endl;
+	if(right!=NULL){
+		right->inorder_visit(level+1,oss);
+	}
+}
+
+void BinaryOPNode::collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total){
+	total++;
+	if (left != NULL) {
+		left->collect_vnodes(nodes,total);
+	}
+	if (right != NULL) {
+		right->collect_vnodes(nodes,total);
+	}
+
+}
+
+void BinaryOPNode::eval_function()
+{
+	assert(left!=NULL && right!=NULL);
+	left->eval_function();
+	right->eval_function();
+	this->calc_eval_function();
+}
+
+void BinaryOPNode::calc_eval_function()
+{
+	double x = NaN_Double;
+	double rx = SV->pop_back();
+	double lx = SV->pop_back();
+	switch (op)
+	{
+	case OP_PLUS:
+		x = lx + rx;
+		break;
+	case OP_MINUS:
+		x = lx - rx;
+		break;
+	case OP_TIMES:
+		x = lx * rx;
+		break;
+	case OP_DIVID:
+		x = lx / rx;
+		break;
+	case OP_POW:
+		x = pow(lx,rx);
+		break;
+	default:
+		cerr<<"op["<<op<<"] not yet implemented!!"<<endl;
+		assert(false);
+		break;
+	}
+	SV->push_back(x);
+}
+
+
+//1. visiting left if not NULL
+//2. then, visiting right if not NULL
+//3. calculating the immediate derivative hu and hv
+void BinaryOPNode::grad_reverse_0()
+{
+	assert(left!=NULL && right != NULL);
+	this->adj = 0;
+	left->grad_reverse_0();
+	right->grad_reverse_0();
+	this->calc_grad_reverse_0();
+}
+
+//right left - right most traversal
+void BinaryOPNode::grad_reverse_1()
+{
+	assert(right!=NULL && left!=NULL);
+	double r_adj = SD->pop_back()*this->adj;
+	right->update_adj(r_adj);
+	double l_adj = SD->pop_back()*this->adj;
+	left->update_adj(l_adj);
+
+	right->grad_reverse_1();
+	left->grad_reverse_1();
+}
+
+void BinaryOPNode::calc_grad_reverse_0()
+{
+	assert(left!=NULL && right != NULL);
+	double l_dh = NaN_Double;
+	double r_dh = NaN_Double;
+	double rx = SV->pop_back();
+	double lx = SV->pop_back();
+	double x = NaN_Double;
+	switch (op)
+	{
+	case OP_PLUS:
+		x = lx + rx;
+		l_dh = 1;
+		r_dh = 1;
+		break;
+	case OP_MINUS:
+		x = lx - rx;
+		l_dh = 1;
+		r_dh = -1;
+		break;
+	case OP_TIMES:
+		x = lx * rx;
+		l_dh = rx;
+		r_dh = lx;
+		break;
+	case OP_DIVID:
+		x = lx / rx;
+		l_dh = 1 / rx;
+		r_dh = -(lx) / pow(rx, 2);
+		break;
+	case OP_POW:
+		if(right->getType()==PNode_Type){
+			x = pow(lx,rx);
+			l_dh = rx*pow(lx,(rx-1));
+			r_dh = 0;
+		}
+		else{
+			assert(lx>0.0); //otherwise log(lx) is not defined in read number
+			x = pow(lx,rx);
+			l_dh = rx*pow(lx,(rx-1));
+			r_dh = pow(lx,rx)*log(lx);  //this is for x1^x2 when x1=0 cause r_dh become +inf, however d(0^x2)/d(x2) = 0
+		}
+		break;
+	default:
+		cerr<<"error op not impl"<<endl;
+		break;
+	}
+	SV->push_back(x);
+	SD->push_back(l_dh);
+	SD->push_back(r_dh);
+}
+
+void BinaryOPNode::hess_reverse_0_init_n_in_arcs()
+{
+	this->left->hess_reverse_0_init_n_in_arcs();
+	this->right->hess_reverse_0_init_n_in_arcs();
+	this->Node::hess_reverse_0_init_n_in_arcs();
+}
+
+void BinaryOPNode::hess_reverse_1_clear_index()
+{
+	this->left->hess_reverse_1_clear_index();
+	this->right->hess_reverse_1_clear_index();
+	this->Node::hess_reverse_1_clear_index();
+}
+
+unsigned int BinaryOPNode::hess_reverse_0()
+{
+	assert(this->left!=NULL && right!=NULL);
+	if(index==0)
+	{
+		unsigned int lindex=0, rindex=0;
+		lindex = left->hess_reverse_0();
+		rindex = right->hess_reverse_0();
+		assert(lindex!=0 && rindex !=0);
+		II->set(lindex);
+		II->set(rindex);
+		double rx,rx_bar,rw,rw_bar;
+		double lx,lx_bar,lw,lw_bar;
+		double x,x_bar,w,w_bar;
+		double r_dh, l_dh;
+		right->hess_reverse_0_get_values(rindex,rx,rx_bar,rw,rw_bar);
+		left->hess_reverse_0_get_values(lindex,lx,lx_bar,lw,lw_bar);
+		switch(op)
+		{
+		case OP_PLUS:
+	//		cout<<"lindex="<<lindex<<"\trindex="<<rindex<<"\tI="<<I<<endl;
+			x = lx + rx;
+	//		cout<<lx<<"\t+"<<rx<<"\t="<<x<<"\t\t"<<toString(0)<<endl;
+			x_bar = 0;
+			l_dh = 1;
+			r_dh = 1;
+			w = lw * l_dh  + rw * r_dh;
+	//		cout<<lw<<"\t+"<<rw<<"\t="<<w<<"\t\t"<<toString(0)<<endl;
+			w_bar = 0;
+			break;
+		case OP_MINUS:
+			x = lx - rx;
+			x_bar = 0;
+			l_dh = 1;
+			r_dh = -1;
+			w = lw * l_dh  + rw * r_dh;
+			w_bar = 0;
+			break;
+		case OP_TIMES:
+			x = lx * rx;
+			x_bar = 0;
+			l_dh = rx;
+			r_dh = lx;
+			w = lw * l_dh  + rw * r_dh;
+			w_bar = 0;
+			break;
+		case OP_DIVID:
+			x = lx / rx;
+			x_bar = 0;
+			l_dh = 1/rx;
+			r_dh = -lx/pow(rx,2);
+			w = lw * l_dh  + rw * r_dh;
+			w_bar = 0;
+			break;
+		case OP_POW:
+			if(right->getType()==PNode_Type)
+			{
+				x = pow(lx,rx);
+				x_bar = 0;
+				l_dh = rx*pow(lx,(rx-1));
+				r_dh = 0;
+				w = lw * l_dh + rw * r_dh;
+				w_bar = 0;
+			}
+			else
+			{
+				assert(lx>0.0); //otherwise log(lx) undefined in real number
+				x = pow(lx,rx);
+				x_bar = 0;
+				l_dh = rx*pow(lx,(rx-1));
+				r_dh = pow(lx,rx)*log(lx);   //log(lx) cause -inf when lx=0;
+				w = lw * l_dh  + rw * r_dh;
+				w_bar = 0;
+			}
+			break;
+		default:
+			cerr<<"op["<<op<<"] not yet implemented!"<<endl;
+			assert(false);
+			break;
+		}
+		TT->set(x);
+		TT->set(x_bar);
+		TT->set(w);
+		TT->set(w_bar);
+		TT->set(l_dh);
+		TT->set(r_dh);
+		assert(TT->index == TT->index);
+		index = TT->index;
+	}
+	return index;
+}
+
+void BinaryOPNode::hess_reverse_0_get_values(unsigned int i,double& x, double& x_bar, double& w, double& w_bar)
+{
+	--i; // skip the r_dh (ie, dh/du)
+	--i; // skip the l_dh (ie. dh/dv)
+	w_bar = TT->get(--i);
+	w = TT->get(--i);
+	x_bar = TT->get(--i);
+	x = TT->get(--i);
+}
+
+void BinaryOPNode::hess_reverse_1(unsigned int i)
+{
+	n_in_arcs--;
+	if(n_in_arcs==0)
+	{
+		assert(right!=NULL && left!=NULL);
+		unsigned int rindex = II->get(--(II->index));
+		unsigned int lindex = II->get(--(II->index));
+	//	cout<<"ri["<<rindex<<"]\tli["<<lindex<<"]\t"<<this->toString(0)<<endl;
+		double r_dh = TT->get(--i);
+		double l_dh = TT->get(--i);
+		double w_bar = TT->get(--i);
+		--i; //skip w
+		double x_bar = TT->get(--i);
+		--i; //skip x
+
+		double lw_bar=0,rw_bar=0;
+		double lw=0,lx=0; left->hess_reverse_1_get_xw(lindex,lw,lx);
+		double rw=0,rx=0; right->hess_reverse_1_get_xw(rindex,rw,rx);
+		switch(op)
+		{
+		case OP_PLUS:
+			assert(l_dh==1);
+			assert(r_dh==1);
+			lw_bar += w_bar*l_dh;
+			rw_bar += w_bar*r_dh;
+			break;
+		case OP_MINUS:
+			assert(l_dh==1);
+			assert(r_dh==-1);
+			lw_bar += w_bar*l_dh;
+			rw_bar += w_bar*r_dh;
+			break;
+		case OP_TIMES:
+			assert(rx == l_dh);
+			assert(lx == r_dh);
+			lw_bar += w_bar*rx;
+			lw_bar += x_bar*lw*0 + x_bar*rw*1;
+			rw_bar += w_bar*lx;
+			rw_bar += x_bar*lw*1 + x_bar*rw*0;
+			break;
+		case OP_DIVID:
+			lw_bar += w_bar*l_dh;
+			lw_bar += x_bar*lw*0 + x_bar*rw*-1/(pow(rx,2));
+			rw_bar += w_bar*r_dh;
+			rw_bar += x_bar*lw*-1/pow(rx,2) + x_bar*rw*2*lx/pow(rx,3);
+			break;
+		case OP_POW:
+			if(right->getType()==PNode_Type){
+				lw_bar += w_bar*l_dh;
+				lw_bar += x_bar*lw*pow(lx,rx-2)*rx*(rx-1) + 0;
+				rw_bar += w_bar*r_dh; assert(r_dh==0.0);
+				rw_bar += 0;
+			}
+			else{
+				assert(lx>0.0); //otherwise log(lx) is not define in Real
+				lw_bar += w_bar*l_dh;
+				lw_bar += x_bar*lw*pow(lx,rx-2)*rx*(rx-1) + x_bar*rw*pow(lx,rx-1)*(rx*log(lx)+1);  //cause log(lx)=-inf when
+				rw_bar += w_bar*r_dh;
+				rw_bar += x_bar*lw*pow(lx,rx-1)*(rx*log(lx)+1) + x_bar*rw*pow(lx,rx)*pow(log(lx),2);
+			}
+			break;
+		default:
+			cerr<<"op["<<op<<"] not yet implemented !"<<endl;
+			assert(false);
+			break;
+		}
+		double rx_bar = x_bar*r_dh;
+		double lx_bar = x_bar*l_dh;
+		right->update_x_bar(rindex,rx_bar);
+		left->update_x_bar(lindex,lx_bar);
+		right->update_w_bar(rindex,rw_bar);
+		left->update_w_bar(lindex,lw_bar);
+
+
+		this->right->hess_reverse_1(rindex);
+		this->left->hess_reverse_1(lindex);
+	}
+}
+void BinaryOPNode::hess_reverse_1_init_x_bar(unsigned int i)
+{
+	TT->at(i-5) = 1;
+}
+void BinaryOPNode::update_x_bar(unsigned int i ,double v)
+{
+	TT->at(i-5) += v;
+}
+void BinaryOPNode::update_w_bar(unsigned int i ,double v)
+{
+	TT->at(i-3) += v;
+}
+void BinaryOPNode::hess_reverse_1_get_xw(unsigned int i,double& w,double& x)
+{
+	w = TT->get(i-4);
+	x = TT->get(i-6);
+}
+void BinaryOPNode::hess_reverse_get_x(unsigned int i,double& x)
+{
+	x = TT->get(i-6);
+}
+
+
+void BinaryOPNode::nonlinearEdges(EdgeSet& edges)
+{
+	for(list<Edge>::iterator it=edges.edges.begin();it!=edges.edges.end();)
+	{
+		Edge e = *it;
+		if(e.a==this || e.b == this){
+			if(e.a == this && e.b == this)
+			{
+				Edge e1(left,left);
+				Edge e2(right,right);
+				Edge e3(left,right);
+				edges.insertEdge(e1);
+				edges.insertEdge(e2);
+				edges.insertEdge(e3);
+			}
+			else
+			{
+				Node* o = e.a==this? e.b: e.a;
+				Edge e1(left,o);
+				Edge e2(right,o);
+				edges.insertEdge(e1);
+				edges.insertEdge(e2);
+			}
+			it = edges.edges.erase(it);
+		}
+		else
+		{
+			it++;
+		}
+	}
+
+	Edge e1(left,right);
+	Edge e2(left,left);
+	Edge e3(right,right);
+	switch(op)
+	{
+	case OP_PLUS:
+	case OP_MINUS:
+		//do nothing for linear operator
+		break;
+	case OP_TIMES:
+		edges.insertEdge(e1);
+		break;
+	case OP_DIVID:
+		edges.insertEdge(e1);
+		edges.insertEdge(e3);
+		break;
+	case OP_POW:
+		edges.insertEdge(e1);
+		edges.insertEdge(e2);
+		edges.insertEdge(e3);
+		break;
+	default:
+		cerr<<"op["<<op<<"] not yet implmented !"<<endl;
+		assert(false);
+		break;
+	}
+	left->nonlinearEdges(edges);
+	right->nonlinearEdges(edges);
+}
+
+#if FORWARD_ENABLED
+
+void BinaryOPNode::hess_forward(unsigned int len, double** ret_vec)
+{
+	double* lvec = NULL;
+	double* rvec = NULL;
+	if(left!=NULL){
+		left->hess_forward(len,&lvec);
+	}
+	if(right!=NULL){
+		right->hess_forward(len,&rvec);
+	}
+
+	*ret_vec = new double[len];
+	hess_forward_calc0(len,lvec,rvec,*ret_vec);
+	//delete lvec, rvec
+	delete[] lvec;
+	delete[] rvec;
+}
+
+void BinaryOPNode::hess_forward_calc0(unsigned int& len, double* lvec, double* rvec, double* ret_vec)
+{
+	double hu = NaN_Double, hv= NaN_Double;
+	double lval = NaN_Double, rval = NaN_Double;
+	double val = NaN_Double;
+	unsigned int index = 0;
+	switch (op)
+	{
+	case OP_PLUS:
+		rval = SV->pop_back();
+		lval = SV->pop_back();
+		val = lval + rval;
+		SV->push_back(val);
+		//calculate the first order derivatives
+		for(unsigned int i=0;i<AutoDiff::num_var;++i)
+		{
+			ret_vec[i] = lvec[i]+rvec[i];
+		}
+		//calculate the second order
+		index = AutoDiff::num_var;
+		for(unsigned int i=0;i<AutoDiff::num_var;++i)
+		{
+			for(unsigned int j=i;j<AutoDiff::num_var;++j){
+				ret_vec[index] = lvec[index] + 0 + rvec[index] + 0;
+				++index;
+			}
+		}
+		assert(index==len);
+		break;
+	case OP_MINUS:
+		rval = SV->pop_back();
+		lval = SV->pop_back();
+		val = lval + rval;
+		SV->push_back(val);
+		//calculate the first order derivatives
+		for(unsigned int i=0;i<AutoDiff::num_var;++i)
+		{
+			ret_vec[i] = lvec[i] - rvec[i];
+		}
+		//calculate the second order
+		index = AutoDiff::num_var;
+		for(unsigned int i=0;i<AutoDiff::num_var;++i)
+		{
+			for(unsigned int j=i;j<AutoDiff::num_var;++j){
+				ret_vec[index] = lvec[index] + 0 - rvec[index] + 0;
+				++index;
+			}
+		}
+		assert(index==len);
+		break;
+	case OP_TIMES:
+		rval = SV->pop_back();
+		lval = SV->pop_back();
+		val = lval * rval;
+		SV->push_back(val);
+		hu = rval;
+		hv = lval;
+		//calculate the first order derivatives
+		for(unsigned int i =0;i<AutoDiff::num_var;++i)
+		{
+			ret_vec[i] = hu*lvec[i] + hv*rvec[i];
+		}
+		//calculate the second order
+		index = AutoDiff::num_var;
+		for(unsigned int i=0;i<AutoDiff::num_var;++i)
+		{
+			for(unsigned int j=i;j<AutoDiff::num_var;++j)
+			{
+				ret_vec[index] = hu * lvec[index] + lvec[i] * rvec[j]+hv * rvec[index] + rvec[i] * lvec[j];
+				++index;
+			}
+		}
+		assert(index==len);
+		break;
+	case OP_POW:
+		rval = SV->pop_back();
+		lval = SV->pop_back();
+		val = pow(lval,rval);
+		SV->push_back(val);
+		if(left->getType()==PNode_Type && right->getType()==PNode_Type)
+		{
+			std::fill_n(ret_vec,len,0);
+		}
+		else
+		{
+			hu = rval*pow(lval,(rval-1));
+			hv = pow(lval,rval)*log(lval);
+			if(left->getType()==PNode_Type)
+			{
+				double coeff = pow(log(lval),2)*pow(lval,rval);
+				//calculate the first order derivatives
+				for(unsigned int i =0;i<AutoDiff::num_var;++i)
+				{
+					ret_vec[i] = hu*lvec[i] + hv*rvec[i];
+				}
+				//calculate the second order
+				index = AutoDiff::num_var;
+				for(unsigned int i=0;i<AutoDiff::num_var;++i)
+				{
+					for(unsigned int j=i;j<AutoDiff::num_var;++j)
+					{
+						ret_vec[index] = 0 + 0 + hv * rvec[index] + rvec[i] * coeff * rvec[j];
+						++index;
+					}
+				}
+			}
+			else if(right->getType()==PNode_Type)
+			{
+				double coeff = rval*(rval-1)*pow(lval,rval-2);
+				//calculate the first order derivatives
+				for(unsigned int i =0;i<AutoDiff::num_var;++i)
+				{
+					ret_vec[i] = hu*lvec[i] + hv*rvec[i];
+				}
+				//calculate the second order
+				index = AutoDiff::num_var;
+				for(unsigned int i=0;i<AutoDiff::num_var;++i)
+				{
+					for(unsigned int j=i;j<AutoDiff::num_var;++j)
+					{
+						ret_vec[index] = hu*lvec[index] + lvec[i] * coeff * lvec[j] + 0 + 0;
+						++index;
+					}
+				}
+			}
+			else
+			{
+				assert(false);
+			}
+		}
+		assert(index==len);
+		break;
+	case OP_SIN: //TODO should move to UnaryOPNode.cpp?
+		assert(left!=NULL&&right==NULL);
+		lval = SV->pop_back();
+		val = sin(lval);
+		SV->push_back(val);
+		hu = cos(lval);
+
+		double coeff;
+		coeff = -val; //=sin(left->val); -- and avoid cross initialisation
+		//calculate the first order derivatives
+		for(unsigned int i =0;i<AutoDiff::num_var;++i)
+		{
+			ret_vec[i] = hu*lvec[i] + 0;
+		}
+		//calculate the second order
+		index = AutoDiff::num_var;
+		for(unsigned int i=0;i<AutoDiff::num_var;++i)
+		{
+			for(unsigned int j=i;j<AutoDiff::num_var;++j)
+			{
+				ret_vec[index] = hu*lvec[index] + lvec[i] * coeff * lvec[j] + 0 + 0;
+				++index;
+			}
+		}
+		assert(index==len);
+		break;
+	default:
+		cerr<<"op["<<op<<"] not yet implemented!";
+		break;
+	}
+}
+#endif
+
+
+string BinaryOPNode::toString(int level){
+	ostringstream oss;
+	string s(level,'\t');
+	oss<<s<<"[BinaryOPNode]("<<op<<")";
+	return oss.str();
+}
+
+} /* namespace AutoDiff */
diff --git a/src/boost/libs/yap/example/autodiff_library/BinaryOPNode.h b/src/boost/libs/yap/example/autodiff_library/BinaryOPNode.h
new file mode 100644
index 00000000..7abc702a
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/BinaryOPNode.h
@@ -0,0 +1,57 @@
+/*
+ * BinaryOPNode.h
+ *
+ *  Created on: 6 Nov 2013
+ *      Author: s0965328
+ */
+
+#ifndef BINARYOPNODE_H_
+#define BINARYOPNODE_H_
+
+#include "OPNode.h"
+
+namespace AutoDiff {
+
+class EdgeSet;
+
+class BinaryOPNode: public OPNode {
+public:
+
+	static OPNode* createBinaryOpNode(OPCODE op, Node* left, Node* right);
+	virtual ~BinaryOPNode();
+
+	void collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total);
+	void eval_function();
+
+	void grad_reverse_0();
+	void grad_reverse_1();
+
+	void hess_forward(unsigned int len, double** ret_vec);
+
+	unsigned int hess_reverse_0();
+	void hess_reverse_0_init_n_in_arcs();
+	void hess_reverse_0_get_values(unsigned int,double&, double&, double&, double&);
+	void hess_reverse_1(unsigned int i);
+	void hess_reverse_1_init_x_bar(unsigned int);
+	void update_x_bar(unsigned int,double);
+	void update_w_bar(unsigned int,double);
+	void hess_reverse_1_get_xw(unsigned int, double&,double&);
+	void hess_reverse_get_x(unsigned int,double& x);
+	void hess_reverse_1_clear_index();
+
+	void nonlinearEdges(EdgeSet& a);
+
+	void inorder_visit(int level,ostream& oss);
+	string toString(int level);
+
+	Node* right;
+
+private:
+	BinaryOPNode(OPCODE op, Node* left, Node* right);
+	void calc_eval_function();
+	void calc_grad_reverse_0();
+	void hess_forward_calc0(unsigned int& len, double* lvec, double* rvec,double* ret_vec);
+};
+
+} /* namespace AutoDiff */
+#endif /* BINARYOPNODE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/Edge.cpp b/src/boost/libs/yap/example/autodiff_library/Edge.cpp
new file mode 100644
index 00000000..4a7183b6
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Edge.cpp
@@ -0,0 +1,54 @@
+/*
+ * Edge.cpp
+ *
+ *  Created on: 12 Nov 2013
+ *      Author: s0965328
+ */
+
+#include "Edge.h"
+#include <iostream>
+#include <sstream>
+
+namespace AutoDiff {
+
+Edge::Edge(Node* a_,Node* b_):a(a_),b(b_) {
+
+}
+
+Edge::~Edge() {
+	a = NULL;
+	b = NULL;
+}
+
+Edge::Edge(const Edge& e)
+{
+	a = e.a;
+	b = e.b;
+}
+
+bool Edge::isEqual(Edge* e)
+{
+	if(e->a == a && e->b == b)
+	{
+		return true;
+	}
+	else if(e->b == a && e->a == b)
+	{
+		return true;
+	}
+	return false;
+}
+
+bool Edge::isEqual(Edge& e)
+{
+	return isEqual(&e);
+}
+
+string Edge::toString()
+{
+	ostringstream oss;
+	oss<<""<<a->toString(0)<<"|"<<a<<" ----- "<<b->toString(0)<<"|"<<b<<""<<endl;
+	return oss.str();
+}
+
+} /* namespace AutoDiff */
diff --git a/src/boost/libs/yap/example/autodiff_library/Edge.h b/src/boost/libs/yap/example/autodiff_library/Edge.h
new file mode 100644
index 00000000..10126ac2
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Edge.h
@@ -0,0 +1,33 @@
+/*
+ * Edge.h
+ *
+ *  Created on: 12 Nov 2013
+ *      Author: s0965328
+ */
+
+#ifndef EDGE_H_
+#define EDGE_H_
+
+
+#include "Node.h"
+
+namespace AutoDiff {
+
+class Edge {
+public:
+	Edge(Node* a, Node* b);
+	Edge(const Edge& e);
+	virtual ~Edge();
+
+	bool isEqual(Edge*);
+	bool isEqual(Edge&);
+	std::string toString();
+
+	Node* a;
+	Node* b;
+
+
+};
+
+} /* namespace AutoDiff */
+#endif /* EDGE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/EdgeSet.cpp b/src/boost/libs/yap/example/autodiff_library/EdgeSet.cpp
new file mode 100644
index 00000000..20e16132
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/EdgeSet.cpp
@@ -0,0 +1,76 @@
+/*
+ * EdgeSet.cpp
+ *
+ *  Created on: 12 Nov 2013
+ *      Author: s0965328
+ */
+
+#include "EdgeSet.h"
+#include "Edge.h"
+#include <sstream>
+
+using namespace std;
+namespace AutoDiff {
+
+EdgeSet::EdgeSet() {
+	// TODO Auto-generated constructor stub
+
+}
+
+EdgeSet::~EdgeSet() {
+	edges.clear();
+}
+
+bool EdgeSet::containsEdge(Edge& e)
+{
+	list<Edge>::iterator it = edges.begin();
+	for(;it!= edges.end();it++){
+		Edge eit = *it;
+		if(eit.isEqual(e))
+		{
+			return true;
+		}
+	}
+	return false;
+}
+
+void EdgeSet::insertEdge(Edge& e) {
+	if(!containsEdge(e)){
+		edges.push_front(e);
+	}
+}
+
+void EdgeSet::clear() {
+	edges.clear();
+}
+
+unsigned int EdgeSet::size(){
+	return edges.size();
+}
+
+unsigned int EdgeSet::numSelfEdges(){
+	unsigned int diag = 0;
+	list<Edge>::iterator it = edges.begin();
+	for(;it!=edges.end();it++)
+	{
+		Edge eit = *it;
+		if(eit.a == eit.b)
+		{
+			diag++;
+		}
+	}
+	return diag;
+}
+
+string EdgeSet::toString()
+{
+	ostringstream oss;
+	list<Edge>::iterator it = edges.begin();
+	for(;it!=edges.end();it++)
+	{
+		oss<<(*it).toString()<<endl;
+	}
+	return oss.str();
+}
+
+} /* namespace AutoDiff */
diff --git a/src/boost/libs/yap/example/autodiff_library/EdgeSet.h b/src/boost/libs/yap/example/autodiff_library/EdgeSet.h
new file mode 100644
index 00000000..26b2be12
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/EdgeSet.h
@@ -0,0 +1,33 @@
+/*
+ * EdgeSet.h
+ *
+ *  Created on: 12 Nov 2013
+ *      Author: s0965328
+ */
+
+#ifndef EDGESET_H_
+#define EDGESET_H_
+
+#include "Edge.h"
+#include <list>
+
+namespace AutoDiff {
+
+class EdgeSet {
+public:
+	EdgeSet();
+	virtual ~EdgeSet();
+
+	void insertEdge(Edge& e);
+	bool containsEdge(Edge& e);
+	unsigned int numSelfEdges();
+	void clear();
+	unsigned int size();
+	std::string toString();
+
+
+	std::list<Edge> edges;
+};
+
+} /* namespace AutoDiff */
+#endif /* EDGESET_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/LICENSE b/src/boost/libs/yap/example/autodiff_library/LICENSE
new file mode 100644
index 00000000..f09534cf
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/LICENSE
@@ -0,0 +1,22 @@
+The MIT License (MIT)
+
+Copyright (c) 2014 fqiang
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
diff --git a/src/boost/libs/yap/example/autodiff_library/Node.cpp b/src/boost/libs/yap/example/autodiff_library/Node.cpp
new file mode 100644
index 00000000..13ed5764
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Node.cpp
@@ -0,0 +1,32 @@
+/*
+ * Node.cpp
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#include "Node.h"
+
+namespace AutoDiff {
+
+unsigned int Node::DEFAULT_INDEX = 0;
+Node::Node():index(Node::DEFAULT_INDEX),n_in_arcs(0){
+}
+
+
+Node::~Node() {
+}
+
+void Node::hess_reverse_0_init_n_in_arcs()
+{
+	n_in_arcs++;
+//	cout<<this->toString(0)<<endl;
+}
+
+void Node::hess_reverse_1_clear_index()
+{
+	index = Node::DEFAULT_INDEX;
+}
+
+}
+
diff --git a/src/boost/libs/yap/example/autodiff_library/Node.h b/src/boost/libs/yap/example/autodiff_library/Node.h
new file mode 100644
index 00000000..f89c2bb1
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Node.h
@@ -0,0 +1,65 @@
+/*
+ * Node.h
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef NODE_H_
+#define NODE_H_
+
+#include <boost/unordered_set.hpp>
+#include "auto_diff_types.h"
+
+using namespace std;
+
+namespace AutoDiff {
+
+class EdgeSet;
+
+class Node {
+public:
+	Node();
+	virtual ~Node();
+
+	virtual void eval_function() = 0;
+	virtual void grad_reverse_0() = 0;
+	virtual void grad_reverse_1_init_adj() = 0;
+	virtual void grad_reverse_1() = 0;
+	virtual void update_adj(double& v) = 0;
+	virtual unsigned int hess_reverse_0() = 0;
+	virtual void hess_reverse_0_init_n_in_arcs();
+	virtual void hess_reverse_0_get_values(unsigned int,double&, double&,double&, double&) = 0;
+	virtual void hess_reverse_1(unsigned int i) = 0;
+	virtual void hess_reverse_1_init_x_bar(unsigned int) = 0;
+	virtual void update_x_bar(unsigned int,double) = 0;
+	virtual void update_w_bar(unsigned int,double) = 0;
+	virtual void hess_reverse_1_get_xw(unsigned int, double&,double&) = 0;
+	virtual void hess_reverse_get_x(unsigned int,double& x)=0;
+	virtual void hess_reverse_1_clear_index();
+	//routing for checking non-zero structures
+	virtual void collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total) = 0;
+	virtual void nonlinearEdges(EdgeSet&) = 0;
+#if FORWARD_ENABLED
+	virtual void hess_forward(unsigned int len, double** ret_vec) = 0;
+#endif
+
+	//other utility methods
+	virtual void inorder_visit( int level,ostream& oss) = 0;
+	virtual string toString(int levl) = 0;
+	virtual TYPE getType() = 0;
+
+
+	//! index on the tape
+	unsigned int index;
+	//! number of incoming arcs
+	//! n_in_arcs in root node equals 1 before evaluation and 0 after evaluation
+	unsigned int n_in_arcs;
+
+	static unsigned int DEFAULT_INDEX;
+
+};
+
+} // end namespace foo
+
+#endif /* NODE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/OPNode.cpp b/src/boost/libs/yap/example/autodiff_library/OPNode.cpp
new file mode 100644
index 00000000..0f647931
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/OPNode.cpp
@@ -0,0 +1,37 @@
+/*
+ * OpNode.cpp
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#include "OPNode.h"
+#include "Stack.h"
+#include "Tape.h"
+#include "PNode.h"
+/***********************************************************
+	     h
+	    / \
+	   u   v  ----- hu  hv represent dh/du dh/dv resepectively.
+	   - - -
+	  x1....xn
+***********************************************************/
+
+namespace AutoDiff{
+
+OPNode::OPNode(OPCODE op, Node* left) : ActNode(), op(op), left(left),val(NaN_Double) {
+}
+
+TYPE OPNode::getType()
+{
+	return OPNode_Type;
+}
+
+OPNode::~OPNode() {
+	if(left->getType()!=VNode_Type)
+	{
+		delete left;
+		this->left = NULL;
+	}
+}
+}
diff --git a/src/boost/libs/yap/example/autodiff_library/OPNode.h b/src/boost/libs/yap/example/autodiff_library/OPNode.h
new file mode 100644
index 00000000..e5ed2da9
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/OPNode.h
@@ -0,0 +1,37 @@
+/*
+ * OpNode.h
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef OPNODE_H_
+#define OPNODE_H_
+
+#include "Node.h"
+#include "ActNode.h"
+
+namespace AutoDiff {
+
+using namespace std;
+
+class OPNode: public ActNode {
+public:
+	OPNode(OPCODE op,Node* left);
+	virtual ~OPNode();
+
+	TYPE getType();
+
+	OPCODE op;
+	Node* left;
+	double val;
+
+
+
+private:
+
+};
+
+}
+
+#endif /* OPNODE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/PNode.cpp b/src/boost/libs/yap/example/autodiff_library/PNode.cpp
new file mode 100644
index 00000000..ee07a1f2
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/PNode.cpp
@@ -0,0 +1,147 @@
+/*
+ * PNode.cpp
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#include "PNode.h"
+#include "Stack.h"
+#include "Tape.h"
+#include "EdgeSet.h"
+
+namespace AutoDiff {
+
+PNode::PNode(double value):pval(value) {
+	assert(!isnan(value));
+}
+
+PNode::~PNode() {
+}
+
+void PNode::inorder_visit(int level,ostream& oss){
+	oss<<this->toString(level)<<endl;
+}
+
+void PNode::collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total)
+{
+	//do not fill this to nodes, as this is a parameter node
+	total++;
+}
+
+void PNode::eval_function()
+{
+	SV->push_back(pval);
+}
+
+string PNode::toString(int level)
+{
+	ostringstream oss;
+	string s(level,'\t');
+	oss<<s<<"[PNode]("<<pval<<")";
+	return oss.str();
+}
+
+void PNode::grad_reverse_0()
+{
+	SV->push_back(pval);
+}
+void PNode::grad_reverse_1_init_adj()
+{
+	//do nothing as PNode does not have adjoint
+}
+
+void PNode::grad_reverse_1()
+{
+	//do nothing
+	//this is a parameter
+}
+
+void PNode::update_adj(double& v)
+{
+	//do nothing
+	//no adj for PNode
+}
+
+unsigned int PNode::hess_reverse_0()
+{
+	if(index==0)
+	{
+		TT->set(pval);
+		assert(TT->index == TT->index);
+		index = TT->index;
+	}
+	return index;
+}
+
+void PNode::hess_reverse_0_get_values(unsigned int i,double& x,double& x_bar,double& w,double& w_bar)
+{
+	x = TT->get(--i);
+	x_bar = 0;
+	w = 0;
+	w_bar = 0;
+}
+
+void PNode::hess_reverse_1(unsigned int i)
+{
+	n_in_arcs--;
+	//leaf node do nothing
+}
+
+void PNode::hess_reverse_1_init_x_bar(unsigned int)
+{
+	//do nothing as Parameter does not have x_bar
+}
+
+void PNode::update_x_bar(unsigned int i ,double v)
+{
+	//do nothing as Parameter does not have x_bar
+}
+void PNode::update_w_bar(unsigned int i ,double v)
+{
+	//do nothing as Parameter does not have w_bar
+}
+void PNode::hess_reverse_1_get_xw(unsigned int i, double& w,double& x)
+{
+	//do nothing as Parameter does not have w
+	x = TT->get(i-1);
+	w = 0;
+}
+void PNode::hess_reverse_get_x(unsigned int i, double& x)
+{
+	x = TT->get(i-1);
+}
+
+void PNode::nonlinearEdges(EdgeSet& edges)
+{
+	for(std::list<Edge>::iterator it=edges.edges.begin();it!=edges.edges.end();)
+	{
+		Edge e = *it;
+		if(e.a == this || e.b == this)
+		{
+			it = edges.edges.erase(it); //erase invalidate the iterator
+		}
+		else
+		{
+			it++;
+		}
+	}
+}
+
+#if FORWARD_ENABLED
+void PNode::hess_forward(unsigned int len, double** ret_vec)
+{
+	//it's a scalar
+	(*ret_vec) = new double[len];
+	std::fill_n(*ret_vec,len,0);
+	SV->push_back(this->pval);
+	assert(SV->size()==1);
+}
+#endif
+
+
+TYPE PNode::getType()
+{
+	return PNode_Type;
+}
+}
diff --git a/src/boost/libs/yap/example/autodiff_library/PNode.h b/src/boost/libs/yap/example/autodiff_library/PNode.h
new file mode 100644
index 00000000..58924d64
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/PNode.h
@@ -0,0 +1,47 @@
+/*
+ * PNode.h
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef PNODE_H_
+#define PNODE_H_
+
+#include "Node.h"
+namespace AutoDiff {
+
+using namespace std;
+class PNode: public Node {
+public:
+	PNode(double value);
+	virtual ~PNode();
+	void collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total);
+	void eval_function();
+	void grad_reverse_0();
+	void grad_reverse_1_init_adj();
+	void grad_reverse_1();
+	void update_adj(double& v);
+	void hess_forward(unsigned int len,double** ret_vec);
+	unsigned int hess_reverse_0();
+	void hess_reverse_0_get_values(unsigned int i,double& x,double& x_bar,double& w,double& w_bar);
+	void hess_reverse_1(unsigned int i);
+	void hess_reverse_1_init_x_bar(unsigned int);
+	void update_x_bar(unsigned int, double);
+	void update_w_bar(unsigned int, double);
+	void hess_reverse_1_get_xw(unsigned int, double&,double&);
+	void hess_reverse_get_x(unsigned int,double& x);
+	
+	void nonlinearEdges(EdgeSet&);
+
+	void inorder_visit(int level,ostream& oss);
+	string toString(int level);
+	TYPE getType();
+
+	double pval;
+
+};
+
+} // end namespace foo
+
+#endif /* PNODE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/README.md b/src/boost/libs/yap/example/autodiff_library/README.md
new file mode 100644
index 00000000..8f375446
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/README.md
@@ -0,0 +1,10 @@
+autodiff_library
+================ <\br>
+
+Automatic Differentiation Library <br>
+1. Very easy access and light weight c++/c library <br>
+2. Build computation graph as a DAG. <br>
+3. Provide function evaluation, reverse gradient, reverse Hessian-vector, forward Hessian computation calls.<br>
+4. State-of-the-art Object-Oriented Design and Implementation in C++ <br>
+5. A modified version of this library is integrated into the Parallel Structured Model Generator(PSMG)--an algebraic modelling langauge for Mathematical Programming.
+See project wiki for more detail.
diff --git a/src/boost/libs/yap/example/autodiff_library/Stack.cpp b/src/boost/libs/yap/example/autodiff_library/Stack.cpp
new file mode 100644
index 00000000..f0e69f84
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Stack.cpp
@@ -0,0 +1,58 @@
+/*
+ * Stack.cpp
+ *
+ *  Created on: 15 Apr 2013
+ *      Author: s0965328
+ */
+
+
+#include <cstddef>
+#include <math.h>
+#include <cassert>
+
+#include "Stack.h"
+
+namespace AutoDiff {
+
+
+Stack* Stack::vals = NULL;
+Stack* Stack::diff = NULL;
+
+Stack::Stack()
+{
+}
+
+Stack::~Stack() {
+	this->clear();
+}
+
+
+double Stack::pop_back()
+{
+	assert(this->lifo.size()!=0);
+	double v = this->lifo.top();
+	lifo.pop();
+	return v;
+}
+void Stack::push_back(double& v)
+{
+	assert(!isnan(v));
+	this->lifo.push(v);
+}
+double& Stack::peek()
+{
+	return this->lifo.top();
+}
+unsigned int Stack::size()
+{
+	return this->lifo.size();
+}
+
+void Stack::clear()
+{
+	while(!this->lifo.empty())
+	{
+		this->lifo.pop();
+	}
+}
+}
diff --git a/src/boost/libs/yap/example/autodiff_library/Stack.h b/src/boost/libs/yap/example/autodiff_library/Stack.h
new file mode 100644
index 00000000..30f5a3f0
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Stack.h
@@ -0,0 +1,38 @@
+/*
+ * Stack.h
+ *
+ *  Created on: 15 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef STACK_H_
+#define STACK_H_
+
+#include <stack>
+
+namespace AutoDiff {
+
+using namespace std;
+#define SV (Stack::vals)
+#define SD (Stack::diff)
+
+class Stack {
+public:
+	Stack();
+	double pop_back();
+	void push_back(double& v);
+	double& peek();
+	unsigned int size();
+	void clear();
+	virtual ~Stack();
+
+	stack<double> lifo;
+
+	static Stack* diff;
+	static Stack* vals;
+
+
+};
+
+}
+#endif /* STACK_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/Tape.cpp b/src/boost/libs/yap/example/autodiff_library/Tape.cpp
new file mode 100644
index 00000000..0a5682f3
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Tape.cpp
@@ -0,0 +1,16 @@
+/*
+ * Tape.cpp
+ *
+ *  Created on: 6 Nov 2013
+ *      Author: s0965328
+ */
+
+
+
+#include "Tape.h"
+
+namespace AutoDiff
+{
+	template<> Tape<unsigned int>* Tape<unsigned int>::indexTape = NULL;
+	template<> Tape<double>* Tape<double>::valueTape = NULL;
+}
diff --git a/src/boost/libs/yap/example/autodiff_library/Tape.h b/src/boost/libs/yap/example/autodiff_library/Tape.h
new file mode 100644
index 00000000..9ca62c4d
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/Tape.h
@@ -0,0 +1,97 @@
+/*
+ * Tape.h
+ *
+ *  Created on: 15 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef TAPE_H_
+#define TAPE_H_
+
+#include <vector>
+#include <string>
+#include <cassert>
+#include <sstream>
+#include <iostream>
+
+
+namespace AutoDiff {
+
+using namespace std;
+#define TT 	(Tape<double>::valueTape)
+#define II 	(Tape<unsigned int>::indexTape)
+
+template<typename T> class Tape {
+public:
+	Tape<T> () : index(0){};
+	T& at(const unsigned int index);
+	const T& get(const unsigned int index);
+	void set(T& v);
+	unsigned int size();
+	void clear();
+	bool empty();
+	string toString();
+	virtual ~Tape();
+
+	vector<T> vals;
+	unsigned int index;
+
+	static Tape<double>* valueTape;
+	static Tape<unsigned int>* indexTape;
+};
+
+
+template<typename T> Tape<T>::~Tape<T>()
+{
+	index = 0;
+	vals.clear();
+}
+
+template<typename T> T& Tape<T>::at(const unsigned int i)
+{
+	assert(this->vals.size()>i);
+	return vals[i];
+}
+template<typename T> const T& Tape<T>::get(const unsigned int i)
+{
+	assert(this->vals.size()>i);
+	return vals[i];
+}
+template <typename T> void Tape<T>::set(T& v)
+{
+	vals.push_back(v);
+	index++;
+}
+
+template<typename T> unsigned int Tape<T>::size()
+{
+	return this->vals.size();
+}
+
+template<typename T> bool Tape<T>::empty()
+{
+	return vals.empty();
+}
+
+template<typename T> void Tape<T>::clear()
+{
+	this->vals.clear();
+	this->index = 0;
+	assert(this->vals.size()==0);
+	assert(this->vals.empty());
+}
+
+template<typename T>  string Tape<T>::toString()
+{
+	assert(vals.size()>=index);
+	ostringstream oss;
+	oss<<"Tape size["<<vals.size()<<"]";
+	for(unsigned int i=0;i<vals.size();i++){
+		if(i%10==0) oss<<endl;
+		oss<<vals[i]<<",";
+	}
+	oss<<endl;
+	return oss.str();
+}
+}
+#endif /* TAPE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/UaryOPNode.cpp b/src/boost/libs/yap/example/autodiff_library/UaryOPNode.cpp
new file mode 100644
index 00000000..73e2711c
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/UaryOPNode.cpp
@@ -0,0 +1,375 @@
+/*
+ * UaryOPNode.cpp
+ *
+ *  Created on: 6 Nov 2013
+ *      Author: s0965328
+ */
+
+#include "UaryOPNode.h"
+#include "BinaryOPNode.h"
+#include "PNode.h"
+#include "Stack.h"
+#include "Tape.h"
+#include "Edge.h"
+#include "EdgeSet.h"
+#include "auto_diff_types.h"
+
+#include <list>
+
+using namespace std;
+
+namespace AutoDiff {
+
+UaryOPNode::UaryOPNode(OPCODE op_, Node* left): OPNode(op_,left) {
+}
+
+OPNode* UaryOPNode::createUnaryOpNode(OPCODE op, Node* left)
+{
+	assert(left!=NULL);
+	OPNode* node = NULL;
+	if(op == OP_SQRT)
+	{
+		double param = 0.5;
+		node = BinaryOPNode::createBinaryOpNode(OP_POW,left,new PNode(param));
+	}
+	else if(op == OP_NEG)
+	{
+		double param = -1;
+		node = BinaryOPNode::createBinaryOpNode(OP_TIMES,left,new PNode(param));
+	}
+	else
+	{
+		node = new UaryOPNode(op,left);
+	}
+	return node;
+}
+
+UaryOPNode::~UaryOPNode() {
+
+}
+
+
+void UaryOPNode::inorder_visit(int level,ostream& oss){
+	if(left!=NULL){
+		left->inorder_visit(level+1,oss);
+	}
+	oss<<this->toString(level)<<endl;
+}
+
+void UaryOPNode::collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total)
+{
+	total++;
+	if(left!=NULL){
+		left->collect_vnodes(nodes,total);
+	}
+}
+
+void UaryOPNode::eval_function()
+{
+	if(left!=NULL){
+		left->eval_function();
+	}
+	this->calc_eval_function();
+}
+
+//1. visiting left if not NULL
+//2. then, visiting right if not NULL
+//3. calculating the immediate derivative hu and hv
+void UaryOPNode::grad_reverse_0(){
+	assert(left!=NULL);
+	this->adj = 0;
+	left->grad_reverse_0();
+	this->calc_grad_reverse_0();
+}
+
+//right left - right most traversal
+void UaryOPNode::grad_reverse_1()
+{
+	assert(left!=NULL);
+	double l_adj = SD->pop_back()*this->adj;
+	left->update_adj(l_adj);
+	left->grad_reverse_1();
+}
+
+void UaryOPNode::calc_grad_reverse_0()
+{
+	assert(left!=NULL);
+	double hu = NaN_Double;
+	double lval = SV->pop_back();
+	double val = NaN_Double;
+	switch (op)
+	{
+	case OP_SIN:
+		val = sin(lval);
+		hu = cos(lval);
+		break;
+	case OP_COS:
+		val = cos(lval);
+		hu = -sin(lval);
+		break;
+	default:
+		cerr<<"error op not impl"<<endl;
+		break;
+	}
+	SV->push_back(val);
+	SD->push_back(hu);
+}
+
+void UaryOPNode::calc_eval_function()
+{
+	double lval = SV->pop_back();
+	double val = NaN_Double;
+	switch (op)
+	{
+	case OP_SIN:
+		assert(left!=NULL);
+		val = sin(lval);
+		break;
+	case OP_COS:
+		assert(left!=NULL);
+		val = cos(lval);
+		break;
+	default:
+		cerr<<"op["<<op<<"] not yet implemented!!"<<endl;
+		assert(false);
+		break;
+	}
+	SV->push_back(val);
+}
+
+void UaryOPNode::hess_reverse_0_init_n_in_arcs()
+{
+	this->left->hess_reverse_0_init_n_in_arcs();
+	this->Node::hess_reverse_0_init_n_in_arcs();
+}
+
+void UaryOPNode::hess_reverse_1_clear_index()
+{
+	this->left->hess_reverse_1_clear_index();
+	this->Node::hess_reverse_1_clear_index();
+}
+
+unsigned int UaryOPNode::hess_reverse_0()
+{
+	assert(left!=NULL);
+	if(index==0)
+	{
+		unsigned int lindex=0;
+		lindex = left->hess_reverse_0();
+		assert(lindex!=0);
+		II->set(lindex);
+		double lx,lx_bar,lw,lw_bar;
+		double x,x_bar,w,w_bar;
+		double l_dh;
+		switch(op)
+		{
+		case OP_SIN:
+			assert(left != NULL);
+			left->hess_reverse_0_get_values(lindex,lx,lx_bar,lw,lw_bar);
+			x = sin(lx);
+			x_bar = 0;
+			l_dh = cos(lx);
+			w = lw*l_dh;
+			w_bar = 0;
+			break;
+		case OP_COS:
+			assert(left!=NULL);
+			left->hess_reverse_0_get_values(lindex,lx,lx_bar,lw,lw_bar);
+			x = cos(lx);
+			x_bar = 0;
+			l_dh = -sin(lx);
+			w = lw*l_dh;
+			w_bar = 0;
+			break;
+		default:
+			cerr<<"op["<<op<<"] not yet implemented!"<<endl;
+			assert(false);
+			break;
+		}
+		TT->set(x);
+		TT->set(x_bar);
+		TT->set(w);
+		TT->set(w_bar);
+		TT->set(l_dh);
+		assert(TT->index == TT->index);
+		index = TT->index;
+	}
+	return index;
+}
+
+void UaryOPNode::hess_reverse_0_get_values(unsigned int i,double& x, double& x_bar, double& w, double& w_bar)
+{
+	--i; // skip the l_dh (ie, dh/du)
+	w_bar = TT->get(--i);
+	w = TT->get(--i);
+	x_bar = TT->get(--i);
+	x = TT->get(--i);
+}
+
+void UaryOPNode::hess_reverse_1(unsigned int i)
+{
+	n_in_arcs--;
+	if(n_in_arcs==0)
+	{
+		double lindex = II->get(--(II->index));
+	//	cout<<"li["<<lindex<<"]\t"<<this->toString(0)<<endl;
+		double l_dh = TT->get(--i);
+		double w_bar = TT->get(--i);
+		--i; //skip w
+		double x_bar = TT->get(--i);
+		--i; //skip x
+	//	cout<<"i["<<i<<"]"<<endl;
+
+		assert(left!=NULL);
+		left->update_x_bar(lindex,x_bar*l_dh);
+		double lw_bar = 0;
+		double lw = 0,lx = 0;
+		left->hess_reverse_1_get_xw(lindex,lw,lx);
+		switch(op)
+		{
+		case OP_SIN:
+			assert(l_dh == cos(lx));
+			lw_bar += w_bar*l_dh;
+			lw_bar += x_bar*lw*(-sin(lx));
+			break;
+		case OP_COS:
+			assert(l_dh == -sin(lx));
+			lw_bar += w_bar*l_dh;
+			lw_bar += x_bar*lw*(-cos(lx));
+			break;
+		default:
+			cerr<<"op["<<op<<"] not yet implemented!"<<endl;
+			break;
+		}
+		left->update_w_bar(lindex,lw_bar);
+		left->hess_reverse_1(lindex);
+	}
+}
+
+void UaryOPNode::hess_reverse_1_init_x_bar(unsigned int i)
+{
+	TT->at(i-4) = 1;
+}
+
+void UaryOPNode::update_x_bar(unsigned int i ,double v)
+{
+	TT->at(i-4) += v;
+}
+void UaryOPNode::update_w_bar(unsigned int i ,double v)
+{
+	TT->at(i-2) += v;
+}
+void UaryOPNode::hess_reverse_1_get_xw(unsigned int i,double& w,double& x)
+{
+	w = TT->get(i-3);
+	x = TT->get(i-5);
+}
+void UaryOPNode::hess_reverse_get_x(unsigned int i, double& x)
+{
+	x = TT->get(i-5);
+}
+
+void UaryOPNode::nonlinearEdges(EdgeSet& edges)
+{
+	for(list<Edge>::iterator it=edges.edges.begin();it!=edges.edges.end();)
+	{
+		Edge& e = *it;
+		if(e.a == this || e.b == this){
+			if(e.a == this && e.b == this)
+			{
+				Edge e1(left,left);
+				edges.insertEdge(e1);
+			}
+			else{
+				Node* o = e.a==this?e.b:e.a;
+				Edge e1(left,o);
+				edges.insertEdge(e1);
+			}
+			it = edges.edges.erase(it);
+		}
+		else
+		{
+			it++;
+		}
+	}
+
+	Edge e1(left,left);
+	switch(op)
+	{
+	case OP_SIN:
+		edges.insertEdge(e1);
+		break;
+	case OP_COS:
+		edges.insertEdge(e1);
+		break;
+	default:
+		cerr<<"op["<<op<<"] is not yet implemented !"<<endl;
+		assert(false);
+		break;
+	}
+	left->nonlinearEdges(edges);
+}
+
+#if FORWARD_ENABLED
+void UaryOPNode::hess_forward(unsigned int len, double** ret_vec)
+{
+	double* lvec = NULL;
+	if(left!=NULL){
+		left->hess_forward(len,&lvec);
+	}
+
+	*ret_vec = new double[len];
+	this->hess_forward_calc0(len,lvec,*ret_vec);
+	delete[] lvec;
+}
+
+void UaryOPNode::hess_forward_calc0(unsigned int& len, double* lvec, double* ret_vec)
+{
+	double hu = NaN_Double;
+	double lval = NaN_Double;
+	double val = NaN_Double;
+	unsigned int index = 0;
+	switch (op)
+	{
+	case OP_SIN:
+		assert(left!=NULL);
+		lval = SV->pop_back();
+		val = sin(lval);
+		SV->push_back(val);
+		hu = cos(lval);
+
+		double coeff;
+		coeff = -val; //=sin(left->val); -- and avoid cross initialisation
+		//calculate the first order derivatives
+		for(unsigned int i =0;i<AutoDiff::num_var;++i)
+		{
+			ret_vec[i] = hu*lvec[i] + 0;
+		}
+		//calculate the second order
+		index = AutoDiff::num_var;
+		for(unsigned int i=0;i<AutoDiff::num_var;++i)
+		{
+			for(unsigned int j=i;j<AutoDiff::num_var;++j)
+			{
+				ret_vec[index] = hu*lvec[index] + lvec[i] * coeff * lvec[j] + 0 + 0;
+				++index;
+			}
+		}
+		assert(index==len);
+		break;
+	default:
+		cerr<<"op["<<op<<"] not yet implemented!";
+		break;
+	}
+}
+#endif
+
+string UaryOPNode::toString(int level)
+{
+	ostringstream oss;
+	string s(level,'\t');
+	oss<<s<<"[UaryOPNode]("<<op<<")";
+	return oss.str();
+}
+
+} /* namespace AutoDiff */
diff --git a/src/boost/libs/yap/example/autodiff_library/UaryOPNode.h b/src/boost/libs/yap/example/autodiff_library/UaryOPNode.h
new file mode 100644
index 00000000..241cd992
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/UaryOPNode.h
@@ -0,0 +1,52 @@
+/*
+ * UaryOPNode.h
+ *
+ *  Created on: 6 Nov 2013
+ *      Author: s0965328
+ */
+
+#ifndef UARYOPNODE_H_
+#define UARYOPNODE_H_
+
+#include "OPNode.h"
+
+namespace AutoDiff {
+
+class UaryOPNode: public OPNode {
+public:
+	static OPNode* createUnaryOpNode(OPCODE op, Node* left);
+	virtual ~UaryOPNode();
+
+	void inorder_visit(int level,ostream& oss);
+	void collect_vnodes(boost::unordered_set<Node*> & nodes,unsigned int& total);
+	void eval_function();
+
+	void grad_reverse_0();
+	void grad_reverse_1();
+#if FORWARD_ENABLED
+	void hess_forward(unsigned int len, double** ret_vec);
+#endif
+	unsigned int hess_reverse_0();
+	void hess_reverse_0_init_n_in_arcs();
+	void hess_reverse_0_get_values(unsigned int i,double& x, double& x_bar, double& w, double& w_bar);
+	void hess_reverse_1(unsigned int i);
+	void hess_reverse_1_init_x_bar(unsigned int);
+	void update_x_bar(unsigned int, double v);
+	void update_w_bar(unsigned int, double v);
+	void hess_reverse_1_get_xw(unsigned int, double&,double&);
+	void hess_reverse_get_x(unsigned int,double& x);
+	void hess_reverse_1_clear_index();
+
+	void nonlinearEdges(EdgeSet&);
+
+	string toString(int level);
+
+private:
+	UaryOPNode(OPCODE op, Node* left);
+	void calc_eval_function();
+	void calc_grad_reverse_0();
+	void hess_forward_calc0(unsigned int& len, double* lvec,double* ret_vec);
+};
+
+} /* namespace AutoDiff */
+#endif /* UARYOPNODE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/VNode.cpp b/src/boost/libs/yap/example/autodiff_library/VNode.cpp
new file mode 100644
index 00000000..285af01f
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/VNode.cpp
@@ -0,0 +1,149 @@
+/*
+ * VNode.cpp
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#include "VNode.h"
+#include "Tape.h"
+#include "Stack.h"
+
+using namespace std;
+
+namespace AutoDiff {
+
+#if FORWARD_ENABLED
+int VNode::DEFAULT_ID = -1;
+#endif
+
+VNode::VNode(double v) : ActNode(), val(v),u(NaN_Double)
+#if FORWARD_ENABLED
+,id(DEFAULT_ID)
+#endif
+{
+}
+
+VNode::~VNode() {
+}
+
+void VNode::collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total)
+{
+	total++;
+	boost::unordered_set<Node*>::iterator it = nodes.find(this);
+	if(it==nodes.end())
+		nodes.insert(this);
+}
+
+void VNode::inorder_visit(int level,ostream& oss)
+{
+	oss<<this->toString(level)<<endl;
+}
+
+void VNode::eval_function()
+{
+	SV->push_back(val);
+}
+
+string VNode::toString(int level)
+{
+	ostringstream oss;
+	string s(level,'\t');
+	oss<<s<<"[VNode](index:"<<index<<",val:"<<val<<",u:"<<u<<") - "<<this;
+	return oss.str();
+}
+
+void VNode::grad_reverse_0()
+{
+	this->adj = 0;
+	SV->push_back(val);
+}
+
+void VNode::grad_reverse_1()
+{
+	//do nothing
+	//this is a leaf node
+}
+
+#if FORWARD_ENABLED
+void VNode::hess_forward(unsigned int len, double** ret_vec)
+{
+	assert(id!=DEFAULT_ID);
+	(*ret_vec) = new double[len];
+	std::fill_n(*ret_vec,len,0);
+	(*ret_vec)[id]=1;
+	SV->push_back(this->val);
+}
+#endif
+
+unsigned int VNode::hess_reverse_0()
+{
+	if(index==0)
+	{//this node is not on tape
+		double nan = NaN_Double;
+		TT->set(val);  	//x_i
+		TT->set(nan);		//x_bar_i
+		TT->set(u);    	//w_i
+		TT->set(nan);    	//w_bar_i
+		index = TT->index;
+	}
+//	cout<<toString(0)<<"  -- "<<index<<endl;
+	return index;
+}
+
+void VNode::hess_reverse_0_get_values(unsigned int i,double& x, double& x_bar, double& w, double& w_bar)
+{
+	w_bar = TT->get(--i);
+	w = TT->get(--i);
+	x_bar = TT->get(--i);
+	x = TT->get(--i);
+}
+
+void VNode::hess_reverse_1(unsigned int i)
+{
+	n_in_arcs--;
+	//leaf node do nothing
+}
+
+void VNode::hess_reverse_1_init_x_bar(unsigned int i)
+{
+	TT->at(i-3) = 1;
+}
+
+void VNode::update_x_bar(unsigned int i ,double v)
+{
+//	cout<<toString(0)<<"  --- "<<__FUNCTION__<<"  v="<<TT->at(i-3)<<"+"<<v<<endl;
+	TT->at(i-3) = isnan(TT->get(i-3))? v: TT->get(i-3) + v;
+}
+
+void VNode::update_w_bar(unsigned int i,double v)
+{
+//	cout<<toString(0)<<"  --- "<<__FUNCTION__<<"  v="<<TT->at(i-1)<<"+"<<v<<endl;
+	TT->at(i-1) = isnan(TT->get(i-1))? v: TT->get(i-1) + v;
+}
+void VNode::hess_reverse_1_get_xw(unsigned int i, double& w,double& x)
+{
+	//cout<<toString(0)<<"  --- "<<__FUNCTION__<<"  w="<<TT->get(i-2)<<"-- "<<"x="<<TT->get(i-4)<<endl;
+	w = TT->get(i-2);
+	x = TT->get(i-4);
+}
+void VNode::hess_reverse_get_x(unsigned int i ,double& x)
+{
+	x = TT->get(i-4);
+}
+
+void VNode::nonlinearEdges(EdgeSet& edges)
+{
+//	for(list<Edge>::iterator it = edges.edges.begin();it!=edges.edges.end();)
+//	{
+//		Edge e=*it;
+//
+//	}
+}
+
+
+TYPE VNode::getType()
+{
+	return VNode_Type;
+}
+}
diff --git a/src/boost/libs/yap/example/autodiff_library/VNode.h b/src/boost/libs/yap/example/autodiff_library/VNode.h
new file mode 100644
index 00000000..d9673d4b
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/VNode.h
@@ -0,0 +1,52 @@
+/*
+ * VNode.h
+ *
+ *  Created on: 8 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef VNODE_H_
+#define VNODE_H_
+
+#include "ActNode.h"
+
+namespace AutoDiff {
+class VNode: public ActNode {
+public:
+	VNode(double v=NaN_Double);
+	virtual ~VNode();
+	void collect_vnodes(boost::unordered_set<Node*>& nodes,unsigned int& total);
+	void eval_function();
+	void grad_reverse_0();
+	void grad_reverse_1();
+	unsigned int hess_reverse_0();
+	void hess_reverse_0_get_values(unsigned int i,double& x, double& x_bar, double& w, double& w_bar);
+	void hess_reverse_1(unsigned int i);
+	void hess_reverse_1_init_x_bar(unsigned int);
+	void update_x_bar(unsigned int,double);
+	void update_w_bar(unsigned int,double);
+	void hess_reverse_1_get_xw(unsigned int,double&,double&);
+	void hess_reverse_get_x(unsigned int,double& x);
+
+	void nonlinearEdges(EdgeSet&);
+
+	void inorder_visit(int level,ostream& oss);
+	string toString(int level);
+	TYPE getType();
+
+#if FORWARD_ENABLED
+	void hess_forward(unsigned int nvar, double** ret_vec);
+	//! used for only forward hessian
+	//! the id has to be assigned starting from 0
+	int id;
+	static int DEFAULT_ID;
+#endif
+	double val;
+	double u;
+
+
+};
+
+} // end namespace foo
+
+#endif /* VNODE_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/auto_diff_types.h b/src/boost/libs/yap/example/autodiff_library/auto_diff_types.h
new file mode 100644
index 00000000..b24c59a1
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/auto_diff_types.h
@@ -0,0 +1,26 @@
+/*
+ * auto_diff_types.h
+ *
+ *  Created on: 18 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef AUTO_DIFF_TYPES_H_
+#define AUTO_DIFF_TYPES_H_
+
+namespace AutoDiff{
+
+#define FORWARD_ENABLED 0
+#if FORWARD_ENABLED
+extern unsigned int num_var;
+#endif
+
+#define NaN_Double std::numeric_limits<double>::quiet_NaN()
+
+typedef enum { OPNode_Type=0, VNode_Type, PNode_Type} TYPE;
+
+
+typedef enum {OP_PLUS=0, OP_MINUS, OP_TIMES, OP_DIVID, OP_SIN, OP_COS, OP_SQRT, OP_POW, OP_NEG} OPCODE;
+
+}
+#endif /* AUTO_DIFF_TYPES_H_ */
diff --git a/src/boost/libs/yap/example/autodiff_library/autodiff.cpp b/src/boost/libs/yap/example/autodiff_library/autodiff.cpp
new file mode 100644
index 00000000..5fae3686
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/autodiff.cpp
@@ -0,0 +1,335 @@
+//============================================================================
+// Name        : autodiff.cpp
+// Author      : 
+// Version     :
+// Copyright   : Your copyright notice
+// Description : Hello World in C++, Ansi-style
+//============================================================================
+
+#include <iostream>
+#include <sstream>
+#include <numeric>
+#include <boost/foreach.hpp>
+#include "autodiff.h"
+#include "Stack.h"
+#include "Tape.h"
+#include "BinaryOPNode.h"
+#include "UaryOPNode.h"
+
+using namespace std;
+
+namespace AutoDiff
+{
+
+#if FORWARD_ENABLED
+
+unsigned int num_var = 0;
+
+void hess_forward(Node* root, unsigned int nvar, double** hess_mat)
+{
+	assert(nvar == num_var);
+	unsigned int len = (nvar+3)*nvar/2;
+	root->hess_forward(len,hess_mat);
+}
+
+#endif
+
+
+PNode* create_param_node(double value){
+	return new PNode(value);
+}
+VNode* create_var_node(double v)
+{
+	return new VNode(v);
+}
+OPNode* create_binary_op_node(OPCODE code, Node* left, Node* right)
+{
+	return BinaryOPNode::createBinaryOpNode(code,left,right);
+}
+OPNode* create_uary_op_node(OPCODE code, Node* left)
+{
+	return UaryOPNode::createUnaryOpNode(code,left);
+}
+double eval_function(Node* root)
+{
+	assert(SD->size()==0);
+	assert(SV->size()==0);
+	root->eval_function();
+	assert(SV->size()==1);
+	double val = SV->pop_back();
+	return val;
+}
+
+double grad_reverse(Node* root,vector<Node*>& vnodes, vector<double>& grad)
+{
+	grad.clear();
+	BOOST_FOREACH(Node* node, vnodes)
+	{
+		assert(node->getType()==VNode_Type);
+		static_cast<VNode*>(node)->adj = NaN_Double;
+	}
+
+	assert(SD->size()==0);
+	root->grad_reverse_0();
+	assert(SV->size()==1);
+	root->grad_reverse_1_init_adj();	
+	root->grad_reverse_1();
+	assert(SD->size()==0);
+	double val = SV->pop_back();
+	assert(SV->size()==0);
+	//int i=0;
+	BOOST_FOREACH(Node* node, vnodes)
+	{
+		assert(node->getType()==VNode_Type);
+		grad.push_back(static_cast<VNode*>(node)->adj);
+		static_cast<VNode*>(node)->adj = NaN_Double;
+	}
+	assert(grad.size()==vnodes.size());
+	//all nodes are VNode and adj == NaN_Double -- this reset adj for this expression tree by root
+	return val;
+}
+
+double grad_reverse(Node* root, vector<Node*>& vnodes, col_compress_matrix_row& rgrad)
+{
+	BOOST_FOREACH(Node* node, vnodes)
+	{
+		assert(node->getType()==VNode_Type);
+		static_cast<VNode*>(node)->adj = NaN_Double;
+	}
+	assert(SD->size()==0);
+	root->grad_reverse_0();
+	assert(SV->size()==1);
+	root->grad_reverse_1_init_adj();
+	root->grad_reverse_1();
+	assert(SD->size()==0);
+	double val = SV->pop_back();
+	assert(SV->size()==0);
+	unsigned int i =0;
+	BOOST_FOREACH(Node* node, vnodes)
+	{
+		assert((node)->getType()==VNode_Type);
+		double diff = static_cast<VNode*>(node)->adj;
+		if(!isnan(diff)){
+			rgrad(i) = diff;
+			static_cast<VNode*>(node)->adj = NaN_Double;
+		}
+		i++;
+	}
+	//all nodes are VNode and adj == NaN_Double -- this reset adj for this expression tree by root
+	assert(i==vnodes.size());
+	return val;
+}
+
+double hess_reverse(Node* root,vector<Node*>& vnodes,vector<double>& dhess)
+{
+	TT->clear();
+	II->clear();
+	assert(TT->empty());
+	assert(II->empty());
+	assert(TT->index==0);
+	assert(II->index==0);
+	dhess.clear();
+
+//	for(vector<Node*>::iterator it=nodes.begin();it!=nodes.end();it++)
+//	{
+//		assert((*it)->getType()==VNode_Type);
+//		(*it)->index = 0;
+//	} //this work complete in hess-reverse_0_init_index
+
+	assert(root->n_in_arcs == 0);
+	root->hess_reverse_0_init_n_in_arcs();
+	assert(root->n_in_arcs == 1);
+	root->hess_reverse_0();
+	double val = NaN_Double;
+	root->hess_reverse_get_x(TT->index,val);
+//	cout<<TT->toString();
+//	cout<<endl;
+//	cout<<II->toString();
+//	cout<<"======================================= hess_reverse_0"<<endl;
+	root->hess_reverse_1_init_x_bar(TT->index);
+	assert(root->n_in_arcs == 1);
+	root->hess_reverse_1(TT->index);
+	assert(root->n_in_arcs == 0);
+	assert(II->index==0);
+//	cout<<TT->toString();
+//	cout<<endl;
+//	cout<<II->toString();
+//	cout<<"======================================= hess_reverse_1"<<endl;
+
+	for(vector<Node*>::iterator it=vnodes.begin();it!=vnodes.end();it++)
+	{
+		assert((*it)->getType()==VNode_Type);
+		dhess.push_back(TT->get((*it)->index-1));
+	}
+
+	TT->clear();
+	II->clear();
+	root->hess_reverse_1_clear_index();
+	return val;
+}
+
+double hess_reverse(Node* root,vector<Node*>& vnodes,col_compress_matrix_col& chess)
+{
+	TT->clear();
+	II->clear();
+	assert(TT->empty());
+	assert(II->empty());
+	assert(TT->index==0);
+	assert(II->index==0);
+
+//	for(vector<Node*>::iterator it=nodes.begin();it!=nodes.end();it++)
+//	{
+//		assert((*it)->getType()==VNode_Type);
+//		(*it)->index = 0;
+//	} //this work complete in hess-reverse_0_init_index
+
+	assert(root->n_in_arcs == 0);
+	//reset node index and n_in_arcs - for the Tape location
+	root->hess_reverse_0_init_n_in_arcs();
+	assert(root->n_in_arcs == 1);
+	root->hess_reverse_0();
+	double val = NaN_Double;
+	root->hess_reverse_get_x(TT->index,val);
+//	cout<<TT->toString();
+//	cout<<endl;
+//	cout<<II->toString();
+//	cout<<"======================================= hess_reverse_0"<<endl;
+	root->hess_reverse_1_init_x_bar(TT->index);
+	assert(root->n_in_arcs == 1);
+	root->hess_reverse_1(TT->index);
+	assert(root->n_in_arcs == 0);
+	assert(II->index==0);
+//	cout<<TT->toString();
+//	cout<<endl;
+//	cout<<II->toString();
+//	cout<<"======================================= hess_reverse_1"<<endl;
+
+	unsigned int i =0;
+	BOOST_FOREACH(Node* node, vnodes)
+	{
+		assert(node->getType() == VNode_Type);
+		//node->index = 0 means this VNode is not in the tree
+		if(node->index!=0)
+		{
+			double hess = TT->get(node->index -1);
+			if(!isnan(hess))
+			{
+				chess(i) = chess(i) + hess;
+			}
+		}
+		i++;
+	}
+	assert(i==vnodes.size());
+	root->hess_reverse_1_clear_index();
+	TT->clear();
+	II->clear();
+	return val;
+}
+
+unsigned int nzGrad(Node* root)
+{
+	unsigned int nzgrad,total = 0;
+	boost::unordered_set<Node*> nodes;
+	root->collect_vnodes(nodes,total);
+	nzgrad = nodes.size();
+	return nzgrad;
+}
+
+/*
+ * number of non-zero gradient in constraint tree root that also belong to vSet
+ */
+unsigned int nzGrad(Node* root, boost::unordered_set<Node*>& vSet)
+{
+	unsigned int nzgrad=0, total=0;
+	boost::unordered_set<Node*> vnodes;
+	root->collect_vnodes(vnodes,total);
+	//cout<<"nzGrad - vnodes size["<<vnodes.size()<<"] -- total node["<<total<<"]"<<endl;
+	for(boost::unordered_set<Node*>::iterator it=vnodes.begin();it!=vnodes.end();it++)
+	{
+		Node* n = *it;
+		if(vSet.find(n) != vSet.end())
+		{
+			nzgrad++;
+		}
+	}
+	return nzgrad;
+}
+
+void nonlinearEdges(Node* root, EdgeSet& edges)
+{
+	root->nonlinearEdges(edges);
+}
+
+unsigned int nzHess(EdgeSet& eSet,boost::unordered_set<Node*>& set1, boost::unordered_set<Node*>& set2)
+{
+	list<Edge>::iterator i = eSet.edges.begin();
+	for(;i!=eSet.edges.end();)
+	{
+		Edge e =*i;
+		Node* a = e.a;
+		Node* b = e.b;
+		if((set1.find(a)!=set1.end() && set2.find(b)!=set2.end())
+			||
+			(set1.find(b)!=set1.end() && set2.find(a)!=set2.end()))
+		{
+			//e is connected between set1 and set2
+			i++;
+		}
+		else
+		{
+			i = eSet.edges.erase(i);
+		}
+	}
+	unsigned int diag=eSet.numSelfEdges();
+	unsigned int nzHess = (eSet.size())*2 - diag;
+	return nzHess;
+}
+
+unsigned int nzHess(EdgeSet& edges)
+{
+	unsigned int diag=edges.numSelfEdges();
+	unsigned int nzHess = (edges.size())*2 - diag;
+	return nzHess;
+}
+
+unsigned int numTotalNodes(Node* root)
+{
+	unsigned int total = 0;
+	boost::unordered_set<Node*> nodes;
+	root->collect_vnodes(nodes,total);
+	return total;
+}
+
+string tree_expr(Node* root)
+{
+	ostringstream oss;
+	oss<<"visiting tree == "<<endl;
+	int level = 0;
+	root->inorder_visit(level,oss);
+	return oss.str();
+}
+
+void print_tree(Node* root)
+{
+	cout<<"visiting tree == "<<endl;
+	int level = 0;
+	root->inorder_visit(level,cout);
+}
+
+void autodiff_setup()
+{
+	Stack::diff = new Stack();
+	Stack::vals = new Stack();
+	Tape<unsigned int>::indexTape = new Tape<unsigned int>();
+	Tape<double>::valueTape = new Tape<double>();
+}
+
+void autodiff_cleanup()
+{
+	delete Stack::diff;
+	delete Stack::vals;
+	delete Tape<unsigned int>::indexTape;
+	delete Tape<double>::valueTape;
+}
+
+} //AutoDiff namespace end
diff --git a/src/boost/libs/yap/example/autodiff_library/autodiff.h b/src/boost/libs/yap/example/autodiff_library/autodiff.h
new file mode 100644
index 00000000..738cb3f9
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_library/autodiff.h
@@ -0,0 +1,114 @@
+/*
+ * autodiff.h
+ *
+ *  Created on: 16 Apr 2013
+ *      Author: s0965328
+ */
+
+#ifndef AUTODIFF_H_
+#define AUTODIFF_H_
+#include <boost/unordered_set.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_sparse.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include "auto_diff_types.h"
+#include "Node.h"
+#include "VNode.h"
+#include "OPNode.h"
+#include "PNode.h"
+#include "ActNode.h"
+#include "EdgeSet.h"
+
+
+/*
+ * + Function and Gradient Evaluation
+ * The tapeless implementation for function and derivative evaluation
+ * Advantage for tapeless:
+ * 			Low memory usage
+ * 			Function evaluation use one stack
+ * 			Gradient evaluation use two stack.
+ * Disadvantage for tapeless:
+ * 			Inefficient if the expression tree have repeated nodes.
+ * 			for example:
+ * 						 root
+ * 						 /  \
+ * 						*    *
+ * 					   / \	/ \
+ * 				      x1 x1 x1 x1
+ * 				      Tapeless implementation will go through all the edges.
+ * 				      ie. adjoint of x will be updated 4 times for the correct
+ * 				      gradient of x1.While the tape implemenation can discovery this
+ * 				      dependence and update adjoint of x1 just twice. The computational
+ * 				      graph (DAG) for a taped implemenation will be look like bellow.
+ * 				       root
+ * 				        /\
+ * 				        *
+ * 				        /\
+ * 				        x1
+ *
+ * + Forward Hessian Evaluation:
+ * This is an inefficient implementation of the forward Hessian method. It will evaluate the diagonal
+ * and upper triangular of the Hessian. The gradient is also evaluation in the same routine. The result
+ * will be returned in an array.
+ * To use this method, one have to provide a len parameter. len = (nvar+3)*nvar/2 where nvar is the number
+ * of independent variables. ie. x_1 x_2 ... x_nvar. And the varaible id need to be a consequent integer start
+ * with 0.
+ * ret_vec will contains len number of doubles. Where the first nvar elements is the gradient vector,
+ * and the rest of (nvar+1)*nvar/2 elements are the upper/lower plus the diagonal part of the Hessian matrix
+ * in row format.
+ * This algorithm is inefficient, because at each nodes, it didn't check the dependency of the independent
+ * variables up to the current node. (or it is hard to do so for this setup). Therefore, it computes a full
+ * for loops over each independent variable (ie. assume they are all dependent), for those independent
+ * variables that are not dependent at the current node, zero will be produced by computation.
+ * By default the forward mode hessian routing is disabled. To enable the forward hessian interface, the
+ * compiler marco FORWARD_ENABLED need to be set equal to 1 in auto_diff_types.h
+ *
+ * + Reverse Hessian*Vector Evaluation:
+ * Simple, building a tape in the forward pass, and a reverse pass will evaluate the Hessian*vector. The implemenation
+ * also discovery the repeated subexpression and use one piece of memory on the tape for the same subexpression. This
+ * allow efficient evaluation, because the repeated subexpression only evaluate once in the forward and reverse pass.
+ * This algorithm can be called n times to compute a full Hessian, where n equals the number of independent
+ * variables.
+ * */
+
+typedef boost::numeric::ublas::compressed_matrix<double,boost::numeric::ublas::column_major,0,std::vector<std::size_t>,std::vector<double> >  col_compress_matrix;
+typedef boost::numeric::ublas::matrix_row<col_compress_matrix > col_compress_matrix_row;
+typedef boost::numeric::ublas::matrix_column<col_compress_matrix  > col_compress_matrix_col;
+
+namespace AutoDiff{
+
+	//node creation methods
+	extern PNode* create_param_node(double value);
+	extern VNode* create_var_node(double v=NaN_Double);
+	extern OPNode* create_uary_op_node(OPCODE code, Node* left);
+	extern OPNode* create_binary_op_node(OPCODE code, Node* left,Node* right);
+
+	//single constraint version
+	extern double eval_function(Node* root);
+	extern unsigned int nzGrad(Node* root);
+	extern double grad_reverse(Node* root,vector<Node*>& nodes, vector<double>& grad);
+	extern unsigned int nzHess(EdgeSet&);
+	extern double hess_reverse(Node* root, vector<Node*>& nodes, vector<double>& dhess);
+
+	//multiple constraints version
+	extern unsigned int nzGrad(Node* root, boost::unordered_set<Node*>& vnodes);
+	extern double grad_reverse(Node* root, vector<Node*>& nodes, col_compress_matrix_row& rgrad);
+	extern unsigned int nzHess(EdgeSet&,boost::unordered_set<Node*>& set1, boost::unordered_set<Node*>& set2);
+	extern double hess_reverse(Node* root, vector<Node*>& nodes, col_compress_matrix_col& chess);
+
+#if FORWARD_ENDABLED
+	//forward methods
+	extern void hess_forward(Node* root, unsigned int nvar, double** hess_mat);
+#endif
+
+	//utiliy methods
+	extern void nonlinearEdges(Node* root, EdgeSet& edges);
+	extern unsigned int numTotalNodes(Node*);
+	extern string tree_expr(Node* root);
+	extern void print_tree(Node* root);
+	extern void autodiff_setup();
+	extern void autodiff_cleanup();
+};
+
+#endif /* AUTODIFF_H_ */