#ifndef LIB_MODELS_ODE_BOOST_TUPLE_HPP_
#define LIB_MODELS_ODE_BOOST_TUPLE_HPP_

//Author : pierre blavy




#include <tuple>
#include <type_traits>




namespace boost {
namespace numeric {
namespace odeint {

//user may specialize (entry point to make boost units works)
template<typename T>
struct std_tuple_algebra_to_double_t{
	static double run(const T&t){return static_cast<double>(t);}
};




template<typename T> double        std_tuple_algebra_to_double(const T&t)     {return std_tuple_algebra_to_double_t<T>::run(t);}
inline               const double& std_tuple_algebra_to_double(const double&t){return t;}


namespace std_tuple_algebra_impl{


	struct Observer_nothing{
		template<typename T, typename U>
		void operator()(const T&t, const U&){}
	};


	//--- norm inf ---
	template<size_t Index>
	struct norm_inf_t{
		template<typename T>
		static void run_r(const T&t, double& write_here){
			norm_inf_t<Index-1>::run_r(t,write_here);
			double d = std_tuple_algebra_to_double(std::get<Index-1>(t));
			if(d<0){d= -d;}
			if(d>write_here){write_here=d;}
		}
	};

	template <>
	struct norm_inf_t<0>{
		template<typename... stuff>
		static void run_r(stuff...){}
	};





    //--- 1 ---
	template <size_t Index>
	struct apply1_t{
		template<typename Op, typename T1>
		static void run_r(Op op, T1 &t1){
			apply1_t<Index-1>::run_r(std::forward<Op>(op),(t1));
			op(std::get<Index-1>(t1));
		}
	};

	template <>
	struct apply1_t<0>{
		template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};

//--- 2 ---
	template <size_t Index>
	struct apply2_t{
		template<typename Op, typename T1, typename T2>
		static void run_r(Op op, T1 &t1, T2 &t2){
			apply2_t<Index-1>::run_r(std::forward<Op>(op),(t1),(t2));
			op(std::get<Index-1>(t1), std::get<Index-1>(t2));
		}
	};

	template <>
	struct apply2_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};




	//--- 3 ---
	template <size_t Index>
	struct apply3_t{
		template<typename Op, typename T1, typename T2, typename T3>
		static void run_r(Op op, T1 &t1, T2 &t2, T3 &t3){
			apply3_t<Index-1>::run_r(std::forward<Op>(op),
					(t1),
					(t2),
					(t3)
			);
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3));
		}
	};

	template <>
	struct apply3_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};


	//--- 4 ---
	template <size_t Index>
	struct apply4_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4
		>
		static void run_r(Op op,
				T1 &t1, T2 &t2, T3 &t3,T4 &t4){
			apply4_t<Index-1>::run_r(std::forward<Op>(op),
					(t1),
					(t2),
					(t3),
					(t4));
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3),std::get<Index-1>(t4)
			);
		}
	};

	template <>
	struct apply4_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};

	//--- 5 ---
	template <size_t Index>
	struct apply5_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4, typename T5
		>
		static void run_r(Op op,
				T1 &t1, T2 &t2, T3 &t3,T4 &t4, T5 &t5){
			apply5_t<Index-1>::run_r(std::forward<Op>(op),
					(t1),
					(t2),
					(t3),
					(t4),
					(t5)
			);
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3),std::get<Index-1>(t4),std::get<Index-1>(t5));
		}
	};

	template <>
	struct apply5_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};

	//--- 6 ---
	template <size_t Index>
	struct apply6_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4, typename T5,
		typename T6
		>
		static void run_r(
				Op op,
				T1 &t1, T2 &t2, T3 &t3,T4 &t4, T5 &t5,
				T6 &t6){
			apply6_t<Index-1>::run_r(std::forward<Op>(op),
					(t1),
					(t2),
					(t3),
					(t4),
					(t5),
					(t6)
			);
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3),std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6)
			);
		}
	};

	template <>
	struct apply6_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};


	//--- 7 ---
	template <size_t Index>
	struct apply7_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4, typename T5,
		typename T6, typename T7
		>
		static void run_r(Op op,
				T1 &t1, T2 &t2, T3 &t3,T4 &t4, T5 &t5,
				T6 &t6, T7 &t7){
			apply7_t<Index-1>::run_r(std::forward<Op>(op),
					(t1),
					(t2),
					(t3),
					(t4),
					(t5),
					(t6),
					(t7)
			);
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3),std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6), std::get<Index-1>(t7)
			);
		}
	};

	template <>
	struct apply7_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};

	//--- 8 ---
	template <size_t Index>
	struct apply8_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4, typename T5,
		typename T6, typename T7, typename T8
		>
		static void run_r(Op op,
				T1 &t1, T2 &t2, T3 &t3,T4 &t4, T5 &t5,
				T6 &t6, T7 &t7, T8 &t8){
			apply8_t<Index-1>::run_r(std::forward<Op>(op),
				(t1),
				(t2),
				(t3),
				(t4),
				(t5),
				(t6),
				(t7),
				(t8)
			);
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3),std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6), std::get<Index-1>(t7), std::get<Index-1>(t8)
			);
		}
	};

	template <>
	struct apply8_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};


	//--- 9 ---
	template <size_t Index>
	struct apply9_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4, typename T5,
		typename T6, typename T7, typename T8, typename T9
		>
		static void run_r(Op op,
				T1 &t1, T2 &t2, T3 &t3, T4 &t4, T5 &t5,
				T6 &t6, T7 &t7, T8 &t8, T9 &t9){
			apply9_t<Index-1>::run_r(std::forward<Op>(op),
					(t1),
					(t2),
					(t3),
					(t4),
					(t5),
					(t6),
					(t7),
					(t8),
					(t9)
			);
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3), std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6), std::get<Index-1>(t7), std::get<Index-1>(t8), std::get<Index-1>(t9)
			);
		}
	};

	template <>
	struct apply9_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};


	//--- 10 ---
	template <size_t Index>
	struct apply10_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4, typename T5,
		typename T6, typename T7, typename T8, typename T9, typename T10
		>
		static void run_r(Op op,
				T1 &t1, T2 &t2, T3 &t3, T4 &t4, T5  &t5,
				T6 &t6, T7 &t7, T8 &t8, T9 &t9, T10 &t10){
			apply10_t<Index-1>::run_r(std::forward<Op>(op),
				(t1),
				(t2),
				(t3),
				(t4),
				(t5),
				(t6),
				(t7),
				(t8),
				(t9),
				(t10)
			);
			op(std::get<Index-1>(t1), std::get<Index-1>(t2), std::get<Index-1>(t3), std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6), std::get<Index-1>(t7), std::get<Index-1>(t8), std::get<Index-1>(t9),std::get<Index-1>(t10)
			);
		}
	};

	template <>
	struct apply10_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};

	//--- 11 ---
	template <size_t Index>
	struct apply11_t{
		template<typename Op,
		typename T1, typename T2, typename T3, typename T4, typename T5,
		typename T6, typename T7, typename T8, typename T9, typename T10,
		typename T11
		>
		static void run_r(Op op,
				T1 &t1  , T2 &t2, T3 &t3, T4 &t4, T5  &t5,
				T6 &t6  , T7 &t7, T8 &t8, T9 &t9, T10 &t10,
				T11 t11){
			apply11_t<Index-1>::run_r(std::forward<Op>(op),
				(t1),
				(t2),
				(t3),
				(t4),
				(t5),
				(t6),
				(t7),
				(t8),
				(t9),
				(t10),
				(t11)
			);
			op(std::get<Index-1>(t1) , std::get<Index-1>(t2) , std::get<Index-1>(t3) , std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6) , std::get<Index-1>(t7) , std::get<Index-1>(t8) , std::get<Index-1>(t9),std::get<Index-1>(t10),
			   std::get<Index-1>(t11)
			);
		}
	};

	template <>
	struct apply11_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};

	//--- 12 ---
	template <size_t Index>
	struct apply12_t{
		template<typename Op,
		typename T1 , typename T2 , typename T3 , typename T4 , typename T5,
		typename T6 , typename T7 , typename T8 , typename T9 , typename T10,
		typename T11, typename T12
		>
		static void run_r(Op op,
				T1 &t1  , T2  &t2 , T3 &t3, T4 &t4, T5  &t5,
				T6 &t6  , T7  &t7 , T8 &t8, T9 &t9, T10 &t10,
				T11 &t11, T12 &t12){
			apply12_t<Index-1>::run_r(std::forward<Op>(op),
				(t1),
				(t2),
				(t3),
				(t4),
				(t5),
				(t6),
				(t7),
				(t8),
				(t9),
				(t10),
				(t11),
				(t12)
			);
			op(std::get<Index-1>(t1) , std::get<Index-1>(t2) , std::get<Index-1>(t3) , std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6) , std::get<Index-1>(t7) , std::get<Index-1>(t8) , std::get<Index-1>(t9),std::get<Index-1>(t10),
			   std::get<Index-1>(t11), std::get<Index-1>(t12)
			);
		}
	};

	template <>
	struct apply12_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};


	//--- 13 ---
	template <size_t Index>
	struct apply13_t{
		template<typename Op,
		typename T1 , typename T2 , typename T3 , typename T4 , typename T5,
		typename T6 , typename T7 , typename T8 , typename T9 , typename T10,
		typename T11, typename T12, typename T13
		>
		static void run_r(Op op,
				T1 &t1  , T2  &t2 , T3 &t3, T4 &t4, T5  &t5,
				T6 &t6  , T7  &t7 , T8 &t8, T9 &t9, T10 &t10,
				T11 &t11, T12 &t12, T13 &t13){
			apply13_t<Index-1>::run_r(std::forward<Op>(op),
				(t1),
				(t2),
				(t3),
				(t4),
				(t5),
				(t6),
				(t7),
				(t8),
				(t9),
				(t10),
				(t11),
				(t12),
				(t13)
			);
			op(std::get<Index-1>(t1) , std::get<Index-1>(t2) , std::get<Index-1>(t3) , std::get<Index-1>(t4),std::get<Index-1>(t5),
			   std::get<Index-1>(t6) , std::get<Index-1>(t7) , std::get<Index-1>(t8) , std::get<Index-1>(t9),std::get<Index-1>(t10),
			   std::get<Index-1>(t11), std::get<Index-1>(t12), std::get<Index-1>(t13)
			);
		}
	};

	template <>
	struct apply13_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};


	//--- 14 ---
	template <size_t Index>
	struct apply14_t{
		template<typename Op,
		typename T1 , typename T2 , typename T3 , typename T4 , typename T5,
		typename T6 , typename T7 , typename T8 , typename T9 , typename T10,
		typename T11, typename T12, typename T13, typename T14
		>
		static void run_r(Op op,
				T1 &t1  , T2  &t2 , T3 &t3, T4 &t4, T5  &t5,
				T6 &t6  , T7  &t7 , T8 &t8, T9 &t9, T10 &t10,
				T11 &t11, T12 &t12, T13 &t13,T14 &t14){
			apply14_t<Index-1>::run_r(std::forward<Op>(op),
				(t1),
				(t2),
				(t3),
				(t4),
				(t5),
				(t6),
				(t7),
				(t8),
				(t9),
				(t10),
				(t11),
				(t12),
				(t13),
				(t14)
			);
			op(std::get<Index-1>(t1) , std::get<Index-1>(t2) , std::get<Index-1>(t3) , std::get<Index-1>(t4) ,std::get<Index-1>(t5),
			   std::get<Index-1>(t6) , std::get<Index-1>(t7) , std::get<Index-1>(t8) , std::get<Index-1>(t9) ,std::get<Index-1>(t10),
			   std::get<Index-1>(t11), std::get<Index-1>(t12), std::get<Index-1>(t13), std::get<Index-1>(t14)
			);
		}
	};

	template <>
	struct apply14_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};

	//--- 15 ---
	template <size_t Index>
	struct apply15_t{
		template<typename Op,
		typename T1 , typename T2 , typename T3 , typename T4 , typename T5,
		typename T6 , typename T7 , typename T8 , typename T9 , typename T10,
		typename T11, typename T12, typename T13, typename T14, typename T15
		>
		static void run_r(Op op,
				T1 &t1  , T2  &t2 , T3 &t3, T4 &t4, T5  &t5,
				T6 &t6  , T7  &t7 , T8 &t8, T9 &t9, T10 &t10,
				T11 &t11, T12 &t12, T13 &t13,T14 &t14, T15 &t15){
			apply15_t<Index-1>::run_r(std::forward<Op>(op),
				(t1),
				(t2),
				(t3),
				(t4),
				(t5),
				(t6),
				(t7),
				(t8),
				(t9),
				(t10),
				(t11),
				(t12),
				(t13),
				(t14),
				(t15)
			);
			op(std::get<Index-1>(t1) , std::get<Index-1>(t2) , std::get<Index-1>(t3) , std::get<Index-1>(t4) ,std::get<Index-1>(t5),
			   std::get<Index-1>(t6) , std::get<Index-1>(t7) , std::get<Index-1>(t8) , std::get<Index-1>(t9) ,std::get<Index-1>(t10),
			   std::get<Index-1>(t11), std::get<Index-1>(t12), std::get<Index-1>(t13), std::get<Index-1>(t14),std::get<Index-1>(t15)
			);
		}
	};

	template <>
	struct apply15_t<0>{
		 template<typename Op, typename... stuff>
		 static void run_r(Op op, stuff...){}
	};


}




struct std_tuple_algebra{

	template<typename... Args>
	static double norm_inf(const std::tuple<Args...> &t){
		double sum=0;
		std_tuple_algebra_impl::norm_inf_t<sizeof...(Args)>::run_r(t,sum);
		return sum;
	}




   //--- 1 ---
	template<typename T1, class Operation>
    static void for_each1(
    	T1 &t1,
   		Operation op )
     {
     	std_tuple_algebra_impl::apply1_t<std::tuple_size<T1>::value>::run_r(std::forward<Operation>(op),(t1));
     }

	//--- 2 ---
    template<
		typename T1, typename T2,
		class Operation
	>
    static void for_each2(
    	T1  &t1 , T2  &t2,
		Operation op )
     {
    	static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each2 for t2");
    	std_tuple_algebra_impl::apply2_t<std::tuple_size<T1>::value>::run_r(std::forward<Operation>(op),(t1),(t2));
     }



    //---3---
    template<
		typename T1, typename T2, typename T3,
		class Operation
	>
    static void for_each3(
    	T1  &t1 , T2  &t2 , T3  &t3 ,
   		Operation op )
     {
    	static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each3 for t2");
    	static_assert(std::tuple_size<T1>::value == std::tuple_size<T3>::value, "size mismatch in for_each3 for t3");
    	std_tuple_algebra_impl::apply3_t<std::tuple_size<T1>::value>::run_r(std::forward<Operation>(op),t1,t2,t3);
     }

    //---4---
    template<
		typename T1, typename T2, typename T3, typename T4,
		class Operation
	>
    static void for_each4(
    	T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 ,
   		Operation op )
     {
    	static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each4 for t2");
    	static_assert(std::tuple_size<T1>::value == std::tuple_size<T3>::value, "size mismatch in for_each4 for t3");
    	static_assert(std::tuple_size<T1>::value == std::tuple_size<T4>::value, "size mismatch in for_each4 for t4");

    	std_tuple_algebra_impl::apply4_t<std::tuple_size<T1>::value>::run_r(std::forward<Operation>(op),(t1),(t2),(t3),(t4));
     }

    //---5---
      template<
  		typename T1, typename T2, typename T3, typename T4,typename T5,
  		class Operation
  	>
      static void for_each5(
    		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
     		Operation op )
       {
      	static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each5 for t2");
      	static_assert(std::tuple_size<T1>::value == std::tuple_size<T3>::value, "size mismatch in for_each5 for t3");
      	static_assert(std::tuple_size<T1>::value == std::tuple_size<T4>::value, "size mismatch in for_each5 for t4");
      	static_assert(std::tuple_size<T1>::value == std::tuple_size<T5>::value, "size mismatch in for_each5 for t5");

      	std_tuple_algebra_impl::apply5_t<std::tuple_size<T1>::value>::run_r(std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5));
       }

      //---6---
	  template<
		typename T1, typename T2, typename T3, typename T4,typename T5,
		typename T6,
		class Operation
	  >
	  static void for_each6(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
	    T6  &t6 ,
		Operation op )
	   {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each6 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3>::value, "size mismatch in for_each6 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4>::value, "size mismatch in for_each6 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5>::value, "size mismatch in for_each6 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6>::value, "size mismatch in for_each6 for t6");

		std_tuple_algebra_impl::apply6_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
			 (t6)
		);
	   }

	  //---7---
	  template<
		typename T1, typename T2, typename T3, typename T4,typename T5,
		typename T6, typename T7,
		class Operation
	  >
	  static void for_each7(
	    T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 ,
		Operation op )
	   {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each7 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3>::value, "size mismatch in for_each7 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4>::value, "size mismatch in for_each7 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5>::value, "size mismatch in for_each7 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6>::value, "size mismatch in for_each7 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7>::value, "size mismatch in for_each7 for t7");

		std_tuple_algebra_impl::apply7_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
				   (t6),(t7)
		);
	   }


	  //---8---
	  template<
		typename T1, typename T2, typename T3, typename T4,typename T5,
		typename T6, typename T7, typename T8,
		class Operation
	  >
	  static void for_each8(
	    T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 ,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each8 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3>::value, "size mismatch in for_each8 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4>::value, "size mismatch in for_each8 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5>::value, "size mismatch in for_each8 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6>::value, "size mismatch in for_each8 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7>::value, "size mismatch in for_each8 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8>::value, "size mismatch in for_each8 for t8");

		std_tuple_algebra_impl::apply8_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
			   (t6),(t7),(t8)
		);
	   }

	  //---9---
	  template<
		typename T1, typename T2, typename T3, typename T4,typename T5,
		typename T6, typename T7, typename T8, typename T9,
		class Operation
	  >
	  static void for_each9(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 , T9  &t9,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2>::value, "size mismatch in for_each9 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3>::value, "size mismatch in for_each9 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4>::value, "size mismatch in for_each9 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5>::value, "size mismatch in for_each9 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6>::value, "size mismatch in for_each9 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7>::value, "size mismatch in for_each9 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8>::value, "size mismatch in for_each9 for t8");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T9>::value, "size mismatch in for_each9 for t9");

		std_tuple_algebra_impl::apply9_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
		     (t6),(t7),(t8),(t9)
		);
	   }

	  //---10---
	  template<
		typename T1, typename T2, typename T3, typename T4,typename T5,
		typename T6, typename T7, typename T8, typename T9,typename T10,
		class Operation
	  >
	  static void for_each10(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 , T9  &t9 , T10 &t10,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2 >::value, "size mismatch in for_each10 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3 >::value, "size mismatch in for_each10 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4 >::value, "size mismatch in for_each10 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5 >::value, "size mismatch in for_each10 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6 >::value, "size mismatch in for_each10 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7 >::value, "size mismatch in for_each10 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8 >::value, "size mismatch in for_each10 for t8");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T9 >::value, "size mismatch in for_each10 for t9");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T10>::value, "size mismatch in for_each10 for t10");

		std_tuple_algebra_impl::apply10_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
			   (t6),(t7),(t8),(t9),(t10)
		);
	   }

	  //---11---
	  template<
		typename T1 , typename T2 , typename T3 , typename T4 ,typename T5 ,
		typename T6 , typename T7 , typename T8 , typename T9 ,typename T10,
		typename T11,
		class Operation
	  >
	  static void for_each11(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 , T9  &t9 , T10 &t10,
		T11 &t11,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2 >::value, "size mismatch in for_each11 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3 >::value, "size mismatch in for_each11 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4 >::value, "size mismatch in for_each11 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5 >::value, "size mismatch in for_each11 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6 >::value, "size mismatch in for_each11 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7 >::value, "size mismatch in for_each11 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8 >::value, "size mismatch in for_each11 for t8");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T9 >::value, "size mismatch in for_each11 for t9");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T10>::value, "size mismatch in for_each11 for t10");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T11>::value, "size mismatch in for_each11 for t11");

		std_tuple_algebra_impl::apply11_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
			   (t6),(t7),(t8),(t9),(t10),
			   (t11)
		);
	   }

	  //---12---
	  template<
		typename T1 , typename T2 , typename T3 , typename T4 ,typename T5 ,
		typename T6 , typename T7 , typename T8 , typename T9 ,typename T10,
		typename T11, typename T12,
		class Operation
	  >
	  static void for_each12(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 , T9  &t9 , T10 &t10,
		T11 &t11, T12 &t12,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2 >::value, "size mismatch in for_each12 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3 >::value, "size mismatch in for_each12 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4 >::value, "size mismatch in for_each12 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5 >::value, "size mismatch in for_each12 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6 >::value, "size mismatch in for_each12 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7 >::value, "size mismatch in for_each12 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8 >::value, "size mismatch in for_each12 for t8");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T9 >::value, "size mismatch in for_each12 for t9");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T10>::value, "size mismatch in for_each12 for t10");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T11>::value, "size mismatch in for_each12 for t11");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T12>::value, "size mismatch in for_each12 for t12");


		std_tuple_algebra_impl::apply12_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
				   (t6),(t7),(t8),(t9),(t10),
				   (t11),  (t12)
		);
	   }


	  //---13---
	  template<
		typename T1 , typename T2 , typename T3 , typename T4 ,typename T5 ,
		typename T6 , typename T7 , typename T8 , typename T9 ,typename T10,
		typename T11, typename T12, typename T13,
		class Operation
	  >
	  static void for_each13(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 , T9  &t9 , T10 &t10,
		T11 &t11, T12 &t12, T13 &t13,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2 >::value, "size mismatch in for_each13 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3 >::value, "size mismatch in for_each13 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4 >::value, "size mismatch in for_each13 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5 >::value, "size mismatch in for_each13 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6 >::value, "size mismatch in for_each13 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7 >::value, "size mismatch in for_each13 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8 >::value, "size mismatch in for_each13 for t8");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T9 >::value, "size mismatch in for_each13 for t9");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T10>::value, "size mismatch in for_each13 for t10");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T11>::value, "size mismatch in for_each13 for t11");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T12>::value, "size mismatch in for_each13 for t12");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T13>::value, "size mismatch in for_each13 for t13");


		std_tuple_algebra_impl::apply13_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
				   (t6),(t7),(t8),(t9),(t10),
				   (t11),(t12),(t13)
		);
	   }

	  //---14---
	  template<
		typename T1 , typename T2 , typename T3 , typename T4 ,typename T5 ,
		typename T6 , typename T7 , typename T8 , typename T9 ,typename T10,
		typename T11, typename T12, typename T13, typename T14,
		class Operation
	  >
	  static void for_each14(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 , T9  &t9 , T10 &t10,
		T11 &t11, T12 &t12, T13 &t13, T14 &t14,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2 >::value, "size mismatch in for_each14 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3 >::value, "size mismatch in for_each14 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4 >::value, "size mismatch in for_each14 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5 >::value, "size mismatch in for_each14 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6 >::value, "size mismatch in for_each14 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7 >::value, "size mismatch in for_each14 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8 >::value, "size mismatch in for_each14 for t8");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T9 >::value, "size mismatch in for_each14 for t9");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T10>::value, "size mismatch in for_each14 for t10");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T11>::value, "size mismatch in for_each14 for t11");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T12>::value, "size mismatch in for_each14 for t12");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T13>::value, "size mismatch in for_each14 for t13");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T14>::value, "size mismatch in for_each14 for t14");


		std_tuple_algebra_impl::apply14_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2),(t3),(t4),(t5),
				  (t6),(t7),(t8),(t9),(t10),
				  (t11),(t12),(t13),(t14)
		);
	   }

	  //---14---
	  template<
		typename T1 , typename T2 , typename T3 , typename T4 ,typename T5 ,
		typename T6 , typename T7 , typename T8 , typename T9 ,typename T10,
		typename T11, typename T12, typename T13, typename T14,typename T15,
		class Operation
	  >
	  static void for_each15(
		T1  &t1 , T2  &t2 , T3  &t3 , T4  &t4 , T5  &t5 ,
		T6  &t6 , T7  &t7 , T8  &t8 , T9  &t9 , T10 &t10,
		T11 &t11, T12 &t12, T13 &t13, T14 &t14, T15 &t15,
		Operation op )
	  {
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T2 >::value, "size mismatch in for_each15 for t2");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T3 >::value, "size mismatch in for_each15 for t3");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T4 >::value, "size mismatch in for_each15 for t4");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T5 >::value, "size mismatch in for_each15 for t5");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T6 >::value, "size mismatch in for_each15 for t6");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T7 >::value, "size mismatch in for_each15 for t7");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T8 >::value, "size mismatch in for_each15 for t8");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T9 >::value, "size mismatch in for_each15 for t9");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T10>::value, "size mismatch in for_each15 for t10");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T11>::value, "size mismatch in for_each15 for t11");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T12>::value, "size mismatch in for_each15 for t12");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T13>::value, "size mismatch in for_each15 for t13");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T14>::value, "size mismatch in for_each15 for t14");
		static_assert(std::tuple_size<T1>::value == std::tuple_size<T15>::value, "size mismatch in for_each15 for t15");


		std_tuple_algebra_impl::apply15_t<std::tuple_size<T1>::value>::run_r(
				std::forward<Operation>(op),(t1),(t2)  ,(t3),(t4),(t5),
			(t6)   ,(t7)  ,(t8),(t9),(t10),
			(t11) ,(t12),(t13),(t14),(t15)
		);
	   }






};

}}}






//---Register algebra in the dispatcher ---
#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
namespace boost {
namespace numeric {
namespace odeint {
	template< typename... Args >
	struct algebra_dispatcher< std::tuple<Args...> >
	{
		typedef std_tuple_algebra algebra_type;
	};

}}}







//default integrate function TODO require more testing
/*
#include <functional>

namespace boost {
namespace numeric {
namespace odeint {




	template< typename...State_args, typename... Speed_args, typename Time_t, typename Oberver_t>
    size_t integrate_const(
		std::function<void(const std::tuple<State_args...>&, std::tuple<Speed_args...>&, const Time_t t)> speed ,
		const std::tuple<State_args...>	&init_state ,
		const Time_t &begin,
		const Time_t &end,
		const Time_t &step,
		Oberver_t    &observer=std_tuple_algebra_impl::Observer_nothing()
	){
	typedef std::tuple<State_args...> State_type;
	typedef std::tuple<Speed_args...> Speed_type;
	typedef runge_kutta_dopri5< State_type , double , Speed_type , Time_t > stepper_type;

	return integrate_const(
			make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) ,
			speed,
			init_state ,
			begin ,
			end ,
			step ,
			observer
	);
    }


	template<
	  typename...State_args,
	  typename... Speed_args,
	  typename Time_t,
	  typename Oberver_t
	>
	size_t integrate(
		std::function<
		  void(
		  const std::tuple<State_args...>&,
		        std::tuple<Speed_args...>&,
		  const Time_t t)
		> speed ,
		const std::tuple<State_args...>	&init_state ,
		const Time_t &begin,
		const Time_t &end,
		const Time_t &step,
		Oberver_t    &observer=std_tuple_algebra_impl::Observer_nothing()
	){
	typedef std::tuple<State_args...> State_type;
	typedef std::tuple<Speed_args...> Speed_type;
	typedef controlled_runge_kutta< runge_kutta_dopri5< State_type , double , Speed_type , Time_t > > stepper_type;

	return integrate_adaptive(
			stepper_type() ,
			speed,
			init_state ,
			begin,
			end,
			step,
			observer
	);
	}


}}}

*/





#endif /* LIB_MODELS_ODE_BOOST_TUPLE_HPP_ */