Adding upstream version 14.2.21.upstream/14.2.21upstream

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

47 files changed, 5503 insertions, 0 deletions

diff --git a/src/boost/libs/yap/example/CMakeLists.txt b/src/boost/libs/yap/example/CMakeLists.txt
new file mode 100644
index 00000000..434ff12b
--- /dev/null
+++ b/src/boost/libs/yap/example/CMakeLists.txt
@@ -0,0 +1,43 @@
+include_directories(${CMAKE_HOME_DIRECTORY})
+
+include(CTest)
+
+enable_testing()
+
+add_custom_target(run_examples COMMAND ${CMAKE_CTEST_COMMAND} -VV -C ${CMAKE_CFG_INTDIR})
+
+macro(add_sample name)
+    add_executable(${name} ${name}.cpp)
+    target_link_libraries(${name} yap)
+    add_test(${name} ${CMAKE_CURRENT_BINARY_DIR}/${name})
+    if (clang_on_linux)
+        target_link_libraries(${name} c++)
+    endif ()
+endmacro()
+
+add_sample(minimal)
+add_sample(hello_world)
+add_sample(hello_world_redux)
+add_sample(calc1)
+add_sample(calc2a)
+add_sample(calc2b)
+add_sample(calc3)
+add_sample(lazy_vector)
+add_sample(tarray)
+add_sample(vec3)
+add_sample(vector)
+add_sample(mixed)
+add_sample(map_assign)
+add_sample(future_group)
+add_sample(transform_terminals)
+add_sample(pipable_algorithms)
+if (constexpr_if_define STREQUAL "-DBOOST_NO_CONSTEXPR_IF=0")
+    add_sample(let)
+    add_sample(self_evaluation)
+endif ()
+
+add_executable(autodiff autodiff_example.cpp)
+target_link_libraries(autodiff yap boost autodiff_library)
+if (clang_on_linux)
+    target_link_libraries(autodiff c++)
+endif ()
diff --git a/src/boost/libs/yap/example/autodiff_example.cpp b/src/boost/libs/yap/example/autodiff_example.cpp
new file mode 100644
index 00000000..597a4094
--- /dev/null
+++ b/src/boost/libs/yap/example/autodiff_example.cpp
@@ -0,0 +1,852 @@
+// 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()
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_ */
diff --git a/src/boost/libs/yap/example/calc1.cpp b/src/boost/libs/yap/example/calc1.cpp
new file mode 100644
index 00000000..20c869ba
--- /dev/null
+++ b/src/boost/libs/yap/example/calc1.cpp
@@ -0,0 +1,27 @@
+// 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)
+//[ calc1
+#include <boost/yap/expression.hpp>
+
+#include <iostream>
+
+
+int main ()
+{
+    using namespace boost::yap::literals;
+
+    // Displays "5"
+    std::cout << evaluate( 1_p + 2.0, 3.0 ) << std::endl;
+
+    // Displays "6"
+    std::cout << evaluate( 1_p * 2_p, 3.0, 2.0 ) << std::endl;
+
+    // Displays "0.5"
+    std::cout << evaluate( (1_p - 2_p) / 2_p, 3.0, 2.0 ) << std::endl;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/calc2a.cpp b/src/boost/libs/yap/example/calc2a.cpp
new file mode 100644
index 00000000..8d28e5bf
--- /dev/null
+++ b/src/boost/libs/yap/example/calc2a.cpp
@@ -0,0 +1,45 @@
+// 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)
+//[ calc2a
+#include <boost/yap/expression.hpp>
+
+#include <iostream>
+
+
+int main ()
+{
+    using namespace boost::yap::literals;
+
+    auto expr_1 = 1_p + 2.0;
+
+    auto expr_1_fn = [expr_1](auto &&... args) {
+        return evaluate(expr_1, args...);
+    };
+
+    auto expr_2 = 1_p * 2_p;
+
+    auto expr_2_fn = [expr_2](auto &&... args) {
+        return evaluate(expr_2, args...);
+    };
+
+    auto expr_3 = (1_p - 2_p) / 2_p;
+
+    auto expr_3_fn = [expr_3](auto &&... args) {
+        return evaluate(expr_3, args...);
+    };
+
+    // Displays "5"
+    std::cout << expr_1_fn(3.0) << std::endl;
+
+    // Displays "6"
+    std::cout << expr_2_fn(3.0, 2.0) << std::endl;
+
+    // Displays "0.5"
+    std::cout << expr_3_fn(3.0, 2.0) << std::endl;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/calc2b.cpp b/src/boost/libs/yap/example/calc2b.cpp
new file mode 100644
index 00000000..cad5b913
--- /dev/null
+++ b/src/boost/libs/yap/example/calc2b.cpp
@@ -0,0 +1,27 @@
+// 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)
+//[ calc2b
+#include <boost/yap/expression.hpp>
+
+#include <iostream>
+
+
+int main ()
+{
+    using namespace boost::yap::literals;
+
+    // Displays "5"
+    std::cout << make_expression_function(1_p + 2.0)(3.0) << std::endl;
+
+    // Displays "6"
+    std::cout << make_expression_function(1_p * 2_p)(3.0, 2.0) << std::endl;
+
+    // Displays "0.5"
+    std::cout << make_expression_function((1_p - 2_p) / 2_p)(3.0, 2.0) << std::endl;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/calc3.cpp b/src/boost/libs/yap/example/calc3.cpp
new file mode 100644
index 00000000..7f0f55b7
--- /dev/null
+++ b/src/boost/libs/yap/example/calc3.cpp
@@ -0,0 +1,100 @@
+// 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)
+//[ calc3
+#include <boost/yap/expression.hpp>
+
+#include <boost/hana/maximum.hpp>
+
+#include <iostream>
+
+
+// Look! A transform!  This one transforms the expression tree into the arity
+// of the expression, based on its placeholders.
+//[ calc3_get_arity_xform
+struct get_arity
+{
+    // Base case 1: Match a placeholder terminal, and return its arity as the
+    // result.
+    template <long long I>
+    boost::hana::llong<I> operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
+                                      boost::yap::placeholder<I>)
+    { return boost::hana::llong_c<I>; }
+
+    // Base case 2: Match any other terminal.  Return 0; non-placeholders do
+    // not contribute to arity.
+    template <typename T>
+    auto operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>, T &&)
+    {
+        using namespace boost::hana::literals;
+        return 0_c;
+    }
+
+    // Recursive case: Match any expression not covered above, and return the
+    // maximum of its children's arities.
+    template <boost::yap::expr_kind Kind, typename... Arg>
+    auto operator() (boost::yap::expr_tag<Kind>, Arg &&... arg)
+    {
+        return boost::hana::maximum(
+            boost::hana::make_tuple(
+                boost::yap::transform(
+                    boost::yap::as_expr(std::forward<Arg>(arg)),
+                    get_arity{}
+                )...
+            )
+        );
+    }
+};
+//]
+
+int main ()
+{
+    using namespace boost::yap::literals;
+
+    // These lambdas wrap our expressions as callables, and allow us to check
+    // the arity of each as we call it.
+
+    auto expr_1 = 1_p + 2.0;
+
+    auto expr_1_fn = [expr_1](auto &&... args) {
+        auto const arity = boost::yap::transform(expr_1, get_arity{});
+        static_assert(arity.value == sizeof...(args), "Called with wrong number of args.");
+        return evaluate(expr_1, args...);
+    };
+
+    auto expr_2 = 1_p * 2_p;
+
+    auto expr_2_fn = [expr_2](auto &&... args) {
+        auto const arity = boost::yap::transform(expr_2, get_arity{});
+        static_assert(arity.value == sizeof...(args), "Called with wrong number of args.");
+        return evaluate(expr_2, args...);
+    };
+
+    auto expr_3 = (1_p - 2_p) / 2_p;
+
+    auto expr_3_fn = [expr_3](auto &&... args) {
+        auto const arity = boost::yap::transform(expr_3, get_arity{});
+        static_assert(arity.value == sizeof...(args), "Called with wrong number of args.");
+        return evaluate(expr_3, args...);
+    };
+
+    // Displays "5"
+    std::cout << expr_1_fn(3.0) << std::endl;
+
+    // Displays "6"
+    std::cout << expr_2_fn(3.0, 2.0) << std::endl;
+
+    // Displays "0.5"
+    std::cout << expr_3_fn(3.0, 2.0) << std::endl;
+
+    // Static-asserts with "Called with wrong number of args."
+    //std::cout << expr_3_fn(3.0) << std::endl;
+
+    // Static-asserts with "Called with wrong number of args."
+    //std::cout << expr_3_fn(3.0, 2.0, 1.0) << std::endl;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/future_group.cpp b/src/boost/libs/yap/example/future_group.cpp
new file mode 100644
index 00000000..15c980f0
--- /dev/null
+++ b/src/boost/libs/yap/example/future_group.cpp
@@ -0,0 +1,129 @@
+// 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)
+//[ future_group
+#include <boost/yap/algorithm.hpp>
+
+#include <boost/hana/concat.hpp>
+
+
+// A custom expression template for future groups.  It supports operators ||
+// and &&.
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct future_expr
+{
+    static boost::yap::expr_kind const kind = Kind;
+
+    future_expr (Tuple && tuple) :
+        elements (std::forward<Tuple &&>(tuple))
+    {}
+
+    Tuple elements;
+
+    // Returns the transformed/flattened expression.
+    auto get () const;
+};
+
+BOOST_YAP_USER_BINARY_OPERATOR(logical_or, future_expr, future_expr)
+BOOST_YAP_USER_BINARY_OPERATOR(logical_and, future_expr, future_expr)
+
+// A special-cased future terminal that matches the semantics from the
+// original Proto example.
+template <typename T>
+struct future :
+    future_expr<boost::yap::expr_kind::terminal, boost::hana::tuple<T>>
+{
+    future (T const & t = T()) :
+        future_expr<boost::yap::expr_kind::terminal, boost::hana::tuple<T>> (boost::hana::tuple<T>{t})
+    {}
+
+    T get () const
+    { return boost::yap::value(*this); }
+};
+
+template <typename T>
+using remove_cv_ref_t = std::remove_cv_t<std::remove_reference_t<T>>;
+
+// A transform that flattens future expressions into a tuple.
+struct future_transform
+{
+    // Transform a terminal into its contained tuple.
+    template <typename T>
+    auto operator() (
+        future_expr<
+            boost::yap::expr_kind::terminal,
+            boost::hana::tuple<T>
+        > const & term
+    ) {
+        return term.elements;
+    }
+
+//[ expr_xform
+    // Transform left || right -> transform(left).
+    template <typename T, typename U>
+    auto operator() (
+        future_expr<
+            boost::yap::expr_kind::logical_or,
+            boost::hana::tuple<T, U>
+        > const & or_expr
+    ) {
+        // Recursively transform the left side, and return the result.
+        // Without the recursion, we might return a terminal expression here
+        // insead of a tuple.
+        return boost::yap::transform(boost::yap::left(or_expr), *this);
+    }
+//]
+
+    // Transform left && right -> concat(transform(left), transform(right)).
+    template <typename T, typename U>
+    auto operator() (
+        future_expr<
+            boost::yap::expr_kind::logical_and,
+            boost::hana::tuple<T, U>
+        > const & and_expr
+    ) {
+        // Recursively transform each side, then combine the resulting tuples
+        // into a single tuple result.
+        return boost::hana::concat(
+            boost::yap::transform(boost::yap::left(and_expr), *this),
+            boost::yap::transform(boost::yap::right(and_expr), *this)
+        );
+    }
+};
+
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+auto future_expr<Kind, Tuple>::get () const
+{ return boost::yap::transform(*this, future_transform{}); }
+
+
+// TEST CASES
+struct A {};
+struct B {};
+struct C {};
+
+// Called "vector" just so the code in main() will match the original Proto
+// example.
+template <typename ...T>
+using vector = boost::hana::tuple<T...>;
+
+int main()
+{
+    future<A> a;
+    future<B> b;
+    future<C> c;
+    future<vector<A,B> > ab;
+
+    // Verify that various future groups have the
+    // correct return types.
+    A                       t0 = a.get();
+    vector<A, B, C>         t1 = (a && b && c).get();
+    vector<A, C>            t2 = ((a || a) && c).get();
+    vector<A, B, C>         t3 = ((a && b || a && b) && c).get();
+    vector<vector<A, B>, C> t4 = ((ab || ab) && c).get();
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/hello_world.cpp b/src/boost/libs/yap/example/hello_world.cpp
new file mode 100644
index 00000000..f509078e
--- /dev/null
+++ b/src/boost/libs/yap/example/hello_world.cpp
@@ -0,0 +1,19 @@
+// 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)
+//[hello_world
+
+#include <boost/yap/expression.hpp>
+
+#include <iostream>
+
+
+int main ()
+{
+    evaluate(boost::yap::make_terminal(std::cout) << "Hello" << ',' << " world!\n");
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/hello_world_redux.cpp b/src/boost/libs/yap/example/hello_world_redux.cpp
new file mode 100644
index 00000000..0bb7eac0
--- /dev/null
+++ b/src/boost/libs/yap/example/hello_world_redux.cpp
@@ -0,0 +1,32 @@
+// 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)
+//[hello_world_redux
+#include <boost/yap/algorithm.hpp>
+
+#include <iostream>
+
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct stream_expr
+{
+    static const boost::yap::expr_kind kind = Kind;
+
+    Tuple elements;
+
+    template <typename T>
+    decltype(auto) operator<< (T && x)
+    { return boost::yap::value(*this) << std::forward<T &&>(x); }
+};
+
+
+int main ()
+{
+    auto cout = boost::yap::make_terminal<stream_expr>(std::cout);
+    cout << "Hello" << ',' << " world!\n";
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/lazy_vector.cpp b/src/boost/libs/yap/example/lazy_vector.cpp
new file mode 100644
index 00000000..dd43ed84
--- /dev/null
+++ b/src/boost/libs/yap/example/lazy_vector.cpp
@@ -0,0 +1,112 @@
+// 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)
+//[ lazy_vector
+#include <boost/yap/expression.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <iostream>
+#include <vector>
+
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct lazy_vector_expr;
+
+
+// This transform turns a terminal of std::vector<double> into a terminal
+// containing the nth double in that vector.  Think of it as turning our
+// expression of vectors into an expression of scalars.
+struct take_nth
+{
+    boost::yap::terminal<lazy_vector_expr, double>
+    operator() (boost::yap::terminal<lazy_vector_expr, std::vector<double>> const & expr);
+
+    std::size_t n;
+};
+
+// A custom expression template that defines lazy + and - operators that
+// produce expressions, and an eager [] operator that returns the nth element
+// of the expression.
+//[ lazy_vector_decl
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct lazy_vector_expr
+{
+    static const boost::yap::expr_kind kind = Kind;
+
+    Tuple elements;
+
+    // Note that this does not return an expression; it is greedily evaluated.
+    auto operator[] (std::size_t n) const;
+};
+
+BOOST_YAP_USER_BINARY_OPERATOR(plus, lazy_vector_expr, lazy_vector_expr)
+BOOST_YAP_USER_BINARY_OPERATOR(minus, lazy_vector_expr, lazy_vector_expr)
+//]
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+auto lazy_vector_expr<Kind, Tuple>::operator[] (std::size_t n) const
+{ return boost::yap::evaluate(boost::yap::transform(*this, take_nth{n})); }
+
+boost::yap::terminal<lazy_vector_expr, double>
+take_nth::operator() (boost::yap::terminal<lazy_vector_expr, std::vector<double>> const & expr)
+{
+    double x = boost::yap::value(expr)[n];
+    // This move is something of a hack; we're forcing Yap to take a copy of x
+    // by using std::move().  The move indicates that the terminal should keep
+    // the value of x (since, being an rvalue, it may be a temporary), rather
+    // than a reference to x.  See the "How Expression Operands Are Treated"
+    // section of the tutorial for details.
+    return boost::yap::make_terminal<lazy_vector_expr, double>(std::move(x));
+}
+
+// In order to define the += operator with the semantics we want, it's
+// convenient to derive a terminal type from a terminal instantiation of
+// lazy_vector_expr.  Note that we could have written a template
+// specialization here instead -- either one would work.  That would of course
+// have required more typing.
+struct lazy_vector :
+    lazy_vector_expr<
+        boost::yap::expr_kind::terminal,
+        boost::hana::tuple<std::vector<double>>
+    >
+{
+    lazy_vector () {}
+
+    explicit lazy_vector (std::vector<double> && vec)
+    { elements = boost::hana::tuple<std::vector<double>>(std::move(vec)); }
+
+    template <boost::yap::expr_kind Kind, typename Tuple>
+    lazy_vector & operator+= (lazy_vector_expr<Kind, Tuple> const & rhs)
+    {
+        std::vector<double> & this_vec = boost::yap::value(*this);
+        for (int i = 0, size = (int)this_vec.size(); i < size; ++i) {
+            this_vec[i] += rhs[i];
+        }
+        return *this;
+    }
+};
+
+int main ()
+{
+    lazy_vector v1{std::vector<double>(4, 1.0)};
+    lazy_vector v2{std::vector<double>(4, 2.0)};
+    lazy_vector v3{std::vector<double>(4, 3.0)};
+
+    double d1 = (v2 + v3)[2];
+    std::cout << d1 << "\n";
+
+    v1 += v2 - v3;
+    std::cout << '{' << v1[0] << ',' << v1[1]
+              << ',' << v1[2] << ',' << v1[3] << '}' << "\n";
+
+    // This expression is disallowed because it does not conform to the
+    // implicit grammar.  operator+= is only defined on terminals, not
+    // arbitrary expressions.
+    // (v2 + v3) += v1;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/let.cpp b/src/boost/libs/yap/example/let.cpp
new file mode 100644
index 00000000..0cac1975
--- /dev/null
+++ b/src/boost/libs/yap/example/let.cpp
@@ -0,0 +1,177 @@
+// 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)
+//[ let
+#include <boost/yap/yap.hpp>
+
+#include <boost/hana/map.hpp>
+#include <boost/hana/at_key.hpp>
+#include <boost/hana/contains.hpp>
+#include <boost/hana/keys.hpp>
+
+#include <vector>
+#include <iostream>
+
+
+// Here, we introduce special let-placeholders, so we can use them along side
+// the normal YAP placeholders without getting them confused.
+template<long long I>
+struct let_placeholder : boost::hana::llong<I>
+{
+};
+
+// Replaces each let-terminal with the expression with which it was
+// initialized in let().  So in 'let(_a = foo)[ _a + 1 ]', this transform will
+// be used on '_a + 1' to replace '_a' with 'foo'.  The map_ member holds the
+// mapping of let-placeholders to their initializers.
+template<typename ExprMap>
+struct let_terminal_transform
+{
+    // This matches only let-placeholders.  For each one matched, we look up
+    // its initializer in map_ and return it.
+    template<long long I>
+    auto operator()(
+        boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
+        let_placeholder<I> i)
+    {
+        // If we have an entry in map_ for this placeholder, return the value
+        // of the entry.  Otherwise, pass i through as a terminal.
+        if constexpr (boost::hana::contains(
+                          decltype(boost::hana::keys(map_))(),
+                          boost::hana::llong_c<I>)) {
+            return map_[boost::hana::llong_c<I>];
+        } else {
+            return boost::yap::make_terminal(i);
+        }
+    }
+
+    ExprMap map_;
+};
+
+// As you can see below, let() is an eager function; this template is used for
+// its return values.  It contains the mapping from let-placeholders to
+// initializer expressions used to transform the expression inside '[]' after
+// a let()'.  It also has an operator[]() member function that takes the
+// expression inside '[]' and returns a version of it with the
+// let-placeholders replaced.
+template<typename ExprMap>
+struct let_result
+{
+    template<typename Expr>
+    auto operator[](Expr && expr)
+    {
+        return boost::yap::transform(
+            std::forward<Expr>(expr), let_terminal_transform<ExprMap>{map_});
+    }
+
+    ExprMap map_;
+};
+
+// Processes the expressions passed to let() one at a time, adding each one to
+// a Hana map of hana::llong<>s to YAP expressions.
+template<typename Map, typename Expr, typename... Exprs>
+auto let_impl(Map && map, Expr && expr, Exprs &&... exprs)
+{
+    static_assert(
+        Expr::kind == boost::yap::expr_kind::assign,
+        "Expressions passed to let() must be of the form placeholder = Expression");
+    if constexpr (sizeof...(Exprs) == 0) {
+        using I = typename std::remove_reference<decltype(
+            boost::yap::value(boost::yap::left(expr)))>::type;
+        auto const i = boost::hana::llong_c<I::value>;
+        using map_t = decltype(boost::hana::insert(
+            map, boost::hana::make_pair(i, boost::yap::right(expr))));
+        return let_result<map_t>{boost::hana::insert(
+            map, boost::hana::make_pair(i, boost::yap::right(expr)))};
+    } else {
+        using I = typename std::remove_reference<decltype(
+            boost::yap::value(boost::yap::left(expr)))>::type;
+        auto const i = boost::hana::llong_c<I::value>;
+        return let_impl(
+            boost::hana::insert(
+                map, boost::hana::make_pair(i, boost::yap::right(expr))),
+            std::forward<Exprs>(exprs)...);
+    }
+}
+
+// Takes N > 0 expressions of the form 'placeholder = expr', and returns an
+// object with an overloaded operator[]().
+template<typename Expr, typename... Exprs>
+auto let(Expr && expr, Exprs &&... exprs)
+{
+    return let_impl(
+        boost::hana::make_map(),
+        std::forward<Expr>(expr),
+        std::forward<Exprs>(exprs)...);
+}
+
+int main()
+{
+    // Some handy terminals -- the _a and _b let-placeholders and std::cout as
+    // a YAP terminal.
+    boost::yap::expression<
+        boost::yap::expr_kind::terminal,
+        boost::hana::tuple<let_placeholder<0>>> const _a;
+    boost::yap::expression<
+        boost::yap::expr_kind::terminal,
+        boost::hana::tuple<let_placeholder<1>>> const _b;
+    auto const cout = boost::yap::make_terminal(std::cout);
+
+    using namespace boost::yap::literals;
+
+    {
+        auto expr = let(_a = 2)[_a + 1];
+        assert(boost::yap::evaluate(expr) == 3);
+    }
+
+    {
+        auto expr = let(_a = 123, _b = 456)[_a + _b];
+        assert(boost::yap::evaluate(expr) == 123 + 456);
+    }
+
+    // This prints out "0 0", because 'i' is passed as an lvalue, so its
+    // decrement is visible outside the let expression.
+    {
+        int i = 1;
+
+        boost::yap::evaluate(let(_a = 1_p)[cout << --_a << ' '], i);
+
+        std::cout << i << std::endl;
+    }
+
+    // Prints "Hello, World" due to let()'s scoping rules.
+    {
+        boost::yap::evaluate(
+            let(_a = 1_p, _b = 2_p)
+            [
+                // _a here is an int: 1
+
+                let(_a = 3_p) // hides the outer _a
+                [
+                    cout << _a << _b // prints "Hello, World"
+                ]
+            ],
+            1, " World", "Hello,"
+        );
+    }
+
+    std::cout << "\n";
+
+    // Due to the macro-substitution style that this example uses, this prints
+    // "3132".  Phoenix's let() prints "312", because it only evaluates '1_p
+    // << 3' once.
+    {
+        boost::yap::evaluate(
+            let(_a = 1_p << 3)
+            [
+                _a << "1", _a << "2"
+            ],
+            std::cout
+        );
+    }
+
+    std::cout << "\n";
+}
+//]
diff --git a/src/boost/libs/yap/example/map_assign.cpp b/src/boost/libs/yap/example/map_assign.cpp
new file mode 100644
index 00000000..d8b304be
--- /dev/null
+++ b/src/boost/libs/yap/example/map_assign.cpp
@@ -0,0 +1,93 @@
+// 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)
+//[ map_assign
+#include <boost/yap/algorithm.hpp>
+
+#include <map>
+#include <iostream>
+
+
+// This transform applies all the call-subexpressions in a map_list_of
+// expression (a nested chain of call operations) as a side effect; the
+// expression returned by the transform is ignored.
+template <typename Key, typename Value, typename Allocator>
+struct map_list_of_transform
+{
+    template <typename Fn, typename Key2, typename Value2>
+    auto operator() (boost::yap::expr_tag<boost::yap::expr_kind::call>,
+                     Fn const & fn, Key2 && key, Value2 && value)
+    {
+        // Recurse into the function subexpression.  Remember, transform()
+        // walks the nodes in an expression tree looking for matches.  Once it
+        // finds a match, it is finished with that matching subtree.  So
+        // without this recursive call, only the top-level call expression is
+        // matched by transform().
+        boost::yap::transform(
+            boost::yap::as_expr<boost::yap::minimal_expr>(fn), *this);
+        map.emplace(
+            std::forward<Key2 &&>(key),
+            std::forward<Value2 &&>(value)
+        );
+        // All we care about are the side effects of this transform, so we can
+        // return any old thing here.
+        return 0;
+    }
+
+    std::map<Key, Value, Allocator> & map;
+};
+
+
+// A custom expression template type for map_list_of expressions.  We only
+// need support for the call operator and an implicit conversion to a
+// std::map.
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct map_list_of_expr
+{
+    static boost::yap::expr_kind const kind = Kind;
+
+    Tuple elements;
+
+    template <typename Key, typename Value, typename Allocator>
+    operator std::map<Key, Value, Allocator> () const
+    {
+        std::map<Key, Value, Allocator> retval;
+        map_list_of_transform<Key, Value, Allocator> transform{retval};
+        boost::yap::transform(*this, transform);
+        return retval;
+    }
+
+    BOOST_YAP_USER_CALL_OPERATOR_N(::map_list_of_expr, 2)
+};
+
+// A tag type for creating the map_list_of function terminal.
+struct map_list_of_tag {};
+
+auto map_list_of = boost::yap::make_terminal<map_list_of_expr>(map_list_of_tag{});
+
+
+int main()
+{
+    // Initialize a map:
+    std::map<std::string, int> op =
+        map_list_of
+            ("<", 1)
+            ("<=",2)
+            (">", 3)
+            (">=",4)
+            ("=", 5)
+            ("<>",6)
+        ;
+
+    std::cout << "\"<\"  --> " << op["<"] << std::endl;
+    std::cout << "\"<=\" --> " << op["<="] << std::endl;
+    std::cout << "\">\"  --> " << op[">"] << std::endl;
+    std::cout << "\">=\" --> " << op[">="] << std::endl;
+    std::cout << "\"=\"  --> " << op["="] << std::endl;
+    std::cout << "\"<>\" --> " << op["<>"] << std::endl;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/minimal.cpp b/src/boost/libs/yap/example/minimal.cpp
new file mode 100644
index 00000000..d79e372f
--- /dev/null
+++ b/src/boost/libs/yap/example/minimal.cpp
@@ -0,0 +1,42 @@
+// 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 <boost/yap/algorithm.hpp>
+
+#include <array>
+#include <iostream>
+
+
+//[ minimal_template
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct minimal_expr
+{
+    static const boost::yap::expr_kind kind = Kind;
+
+    Tuple elements;
+};
+//]
+
+
+int main()
+{
+//[ minimal_template_manual_construction
+    auto left = boost::yap::make_terminal<minimal_expr>(1);
+    auto right = boost::yap::make_terminal<minimal_expr>(41);
+
+    auto expr = boost::yap::make_expression<
+        minimal_expr,
+        boost::yap::expr_kind::plus
+    >(left, right);
+//]
+
+//[ minimal_template_evaluation
+    auto result = boost::yap::evaluate(expr);
+
+    std::cout << result << "\n"; // prints "42"
+//]
+
+    return 0;
+}
diff --git a/src/boost/libs/yap/example/mixed.cpp b/src/boost/libs/yap/example/mixed.cpp
new file mode 100644
index 00000000..e7dee54c
--- /dev/null
+++ b/src/boost/libs/yap/example/mixed.cpp
@@ -0,0 +1,213 @@
+// 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)
+//[ mixed
+#include <boost/yap/yap.hpp>
+
+#include <complex>
+#include <list>
+#include <vector>
+#include <iostream>
+
+
+// This wrapper makes the pattern matching in transforms below (like deref and
+// incr) a lot easier to write.
+template <typename Iter>
+struct iter_wrapper
+{
+    Iter it;
+};
+
+template <typename Iter>
+auto make_iter_wrapper (Iter it)
+{ return iter_wrapper<Iter>{it}; }
+
+
+// A container -> wrapped-begin transform.
+struct begin
+{
+    template <typename Cont>
+    auto operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
+                     Cont const & cont)
+        -> decltype(boost::yap::make_terminal(make_iter_wrapper(cont.begin())))
+    { return boost::yap::make_terminal(make_iter_wrapper(cont.begin())); }
+};
+
+// A wrapped-iterator -> dereferenced value transform.
+struct deref
+{
+    template <typename Iter>
+    auto operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
+                     iter_wrapper<Iter> wrapper)
+        -> decltype(boost::yap::make_terminal(*wrapper.it))
+    { return boost::yap::make_terminal(*wrapper.it); }
+};
+
+// A wrapped-iterator increment transform, using side effects.
+struct incr
+{
+    template <typename Iter>
+    auto operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
+                     iter_wrapper<Iter> & wrapper)
+        -> decltype(boost::yap::make_terminal(wrapper.it))
+    {
+        ++wrapper.it;
+        // Since this transform is valuable for its side effects, and thus the
+        // result of the transform is ignored, we could return anything here.
+        return boost::yap::make_terminal(wrapper.it);
+    }
+};
+
+
+// The implementation of elementwise evaluation of expressions of sequences;
+// all the later operations use this one.
+template <
+    template <class, class> class Cont,
+    typename T,
+    typename A,
+    typename Expr,
+    typename Op
+>
+Cont<T, A> & op_assign (Cont<T, A> & cont, Expr const & e, Op && op)
+{
+    decltype(auto) expr = boost::yap::as_expr(e);
+    // Transform the expression of sequences into an expression of
+    // begin-iterators.
+    auto expr2 = boost::yap::transform(boost::yap::as_expr(expr), begin{});
+    for (auto && x : cont) {
+        // Transform the expression of iterators into an expression of
+        // pointed-to-values, evaluate the resulting expression, and call op()
+        // with the result of the evaluation.
+        op(x, boost::yap::evaluate(boost::yap::transform(expr2, deref{})));
+        // Transform the expression of iterators into an ignored value; as a
+        // side effect, increment the iterators in the expression.
+        boost::yap::transform(expr2, incr{});
+    }
+    return cont;
+}
+
+template <
+    template <class, class> class Cont,
+    typename T,
+    typename A,
+    typename Expr
+>
+Cont<T, A> & assign (Cont<T, A> & cont, Expr const & expr)
+{
+    return op_assign(cont, expr, [](auto & cont_value, auto && expr_value) {
+        cont_value = std::forward<decltype(expr_value)>(expr_value);
+    });
+}
+
+template <
+    template <class, class> class Cont,
+    typename T,
+    typename A,
+    typename Expr
+>
+Cont<T, A> & operator+= (Cont<T, A> & cont, Expr const & expr)
+{
+    return op_assign(cont, expr, [](auto & cont_value, auto && expr_value) {
+        cont_value += std::forward<decltype(expr_value)>(expr_value);
+    });
+}
+
+template <
+    template <class, class> class Cont,
+    typename T,
+    typename A,
+    typename Expr
+>
+Cont<T, A> & operator-= (Cont<T, A> & cont, Expr const & expr)
+{
+    return op_assign(cont, expr, [](auto & cont_value, auto && expr_value) {
+        cont_value -= std::forward<decltype(expr_value)>(expr_value);
+    });
+}
+
+// A type trait that identifies std::vectors and std::lists.
+template <typename T>
+struct is_mixed : std::false_type {};
+
+template <typename T, typename A>
+struct is_mixed<std::vector<T, A>> : std::true_type {};
+
+template <typename T, typename A>
+struct is_mixed<std::list<T, A>> : std::true_type {};
+
+// Define expression-producing operators over std::vectors and std::lists.
+BOOST_YAP_USER_UDT_UNARY_OPERATOR(negate, boost::yap::expression, is_mixed); // -
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(multiplies, boost::yap::expression, is_mixed); // *
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(divides, boost::yap::expression, is_mixed); // /
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(modulus, boost::yap::expression, is_mixed); // %
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(plus, boost::yap::expression, is_mixed); // +
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(minus, boost::yap::expression, is_mixed); // -
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(less, boost::yap::expression, is_mixed); // <
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(greater, boost::yap::expression, is_mixed); // >
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(less_equal, boost::yap::expression, is_mixed); // <=
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(greater_equal, boost::yap::expression, is_mixed); // >=
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(equal_to, boost::yap::expression, is_mixed); // ==
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(not_equal_to, boost::yap::expression, is_mixed); // !=
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(logical_or, boost::yap::expression, is_mixed); // ||
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(logical_and, boost::yap::expression, is_mixed); // &&
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(bitwise_and, boost::yap::expression, is_mixed); // &
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(bitwise_or, boost::yap::expression, is_mixed); // |
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(bitwise_xor, boost::yap::expression, is_mixed); // ^
+
+// Define a type that can resolve to any overload of std::sin().
+struct sin_t
+{
+    template<typename T>
+    T operator()(T x)
+    {
+        return std::sin(x);
+    }
+};
+
+int main()
+{
+    int n = 10;
+    std::vector<int> a,b,c,d;
+    std::list<double> e;
+    std::list<std::complex<double>> f;
+
+    int i;
+    for(i = 0;i < n; ++i)
+    {
+        a.push_back(i);
+        b.push_back(2*i);
+        c.push_back(3*i);
+        d.push_back(i);
+        e.push_back(0.0);
+        f.push_back(std::complex<double>(1.0, 1.0));
+    }
+
+    assign(b, 2);
+    assign(d, a + b * c);
+    a += if_else(d < 30, b, c);
+
+    assign(e, c);
+    e += e - 4 / (c + 1);
+
+    auto sin = boost::yap::make_terminal(sin_t{});
+    f -= sin(0.1 * e * std::complex<double>(0.2, 1.2));
+
+    std::list<double>::const_iterator ei = e.begin();
+    std::list<std::complex<double>>::const_iterator fi = f.begin();
+    for (i = 0; i < n; ++i)
+    {
+        std::cout
+            << "a(" << i << ") = " << a[i]
+            << " b(" << i << ") = " << b[i]
+            << " c(" << i << ") = " << c[i]
+            << " d(" << i << ") = " << d[i]
+            << " e(" << i << ") = " << *ei++
+            << " f(" << i << ") = " << *fi++
+            << std::endl;
+    }
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/pipable_algorithms.cpp b/src/boost/libs/yap/example/pipable_algorithms.cpp
new file mode 100644
index 00000000..07c71f54
--- /dev/null
+++ b/src/boost/libs/yap/example/pipable_algorithms.cpp
@@ -0,0 +1,147 @@
+// 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)
+//[ pipable_algorithms
+#include <boost/yap/algorithm.hpp>
+
+#include <algorithm>
+#include <vector>
+
+
+//[ pipable_algorithms_and_simple_range
+// An enum to represent all the standard algorithms we want to adapt.  In this
+// example, we only care about std::sort() and std::unique().
+enum class algorithm_t { sort, unique };
+
+// Just about the simplest range template you could construct.  Nothing fancy.
+template<typename Iter>
+struct simple_range
+{
+    Iter begin() const { return first_; }
+    Iter end() const { return last_; }
+
+    Iter first_;
+    Iter last_;
+};
+
+// This simply calls the standard algorithm that corresponds to "a".  This
+// certainly won't work for all the algorithms, but it works for many of them
+// that just take a begin/end pair.
+template<typename Range>
+auto call_algorithm(algorithm_t a, Range & r)
+{
+    simple_range<decltype(r.begin())> retval{r.begin(), r.end()};
+    if (a == algorithm_t::sort) {
+        std::sort(r.begin(), r.end());
+    } else if (a == algorithm_t::unique) {
+        retval.last_ = std::unique(r.begin(), r.end());
+    }
+    return retval;
+}
+//]
+
+// This is the transform that evaluates our piped expressions.  It returns a
+// simple_range<>, not a Yap expression.
+struct algorithm_eval
+{
+    // A pipe should always have a Yap expression on the left and an
+    // algorithm_t terminal on the right.
+    template<typename LExpr>
+    auto operator()(
+        boost::yap::expr_tag<boost::yap::expr_kind::bitwise_or>,
+        LExpr && left,
+        algorithm_t right)
+    {
+        // Recursively evaluate the left ...
+        auto const left_result =
+            boost::yap::transform(std::forward<LExpr>(left), *this);
+        // ... and use the result to call the function on the right.
+        return call_algorithm(right, left_result);
+    }
+
+    // A call operation is evaluated directly.  Note that the range parameter
+    // is taken as an lvalue reference, to prevent binding to a temporary and
+    // taking dangling references to its begin and end.  We let the compiler
+    // deduce whether the lvalue reference is const.
+//[ tag_xform
+    template<typename Range>
+    auto operator()(
+        boost::yap::expr_tag<boost::yap::expr_kind::call>,
+        algorithm_t a,
+        Range & r)
+    {
+        return call_algorithm(a, r);
+    }
+//]
+};
+
+// This is the expression template we use for the general case of a pipable
+// algorithm expression.  Terminals are handled separately.
+template<boost::yap::expr_kind Kind, typename Tuple>
+struct algorithm_expr
+{
+    static boost::yap::expr_kind const kind = Kind;
+
+    Tuple elements;
+
+    // This is how we get the nice initializer semantics we see in the example
+    // usage below.  This is a bit limited though, because we always create a
+    // temporary.  It might therefore be better just to create algorithm_expr
+    // expressions, call yap::evaluate(), and then use the sequence containers
+    // assign() member function to do the actual assignment.
+    template<typename Assignee>
+    operator Assignee() const
+    {
+        // Exercise left for the reader: static_assert() that Assignee is some
+        // sort of container type.
+        auto const range = boost::yap::transform(*this, algorithm_eval{});
+        return Assignee(range.begin(), range.end());
+    }
+};
+
+
+// This is a bit loose, because it allows us to write "sort(v) | unique(u)" or
+// similar.  It works fine for this example, but in production code you may
+// want to write out the functions generated by this macro, and add SFINAE or
+// concepts constraints on the right template parameter.
+BOOST_YAP_USER_BINARY_OPERATOR(bitwise_or, algorithm_expr, algorithm_expr)
+
+// It's useful to specially handle terminals, because we want a different set
+// of operations to apply to them.  We don't want "sort(x)(y)" to be
+// well-formed, for instance, or "sort | unique" either.
+struct algorithm
+{
+    static boost::yap::expr_kind const kind = boost::yap::expr_kind::terminal;
+
+    boost::hana::tuple<algorithm_t> elements;
+
+    BOOST_YAP_USER_CALL_OPERATOR_N(::algorithm_expr, 1)
+};
+
+// Here are ready-made Yap terminals, one for each algorithm enumerated in
+// algorithm_t.
+constexpr algorithm sort{{algorithm_t::sort}};
+constexpr algorithm unique{{algorithm_t::unique}};
+
+int main()
+{
+    {
+//[ typical_sort_unique_usage
+        std::vector<int> v1 = {0, 2, 2, 7, 1, 3, 8};
+        std::sort(v1.begin(), v1.end());
+        auto it = std::unique(v1.begin(), v1.end());
+        std::vector<int> const v2(v1.begin(), it);
+        assert(v2 == std::vector<int>({0, 1, 2, 3, 7, 8}));
+//]
+    }
+    {
+//[ pipable_sort_unique_usage
+        std::vector<int> v1 = {0, 2, 2, 7, 1, 3, 8};
+        std::vector<int> const v2 = sort(v1) | unique;
+        assert(v2 == std::vector<int>({0, 1, 2, 3, 7, 8}));
+//]
+    }
+}
+//]
diff --git a/src/boost/libs/yap/example/self_evaluation.cpp b/src/boost/libs/yap/example/self_evaluation.cpp
new file mode 100644
index 00000000..d5c7323d
--- /dev/null
+++ b/src/boost/libs/yap/example/self_evaluation.cpp
@@ -0,0 +1,237 @@
+// 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)
+//[ self_evaluation
+#include <boost/yap/expression.hpp>
+
+#include <boost/optional.hpp>
+#include <boost/hana/fold.hpp>
+#include <boost/hana/maximum.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <iostream>
+#include <vector>
+
+
+// A super-basic matrix type, and a few associated operations.
+struct matrix
+{
+    matrix() : values_(), rows_(0), cols_(0) {}
+
+    matrix(int rows, int cols) : values_(rows * cols), rows_(rows), cols_(cols)
+    {
+        assert(0 < rows);
+        assert(0 < cols);
+    }
+
+    int rows() const { return rows_; }
+    int cols() const { return cols_; }
+
+    double operator()(int r, int c) const
+    { return values_[r * cols_ + c]; }
+    double & operator()(int r, int c)
+    { return values_[r * cols_ + c]; }
+
+private:
+    std::vector<double> values_;
+    int rows_;
+    int cols_;
+};
+
+matrix operator*(matrix const & lhs, double x)
+{
+    matrix retval = lhs;
+    for (int i = 0; i < retval.rows(); ++i) {
+        for (int j = 0; j < retval.cols(); ++j) {
+            retval(i, j) *= x;
+        }
+    }
+    return retval;
+}
+matrix operator*(double x, matrix const & lhs) { return lhs * x; }
+
+matrix operator+(matrix const & lhs, matrix const & rhs)
+{
+    assert(lhs.rows() == rhs.rows());
+    assert(lhs.cols() == rhs.cols());
+    matrix retval = lhs;
+    for (int i = 0; i < retval.rows(); ++i) {
+        for (int j = 0; j < retval.cols(); ++j) {
+            retval(i, j) += rhs(i, j);
+        }
+    }
+    return retval;
+}
+
+// daxpy() means Double-precision AX Plus Y.  This crazy name comes from BLAS.
+// It is more efficient than a naive implementation, because it does not
+// create temporaries.  The covnention of using Y as an out-parameter comes
+// from FORTRAN BLAS.
+matrix & daxpy(double a, matrix const & x, matrix & y)
+{
+    assert(x.rows() == y.rows());
+    assert(x.cols() == y.cols());
+    for (int i = 0; i < y.rows(); ++i) {
+        for (int j = 0; j < y.cols(); ++j) {
+            y(i, j) += a * x(i, j);
+        }
+    }
+    return y;
+}
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct self_evaluating_expr;
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+auto evaluate_matrix_expr(self_evaluating_expr<Kind, Tuple> const & expr);
+
+// This is the primary template for our expression template.  If you assign a
+// self_evaluating_expr to a matrix, its conversion operator transforms and
+// evaluates the expression with a call to evaluate_matrix_expr().
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct self_evaluating_expr
+{
+    operator auto() const;
+
+    static const boost::yap::expr_kind kind = Kind;
+
+    Tuple elements;
+};
+
+// This is a specialization of our expression template for assignment
+// expressions.  The destructor transforms and evaluates via a call to
+// evaluate_matrix_expr(), and then assigns the result to the variable on the
+// left side of the assignment.
+//
+// In a production implementation, you'd need to have specializations for
+// plus_assign, minus_assign, etc.
+template <typename Tuple>
+struct self_evaluating_expr<boost::yap::expr_kind::assign, Tuple>
+{
+    ~self_evaluating_expr();
+
+    static const boost::yap::expr_kind kind = boost::yap::expr_kind::assign;
+
+    Tuple elements;
+};
+
+struct use_daxpy
+{
+    // A plus-expression, which may be of the form double * matrix + matrix,
+    // or may be something else.  Since our daxpy() above requires a mutable
+    // "y", we only need to match a mutable lvalue matrix reference here.
+    template <typename Tuple>
+    auto operator()(
+        boost::yap::expr_tag<boost::yap::expr_kind::plus>,
+        self_evaluating_expr<boost::yap::expr_kind::multiplies, Tuple> const & expr,
+        matrix & m)
+    {
+        // Here, we transform the left-hand side into a pair if it's the
+        // double * matrix operation we're looking for.  Otherwise, we just
+        // get a copy of the left side expression.
+        //
+        // Note that this is a bit of a cheat, done for clarity.  If we pass a
+        // larger expression that happens to contain a double * matrix
+        // subexpression, that subexpression will be transformed into a tuple!
+        // In production code, this transform should probably only be
+        // performed on an expression with all terminal members.
+        auto lhs = boost::yap::transform(
+            expr,
+            [](boost::yap::expr_tag<boost::yap::expr_kind::multiplies>,
+               double d,
+               matrix const & m) {
+                return std::pair<double, matrix const &>(d, m);
+            });
+
+        // If we got back a copy of expr above, just re-construct the
+        // expression this function mathes; in other words, do not effectively
+        // transform anything.  Otherwise, replace the expression matched by
+        // this function with a call to daxpy().
+        if constexpr (boost::yap::is_expr<decltype(lhs)>::value) {
+            return expr + m;
+        } else {
+            return boost::yap::make_terminal(daxpy)(lhs.first, lhs.second, m);
+        }
+    }
+};
+
+
+// This is the heart of what self_evaluating_expr does.  If we had other
+// optimizations/transformations we wanted to do, we'd put them in this
+// function, either before or after the use_daxpy transformation.
+template <boost::yap::expr_kind Kind, typename Tuple>
+auto evaluate_matrix_expr(self_evaluating_expr<Kind, Tuple> const & expr)
+{
+    auto daxpy_form = boost::yap::transform(expr, use_daxpy{});
+    return boost::yap::evaluate(daxpy_form);
+}
+
+template<boost::yap::expr_kind Kind, typename Tuple>
+self_evaluating_expr<Kind, Tuple>::operator auto() const
+{
+    return evaluate_matrix_expr(*this);
+}
+
+template<typename Tuple>
+self_evaluating_expr<boost::yap::expr_kind::assign, Tuple>::
+    ~self_evaluating_expr()
+{
+    using namespace boost::hana::literals;
+    boost::yap::evaluate(elements[0_c]) = evaluate_matrix_expr(elements[1_c]);
+}
+
+// In order to define the = operator with the semantics we want, it's
+// convenient to derive a terminal type from a terminal instantiation of
+// self_evaluating_expr.  Note that we could have written a template
+// specialization here instead -- either one would work.  That would of course
+// have required more typing.
+struct self_evaluating :
+    self_evaluating_expr<
+        boost::yap::expr_kind::terminal,
+        boost::hana::tuple<matrix>
+    >
+{
+    self_evaluating() {}
+
+    explicit self_evaluating(matrix m)
+    { elements = boost::hana::tuple<matrix>(std::move(m)); }
+
+    BOOST_YAP_USER_ASSIGN_OPERATOR(self_evaluating_expr, ::self_evaluating_expr);
+};
+
+BOOST_YAP_USER_BINARY_OPERATOR(plus, self_evaluating_expr, self_evaluating_expr)
+BOOST_YAP_USER_BINARY_OPERATOR(minus, self_evaluating_expr, self_evaluating_expr)
+BOOST_YAP_USER_BINARY_OPERATOR(multiplies, self_evaluating_expr, self_evaluating_expr)
+
+
+int main()
+{
+    matrix identity(2, 2);
+    identity(0, 0) = 1.0;
+    identity(1, 1) = 1.0;
+
+    // These are YAP-ified terminal expressions.
+    self_evaluating m1(identity);
+    self_evaluating m2(identity);
+    self_evaluating m3(identity);
+
+    // This transforms the YAP expression to use daxpy(), so it creates no
+    // temporaries.  The transform happens in the destructor of the
+    // assignment-expression specialization of self_evaluating_expr.
+    m1 = 3.0 * m2 + m3;
+
+    // Same as above, except that it uses the matrix conversion operator on
+    // the self_evaluating_expr primary template, because here we're assigning
+    // a YAP expression to a non-YAP-ified matrix.
+    matrix m_result_1 = 3.0 * m2 + m3;
+
+    // Creates temporaries and does not use daxpy(), because the A * X + Y
+    // pattern does not occur within the expression.
+    matrix m_result_2 = 3.0 * m2;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/tarray.cpp b/src/boost/libs/yap/example/tarray.cpp
new file mode 100644
index 00000000..6aa4e109
--- /dev/null
+++ b/src/boost/libs/yap/example/tarray.cpp
@@ -0,0 +1,186 @@
+// 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)
+//[ tarray
+#include <boost/yap/algorithm.hpp>
+#include <boost/yap/print.hpp>
+
+#include <array>
+#include <iostream>
+
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct tarray_expr;
+
+
+struct take_nth
+{
+    boost::yap::terminal<tarray_expr, int>
+    operator() (boost::yap::terminal<tarray_expr, std::array<int, 3>> const & expr);
+
+    std::size_t n;
+};
+
+// Another custom expression template.  In this case, we static_assert() that
+// it only gets instantiated with terminals with pre-approved value types.
+template <boost::yap::expr_kind Kind, typename Tuple>
+struct tarray_expr
+{
+    // Make sure that, if this expression is a terminal, its value is one we
+    // want to support.  Note that the presence of expr_kind::expr_ref makes
+    // life slightly more difficult; we have to account for int const & and
+    // int & as well as int.
+    static_assert(
+        Kind != boost::yap::expr_kind::terminal ||
+        std::is_same<Tuple, boost::hana::tuple<int const &>>{} ||
+        std::is_same<Tuple, boost::hana::tuple<int &>>{} ||
+        std::is_same<Tuple, boost::hana::tuple<int>>{} ||
+        std::is_same<Tuple, boost::hana::tuple<std::array<int, 3>>>{},
+        "tarray_expr instantiated with an unsupported terminal type."
+    );
+
+    static const boost::yap::expr_kind kind = Kind;
+
+    Tuple elements;
+
+    int operator[] (std::size_t n) const
+    { return boost::yap::evaluate(boost::yap::transform(*this, take_nth{n})); }
+};
+
+// Define operators +, -, *, and /.
+BOOST_YAP_USER_BINARY_OPERATOR(plus, tarray_expr, tarray_expr)
+BOOST_YAP_USER_BINARY_OPERATOR(minus, tarray_expr, tarray_expr)
+BOOST_YAP_USER_BINARY_OPERATOR(multiplies, tarray_expr, tarray_expr)
+BOOST_YAP_USER_BINARY_OPERATOR(divides, tarray_expr, tarray_expr)
+
+
+boost::yap::terminal<tarray_expr, int>
+take_nth::operator() (boost::yap::terminal<tarray_expr, std::array<int, 3>> const & expr)
+{
+    int x = boost::yap::value(expr)[n];
+    // Again, this is the move hack to get x into the resulting terminal as a
+    // value instead of a reference.
+    return boost::yap::make_terminal<tarray_expr>(std::move(x));
+}
+
+
+// Stream-out operators for the two kinds of terminals we support.
+
+std::ostream & operator<< (std::ostream & os, boost::yap::terminal<tarray_expr, int> expr)
+{ return os << '{' << boost::yap::value(expr) << '}'; }
+
+std::ostream & operator<< (std::ostream & os, boost::yap::terminal<tarray_expr, std::array<int, 3>> expr)
+{
+    std::array<int, 3> const & a = boost::yap::value(expr);
+    return os << '{' << a[0] << ", " << a[1] << ", " << a[2] << '}';
+}
+
+// Stream-out operators for general expressions.  Note that we have to treat
+// the reference case separately; this also could have been done using
+// constexpr if in a single function template.
+
+template <typename Tuple>
+std::ostream & operator<< (std::ostream & os, tarray_expr<boost::yap::expr_kind::expr_ref, Tuple> const & expr)
+{ return os << boost::yap::deref(expr); }
+
+template <boost::yap::expr_kind Kind, typename Tuple>
+std::ostream & operator<< (std::ostream & os, tarray_expr<Kind, Tuple> const & expr)
+{
+    if (Kind == boost::yap::expr_kind::plus || Kind == boost::yap::expr_kind::minus)
+        os << '(';
+    os << boost::yap::left(expr) << " " << op_string(Kind) << " " << boost::yap::right(expr);
+    if (Kind == boost::yap::expr_kind::plus || Kind == boost::yap::expr_kind::minus)
+        os << ')';
+    return os;
+}
+
+
+// Since we want different behavior on terminals than on other kinds of
+// expressions, we create a custom type that does so.
+struct tarray :
+    tarray_expr<
+        boost::yap::expr_kind::terminal,
+        boost::hana::tuple<std::array<int, 3>>
+    >
+{
+    explicit tarray (int i = 0, int j = 0, int k = 0)
+    {
+        (*this)[0] = i;
+        (*this)[1] = j;
+        (*this)[2] = k;
+    }
+
+    explicit tarray (std::array<int, 3> a)
+    {
+        (*this)[0] = a[0];
+        (*this)[1] = a[1];
+        (*this)[2] = a[2];
+    }
+
+    int & operator[] (std::ptrdiff_t i)
+    { return boost::yap::value(*this)[i]; }
+
+    int const & operator[] (std::ptrdiff_t i) const
+    { return boost::yap::value(*this)[i]; }
+
+    template <typename T>
+    tarray & operator= (T const & t)
+    {
+        // We use as_expr() here to make sure that the value passed to
+        // assign() is an expression.  as_expr() simply forwards expressions
+        // through, and wraps non-expressions as terminals.
+        return assign(boost::yap::as_expr< ::tarray_expr>(t));
+    }
+
+    template <typename Expr>
+    tarray & printAssign (Expr const & expr)
+    {
+        *this = expr;
+        std::cout << *this << " = " << expr << std::endl;
+        return *this;
+    }
+
+private:
+    template <typename Expr>
+    tarray & assign (Expr const & expr)
+    {
+        (*this)[0] = expr[0];
+        (*this)[1] = expr[1];
+        (*this)[2] = expr[2];
+        return *this;
+    }
+};
+
+
+int main()
+{
+    tarray a(3,1,2);
+
+    tarray b;
+
+    std::cout << a << std::endl;
+    std::cout << b << std::endl;
+
+    b[0] = 7; b[1] = 33; b[2] = -99;
+
+    tarray c(a);
+
+    std::cout << c << std::endl;
+
+    a = 0;
+
+    std::cout << a << std::endl;
+    std::cout << b << std::endl;
+    std::cout << c << std::endl;
+
+    a = b + c;
+
+    std::cout << a << std::endl;
+
+    a.printAssign(b+c*(b + 3*c));
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/transform_terminals.cpp b/src/boost/libs/yap/example/transform_terminals.cpp
new file mode 100644
index 00000000..2bab93d6
--- /dev/null
+++ b/src/boost/libs/yap/example/transform_terminals.cpp
@@ -0,0 +1,72 @@
+// 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)
+//[ transform_terminals
+#include <boost/yap/yap.hpp>
+
+
+struct iota_terminal_transform
+{
+    // Base case. Note that we're not treating placeholders specially for this
+    // example (they're easy to special-case if necessary).
+    template<typename T>
+    auto operator()(boost::yap::expr_tag<boost::yap::expr_kind::terminal>, T && t)
+    {
+        // Like the std::iota() algorithm, we create replacement int terminals
+        // with the values index_, index_ + 1, index_ + 2, etc.
+        return boost::yap::make_terminal(index_++);
+    }
+
+    // Recursive case: Match any call expression.
+    template<typename CallableExpr, typename... Arg>
+    auto operator()(boost::yap::expr_tag<boost::yap::expr_kind::call>,
+                    CallableExpr callable, Arg &&... arg)
+    {
+        // Even though the callable in a call expression is technically a
+        // terminal, it doesn't make a lot of sense to replace it with an int,
+        // so we'll only transform the argument subexpressions.
+        return boost::yap::make_expression<boost::yap::expr_kind::call>(
+            boost::yap::as_expr(callable),
+            boost::yap::transform(boost::yap::as_expr(arg), *this)...);
+    }
+
+    int index_;
+};
+
+int sum(int a, int b) { return a + b; }
+
+int main()
+{
+    {
+        // This simple sum(8, 8) expression requires both overloads of
+        // iota_terminal_transform.
+        auto expr = boost::yap::make_terminal(sum)(8, 8);
+        assert(evaluate(expr) == 16);
+
+        auto iota_expr = boost::yap::transform(expr, iota_terminal_transform{1});
+        assert(evaluate(iota_expr) == 3);
+    }
+
+    {
+        // This expression requires only the terminal case of
+        // iota_terminal_transform.
+        auto expr = -(boost::yap::make_terminal(8) + 8);
+        assert(evaluate(expr) == -16);
+
+        auto iota_expr = boost::yap::transform(expr, iota_terminal_transform{0});
+        assert(evaluate(iota_expr) == -1);
+    }
+
+    {
+        // Like the first expression above, this expression requires both
+        // overloads of iota_terminal_transform.
+        auto expr = boost::yap::make_terminal(sum)(-(boost::yap::make_terminal(8) + 8), 0);
+        assert(evaluate(expr) == -16);
+
+        auto iota_expr = boost::yap::transform(expr, iota_terminal_transform{0});
+        assert(evaluate(iota_expr) == -3);
+    }
+}
+//]
diff --git a/src/boost/libs/yap/example/vec3.cpp b/src/boost/libs/yap/example/vec3.cpp
new file mode 100644
index 00000000..5a458dda
--- /dev/null
+++ b/src/boost/libs/yap/example/vec3.cpp
@@ -0,0 +1,129 @@
+// 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)
+//[ vec3
+#include <boost/yap/yap.hpp>
+
+#include <array>
+#include <iostream>
+
+
+struct take_nth
+{
+    auto operator() (boost::yap::terminal<boost::yap::expression, std::array<int, 3>> const & expr)
+    {
+        int x = boost::yap::value(expr)[n];
+        // The move forces the terminal to store the value of x, not a
+        // reference.
+        return boost::yap::make_terminal(std::move(x));
+    }
+
+    std::size_t n;
+};
+
+// Since this example doesn't constrain the operators defined on its
+// expressions, we can just use boost::yap::expression<> as the expression
+// template.
+using vec3_terminal = boost::yap::expression<
+    boost::yap::expr_kind::terminal,
+    boost::hana::tuple<std::array<int, 3>>
+>;
+
+// Customize the terminal type we use by adding index and assignment
+// operations.
+struct vec3 : vec3_terminal
+{
+    explicit vec3 (int i = 0, int j = 0, int k = 0)
+    {
+        (*this)[0] = i;
+        (*this)[1] = j;
+        (*this)[2] = k;
+    }
+
+    explicit vec3 (std::array<int, 3> a)
+    {
+        (*this)[0] = a[0];
+        (*this)[1] = a[1];
+        (*this)[2] = a[2];
+    }
+
+    int & operator[] (std::ptrdiff_t i)
+    { return boost::yap::value(*this)[i]; }
+
+    int const & operator[] (std::ptrdiff_t i) const
+    { return boost::yap::value(*this)[i]; }
+
+    template <typename T>
+    vec3 & operator= (T const & t)
+    {
+        decltype(auto) expr = boost::yap::as_expr(t);
+        (*this)[0] = boost::yap::evaluate(boost::yap::transform(expr, take_nth{0}));
+        (*this)[1] = boost::yap::evaluate(boost::yap::transform(expr, take_nth{1}));
+        (*this)[2] = boost::yap::evaluate(boost::yap::transform(expr, take_nth{2}));
+        return *this;
+    }
+
+    void print() const
+    {
+        std::cout << '{' << (*this)[0]
+                  << ", " << (*this)[1]
+                  << ", " << (*this)[2]
+                  << '}' << std::endl;
+    }
+};
+
+// This is a stateful transform that keeps a running count of the terminals it
+// has seen.
+struct count_leaves_impl
+{
+    auto operator() (vec3_terminal const & expr)
+    {
+        value += 1;
+        return expr;
+    }
+
+    int value = 0;
+};
+
+template <typename Expr>
+int count_leaves (Expr const & expr)
+{
+    count_leaves_impl impl;
+    boost::yap::transform(boost::yap::as_expr(expr), impl);
+    return impl.value;
+}
+
+
+int main()
+{
+    vec3 a, b, c;
+
+    c = 4;
+
+    b[0] = -1;
+    b[1] = -2;
+    b[2] = -3;
+
+    a = b + c;
+
+    a.print();
+
+    vec3 d;
+    auto expr1 = b + c;
+    d = expr1;
+    d.print();
+
+    int num = count_leaves(expr1);
+    std::cout << num << std::endl;
+
+    num = count_leaves(b + 3 * c);
+    std::cout << num << std::endl;
+
+    num = count_leaves(b + c * d);
+    std::cout << num << std::endl;
+
+    return 0;
+}
+//]
diff --git a/src/boost/libs/yap/example/vector.cpp b/src/boost/libs/yap/example/vector.cpp
new file mode 100644
index 00000000..e6bd185b
--- /dev/null
+++ b/src/boost/libs/yap/example/vector.cpp
@@ -0,0 +1,145 @@
+// 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)
+//[ vector
+#include <boost/yap/yap.hpp>
+
+#include <vector>
+#include <iostream>
+
+
+//[ vector_take_nth_xform
+struct take_nth
+{
+    template <typename T>
+    auto operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
+                     std::vector<T> const & vec)
+    { return boost::yap::make_terminal(vec[n]); }
+
+    std::size_t n;
+};
+//]
+
+// A stateful transform that records whether all the std::vector<> terminals
+// it has seen are equal to the given size.
+struct equal_sizes_impl
+{
+    template <typename T>
+    auto operator() (boost::yap::expr_tag<boost::yap::expr_kind::terminal>,
+                     std::vector<T> const & vec)
+    {
+        auto const expr_size = vec.size();
+        if (expr_size != size)
+            value = false;
+        return 0;
+    }
+
+    std::size_t const size;
+    bool value;
+};
+
+template <typename Expr>
+bool equal_sizes (std::size_t size, Expr const & expr)
+{
+    equal_sizes_impl impl{size, true};
+    boost::yap::transform(boost::yap::as_expr(expr), impl);
+    return impl.value;
+}
+
+
+// Assigns some expression e to the given vector by evaluating e elementwise,
+// to avoid temporaries and allocations.
+template <typename T, typename Expr>
+std::vector<T> & assign (std::vector<T> & vec, Expr const & e)
+{
+    decltype(auto) expr = boost::yap::as_expr(e);
+    assert(equal_sizes(vec.size(), expr));
+    for (std::size_t i = 0, size = vec.size(); i < size; ++i) {
+        vec[i] = boost::yap::evaluate(
+            boost::yap::transform(boost::yap::as_expr(expr), take_nth{i}));
+    }
+    return vec;
+}
+
+// As assign() above, just using +=.
+template <typename T, typename Expr>
+std::vector<T> & operator+= (std::vector<T> & vec, Expr const & e)
+{
+    decltype(auto) expr = boost::yap::as_expr(e);
+    assert(equal_sizes(vec.size(), expr));
+    for (std::size_t i = 0, size = vec.size(); i < size; ++i) {
+        vec[i] += boost::yap::evaluate(
+            boost::yap::transform(boost::yap::as_expr(expr), take_nth{i}));
+    }
+    return vec;
+}
+
+// Define a type trait that identifies std::vectors.
+template <typename T>
+struct is_vector : std::false_type {};
+
+template <typename T, typename A>
+struct is_vector<std::vector<T, A>> : std::true_type {};
+
+// Define all the expression-returning numeric operators we need.  Each will
+// accept any std::vector<> as any of its arguments, and then any value in the
+// remaining argument, if any -- some of the operators below are unary.
+BOOST_YAP_USER_UDT_UNARY_OPERATOR(negate, boost::yap::expression, is_vector); // -
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(multiplies, boost::yap::expression, is_vector); // *
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(divides, boost::yap::expression, is_vector); // /
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(modulus, boost::yap::expression, is_vector); // %
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(plus, boost::yap::expression, is_vector); // +
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(minus, boost::yap::expression, is_vector); // -
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(less, boost::yap::expression, is_vector); // <
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(greater, boost::yap::expression, is_vector); // >
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(less_equal, boost::yap::expression, is_vector); // <=
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(greater_equal, boost::yap::expression, is_vector); // >=
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(equal_to, boost::yap::expression, is_vector); // ==
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(not_equal_to, boost::yap::expression, is_vector); // !=
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(logical_or, boost::yap::expression, is_vector); // ||
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(logical_and, boost::yap::expression, is_vector); // &&
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(bitwise_and, boost::yap::expression, is_vector); // &
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(bitwise_or, boost::yap::expression, is_vector); // |
+BOOST_YAP_USER_UDT_ANY_BINARY_OPERATOR(bitwise_xor, boost::yap::expression, is_vector); // ^
+
+int main()
+{
+    int i;
+    int const n = 10;
+    std::vector<int> a,b,c,d;
+    std::vector<double> e(n);
+
+    for (i = 0; i < n; ++i)
+    {
+        a.push_back(i);
+        b.push_back(2*i);
+        c.push_back(3*i);
+        d.push_back(i);
+    }
+
+    // After this point, no allocations occur.
+
+    assign(b, 2);
+    assign(d, a + b * c);
+
+    a += if_else(d < 30, b, c);
+
+    assign(e, c);
+    e += e - 4 / (c + 1);
+
+    for (i = 0; i < n; ++i)
+    {
+        std::cout
+            << " a(" << i << ") = " << a[i]
+            << " b(" << i << ") = " << b[i]
+            << " c(" << i << ") = " << c[i]
+            << " d(" << i << ") = " << d[i]
+            << " e(" << i << ") = " << e[i]
+            << std::endl;
+    }
+
+    return 0;
+}
+//]