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


//serie odot(serie &,serie &);
/////////////////////////////////////////////////////////////////////
/* Example : Duality and interval analysis over idempotent semiring
Thomas Brunsch, Laurent Hardouin, Carlos Andrey Maia, Joerg Raisch
Submited to : Linear Algebra And its Applications

The goal of this programm is to compute with interval analysis the
projection of a degenerated interval X=\underline{X}=\overline{X}
in the solution set of the equation
A\otimes X <= X <= B \odot X

where A and B are matrices of interval whith bounds in Zmax of gamma
Each bound of entries of B are assumed to be monomial

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

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


int main()
{
try
{
    int N=4;
    smatrix A(N,N);
    smatrix B(N,N);
    smatrix A_Star(N,N);
    smatrix B_Star(N,N);
    smatrix Z;

    smatrix X0(N,1);
    smatrix X1(N,1);
smatrix X2(N,1);

    int i,j;
   srand(time(NULL));
  /*for(i=0;i<A_low.getrow();i++)
  {
      for(j=0;j<A_low.getcol();j++)
      {
          B_low(i,j)=gd(_infinity,infinity);
           A_low(i,j)=gd(_infinity,infinity);


      }
*/

for(i=0;i<N;i++)
{
    X0(i,0)=gd(rand()%10,rand()%10);
}

  A(0,2)=gd(0,5);
  A(2,1)=gd(5,0);
  A(1,0)=gd(0,7);
  A(2,3)=gd(0,8);
  A(3,2)=gd(2,1);
  // here it is not fair, I choose B such that a solution exists, but it is to show the efficiency of the algorithm
Z=otimes(A,X0);
B=rfrac(Z,X0);
 /* here I define a B matrix

 B(0,2)=gd(0,5);
  B(2,1)=gd(3,0);
  B(1,0)=gd(0,9);
  B(2,3)=gd(0,12);
  B(3,2)=gd(2,1);
*/




 // A_Star=star(A);
 // B_Star=star(B);

system("pause");
cout<<"A\n"<<A<<endl;
//cout<<"A_Star\n"<<A_Star<<endl;
cout<<"B\n"<<B<<endl;
//out<<"B_Star\n"<<B_Star<<endl;
//A=A_Star;
//B=B_Star;
// I want to solve A_up X= A_low X with X<= X0
for (i=0;i<N;i++) // here I choose x0 equal to e
{
    X0(i,0)=gd(0,0);
}
for(i=0;i<N;i++) // here it is a random X0
{
    X0(i,0)=gd(rand()%10,rand()%10);
}
smatrix Xinit=X0;

// xn+1= B\Axn wedge  A\Bxn wedge xn
X1=otimes(A,X0);
X1=lfrac(X1,B);
X2=otimes(B,X0);
X2=lfrac(X2,A);
X1=inf(X1,X2);
X1=inf(X1,X0);
int iterator=0;
while(!(X0==X1)&& iterator<10)
         { iterator++;
             X0=X1;
             X1=otimes(A,X0);
X1=lfrac(X1,B);
X2=otimes(B,X0);
X2=lfrac(X2,A);
X1=inf(X1,X2);
X1=inf(X1,X0);
cout<<"X1\n"<<X1<<endl; system("pause");
         }
smatrix X3;
cout<<"X1\n"<<X1<<endl;
cout<<"Xinit\n"<<Xinit<<endl;
X2=otimes(A,X1);
X3=otimes(B,X1);
if(X2==X3)
{
     cout<<" Yes it is a solution "<<endl;
}
/*poly p,q;
gd r;
for(i=0;i<N;i++)
 X(i,0).init(p.init(0,20),q.init(0,20),r.init(4,26));

 //X=otimes(X,AdualBupStar);

 //X=prcaus(X);

cout<<"X\n"<<X<<endl;


smatrix aff;


Px_up=lfrac(X,AdualBupStar);

Px_low=lfrac(X, AdualBlowStar);
Px_low=inf(Px_low,Px_up);

cout<<"PX_low\n"<<Px_low<<endl;
cout<<"PX_up\n"<<Px_up<<endl;


if (otimes(A_lowStar,Px_low)==Px_low) cout<<"OK A_lowStar PX == PX"<<endl;
if (otimes(A_upStar,Px_up)==Px_up) cout<<"OK A_upStar PX_up == PX_up"<<endl;

if (odot(B_lowDualStar ,Px_low)==Px_low) cout<<"OK B_lowDualStar odot PX == PX"<<endl;
else
{
    cout<<"B_lowDualStar "<<B_lowDualStar<<endl;
    aff=odot(B_lowDualStar ,Px_low);
     cout<<"B_lowDualStar odot PX_low = "<<aff <<endl;

}

if (odot(B_upDualStar ,Px_up)==Px_up) cout<<"OK B_upDualStar odot PX_up == PX_up"<<endl;
else
{
    cout<<"B_upDualStar "<<B_upDualStar<<endl;
    aff=odot(B_upDualStar ,Px_up);
     cout<<"B_upDualStar odot PX_up = "<<aff <<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);
 }
}


