#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <kv/autodif.hpp>

namespace ub = boost::numeric::ublas;

template <class T> ub::vector<T> f(const ub::vector<T>& x) {
	ub::vector<T> y(2);

	y(0) = 2. * x(0) * x(0) * x(1) - 1.;
	y(1) = x(0) + 0.5 * x(1) * x(1) - 2.;

	return y;
}

int main()
{
	ub::vector<double> v1, v2;
	ub::vector< kv::autodif<double> > va1, va2;
	ub::matrix<double> m;

	v1.resize(2);
	v1(0) = 5.; v1(1) = 6.;

	va1 = kv::autodif<double>::init(v1);

	va2 = f(va1);

	kv::autodif<double>::split(va2, v2, m);
	
	std::cout << v2 << "\n"; // f(5, 6)
	std::cout << m << "\n"; // Jacobian matrix
}