function S=sudoku_gurobi(X)
% SEDOKU_GUROBI Solve a sedoku puzzle using Gurobi for an Binary LP
% 
% S=sudoku_gurobi(X)
%   X has a set of initial conditions where anything ~= 1 ... 9 is an 
%   unknown
% S is the solved sudoku puzzle

assert(all(size(X) == [9 9]))

% create variables
nvars = 9*9*9;
vtypes = char('B'*ones(1,nvars));

% index the variables for easier access
vars = zeros(9,9,9);
ind = 1;
for i=1:9
    for j=1:9
        for v=1:9
            vars(i,j,v) = ind;
            ind = ind + 1;
        end
    end
end

% 
obj = ones(1,nvars); % no objective, just constraints


% add known cells as constraints
A = sparse([])
b = [];
for i=1:9
    for j=1:9
        if X(i,j) < 1 || X(i,j) > 9
            % this is a free variable, ignore
        else
            % this is a constraint, add it
            cons = zeros(nvars,1);
            cons(vars(i,j,X(i,j))) = 1;
            A(:,end+1) = sparse(cons);
            b(end+1) = 1;
        end
    end
end

%
% add sudoku constraints
%

% each cell only has one value
for i=1:9
    for j=1:9
        cons = zeros(nvars,1);
        for v=1:9
            cons(vars(i,j,v)) = 1;
        end
        A(:,end+1) = sparse(cons);
        b(:,end+1) = 1;
    end
end


% each row has 1-9
for i=1:9
    for v=1:9
        cons = zeros(nvars,1);
        % the value v only appears in each row once
        for j=1:9
            cons(vars(i,j,v)) = 1;
        end
        A(:,end+1) = sparse(cons);
        b(:,end+1) = 1;
    end
end    

% each column has 1-9
for j=1:9
    for v=1:9
        cons = zeros(nvars,1);
        % the value v only appears in each row once
        for i=1:9
            cons(vars(i,j,v)) = 1;
        end
        A(:,end+1) = sparse(cons);
        b(:,end+1) = 1;
    end
end    

% subgrid has 1-9
for v=1:9
    for i0=0:2
        for j0=0:2
            cons = zeros(nvars,1);
            for i=1:3
                for j=1:3
                    cons(vars(i+3*i0,j+3*j0,v)) = 1;
                end
            end
            A(:,end+1) = sparse(cons);
            b(:,end+1) = 1;
        end
    end
end

% setup Gurobi call
objtype = 1; % 1 - minimize, -1 - maximize
A = A'
lb = [];
ub = [];
contypes = char('='*ones(1,size(A,1)));
clear opts
opts.IterationLimit = 1000;
opts.FeasibilityTol = 1e-6;
opts.IntFeasTol = 1e-5;
opts.OptimalityTol = 1e-6;
opts.Method = 1; % 0 - primal, 1 - dual
opts.Presolve = -1; % -1 - auto, 0 - no, 1 - conserv, 2 - aggressive
opts.Display = 1;
opts.LogFile = 'test_gurobi_mex_MIP.log';
opts.WriteToFile = 'test_gurobi_mex_MIP.mps';
[x,val,exitflag,output] = gurobi_mex(obj,objtype,A,b,contypes,lb,ub,vtypes,opts);

% assign solution
S = zeros(9);
for i=1:9
    for j=1:9
        for v=1:9
            if x(vars(i,j,v)) > 0
                S(i,j) = v;
            end
        end
    end
end