function [solution] = solver(puzzle,list)
rand('seed',1);
list=list(:,end:-1:1);
global h iX iN col_selection blockmap rowmap colmap;
h = reshape(1:9,3,3);
g = ceil((1:9)/3);
h = uint8((h(g,g)-1)*9 + repmat(h,3));
h1 = reshape(1:81,9,9);
iX = uint8([h h1 h1']);
iN = ones(81,81)*3;
for rcgi=1:27
for ci=1:9
for ci2=ci+1:9
index1=iX(ci,rcgi);index2=iX(ci2,rcgi);
iN(index1,index2)=iN(index1,index2)-1;
end
end
end
for ci=1:81
for ci2=ci+1:81
iN(ci2,ci)=iN(ci,ci2);
end
end
col_selection = [[4 5 6 7 8 9];[7 8 9 7 8 9]];
blockmap=kron(reshape(1:9,3,3),ones(3));
rowmap = repmat((1:9)',1,9);
colmap = repmat((1:9),9,1);
[solution,s] =solverC(puzzle,list);
if floor(s) == 19; [sol2,s2]=solverC(puzzle(:,end:-1:1),list); if s2<s;solution=sol2(:,end:-1:1);s=s2;end;end;
if s > 200; [sol2,s2]=solverC(puzzle(end:-1:1,:),list); if s2<s;solution=sol2(end:-1:1,:);s=s2;end;end;
[solution1,s2]=improver3(puzzle,list,solution); if s2<s;solution=solution1; s=s2;end;
if s>100 [solution1,s2]=improver3(puzzle,list,solution); if s2<s;solution=solution1; s=s2;end;end;
flag=1;while flag; [solution,flag]=improver4(puzzle,list,solution);end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [solution,s] = solverC(puzzle,list)
NUMREPS = 25000;
NUMTRIES1 = 200;
NUMTRIES2 = 200;
NUMTRIES3 = 300;
global h iX iN col_selection blockmap rowmap colmap;
free = puzzle==0;
q = uint8(find(free)');
missing = numel(q);
lq = missing;
list_count = numel(list);
listO=list;
NUMRUNS = (list_count / 10)+10;
% Calculate the average of all values available
total_puzzle = sum(sum(puzzle));
total_all = sum(list) + total_puzzle;
total_count = (81 - missing) + list_count;
average = total_all / total_count;
target = average*9;
bestscore2 = 1e10;
bestscore = 1e10;
reruns=0;
index1v = ceil(rand(1,NUMREPS)*lq);
index1v = q(index1v);
index2v = mod((index1v+uint8(floor(rand(1,NUMREPS)*(lq-1)))),lq)+1;
index2v = q(index2v);
while reruns<3
for run = 1:NUMRUNS
% Fill in values to optimize the rows
list = listO;
list_left = list_count;
solution = puzzle;
% Fill in values to optimize the rows
for r = 1:9
% Fill in values for all empty spaces in this row
for c = 1:9
if free(r,c)
list_index = ceil(rand*list_left);
solution(r,c) = list(list_index);
list(list_index)=[];
list_left = list_left-1;
end
end
list(end+1:list_count) = 1e10;
if list_left
% Get row as close as possible to
% target value
row_total = sum(solution(r,:)) - target;
est_error = abs(row_total);
tq = r + 9*(find(free(r,:))-1);
indv = tq(ceil(rand(NUMTRIES1,1)*numel(tq)));
for tries = 1:NUMTRIES1
ind = indv(tries);
list_index = ceil(rand*list_left);
list_value = list(list_index);
delta = list_value - solution(ind);
new_est_error = abs(row_total + delta);
if (new_est_error < est_error)
row_total = row_total + delta;
est_error = new_est_error;
list(list_index) = solution(ind);
solution(ind) = list_value;
end
end
end
end
target = sum(sum(solution))/9;
% Go about optimizing blocks without changing rows
for rb = 1:3
scale_row = (rb-1)*3;
for cb = 1:2
scale_col = (cb-1)*3;
cols_all_outside = col_selection(cb,:);
block = sum(sum(solution(scale_row+1:rb*3,scale_col+1:cb*3)));
for tries = 1:NUMTRIES2
row = scale_row + ceil(rand*3);
col_inside = scale_col + ceil(rand*3);
if free(row,col_inside)
col_outside = cols_all_outside(ceil(rand*6));
if free(row,col_outside)
val_inside = solution(row, col_inside);
val_outside = solution(row, col_outside);
new_block = block - val_inside + val_outside;
if abs(new_block - target) < abs(block - target)
solution(row,col_outside) = val_inside;
solution(row,col_inside) = val_outside;
block = new_block;
end
end
end
end
end
end
%Solve the columns without changing rows or blocks
col_total = sum(solution) - target;
for tries = 1:NUMTRIES3
el = q(ceil(rand*missing));
col = colmap(el);
row = el - 9*(col-1);
col_swap = ceil(rand*3) + (ceil(col/3)-1)*3;
if free(row,col_swap)
value_swap = solution(row,col_swap);
error = abs(col_total(col));
value = solution(row,col);
col_swap_total = col_total(col_swap);
col_total_new = col_total(col) + value_swap - value;
col_swap_total_new = col_swap_total + value - value_swap;
error_total_new = abs(col_total_new) + abs(col_swap_total_new);
if (error_total_new < error + abs(col_swap_total))
solution(el) = value_swap;
solution(row, col_swap) = value;
col_total(col_swap) = col_swap_total_new;
col_total(col) = col_total_new;
end
end
end
sums = sum(solution(iX));
score = sum(abs(sum(sums)/27-sums));
if score < bestscore
bestinit = solution;
bestscore = score;
bestlist = list;
end
solution = bestinit;
list = bestlist;
sTemp = sum(solution);
target = sum(sTemp)/9;target3=target*3;
rowtotals = sum(solution,2);
coltotals = sTemp;
blocktotals = sum(solution(h));
if list_left ~= 0
rep = 0;
for repl = 1:5000 % blimey! loop unroll! (drumroll?)
% ----- loop 0, do the first thing (loops 0-3 are the same)
rep = rep + 1;
index1 = index1v(rep); index2 = index2v(rep);
index1r = rowmap(index1); index2r = rowmap(index2);
index1c = colmap(index1); index2c = colmap(index2);
value1 = solution(index1); value2 = solution(index2);
block1i = blockmap(index1);block2i = blockmap(index2);
valdiff = value1-value2;
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
subtotal2_old = rowtotals(index2r) + coltotals(index2c) + blocktotals(block2i);
subscoret_old = abs(subtotal1_old - target3) + abs(subtotal2_old - target3);
subscoret_new = abs(subtotal1_old - target3-iN(index1,index2)*valdiff) + abs(subtotal2_old - target3+iN(index1,index2)*valdiff);
if (subscoret_new - subscoret_old) < (rand/6)
solution(index1) = value2;
solution(index2) = value1;
% dont bother to compare if the rows/cols are the same..
rowtotals(index1r) = rowtotals(index1r) - valdiff;
rowtotals(index2r) = rowtotals(index2r) + valdiff;
coltotals(index1c) = coltotals(index1c) - valdiff;
coltotals(index2c) = coltotals(index2c) + valdiff;
blocktotals(block1i) = blocktotals(block1i) - valdiff;
blocktotals(block2i) = blocktotals(block2i) + valdiff;
end
% ----- loop 1, do the first thing
rep = rep + 1;
index1 = index1v(rep); index2 = index2v(rep);
index1r = rowmap(index1); index2r = rowmap(index2);
index1c = colmap(index1); index2c = colmap(index2);
value1 = solution(index1); value2 = solution(index2);
block1i = blockmap(index1);block2i = blockmap(index2);
valdiff = value1-value2;
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
subtotal2_old = rowtotals(index2r) + coltotals(index2c) + blocktotals(block2i);
subscoret_old = abs(subtotal1_old - target3) + abs(subtotal2_old - target3);
subscoret_new = abs(subtotal1_old - target3-iN(index1,index2)*valdiff) + abs(subtotal2_old - target3+iN(index1,index2)*valdiff);
if (subscoret_new - subscoret_old) < (rand/6)
solution(index1) = value2;
solution(index2) = value1;
% dont bother to compare if the rows/cols are the same..
rowtotals(index1r) = rowtotals(index1r) - valdiff;
rowtotals(index2r) = rowtotals(index2r) + valdiff;
coltotals(index1c) = coltotals(index1c) - valdiff;
coltotals(index2c) = coltotals(index2c) + valdiff;
blocktotals(block1i) = blocktotals(block1i) - valdiff;
blocktotals(block2i) = blocktotals(block2i) + valdiff;
end
% ----- loop 2, do the first thing
rep = rep + 1;
index1 = index1v(rep); index2 = index2v(rep);
index1r = rowmap(index1); index2r = rowmap(index2);
index1c = colmap(index1); index2c = colmap(index2);
value1 = solution(index1); value2 = solution(index2);
block1i = blockmap(index1);block2i = blockmap(index2);
valdiff = value1-value2;
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
subtotal2_old = rowtotals(index2r) + coltotals(index2c) + blocktotals(block2i);
subscoret_old = abs(subtotal1_old - target3) + abs(subtotal2_old - target3);
subscoret_new = abs(subtotal1_old - target3-iN(index1,index2)*valdiff) + abs(subtotal2_old - target3+iN(index1,index2)*valdiff);
if (subscoret_new - subscoret_old) < (rand/6)
solution(index1) = value2;
solution(index2) = value1;
% dont bother to compare if the rows/cols are the same..
rowtotals(index1r) = rowtotals(index1r) - valdiff;
rowtotals(index2r) = rowtotals(index2r) + valdiff;
coltotals(index1c) = coltotals(index1c) - valdiff;
coltotals(index2c) = coltotals(index2c) + valdiff;
blocktotals(block1i) = blocktotals(block1i) - valdiff;
blocktotals(block2i) = blocktotals(block2i) + valdiff;
end
% ----- loop 3, do the first thing
rep = rep + 1;
index1 = index1v(rep); index2 = index2v(rep);
index1r = rowmap(index1); index2r = rowmap(index2);
index1c = colmap(index1); index2c = colmap(index2);
value1 = solution(index1); value2 = solution(index2);
block1i = blockmap(index1);block2i = blockmap(index2);
valdiff = value1-value2;
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
subtotal2_old = rowtotals(index2r) + coltotals(index2c) + blocktotals(block2i);
subscoret_old = abs(subtotal1_old - target3) + abs(subtotal2_old - target3);
subscoret_new = abs(subtotal1_old - target3-iN(index1,index2)*valdiff) + abs(subtotal2_old - target3+iN(index1,index2)*valdiff);
if (subscoret_new - subscoret_old) < (rand/6)
solution(index1) = value2;
solution(index2) = value1;
% dont bother to compare if the rows/cols are the same..
rowtotals(index1r) = rowtotals(index1r) - valdiff;
rowtotals(index2r) = rowtotals(index2r) + valdiff;
coltotals(index1c) = coltotals(index1c) - valdiff;
coltotals(index2c) = coltotals(index2c) + valdiff;
blocktotals(block1i) = blocktotals(block1i) - valdiff;
blocktotals(block2i) = blocktotals(block2i) + valdiff;
end
%------ loop 4, do the other thing
rep = rep+1;
index1 = index1v(rep);
index1r = rowmap(index1);
index1c = colmap(index1);
value1 = solution(index1);
block1i = blockmap(index1);
list_index = ceil(rand*list_left);
replace = list(list_index);
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
newtarg = target + (replace-value1)/9;
score_diff = abs(subtotal1_old+3*(replace-value1-newtarg)) - abs(subtotal1_old - target3);
if score_diff < 0.93*rand
solution(index1) = replace;
list(list_index) = value1;
delta=replace-value1;
rowtotals(index1r) = rowtotals(index1r) + delta;
coltotals(index1c) = coltotals(index1c) + delta;
blocktotals(block1i) = blocktotals(block1i) + delta;
target = newtarg;target3=target*3;
end
% ---- end loop
end
else % (list_left == 0)
% time for fresh random numbers...
index1v = ceil(rand(NUMREPS,1)*lq);
index1v = q(index1v);
for rep = 1:NUMREPS
% --- this is same as the "first thing" above
index1 = index1v(rep); index2 = index2v(rep);
index1r = rowmap(index1); index2r = rowmap(index2);
index1c = colmap(index1); index2c = colmap(index2);
value1 = solution(index1); value2 = solution(index2);
block1i = blockmap(index1);block2i = blockmap(index2);
valdiff = value1-value2;
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
subtotal2_old = rowtotals(index2r) + coltotals(index2c) + blocktotals(block2i);
subscoret_old = abs(subtotal1_old - target3) + abs(subtotal2_old - target3);
subscoret_new = abs(subtotal1_old - target3-iN(index1,index2)*valdiff) + abs(subtotal2_old - target3+iN(index1,index2)*valdiff);
if (subscoret_new - subscoret_old) < (rand/6)
solution(index1) = value2;
solution(index2) = value1;
% dont bother to compare if the rows/cols are the same..
rowtotals(index1r) = rowtotals(index1r) - valdiff;
rowtotals(index2r) = rowtotals(index2r) + valdiff;
coltotals(index1c) = coltotals(index1c) - valdiff;
coltotals(index2c) = coltotals(index2c) + valdiff;
blocktotals(block1i) = blocktotals(block1i) - valdiff;
blocktotals(block2i) = blocktotals(block2i) + valdiff;
end
end
end
% --- end the inner loop
sums = sum(solution(iX));
score = sum(abs(sum(sums)/27-sums));
if score < bestscore2
bestsol = solution;
bestscore2 = score;
end
if bestscore2 < 5
break
end
end
if ( reruns==0 & bestscore2>9 )
reruns=reruns+1;
elseif ( reruns==1 & bestscore2>25 )
reruns=reruns+1;
else
reruns=1e8;
end
end
solution = bestsol;
% --- prepare for final improvers
n=ceil(bestscore2);
X = solution(iX);
sX = sum(X);
mX = sum(sX)/27;
XX=zeros(81,3);
for i=q
XX(i,:)=find(sum(iX==i));
end
bv=1;
for i=1:n
bVerb=0;
for iEl=q
I=XX(iEl,:);
s1=sum(sX(I)-mX);
if (s1>0&&sum(sX(I)>mX)<2)||(s1<0&&sum(sX(I)<mX)<2)
continue
end
x=solution(iEl);
dx=solution(q)-x;
if s1>0
j=find(dx<0.8&dx+s1*2>0&dx~=0);
else
j=find(dx>-0.8&dx+s1*2<0&dx~=0);
end
if ~isempty(j)
S0=sX-mX;
s=sum(abs(S0));
m=0;
for k=1:numel(j)
S=S0;
S(I)=S(I)+dx(j(k));
J=XX(q(j(k)),:);
S(J)=S(J)-dx(j(k));
T=sum(abs(S));
if T<s
m=k;
s=T;
end
end
if m
j=q(j(m));
solution(iEl)=solution(j);
solution(j)=x;
X=solution(iX);
sX=sum(X);
bVerb=1;
bv=1;
end
end
end
if ~bVerb
break
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [out,vprev]=improver3(in,list,solution)
%Initialize
idx=find(in==0); a=solution(idx);
a=[a ; list(~ismembc(list,sort(a)))];
N1=numel(idx);
% Initialise the square indices
global h;
isquare=h';
square=kron(reshape(1:9,3,3),ones(3));
% Initialise the solution
[x,y]=find(in==0);
out=solution;
N=numel(a)-N1;
% Calculate bounds on the block values
rowsum=sum(out,2);
optsum=sum(rowsum)/9;
colsum=sum(out,1)';
sqsum=sum(out(isquare),2);
S=[rowsum(x) colsum(y) sqsum(square(idx))]-optsum;
vprev=Inf;
ii=1; cnt=1;
while 1 % Approximate the change for swapping with each weight
yy=ii+1;
tmp=S(zeros(numel(a)-ii,1)+ii,:);
tmp=tmp+a(yy:end,[1 1 1])-a(ii);
tmp=sum(abs(tmp),2)-sum(abs(S(ii,:)));
tmp2=S(yy:end,:)-a(yy:numel(idx),[1 1 1])+a(ii);
tmp2=[sum(abs(tmp2),2)-sum(abs(S(yy:end,:)),2) ; zeros(N,1)];
tmp=tmp+tmp2;
[tmp,idx3]=min(tmp);
if (tmp<0) % Replace one list value by another
zz = idx3+ii;
tmp=a(zz); a(zz)=a(ii); a(ii)=tmp;
out(x(ii),y(ii))=a(ii);
if (zz<=size(S,1)) % Swaps the two currently used list values
out(x(zz),y(zz))=a(zz);
else % Re-calculate the optimal sum
optsum=optsum+(a(ii)-a(zz))/9;
end % Re-calculate the various sums
rowsum=sum(out,2);
colsum=sum(out,1)';
sqsum=sum(out(isquare),2);
S=[rowsum(x) colsum(y) sqsum(square(idx))]-optsum;
end
if (ii==1)
sums = [rowsum ; colsum ; sqsum];
v = sum(abs(sum(sums)/numel(sums)-sums));
if (v>=vprev)
break
end
vprev=v;
end
ii=yy;
if (ii==size(S,1));ii=1;cnt=cnt+1;if (cnt==4);break;end;end;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All of the previous algorithms have focussed on pairwise switching
% of weights. The main loop is a fast random method for choosing the
% blocks to be switched. What follows that (formerly improver2) is a
% slower method comparing the scores before and after each pair of blocks
% already in the matrix have been swapped. Then improver3 is a faster
% method of doing the same thing, but only approximates the change in
% score, and also allows swaps with unassigned elements in the list.
% This improver looks at swapping three elements at a time (although
% not exhaustively), and can trivially be altered for more than three
% elements (increasing to more than 3 would probably minimally affect
% the score though).
function [out,flag]=improver4(in,list,solution)
global iX;
%Initialize
idx=find(in==0); a=solution(idx); N=numel(a);
a=[a ; list(~ismembc(list,sort(a)))];
sums = sum(solution(iX));
bsc = sum(abs(sum(sums)/27-sums)); bsc2=bsc;
ba = a;
[b,idx2]=sort(a);
zok=length(b)-2;
for ii=1:zok
tmp=a(idx2(ii:ii+2));
a=ba;
a(idx2(ii:ii+2))=[tmp(2:3) ; tmp(1)];
solution(idx)=a(1:N);
sums = sum(solution(iX));
score = sum(abs(sum(sums)/27-sums));
if (score<bsc) bsc=score; ba=a; continue; end;
tmp=[tmp(2:3) ; tmp(1)];
a(idx2(ii:ii+2))=[tmp(2:3) ; tmp(1)];
solution(idx)=a(1:N);
sums = sum(solution(iX));
score = sum(abs(sum(sums)/27-sums));
if (score<bsc) bsc=score; ba=a; continue; end;
a(idx2(ii:ii+2))=tmp([3;1;2]);
end;
flag=(bsc<bsc2);
out=in; out(idx)=ba(1:N);
|
