wolffd@0: function [C,T]=hungarian(A)
wolffd@0: %HUNGARIAN Solve the Assignment problem using the Hungarian method.
wolffd@0: %
wolffd@0: %[C,T]=hungarian(A)
wolffd@0: %A - a square cost matrix.
wolffd@0: %C - the optimal assignment.
wolffd@0: %T - the cost of the optimal assignment.
wolffd@0: 
wolffd@0: % Adapted from the FORTRAN IV code in Carpaneto and Toth, "Algorithm 548:
wolffd@0: % Solution of the assignment problem [H]", ACM Transactions on
wolffd@0: % Mathematical Software, 6(1):104-111, 1980.
wolffd@0: 
wolffd@0: % v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.
wolffd@0: %                 Department of Computing Science, Umeå University,
wolffd@0: %                 Sweden. 
wolffd@0: %                 All standard disclaimers apply.
wolffd@0: 
wolffd@0: % A substantial effort was put into this code. If you use it for a
wolffd@0: % publication or otherwise, please include an acknowledgement or at least
wolffd@0: % notify me by email. /Niclas
wolffd@0: 
wolffd@0: [m,n]=size(A);
wolffd@0: 
wolffd@0: if (m~=n)
wolffd@0:     error('HUNGARIAN: Cost matrix must be square!');
wolffd@0: end
wolffd@0: 
wolffd@0: % Save original cost matrix.
wolffd@0: orig=A;
wolffd@0: 
wolffd@0: % Reduce matrix.
wolffd@0: A=hminired(A);
wolffd@0: 
wolffd@0: % Do an initial assignment.
wolffd@0: [A,C,U]=hminiass(A);
wolffd@0: 
wolffd@0: % Repeat while we have unassigned rows.
wolffd@0: while (U(n+1))
wolffd@0:     % Start with no path, no unchecked zeros, and no unexplored rows.
wolffd@0:     LR=zeros(1,n);
wolffd@0:     LC=zeros(1,n);
wolffd@0:     CH=zeros(1,n);
wolffd@0:     RH=[zeros(1,n) -1];
wolffd@0:     
wolffd@0:     % No labelled columns.
wolffd@0:     SLC=[];
wolffd@0:     
wolffd@0:     % Start path in first unassigned row.
wolffd@0:     r=U(n+1);
wolffd@0:     % Mark row with end-of-path label.
wolffd@0:     LR(r)=-1;
wolffd@0:     % Insert row first in labelled row set.
wolffd@0:     SLR=r;
wolffd@0:     
wolffd@0:     % Repeat until we manage to find an assignable zero.
wolffd@0:     while (1)
wolffd@0:         % If there are free zeros in row r
wolffd@0:         if (A(r,n+1)~=0)
wolffd@0:             % ...get column of first free zero.
wolffd@0:             l=-A(r,n+1);
wolffd@0:             
wolffd@0:             % If there are more free zeros in row r and row r in not
wolffd@0:             % yet marked as unexplored..
wolffd@0:             if (A(r,l)~=0 & RH(r)==0)
wolffd@0:                 % Insert row r first in unexplored list.
wolffd@0:                 RH(r)=RH(n+1);
wolffd@0:                 RH(n+1)=r;
wolffd@0:                 
wolffd@0:                 % Mark in which column the next unexplored zero in this row
wolffd@0:                 % is.
wolffd@0:                 CH(r)=-A(r,l);
wolffd@0:             end
wolffd@0:         else
wolffd@0:             % If all rows are explored..
wolffd@0:             if (RH(n+1)<=0)
wolffd@0:                 % Reduce matrix.
wolffd@0:                 [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR);
wolffd@0:             end
wolffd@0:             
wolffd@0:             % Re-start with first unexplored row.
wolffd@0:             r=RH(n+1);
wolffd@0:             % Get column of next free zero in row r.
wolffd@0:             l=CH(r);
wolffd@0:             % Advance "column of next free zero".
wolffd@0:             CH(r)=-A(r,l);
wolffd@0:             % If this zero is last in the list..
wolffd@0:             if (A(r,l)==0)
wolffd@0:                 % ...remove row r from unexplored list.
wolffd@0:                 RH(n+1)=RH(r);
wolffd@0:                 RH(r)=0;
wolffd@0:             end
wolffd@0:         end
wolffd@0:         
wolffd@0:         % While the column l is labelled, i.e. in path.
wolffd@0:         while (LC(l)~=0)
wolffd@0:             % If row r is explored..
wolffd@0:             if (RH(r)==0)
wolffd@0:                 % If all rows are explored..
wolffd@0:                 if (RH(n+1)<=0)
wolffd@0:                     % Reduce cost matrix.
wolffd@0:                     [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR);
wolffd@0:                 end
wolffd@0:                 
wolffd@0:                 % Re-start with first unexplored row.
wolffd@0:                 r=RH(n+1);
wolffd@0:             end
wolffd@0:             
wolffd@0:             % Get column of next free zero in row r.
wolffd@0:             l=CH(r);
wolffd@0:             
wolffd@0:             % Advance "column of next free zero".
wolffd@0:             CH(r)=-A(r,l);
wolffd@0:             
wolffd@0:             % If this zero is last in list..
wolffd@0:             if(A(r,l)==0)
wolffd@0:                 % ...remove row r from unexplored list.
wolffd@0:                 RH(n+1)=RH(r);
wolffd@0:                 RH(r)=0;
wolffd@0:             end
wolffd@0:         end
wolffd@0:         
wolffd@0:         % If the column found is unassigned..
wolffd@0:         if (C(l)==0)
wolffd@0:             % Flip all zeros along the path in LR,LC.
wolffd@0:             [A,C,U]=hmflip(A,C,LC,LR,U,l,r);
wolffd@0:             % ...and exit to continue with next unassigned row.
wolffd@0:             break;
wolffd@0:         else
wolffd@0:             % ...else add zero to path.
wolffd@0:             
wolffd@0:             % Label column l with row r.
wolffd@0:             LC(l)=r;
wolffd@0:             
wolffd@0:             % Add l to the set of labelled columns.
wolffd@0:             SLC=[SLC l];
wolffd@0:             
wolffd@0:             % Continue with the row assigned to column l.
wolffd@0:             r=C(l);
wolffd@0:             
wolffd@0:             % Label row r with column l.
wolffd@0:             LR(r)=l;
wolffd@0:             
wolffd@0:             % Add r to the set of labelled rows.
wolffd@0:             SLR=[SLR r];
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: % Calculate the total cost.
wolffd@0: T=sum(orig(logical(sparse(C,1:size(orig,2),1))));
wolffd@0: 
wolffd@0: 
wolffd@0: function A=hminired(A)
wolffd@0: %HMINIRED Initial reduction of cost matrix for the Hungarian method.
wolffd@0: %
wolffd@0: %B=assredin(A)
wolffd@0: %A - the unreduced cost matris.
wolffd@0: %B - the reduced cost matrix with linked zeros in each row.
wolffd@0: 
wolffd@0: % v1.0  96-06-13. Niclas Borlin, niclas@cs.umu.se.
wolffd@0: 
wolffd@0: [m,n]=size(A);
wolffd@0: 
wolffd@0: % Subtract column-minimum values from each column.
wolffd@0: colMin=min(A);
wolffd@0: A=A-colMin(ones(n,1),:);
wolffd@0: 
wolffd@0: % Subtract row-minimum values from each row.
wolffd@0: rowMin=min(A')';
wolffd@0: A=A-rowMin(:,ones(1,n));
wolffd@0: 
wolffd@0: % Get positions of all zeros.
wolffd@0: [i,j]=find(A==0);
wolffd@0: 
wolffd@0: % Extend A to give room for row zero list header column.
wolffd@0: A(1,n+1)=0;
wolffd@0: for k=1:n
wolffd@0:     % Get all column in this row. 
wolffd@0:     cols=j(k==i)';
wolffd@0:     % Insert pointers in matrix.
wolffd@0:     A(k,[n+1 cols])=[-cols 0];
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function [A,C,U]=hminiass(A)
wolffd@0: %HMINIASS Initial assignment of the Hungarian method.
wolffd@0: %
wolffd@0: %[B,C,U]=hminiass(A)
wolffd@0: %A - the reduced cost matrix.
wolffd@0: %B - the reduced cost matrix, with assigned zeros removed from lists.
wolffd@0: %C - a vector. C(J)=I means row I is assigned to column J,
wolffd@0: %              i.e. there is an assigned zero in position I,J.
wolffd@0: %U - a vector with a linked list of unassigned rows.
wolffd@0: 
wolffd@0: % v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.
wolffd@0: 
wolffd@0: [n,np1]=size(A);
wolffd@0: 
wolffd@0: % Initalize return vectors.
wolffd@0: C=zeros(1,n);
wolffd@0: U=zeros(1,n+1);
wolffd@0: 
wolffd@0: % Initialize last/next zero "pointers".
wolffd@0: LZ=zeros(1,n);
wolffd@0: NZ=zeros(1,n);
wolffd@0: 
wolffd@0: for i=1:n
wolffd@0:     % Set j to first unassigned zero in row i.
wolffd@0: 	lj=n+1;
wolffd@0: 	j=-A(i,lj);
wolffd@0: 
wolffd@0:     % Repeat until we have no more zeros (j==0) or we find a zero
wolffd@0: 	% in an unassigned column (c(j)==0).
wolffd@0:     
wolffd@0: 	while (C(j)~=0)
wolffd@0: 		% Advance lj and j in zero list.
wolffd@0: 		lj=j;
wolffd@0: 		j=-A(i,lj);
wolffd@0: 	
wolffd@0: 		% Stop if we hit end of list.
wolffd@0: 		if (j==0)
wolffd@0: 			break;
wolffd@0: 		end
wolffd@0: 	end
wolffd@0: 
wolffd@0: 	if (j~=0)
wolffd@0: 		% We found a zero in an unassigned column.
wolffd@0: 		
wolffd@0: 		% Assign row i to column j.
wolffd@0: 		C(j)=i;
wolffd@0: 		
wolffd@0: 		% Remove A(i,j) from unassigned zero list.
wolffd@0: 		A(i,lj)=A(i,j);
wolffd@0: 
wolffd@0: 		% Update next/last unassigned zero pointers.
wolffd@0: 		NZ(i)=-A(i,j);
wolffd@0: 		LZ(i)=lj;
wolffd@0: 
wolffd@0: 		% Indicate A(i,j) is an assigned zero.
wolffd@0: 		A(i,j)=0;
wolffd@0: 	else
wolffd@0: 		% We found no zero in an unassigned column.
wolffd@0: 
wolffd@0: 		% Check all zeros in this row.
wolffd@0: 
wolffd@0: 		lj=n+1;
wolffd@0: 		j=-A(i,lj);
wolffd@0: 		
wolffd@0: 		% Check all zeros in this row for a suitable zero in another row.
wolffd@0: 		while (j~=0)
wolffd@0: 			% Check the in the row assigned to this column.
wolffd@0: 			r=C(j);
wolffd@0: 			
wolffd@0: 			% Pick up last/next pointers.
wolffd@0: 			lm=LZ(r);
wolffd@0: 			m=NZ(r);
wolffd@0: 			
wolffd@0: 			% Check all unchecked zeros in free list of this row.
wolffd@0: 			while (m~=0)
wolffd@0: 				% Stop if we find an unassigned column.
wolffd@0: 				if (C(m)==0)
wolffd@0: 					break;
wolffd@0: 				end
wolffd@0: 				
wolffd@0: 				% Advance one step in list.
wolffd@0: 				lm=m;
wolffd@0: 				m=-A(r,lm);
wolffd@0: 			end
wolffd@0: 			
wolffd@0: 			if (m==0)
wolffd@0: 				% We failed on row r. Continue with next zero on row i.
wolffd@0: 				lj=j;
wolffd@0: 				j=-A(i,lj);
wolffd@0: 			else
wolffd@0: 				% We found a zero in an unassigned column.
wolffd@0: 			
wolffd@0: 				% Replace zero at (r,m) in unassigned list with zero at (r,j)
wolffd@0: 				A(r,lm)=-j;
wolffd@0: 				A(r,j)=A(r,m);
wolffd@0: 			
wolffd@0: 				% Update last/next pointers in row r.
wolffd@0: 				NZ(r)=-A(r,m);
wolffd@0: 				LZ(r)=j;
wolffd@0: 			
wolffd@0: 				% Mark A(r,m) as an assigned zero in the matrix . . .
wolffd@0: 				A(r,m)=0;
wolffd@0: 			
wolffd@0: 				% ...and in the assignment vector.
wolffd@0: 				C(m)=r;
wolffd@0: 			
wolffd@0: 				% Remove A(i,j) from unassigned list.
wolffd@0: 				A(i,lj)=A(i,j);
wolffd@0: 			
wolffd@0: 				% Update last/next pointers in row r.
wolffd@0: 				NZ(i)=-A(i,j);
wolffd@0: 				LZ(i)=lj;
wolffd@0: 			
wolffd@0: 				% Mark A(r,m) as an assigned zero in the matrix . . .
wolffd@0: 				A(i,j)=0;
wolffd@0: 			
wolffd@0: 				% ...and in the assignment vector.
wolffd@0: 				C(j)=i;
wolffd@0: 				
wolffd@0: 				% Stop search.
wolffd@0: 				break;
wolffd@0: 			end
wolffd@0: 		end
wolffd@0: 	end
wolffd@0: end
wolffd@0: 
wolffd@0: % Create vector with list of unassigned rows.
wolffd@0: 
wolffd@0: % Mark all rows have assignment.
wolffd@0: r=zeros(1,n);
wolffd@0: rows=C(C~=0);
wolffd@0: r(rows)=rows;
wolffd@0: empty=find(r==0);
wolffd@0: 
wolffd@0: % Create vector with linked list of unassigned rows.
wolffd@0: U=zeros(1,n+1);
wolffd@0: U([n+1 empty])=[empty 0];
wolffd@0: 
wolffd@0: 
wolffd@0: function [A,C,U]=hmflip(A,C,LC,LR,U,l,r)
wolffd@0: %HMFLIP Flip assignment state of all zeros along a path.
wolffd@0: %
wolffd@0: %[A,C,U]=hmflip(A,C,LC,LR,U,l,r)
wolffd@0: %Input:
wolffd@0: %A   - the cost matrix.
wolffd@0: %C   - the assignment vector.
wolffd@0: %LC  - the column label vector.
wolffd@0: %LR  - the row label vector.
wolffd@0: %U   - the 
wolffd@0: %r,l - position of last zero in path.
wolffd@0: %Output:
wolffd@0: %A   - updated cost matrix.
wolffd@0: %C   - updated assignment vector.
wolffd@0: %U   - updated unassigned row list vector.
wolffd@0: 
wolffd@0: % v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.
wolffd@0: 
wolffd@0: n=size(A,1);
wolffd@0: 
wolffd@0: while (1)
wolffd@0:     % Move assignment in column l to row r.
wolffd@0:     C(l)=r;
wolffd@0:     
wolffd@0:     % Find zero to be removed from zero list..
wolffd@0:     
wolffd@0:     % Find zero before this.
wolffd@0:     m=find(A(r,:)==-l);
wolffd@0:     
wolffd@0:     % Link past this zero.
wolffd@0:     A(r,m)=A(r,l);
wolffd@0:     
wolffd@0:     A(r,l)=0;
wolffd@0:     
wolffd@0:     % If this was the first zero of the path..
wolffd@0:     if (LR(r)<0)
wolffd@0:         ...remove row from unassigned row list and return.
wolffd@0:         U(n+1)=U(r);
wolffd@0:         U(r)=0;
wolffd@0:         return;
wolffd@0:     else
wolffd@0:         
wolffd@0:         % Move back in this row along the path and get column of next zero.
wolffd@0:         l=LR(r);
wolffd@0:         
wolffd@0:         % Insert zero at (r,l) first in zero list.
wolffd@0:         A(r,l)=A(r,n+1);
wolffd@0:         A(r,n+1)=-l;
wolffd@0:         
wolffd@0:         % Continue back along the column to get row of next zero in path.
wolffd@0:         r=LC(l);
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR)
wolffd@0: %HMREDUCE Reduce parts of cost matrix in the Hungerian method.
wolffd@0: %
wolffd@0: %[A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR)
wolffd@0: %Input:
wolffd@0: %A   - Cost matrix.
wolffd@0: %CH  - vector of column of 'next zeros' in each row.
wolffd@0: %RH  - vector with list of unexplored rows.
wolffd@0: %LC  - column labels.
wolffd@0: %RC  - row labels.
wolffd@0: %SLC - set of column labels.
wolffd@0: %SLR - set of row labels.
wolffd@0: %
wolffd@0: %Output:
wolffd@0: %A   - Reduced cost matrix.
wolffd@0: %CH  - Updated vector of 'next zeros' in each row.
wolffd@0: %RH  - Updated vector of unexplored rows.
wolffd@0: 
wolffd@0: % v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.
wolffd@0: 
wolffd@0: n=size(A,1);
wolffd@0: 
wolffd@0: % Find which rows are covered, i.e. unlabelled.
wolffd@0: coveredRows=LR==0;
wolffd@0: 
wolffd@0: % Find which columns are covered, i.e. labelled.
wolffd@0: coveredCols=LC~=0;
wolffd@0: 
wolffd@0: r=find(~coveredRows);
wolffd@0: c=find(~coveredCols);
wolffd@0: 
wolffd@0: % Get minimum of uncovered elements.
wolffd@0: m=min(min(A(r,c)));
wolffd@0: 
wolffd@0: % Subtract minimum from all uncovered elements.
wolffd@0: A(r,c)=A(r,c)-m;
wolffd@0: 
wolffd@0: % Check all uncovered columns..
wolffd@0: for j=c
wolffd@0:     % ...and uncovered rows in path order..
wolffd@0:     for i=SLR
wolffd@0:         % If this is a (new) zero..
wolffd@0:         if (A(i,j)==0)
wolffd@0:             % If the row is not in unexplored list..
wolffd@0:             if (RH(i)==0)
wolffd@0:                 % ...insert it first in unexplored list.
wolffd@0:                 RH(i)=RH(n+1);
wolffd@0:                 RH(n+1)=i;
wolffd@0:                 % Mark this zero as "next free" in this row.
wolffd@0:                 CH(i)=j;
wolffd@0:             end
wolffd@0:             % Find last unassigned zero on row I.
wolffd@0:             row=A(i,:);
wolffd@0:             colsInList=-row(row<0);
wolffd@0:             if (length(colsInList)==0)
wolffd@0:                 % No zeros in the list.
wolffd@0:                 l=n+1;
wolffd@0:             else
wolffd@0:                 l=colsInList(row(colsInList)==0);
wolffd@0:             end
wolffd@0:             % Append this zero to end of list.
wolffd@0:             A(i,l)=-j;
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: % Add minimum to all doubly covered elements.
wolffd@0: r=find(coveredRows);
wolffd@0: c=find(coveredCols);
wolffd@0: 
wolffd@0: % Take care of the zeros we will remove.
wolffd@0: [i,j]=find(A(r,c)<=0);
wolffd@0: 
wolffd@0: i=r(i);
wolffd@0: j=c(j);
wolffd@0: 
wolffd@0: for k=1:length(i)
wolffd@0:     % Find zero before this in this row.
wolffd@0:     lj=find(A(i(k),:)==-j(k));
wolffd@0:     % Link past it.
wolffd@0:     A(i(k),lj)=A(i(k),j(k));
wolffd@0:     % Mark it as assigned.
wolffd@0:     A(i(k),j(k))=0;
wolffd@0: end
wolffd@0: 
wolffd@0: A(r,c)=A(r,c)+m;