%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   A program for finding a solution to non homogeneous system Ax + c = Bx + d in max algebra where A & B matrices, c & d,vectors   %
%   Alternating method is used in finding the solution, if it exist or indicate that no such solution exist.                        % 
%   User have the choice of  input entry manually or by random                                                                      % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function nonhomogeneous; 
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 + c = Bx + d in max algebra')
disp('where A & B are given matrices, c & d are vectors 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);
    disp('Note: The vectors c & d must have the same number of rows as A & B')
    disp('...Press any key to continue...')
    pause
    c = input('Enter the vector c: ')
    d = input('Enter the vector d: ')
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
        
        Rc = rand(rowAB,1);
        for j = 1 : rowAB
            Rc(j) = round(Rc(j) * range) + Lbound;
        end
        c = Rc
        
        Rd = rand(rowAB,1);
        for j = 1 : rowAB
            Rd(j) = round(Rd(j) * range)+ Lbound;
        end
        d = Rd                
               
    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
        
        Rc = rand(rowAB,1);
        Rc = round(Rc * range);
        c = Rc
        
        Rd = rand(rowAB,1);
        Rd = round(Rd * range);
        d = Rd
    end    
    disp('... Press any key to continue ...');
    pause
    
end

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

A = [A c];
B = [B d];

[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

for i = 1 : m
    if c(i) < d(i)
        temp = A(i,:);
        A(i,:) = B(i,:);
        B(i,:)  = temp;
    end
end

%disp('After the arrangement the matrices are:')
%A
%B

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); 
    x = [x 0];
else
    x = zeros(1,colAB)
    x = [x 0];
end
xi = x;


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

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
xi = xo;

I = rand(m,m);

for s = 1 : m
    for t = 1 :m
        if s == t
            I(s,t) = 0;
        else
            I(s,t) = -inf;
        end
    end
end

I;
p = [A;B];
q = [I;I];

sp = -p';
sq = -q';

[mm,nn] = size(p);
[mm,kk] = size(q);
            
while r ~= -1
    
    px = [];      % Maximum of px1 for each i
    px1 = [];     % Sum of ith row of matrix p and vector x
    qy = [];      % Maximum of qy1 for each i
    qy1 = [];     % Sum of ith row of matrix q and vector y
    kx = [];
    ky = [];
    x1 = [];      % Sum of ith row of matrix sp and px 
    w = [];
    wy = [];
    y1 = [];      % Sum of ith row of matrix sq and px
    y = [];       % minimum of y1 for each i
    z = [];
    zy = [];
            
    for i = 1 : mm
        px1 = p(i,:) + x;
        px(i) = max(px1);
    end
    
    px;
    
    for i = 1 : kk
        y1 = sq(i,:) + px;
        y(i) = min(y1);
    end 
    
    y;
    
    yn=y;
    
    for i = 1 : mm       
        qy1 = q(i,:) + y;
        qy(i) = max(qy1);
    end
    
    qy;
    
    if isequal(px,qy)
        disp('The problem has solution. A solution is: ');
        r=-1;
        for j = 1 : n-1
            xfinal(j) = x(j)-x(n);
        end
        disp(xfinal);
        return;
        
    else
        
        for i = 1 : nn
            x1 = sp(i,:) + qy;
            x(i) = min(x1);
        end
        
        x; 
        
        xn = x;
        
        r = r + 1;
        
        if step == 1
            
            yo = y;
        end
        
        yo;
        
        for i = 1 : nn
            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 : kk
            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






