function solved = sudoku_solve(default)
%Solve a sudoku puzzle.  The puzzle is specified as a 9x9 matrix, with 0
%for a blank space, and numbers specifying the initial conditions of the
%puzzle.  We solve the puzzle using 3 methods:
%     1) assign values to squares that can take only one value, then
%        eliminate that possibility from the row, column and region
%     2) eliminate in reverse.  see which values in a row/column/region can
%        only be fixed to a certain place and eliminate those possibilities
%        in other squares
%     3) look for pairs/triples that force numbers into one of two or 3
%        squares, even though we can't decide which.  that helps eliminate
%        possible values from a row/col/region, even though it doesn't
%        necessarily lead to an assignment, but helps eliminate values for
%        future iterations.
% We stop when either all squares are filled, or the board doesn't change
% twice in a row.  When the latter case happens, it means that we can only
% assign values at random and hope that we find a consistent (although not
% necessarily unique) solution.  we could do that using A* searching
% algorithm to try different combinations, but that will be left for
% another release.  In this case, we will fail gracefully (e.g. stop
% running) leaving some squares unfilled.  You can see an examples of this
% by uncommenting the 12/18 puzzle.
%
% The sudoku puzzle will be solved in 1-1.5 seconds.  We test to make sure
% the solution is a valid one.
%
% Code by Rob Turetsky rob@ee.columbia.edu, who allowed this to frustrate
% me instead of working on stuff that i should have been.  I did not
% optimize the code in any way, so feel free to modify it on your own.

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

%test puzzles - you can uncomment one at a time or pass in your own

% if ~exist('default')
%     default = [ 0 0 0 0 8 5 0 0 0; ... 
%                 0 0 4 0 0 0 0 0 7; ...
%                 9 0 6 4 0 0 0 0 0; ...
%                 0 0 7 0 5 0 0 6 3; ...
%                 4 0 2 0 0 0 7 0 9; ...
%                 6 3 0 0 2 0 4 0 0; ...
%                 0 0 0 0 0 7 8 0 5; ...
%                 8 0 0 0 0 0 6 0 0; ...
%                 0 0 0 9 3 0 0 0 0];
% end

if ~exist('default')
    default = [ 9 0 2 0 0 0 6 0 0; ... 
                0 0 6 8 2 3 0 0 0; ...
                1 0 0 0 0 0 2 0 0; ...
                5 0 7 0 0 0 0 0 6; ...
                8 0 0 4 0 5 0 0 2; ...
                2 0 0 0 0 0 1 0 9; ...
                0 0 5 0 0 0 0 0 4; ...
                0 0 0 7 6 2 8 0 0; ...
                0 0 8 0 0 0 9 0 3];
end

%12/18/2005  http://www.msnbc.com/comics/games/sudoku.asp?sFile=sudoc051218
% compy has no problems with this one
% if ~exist('default')
%     default = [ 2 0 0 0 0 1 7 0 9; ...
%                 0 6 0 0 0 0 8 0 0; ...
%                 0 0 0 0 3 7 4 0 0; ...
%                 0 0 0 0 0 0 1 3 0; ...
%                 0 0 0 2 9 8 0 0 0; ...
%                 0 5 4 0 0 0 0 0 0; ...
%                 0 0 8 5 6 0 0 0 0; ...
%                 0 0 6 0 0 0 0 7 0; ...
%                 9 0 3 7 0 0 0 0 1];
% end

%12/11/2005   http://www.msnbc.com/comics/games/sudoku.asp?sFile=sudoc051211
% if ~exist('default')
%     default = [ 0 0 2 7 0 9 3 0 0; ...
%                 6 0 0 0 0 0 0 0 0; ...
%                 3 0 8 0 0 0 0 5 0; ...
%                 1 0 0 0 3 0 0 0 2; ...
%                 9 0 0 0 5 0 0 0 7; ...
%                 2 0 0 0 4 0 0 0 6; ...
%                 0 4 0 0 0 0 8 0 9; ...
%                 0 0 0 0 0 0 0 0 3; ...
%                 0 0 9 3 0 2 1 0 0];
% end

%12/10/2005  http://www.msnbc.com/comics/games/sudoku.asp?sFile=sudoc051210
% this one cannot be solved with logic (i think) - you need to randomly
% fill in pieces.  ouch.  now, we fail gracefully, but in theory i could
% add in an A* search to test possibilities, but i don't want to do that
% if ~exist('default')
%     default = [ 0 0 0 9 4 6 0 0 0; ...
%                 7 0 0 5 0 2 0 3 6; ...
%                 6 0 0 0 0 0 0 5 9; ...
%                 0 0 0 0 0 0 0 0 1; ...
%                 9 0 0 3 0 8 0 0 7; ...
%                 2 0 0 0 0 0 0 0 0; ...
%                 8 9 0 0 0 0 0 0 5; ...
%                 4 6 0 8 0 9 0 0 3; ...
%                 0 0 0 1 7 4 0 0 0];
% end

%global
global numstr;
numstr = '123456789';

%create the board
board = repmat(struct('p1',1,'p2',1,'p3',1,'p4',1,'p5',1,'p6',1,'p7',1,'p8',1,'p9',1,'set',0), 9, 9);
last_board = board;

%add in the defaults
'laying out initial conditions'
for i=1:9
    for j=1:9
        if default(i,j) ~= 0
            board = set_num(board, i, j, default(i,j));
        end
    end
end

'scanning the board'
hit = 1; iter = 1;
while hit > 0 & board_changes(board, last_board)

    last_board = board;
    
    %scan board for singles - the easiest method for solving
    hit = 0;
    for i=1:9
        for j=1:9
            if board(i,j).set == 0
                num = is_ready(board(i,j));
                if num > 0
                    ['found ' numstr(num) ' at (' numstr(i) ', ' numstr(j) ')']
                    board = set_num(board, i, j, num);
                    hit = hit + 1;
                end
            end
        end        
    end    
    
    full_hit = 0;
    %try to do this in reverse - a little trickier
    if hit == 0
        ['trying to eliminate in reverse']
        for r=1:27
            mask = get_mask(r);
            miss = find(missing(board, mask));
            count = zeros(1,length(miss));
            place = zeros(1,length(miss));
            [mu,mv] = find(empty(board,mask));
            for i=1:length(miss)
                pv = possible_vector(board(mu(i),mv(i)));
                for j=1:length(miss)
                    if pv(miss(j))
                        count(j) = count(j)+1;
                        place(j) = i;
                    end
                end
            end
            for i=1:length(miss)
                if count(i) == 1
                    board = set_num(board, mu(place(i)), mv(place(i)), miss(i));
                    full_hit = full_hit + 1;
                end                
            end
        end
    end
    
    back_hit = 0;
    if hit == 0 & full_hit == 0
        %lets try 2-5 at a time - the hardest way, does not necessarily
        %yield eliminations.  if so they'll show up on next iteration
        ['trying to eliminate groups']
        for r=1:27
            mask = empty(board, get_mask(r));
            [a,b] = find(mask);
            sm = zeros(length(a), length(a));
            for i=1:length(a)
                for j=setdiff(1:length(a),i)
                    sm(i,j) = is_equal(board(a(i),b(i)), board(a(j),b(j)));
                end
            end
            if any(any(sm))
                for i=2:5
                    [u,v] = find(sm == i);
                    if length(u) == i
                        pv = find(possible_vector(board(a(u(1)),b(v(1)))));
                        board = cancel(board, mask, u, pv);
                        back_hit = back_hit + 1;
                    elseif sum(sum(sm == i)) > i
                        %we have multiples
                        group = zeros(1,length(u)); maxgroup = 0;
                        for j=1:length(u)-1
                            if group(j) == 0
                                maxgroup = maxgroup + 1;                                
                                group(j) = maxgroup;
                            end
                            for k=j+1:length(u)
                                if is_equal(board(a(u(j)),b(v(j))), board(a(u(k)),b(v(k))))
                                    group(k) = group(j);
                                end
                            end                            
                        end
                        for j=1:maxgroup
                            ug = u(group == group(j)); vg = v(group == group(j));
                            pv = find(possible_vector(board(a(ug(1)),b(vg(1)))));
                            board = cancel(board, mask, ug, pv);
                            back_hit = back_hit + 1;
                        end                        
                    end
                end
            end
        end
    end
    
    if full_hit > 0
        hit = full_hit
    elseif back_hit > 0
        hit = back_hit;
    end
    
    ['found ' num2str(hit) ' squares on iteration ' num2str(iter)]
    iter = iter + 1;
end

solved = reshape([board.set], 9, 9);
if sum(sum(solved == 0)) == 0 & validate(board)
    'We solved it correctly'
else
    'There was an error!'
end
display_board(board);
clear numstr;  %remove global variable

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

function mask = get_mask(i)
mask = zeros(9,9);
modk = mod(i,9); modk(modk == 0) = 9;
if i > 0 & i <= 9
    %the row
    mask(i,:) = 1;
elseif i >= 10 & i <= 18
    %the column
    mask(:,modk) = 1;
elseif i >= 19 & i <= 27
    mask = region_mask(modk);
end

function board = cancel(board, mask, c, n)
global numstr;
[a,b] = find(mask);
for i=setdiff(1:length(a), c)
    for j=1:length(n)
        board(a(i),b(i)).(['p' numstr(n(j))]) = 0;
    end
end

function num = is_ready(square)
global numstr;
num = 0;
if square.set == 0
    hit = 0;
    for i=1:9
        if square.(['p' numstr(i)]) == 1
            num = i;
            hit = hit + 1;
        end
    end
    if hit > 1
        num = 0;
    end    
end

function e = is_equal(sq1, sq2)
e = 0;
if sq1.set == 0 & sq2.set == 0
    v1 = possible_vector(sq1);
    v2 = possible_vector(sq2);
    if (length(find(v1)) == length(find(v2))) & (sum(v1 & v2) == length(find(v1)))
        e = length(find(v1));
    end        
end

function emask = empty(board, mask)
emask = zeros(size(mask));
[i,j] = find(mask);
for k=1:length(i)
    emask(i(k),j(k)) = ~board(i(k),j(k)).set;
end

function m = missing(board, mask)
m = ones(1,9);
[a,b] = find(mask);
for i=1:length(a)
    if board(a(i),b(i)).set
        m(board(a(i),b(i)).set) = 0;
    end
end


function v = possible_vector(square)
%return a binary vector of possible numbers
global numstr;
v = zeros(1,9);
for i=1:9
    v(i) = square.(['p' numstr(i)]);
end

function board = set_num(board, x, y, num)
%set x + y to num, get rid of other regions
global numstr;
board(x,y).set = num;
board(x,y).(['p' numstr(num)]) = 1;
for i=setdiff(1:9,num)
    board(x,y).(['p' numstr(i)]) = 0;
end
    
for i=setdiff(1:9,x)
    board(i,y).(['p' numstr(num)]) = 0;
end

for j=setdiff(1:9,y)
    board(x,j).(['p' numstr(num)]) = 0;
end
    
r = region(x,y);
r(x,y) = 0;
[a,b] = find(r);
for i=1:length(a)
    board(a(i),b(i)).(['p' numstr(num)]) = 0;
end

function r = region(i,j)
r = zeros(9,9);
if i >= 1 & i <= 3
    if j >= 1 & j <= 3
        r(1:3,1:3) = 1;
    elseif j >= 4 & j <= 6
        r(1:3,4:6) = 1;
    else
        r(1:3,7:9) = 1;
    end
elseif i >= 4 & i <= 6
    if j >= 1 & j <= 3
        r(4:6,1:3) = 1;
    elseif j >= 4 & j <= 6
        r(4:6,4:6) = 1;
    else
        r(4:6,7:9) = 1;
    end
else
    if j >= 1 & j <= 3
        r(7:9,1:3) = 1;
    elseif j >= 4 & j <= 6
        r(7:9,4:6) = 1;
    else
        r(7:9,7:9) = 1;
    end
end

function emask = empty(board,mask)
emask = zeros(9,9);
[i,j] = find(mask);
for i=1:length(i)
    if board(i,j).set == 0
        emask(i,j) = 1;
    end
end

function display_board(board)
s = repmat('-',1,19);
for i=1:9
    si = '|';
    for j=1:9
        si = [si num2str(board(i,j).set)];
        if j == 3 | j == 6 | j == 9
            si = [si '|'];
        else 
            si = [si ' '];
        end
    end
    s = [s; si];
    if i == 3 | i == 6
        s = [s; '|' repmat('-',1,17) '|'];
    end
end
s = [s; repmat('-',1,19)];
disp(s)

function m = region_mask(i)
m = zeros(9,9);
if i == 1
    m(1:3,1:3) = 1;
elseif i == 2
    m(1:3,4:6) = 1;
elseif i == 3
    m(1:3,7:9) = 1;
elseif i == 4
    m(4:6,1:3) = 1;
elseif i == 5
    m(4:6,4:6) = 1;
elseif i == 6
    m(4:6,7:9) = 1;
elseif i == 7
    m(7:9,1:3) = 1;
elseif i == 8
    m(7:9,4:6) = 1;
elseif i == 9
    m(7:9,7:9) = 1;
end

function c = board_changes(b1, b2)
c = 0;
for i=1:9
    for j=1:9
        v1 = [possible_vector(b1(i,j)) b1(i,j).set];
        v2 = [possible_vector(b2(i,j)) b2(i,j).set];
        if ~all(v1 == v2)
            c = 1;
        end
    end
end

function v = validate(board)
if sum([board.set] == 0) == 0
    v = 1;
else 
    v = -1;
end

for i=1:27
    mask = get_mask(i);
    set = [board(find(mask)).set];
    set = setdiff(set,0);
    if length(unique(set)) ~= length(set)
        v = 0;
    end
end