
function [sol]=  maximisation2(A1,B1,c1,d1,xfinal1,epsilon1,decimal1)

disp('The program now finds an optimal solution to the problem.')
disp('...Press any key to continue...')
pause

epsilon = epsilon1;
decimal = decimal1;
epsilon=(10^(-decimal))*(fix(epsilon*10^(decimal)));

AA = A1;
BB = B1;

[rAA,cAA]=size(AA);

for i = 1 : rAA
    for j = 1 : cAA - 1
        A(i,j) = AA(i,j);
    end
end
[rBB,cBB]=size(BB);
for i = 1 : rBB
    for j = 1 : cBB - 1
        B(i,j) = BB(i,j);
    end
end
c = AA(:,cAA);
d = BB(:,cBB);

[m,n] = size(A);  % m,n are the no. of rows and columns of matrix A respectively
[m,n] = size(B);  % m,k are the no. of rows and columns of matrix B respectively

xnot=xfinal1;
disp('The program will now check for finiteness (or otherwise) of an optimal solution to the problem.')
disp('... Press any key to continue ...')
pause
%%%%%%%%%%Checking for finiteness (or otherwise) of an optimal solution.%%%%%%%%%%%%%%%%
A;
B;
I = rand(m,n);
for k = 1 : m
    for l = 1 : n
        if k == l
            I(k,l) = 0;
        else
            I(k,l) = -inf;
        end
    end
end
AC = [A;B];
BC = [I;I];
[rAC,cAC] = size(AC);
[rBC,cBC] = size(BC);
x=zeros(1,cAC);
sAC = -AC';          % Negation and transposition of AConvert
sBC = -BC' ;         % Negation and transposition of BConvert
r = 0;             % A non-negative integer and act as the condition of the (while) loop
step = 1;          % count the number of step and act as the condition of the (while)loop
xn=[];            % Any x generated in the program depending on r, called x new
yo=[];            % y generated at the first step with r=0, called y old
yn=[];            % Any y generated in the program depending on r,called y new

xo=x;              % To keep the initial vector x, as xo


while r ~= -1
    
    ACx = [];      % Maximum of ax1 for each i
    ACx1 = [];     % Sum of ith row of matrix A and vector x
    BCy = [];      % Maximum of by1 for each i
    BCy1 = [];     % Sum of ith row of matrix B and vector y
    kx = [];
    ky = [];
    x1 = [];      % Sum of ith row of matrix sA and Ax 
    w = [];
    wy = [];
    y1 = [];      % Sum of ith row of matrix sB and Bx
    y = [];       % minimum of y1 for each i
    z = [];
    zy = [];
            
    for i = 1 : rAC
        ACx1 = AC(i,:) + x;
        ACx(i) = max(ACx1);
    end
    
    ACx;
    
    for i = 1 : cBC
        y1 = sBC(i,:) + ACx;
        y(i) = min(y1);
    end 
    
    y;
    
    yn=y;
    
    for i = 1 : rBC       
        BCy1 = BC(i,:) + y;
        BCy(i) = max(BCy1);
    end
    
    BCy;
    
    if isequal(ACx,BCy)
        disp('The optimal solution is unbounded from above. Therefore the optimal value is +inf.');
        sol=0;
        return
        
    else
        
        for i = 1 : cAC
            x1 = sAC(i,:) + BCy;
            x(i) = min(x1);
        end
        
        x; 
        
        xn=x;
        
        r = r + 1;
        
        if step == 1
            
            yo=y;
        end
        
        yo;
        
        for i = 1 : cAC
            if xn(i) < xo(i)
                t = 1;
            else
                t = 0;
            end
            w = union(t,kx);
            z = union(z,w);
        end
        
        if all(z)
            disp('The problem has a finite optimal solution.')
            disp('...Press any key to find the solution...')
            pause
            r=-1;
        end
        
        if r~=-1
            for i = 1 : cBC
                if yn(i) < yo(i)
                    ty = 1;
                else
                    ty = 0;
                end
                wy = union(ty,ky);
                zy = union(zy,wy);
            end
        
            if all(zy)
                disp('The problem has a finite optimal solution.')
                disp('...Press any key to find the solution...')
                pause
                r=-1;
            end
        end
        
        r;
        step = step+1;
        
    end
end
%%%End of checking for finiteness of optimal solution to the problem %%%%%%%%

disp('Do you want to enter a vector h of coefficients of the objective function? ')
disp('The answer ''no'' means the coefficients will remain zero.');
vopt = input('Press 1 to enter, 0 otherwise: ');
if vopt
    dec = input('Press 1 to enter manually, 0 by random: ');
    if dec
        disp('Note: The vector h of the objective function coefficients must have the same dimension as x0.');
        disp('... Press any key to continue ...');
        pause;
        h = input('Enter the vector h:');
        disp(h);
    else
        lambda = input('Enter (Positive integer) range for elements of the vector (e.g 5 for the range [-5 5]): ');
        rh = rand(1,n);
        for j = 1 : n
            c2 = rand();
            if c2 < 0.5
                s = -1;
            else
                s = 1;
            end
            rh(j) = s * round(rh(j) * lambda);
        end
        h = rh        
    end
    
    y = xfinal1+h;
    diagh = rand(n);
    diaghinv = rand(n);
    
    for i = 1:n
        for j = 1:n
            if i==j
                diagh(i,j) = h(j);
                diaghinv(i,j) = -h(j);
            else
                diagh(i,j) = -inf;
                diaghinv(i,j) = -inf;
            end
        end
    end
    diagh;
    diaghinv;
    for i = 1:m
        for j = 1:n
            A(i,j)=max(A(i,:)+diaghinv(:,j)');
            B(i,j)=max(B(i,:)+diaghinv(:,j)');
            xnot=y;
        end
    end
    disp('Since the vector of coefficients of the objective function is not zero, ')
    disp('the problem has been transformed to the same problem for the following matrices and vectors. ');
    disp('Note: The sides of the equations have been swapped (if necessary) to make cc >= dd');
    AA=A
    BB=B
    cc=c
    dd=d
    disp('The objective fuction value at x0 is :');
    disp(max(xnot));
else
    h=0;  
end
disp('... Press any key to continue ...');
pause
count=1;
for i=1:m
    AN(i,:)=c(i)-A(i,:);
end
for i=1:m
    ANmax(count)=max(AN(i,:));
    count = count + 1;
end
U=max(ANmax);
L=max(xnot);
m=U;
xhat =  xnot;
D = epsilon + 1;

while D > epsilon
    z=m;
    Zero=zeros(1,n);
    One=ones(1,n);
    mOne=-One;
    P=[A c;Zero z-1];
    Q=[B d;mOne z];
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    S =P; 
    T = Q;
    [rS,cS] = size(S);  % rS,cS are the no. of rows and columns of matrix S respectively
    [rT,cT] = size(T);  % rT,cT are the no. of rows and columns of matrix T respectively
    I = rand(rS,rS);
    
    for s = 1 : rS
        for t = 1 :rS
            if s == t
                I(s,t) = 0;
            else
                I(s,t) = -inf;
            end
        end
    end
    
    I;
    PP = [S;T]; 
    QQ = [I;I]; 
    kut=1;

    [rpp,cpp] = size(PP);  % rp,cp are the no. of rows and columns of matrix P respectively
    [rqq,cqq] = size(QQ);  % rq,cq are the no. of rows and columns of matrix Q respectively

    w = zeros(1,cpp);

    wi = w;
    
    sPP = -PP';          % Negation and transposition of P
    sQQ = -QQ';          % Negation and transposition of Q

    rr = 0;             % A non-negative integer and act as the condition of the (while) loop
    step = 1;          % count the number of step and act as the condition of the (while)loop
    xn=[];            % Any x generated in the program depending on r, called x new
    yo=[];            % y generated at the first step with r=0, called y old
    yn=[];            % Any y generated in the program depending on r,called y new

    wo=w;              
    wi = wo;

    while rr ~= -1
        PPx = [];      % Maximum of ax1 for each i
        PPx1 = [];     % Sum of ith row of matrix P and vector x
        PPy = [];      % Maximum of Py1 for each i
        QQy1 = [];     % Sum of ith row of matrix Q and vector y
        kx = [];
        ky = [];
        w1 = [];      % Sum of ith row of matrix sa and ax 
        ww = [];
        wy = [];
        y1 = [];      % Sum of ith row of matrix sb and ax
        y = [];       % minimum of y1 for each i
        zz = [];
        zy = [];
            
        for i = 1 : rpp
            PPx1 = PP(i,:) + w;
            PPx(i) = max(PPx1);
        end
        PPx;
    
        for i = 1 : cqq
            y1 = sQQ(i,:) + PPx;
            y(i) = min(y1);
        end
        
        y;
        yn=y;
        
        for i = 1 : rqq
            QQy1 = QQ(i,:) + y;
            QQy(i) = max(QQy1);
        end
        
        if isequal(PPx,QQy)
            rr = -1;
            Existence =1;
            for j = 1 : cpp-1
                xhat1(kut)=w(j)-w(cpp);
                kut=kut+1;
            end
            xhat = xhat1;
            rr=-1;
        else
            for i = 1 : cpp
                w1 = sPP(i,:) + QQy;
                w(i) = min(w1);
            end
            w;
            wn=w;
            rr = rr + 1;
            if step == 1
                yo=y;
                yi = yo;
            end
            
            yo;
            
            for i = 1 : cpp
                if wn(i) < wo(i)
                    t = 1;
                else
                    t = 0;
                end
                ww = union(t,kx);
                zz = union(zz,ww);
            end
            
            if all(zz)
                Existence =0;
                rr=-1;
            end
            
            if rr ~=-1
                for i = 1 : cqq
                    if yn(i) < yo(i)
                        ty = 1;
                    else
                    ty = 0;
                    end
                    wy = union(ty,ky);
                    zy = union(zy,wy);
                end
                if all(zy)
                    Existence = 0;
                    rr=-1;
                end
                rr;  
                step = step+1;
            end
        end
    end

   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if Existence
        L=m;
    else
        U=m;        
    end
    m=(L+U)/2;
    if ~isequal(decimal,4)
        m = (10^(-decimal))*(round(m*10^(decimal)));
    end
    D = abs(U-L);
    D=(10^(-decimal))*(fix(D*10^(decimal)));       
end
disp('An optimal solution to hx --> max subject to Ax + c = Bx + d (all expressions are in max algebra) is :');
osol = xhat - h;
disp(osol);
disp('The optimal value is:')
fmax=L;
ovalue1 = osol + h;
ovalue = max(ovalue1);
disp(ovalue);
sol=0;
end



