%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   A program for finding a solution to the system Ax = By in max algebra where A & B, matrices.       %
%   Alternating Method is used in finding a solution for the system.                                   % 
%   The user have a choice of  input entry manually or by random.                                      % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function homogeneous1 

disp('              THIS PROGRAM WAS DESIGNED AND WRITTEN BY     ');
disp('');
disp('                    ***** ABDULHADI AMINU *****    ');
disp('                    *****      710779     *****    ');
disp('                      AS PART OF HIS PHD WORK            ');
disp('                      UNDER THE SUPERVISION OF    ');
disp('                   ***** DR PETER BUTKOVIC *****'); 
disp('')
disp('                       SCHOOL OF MATHEMATICS           ');
disp('                   THE UNIVERSITY OF BIRMINGHAM        ');
disp('                             EDGBASTON                   ')
disp('                        BIRMINGHAM  B15 2TT            ' );
disp('                          UNITED KINGDOM   ');


disp('...Press any key to continue...')
pause

disp('This program is designed to find a solution to Ax = By in max algebra ')
disp('where A & B are given matrices or to indicate that no such solution exists.')
disp('You have the choice of entering matrices manually (row by row) or to generate them randomly.')
disp('Please use square brackets to indicate the begining or end of a matrix or a vector.')
disp('Use semicolons to separate the rows and spaces to separate the individual entries.')
disp('... Press any key to continue ...')
pause
entryopt = input('Press 1 to enter manually, 0 by random: ');

if entryopt
    A = input('Enter the matrix A: ');
    disp(A);
    disp('Note: The number of rows for A & B must be the same.')
    disp('... Press any key to continue ...')
    pause
    B = input('Enter the matrix B: ');
    disp(B);
    [rowAB,colA]=size(A);
    [rowAB,colB]=size(B);
    
else
    rowAB = input('Enter the number of rows for matrices A and B: ');
    colA = input('Enter the number of columns for matrices A: ');
    colB = input('Enter the number of columns for matrices B: ');
    disp('Note: Default range of entries for matrices is [0,100].');
    
    rangeopt = input('Press 1 to enter own range, 0 otherwise: ');
    
    if rangeopt
        Ubound = input('Enter an upper bound for the matrix entries: ');
        Lbound = input('Enter a lower bound for the matrix entries: ');
        range = Ubound - Lbound;
        
        RA = rand(rowAB,colA);
        for i = 1 : rowAB
            for j = 1 : colA
                RA(i,j) = round(RA(i,j) * range) + Lbound;
            end
        end
        A = RA;
        
        RB = rand(rowAB,colB);
        for i = 1 : rowAB
            for j = 1 : colB
                RB(i,j) = round(RB(i,j) * range) + Lbound;
            end
        end
        B = RB;
        
        A
        B
        
    else
        range = 100;
        
        RA = rand(rowAB,colA);
        RA = round(RA * range);
        A = RA;
        
        RB = rand(rowAB,colB);
        RB = round(RB * range);
        B = RB;
        
        A
        B
    end       
    
    disp('... Press any key to continue ...');
    pause
end

if isempty(A) | isempty(B)
    disp('The matrices can not be empty.')
    return
end

disp('Note: The starting vector is x=0.')

decision = input('Press 1 if you want a different starting vector, 0 otherwise: ');

if decision
    disp('Note: The number of columns in matrix A must be the same as the dimension of x.')
    disp('... Press any key to continue ...');
    pause
    x = input('Enter the starting vector x: ');
    disp(x);    
else
    x = zeros(1,colA)
end
xi = x;

disp('... Press any key to run the program ...')
pause

sA = -A';          % Negation and transposition of A
sB = -B' ;         % Negation and transposition of B

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

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
    
    Ax = [];      % Maximum of ax1 for each i
    Ax1 = [];     % Sum of ith row of matrix A and vector x
    By = [];      % Maximum of by1 for each i
    By1 = [];     % 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 : m
        Ax1 = A(i,:) + x;
        Ax(i) = max(Ax1);
    end
    
    Ax;
    
    for i = 1 : k
        y1 = sB(i,:) + Ax;
        y(i) = min(y1);
    end 
    
    y;
    
    yn=y;
    
    for i = 1 : m       
        By1 = B(i,:) + y;
        By(i) = max(By1);
    end
    
    By;
    
    if isequal(Ax,By)
        disp('The problem has a solution. The solution is: ');
        r;
        x
        y
        
        r = -1;
        
    else
        
        for i = 1 : n
            x1 = sA(i,:) + By;
            x(i) = min(x1);
        end
        
        x; 
        
        xn=x;
        
        r = r + 1;
        
        if step == 1
            
            yo=y;
        end
        
        yo;
        
        for i = 1 : n
            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 no solution.')
            return
        end
        
        for i = 1 : k
            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 no solution.')
            return
        end
        
        r;
        step = step+1;
        
    end
end



