#include "ibex.h"

#include <vector>
#include <queue>

#define DECOMPOSE 1

// Macro variables to set the values of the maximal length of a
// pattern (NB_K) and the maximal number of splittings of the
// state-space (NB_D)
#define NB_K 4
#define NB_D 1
#define NB_M 16 //number of modes 
#define NB_scale 4 //scaling of the original system
#define SW_period 5.0
#define CS_period 2.5

using namespace ibex;

// ************************************************************
// MINIMATOR FUNCTIONS ****************************************
// ************************************************************


// Data type for sampled switched systems
typedef struct  {
	double period;
	//std::vector<Function*> dynamics;
} sampledSwitchedSystem;

typedef struct {
	IntervalVector *Yinit;
	IntervalVector *Ycurrent;
	std::vector<int> pattern;  
} node;


// A process to compute the successor of binary word of a given length 
bool nextPattern(std::vector<int>& pattern) {
	for (int i = pattern.size() - 1; i >= 0; i-- ){
		pattern[i] += 1 ;
		if ( pattern[i] > NB_M-1 ){
			pattern[i] = 0 ;
			continue ;
		}
		return true ;
	}
	return false;
}


// Binary strings generation function
void printTheArray(int arr[], int n, int res[][NB_scale], int * cpt) 
{ 
	for (int i = 0; i < n; i++) { 
		res[*cpt][i] = arr[i];
	       	cout << arr[i] << " "; 
	} 
	(*cpt)++;
	cout << endl; 
} 

// Binary strings generation function (2)
void generateAllBinaryStrings(int n, int arr[], int i, int res[][NB_scale], int * cpt) 
{ 
	if (i == n) { 
		printTheArray(arr, n, res, cpt); 
		cpt++;
		return; 
		} 
	arr[i] = 0; 
	generateAllBinaryStrings(n, arr, i + 1, res, cpt); 
	arr[i] = 1; 
	generateAllBinaryStrings(n, arr, i + 1, res, cpt); 
} 


// Function that finds a control pattern for box W
bool findPattern2 (const sampledSwitchedSystem& sys, const IntervalVector W,
		  const IntervalVector R1, const IntervalVector R2, const IntervalVector R3, const IntervalVector R4, const IntervalVector B,  const IntervalVector S1, const IntervalVector S2, const IntervalVector S3, const IntervalVector S4, unsigned int k,
		  std::vector<int>&  res_pattern, void * place_function1, void * place_ivp_ode1, void * place_simulation1, void * place_function2, void * place_ivp_ode2, void * place_simulation2, void * place_function3, void * place_ivp_ode3, void * place_simulation3, void * place_function4, void * place_ivp_ode4, void * place_simulation4) {

	node node_init;
	node_init.Yinit = new IntervalVector(W);
	node_init.Ycurrent = new IntervalVector(W);
	node_init.pattern = std::vector<int>();
	  
	std::list<node> list_node;
	  
	list_node.push_back(node_init);
	  
	int nb_pattern = NB_M;
	  
	std::cout << "W : " << W << std::endl;

	// Model parameters 
	const double T_e = 10.0;
	const double a_12 = 0.05;
	const double a_21 = 0.05;
	const double a_e1 = 0.005;
	const double a_e2 = 0.0033;
	const double a_f =  0.0083;
	const double T_f = 50;

	int m = NB_scale; 
	int arr[m]; 
	int test = pow(2.0,m);
	int res[test][NB_scale];
	int cpt = 0;
	generateAllBinaryStrings(m, arr, 0, res, &cpt); 

	bool fin_R1;
	bool stay_S1;
	bool out_S1;

	bool fin_R2;
	bool stay_S2;
	bool out_S2;

	bool fin_R3;
	bool stay_S3;
	bool out_S3;

	bool fin_R4;
	bool stay_S4;
	bool out_S4;

	double Tn = 0.0;

	// Cosimulation parameters 
	double H = CS_period;      
	int nb_cosim_period = SW_period/CS_period;

	// Declaration of local variables 
	Variable y1(3);
	Variable y2(3);
	Variable y3(3);
	Variable y4(3);

	IntervalVector picard(8);
	IntervalVector yin1_test(3);
	IntervalVector yin2_test(3);
	IntervalVector yin3_test(3);
	IntervalVector yin4_test(3);
	IntervalVector yinit1(3);
	IntervalVector yinit2(3);
	IntervalVector yinit3(3);
	IntervalVector yinit4(3);
	IntervalVector ypic1(3);
	IntervalVector ypic2(3);
	IntervalVector ypic3(3);
	IntervalVector ypic4(3);
	IntervalVector yin1(1);
	IntervalVector yin2(2);
	IntervalVector yin3(2);
	IntervalVector yin4(1);


	while(!list_node.empty())
	{
	    
		node node_current = list_node.front();
		list_node.pop_front();
	    
		for (int i_mode=0; i_mode < nb_pattern; i_mode++)
		{
	      
			if (node_current.pattern.empty())
				std::cout << "current pattern : " << i_mode << std::endl;
			else
				std::cout << "current pattern : " << node_current.pattern << " + " << i_mode << std::endl;
		      
			IntervalVector temp_init = *node_current.Ycurrent;

			// Initialisation of local variables
			yinit1[0] = temp_init[0];
			yinit1[1] = temp_init[1];
			yinit1[2] = Interval(0.0);

			yinit2[0] = temp_init[2];
			yinit2[1] = temp_init[3];
			yinit2[2] = Interval(0.0);

			yinit3[0] = temp_init[4];
			yinit3[1] = temp_init[5];
			yinit3[2] = Interval(0.0);

			yinit4[0] = temp_init[6];
			yinit4[1] = temp_init[7];
			yinit4[2] = Interval(0.0);
		      
			for (int i=0; i<nb_cosim_period; i++)
			{
				int test = 0;

				yin1_test[0] = temp_init[0];
				yin1_test[1] = temp_init[1];
				yin1_test[2] = Interval(Tn,Tn+H);

				yin2_test[0] = temp_init[2];
				yin2_test[1] = temp_init[3];
				yin2_test[2] = Interval(Tn,Tn+H);

				yin3_test[0] = temp_init[4];
				yin3_test[1] = temp_init[5];
				yin3_test[2] = Interval(Tn,Tn+H);

				yin4_test[0] = temp_init[6];
				yin4_test[1] = temp_init[7];
				yin4_test[2] = Interval(Tn,Tn+H);


				ypic1[0] = Interval::ALL_REALS;
				ypic1[1] = Interval::ALL_REALS;
				ypic1[2] = Interval::ALL_REALS;
			
				ypic2[0] = Interval::ALL_REALS;
				ypic2[1] = Interval::ALL_REALS;
				ypic2[2] = Interval::ALL_REALS;

				ypic3[0] = Interval::ALL_REALS;
				ypic3[1] = Interval::ALL_REALS;
				ypic3[2] = Interval::ALL_REALS;

				ypic4[0] = Interval::ALL_REALS;
				ypic4[1] = Interval::ALL_REALS;
				ypic4[2] = Interval::ALL_REALS;

				yin1[0] = yin2_test[0]; 

				yin2[0] = yin1_test[1];
				yin2[1] = yin3_test[0];

				yin3[0] = yin2_test[1];
				yin3[1] = yin4_test[0];

				yin4[0] = yin3_test[1];


				// Cross-picard computation
				while((((!ypic2.is_subset(yin2_test)) && (!ypic1.is_subset(yin1_test)) && (!ypic3.is_subset(yin3_test)) && (!ypic4.is_subset(yin4_test))) || (test<2)) && (test<100))
				{
					yin1_test = ypic1;
					yin2_test = ypic2;
					yin3_test = ypic3;
					yin4_test = ypic4;
					test++;
			      


					Function * ypic11 = new(place_function1) Function (y1, Return( (-1)*(a_21+a_e1+res[i_mode][0]*a_f)*y1[0]+a_21*y1[1]+a_e1*T_e+res[i_mode][0]*a_f*T_f,
									a_12*y1[0]+(-1)*(a_12+a_21+a_e2)*y1[1]+a_21*yin1[0]+a_e2*T_e,
							Interval(1.0) ));



					ivp_ode * problem_pic1 = new(place_ivp_ode1) ivp_ode(*ypic11,i*H,yinit1,SYMBOLIC);
					simulation * simu_pic1 = new(place_simulation1) simulation(problem_pic1,(i+1)*H,HEUN,1e-3);
					simu_pic1->setHmax(CS_period);
					simu_pic1->setHmin(CS_period);
					solution_j_rk4 u_j_pic1(Affine2Vector(yinit1,true),i*H,H, problem_pic1,CS_period,CS_period);
					IntervalVector test1 = u_j_pic1.picard_tayl(yinit1,problem_pic1,2);
					ypic1 = test1;



					Function * ypic22 = new(place_function2) Function (y2, Return(a_12*yin2[0]+(-1)*(a_12+a_21+a_e1+res[i_mode][1]*a_f)*y2[0]+a_21*y2[1]+a_e1*T_e+res[i_mode][1]*a_f*T_f,
									a_12*y2[0]+(-1)*(a_12+a_21+a_e2)*y2[1]+a_21*yin2[1]+a_e2*T_e,
							Interval(1.0) ));
				
					ivp_ode * problem_pic2 = new(place_ivp_ode2) ivp_ode(*ypic22,i*H,yinit2,SYMBOLIC);
					simulation * simu_pic2 = new(place_simulation2) simulation(problem_pic2,(i+1)*H,HEUN,1e-3);
					simu_pic2->setHmax(CS_period);
					simu_pic2->setHmin(CS_period);
					solution_j_rk4 u_j_pic2(Affine2Vector(yinit2,true),i*H,H, problem_pic2,CS_period,CS_period);
					IntervalVector test2 = u_j_pic2.picard_tayl(yinit2,problem_pic2,2);
		   			ypic2 = test2;	



					Function * ypic33 = new(place_function3) Function (y3, Return( a_12*yin3[0]+(-1)*(a_12+a_21+a_e1+res[i_mode][2]*a_f)*y3[0]+a_21*y3[1]+a_e1*T_e+res[i_mode][2]*a_f*T_f,
									a_12*y3[0]+(-1)*(a_12+a_21+a_e2)*y3[1]+a_21*yin3[1]+a_e2*T_e,
							Interval(1.0) )) ;

					ivp_ode * problem_pic3 = new(place_ivp_ode3) ivp_ode(*ypic33,i*H,yinit3,SYMBOLIC);
					simulation * simu_pic3 = new(place_simulation3) simulation(problem_pic3,(i+1)*H,HEUN,1e-3);
					simu_pic3->setHmax(CS_period);
					simu_pic3->setHmin(CS_period);
					solution_j_rk4 u_j_pic3(Affine2Vector(yinit3,true),i*H,H, problem_pic3,CS_period,CS_period);
					IntervalVector test3 = u_j_pic3.picard_tayl(yinit3,problem_pic3,2);
			   		ypic3 = test3;	


					Function * ypic44 = new(place_function4) Function (y4, Return( a_12*yin4[0]+(-1)*(a_12+a_21+a_e1+res[i_mode][3]*a_f)*y4[0]+a_21*y4[1]+a_e1*T_e+res[i_mode][3]*a_f*T_f,
										a_12*y4[0]+(-1)*(a_12+a_e2)*y4[1]+a_e2*T_e,
								Interval(1.0) )) ;

					ivp_ode * problem_pic4 = new(place_ivp_ode4) ivp_ode(*ypic44,i*H,yinit4,SYMBOLIC);
					simulation * simu_pic4 = new(place_simulation4) simulation(problem_pic4,(i+1)*H,HEUN,1e-3);
					simu_pic4->setHmax(CS_period);
					simu_pic4->setHmin(CS_period);
					solution_j_rk4 u_j_pic4(Affine2Vector(yinit4,true),i*H,H, problem_pic4,CS_period,CS_period);
					IntervalVector test4 = u_j_pic4.picard_tayl(yinit4,problem_pic4,2);
			   		ypic4 = test4;	
						

					ypic11->~Function();
					ypic22->~Function();
					ypic33->~Function();
					ypic44->~Function();
					problem_pic1->~ivp_ode();
					problem_pic2->~ivp_ode();
					problem_pic3->~ivp_ode();
					problem_pic4->~ivp_ode();
					simu_pic1->~simulation();
					simu_pic2->~simulation();
					simu_pic3->~simulation();
					simu_pic4->~simulation();

					yin1[0] = ypic2[0]; 

					yin2[0] = ypic1[1];
					yin2[1] = ypic3[0];

					yin3[0] = ypic2[1];
					yin3[1] = ypic4[0];

					yin4[0] = ypic3[1];

					picard[0] = ypic1[0];
					picard[1] = ypic1[1];
					picard[2] = ypic2[0];
					picard[3] = ypic2[1];
					picard[4] = ypic3[0];
					picard[5] = ypic3[1];
					picard[6] = ypic4[0];
					picard[7] = ypic4[1];

				}


				// Cosimulation step
				{
					yin1[0] = picard[2]; 

					yin2[0] = picard[1];
					yin2[1] = picard[4];

					yin3[0] = picard[3];
					yin3[1] = picard[6];

					yin4[0] = picard[5];

					#pragma omp parallel
					{
						#pragma omp sections
						{
							{

								Function * ydot1 = new(place_function1) Function (y1, Return( (-1)*(a_21+a_e1+res[i_mode][0]*a_f)*y1[0]+a_21*y1[1]+a_e1*T_e+res[i_mode][0]*a_f*T_f,
													a_12*y1[0]+(-1)*(a_12+a_21+a_e2)*y1[1]+a_21*yin1[0]+a_e2*T_e,
											Interval(1.0) ));


								ivp_ode * problem1 = new(place_ivp_ode1) ivp_ode(*ydot1,i*H,yinit1,SYMBOLIC);
								simulation * simu1 = new(place_simulation1) simulation(problem1,(i+1)*H,HEUN,1e-4);

								simu1->setHmax(0.9);
								simu1->setHmin(1e-6);

								simu1->run_simulation();
								yinit1 = simu1->get_last();

								fin_R1 = simu1->finished_in(R1);
								stay_S1 = simu1->stayed_in(S1);
								out_S1 = simu1->go_out(S1);

								ydot1->~Function();
								problem1->~ivp_ode();
								simu1->~simulation();

							}
							
							#pragma omp section
							{

								Function * ydot2 = new(place_function2) Function (y2, Return(a_12*yin2[0]+(-1)*(a_12+a_21+a_e1+res[i_mode][1]*a_f)*y2[0]+a_21*y2[1]+a_e1*T_e+res[i_mode][1]*a_f*T_f,
												a_12*y2[0]+(-1)*(a_12+a_21+a_e2)*y2[1]+a_21*yin2[1]+a_e2*T_e,
										Interval(1.0) ));

								ivp_ode * problem2 = new(place_ivp_ode2) ivp_ode(*ydot2,i*H,yinit2,SYMBOLIC);
								simulation * simu2 = new(place_simulation2) simulation(problem2,(i+1)*H,HEUN,1e-4);
								simu2->run_simulation();
								yinit2 = simu2->get_last();

								fin_R2 = simu2->finished_in(R2);
								stay_S2 = simu2->stayed_in(S2);
								out_S2 = simu2->go_out(S2);

								ydot2->~Function();
								problem2->~ivp_ode();
								simu2->~simulation();

							}

							#pragma omp section
							{

								Function * ydot3 = new(place_function3) Function (y3, Return( a_12*yin3[0]+(-1)*(a_12+a_21+a_e1+res[i_mode][2]*a_f)*y3[0]+a_21*y3[1]+a_e1*T_e+res[i_mode][2]*a_f*T_f,
												a_12*y3[0]+(-1)*(a_12+a_21+a_e2)*y3[1]+a_21*yin3[1]+a_e2*T_e,
										Interval(1.0) )) ;

								ivp_ode * problem3 = new(place_ivp_ode3) ivp_ode(*ydot3,i*H,yinit3,SYMBOLIC);
								simulation * simu3 = new(place_simulation3) simulation(problem3,(i+1)*H,HEUN,1e-4);
								simu3->run_simulation();
								yinit3 = simu3->get_last();

								fin_R3 = simu3->finished_in(R3);
								stay_S3 = simu3->stayed_in(S3);
								out_S3 = simu3->go_out(S3);

								ydot3->~Function();
								problem3->~ivp_ode();
								simu3->~simulation();

							}

							#pragma omp section
							{

								Function * ydot4 = new(place_function4) Function (y4, Return( a_12*yin4[0]+(-1)*(a_12+a_21+a_e1+res[i_mode][3]*a_f)*y4[0]+a_21*y4[1]+a_e1*T_e+res[i_mode][3]*a_f*T_f,
												a_12*y4[0]+(-1)*(a_12+a_e2)*y4[1]+a_e2*T_e,
										Interval(1.0) )) ;

								ivp_ode * problem4 = new(place_ivp_ode4) ivp_ode(*ydot4,i*H,yinit4,SYMBOLIC);
								simulation * simu4 = new(place_simulation4) simulation(problem4,(i+1)*H,HEUN,1e-4);
								simu4->run_simulation();
								yinit4 = simu4->get_last();

								fin_R4 = simu4->finished_in(R4);
								stay_S4 = simu4->stayed_in(S4);
								out_S4 = simu4->go_out(S4);

								ydot4->~Function();
								problem4->~ivp_ode();
								simu4->~simulation();

							}
						}
					}
				}	
			}

	      
			if (fin_R1 && fin_R2 && fin_R3 && fin_R4 && stay_S1 && stay_S2 && stay_S3 && stay_S4)
			{
				std::cout << "sol found !" << std::endl;
				node_current.pattern.push_back(i_mode);

				res_pattern = node_current.pattern;

				return true;
			}
			else 
			{
				if (out_S1 || out_S2 || out_S3 || out_S4) 
				{
					std::cout << "Wrong direction !" << std::endl;
				}
				else 
				{
					if (stay_S1 && stay_S2 && stay_S3 && stay_S4)
				  	{
				   		if (node_current.pattern.size() < k)
				    		{
			     			std::cout << "Increment of pattern !" << std::endl;	      

			    			node new_node;
			      			new_node.Yinit = new IntervalVector(*node_current.Yinit);

					      	IntervalVector temp(8);
					     	temp[0] = yinit1[0];
					  	temp[1] = yinit1[1];

						temp[2] = yinit2[0];
						temp[3] = yinit2[1];

						temp[4] = yinit3[0];
						temp[5] = yinit3[1];

						temp[6] = yinit4[0];
						temp[7] = yinit4[1];

						new_node.Ycurrent = new IntervalVector(temp);
						      
						std::vector<int> new_pattern (node_current.pattern);
						new_pattern.push_back(i_mode);
						new_node.pattern = new_pattern;
						list_node.push_back(new_node);	      

			    			}
			  		}
				}
			}
		}
	}
	return false;
}



// The main algorithm of minimator
bool decompose (const sampledSwitchedSystem& sys, const IntervalVector W,
		const IntervalVector R, const IntervalVector R1, const IntervalVector R2, const IntervalVector R3, const IntervalVector R4, const IntervalVector B, const IntervalVector S, const IntervalVector S1, const IntervalVector S2, const IntervalVector S3, const IntervalVector S4, unsigned int k, unsigned int d, std::vector<std::pair<IntervalVector, std::vector<int> > >& result, void * place_function1, void * place_ivp_ode1, void * place_simulation1, void * place_function2, void * place_ivp_ode2, void * place_simulation2, void * place_function3, void * place_ivp_ode3, void * place_simulation3, void * place_function4, void * place_ivp_ode4, void * place_simulation4) {
  std::vector<int> res_pattern;

	  // q is a queue of sub-state-space to explore to prove the invariance
	std::queue<IntervalVector> q;
	q.push (W);

	unsigned int nbStep = d;

	while (!q.empty() && nbStep > 0) {
		IntervalVector current = q.front ();
		q.pop();
		bool flag = findPattern2 (sys, current, R1, R2, R3, R4, B, S1, S2, S3, S4, k, res_pattern, place_function1, place_ivp_ode1, place_simulation1, place_function2, place_ivp_ode2, place_simulation2, place_function3, place_ivp_ode3, place_simulation3, place_function4, place_ivp_ode4, place_simulation4);
		if (flag) {
			// Validated result, we keep it in the vector of results
			result.push_back (std::pair< IntervalVector, std::vector<int> > (current, res_pattern));
	    	}
	    	else {
			LargestFirst bbb(0.001,0.5);
			std::pair<IntervalVector,IntervalVector> p = bbb.bisect(current);
			q.push (p.first);
			q.push (p.second);
			nbStep--;
	    	}
	}

	if (nbStep <= 0) {
		return false;
	}
	else {
		return true;
	}
}

// ********************************************************
// ********************************************************
// ********************************************************


double sommeVector (const Vector &v) {
	double res = 0.0;
	for (int i = 0; i < v.size(); i++) {
		res += v[i];
	}
	return res;
}



int main(){
	int m = NB_scale; 
	int arr[m]; 
	int test = pow(2.0,m);
	int res[test][NB_scale];
	int cpt = 0;
	generateAllBinaryStrings(m, arr, 0, res, &cpt); 

	std::cout << res[0][0] << res[0][1] << res[0][2] << res[0][3] << std::endl;
	std::cout << res[test-1][0] << res[test-1][1] << res[test-1][2] << res[test-1][3] << std::endl;

	const double T_e = 10.0;
	const double a_12 = 0.05;
	const double a_21 = 0.05;
	const double a_e1 = 0.005;
	const double a_e2 = 0.0033;
	const double a_f =  0.0083;
	const double T_f = 50;

	const int n = 8;//2*NB_scale;
	Variable y(n);


	sampledSwitchedSystem sys;
	sys.period = SW_period;

	IntervalVector R(8);           // Recurrence set for global system
	R[0]= Interval(19,21);
	R[1]= Interval(19,21);
	R[2]= Interval(19,21);
	R[3]= Interval(19,21);
	R[4]= Interval(19,21);
	R[5]= Interval(19,21);
	R[6]= Interval(19,21);
	R[7]= Interval(19,21);


	IntervalVector R1(3);           // Recurrence set for system 1
	R1[0]= Interval(19,21);
	R1[1]= Interval(19,21);
	R1[2]= Interval(-1000,1000);

	IntervalVector R2(3);           // Recurrence set for system 2
	R2[0]= Interval(19,21);
	R2[1]= Interval(19,21);
	R2[2]= Interval(-1000,1000);

	IntervalVector R3(3);           // Recurrence set for system 3
	R3[0]= Interval(19,21);
	R3[1]= Interval(19,21);
	R3[2]= Interval(-1000,1000);

	IntervalVector R4(3);           // Recurrence set for system 4
	R4[0]= Interval(19,21);
	R4[1]= Interval(19,21);
	R4[2]= Interval(-1000,1000);

	#ifdef DECOMPOSE

	IntervalVector W(8);              // Starting set
	W[0]= Interval(19,21);
	W[1]= Interval(19,21);
	W[2]= Interval(19,21);
	W[3]= Interval(19,21);
	W[4]= Interval(19,21);
	W[5]= Interval(19,21);
	W[6]= Interval(19,21);
	W[7]= Interval(19,21);

	IntervalVector B(8);               // Forbidden/avoid set
	B[0] = Interval(-10.0,-9.01);
	B[1] = Interval(-10.0,-9.01);
	B[2] = Interval(-10.0,-9.01);
	B[3] = Interval(-10.0,-9.01);
	B[4] = Interval(-10.0,-9.01);
	B[5] = Interval(-10.0,-9.01);
	B[6] = Interval(-10.0,-9.01);
	B[7] = Interval(-10.0,-9.01);

	IntervalVector S(8);		// Safety set for global system
	S[0] = Interval(18, 22);
	S[1] = Interval(18, 22);
	S[2] = Interval(18, 22);
	S[3] = Interval(18, 22);
	S[4] = Interval(18, 22);
	S[5] = Interval(18, 22);
	S[6] = Interval(18, 22);
	S[7] = Interval(18, 22);

	IntervalVector S1(3);          // Recurrence set for system 1
	S1[0] = Interval(18, 22);
	S1[1] = Interval(18, 22);
	S1[2] = Interval(-2000,2000);

	IntervalVector S2(3);		// Recurrence set for system 2
	S2[0] = Interval(18, 22);
	S2[1] = Interval(18, 22);
	S2[2] = Interval(-2000,2000);

	IntervalVector S3(3);		// Recurrence set for system 3
	S3[0] = Interval(18, 22);
	S3[1] = Interval(18, 22);
	S3[2] = Interval(-2000,2000);

	IntervalVector S4(3);		// Recurrence set for system 4
	S4[0] = Interval(18, 22);
	S4[1] = Interval(18, 22);
	S4[2] = Interval(-2000,2000);

	std::list<IntervalVector> list_W;
	list_W.push_back(W);


	char memory_function1[sizeof(Function)];
	void* place_function1 = memory_function1;
	char memory_ivp_ode1[sizeof(ivp_ode)];
	void* place_ivp_ode1 = memory_ivp_ode1;
	char memory_simulation1[sizeof(simulation)];
	void* place_simulation1 = memory_simulation1;

	char memory_function2[sizeof(Function)];
	void* place_function2 = memory_function2;
	char memory_ivp_ode2[sizeof(ivp_ode)];
	void* place_ivp_ode2 = memory_ivp_ode2;
	char memory_simulation2[sizeof(simulation)];
	void* place_simulation2 = memory_simulation2;

	char memory_function3[sizeof(Function)];
	void* place_function3 = memory_function3;
	char memory_ivp_ode3[sizeof(ivp_ode)];
	void* place_ivp_ode3 = memory_ivp_ode3;
	char memory_simulation3[sizeof(simulation)];
	void* place_simulation3 = memory_simulation3;

	char memory_function4[sizeof(Function)];
	void* place_function4 = memory_function4;
	char memory_ivp_ode4[sizeof(ivp_ode)];
	void* place_ivp_ode4 = memory_ivp_ode4;
	char memory_simulation4[sizeof(simulation)];
	void* place_simulation4 = memory_simulation4;

	// File for solutions
	std::ofstream file("./result_minimator_8_rooms_cosim.txt", std::ios::out | std::ios::trunc);

	while(!list_W.empty())
	{
		IntervalVector current_W = list_W.front();
		list_W.pop_front();
	
		// Pre-bisection of the state-space 
		if (current_W.diam().max() > 1.0)
		{
			LargestFirst bbb(0.1,0.5);
			std::pair<IntervalVector,IntervalVector> p = bbb.bisect(current_W);
			list_W.push_back(p.first);
			list_W.push_back(p.second);
			continue;
		}


	std::cout << "taille list_W : " << list_W.size() << std::endl;

	std::vector< std::pair<IntervalVector, std::vector<int> > > result;
	unsigned int k = NB_K;
	unsigned int d = NB_D;

	bool flag = decompose(sys, current_W, R, R1, R2, R3, R4, B, S, S1, S2, S3, S4, k, d, result, place_function1, place_ivp_ode1, place_simulation1, place_function2, place_ivp_ode2, place_simulation2, place_function3, place_ivp_ode3, place_simulation3, place_function4, place_ivp_ode4, place_simulation4);

	if (result.empty()) {
		std::cerr << "No solution with k = " << k << " and d = " << d << std::endl;
		file << current_W << " : no sol" << std::endl;
	}
	else
	{
		if (flag) {
		std::cerr << "Complete result -> PROOF" << std::endl;
		}
		else {
		std::cerr << "Incomplete result" << std::endl;
		}
	  
		// Display result
		std::vector< std::pair<IntervalVector, std::vector<int> > >::const_iterator it = result.begin();
		for (; it != result.end(); it++) {
			file << "zon = str2zon2(\"" ;
			file << it->first << "\"); \n if is_in_zonotope(X,zon) \t \n pattern = [";

			std::vector<int>::const_iterator it_pat = (it->second).begin();
			for (; it_pat != (it->second).end(); it_pat++) {
				file << *it_pat << " ";
			}
			file << "]; \n endif" << std::endl;
			}
		}
	}
	file.close();
	#endif

	return 0;
}