% taille = 4;
% connus = [1,2,2; 1,3,4; 2,4,2; 3,1,3; 4,2,1; 4,3,3];

taille = 9;
connus = [1,2,2;
    1,5,3;
    1,8,4;
    2,1,6;
    2,9,3;
    3,3,4;
    3,7,5;
    4,4,8;
    4,6,6;
    5,1,8;
    5,5,1;
    5,9,6;
    6,4,7;
    6,6,5;
    7,3,7;
    7,7,6;
    8,1,4;
    8,9,8;
    9,2,3;
    9,5,4;
    9,8,2];

nbvariables = taille^3;
nbconditions = (taille^(2))*4+length(connus);

Aeq = zeros(nbconditions,nbvariables);
beq = ones(nbconditions,1);

condrow = 1;
% =========== lignes ===========
for i = 1:taille
    for k = 1:taille
        for j = 1:taille
            Aeq(condrow, (i-1)*taille*taille+(j-1)*taille+k) = 1;
        end
        condrow = condrow + 1;
    end
end
% =========== colonnes ===========
for k = 1:taille
    for j = 1:taille
        for i = 1:taille
            Aeq(condrow, (i-1)*taille*taille+(j-1)*taille+k) = 1;
        end
        condrow = condrow + 1;
    end
end
% =========== carrés ===========
for k = 1:taille
	n = sqrt(taille);
    for u = 0:n-1
        p = u*n+1;
        for v = 0:n-1
            q = v*n+1;
            for i = p:p+n-1
                for j = q:q+n-1
                    Aeq(condrow, (i-1)*taille*taille+(j-1)*taille+k) = 1;
                end
            end
            condrow = condrow + 1;
        end
    end
end

% unicité chiffres
for i = 1:taille
    for j = 1:taille
        for k = 1:taille
            Aeq(condrow, (i-1)*taille*taille+(j-1)*taille+k) = 1;
        end
        condrow = condrow + 1;
    end
end

% Nombres connus
for n = 1:size(connus)
    i = connus(n,1);
    j = connus(n,2);
    k = connus(n,3);
    Aeq(condrow, (i-1)*taille*taille+(j-1)*taille+k) = 1;
    condrow = condrow + 1;
end

% matrice
lb = zeros(nbvariables,1);
ub = ones(nbvariables,1);

f = zeros(nbvariables,1);
intcon = nbvariables;

A = [];
B = [];

x = intlinprog(f,intcon,A,B,Aeq,beq,lb,ub);

grid = zeros(taille);

for i = 1:taille
    for j = 1:taille
        for k = 1:taille
            index = (i-1)*taille*taille + (j-1)*taille + k;
            if int8(x(index)) == 1
                grid(i,j) = k;
                break
            end
        end
    end
end

grid