Euler 345th problem solution

euler345.m

function [Opt,Value]=euler345(InMatrix)
%Solves the Project Euler 345th Problem
%
%[Opt,Value]=euler(InMatrix)
%InMatrix - a square Opt matrix.
%Opt - the optimal assignment.
%Value - the Opt of the optimal assignment.

[m,n]=size(InMatrix);

if (m~=n)
    error('euler: Opt matrix must be square');
end

% Save original Opt matrix.
orig=InMatrix;

% Reduce matrix.
InMatrix=hminired(InMatrix);

% Do an initial assignment.
[InMatrix,Opt,U]=hminiass(InMatrix);

% 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 (InMatrix(r,n+1)~=0)
            % ...get column of first free zero.
            l=-InMatrix(r,n+1);

            % If there are more free zeros in row r and row r in not
            % yet marked as unexplored..
            if (InMatrix(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)=-InMatrix(r,l);
            end
        else
            % If all rows are explored..
            if (RH(n+1)<=0)
                % Reduce matrix.
                [InMatrix,CH,RH]=hmreduce(InMatrix,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)=-InMatrix(r,l);
            % If this zero is last in the list..
            if (InMatrix(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 Opt matrix.
                    [InMatrix,CH,RH]=hmreduce(InMatrix,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)=-InMatrix(r,l);

            % If this zero is last in list..
            if(InMatrix(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 (Opt(l)==0)
            % Flip all zeros along the path in LR,LC.
            [InMatrix,Opt,U]=hmflip(InMatrix,Opt,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=Opt(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 Opt.
Value=sum(orig(logical(sparse(Opt,1:size(orig,2),1))));


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

[~,n]=size(InMatrix);

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

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

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

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


function [InMatrix,Opt,U]=hminiass(InMatrix)
%HMINIASS Initial assignment of the euler method.
%
%[B,Opt,U]=hminiass(InMatrix)
%InMatrix - the reduced Opt matrix.
%B - the reduced Opt matrix, with assigned zeros removed from lists.
%Opt - a vector. Opt(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.

[n,np1]=size(InMatrix);

% Initalize return vectors.
Opt=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=-InMatrix(i,lj);

    % Repeat until we have no more zeros (j==0) or we find a zero
        % in an unassigned column (c(j)==0).

        while (Opt(j)~=0)
                % Advance lj and j in zero list.
                lj=j;
                j=-InMatrix(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.
                Opt(j)=i;

                % Remove InMatrix(i,j) from unassigned zero list.
                InMatrix(i,lj)=InMatrix(i,j);

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

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

                % Check all zeros in this row.

                lj=n+1;
                j=-InMatrix(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=Opt(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 (Opt(m)==0)
                                        break;
                                end

                                % Advance one step in list.
                                lm=m;
                                m=-InMatrix(r,lm);
                        end

                        if (m==0)
                                % We failed on row r. Continue with next zero on row i.
                                lj=j;
                                j=-InMatrix(i,lj);
                        else
                                % We found a zero in an unassigned column.

                                % Replace zero at (r,m) in unassigned list with zero at (r,j)
                                InMatrix(r,lm)=-j;
                                InMatrix(r,j)=InMatrix(r,m);

                                % Update last/next pointers in row r.
                                NZ(r)=-InMatrix(r,m);
                                LZ(r)=j;

                                % Mark InMatrix(r,m) as an assigned zero in the matrix . . .
                                InMatrix(r,m)=0;

                                % ...and in the assignment vector.
                                Opt(m)=r;

                                % Remove InMatrix(i,j) from unassigned list.
                                InMatrix(i,lj)=InMatrix(i,j);

                                % Update last/next pointers in row r.
                                NZ(i)=-InMatrix(i,j);
                                LZ(i)=lj;

                                % Mark InMatrix(r,m) as an assigned zero in the matrix . . .
                                InMatrix(i,j)=0;

                                % ...and in the assignment vector.
                                Opt(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=Opt(Opt~=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 [InMatrix,Opt,U]=hmflip(InMatrix,Opt,LC,LR,U,l,r)
%HMFLIP Flip assignment state of all zeros along a path.
%
%[InMatrix,Opt,U]=hmflip(InMatrix,Opt,LC,LR,U,l,r)
%Input:
%InMatrix   - the Opt matrix.
%Opt   - the assignment vector.
%LC  - the column label vector.
%LR  - the row label vector.
%U   - the
%r,l - position of last zero in path.
%Output:
%InMatrix   - updated Opt matrix.
%Opt   - updated assignment vector.
%U   - updated unassigned row list vector.

n=size(InMatrix,1);

while (1)
    % Move assignment in column l to row r.
    Opt(l)=r;

    % Find zero to be removed from zero list..

    % Find zero before this.
    m= InMatrix(r,:)==-l;

    % Link past this zero.
    InMatrix(r,m)=InMatrix(r,l);

    InMatrix(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.
        InMatrix(r,l)=InMatrix(r,n+1);
        InMatrix(r,n+1)=-l;

        % Continue back along the column to get row of next zero in path.
        r=LC(l);
    end
end


function [InMatrix,CH,RH]=hmreduce(InMatrix,CH,RH,LC,LR,~,SLR)
%HMREDUCE Reduce parts of Opt matrix in the Hungerian method.
%
%[InMatrix,CH,RH]=hmreduce(InMatrix,CH,RH,LC,LR,SLC,SLR)
%Input:
%InMatrix   - 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:
%InMatrix   - Reduced Opt matrix.
%CH  - Updated vector of 'next zeros' in each row.
%RH  - Updated vector of unexplored rows.

n=size(InMatrix,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(InMatrix(r,c)));

% Subtract minimum from all uncovered elements.
InMatrix(r,c)=InMatrix(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 (InMatrix(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=InMatrix(i,:);
            colsInList=-row(row<0);
            if (isempty(colsInList))
                % No zeros in the list.
                l=n+1;
            else
                l=colsInList(row(colsInList)==0);
            end
            % Append this zero to end of list.
            InMatrix(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(InMatrix(r,c)<=0);

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

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

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