%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   A program for finding a solution to the homogeneous system Ax = Bx in max algebra where A & B are matrices.       %
%   Stepping stone method is used in finding a solution to the system.                                                % 
%   The user have a choice of  input entry manually or by random.                                                     % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function homogeneous2 
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 = Bx 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 and columns for A & B must be the same.')
    disp('... Press any key to continue ...')
    pause
    B = input('Enter the matrix B: ');
    disp(B);
    [rowAB,colAB]=size(A);
    [rowAB,colAB]=size(B);
else
    rowAB = input('Enter the number of rows for matrices A and B : ');
    colAB = input('Enter the number of columns for matrices A and 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,colAB);
        for i = 1 : rowAB
            for j = 1 : colAB
                RA(i,j) = round(RA(i,j) * range) + Lbound;
            end
        end
        A = RA;
        
        RB = rand(rowAB,colAB);
        for i = 1 : rowAB
            for j = 1 : colAB
                RB(i,j) = round(RB(i,j) * range) + Lbound;
            end
        end
        B = RB;
                
        A
        B
        
    else
        
        range = 100;
        
        RA = rand(rowAB,colAB);
        RA = round(RA * range);
        A = RA;
        
        RB = rand(rowAB,colAB);
        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 zero vector.')

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

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

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

iterations = 1;  % count the number of iterations and act as the condition of the (while)loop

while iterations ~= 0
    [m,n] = size(A);  % m,n are the no. of rows and columns of matrix A respectively
    [m,n] = size(B);  % m,n are the no. of rows and columns of matrix B respectively
    N = 1 : n;        % Set of integers from 1 to n
    S = 1 : m;        % Set of integers from 1 to m 
    
    % Begining of the counters to help in getting members of some sets
    c = 1;
    cc = 1;
    e = 1;
    ee = 1;
    e3 = 1;
    ee3 = 1;    
    ja = 1; 
    jj = 1;
    jk = 1;
    jl = 1;
    jm = 1;
    jn = 1;
    jo = 1;
    k = 1;
    l = 1;
    v = 1;
    % End of the counters 
    
    Ay = [];          % Same as Py
    aa = [];          % A vector, which contains maximum of p for each row 
    al = [];          % element of p for each i in L1 and j in NP
    alpha = [];       % maximum element of p for each i in L1 and j in NP
    bb = [];          % A vector, which contains maximum of p for each row
    beta = [];        % maximum element of q for each i in L2 and j in NP
    beta3 = [];       % beta for finding t3
    bt = [];          % element of q for each i in L2 and j in NP
    E = [];           % set containing indices of rows s.t aa = bb
    E1 = [];          % set of indices i in E  s.t ff in Py and gg not in Py
    E2 = [];          % set of indices i in E  s.t ff not in Py and gg in Py
    f = [];           % set of indices(columns)which gives the maximun on the RHS
    ff = [];          % same as f
    ffp = [];         % set of indices j(columns)in ff but not in Py
    ffs=[];           % set formed to detect the elements of ff in (not in) Py 
    fs=[];            % set formed to detect the elements of f in (not in) Py 
    g = [];           % set of indices(columns)which gives the maximun on the LHS
    gg = [];          % same as g
    ggp = [];         % set of indices j(columns)in gg but not in Py
    ggs=[];           % set formed to detect the elements of gg in (not in) Py 
    gs=[];            % set formed to detect the elements of g in (not in) Py 
    H = [];           % set containing indices of rows s.t aa > bb
    kl = [];          % set formed to help in getting the elements of ugfp
    L1 = [];          % set of indices( no. of rows) s.t f is subset of Py
    L2 = [];          % set of indices( no. of rows) s.t g is subset of Py
    L3 = [];          % set of indices( no. of rows) in H s.t g is not a subset of Py
    LL11 = [];        % Formed to help in getting the elements of L1
    LL12 = [];        % Formed to help in getting the elements of L1
    LL21 = [];        % Formed to help in getting the elements of L2
    LL22 = [];        % Formed to help in getting the elements of L2
    NP = [];          % Elements in N not in Py
    Py = [];          % set containing indices j (of x(j))to be decreased
    p = [];           % A matrix, which contains the sum of row i in a and y
    q = [];           % A matrix, which contains the sum of row i in b and y
    repeat = 1;       % A condition of the (while)loop
    t = [];           % minimum of t1, t2 and t3
    t1 = [];          % minimum value of the difference btw aa and alpha for  each i in L1
    t2 = [];          % minimum value of the difference btw bb and beta for each i in L2
    t3 = [];          % minimum value of the difference btw aa and beta3 for each i in L3
    temp1 = [];       % Temporary matrix to help in swapping
    temp2 = [];       % Temporary matrix to help in swapping
    temp3 = [];       % Temporary matrix to help in swapping
    tt1 = [];         % Formed to help in getting the elements of t1
    tt2 = [];         % Formed to help in getting the elements of t2
    tt3 = [];         % Formed to help in getting the elements of t3
    tk = [];          % row vector of t1,t2, and t3
    ugfp = [];        % The overall union of ffp and ggp for each i
    ugfp1 = [];       % union of ffp and ggp for each i
    w = [];           % set formed to help in getting the elements of Ay
    
    
    
    for i = 1 : m
        y = x;
        p(i,:) = (A(i,:) + y);
        aa(i) = max( p(i,:) );
        q(i,:) = (B(i,:) + y);
        bb(i )= max(q(i,:));
        
        if aa(i) < bb(i)
            temp1 = q(i,:);
            q(i,:) = p(i,:);
            p(i,:) = temp1;
            temp2 = bb(i);
            bb(i) = aa(i);
            aa(i) = temp2;
            temp3 = B(i,:);
            B(i,:) = A(i,:);
            A(i,:) = temp3;
        elseif aa(i) == bb(i)
            p(i,:) = p(i,:);
            q(i,:) = q(i,:);
            A(i,:) = A(i,:);
            B(i,:) = B(i,:);
        end
    end
    
   % aa
    %bb
    %p
    %q
    %A
    %B
    
    for i = 1 : m
        if aa(i) > bb(i)
            H(ja) = i;
            ja = ja + 1;
        end
    end
    
    %H
    E = setdiff(S,H);
    
    if isempty(H)
        disp('The problem has a solution. The solution is: ');
        disp(y);
        return;
    end
    
    for i = 1 : m
        if ismember(i,H)
            f = find(aa(i) == p(i,:)); 
            Ay = union(f,w);
        end
        w = Ay;
    end
    
    Py = Ay;
    
    while repeat ~= 0
        E1 = [];
        E2 = [];
        k = 1;
        l = 1;
        for i = 1 : m
            if ismember(i,E)
                ff = find(aa(i) == p(i,:));
                gg = find(bb(i) == q(i,:));
                
                ffs=ismember(ff,Py);
                ggs=ismember(gg,Py);
                if all(ffs) & ~all(ggs)
                    E1(k) = i;
                    k = k + 1;
                    ggp = gg;
                end
                if ~all(ffs) & all(ggs)
                    E2(l) = i;
                    ffp = ff;
                    l = l + 1;
                end
                ugfp1 = union(ggp,ffp);
                ugfp = union(ugfp1,kl);
                kl = ugfp;
            end
        end
        
        if isempty(union(E1,E2))
            %disp('Py is the set of indices of variable to be decreased');
            %disp(Py);
            repeat = 0;
        else
            Py = union(Py,ugfp);
            if isequal(Py,N)
                disp('The problem has no solution.');
                return;
            end
            repeat = repeat + 1;
        end
    end
    
    %Py
    
    if isequal(Py,N)
        disp('The problem has no solution.');
        return;
    end
    
    for i = 1 : m
        if ismember(i,S)
            if ismember(i,H)
                f = find(aa(i) == p(i,:));
                fs = ismember(f,Py);
                
                if all(fs)
                    LL11(jj) = i;
                    jj = jj + 1;
                end
                g = find(bb(i) == q(i,:));
                gs = ismember(g,Py);
               
                if all(gs)
                    LL21(jl) = i;
                    jl = jl + 1;
                end
            end
            if ismember(i,E) 
                ff = find(aa(i) == p(i,:));
                ffs = ismember(ff,Py);
                
                if all(ffs)
                    LL12(jk) = i;
                    jk = jk + 1;
                end
                gg = find(bb(i) == q(i,:));
                ggs = ismember(gg,Py);
                
                if all(ggs)
                    LL22(jm) = i;
                    jm = jm + 1;
                end
            end
            L1 = union(LL11,LL12);
            L2 = union(LL21,LL22);
        end
    end
    
    for i = 1 : m
        if ismember(i,H)
            g = find(bb(i) == q(i,:));
            gs = ismember(g,Py);
            
            if ~all(gs)
                L3(jn) = i;
                jn = jn + 1;
            end
        end
    end
    
     NP = setdiff(N,Py);
    
    for i = 1 : m
        al = [];      % we want different al for each i
        for j = 1 : n
            if ismember(i,L1)
                if ismember(j,NP)
                    al(c) = A(i,j) + x(j);
                    c = c + 1;
                end
            end
        end
        alpha = max(al);
        c = 1;        % for each i we want a different al
        if ismember(i,L1)
            tt1(e) = aa(i) - alpha;
            e = e + 1;
        end
    end
    
    if isempty(tt1)
        %disp('t1 is empty');
        t1 = inf;
    else
        t1 = min(tt1);
    end
    
    for i = 1 : m
        bt = [];       % we want different bt for each i
        for j = 1 : n
            if ismember(i,L2)
                if ismember(j,NP)
                    bt(cc) = B(i,j) + x(j);
                    cc = cc + 1;
                end
            end
        end
        beta = max(bt);
        cc = 1;         % for each i we want a different bt
        if ismember(i,L2)
            tt2(ee) = bb(i) - beta;
            ee = ee + 1;
        end
    end
    
    if isempty(tt2)
        %disp('t2 is empty');
        t2 = inf;
    else
        t2 = min(tt2);
    end
    
    for i = 1 : m
        bt3 = [];       % we want different bt3 for each i
        for j = 1 : n
            if ismember(i,L3)
                if ismember(j,NP)
                    bt3(e3) = B(i,j) + x(j);
                    e3 = e3 + 1;
                end
            end
        end
        beta3 = max(bt3);
        e3 = 1;         % for each i we want a different bt3
        if ismember(i,L3)
            tt3(ee3) = aa(i)- beta3;
            ee3 = ee3 + 1;
        end
    end
    
    if isempty(tt3)
        %disp('t3 is empty');
        t3 = inf;
    else
        t3 = min(tt3);
    end
    
    tk = [t1, t2, t3];
    
    t=min(tk);
    
    for j = 1 : n
        if ismember(j,Py)
            y(j) = y(j) - t;
        else
            y(j) = y(j);
        end
    end
    
    y;
    
    if (~isequal(aa,bb))
        iterations = iterations + 1;
        x = y;
    end
end    
        
         
                 
       
        
           
        
        
        

    
    
    
    
