clear
/////////////////////////////////////////////////////////////////////////////////////////////
function N=Twin(M)
    N=maxplus(-plustimes(full(M)));
endfunction
/////////////////////////////////////////////////////////////////////////////////////////////

function X=mpm(Z,Y)
// min plus multiplication of two max-plus matrices

[liny,coly]=size(Y);
[linz,colz]=size(Z);
Y=plustimes(Y);Z=-plustimes(Z);

X=maxplus(Z)*maxplus(-Y);
X=plustimes(X);

X=maxplus(-X);


endfunction
/////////////////////////////////////////////////////////////////////////////////////////////

function PsiPlus1=PsiI(alpha,M,A,B,C,D,x)

// x is the precedent psi
// C and D are the matrices for the cone representation of the AB-invariant module
// A and B descriB the system, x(k+1)=Ax(k)+Bu(k), A*mu is the column j of A when mu is a vector 
// made of %0 except at line j.
// One apply this algorithm for each column of the friend F choosing the right mu
// The algorithm stops when PsiPlus1 is equal to precedent PsiPlus1. This function descriBd
// here is one step of the algorithm.

    Bta=mpm(mpm(maxplus(-plustimes(C*B)'),(D*B*x*M+D*A*M)),maxplus(-plustimes(M)')); //mpm is defined above
    Gamma=mpm(mpm(maxplus(-plustimes(D*B)'),(C*B*x*M+C*A*M)),maxplus(-plustimes(M)'));;
    PsiPlus1=maxplus(-plustimes(alpha))+maxplus(-plustimes(Bta))+maxplus(-plustimes(Gamma));
    PsiPlus1=maxplus(-plustimes(PsiPlus1));
    
endfunction
/////////////////////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////////////////////
function B=Normalizep(A)
    //Do not work yet, well, maybe it works now, try it!
    [nlin,ncol]=size(A);
    A=Normalize(A);
    B=%zeros(nlin,ncol);
    for  i=1:ncol
       r=Twin(mpm(%ones(1,nlin),(A(:,i)))); 
       B(:,i)=r*A(:,i); 
 end;

endfunction


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

function B=Normalize(A)
    [nlin,ncol]=size(A);
    B=%zeros(nlin,ncol);
    for  i=1:ncol
       r=mpm(Twin(%ones(1,nlin)*A(:,i))',%1); 
       B(:,i)=r*A(:,i); 
   end;
endfunction

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

A=%zeros(8,8);
A(1,2)=100;
A(1,5)=280;
A(2,2)=115;
A(2,5)=295;
A(3,2)=5;
A(3,5)=240;
A(4,2)=20;
A(4,5)=255;
A(5,4)=0;
A(6,2)=0;
A(7,4)=0;
A(8,7)=0;


disp(full(A),'A = ');
B=%zeros(8,1);
B(2,1)=0;
disp(full(B), 'B = ');

[linB,colB]=size(B); p=colB; n=linB;// To B used further.

Cx=%zeros(8,8);
Cx(8,3)=-250;
Cx(6,1)=-110;
Cx=full(Cx);


K=star(Cx);

disp(full(K),' The specifications are x should remain in Im K, with K = ');
Mj=K;j=0;
[LinM,ColM]=size(Mj);

rep=1;change=2;
//while rep==1 do
for k=1:4 do

Mi=weakbasis(full(Mj));
[Cpt,Dpt]=mpkernel(Mi'); Pi=Cpt';Qi=Dpt';[LinP2,ColP2]=size(Pi);
Hei=[Pi*A Pi*B;Pi %zeros(LinP2,p)];Qei=[Qi*A Qi*B;Qi %zeros(LinP2,p)];
[LinH,ColH]=size(Hei);
Hei=full(Hei);
Qei=full(Qei);
Mj=mpsolve(Hei,Qei);
//  Mj=VariasLineasAxLeqBx(Hei,Qei,%eye(ColH,ColH));
//  Mj=VariasLineasAxLeqBx(Qei,Hei,Mj);
Mj=weakbasis(full(Mj));
Uj=Mj(LinM+1:ColH,:);Mj=Mj(1:LinM,:);

disp(includespan(full(Mj),full(A*Mj)),'Is M A-invariant ?');
disp(includespan(full(Mj),full(A*Mj+B*Uj)),'Is M AB-invariant ?');
j=j+1;
disp(full(Normalizep(weakbasis(full(Mj)))),j,'M =');
Mj=weakbasis(full(Mj));
[LinMj,ColMj]=size(Mj);

//change=x_choose(['Remove one column','Do not change anything'],'Do you want to ... ?');
//if change==1 then 
//      cc=evstr(x_dialog('Give the column numBr to remove :',['1']));
//      disp(cc,'remove column:');
//      Mj(:,cc)=%0;
//      Mj=weakbasis(Mj);
//end //end if
//rep=x_choose(['Continue the sequence','Stop'],'Do you want to ... ?');
end
M=Mj;


/////////////////////////////////////////////////////////////////////////////////////////////
F=Twin(%zeros(colB,linB)); //Twin is defined above, just to have -x

[LinM,ColM]=size(M);
[Pt,Qt]=mpkernel(M');
P=Pt';Q=Qt';
for i=1:ColM
Alpha=#(295)*%ones(p,n);Psi=Alpha;//Alpha=295 and above is ok but 10.99 is not big enougth; Psi is defined above
while (or(PsiI(Alpha,M(:,i),A,B,P,Q,Psi)~=Psi)) do
Psi=PsiI(Alpha,M(:,i),A,B,P,Q,Psi);
end;
F=Twin(Twin(F)+Twin(Psi));
end;
disp(full(F),'F=');
disp(and(P*A*M<=(Q*A*M)+Q*B*F*M),'Is P*A*M<=(Q*A*M)+Q*B*F*M ?');
disp(and(Q*A*M<=(P*A*M)+P*B*F*M),'Is Q*A*M<=(P*A*M)+P*B*F*M ?');
disp(includespan(full(M),full((A+B*F)*M)),'Is M A+BF invariant ?');