# MATLAB: Knuth’s Algorithm X (Dancing Links; Exact Cover)

I'm trying to write a Matlab version of Knuth's "Algorithm X." This algorithm makes use of the Dancing Links routine and is a recursive, nondeterministic, depth-first, backtracking algorithm. I am using this algorithm to solve a simple application of the exact cover problem. For more information:
Here is pseudocode for Algorithm X:
`` 1. If matrix A is empty, the problem is solved; terminate successfully. 2. Otherwise: choose a column, c, of A with the least number of 1's. 3. Choose a row, r, that contains a 1 in column c (e.g., A(r,c) = 1). 4. Add r to the partial solution. 5. For each column j such that A(r,j) = 1:         for each row i such that A(i,j) = 1:              delete row i from matrix A         delete column j from matrix A 6. Repeat this algorithm recursively on the reduced matrix A.``
Following the Wiki-page, here is a working example:
``A = zeros(1,7); B = A; C = B; D = C; E = D; F = E;A([1 4 7]) = 1; B([1 4]) = 1; C([4 5 7]) = 1; D([3 5 6]) = 1; E([2 3 6 7]) = 1; F([2 7]) = 1; M = [A; B; C; D; E; F]; % the binary 'A' matrix``
Knuth's algorithm should return the following exact cover of M: {B, D, F}; i.e., rows 2, 4, and 6 of the original (unreduced) matrix collectively contribute a 1 in each column of M (exactly once).
Here is what I have so far:
``function sol = dlx(X,sol)% if isempty(X)%     return% end% Find the first column of X having the lowest number of 1s in any positionc = find(sum(X,1) == min(sum(X,1)), 1);% If column c is entirely zero, terminate unsuccessfully% if (sum(X(:,c)) == 0)%     sol = [];%     return% end% Find the rows that have a 1 in column cr = find(X(:,c) == 1);% if isempty(r)%     sol = [];%     return% end% if isempty(r)%     sol = sol(1:end-1);%     return% end% Loop over each row that has a 1 in column cfor rr = 1:length(r)      % include row r in the partial solution      sol = [sol r(rr)];      % find the columns in which A(r,j) = 1      J = find(X(r(rr),:) == 1);      % initialize the reduced matrix      Xred = X;      I = [];      for jj = 1:length(J)          I = [I; find(Xred(:,J(jj)) == 1)];      end      I = unique(I);      Xred(I,:) = []; % delete the I rows      Xred(:,J) = []; % delete the J columns      % repeat this algorithm recursively on the reduced matrix Xred.      sol = dlx(Xred, sol);endend``
For the above example,
``sol = dlx(X, [])sol =       1     2     1     1``
Except that my algorithm should return sol = [2 4 6].
Here is another simple example that also bugs my code:
``A = [1 1 0; 0 1 1];sol = dlx(A, [])``
sol =
``     1``
In this case, the there is no exact cover of the matrix A, and so my algorithm should return an empty solution (e.g., sol = []).
As you can see, I am very close. I believe that I have successfully written the recursion of reducing the matrix, but I am seeking help on the following:
1. Looking at the pseudocode, should the "Choose a row, r, that contains a 1 in column c" be implemented as a for-loop in my code? That is, should I loop over all such rows?
2. I'm having trouble formulating the final solution from the partial solution (e.g., so that my function returns sol = [2 4 6]). During the recursion, the partial solution should 'backtrack' when the partial solution turns out not to be valid.
3. Where should I place the 'if A is empty' condition from the pseudocode? You'll see that I've tried putting this condition at the beginning of the function, and also at the 'isempty(r)' line.
4. I'm also open to suggestions for making my code more efficient; i.e., should/can I replace the 'find' or 'for-loops' with better code?
Can any one offer me suggestions to my questions? Thanks in advance!

``function [sol, success] = dlx(A)    sol = recursion(A, 1:size(A, 1));endfunction [sol] = recursion(A, rownum)    sol = {};   %empty cell array indicates failure       %step 1    if isempty(A)        %success        sol = {[]};        return;    end    %step 2    [~, c] = min(sum(A, 1));    if ~any(A(:, c))        %failure        return;    end    %step 3    possiblerows = find(A(:, c));    for tryrow = possiblerows.'        %step 4 for later, only if success        %step 5        deletecols = A(tryrow, :) == 1;        deleterows = any(A(:, deletecols) == 1, 2);        reducedA = A(:,~deletecols);        reducedA = reducedA(~deleterows, :);        reducedrows = rownum(~deleterows);        %step 6        [subsol] = recursion(reducedA, reducedrows);        if ~isempty(subsol)            %step 4             %add current row to the all the row combinations returned by the recursion            sol = [sol; cellfun(@(ss) [rownum(tryrow), ss], subsol, 'UniformOutput', false)]; %#ok<AGROW>        end    end    %if all the tryrow ended in failure then sol is still an empty cell array which indicates failure on returnend``