Winner the cyclist (CT 1)

2005-11-09 09:00:00 UTC

# anneal

by kup

Status: Failed
Results:

kup
09 Nov 2005
simplified simulated annealing approach with the getStartCombo().
Code
```function solution = solver(originalpuzzle,originalnumbers)
% SOLUTION = SOLVER(PUZZLE,LIST) Sudoku solver

% parameters
numcycles = 5;
numiterations = 4000;
level = 1;
damp = 0.98;

changecount = 1;
TAB = sprintf('\t');
NL = sprintf('\n');

% definitions for score calculations
h = reshape(1:9,3,3);
g = ceil((1:9)/3);
h = (h(g,g)-1)*9 + repmat(h,3);
z = cat(2,cat(1,ones(3),2*ones(3),3*ones(3)),3+cat(1,ones(3),2*ones(3),3*ones(3)),6+cat(1,ones(3),2*ones(3),3*ones(3)));

k = 1;
m = size(originalpuzzle,1);
for i=1:m
for j=1:size(originalpuzzle,2)
if originalpuzzle(j,i) == 0
index(k,1) = (i-1)*m+j;     % linear index
index(k,2) = 9 + i;         % row
index(k,3) = 18 + j;        % column
index(k,4) = z(j,i);        % square no.
k = k+1;
end
end
end
numindex = size(index,1);
numnumbers = numel(originalnumbers);

for cycle = 1:numcycles
numbers = originalnumbers;
[aux puzzle score] = getStartCombo(originalpuzzle, numbers, h);

for j=1:numnumbers
for k = 1:numindex
if aux(k) == numbers(j)
numbers(j) = 0;
end
end
end
k = numindex+1;
for i=1:numnumbers
if numbers(i) ~= 0
aux(k) = numbers(i);
k = k+1;
end
end
numbers = aux;

% calculate initial score
sums  = [sum(puzzle(h)) sum(puzzle) sum(puzzle,2)'];
oldscore = sum(abs(mean(sums)-sums));
if cycle == 1
solutionscore = oldscore;
solution = puzzle;
end

n=1;
while n < numiterations
n = n+1;
% swap numbers
i1 = floor(numnumbers * rand)+1;
i2 = floor(numnumbers * rand)+1;
y1 = numbers(i1);
y2 = numbers(i2);
numbers(i1) = y2;
numbers(i2) = y1;

% incremental score calculation
delta = y2-y1;
if i1 <= numindex
sums(index(i1,4)) = sums(index(i1,4)) + delta;
sums(index(i1,2)) = sums(index(i1,2)) + delta;
sums(index(i1,3)) = sums(index(i1,3)) + delta;
end
if i2 <= numindex
sums(index(i2,4)) = sums(index(i2,4)) - delta;
sums(index(i2,2)) = sums(index(i2,2)) - delta;
sums(index(i2,3)) = sums(index(i2,3)) - delta;
end
score = sum(abs(mean(sums)-sums));

if score < oldscore | rand < level
oldscore = score;
changecount = 1;
else
changecount = changecount+1;
numbers(i1) = y1;
numbers(i2) = y2;
if i1 <= numindex
sums(index(i1,4)) = sums(index(i1,4)) - delta;
sums(index(i1,2)) = sums(index(i1,2)) - delta;
sums(index(i1,3)) = sums(index(i1,3)) - delta;
end
if i2 <= numindex
sums(index(i2,4)) = sums(index(i2,4)) + delta;
sums(index(i2,2)) = sums(index(i2,2)) + delta;
sums(index(i2,3)) = sums(index(i2,3)) + delta;
end
end
level = level * damp;
end % iterations

if solutionscore >= oldscore
for k=1:numindex
puzzle(index(k,1)) = numbers(k);
end
solution = puzzle;
solutionscore = oldscore;
end

end % cycle

% #########################################################################
function [bestx A, G] = getStartCombo(A,list,h)

% nloop, greedyfactor
p = [10, 1.5];

bestcost = inf;
bestx = [];

% generate system matrix
R = repmat(eye(9),1,9); % row sum
C = kron(eye(9),ones(1,9)); % column sum
% region sum
B = repmat(kron(eye(3),ones(1,3)),1,3);
B = blkdiag(B,B,B);

% system matrix
MF = [R; C; B] - 1/9;
b = -MF*A(:);

% pick the free variables
M = MF(:,A(:)==0);

% try to solve
% M*x = b
% using as x only pertubations of a given list

% solve the relaxed problem
% (bias so that |x_acc - mean(list)| is minimized in 2-norm)
pinvM = pinv(M);
x_acc = pinvM*(b - sum(M,2)*median(list)) + median(list);

% try to find feasible solution close to the x_acc
x = fittosol(x_acc, list);
f = b - M*x;

for i = 1:p(1)

%Gauss-Newton
dx = p(2)*pinvM*f;
x = fittosol(x+dx, list);
f = b-M*x;

if sum(abs(f)) < bestcost
bestx = x;
bestcost = sum(abs(f));
if bestcost==0 % optimum
break
end
end

end

% substitute solution to A
A(A(:)==0) = bestx;

sums  = [sum(A(h)) sum(A) sum(A,2)'];
mu_sums = sum(sums)/length(sums);
G = sum(abs(mu_sums-sums));

%%%%%
% return y such that y is a (partial) permutation of list
% and |y-x| is minimized
% (assume list is sorted)

function x = fittosol(x, list)

extra = length(list)-length(x);
for i=1:extra
% random element introduced here!
list(ceil(length(list)*rand)) = [];
end
[sortedx, si] = sort(x);
x(si) = list;

```