#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
*/
using namespace std;
#include "..\src\gd.cpp"
#include "..\src\poly.cpp"
#include "..\src\serie.cpp"
#include "..\src\smatrix.cpp"
#include "..\src\tools.cpp"

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


int main()
{
try
{
    serie s1,s2,s3,sinf,ssup;
    poly p1;

    poly q1;
    poly p2;
    poly q2;
    gd r1;

    gd r2;



	smatrix A0(6,6);	// matrice d'état
	smatrix A1(6,6);	// matrice d'état
	smatrix A(6,6);	// matrice d'état
    smatrix B0(6,2);	// matrice d'état
    smatrix B(6,2);	// matrice d'état
    smatrix E(6,6);
    smatrix Id(6,6);
    smatrix ES;
    smatrix AES;
    smatrix G0,F0,F1,G1,G2,F2,tampon,tampon1;
    smatrix D0,D1;
    smatrix t;
    int n=0,exit=0;
	int i,j;


 p1.init(0,2)(4,3)(5,8);
    q1.init(8,10);
    r1.init(2,2);
    p2.init(1,1)(2,4)(4,5)(6,7);
    q2.init(7,9);
r2.init(2,2);
    s1.init(p1,q1,r1);
    s2.init(p2,q2,r2);

  cout<<"s1 "<<s1<<endl;
  cout<<"s2 "<<s2<<endl;
  sinf=inf(s1,s2);
  ssup=oplus(s1,s2);
  s3=frac(sinf,ssup);
  cout<<"sinf/sup"<<s3<<endl;

//s3=\gamma \delta^{-2} \oplus \gamma^2 \delta^{-1} \oplus \gamma^3  \oplus \gamma^4 \delta \oplus \gamma^5 \delta^3 (\gamma \delta)^*

  for(i=0;i<6;i++)
  Id(i,i)=gd(0,0);

	A0(1,0)=gd(0,10);//(1,4)
	A0(3,2)=gd(0,35);//(1,2)

	A0(4,1)=gd(0,2);

	A0(4,3)=gd(0,2);

	A0(5,4)=gd(0,50);
	A1(0,1)=gd(0,0);
	A1(2,3)=gd(0,0);
	A1(4,5)=gd(0,0);
    for(i=0;i<6;i++)
	A1(i,i)=gd(0,0);

   B0(0,0)=gd(0,6);
   B0(2,1)=gd(0,12);


E(3,1)=gd(0,0);
E(1,3)=gd(0,0);
E(0,5)=gd(0,-65);
E(2,5)=gd(0,-90);
E=oplus(E,A0);

ES=star(E);


	A=star(A0);
	B=otimes(A,B0);
	A=otimes(A,A1);

  cout<<" A0*A1 "<<A<<endl;
  cout<<" A0*B0 "<<B<<endl;
 // cout<<" ES "<<ES<<endl;
  AES=otimes(A,ES);

  G0=lfrac(AES,ES); //G0= E*\AE*
  t(0,0)=gd(0,100);
  D0=Id;
    for (i=1;i<6;i++)
    D0(i,i)=otimes(t(0,0),D0(i,i));
  G0=otimes(G0,D0);
  if(AES==otimes(ES,G0)) cout<<"OK"<<endl;
  else {
      cout<<"NOK"<<endl;
      cout<<"AES "<<AES<<endl;
      tampon=otimes(ES,G0);
      for(i=0;i<6;i++)
      for(j=0;j<6;j++)

       cout<<"E* (E* \ AE*) "<<tampon(i,j)<<" AE*"<<AES(i,j)<<endl;
  getchar();
  }




  tampon=otimes(ES,G0); //   E*G
  D0=rfrac(tampon,AES);
  t(0,0)=D0(0,0);
    for (i=1;i<6;i++)
    t(0,0)=inf(t(0,0),D0(i,i));
    D0=Id;
    for (i=1;i<6;i++)
    D0(i,i)=otimes(t(0,0),D0(i,i));


  tampon=lfrac(tampon,B); // B\E*G
  F0=rfrac(tampon,ES); //  B\E*G/E*
  F0=prcaus(F0);
  exit=0;
  n=0;
  do
  {
    tampon=otimes(D0,AES);
    tampon1=otimes(F0,ES);
    tampon1=otimes(B,tampon1);
    tampon=oplus(tampon1,tampon);
    G1=lfrac(tampon1,ES);

    tampon=otimes(ES,G1); //   E*G
    D1=rfrac(tampon,AES);
    t(0,0)=D1(0,0);
    for (i=1;i<6;i++)
    t(0,0)=inf(t(0,0),D1(i,i));
    D1=Id;
    for (i=0;i<6;i++)
    D1(i,i)=otimes(t(0,0),D1(i,i));

    tampon=lfrac(tampon,B); // B\E*G
    F1=rfrac(tampon,ES);
    F1=prcaus(F1);

  if(G1==G0 && F1==F0 && D1==D0)
  { exit=1;}
  else
  {  n++;
      D0=D1;
      G0=G1;
      F0=F1;
  }
  }
  while(exit!=1 && n<=10 );

 if (n>=10)
 cout<<"Pas bon "<<endl;
 else
 {
cout<<"n"<<n<<endl;
 cout<<"Feedback F="<<F1<<endl;
 cout<<"Matrix   D="<<D1<<endl;
 cout<<"Matrix   G="<<G1<<endl;

 }

 tampon=otimes(D1,A);
 tampon1=otimes(B,F1);
 tampon=oplus(tampon,tampon1);
 tampon=otimes(tampon,ES);
 tampon1=otimes(ES,G1);

 if (tampon==tampon1)
 { cout<<"On a bien égalité "<<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);
 }
}
