#ifndef _WIN32
#include "../include/lminmaxgd.h"
#else
#include "..\include\lminmaxgd.h"
using namespace std;
#endif
#include <time.h>


//serie odot(serie &,serie &);
/////////////////////////////////////////////////////////////////////
/*

*/
#include "..\src\gd.cpp"
#include "..\src\poly.cpp"
#include "..\src\serie.cpp"
#include "..\src\smatrix.cpp"
#include "..\src\tools.cpp"

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


int main()
{
try
{
/*
	F = smatrix(4,8);

		F(1,1)=gd(1, 14);
		F(1,2)=gd(1 ,14);
		F(1,3)=gd(1 ,14);
		F(1,4)=gd(1 ,13);
		F(1,5)=gd(1 ,14);
		F(1,6)=gd(1 ,14);
		F(1,7)=gd(1 ,14);
		F(1,8)=gd(1 ,14);
		F(2,1)=gd(1 ,11);
		F(2,2)=gd(1 ,14);
		F(2,3)=gd(1 ,11);
		F(2,4)=gd(1 ,10);
		F(2,5)=gd(1 ,14);
		F(2,6)=gd(1 ,14);
		F(2,7)=gd(1 ,14);
		F(2,8)=gd(1 ,14);
		F(3,1)=gd(1 ,14);
		F(3,2)=gd(1 ,14);
		F(3,3)=gd(1 ,14);
		F(3,4)=gd(1 ,13);
		F(3,5)=gd(1 ,14);
		F(3,6)=gd(1 ,14);
		F(3,7)=gd(1 ,14);
		F(3,8)=gd(1 ,14);
		F(4,1)=gd(1 ,14);
		F(4,2)=gd(1 ,14);
		F(4,3)=gd(1 ,14);
		F(4,4)=gd(1 ,14);
		F(4,5)=gd(1 ,14);
		F(4,6)=gd(1 ,14);
		F(4,7)=gd(1 ,14);
		F(4,8)=gd(1 ,14);







	E = smatrix(8,8);

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

	E(1,6)=gd(0, t2);

	E(2,6)=gd(0, 0);

	E(2,7)=gd(0, t3);

	E(2,8)=gd(0, t4);

	E(3,5)=gd(0, t1);

	E(3,7)=gd(0, t3);

	E(3,8)=gd(0, t4);

	E(4,5)=gd(0, t1);

	E(4,7)=gd(0, t3);

	E(4,8)=gd(0, 0);

	E(5,1)=gd(0, -L);

	E(5,3)=gd(0, (-M)+(-t1));

	E(5,4)=gd(0, (-M)+(-t1));

	E(6,1)=gd(0, (-M)+(-t2));

	E(6,2)=gd(0, -L);

	E(7,2)=gd(0, (-M)+(-t3));

	E(7,3)=gd(0, max((-L),(-M)+(-t3)));

	E(7,4)=gd(0 (-M)+(-t3));

	E(8,2)=gd(0 (-M)+(-t4));

	E(8,3)=gd(0 (-M)+(-t4));

	E(8,4)=gd(0 -L);


*/










    int N=8;
    smatrix A(N,N);
    smatrix E(N,N);
    smatrix B(N,4);

    smatrix DualIdentity(N,N);
    smatrix X0(N,1);


    smatrix X(N,1);


    int i,j;

    int t1 = 14;
    int t2 = 17;
    int t3 = 11;
    int t4 = 9;
    int L = 15;  //15
    int M = 4;

   srand(time(NULL));
  for(i=0;i<E.getrow();i++)
  {
      for(j=0;j<E.getcol();j++)
      {
          E(i,j)=gd(_infinity,infinity);

//          A_low(i,j)=gd(rand()%10,rand()%10);
//          A_up(i,j)=gd(rand()%10,rand()%10);
////          A_up(i,j)=otimes(A_low(i,j),A_up(i,j));
//          A_up(i,j)=oplus(A_low(i,j),A_up(i,j));
//
//          B_up(i,j)=gd(0,rand()%10);
//          B_low(i,j)=gd(0,rand()%10);
//
//          B_low(i,j)=inf(B_low(i,j),B_up(i,j));


      }


  }

smatrix  Fref(4,8);

		Fref(1-1,1-1)=gd(1, 14);
		Fref(1-1,2-1)=gd(1 ,14);
		Fref(1-1,3-1)=gd(1 ,14);
		Fref(1-1,4-1)=gd(1 ,13);
		Fref(1-1,5-1)=gd(1 ,14);
		Fref(1-1,6-1)=gd(1 ,14);
		Fref(1-1,7-1)=gd(1 ,14);
		Fref(1-1,8-1)=gd(1 ,14);
		Fref(2-1,1-1)=gd(1 ,11);
		Fref(2-1,2-1)=gd(1 ,14);
		Fref(2-1,3-1)=gd(1 ,11);
		Fref(2-1,4-1)=gd(1 ,10);
		Fref(2-1,5-1)=gd(1 ,14);
		Fref(2-1,6-1)=gd(1 ,14);
		Fref(2-1,7-1)=gd(1 ,14);
		Fref(2-1,8-1)=gd(1 ,14);
		Fref(3-1,1-1)=gd(1 ,14);
		Fref(3-1,2-1)=gd(1 ,14);
		Fref(3-1,3-1)=gd(1 ,14);
		Fref(3-1,4-1)=gd(1 ,13);
		Fref(3-1,5-1)=gd(1 ,14);
		Fref(3-1,6-1)=gd(1 ,14);
		Fref(3-1,7-1)=gd(1 ,14);
		Fref(3-1,8-1)=gd(1 ,14);
		Fref(4-1,1-1)=gd(1 ,14);
		Fref(4-1,2-1)=gd(1 ,14);
		Fref(4-1,3-1)=gd(1 ,14);
		Fref(4-1,4-1)=gd(1 ,14);
		Fref(4-1,5-1)=gd(1 ,14);
		Fref(4-1,6-1)=gd(1 ,14);
		Fref(4-1,7-1)=gd(1 ,14);
		Fref(4-1,8-1)=gd(1 ,14);

int n=10;

A(0,1)=gd(1,t2);
A(2,2)=gd(1,t3);
A(1,3)=gd(1,t4);
A(2,0)=gd(1,t1);
A(2,2)=gd(1,t3);
A(2,3)=gd(1,t4);
A(3,0)=gd(1,t1);
A(3,2)=gd(1,t3);
A(4,0)=gd(1,0);
A(5,1)=gd(1,0);
A(6,2)=gd(1,0);
A(7,3)=gd(1,0);

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



X0(0,0)=gd(0, 17);
		X0(1,0)=gd(0,15);
		X0(2,0)=gd(0 ,18);
		X0(3,0)=gd(0 ,19);
		X0(4,0)=gd(0 ,4);
		X0(5,0)=gd(0 ,0);
		X0(6,0)=gd(0 ,4);
		X0(7,0)=gd(0 ,5);
smatrix Eref(N,N);
int LM;
Eref(1-1,5-1)=gd(0, 0);
  Eref(1-1,6-1)=gd(0, t2);
  Eref(2-1,6-1)=gd(0 ,0);
  Eref(2-1,7-1)=gd(0 ,t3);
  Eref(2-1,8-1)=gd(0 ,t4);
  Eref(3-1,5-1)=gd(0 ,t1);
  Eref(3-1,7-1)=gd(0 ,t3);
  Eref(3-1,8-1)=gd(0 ,t4);
  Eref(4-1,5-1)=gd(0 ,t1);
  Eref(4-1,7-1)=gd(0 ,t3);
  Eref(4-1,8-1)=gd(0 ,0);
  Eref(5-1,1-1)=gd(0 ,-L);
  Eref(5-1,3-1)=gd(0 ,(-M)+(-t1));
  Eref(5-1,4-1)=gd(0 ,(-M)+(-t1));
  Eref(6-1,1-1)=gd(0 ,(-M)+(-t2));
  Eref(6-1,2-1)=gd(0 ,-L);
  Eref(7-1,2-1)=gd(0, (-M)+(-t3));
  if (-L>=(-M-t3)) LM=-L;
    else LM=-M-t3;
  Eref(7-1,3-1)=gd(0 ,LM);
  Eref(7-1,4-1)=gd(0, (-M)+(-t3));
  Eref(8-1,2-1)=gd(0, (-M)+(-t4));
  Eref(8-1,3-1)=gd(0, (-M)+(-t4));
  Eref(8-1,4-1)=gd(0 ,-L);

/*E(1,5)=gd(0, 0);
  E(1,6)=gd(0, t2);
  E(2,6)=gd(0 ,0);
  E(2,7)=gd(0 ,t3);
  E(2,8)=gd(0 ,t4);
  E(3,5)=gd(0 ,t1);
  E(3,7)=gd(0 ,t3);
  E(3,8)=gd(0 ,t4);
  E(4,5)=gd(0 ,t1);
  E(4,7)=gd(0 ,t3);
  E(4,8)=gd(0 ,0);
  E(5,1)=gd(0 ,-L);
  E(5,3)=gd(0 ,(-M)+(-t1));
  E(5,4)=gd(0 ,(-M)+(-t1));
  E(6,1)=gd(0 ,(-M)+(-t2));
  E(6,2)=gd(0 ,-L);
  E(7,2)=gd(0, (-M)+(-t3));
  if (-L>=(-M-t3)) LM=-L;
    else LM=-M-t3;
  E(7,3)=gd(0 ,LM);
  E(7,4)=gd(0, (-M)+(-t3));
  E(8,2)=gd(0, (-M)+(-t4));
  E(8,3)=gd(0, (-M)+(-t4));
  E(8,4)=gd(0 ,-L);
*/

  E(4,0)=gd(0, 0);
  E(5,0)=gd(0, -t2);
  E(5,1)=gd(0 ,0);
  E(6,1)=gd(0 ,-t3);
  E(7,1)=gd(0 ,-t4);
  E(4,2)=gd(0 ,-t1);
  E(6,2)=gd(0 ,-t3);
  E(7,2)=gd(0 ,-t4);
  E(4,3)=gd(0 ,-t1);
  E(6,3)=gd(0 ,-t3);
  E(7,3)=gd(0 ,0);
  E(0,4)=gd(0 ,L);
  E(2,4)=gd(0,-((-M)+(-t1)));
  E(3,4)=gd(0 ,-((-M)+(-t1)));
  E(0,5)=gd(0 ,-((-M)+(-t2)));
  E(1,5)=gd(0 ,L);
  E(1,6)=gd(0,-( (-M)+(-t3)));

  if (-L>=(-M-t3)) LM=-L;
    else LM=-M-t3;
  E(2,6)=gd(0 ,-LM);
  E(3,6)=gd(0, -((-M)+(-t3)));
  E(1,7)=gd(0, -((-M)+(-t4)));
  E(2,7)=gd(0, -((-M)+(-t4)));
  E(3,7)=gd(0 ,L);

  smatrix AStar=star(A);


cout<<"AStar\n"<<AStar<<endl;
poly pepsilon=epsilon;
poly qe=gd(0,0);
gd r(1,14);
smatrix GDIdentity(N,N);
   for(i=0;i<DualIdentity.getrow();i++)
  {
      for(j=0;j<DualIdentity.getcol();j++)
      {
    if (i==j) {
        DualIdentity(i,j)=gd(0,0);
        GDIdentity(i,j)=serie(pepsilon,qe,r);}
    else DualIdentity(i,j)=gd(_infinity,infinity);


      }
  }



  smatrix Temp;
  smatrix En=E;
  smatrix EdualStar=DualIdentity;
  Temp=inf(En,EdualStar);

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

//smatrix s1(1,1);
//smatrix s2(1,1);
//smatrix AStardualIdentity=Duallfrac(AStar,DualIdentity);
//if (AStardualIdentity==AStar) cout<<" OUFOUF"<<endl;
//else
//{for (i=0;i<DualIdentity.getrow();i++)
//{
//
//
//for (j=0;j<DualIdentity.getcol();j++)
//{ s2(0,0)=AStar(i,j);
//  s1(0,0)=AStardualIdentity(i,j);
//    if(s2==s1) cout<<"ok\n";
//    else
//    {
//
//
//    cout<<"AStar"<<i<<","<<j<<"]"<<s2(0,0)<<endl;
//    cout<<"AStardualidentity["<<i<<","<<j<<"]"<<s1(0,0)<<endl;
//    }
//
//}
//}
//}
//getchar();


smatrix AdualE=Duallfrac(AStar,EdualStar);
cout<<"AdualE\n"<<AdualE<<endl;
//
//for (i=0;i<AdualE.getrow();i++)
//{
//
//
//for (j=0;j<AdualE.getcol();j++)
//{ s2(0,0)=oplus(AdualE(i,j),AStar(i,j));
//  s1(0,0)=AdualE(i,j);
//    if(s2==s1) cout<<"ok\n";
//    else
//    {
//
//
//    cout<<"AdualE["<<i<<","<<j<<"]"<<AdualE(i,j)<<endl;
//    cout<<"AStar["<<i<<","<<j<<"]"<<AStar(i,j)<<endl;
//    }
//
//}
//}
if (AdualE==oplus(AdualE,AStar))
cout<<"condition fulfilled AdualE is greater than AStar"<<endl;
getchar();


smatrix AdualEStar=star(AdualE);
cout<<" AdualEstar : the projection matrix "<<AdualEStar<<endl;
smatrix AB=otimes(AStar,B);
cout <<"A^ast B"<<AB<<endl;

smatrix Xref=otimes(AdualEStar,X0);
 Xref=otimes(GDIdentity,Xref); // a reference model, he has to be greater than otimes(AdualEStar,X0)

//Xopt=AdualEstarB;
//Xopt=AdualEStar;
if( Xref==lfrac(Xref, AdualEStar))
{
cout<<"Xopt is already in the intersection\n";
 getchar();
}
smatrix Xopt=lfrac(Xref, AdualEStar); //projection in the set of solution
cout<<"Xopt \n"<<Xopt;

smatrix tp=otimes(AStar,Xopt);

if( Xopt==tp )
{
 cout<<"Xopt is in imAstar  ";getchar();
 if (Xopt==odot(EdualStar,Xopt)) cout<<"and in imEdualStar \n";getchar();
}

smatrix F=rfrac(Xopt,Xopt);
if( F==otimes(AStar,F) && F==odot(EdualStar,F))
{
 cout<<"Xopt/Xopt is in imAstar and in imEdualStar \n";getchar();
}

if (F==Xopt) {cout<<"Furthermore Xopt==Xopt/Xopt\n";getchar();}

//F=rfrac(Xopt,X0);

//F=lfrac(F,AStar);
F=lfrac(F,B);
//F=rfrac(F,AStar);
F=prcaus(F);
cout <<"F"<<F<<endl;

smatrix BF=otimes(B,F);
smatrix ABF=oplus(A,BF);
smatrix ABFStar=star(ABF);

smatrix ABFStarXopt=otimes(ABFStar,Xopt);

smatrix ABFStarX0=otimes(ABFStar,X0);

smatrix EdualStarABFStarX0=odot(EdualStar,ABFStarX0);
if (ABFStarX0==EdualStarABFStarX0) {cout<<"The state is in the good set (intersection of two images) \n"; getchar();}
else{cout<<" there is a problem ???\n";
//for (i=0;i<ABFStar.getrow();i++)
//{
//smatrix s2(1,1);
//smatrix s1(1,1);
//
//for (j=0;j<ABFStar.getcol();j++)
//{ s2(0,0)=ABFStar(i,j);
//  s1(0,0)=AdualEStar(i,j);
//    if(s2==oplus(s1,s2)) cout<<"ok\n";
//    else
//    {
//
//
//    cout<<"ABFStar(["<<i<<","<<j<<"]"<<ABFStar(i,j)<<endl;
//    cout<<"AdualEstar["<<i<<","<<j<<"]"<<AdualEStar(i,j)<<endl;
//    }
//
//}
//}
}
// Here jsut to check with the Katz Matrix, it is the same thant the previous state,
smatrix Erefstar=star(Eref);
smatrix ErefstarABFStarX0=otimes(Erefstar,ABFStarX0);

if (ABFStarX0==ErefstarABFStarX0)
{cout<<"OK with the Katz point of view"<<endl;
getchar();
}


if(ABFStarXopt==Xopt)
cout<<" For Xopt it is OK\n";
getchar();
smatrix ABFStarB=otimes(ABFStar,B);

ABFStarX0=otimes(ABFStar,X0);
smatrix AdualEStarX0=otimes(AdualEStar,X0);
if (ABFStarX0==oplus(ABFStarX0,AdualEStarX0))
cout<<"CNS Ok (Ouerghi PhD proposition 4.9)\n";
else
{ cout<<" CNS NOK\n";
//    for (i=0;i<ABFStarX0.getrow();i++)
//{
//smatrix s2(1,1);
//smatrix s1(1,1);
//
//for (j=0;j<ABFStarX0.getcol();j++)
//{ s2(0,0)=ABFStarX0(i,j);
//  s1(0,0)=AdualEStarX0(i,j);
//    if(s2==oplus(s2,s1)) cout<<"ok\n";
//    else
//    {
//
//
//    cout<<"ABFStarX0(["<<i<<","<<j<<"]"<<ABFStarX0(i,j)<<endl;
//    cout<<"AdualEstarX0["<<i<<","<<j<<"]"<<AdualEStarX0(i,j)<<endl;
//    }
//
//}
//}


}

smatrix FOuerghi=F; // on sauvegarde pour comapraison
smatrix ABFStarOuerghi=ABFStar;
smatrix ABFStarX0Ouerghi=ABFStarX0;
getchar();
cout<<"Now we compute with the feedback given by Katz";
getchar();
F=Fref;
BF=otimes(B,F);
ABF=oplus(A,BF);
ABFStar=star(ABF);

 ABFStarX0=otimes(ABFStar,X0);
 Erefstar=star(Eref);
ErefstarABFStarX0=otimes(AdualEStar,ABFStarX0);
if (ABFStarX0==ErefstarABFStarX0)
{cout<<"The goal is achieved with the Katz feedback "<<endl;
getchar();
}

smatrix ErefStarABFStar=otimes(AdualEStar,ABFStar);

AdualEStarX0=otimes(AdualEStar,X0);

smatrix s2(1,1);
smatrix s1(1,1);
ABFStarB=otimes(ABFStar,B);
if((ABFStarX0==oplus(ABFStarX0,AdualEStarX0)))
{
cout<<"And CNS Ouerghi prop.4.9 is fulfilled\n";
getchar();
}

//{
// for (i=0;i<ABFStar.getrow();i++)
//{
//
//for (j=0;j<ABFStarX0.getcol();j++)
//{ s2(0,0)=ABFStarX0(i,j);
//  s1(0,0)=AdualEStarX0(i,j);
//    if(s2==oplus(s1,s2)) cout<<"ok \n";
//    else
//    {
//
//
//    cout<<"ABFStarX0(["<<i<<","<<j<<"]"<<ABFStarX0(i,j)<<endl;
//    cout<<"AdualEstarX0["<<i<<","<<j<<"]"<<AdualEStarX0(i,j)<<endl;
//    }
//
//}
//}
//}

cout<<"Comparison of the two feedback"<<endl;
if(FOuerghi==oplus(FOuerghi,Fref)) cout<<"Ouerghi Feedback is greater \n";
else{
 for (i=0;i<F.getrow();i++)
{
for (j=0;j<F.getcol();j++)
{ s2(0,0)=F(i,j);
  s1(0,0)=FOuerghi(i,j);
    if(s1==oplus(s1,s2)) cout<<"Ouerghi greater \n";
    else
    {


    cout<<"FKatz(["<<i<<","<<j<<"]"<<Fref(i,j)<<endl;
    cout<<"FOuerghi["<<i<<","<<j<<"]"<<FOuerghi(i,j)<<endl;
    }

}
}
}

cout<<"Comparison of the two (A+BF)^ast"<<endl;
if(ABFStarOuerghi==oplus(ABFStarOuerghi,ABFStar)) cout<<"Ouerghi transfer is greater \n";
else{
 for (i=0;i<ABFStar.getrow();i++)
{
for (j=0;j<ABFStar.getcol();j++)
{ s2(0,0)=ABFStar(i,j);
  s1(0,0)=ABFStarOuerghi(i,j);
    if(s1==oplus(s1,s2)) cout<<"Ouerghi is greater \n";
    else
    {


    cout<<"ABFStarKatz(["<<i<<","<<j<<"]"<<ABFStar(i,j)<<endl;
    cout<<"ABFStarOuerghi["<<i<<","<<j<<"]"<<ABFStarOuerghi(i,j)<<endl;
    }

}
}
}

cout<<"Comparison of the two (A+BF)^astX0"<<endl;
if(ABFStarX0Ouerghi==oplus(ABFStarX0Ouerghi,ABFStarX0)) cout<<"Ouerghi (A+BF)^astX0 is greater \n";
else{
 for (i=0;i<ABFStarX0.getrow();i++)
{
for (j=0;j<ABFStarX0.getcol();j++)
{ s2(0,0)=ABFStarX0(i,j);
  s1(0,0)=ABFStarX0Ouerghi(i,j);
    if(s1==oplus(s1,s2)) cout<<"Ouerghi is greater \n";
    else
    {


    cout<<"ABFStarX0Katz(["<<i<<","<<j<<"]"<<ABFStarX0(i,j)<<endl;
    cout<<"ABFStarX0Ouerghi["<<i<<","<<j<<"]"<<ABFStarX0Ouerghi(i,j)<<endl;
    }

}
}
}
}


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


