#ifndef _WIN32
#include "../include/lminmaxgd.h"
#else
#include "..\include\lminmaxgd.h"
using namespace std;
#endif
/////////////////////////////////////////////////////////////////////
/* Example : Observer Synthesis, Example taken from paper
[1] arXiv:0901.2915v1 (authors DiLoreto, Gaubert Katz, Loiseau
Synthesis down thanks to results in
[2] "Observer Design for max,+ linear systems
(IEEE TAC, 2010, Hardouin, Maia, Cottenceau, Lhommeau)

The goal of this programm is to compute with interval analysis the
observer matrix L =[L_inf L_sup]
and to compare with the one given in [2]

// you must link with a library which contains all source files given in folder src
// gd.cpp poly.cpp serie.cpp smatrix.cpp tools.cpp interf.cpp
// or add these files in your project
// or alternatively you can include the folowing files
*/
#include "..\src\gd.cpp"
#include "..\src\poly.cpp"
#include "..\src\serie.cpp"
#include "..\src\smatrix.cpp"
#include "..\src\tools.cpp"

/////////////////////////////////////////////////////////////////////


int main()
{
try
{

	smatrix A(9,9);	// matrice d'état
    smatrix A_inf(9,9);	// matrice d'état
    smatrix A_sup(9,9);	// matrice d'état
	smatrix L,L1,L2;
	smatrix L_inf,L1_inf,L2_inf;
	smatrix L_sup,L1_sup,L2_sup;
	int i,j;

	smatrix CA,AB,CAB,AR,CAR;

	smatrix CA_sup,AB_sup,CAB_sup,AR_sup,CAR_sup;
	smatrix CA_inf,AB_inf,CAB_inf,AR_inf,CAR_inf;


	A(0,2)=gd(1,4);//(1,4)
	A(0,6)=gd(1,2);//(1,2)

	A(1,0)=gd(1,1);
	A(1,7)=gd(1,3);

	A(2,1)=gd(1,5);
	A(2,8)=gd(1,1);

	A(3,0)=gd(1,4);
	A(3,5)=gd(1,3);

	A(4,1)=gd(1,3);
	A(4,3)=gd(1,1);

	A(5,2)=gd(1,5);
	A(5,4)=gd(1,4);

	A(6,3)=gd(1,4);
	A(6,8)=gd(1,3);

	A(7,6)=gd(1,5);
	A(7,4)=gd(1,3);

    A(8,7)=gd(1,4);
	A(8,5)=gd(1,2);

    A_inf=A;
    A_sup=A;
    A_sup(1,0)=gd(1,7);
    A_sup(1,7)=gd(1,3);
    A_sup(4,3)=gd(1,3);
    A_sup(4,1)=gd(1,5);

	smatrix C(9,9);	// output matrix only coefficient on diagonl equal to e or epsilon
    smatrix C_inf(9,9);	// output matrix only coefficient on diagonl equal to e or epsilon
    smatrix C_sup(9,9);	// output matrix only coefficient on diagonl equal to e or epsilon
   // C(0,0)=gd(0,0);
    //C(1,1)=gd(0,0);
 	C(2,2)=gd(0,0); // in paper [1] the first sensor is her
	//C(3,3)=gd(0,0);
    //C(4,4)=gd(0,0);
	C(5,5)=gd(0,0); // in paper [1] the second sensor is here
    //C(6,6)=gd(0,0);
	C(7,7)=gd(0,0); // in paper [1] the third sensor is here
    //C(8,8)=gd(0,0);
	/* */
// test we can change the position of sensors to achieve the equality, you have to comment the previous three
// line and to replace by the following one : in this case all the sate are perfectly estimated
// since the sufficient condition of proposition 4 of paper [2], see reference in the head, is fulfilled
/*	C(1,1)=gd(0,0);
	C(4,4)=gd(0,0);
	C(7,7)=gd(0,0);
*/
    C_inf=C;
    C_sup=C;

    smatrix B(9,1);	// matrice entrée et sortie, epsilon everywhere
	smatrix R(9,3);
	smatrix B_inf(9,1);	// matrice entrée et sortie, epsilon everywhere
	smatrix R_inf(9,3);
	smatrix B_sup(9,1);	// matrice entrée et sortie, epsilon everywhere
	smatrix R_sup(9,3);


	R(1,0)=gd(0,0);
	R(4,0)=gd(0,0);
    R(1,0)=gd(0,0);
    R_inf=R;
    R_sup=R;
	smatrix As=star(A);
	smatrix As_inf=star(A);
	smatrix As_sup=star(A);

	CA=otimes(C,As);  // CA*
    CA_inf=otimes(C,As);  // CA*
    CA_sup=otimes(C,As);  // CA*

    AB=otimes(As,B);
    AB_inf=otimes(As_inf,B_inf);
    AB_sup=otimes(As_sup,B_sup);
    CAB=otimes(C,AB);
    CAB_inf=otimes(C_inf,AB_inf);
    CAB_sup=otimes(C_sup,AB_sup);
    AR=otimes(As,R);
    AR_inf=otimes(As_inf,R_inf);
    AR_sup=otimes(As_sup,R_sup);
    CAR=otimes(C,AR);
    CAR_inf=otimes(C_inf,AR_inf);
    CAR_sup=otimes(C_sup,AR_sup);

	L1=rfrac(AB,CAB);
	L2=rfrac(AR,CAR);
	L=inf(L1,L2);
	L=prcaus(L);

	L1_inf=rfrac(AB_inf,CAB_inf);
	L2_inf=rfrac(AR_inf,CAR_inf);
	L_inf=inf(L1_inf,L2_inf);
	L_inf=prcaus(L_inf);

    L1_sup=rfrac(AB_sup,CAB_sup);
	L2_sup=rfrac(AR_sup,CAR_sup);
	L_sup=inf(L1_sup,L2_sup);
	L_sup=prcaus(L_sup);

    L_inf=inf(L_inf,L_sup); // projection in the set of intervals
    if(L_inf==L)
    cout<<"it is the same"<<endl; // that means the upper bounds of the intervals have no sense
                                // on this example, the only point which is take into account itis the lower bound




    smatrix X;
    smatrix Xestim;
    smatrix ALCA,LC,LCA,LCs,ALCAR;
    smatrix ALCA_sup,LC_sup,LCA_sup,LCs_sup,ALCAR_sup;
    X=AR;
    cout<<" X = AR"<<X<<endl;
    LC=otimes(L,C);
    LCA=otimes(LC,As);
    LCs=oplus(LC,A);
    LCs=star(LCs);
    ALCA=otimes(LCs,LCA);
    ALCAR=otimes(ALCA,R);

    LC_sup=otimes(L_sup,C_sup);
    LCA_sup=otimes(LC_sup,As_sup);
    LCs_sup=oplus(LC_sup,A_sup);
    LCs_sup=star(LCs_sup);
    ALCA_sup=otimes(LCs_sup,LCA_sup);
    ALCAR_sup=otimes(ALCA_sup,R_sup);

    if (ALCAR==AR) cout<<" it is OK"<<endl;


	cout << "Xestim :" << ALCAR << endl;
    cout << "X :" << AR << endl;
   for(i=0;i<9;i++)
   { for (j=0;j<3;j++)
   {
    if (ALCAR(i,j)==AR(i,j))
    cout<<"equality between x"<<i+1<<"and xestimate"<<i+1<<endl; // that means x and xestimate will be equal whatever be the disturbance
    else {
        cout<<"NO equality for transition labelled :"<<i+1<<endl; // for these transition equality is not achieved but the best is done
        cout<<"It is due to  difference between AR("<<i+1<<","<<j+1<<") "<<AR(i,j)<<endl;
        cout<<"and ALCAR("<<i+1<<","<<j+1<<") "<<ALCAR(i,j)<<endl;
    }

    if (ALCAR_sup(i,j)==AR_sup(i,j))
    cout<<"equality between x"<<i+1<<"and xestimate_sup"<<i+1<<endl; // that means x and xestimate will be equal whatever be the disturbance
    else {
        cout<<"NO equality for transition labelled :"<<i+1<<endl; // for these transition equality is not achieved but the best is done
        cout<<"It is due to  difference between AR("<<i+1<<","<<j+1<<") "<<AR(i,j)<<endl;
        cout<<"and ALCAR_sup("<<i+1<<","<<j+1<<") "<<ALCAR(i,j)<<endl;
    }

   }
   }
       if (ALCAR==AR) cout<<" it is OK"<<endl;

//cout<<"CA "<<CA<<endl;
//cout<<"L "<<L<<endl;
	return(0);
}

 catch(mem_limite l)
 {
	 cout<<"Exception : too many coefficents in polynom "<<l.memoire<<endl;
	 return(1);
 }

 catch(taille_incorrecte obj)
 { // 0 : r non causal
   // 1 : tentative d'accès à un element d'une matrice avec un indice incorrect
   // 2 : matrice de taille incompatible pour oplus, inf, otimes, rfrac, lfrac
   // 3 : etoile de matrice carrée uniquement
	 cout<<"Exception  "<<obj.erreur<<endl;
	 return(1);
 }
}
