#ifndef _WIN32
#include "../include/lminmaxgd.h"
#else
#include "..\include\lminmaxgd.h"
using namespace std;
#endif
/////////////////////////////////////////////////////////////////////
/* Example : transport avec ZeCury
*/
#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
    int t1=10;
    int t2=10;
    int t3=10;
    int T2= 15;
    int g1=5;

    int r1=6;
    int R1= 14;

    int g2=5;
    int r2=7; //bien
    //int r2=15;
    int R2=10;

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


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

    A(2,0)=gd(0,r1);
	A(3,1)=gd(0,r2);


	C(0,1)=gd(0,t3);
	B(0,2)=gd(0,t1);
B(2,0)=gd(0,0);
B(3,1)=gd(0,0);
smatrix AStar=star(A);
cout<<"AStar"<<AStar; getchar();
smatrix DualIdentity(4,4);
smatrix E(4,4);
   for(i=0;i<DualIdentity.getrow();i++)
  {
      for(j=0;j<DualIdentity.getcol();j++)
      {
    if (i==j) DualIdentity(i,j)=e;
    else DualIdentity(i,j)=Top;
 E(i,j)=Top;

      }
  }


E(1,0)=gd(0,T2);
E(2,0)=gd(0,R1);
E(3,1)=gd(0,R2);
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"<<AdualE<<endl;

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

smatrix AB=otimes(AStar,B);
cout <<"AB"<<AB<<endl;
smatrix AdualEStarB=otimes(AdualEStar,B);
cout <<"AdualEStarB"<<AdualEStarB<<endl;getchar();

smatrix Parallel(4,4);
poly p,q;
gd r;
p.init(infinity,_infinity);
q.init(0,0);
r.init(0,0);
 Parallel(0,0)=serie(p,q,r);
 Parallel(1,1)=serie(p,q,r);
 Parallel(2,2)=serie(p,q,r);
 Parallel(3,3)=serie(p,q,r);

smatrix X0=otimes(Parallel,AdualEStarB);
X0=otimes(AdualEStar,X0);

//smatrix X0=otimes(Parallel,AdualEStar);
cout<<"X0"<<X0<<endl;getchar();
 smatrix Gref=otimes(C,X0);
 //Gref=otimes(Gref,B);
 X0=lfrac(Gref,C);
 //X0=otimes(X0,)

smatrix Xopt=lfrac(X0, AdualEStar);

cout<<"Xopt"<<Xopt<<endl;getchar();
smatrix F=rfrac(Xopt,Xopt);
//cout <<"F"<<F<<endl;
F=lfrac(F,AB);
cout <<"F"<<F<<endl;getchar();
// on verifie que (A+BF)*X0=X0 avec X0 \in im AdualEStar ou bien \in im ECARLOSStar
    smatrix ECARLOS(4,4);
    ECARLOS(0,1)=gd(0,-T2);
    ECARLOS(0,2)=gd(0,-R1);
    ECARLOS(1,0)=gd(0,t2);
    ECARLOS(1,3)=gd(0,-R2);
    ECARLOS(2,0)=gd(0,r1);
    ECARLOS(3,1)=gd(0,r2);

    smatrix ECARLOSStar=star(ECARLOS);
cout<<"ECarlosStar"<<ECARLOSStar;
cout<<"EdualStar"<<EdualStar<<endl;getchar();
    smatrix ABF=otimes(B,F);
    ABF=otimes(AStar,ABF);
    cout<<"ABF"<<ABF<<endl;
    smatrix ABFS=star(ABF);
//     ABFS=otimes(ABFS,AStar);
    cout<<"ABFStar "<<ABFS<<endl;getchar();
// une Ci compatible
smatrix alpha(4,1);
q.init(0,0);
r.init(0,0);
 alpha(0,0)=serie(p,q,r);
 alpha(1,0)=serie(p,q,r);
 alpha(2,0)=serie(p,q,r);
alpha(3,0)=serie(p,q,r);
alpha=otimes(AdualEStar,alpha);
cout<<"alpha in imE*"<<alpha<<endl;
r.init(0,0);
q.init(0,0)(1,12)(2,24)(3,36)(4,infinity);
 alpha(0,0)=serie(p,q,r);
 q.init(0,10)(1,22)(2,34)(3,46)(4,infinity);
 alpha(1,0)=serie(p,q,r);
 q.init(0,7)(1,19)(2,31)(3,42)(4,infinity);

 alpha(2,0)=serie(p,q,r);
 q.init(0,17)(1,29)(2,41)(3,53)(4,infinity);
alpha(3,0)=serie(p,q,r);
smatrix  test=otimes(AdualEStar,alpha);
if (test==alpha) cout<<"Yes alpha in im E* !!!!\n"<<endl;
else cout<<"No, in im E*"<<endl;
getchar();
    smatrix ABFSB=otimes(ABFS,AdualEStar);
    ABFSB=otimes(ABFS,alpha); // un CI particulière
    smatrix ASCB=otimes(AdualEStar,ABFSB);
test=oplus(ASCB,ABFSB);
if (ASCB==ABFSB) cout<<"Ok Dynamic Feedback Non Causal"<<endl;
else cout<<"Nok feedback Dynamic Non Causal"<<endl;
cout<<"alpha"<<alpha;
cout<<"(A+BF)alpha"<<ABFSB;
getchar();


/** Il faut extraire les C.I compatibles **/


q.init(3,20);
r.init(0,0);
 Parallel(0,0)=serie(p,q,r);
 Parallel(1,1)=serie(p,q,r);
 Parallel(2,2)=serie(p,q,r);
Parallel(3,3)=serie(p,q,r);
smatrix Fc=otimes(F,Parallel);
Fc=prcaus(F);
    smatrix ABFc=otimes(B,Fc);
    ABFc=otimes(AStar,ABFc);
    cout<<"ABFc"<<ABFc<<endl;
    smatrix ABFcS=star(ABFc);
     //ABFcS=otimes(ABFcS,AStar);
    cout<<"ABFStar "<<ABFcS<<endl;getchar();


    smatrix ABFcSB=otimes(ABFcS,AdualEStar);
        ABFcSB=otimes(ABFcS,alpha); // un CI particulière
    smatrix ASCBc=otimes(AdualEStar,ABFcSB);
      test=oplus(ASCBc,ABFcSB);
if (ASCBc==ABFcSB) cout<<"Ok Dynamic Feedback  Causal"<<endl;
else cout<<"Nok feedback Dynamic Causal"<<endl;
cout<<"(A+BF)alpha"<<ABFcSB;getchar();


/**
smatrix P=rfrac(ABFSB,ABFSB);
cout<<"P"<<P<<endl;getchar();
q.init(0,20);
r.init(0,0);
 Parallel(0,0)=serie(p,q,r);
 Parallel(1,1)=serie(p,q,r);
 Parallel(2,2)=serie(p,q,r);
Parallel(3,3)=serie(p,q,r);

P=otimes(P,Parallel);
cout<<"P"<<P<<endl;getchar();
Gref=otimes(ABFSB,P);
Xopt=rfrac(Gref,Gref);
F=lfrac(Xopt,B);

cout<<"F"<<F;getchar();
****/
/******* METHODE CARLOS ***************************/
smatrix FCarlos(3,4);
for (i=0;i<3;i++)
for(j=0;j<4;j++)
 {
     FCarlos(0,j)=gd(1,21);
FCarlos(1,j)=gd(1,27);
FCarlos(2,j)=gd(1,0);

 }
 // un autre
 FCarlos(0,0)=gd(3,29);
 FCarlos(0,1)=gd(3,24);
 FCarlos(0,2)=gd(3,35);
 FCarlos(0,3)=gd(3,25);
 FCarlos(1,0)=gd(3,35);
 FCarlos(1,1)=gd(3,30);
 FCarlos(1,2)=gd(3,41);
 FCarlos(1,3)=gd(3,31);
 FCarlos(2,0)=gd(3,8);
 FCarlos(2,1)=gd(3,3);
 FCarlos(2,2)=gd(3,14);
 FCarlos(2,3)=gd(3,4);
//F=prcaus(F);
     cout<<"FCarlos"<<FCarlos<<endl;
    smatrix ABFCarlos=otimes(B,FCarlos);
    ABFCarlos=oplus(A,ABFCarlos);
    smatrix ABFSCarlos=star(ABFCarlos);
    smatrix ESCARLOS=star(ECARLOS);

    smatrix ABFSBCarlos=otimes(ABFSCarlos,AdualEStar);
    smatrix ASCBCarlos=otimes(AdualEStar,ABFSBCarlos);

if (ASCBCarlos==ABFSBCarlos) cout<<"Ok feedback Static"<<endl;

else cout<<"Nok feedback Static"<<endl; getchar();




ABFSBCarlos=otimes(ABFSCarlos,B);
cout<<"(A+BF)BCARLOS"<<ABFSBCarlos;
smatrix X=lfrac(ABFSBCarlos,AB);
X=rfrac(X,X);
smatrix F2=rfrac(X,AB);
F2=prcaus(F2);
cout<<"F2"<<F2;

smatrix ABF2=otimes(B,F2);
ABF2=oplus(A,ABF2);
smatrix ABF2S=star(ABF2);

 ESCARLOS=star(ECARLOS);
smatrix ABF2SB=otimes(ABF2S,AdualEStar);
smatrix ASCB2=otimes(AdualEStar,ABF2SB);

cout<<"F2"<<F2;
cout<<"FCarlos"<<FCarlos;getchar();
if (ASCB2==ABF2SB) cout<<"Ok feedback F2"<<endl;

else cout<<"Nok feedback F2"<<endl; 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);
 }
}
