#ifndef _WIN32
#include "../include/lminmaxgd.h"
#else
#include "..\include\lminmaxgd.h"
using namespace std;
#endif
/////////////////////////////////////////////////////////////////////
/* Example : transport system paper ZeCury
Maia C.A., Hardouin L., Cury J.E.R.
"Some results on the feedback control of max-plus linear systems under state constrains",
Conference on Decision and Control, CDC 2013, Firenze.
*/
#include "..\src\gd.cpp"
#include "..\src\poly.cpp"
#include "..\src\serie.cpp"
#include "..\src\smatrix.cpp"
#include "..\src\tools.cpp"

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


int main()
{
try
{
int i,j;
	smatrix A(4,4);	// matrice d'état


  smatrix C(2,4);
  smatrix B(4,3);


	A(0,2)=gd(1,1);//
	A(1,3)=gd(2,3);
	A(2,0)=gd(0,2);
	A(2,3)=gd(0,1);
	A(3,1)=gd(0,3);


	C(0,2)=gd(0,3);
	C(1,3)=gd(0,1);

	B(0,0)=gd(0,1);
	B(0,1)=gd(0,5);

    B(1,1)=gd(0,7);
    B(1,2)=gd(0,5);


   smatrix AStar=star(A);


cout<<"AStar"<<AStar; getchar();
smatrix DualIdentity(4,4);
smatrix E(4,4);  // The matrix of constraint here
   for(i=0;i<DualIdentity.getrow();i++)
  {
      for(j=0;j<DualIdentity.getcol();j++)
      {
    if (i==j) DualIdentity(i,j)=e; // identity
    else DualIdentity(i,j)=Top; // top element (-oo, +oo )

    E(i,j)=Top;

      }
  }

// contraint definition,
E(2,0)=gd(0,4); //  no more than 4 time units between x1 and x3
E(2,1)=gd(0,6); // no more than 6 time units between x2 and x3


smatrix EdualStar=DualIdentity;
smatrix En=E;
  smatrix Temp=inf(En,EdualStar);
  cout<<"Temp"<<Temp<<endl;
cout<<"E_DualStar\n"<<EdualStar<<endl;
  while(!(EdualStar==Temp))
         {
             EdualStar=Temp;
             En=odot(En,E);
             Temp=inf(Temp,En);
         }
//cout<<"E\n"<<E<<endl;
cout<<"E_DualStar\n"<<EdualStar<<endl;
getchar();


smatrix AdualE=Duallfrac(AStar,EdualStar);
cout<<" AdualE \n"<<AdualE<<endl;
getchar();

smatrix AdualEStar=star(AdualE);
cout<<" AdualEstar \n"<<AdualEStar<<endl;



smatrix AB=otimes(AStar,B);
cout <<"AB \n"<<AB<<endl;
smatrix CAB=otimes(C,AB);
cout <<"CAB \n"<<CAB<<endl;
smatrix CA=otimes(C,AStar);
cout <<"CA \n"<<CA<<endl;
smatrix AdualEStarB=otimes(AdualEStar,B);
cout <<"AdualEStarB \n"<<AdualEStarB<<endl;getchar();

smatrix CAdualEStarB=otimes(C,AdualEStarB);
smatrix Gref=CAdualEStarB;





smatrix P;
smatrix test1,test2,x1,x2,x3,x4,x;
int n=0;

P=lfrac(Gref,CAB);
x4=P;
do
{
 x=x4;

 test1=otimes(AdualEStarB,x);
 test2=otimes(AB,x);
 x1= lfrac(test1,AB);
 x2= lfrac(test2,AdualEStarB);

 x2=inf(x2,x1);
 x3=prcaus(x);

 x4=inf(x3,x2);




 cout<<"n"<<++n<<endl;

 }
while((!(x==x4)) &&n<10);

if (n<10)
{
     cout<<" it is ok \n"<<endl;
}
else
{
     cout<<" it is Nok \n"<<endl;
}

P=x;
cout <<"P"<<P<<endl;getchar();
/*
P(0,0)=gd(3,9);
P(1,0)=gd(4,2);
P(2,0)=gd(3,1);

P(0,1)=gd(0,8);
P(1,1)=gd(1,1);
gd e=gd(0,0);
P(1,1)=oplus(P(1,1),e);

gd r;
r=gd(0,2);
P(2,1)=gd(1,3);
P(2,1)=oplus(P(2,1),r);

P(0,2)=gd(0,6);
r=gd(2,4);
P(1,2)=gd(3,5);
P(1,2)=oplus(P(1,2),r);

P(2,2)=P(1,1);
cout <<"P"<<P<<endl;getchar();
*/
gd r;
// on verifie que (A+BF)*BP=E*A*BP


    smatrix ABP=otimes(AB,P);
    smatrix CABP=otimes(C,ABP);

    cout<<"ABP"<<ABP<<endl;


    smatrix test3=otimes(AdualEStar,ABP);
    if (test3==ABP)
        cout<<" P est causal et ABP in \in E* \n P= \n"<<P<<endl;
        else cout<<" NOK  P est non causal et ABP in \in E* "<<endl;


    getchar();
    smatrix F;
    F=lfrac(CABP,CABP);
    F=rfrac(F,ABP);
    F=prcaus(F);
    cout <<"F \n"<<F<<endl;
    F(0,0)=epsilon;
    F(1,0)=epsilon;
    F(2,0)=epsilon;

    F(0,1)=epsilon;
    F(1,1)=epsilon;
    F(2,1)=epsilon;

    F(0,2)=gd(1,0);
    F(1,2)=gd(5,2);
    F(2,2)=gd(4,1);

    F(0,3)=gd(1,1);
    F(1,3)=gd(4,2);
    r=gd(5,3);
    F(1,3)=oplus(F(1,3),r);

    F(2,3)=gd(4,4);
    r=gd(5,5);
    F(2,3)=oplus(F(2,3),r);
cout <<"F \n"<<F<<endl; getchar();
    smatrix BP=otimes(B,P);
    smatrix BPF=otimes(BP,F);
    smatrix ABF=oplus(A,BPF);

    smatrix ABFS=star(ABF);
    smatrix CABFS=otimes(C,ABFS);
    smatrix CABFSBP=otimes(CABFS,BP);
    smatrix ABFSBP=otimes(ABFS,BP);
    if(ABFSBP==ABP)
    cout<<"F est causal "<<F<<endl;
    getchar();
    getchar();

	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);
 }
}
