%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   A program for solving max-linear programs with two-sided constraints (maximisation case).          %  
%   Alternating method is used in finding feasible and optimal solutions (if it exist).                % 
%   The Program is designed as 2 in 1.                                                                 %                                          
%   The current find a feasible solution (if it exist) and display 'No solution' otherwise.            % 
%   If a feasible solution exists maximisation2 will be engaged.                                       %              
%   The second (maximisation2) will then find an optimal solution to the problem.                      %                   
%   The user have a choice of  input entry manually or by random.                                      % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function maximisation
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 an optimal solution for the max-linear programs')
disp('(maximisation) with real\integer coefficients and non-homogeneous constraints,')
disp('Ax + c = Bx + d, where A & B are given matrices, c & d are given vectors.')
disp('The program checks the feasibility before the optimisation begins.')
disp('The optimisation part of the program works with objective function whose all ')
disp('coefficients are zero by default. There will be a possibility to change these once ')
disp('the feasibility is confirmed.')
disp('You have the choice of entering the matrices and vectors manually (row by row) or')
disp('to generate them randomly.')
disp('... Press any key to continue ...')
pause
disp('Please use square brackets to indicate the beginning 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
disp('Note: The maximum number of decimal places for all entries is 4.')
disp('      The default number of decimal places is 4.');
disp('      For problems with integer entries, please enter 0 as the number of decimal places.')
disp('... Press any key to continue ...')
pause
decimal = input('Enter a nonnegative number, the number of decimal places for all entries: ');
if isempty(decimal)
    disp('The number of decimal places cannot be empty.');
    return
end
if decimal > 4
    disp('The number of decimal places cannot be more than 4.');
    return
end
if decimal < 0
    disp('The number of decimal places cannot be negative.');
    return
end
if decimal > 0
    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,10].');
    
        rangeopt = input('Press 1 to enter own range, 0 otherwise: ');
        if rangeopt
            disp('Note: The range of entries for matrix will now be [0,range].')
            range = input('Enter a range for the matrix entries: ');
            RA = rand(rowAB,colAB);
            for i = 1 : rowAB
                for j = 1 : colAB
                    R1=rand();
                    if R1 < 0.5
                        RA(i,j) = RA(i,j);
                    else
                        RA(i,j) = RA(i,j) * range;
                    end
                end
            end
            A = RA;
            RB = rand(rowAB,colAB);
            for i = 1 : rowAB
                for j = 1 : colAB
                    R2=rand();
                    if R2 < 0.5
                        RB(i,j) = RB(i,j);
                    else
                        RB(i,j) = RB(i,j) * range;
                    end
                end
            end
            B = RB;
            A
            B
            Rc = rand(rowAB,1);
            for j = 1 : rowAB
                c1 = rand();
                if c1 < 0.5
                    Rc(j) = Rc(j);
                else
                    Rc(j) = Rc(j) * range;
                end
            end
            c = Rc
            Rd = rand(rowAB,1);
            for j = 1 : rowAB
                d1 = rand();
                if d1 < 0.5
                    Rd(j) = Rd(j);
                else
                    Rd(j) = Rd(j) * range;
                end
            end
            d = Rd
        else
            range = 10;
            RA = rand(rowAB,colAB);
            for i = 1 : rowAB
                for j = 1 : colAB
                    R1=rand();
                    if R1 < 0.5
                        RA(i,j) = RA(i,j);
                    else
                        RA(i,j) = RA(i,j) * range;
                    end
                end
            end
            A = RA;
            RB = rand(rowAB,colAB);
            for i = 1 : rowAB
                for j = 1 : colAB
                    R2=rand();
                    if R2 < 0.5
                        RB(i,j) = RB(i,j);
                    else
                        RB(i,j) = RB(i,j) * range;
                    end
                end
            end
            B = RB;
            Rc = rand(rowAB,1);
            for j = 1 : rowAB
                c1 = rand();
                if c1 < 0.5
                    Rc(j) = Rc(j);
                else
                    Rc(j) = Rc(j) * range;
                end 
            end
            c = Rc;
            Rd = rand(rowAB,1);
            for j = 1 : rowAB
                d1 = rand();
                if d1 < 0.5
                    Rd(j) = Rd(j);
                else
                    Rd(j) = Rd(j) * range;
                end 
            end
            d = Rd;
        end
        disp('... Press any key to continue ...');
        pause
    end
    power = 10^decimal;
    mpower= 10^(-decimal);
    epsilon = 10^(-decimal)
    if ~isequal(decimal,4)
        A =(mpower)*(fix(A*power));
        B =(mpower)*(fix(B*power));
        c = (mpower)*(fix(c*power));
        d = (mpower)*(fix(d*power));
    end
else %(for the if decimal >0)
    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: ')
        power = 10^decimal;
        mpower= 10^(-decimal);
        A =(mpower)*(fix(A*power));
        B =(mpower)*(fix(B*power));
        c = (mpower)*(fix(c*power));
        d = (mpower)*(fix(d*power));
    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;
            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;
            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
    epsilon =1
end

    %end of integer
A
B
c
d
%epsilon=0.3
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('...Press any key to continue');
%pause;
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('A feasible solution is (x0):');
        r=-1;
        for j = 1 : n-1
            xfinal(j) = x(j)-x(n);
        end
        %disp(xfinal);
        xzero = xfinal;
        disp(xzero); 
        disp('The objective function value is :');
        disp(max(xfinal));
        cequald = 0;
        
        [sol] = maximisation2(A,B,c,d,xfinal,epsilon,decimal);
        
        if cequald
            return;
        end
        
        if sol == 0
           return  ;
        end
    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('No feasible solution exists for the problem.');
            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('No feasible solution exists for the problem.');
            return
        end
        r;   
        step = step+1;  
       
    end
end






