Winner the cyclist (CT 1)

Finish 2005-11-09 09:00:00 UTC

anneal

by kup

Status: Failed
Results:

Comments
kup
09 Nov 2005
simplified simulated annealing approach with the getStartCombo().
Please login or create a profile.
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;