#include <kv/poincaremap.hpp>
#include <kv/kraw-approx.hpp>
#include <kv/ode-maffine.hpp>

namespace ub = boost::numeric::ublas;

typedef kv::interval<double> itv;

struct VDP {
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x, T t){
		ub::vector<T> y(2);

		y(0) = x(1);
		y(1) = 1. * (1. - x(0)*x(0)) * x(1) - x(0);

		return y;
	}
};

struct VDPPoincareSection {
	template <class T> T operator() (const ub::vector<T>& x){
		T y;

		y = x(0) - 0.;

		return y;
	}
};


int main()
{
	ub::vector<double> x;
	ub::vector<itv> ix;
	int i;
	bool r;

	std::cout.precision(17);

	VDP f;
	VDPPoincareSection ps;
	kv::PoincareMap<VDP,VDPPoincareSection,double> po(f, ps, (itv)0.);

	x.resize(3);
	x(0) = 0.; 
	x(1) = 2.;
	x(2) = 6.;

	r = kv::krawczyk_approx(po, x, ix, 5, 0);
	if (!r) return 0;

	std::list< kv::interval<double> > time_list;
	std::list< ub::vector< kv::interval<double> > > value_list;
	ub::vector<itv> ix2;
	itv end;

	ix2.resize(2);
	ix2(0) = ix(0);
	ix2(1) = ix(1);
	end = ix(2);

	r = kv::odelong_maffine(f, ix2, itv(0.), end, kv::ode_param<double>(), kv::ode_callback_dense_list<double>(itv(0.), itv(0.05), time_list, value_list)); 

	time_list.push_back(end);
	value_list.push_back(ix2);

	std::list< kv::interval<double> >::iterator pt;
	std::list< ub::vector< kv::interval<double> > >::iterator pv;
	pt = time_list.begin();
	pv = value_list.begin();

	while (pt != time_list.end()) {
		std::cout << mid(*pt) << " ";
		for (i=0; i<2; i++) {
			std::cout << mid((*pv)(i)) << " ";
		}
		std::cout << "\n";
		pt++; pv++;
	}
}