Contents

% SUDOKU_SOLVE   Will solve sudoku table of arbitrary size.
%   [ solution, stepCount, rating ] = SUDOKU_SOLVE( A )
%   [ solution, stepCount, rating ] = SUDOKU_SOLVE( A, C )
%   [ solution, stepCount, rating ] = SUDOKU_SOLVE( A, C, maxRating )
%   [ solution, stepCount, rating, solutionCount ] = SUDOKU_SOLVE( A )
%   [ solution, stepCount, rating, solutionCount ] = SUDOKU_SOLVE( A, C )
%   [ solution, stepCount, rating, solutionCount ] = SUDOKU_SOLVE( A, C, maxRating )
%
% Syntax:  [ solution, stepCount, rating, solutionCount ] = sudoku_solve( A, C, maxRating )
%
% Inputs:
%    A - Sudoku table to be solved, must be square with edge length that
%    has integer square root. Smalest valid sudoku table is 4x4.
%    C - Index cache table generated by generate_influence_indices function
%    with input parameter being dimension of one block.
%    maxRating - Rating that if reached terminates execution. This rating
%    more or less correlates with execution time requirements. Default
%    value is set to Infinity as search is unbounded.
%
% Outputs:
%    solution - Table with solution or input in case there is none.
%    stepCount - Number of steps taken to reach the result. Step count does
%    not include backtrack operations.
%    rating - Difficulty rating of give sudoku as positive value, higher is
%    more difficult. In case it is -1, given sudoku has no solution. This
%    rating includes backtrack operations. In case it is -2, solver
%    terminated on maxRating.
%    solutionCount - If defined, solver will not terminate on first found
%    solution but will look for all solutions instead. Only first found
%    solution will be returned.
%
% Example 1:
%    load('testData', 'data');
%    solution = sudoku_solve( reshape(data(1,:,:),[9,9]) )
%
% Example 2:
%    load('testData', 'data');
%    C = generate_influence_indices(3);
%    solution = sudoku_solve( reshape(data(1,:,:),[9,9]), C )
%    solution = sudoku_solve( reshape(data(2,:,:),[9,9]), C, 1000 )
%
% Example 3:
%    load('testData', 'data');
%    C = generate_influence_indices(3);
%    [ ~, ~, ~, solutionCount ] = sudoku_solve( reshape(data(1,:,:),[9,9]), C )
%
% Other m-files required: generate_influence_indices.m
% Subfunctions: none
% MAT-files required: none
%
% See also: GENERATE_INFLUENCE_INDICES,  TEST_DATA

Sudoku solver

Function sudoku_solve solves given sudoku table if possible.

Sudoku puzzle

Consider a square table of size $N$ where $N>4$ and $M = \sqrt{N}$ is integer. This table is further divided into square blocks of size $M$. we will now define term \textit{house}: all cells of the table that lie within same row, column, or block belong to the same house. Table is filled with at least 17 numbers from range $[1,M]$ where no two cells with equal value belong to the same house. This table is called sudoku table.

Solver is given sudoku table with task to fill remaining cells with numbers from range $[1,M]$ while not violating initial rule that two cells with same value may not be part of same house.

There is unwritten rule to sudoku tables that says they have to have exactly one correct solution.

Sudoku as CSP

Sudoku puzzle can be, as many others, expressed as Constraint Satisfaction Problem. That means we can map whole puzzle to set of variables that are connected with set of constraints, and each variable has domain of values that can be assigned to it. Domain of each variable depends on other variables and constraints between them. When solving a CSP, we look for such assignment of values to variables that does not violate any constraints and all assigned values are from domains of respective variables.

In our case, each cell in sudoku table is a variable with domain containing values 1 through edge size of the table, while each variable as constraint it must not contain same value as any other variable on same row, column, and in same block.

CSP solver

Process of solving of above defined sudoku CSP follows simple algorithm. First, most constrained variable is chosen, based on current domain sizes, then values from its domain are assigned to it and same step is called with each. After each assignment is done check whether domain of any unassigned variable was reduced to zero, in which case this branch is terminated and we continue with next available value from the domain. In case all values from the domain proved not to be feasible, previous step is reverted.

In case we revert root step, given sudoku has no solution, as the most constrained variable has no assignment that leads to valid table.

After each assignemnt is done check whether maximum amount of variables was assigned, and if so, we terminate with correct solution. This means we find only one solution in case there is more of them.

Iteration through domain of a variable as implemented heuristic that takes the least constraining value first. Least constraining value is such, that reduces domains of variables it has influence over the least. Choosing such value will increase number of possibilities we may have in future increasing the chance one of them will be feasible.

Cached influence indices

When modifying domains of influenced variables after assignment to a variable, we need to select these domains based on dynamic row, column and block selector. To remove number of computations, all possible influence index configurations are precalculated and cached, so we can check all influenced variables with simple [row, col] query on the cache.

Since it is intended to run large number of solve calls in a sequence, this cache can be calculated once and handed to this solver. This approach not only avoids large number of calculations, but also frequent allocations and deallocations.

Complexity rating

When evaluating complexity of sudoku puzzle, we need to consider steps taken by solver to solve it with least amount of backtracking. There are various different ways to rate complexity of sudoku puzzles, mainly based on weighting of different techniques (difficulty to spot) and aggregating this weight as they are applied. (see https://www.fi.muni.cz/~xpelanek/publications/sudoku-arxiv.pdf)

Since implementation of this solver has limited amount of heuristics, rating is calculated from number of steps taken with weight of selection of least constrained (the sorting process) and then penalization for backtrack.

Possible improvements

This implementation, naturally, does not contain all possible optimizations. Here are suggested few improvements that could be done to this function.

Hidden singles assignment

As we know, each block, row, and column has to contain on of each value. That means if a row, block, or column in page for specific value contains only one valid variable, we have to assign that value to that variable no matter the actual size of the domain for that variable. This way can be avoided unnecessary branching and we gain free domain reduction for other variables.

It would be rather expensive to check for each unassigned variable all its values for their uniqueness in tree different spaces. That can be solved by implementation of row, column, and block sums, where we maintain for each value and row/column/block sum of available variables within. This way we only check whether any of these sum matrices contains number 1, and if so, we already know row/column and value of hidden single variable, making calculation of remaining coordinate trivial.

Maintenance of these sum structures is not trivial and requires a lot of thought, mainly because with each assignment is changed domain of assigned variable, sums of all values for that variable, and domains of variables with assigned value available on the row, column and block.

Naked pair domain reduction

When two variables within a block, row, or column have domain size equal to two, and values within their domains are the same, we can dafely assume one has to have have one value, and the other has to have have the other value. That means no other variable within that block, row, or column can have these values assigned, allowing us to remove them from their domains, possibly exposing other clear assignments.

This improvement can be extended up to four variables with no extra computational cost, as we always compare domain arrays as a whole. Checking all pairs, triplets and quadruplets of all variables would be rather demanding, but we can limit our search only to those variables that changed their domain in last assignment.

Enhanced influence index caching

When using influence index cache, we first obtain index within a page to adjust domain size for assigned variable, and then we calculate indices within whole domain matrix by adding multiple of total number of variables. This can be avoided by precalculation of the cache not only for each position within the page, but also for each page.

Allocation reduction

When going up and down the stack during recursion, large number of variables, matrices, and arrays is allocated, and quite unnecessarily. Since we know maximum depth of the recursion at the time of invocation, we can preallocate all these structures before recursion, and then simply address appropriate stack 'context' with recursion depth.

Recursion removal

Recursion can be removed once previous optimization is implemented, as we only increment/decrement stack page pointer within a while loop, exitting either when all is assigned, or stack pointer is decremented to zero. This improvement more than anything else eleviates the restriction on recursion depth count, which can occour on older matlab versions and higher sudoku sizes (sudoku with block size 5 has edge 25, and 625 variables to have assigned, which means 625 recursion steps, 125 above default limit).

Head of function

First argument is mandatory while second serves only as optimization to allow more efficient mass solving when all tables have the same size.

function [ solution, stepCount, rating, solutionCount ] = sudoku_solve( A, C, maxRating )

Check validity of inputs. Expected at least one argument that is square 2D matrix of size larger than 4x4. This matrix has to have integer square root of its edge and contain only integers from range [0, edge_size]. More than two areguments will not be accepted. Second argument is checked only on size, as that is telltale sign that generate_influence_indices call was done with incompatible argument.

if nargin < 1
    return;
end
if isempty(A) || size(A,1) ~= size(A,2)  || length(size(A))~=2
    error('Input table has to be non-empty and square!');
end
if size(A,1)<4
    error('Valid assignment has to be at least 4x4!');
end
if sqrt(size(A,1))*sqrt(size(A,1)) ~= size(A,1)
    error('Edge of the table must have integer square root!');
end
if ~isequal(fix(A),A)
    error('Given table must contain only integer numbers!');
end
if min(min(A)) < 0 || max(max(A)) > size(A,1)
    error('Values in given table must be in range [0, size(A,1)]!');
end
if nargin > 3
    error('Too many input arguments!');
end
if nargin == 3 && maxRating < 0
    error('Maximum rating must be above zero for at least direct solution!');
end
if nargin < 3
    maxRating = Inf;
end

Initialize basic variables and constants.

solution = A;               % Solution matrix for A.
stepCount = 0;              % Number of processing steps it took to solve A.
d = sqrt(size(A,1));        % Size of one block.
r = d*d;                    % Range of values (also edge size of the board).
lr = r+1;                   % Range increased by one as marker for assigned.
t = r*r;                    % Total number of variables in the system.
rating = sum(sum(A>0))-t;   % Difficulty rating of given assignment.
asn = sum(sum(A>0));        % Number of assigned values from assignment.

Initialize structure for solution counting. Count solution flag will tell the algorithm to not terminate on found solution, but instead just increment solution count and record last found solution.

countSolutions = (nargout == 4);
solutionCount = 0;
firstSolution = [];
firstRating = -1;
firstStepCount = 0;

Generate relation indices for all positions in the grid. These indices point at all variables that are influenced by variable of lookup index.

if ~exist('C','var')
    C = generate_influence_indices( d );
elseif ~isequal(size(C),[d*d+2*(r-d),r,r])
    error('Received relation index cache of incorrect dimensions!');
end
[domains, domainSizes] = sudoku_domains( A, C );

Start the recursion within try catch block, as success is reported via exception. In case no exception was thrown, we searched whole tree of possibilities and there is no solution, or solution counting is on.

try
    step();

    if countSolutions && solutionCount>0
        solution = firstSolution;
        rating = firstRating;
        stepCount = firstStepCount;
        return;
    end

    rating = -1;
    return;
catch ME
    if (strcmp(ME.identifier,'Solver:found'))
        return;
    elseif (strcmp(ME.identifier,'Solver:maxOut'))
        rating = -2;
        return;
    end

    % Rethrow caught exception if it is not success report or maxOut report,
    % should not happen once everything is tested and implemented.
    rethrow(ME);

end

One solve step called in recursion. Will pick variable for assignment, sequentially assign values from its domain to it and call itself recursively in between.

    function step()
        stepCount = stepCount + 1;

        % Find what size is the most restricted domain so we can query its
        % position.
        domainSize = min(min(domainSizes));

        % Assignment of value that is the only option is lower than when
        % player has to choose from multiple values.
        rating = rating + domainSize;

        % Check maximum rating end exit early if exceeded.
        if rating > maxRating
            throw(MException('Solver:maxOut','Found a solution!'))
        end

        % Find Most Constrained variable, witch is the one with lowest
        % number of options in its domain. We need only one, and which is
        % arbitrary.
        [mcR, mcC] = find(domainSizes == domainSize, 1);

        % Take domain of picked variable.
        domain = find(domains(mcR, mcC, :));

        % Sort domain by constraining factor to pick the least constraining.
        toCheck = cell(length(domain), 4);
        for index = 1:length(domain)
            value = domain(index);
            rel = C(:, mcR, mcC);
            ser = domains(:, :, value);

            % Get indices to all variables that should have value removed
            % from their domains.
            toCheck{index,1} = value;
            toCheck{index,2} = rel(logical(ser(rel)));
            toCheck{index,3} = toCheck{index,2}+(value-1)*t;
            toCheck{index,4} = length(toCheck{index,3});
        end
        toCheck = sortrows(toCheck, 4);

        % Remove value from domain of changed variable.
        domains(mcR, mcC, domain) = 0;
        domainSizes(mcR, mcC) = lr;

        % Increase number of assigned values as new value will be assigned
        % for following solveStep calls.
        asn = asn + 1;

        % Iterate the domain.
        for index = 1:length(domain)
            value = toCheck{index,1};
            pos = toCheck{index,2};
            posG = toCheck{index,3};

            % Assign value to the solution.
            solution(mcR, mcC) = value;

            % Reduce influenced domains.
            domainSizes(pos) = domainSizes(pos) - 1;
            domains(posG) = 0;

            % Check success conditions (throw exception). If count of
            % assigned values is equal to total numbers to be assigned, we
            % have solution and no more checks are needed.
            if asn == t
                if countSolutions
                    if solutionCount == 0
                        firstSolution = solution;
                        firstRating = rating;
                        firstStepCount = stepCount;
                    end
                    solutionCount = solutionCount + 1;
                else
                    throw(MException('Solver:found','Found a solution!'));
                end
            else
                % Check fail conditions (all domain sizes are non-zero, which
                % means there is either something to assign, thus number of
                % avalable numbers is present, or something already assigned,
                % and value is above maximum number of values in a domain = lr)
                if all(all(domainSizes(pos)))
                    step();
                end
            end

            % Undo reduction of influenced domains.
            domainSizes(pos) = domainSizes(pos) + 1;
            domains(posG) = 1;

            % Player is punished when he has to revert too many reductions.
            rating = rating + numel(pos);

            % Check maximum rating end exit early if exceeded.
            if rating > maxRating
                throw(MException('Solver:maxOut','Found a solution!'));
            end
        end

        % Undo assignment of the value.
        solution(mcR, mcC) = 0;

        % Decrease number of assigned values as we are backtracking from
        % assignment to this variable.
        asn = asn -1;

        % Restore availability of value in domain for assigned variable.
        domains(mcR, mcC, domain) = 1;
        domainSizes(mcR, mcC) = domainSize;
    end
end