

#ifndef _WIN32
#include "../include/lminmaxgd.h"
#else
#include "..\include\lminmaxgd.h"

#endif
#include <time.h>
using namespace std;
#include <iostream>
#include <time.h>
#include "ominus.cpp"
#include "expandorigin.cpp"
//#include "..\src\gd.cpp"
//#include "..\src\poly.cpp"
//#include "..\src\serie.cpp"
//#include "..\src\smatrix.cpp"
//#include "..\src\tools.cpp"
#include "gd.cpp"
#include "poly.cpp"
#include "serie.cpp"
#include "smatrix.cpp"
#include "tools.cpp"


double inline min(double a,double b)
{
		if(a>b) return b;
		else return a;
}

poly rresid( serie s, long T){
poly temp;
gd m, cr, tempmon;
int crossi=1000;
int crossj=1000;

for (int i=0; i<s.getp().getn(); i++)
{
    m = s.getp().getpol(i);
    if(m.getd()<= T)
    {
        temp = oplus(temp,m);
    }

       if(m.getd() > T)
    {
        crossi= min(i,crossi);
     //   std::cout << " crossi:" << std::endl << crossi<< std::endl;
    }
    if(crossi < 1000)
     {
     cr (s.getp().getpol(crossi).getg(),infinity);
         temp = oplus(temp,cr);
     }
}

for (int j=0; j<=s.getq().getn(); j++)
{
    m = s.getq().getpol(j);
    if(m.getd()<= T)
    {
    temp = oplus(temp,m);
    }
//       if(m.getd() > T)
//    {
//        crossj= min(j,crossj);
//  //      std::cout << " crossj:" << std::endl << crossj<< std::endl;
//    }

}
//std::cout << " rresid(u)" << std::endl << temp<< std::endl;
return temp;


}

poly ru( serie s, long T){
poly temp;
gd m, cr, tempmon;
int crossi=1000;
int crossj=1000;

for (int i=0; i<s.getp().getn(); i++)
{
    m = s.getp().getpol(i);
    if(m.getd()<= T)
    {
        temp = oplus(temp,m);
    }

       if(m.getd() > T)
    {
        crossi= min(i,crossi);
       // std::cout << " crossi:" << std::endl << crossi<< std::endl;
    }
    if(crossi < 1000)
     {
    cr (s.getp().getpol(crossi).getg(),T);
    //tempmon (infinity,infinity);
//cr=oplus(cr,tempmon);
     //  cr (infinity,T);

         temp = oplus(temp,cr);
     }
}

for (int j=0; j<=s.getq().getn(); j++)
{
    m = s.getq().getpol(j);
    if(m.getd()<= T)
    {
    temp = oplus(temp,m);
    }
       if(m.getd() > T)
    {
        crossj= min(j,crossj);
       // std::cout << " crossj:" << std::endl << crossj<< std::endl;
    }

}
// std::cout << " r(u)" << std::endl << temp<< std::endl;
return temp;


}

/////////////////////////////////////////////////////////////////////
/* Example :  Hardouin L., Cottenceau B., Lagrange S., Le Corronc E.

Performance a&nalysis of linear systems over idempotent semiring with additive inputs
Wodes 2008, Goteborg, Sweden
	*/
/////////////////////////////////////////////////////////////////////



int main()
{
    smatrix z1o,z2o, z1p, z2p,y1underline;
    long T=28;
	gd temp_gd, h1, h2;
	poly temp,SR_tf1,SR_tf2,temp1,temp2,x11,x12,x21,x22;
	serie temp_serie,temp_serie2,s1,s,serie,u1opts, u2opts;//,s3,s4;
	smatrix transfunctionH1,transfunctionH2,temp_matrix;
	int time_length;
	int index;
	token_num* token_info;
	token_num* window_info;
	int window_length;
    smatrix A1(2,2),B1(2,1),C1(1,2),AS1(2,2);
    smatrix A2(2,2),B2(2,1),C2(1,2),AS2(2,2);
    smatrix u1_star,x1_star;
    smatrix u2_star,x2_star,u1p_star,u2rcutm;

    smatrix z1(1,1);//,s3(1,1), s4(1,1);
    smatrix z2(1,1);
/////////////////////////////////////////////////////////////

    SR_tf1.init(1,1);
//////////
    h1.init(0,3);
    h2.init(0,4);


    ///* Test of Laurent ru and rresidu */
    temp1.init(0,10)(10,11)(20,30)(30,45)(6,infinity);
    temp.init(40,10)(61,11)(62,30)(73,45)(6,infinity);
    z1o(0,0).init(temp1,temp,e);
    temp1=ru(z1o(0,0),30);
    temp2=rresid(z1o(0,0),30);
    cout<<"z1o"<<z1o<<endl;
    cout<<"temp1"<<z1o<<endl;
    cout<<"temp2"<<z1o<<endl;
///* End test laurent */






     temp.init(0,10)(1,11)(2,30)(3,45)(6,infinity);
     z1(0,0).init(epsilon,temp,e);

////    rresid( z1(0,0), 25);
////    ru( z1(0,0), 25);
    temp.init(0,32)(1,38)(2,44)(5,infinity);
    z2(0,0).init(epsilon,temp,e);

    A1(0,1).init(temp.init(1,1),epsilon,e);
    A1(1,0).init(temp.init(0,3),epsilon,e);


    A2(0,1).init(temp.init(1,1),epsilon,e);
    A2(1,0).init(temp.init(0,4),epsilon,e);



      B1(0,0).init(epsilon,e,e);
      B2(0,0).init(epsilon,e,e);

     C1(0,1).init(epsilon,e,e);
    C2(0,1).init(epsilon,e,e);
////////////////////////////////////////////////////////////////
     AS1 = star(A1);
     AS2 = star(A2);
    temp_matrix = otimes(C1,AS1);
    transfunctionH1 = otimes(temp_matrix,B1);
    std::cout << " H1 :" << std::endl << transfunctionH1<< std::endl;
        u1_star = lfrac(z1,transfunctionH1);
      //  u1p_star = lfrac(z1p,transfunctionH1);
    temp_matrix = otimes(C2,AS2);
    transfunctionH2 = otimes(temp_matrix,B2);
        std::cout << " H2 :" << std::endl << transfunctionH2<< std::endl;

    u2_star = lfrac(z2,transfunctionH2);
    smatrix ASB1;
    ASB1 = otimes(AS1,B1);
    x1_star = otimes(ASB1,u1_star);
        std::cout<<"u1* : "<<std::endl<<u1_star<<std::endl;
        std::cout<<"x11* : "<<std::endl<<x1_star(0,0)<<std::endl;


    smatrix ASB2;
    ASB2 = otimes(AS2,B2);
    x2_star = otimes(ASB2,u2_star);
    std::cout<<"x2 star : "<<std::endl<<x2_star<<std::endl;
//
/////////////////////////////////////////////////////////////////
//
 poly z1poly,z2poly,z1pp;
 poly x21starpoly,x22starpoly,x21poly,x22poly;
 poly x11poly, x12poly,x11starpoly,x12starpoly;
 poly betah1,b3,G2;



//
// //   x23poly=frac(z2poly,H2);

// transformation as poly in order to use odot , odotfrac and odotflta
  x21starpoly=oplus(x2_star(0,0).getp(),x2_star(0,0).getq());
   x22starpoly=oplus(x2_star(1,0).getp(),x2_star(1,0).getq());
 x21poly=x21starpoly;
  x22poly=x22starpoly;
//
//
      x11starpoly=oplus(x1_star(0,0).getp(),x1_star(0,0).getq());
      x12starpoly=oplus(x1_star(1,0).getp(),x1_star(1,0).getq());
      x11poly=x11starpoly;
       x12poly=x12starpoly;


/** HERE COMPUTATION OF  optimal control for SYSTEM 2 **/
int i=0;
//
int step=1;
//
//  //   x4=z2;
//     //
//
      do{
   i++;

//cout << "i"<< i <<endl;

betah1= x21poly; // prevoius value of x21
b3=odot(x21poly,x11starpoly);
b3=frac(b3,SR_tf1); // Sr_tf1 is beta_1 actually

b3=fracodotsharp(b3,x12starpoly);

//
G2=oplus(A2(1,0).getp(),A2(1,0).getq());

b3=frac(b3,G2);


betah1=inf(x21poly,b3);
//cout << "betah1"<< betah1 <<endl;
if(!(x21poly==betah1)){
x21poly=betah1;
step ++;
}
else
step=21;


} while ( step <20);

if (step==21) cout<<" x21* "<<x21poly<<endl;
    else cout<<"the optimal x21* doesn't exist"<<endl;
getchar();
/** HERE COMPUTATION OF  optimal state and output for the two systems **/


u1_star(0,0).init(x11poly,x11poly,e);
u2_star(0,0).init(x21poly,x21poly,e);
x1_star = otimes(ASB1,u1_star);
x2_star = otimes(ASB2,u2_star);
cout<<"x1_star"<<x1_star<<endl;
cout<<"x2_star"<<x2_star<<endl;
getchar();

    x11poly=oplus(x1_star(0,0).getp(),x1_star(0,0).getq());
    x12poly=oplus(x1_star(1,0).getp(),x1_star(1,0).getq());
    x21poly=oplus(x2_star(0,0).getp(),x2_star(0,0).getq());
    x22poly=oplus(x2_star(1,0).getp(),x2_star(1,0).getq());


    std::cout<<"SR_tf1  "<<SR_tf1<<std::endl;
    poly p1,p2,p3;
     //
     p2=x21poly;
   do{
   p3=p2;
   G2=oplus(A2(1,0).getp(),A2(1,0).getq());
       //s.init(p3,p3,e);
//s1=otimes(AS2(1,0),s);
x22poly=otimes(G2,p2);

//x22poly=oplus( s1.getp(),s1.getq());
p1=odot(x12poly,x22poly);
p1=otimes(SR_tf1,p1);



p2=fracodotflat(p1,x11poly); // is only working for polynomial until now
cout<<"p2"<<p2<<endl;
p2=oplus(p2,x21poly);
cout<<"p2"<<p2<<endl;


getchar();
   }
while(!(p2==p3));
u2_star(0,0).init(p3,p3,e);
x2_star=otimes(ASB2,u2_star);

cout<<" u2star "<<u2_star(0,0)<<endl;
cout<<" ASb2"<<ASB2<<endl;

if(p3==x21poly) cout<<"ok"<<endl;





 std::cout<<"X1 : "<<std::endl<<x1_star<<std::endl;
 std::cout<<"X2 : "<<std::endl<<x2_star<<std::endl;


T=18; /// change of the reference input of z1, for example.

///* Computation of the smallest output y1 and y2 we can achieve *///
/*
poly u1_cut,u1_min,u2_cut,u2_min,u1_new;
smatrix u2(1,1);
u1_cut=rresid(x11poly,T);
u1_min=ru(x11poly,T);
u2_cut=rresid(x21poly,T);
u2_min=ru(x21poly,T);
u2(0,0).init(u2_cut,u2_cut,e);
x22=otimes(ASB2,u2);
u1_new=u1_min;

 do{
       temp=u1_new;
temp1=otimes(AS1(1,0),temp);
//cout<<"temp1"<<temp1<<endl;
temp1=odot(x22poly,temp1);

temp1=otimes(beta(0,0),temp1);
cout<<"temp1"<<temp1<<endl;
p1=oplus(temp1.getp(),temp1.getq());
p2=oplus(x1_1.getp(),x1_1.getq());
p1=fracodotflat(p1,p2); // is only working for polynomial until now
temp1.init(p1,p1,e);
cout<<"temp1"<<temp1<<endl;
x1_2=oplus(temp1,su2(0,0));

cout<<" x1_2 "<<x1_2;
getchar();
   }
while(!(temp==x1_2));

cout<<" temp "<<temp;

*/

////poly u2star;
////u2star=frac(y2star,H2);
//
//if (step==22)
//{
//     std::cout<<"x11* before changing the ref z is ok"<<std::endl<<x11poly<<std::endl;
//
//}
//
//u1opt=x11poly;
//u2opt=x21poly;
//
//i=0;
//step=1;
//smatrix u1scutm;
//poly u1p,u2p;
//u1opts.init(epsilon,u1opt,e);
//u2opts.init(epsilon,u2opt,e);
//u1scut=ru(u1opts,T);
//u2rcut=rresid(u2opts,T);
//u1rcut=rresid(u1opts,T);
//u2rcutm(0,0).init(epsilon,u1rcut,e);
////temp_matrix=otimes(u1scutm,transfunctionH1);
//
//temp_matrix=otimes(AS2,B2);
//temp_matrix= otimes(temp_matrix,u2rcutm);
//
//x21=oplus(temp_matrix(0,0).getp(),temp_matrix(0,0).getq());
//x22=oplus(temp_matrix(1,0).getp(),temp_matrix(1,0).getq());
//x11=x11poly;
//i=0;
////x11.init(0,0)(1,5)(2,7)(3,10)(4,12)(5,15)(6,infinity);
////x12=otimes(h1,x11);
//step=1;
////x21=oplus(x2_star(0,0).getp(),x2_star(0,0).getq());
// //x21=epsilon; //*****
//
//     //
//
//     do{
//   i++;
//
//
//b3= x11;// prevoius value of x21
//temp1=otimes(h1,x11);
////temp2=otimes(h1,x11);
////cout << " H1 . x11"<< temp2 <<endl;
//temp=odot_mark2(temp1,x22);
//temp=otimes(SR_tf1,temp);
//temp=odotflat(temp,x21);
//temp=oplus(temp,u1scut);
//
//cout << "temp"<< temp <<endl;
//
//if(!(temp==b3)){
//x11=temp;
//}
//else
//step=21;
//
//step ++;
//} while ( step <20);
//
//if (step==22)
//{
//     std::cout<<"x11 is ok"<<std::endl<<x11<<std::endl;
//          std::cout<<"num of itterations"<<std::endl<<i<<std::endl;
//
//
//}
//x12=otimes(h1,x11);
//
//y1underline(0,0).init(epsilon,x12,e);
//
//
//
//
//
//
//
//temp_matrix=oplus(y1underline,z1p);
//
//temp_matrix=lfrac(temp_matrix,transfunctionH1);
//u1p=oplus(temp_matrix(0,0).getp(),temp_matrix(0,0).getq());
////
////
////x21cut=rresid(x2_star(0,0),T);
////// u1_star = lfrac(z1,transfunctionH1);
////
////
//
//
//
//   do{
//   i++;
//
//
//betah1= u1p; // prevoius value of x22
//b3=odot(u2rcut,u1p);
//b3=frac(b3,SR_tf1);
//temp=otimes(h2,u2rcut);
//b3=fracodotsharp(b3,temp);
////cout << "b3"<< b3 <<endl;
//
//b3=frac(b3,G1);
//betah1=inf(betah1,b3);
//betah1=inf(betah1,u1rcut);
////betah1=oplus(betah1,u1scut);
//if(!(u1p==betah1)){
//u1p=betah1;
//}
//else
//step=41;
//
//step ++;
//} while ( step <40);
//if (step==42)
//{
//     std::cout<<"x11 after change is ok"<<std::endl<<u1p<<std::endl;
//
//}
//


}
