function f = heuristic_jrp



clear all;

choice = 3 ;

while choice ~= 0
    choice = menu('"THE JOB ROTATION PROBLEM" (MAXIMIZATION)'     , ...
             '0. Stop'        , ...
             '1. Enter the matrix by hand'  , ...
	         '2. Random generation of the matrix  ') ;


    switch(choice)
        
        case 1
            choice = 0 ;
               
        case 2
            choice = 1 ;
            A = input('\n\n\n - Enter your matrix in the Matlab format - \n ex : [a b\n       c d]     \n\nA = ');
            [n,m] = size(A);
            k = input('\n\n\n - Enter your choice for k: ');
            [li,co] = size(A) ;

            if ( (n < k) | (n < 1) | (k < 1) )
                error('Error on n or k')
            elseif ( li ~= co)
                error('This is not a square matrix') 
            else       
                sprintf('\n\n\n  We consider the following matrix : ')
                A
                [f1,f2,f3,f4] = heuristic4(A,k) ;
                sprintf('\n\n  Indices of the principal submatrix : %s ',mat2str(f1)) 
                sprintf('\n\n  The principal submatrix : ')
                f2
                sprintf ('\n\n  Optimal solution to the assignment problem for the principal submatrix : %s ',mat2str(f3)) 
                sprintf ('\n - Optimal or suboptimal value for the job rotation problem : %d' , f4)

            end
            

  

        case 3
            choice = 2 ;
            choice2 = menu('Choose your command'     , ...
             '1. Random entries of the matrix : integers from [0,100]'  , ...
	         '2. Choice of the lower and the upper bound  ') ;
            
           switch(choice2)
               
                case 1
                    n = input('\n\n\n - Enter n : ');
                    k = input('\n\n\n - Enter k : ');
                    A = floor(100*rand(n)); 
                    [li,co] = size(A) ;

                    if ( (n < k) | (n < 1) | (k < 1) )
                        error('Error on n or k')
                    elseif ( li ~= co)
                        error('This is not a square matrix') 
                    else       
                        sprintf('\n\n\n  We consider the following matrix : ')
                        A
                        [f1,f2,f3,f4] = heuristic4(A,k) ;
                        sprintf('\n\n  Indices of the principal submatrix : %s ',mat2str(f1)) 
                        sprintf('\n\n  The principal submatrix : ')
                        f2
                        sprintf ('\n\n  Optimal solution to the assignment problem for the principal submatrix : %s ',mat2str(f3)) 
                        sprintf ('\n - Optimal or suboptimal value for the job rotation problem : %d' , f4)

                    end
                    
                case 2
                    n = input('\n\n\n - Enter n : ');
                    k = input('\n\n\n - Enter k : ');
                    l = input('\n\n\n - Enter your choice for lower bound:  l = ');
                    u = input('\n\n\n - Enter your choice for upper bound:  u = ');
                    sprintf('\n\n\n - We consider the following matrix - ') 
                    A = floor(l*ones(n) + (u-l)*rand(n)); 
                    [li,co] = size(A) ;

                    if ( (n < k) | (n < 1) | (k < 1) )
                        error('Error on n or k')
                    elseif ( li ~= co)
                        error('This is not a square matrix') 
                    else       
                        sprintf('\n\n\n  We consider the following matrix : ')
                        A
                        [f1,f2,f3,f4] = heuristic4(A,k) ;
                        sprintf('\n\n  Indices of the principal submatrix : %s ',mat2str(f1)) 
                        sprintf('\n\n  The principal submatrix is : ')
                        f2
                        sprintf ('\n\n  Optimal solution to the assignment problem for the principal submatrix : %s ',mat2str(f3)) 
                        sprintf ('\n - Optimal or suboptimal value for the job rotation problem : %d' , f4)

                    end
                end
            
            
            
            
    end;
    
end;
 

%-----------------------------------------------------------------------------------------------------------


%compute the heuristic 2
function [f1,f2,f3] = heuristic2(A,V,k)
% A matrix n*n
% V vector of size k

[m,n] = size(A) ;
       
S1 = submatrix(A,V,k) ;
[iz,T] = hungarian(-S1') ;
perm = iz ;
T = -T ;
c2=T ;
W = V ;
cpt =  0;

compW = complementary_general(W,n,k) ;
 
while (cpt < k) 
    indice = W(k-cpt) ;
    indice_comp = 0 ;
    X=W ;
    for i = 1:(n-k)
        X(k-cpt) = compW(i) ; 
        S = submatrix(A,X,k) ;
        [iz, c] = hungarian(-S') ; 
        c = -c ;
        if (c > c2)
            c2 = c ;
            indice = compW(i) ;
            indice_comp = i ;
            perm = iz ;
        end
    end
    if (W(k-cpt)~=indice)
        compW(indice_comp)=W(k-cpt);
        W(k-cpt) = indice ;
    end
    cpt = cpt + 1 ;
end
f1 = W ;
f2 = perm ;
f3 = c2 ;


%------------------------------------------------------------------------------------------


% compute heuristic 4
function [f1,f2,f3,f4] = heuristic4(A,k)

[n,m] = size(A) ;

[V1,V2] = heuristic4_temp(A,k) ;

[r1,r2,r3] = heuristic2(A,V1,k) ;
[r4,r5,r6] = heuristic2(A,V2,k) ;


if (r3 > r6)
    V = order(r1) ;
    S = submatrix(A,V,k) ;
    [iz, c] = hungarian(-S') ; 
    c = -c ;

else
    V = order(r4) ;
    S = submatrix(A,V,k) ;
    [iz, c] = hungarian(-S') ; 
    c = -c ;

end

if (c == max([r3 r6]))
    sprintf('hello--------------------------------------------')
end


f1 = V ;
f2 = S ;
f3 = iz ;
f4 = c ;

%----------------------------------------------------------------------------------------------

% help function for heuristic 4
function [f1,f2] = heuristic4_temp(A,k)

% We transform A of size n in B of size 2*n-k such that
% B = [   A      0...0
%                  .
%                  .
%      0...0  -inf...-inf   ]

[n,m] = size(A) ;

B = zeros(2*n-k,2*n-k) ;
B(1:n,1:n) = A ;
B(n+1:2*n-k,n+1:2*n-k) = -1e+8 ;


[C,T] = hungarian(-B') ;


cpt = 1 ;

for i = 1:n
    if C(i) < (n+1)
        lig(cpt) = i ;
        col(cpt) = C(i) ;
        cpt = cpt + 1 ;
    end
end

f1=lig ;
f2=col ;



%--------------------------------------------------------------------------------------------


function [C,T]=hungarian(A)
%HUNGARIAN Solve the Assignment problem using the Hungarian method.
%
%[C,T]=hungarian(A)
%A - a square cost matrix.
%C - the optimal assignment.
%T - the cost of the optimal assignment.

% Adapted from the FORTRAN IV code in Carpaneto and Toth, "Algorithm 548:
% Solution of the assignment problem [H]", ACM Transactions on
% Mathematical Software, 6(1):104-111, 1980.

% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.
%                 Department of Computing Science, Umeå University,
%                 Sweden. 
%                 All standard disclaimers apply.

% A substantial effort was put into this code. If you use it for a
% publication or otherwise, please include an acknowledgement or at least
% notify me by email. /Niclas

[m,n]=size(A);

if (m~=n)
    error('HUNGARIAN: The cost matrix must be square');
end

% Save original cost matrix.
orig=A;

% Reduce matrix.
A=hminired(A);

% Do an initial assignment.
[A,C,U]=hminiass(A);

% Repeat while we have unassigned rows.
while (U(n+1))
    % Start with no path, no unchecked zeros, and no unexplored rows.
    LR=zeros(1,n);
    LC=zeros(1,n);
    CH=zeros(1,n);
    RH=[zeros(1,n) -1];
    
    % No labelled columns.
    SLC=[];
    
    % Start path in first unassigned row.
    r=U(n+1);
    % Mark row with end-of-path label.
    LR(r)=-1;
    % Insert row first in labelled row set.
    SLR=r;
    
    % Repeat until we manage to find an assignable zero.
    while (1)
        % If there are free zeros in row r
        if (A(r,n+1)~=0)
            % ...get column of first free zero.
            l=-A(r,n+1);
            
            % If there are more free zeros in row r and row r in not
            % yet marked as unexplored..
            if (A(r,l)~=0 & RH(r)==0)
                % Insert row r first in unexplored list.
                RH(r)=RH(n+1);
                RH(n+1)=r;
                
                % Mark in which column the next unexplored zero in this row
                % is.
                CH(r)=-A(r,l);
            end
        else
            % If all rows are explored..
            if (RH(n+1)<=0)
                % Reduce matrix.
                [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR);
            end
            
            % Re-start with first unexplored row.
            r=RH(n+1);
            % Get column of next free zero in row r.
            l=CH(r);
            % Advance "column of next free zero".
            CH(r)=-A(r,l);
            % If this zero is last in the list..
            if (A(r,l)==0)
                % ...remove row r from unexplored list.
                RH(n+1)=RH(r);
                RH(r)=0;
            end
        end
        
        % While the column l is labelled, i.e. in path.
        while (LC(l)~=0)
            % If row r is explored..
            if (RH(r)==0)
                % If all rows are explored..
                if (RH(n+1)<=0)
                    % Reduce cost matrix.
                    [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR);
                end
                
                % Re-start with first unexplored row.
                r=RH(n+1);
            end
            
            % Get column of next free zero in row r.
            l=CH(r);
            
            % Advance "column of next free zero".
            CH(r)=-A(r,l);
            
            % If this zero is last in list..
            if(A(r,l)==0)
                % ...remove row r from unexplored list.
                RH(n+1)=RH(r);
                RH(r)=0;
            end
        end
        
        % If the column found is unassigned..
        if (C(l)==0)
            % Flip all zeros along the path in LR,LC.
            [A,C,U]=hmflip(A,C,LC,LR,U,l,r);
            % ...and exit to continue with next unassigned row.
            break;
        else
            % ...else add zero to path.
            
            % Label column l with row r.
            LC(l)=r;
            
            % Add l to the set of labelled columns.
            SLC=[SLC l];
            
            % Continue with the row assigned to column l.
            r=C(l);
            
            % Label row r with column l.
            LR(r)=l;
            
            % Add r to the set of labelled rows.
            SLR=[SLR r];
        end
    end
end

% Calculate the total cost.
T=sum(orig(logical(sparse(C,1:size(orig,2),1))));


function A=hminired(A)
%HMINIRED Initial reduction of cost matrix for the Hungarian method.
%
%B=assredin(A)
%A - the unreduced cost matris.
%B - the reduced cost matrix with linked zeros in each row.

% v1.0  96-06-13. Niclas Borlin, niclas@cs.umu.se.

[m,n]=size(A);

% Subtract column-minimum values from each column.
colMin=min(A);
A=A-colMin(ones(n,1),:);

% Subtract row-minimum values from each row.
rowMin=min(A')';
A=A-rowMin(:,ones(1,n));

% Get positions of all zeros.
[i,j]=find(A==0);

% Extend A to give room for row zero list header column.
A(1,n+1)=0;
for k=1:n
    % Get all column in this row. 
    cols=j(k==i)';
    % Insert pointers in matrix.
    A(k,[n+1 cols])=[-cols 0];
end


function [A,C,U]=hminiass(A)
%HMINIASS Initial assignment of the Hungarian method.
%
%[B,C,U]=hminiass(A)
%A - the reduced cost matrix.
%B - the reduced cost matrix, with assigned zeros removed from lists.
%C - a vector. C(J)=I means row I is assigned to column J,
%              i.e. there is an assigned zero in position I,J.
%U - a vector with a linked list of unassigned rows.

% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.

[n,np1]=size(A);

% Initalize return vectors.
C=zeros(1,n);
U=zeros(1,n+1);

% Initialize last/next zero "pointers".
LZ=zeros(1,n);
NZ=zeros(1,n);

for i=1:n
    % Set j to first unassigned zero in row i.
	lj=n+1;
	j=-A(i,lj);

    % Repeat until we have no more zeros (j==0) or we find a zero
	% in an unassigned column (c(j)==0).
    
	while (C(j)~=0)
		% Advance lj and j in zero list.
		lj=j;
		j=-A(i,lj);
	
		% Stop if we hit end of list.
		if (j==0)
			break;
		end
	end

	if (j~=0)
		% We found a zero in an unassigned column.
		
		% Assign row i to column j.
		C(j)=i;
		
		% Remove A(i,j) from unassigned zero list.
		A(i,lj)=A(i,j);

		% Update next/last unassigned zero pointers.
		NZ(i)=-A(i,j);
		LZ(i)=lj;

		% Indicate A(i,j) is an assigned zero.
		A(i,j)=0;
	else
		% We found no zero in an unassigned column.

		% Check all zeros in this row.

		lj=n+1;
		j=-A(i,lj);
		
		% Check all zeros in this row for a suitable zero in another row.
		while (j~=0)
			% Check the in the row assigned to this column.
			r=C(j);
			
			% Pick up last/next pointers.
			lm=LZ(r);
			m=NZ(r);
			
			% Check all unchecked zeros in free list of this row.
			while (m~=0)
				% Stop if we find an unassigned column.
				if (C(m)==0)
					break;
				end
				
				% Advance one step in list.
				lm=m;
				m=-A(r,lm);
			end
			
			if (m==0)
				% We failed on row r. Continue with next zero on row i.
				lj=j;
				j=-A(i,lj);
			else
				% We found a zero in an unassigned column.
			
				% Replace zero at (r,m) in unassigned list with zero at (r,j)
				A(r,lm)=-j;
				A(r,j)=A(r,m);
			
				% Update last/next pointers in row r.
				NZ(r)=-A(r,m);
				LZ(r)=j;
			
				% Mark A(r,m) as an assigned zero in the matrix . . .
				A(r,m)=0;
			
				% ...and in the assignment vector.
				C(m)=r;
			
				% Remove A(i,j) from unassigned list.
				A(i,lj)=A(i,j);
			
				% Update last/next pointers in row r.
				NZ(i)=-A(i,j);
				LZ(i)=lj;
			
				% Mark A(r,m) as an assigned zero in the matrix . . .
				A(i,j)=0;
			
				% ...and in the assignment vector.
				C(j)=i;
				
				% Stop search.
				break;
			end
		end
	end
end

% Create vector with list of unassigned rows.

% Mark all rows have assignment.
r=zeros(1,n);
rows=C(C~=0);
r(rows)=rows;
empty=find(r==0);

% Create vector with linked list of unassigned rows.
U=zeros(1,n+1);
U([n+1 empty])=[empty 0];


function [A,C,U]=hmflip(A,C,LC,LR,U,l,r)
%HMFLIP Flip assignment state of all zeros along a path.
%
%[A,C,U]=hmflip(A,C,LC,LR,U,l,r)
%Input:
%A   - the cost matrix.
%C   - the assignment vector.
%LC  - the column label vector.
%LR  - the row label vector.
%U   - the 
%r,l - position of last zero in path.
%Output:
%A   - updated cost matrix.
%C   - updated assignment vector.
%U   - updated unassigned row list vector.

% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.

n=size(A,1);

while (1)
    % Move assignment in column l to row r.
    C(l)=r;
    
    % Find zero to be removed from zero list..
    
    % Find zero before this.
    m=find(A(r,:)==-l);
    
    % Link past this zero.
    A(r,m)=A(r,l);
    
    A(r,l)=0;
    
    % If this was the first zero of the path..
    if (LR(r)<0)
        ...remove row from unassigned row list and return.
        U(n+1)=U(r);
        U(r)=0;
        return;
    else
        
        % Move back in this row along the path and get column of next zero.
        l=LR(r);
        
        % Insert zero at (r,l) first in zero list.
        A(r,l)=A(r,n+1);
        A(r,n+1)=-l;
        
        % Continue back along the column to get row of next zero in path.
        r=LC(l);
    end
end


function [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR)
%HMREDUCE Reduce parts of cost matrix in the Hungerian method.
%
%[A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR)
%Input:
%A   - Cost matrix.
%CH  - vector of column of 'next zeros' in each row.
%RH  - vector with list of unexplored rows.
%LC  - column labels.
%RC  - row labels.
%SLC - set of column labels.
%SLR - set of row labels.
%
%Output:
%A   - Reduced cost matrix.
%CH  - Updated vector of 'next zeros' in each row.
%RH  - Updated vector of unexplored rows.

% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.

n=size(A,1);

% Find which rows are covered, i.e. unlabelled.
coveredRows=LR==0;

% Find which columns are covered, i.e. labelled.
coveredCols=LC~=0;

r=find(~coveredRows);
c=find(~coveredCols);

% Get minimum of uncovered elements.
m=min(min(A(r,c)));

% Subtract minimum from all uncovered elements.
A(r,c)=A(r,c)-m;

% Check all uncovered columns..
for j=c
    % ...and uncovered rows in path order..
    for i=SLR
        % If this is a (new) zero..
        if (A(i,j)==0)
            % If the row is not in unexplored list..
            if (RH(i)==0)
                % ...insert it first in unexplored list.
                RH(i)=RH(n+1);
                RH(n+1)=i;
                % Mark this zero as "next free" in this row.
                CH(i)=j;
            end
            % Find last unassigned zero on row I.
            row=A(i,:);
            colsInList=-row(row<0);
            if (length(colsInList)==0)
                % No zeros in the list.
                l=n+1;
            else
                l=colsInList(row(colsInList)==0);
            end
            % Append this zero to end of list.
            A(i,l)=-j;
        end
    end
end

% Add minimum to all doubly covered elements.
r=find(coveredRows);
c=find(coveredCols);

% Take care of the zeros we will remove.
[i,j]=find(A(r,c)<=0);

i=r(i);
j=c(j);

for k=1:length(i)
    % Find zero before this in this row.
    lj=find(A(i(k),:)==-j(k));
    % Link past it.
    A(i(k),lj)=A(i(k),j(k));
    % Mark it as assigned.
    A(i(k),j(k))=0;
end

A(r,c)=A(r,c)+m;



%----------------------------------------------------------------------------------------

function f = submatrix(A,V,k)

% give the submatrix of A with the coefficient of V

B =  zeros(k,k) ;

for i = 1:k
    for j = 1:k
        B(i,j) = A(V(i),V(j)) ;
    end
end
 
f = B ;



%--------------------------------------------------------------------------------------


function f = complementary_general(V,n,k)

% we compute the complementary of V: compV
compV = ones(n-k,1) ;

j=1;
for i=1:n
    if (belong(i,V)==0)
        compV(j)=i ;
        j=j+1;
    end
end
f=compV ;


%------------------------------------------------------------------------------------

%test if the element e belongs to the vector V 
function f = belong(e,V)

n = length(V) ;

answer = 0 ;
cpt = 0 ;

while ( (cpt < n) & (answer == 0) )
    cpt = cpt + 1 ;
    if (V(cpt) == e)
        answer = 1 ;
    end
end

f = answer ;
    

%-----------------------------------------------------------------------------------



function f = order(V)


if (length(V)==1)
    W = V ;
else
    W = [V(1)] ;
    for (i = 2:length(V))
        [T,ansver,B,E] = ordertemp(V(i),W,[],[]) ;
        while (ansver == 0)
            [T1,ansver,B1,E1] = ordertemp(V(i),T,B,E) ;
            T = T1 ;
            B = B1 ;
            E = E1 ;
        end
        W = [B T E] ;
    end
end

f = W ;


%------------------------------------------------------------------------------


function [f1,f2,f3,f4] = ordertemp(e,V,beg,ends)

n = length(V) ;

if (n==0)
     ansver = 1 ;
     W = [e] ;
     
elseif (n==1)
    ansver = 1 ;
    if e > V(1)
        W = [V e] ;
    else
        W = [e V] ;
    end

else

    if e < V(1)
        W = [e V] ;
        ansver = 1 ;
    
    elseif e > V(n)
        W = [V e] ;
        ansver = 1 ;
    else
        W = V(2:n-1) ;
        ansver = 0 ;
        beg = [beg V(1)] ;
        ends = [V(n) ends] ;
    end
end

f1 = W ;
f2 = ansver ;
f3 = beg ;
f4 = ends ;




