%Based on Hoffa
function solution = solver(puzzle,list)
% SOLUTION = SOLVER(PUZZLE,LIST) Sudoku solver
NUMREPS=10000;
NUMRUNS=11;
h = reshape(1:9,3,3);
g = ceil((1:9)/3);
h = (h(g,g)-1)*9 + repmat(h,3);
blocks = ceil(h/9);
h1=reshape(1:81,9,9);
iX=[h h1 h1'];
rand('seed',1);
q = find(~puzzle);
missing = numel(q);
listO=list;
solution = puzzle;
% Calculate the average of all values available
total_puzzle = sum(sum(puzzle));
total_all = sum(list) + total_puzzle;
total_count = (81 - missing) + numel(list);
average = total_all / total_count;
target = average*9;
bestscore2 = 1e10;
bestscore=1e10;
block_errors = zeros(1,9);
for run= 1:NUMRUNS
% Fill in values to optimize the rows
list=listO;
list_left = numel(list);
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 puzzle(r,c) == 0
list_index = ceil(rand*list_left);
list_value = list(list_index);
list(list_index:list_left-1) = list(list_index+1:list_left);
list(list_left) = 1e10;
list_left = list_left-1;
solution(r,c) = list_value;
end;
end;
if list_left
% Get row as close as possible to target value
row_total = sum(solution(r,:)) - target;
est_error = abs(row_total);
tq=find(puzzle(r,:)==0);
colv=tq(ceil(rand(300,1)*length(tq)));
for tries = 1:300
col = colv(tries);
list_index = ceil(rand*list_left);
list_value = list(list_index);
delta=list_value - solution(r,col);
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(r,col);
solution(r,col) = list_value;
end;
end;
end;
end;
target = sum(sum(solution))/9;
col_selection = [[4 5 6 7 8 9];[7 8 9 7 8 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:200
row = scale_row + ceil(rand*3);
col_inside = scale_col + ceil(rand*3);
if (puzzle(row,col_inside) == 0)
val_inside = solution(row,col_inside);
col_index = ceil(rand*6);
col_outside = cols_all_outside(col_index);
if (puzzle(row,col_outside) == 0)
val_outside = solution(row,col_outside);
error = abs(block - target);
error_new = abs(block + val_outside - val_inside - target);
if (error_new < error)
solution(row,col_outside) = val_inside;
solution(row,col_inside) = val_outside;
block = block - val_inside + val_outside;
end;
end;
end;
end;
end;
end;
index = 1;
for rb = 1:3
scale_row = (rb-1)*3;
for cb = 1:3
scale_col = (cb-1)*3;
block = sum(sum(solution(scale_row+1:rb*3,scale_col+1:cb*3)));
block_errors(index) = block - target;
index = index + 1;
end;
end;
%Solve the columns without changing rows or blocks
col_total = sum(solution)-target;
for tries = 1:1000
col = ceil(rand*9);
error = abs(col_total(col));
row = ceil(rand*9);
value = solution(row,col);
if (puzzle(row,col) == 0)
col_swap = ceil(rand*3) + (ceil(col/3)-1)*3;
value_swap = solution(row,col_swap);
if (puzzle(row,col_swap) == 0)
col_swap_total = col_total(col_swap);
error_swap = abs(col_swap_total);
error_total = error + error_swap;
col_total_new = col_total(col) + value_swap - value;
error_new = abs(col_total_new);
col_swap_total_new = col_swap_total + value - value_swap;
error_swap_new = abs(col_swap_total_new);
error_total_new = error_new + error_swap_new;
if (error_total_new < error_total)
solution(row,col) = value_swap;
solution(row,col_swap) = value;
col_total(col_swap)=col_swap_total_new;
col_total(col)=col_total_new;
end;
end;
end;
end;
sums = sum(solution(iX));
score = sum(abs(sum(sums)/27-sums));
if (score<bestscore)
bestinit=solution;
bestscore=score;
bestlist=list;
end
lq=length(q);
index1v = ceil(rand(NUMREPS,1)*lq);
index2v = mod((index1v+floor(rand(NUMREPS,1)*(lq-1))),lq)+1;
index1v = q(index1v);
index2v = q(index2v);
solution=bestinit;
list=bestlist;
target = sum(sum(solution))/9;
rowtotals=sum(solution,2);
coltotals=sum(solution);
blocktotals=sum(solution(h));
for rep= 1:NUMREPS
index1 = index1v(rep);
index2 = index2v(rep);
index1r = mod(index1-1,9)+1;
index2r = mod(index2-1,9)+1;
index1c = ceil(index1/9);
index2c = ceil(index2/9);
if (rand > 0.8) & (list_left ~= 0)
value1 = solution(index1r,index1c);
block1i = blocks(index1r,index1c);
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
score1_old = abs(subtotal1_old-3*target);
list_index = ceil(rand*list_left);
replace = list(list_index);
score1_new = abs(subtotal1_old+3*replace-3*value1-3*target);
if (score1_new - score1_old) < rand
solution(index1r,index1c) = replace;
list(list_index) = value1;
rowtotals(index1r)=rowtotals(index1r)+replace-value1;
coltotals(index1c)=coltotals(index1c)+replace-value1;
blocktotals(block1i)=blocktotals(block1i)+replace-value1;
target = sum(sum(solution))/9;
end;
else
value1 = solution(index1r,index1c);
value2 = solution(index2r,index2c);
block1i = blocks(index1r,index1c);
block2i = blocks(index2r,index2c);
subtotal1_old = rowtotals(index1r) + coltotals(index1c) + blocktotals(block1i);
subtotal2_old = rowtotals(index2r) + coltotals(index2c) + blocktotals(block2i);
subscore1_old = abs(subtotal1_old-3*target);
subscore2_old = abs(subtotal2_old-3*target);
subscoret_old = subscore1_old + subscore2_old;
delta1 = 3*value1;
delta2 = 3*value2;
subtotal1_new = subtotal1_old - delta1 + delta2;
subtotal2_new = subtotal2_old - delta2 + delta1;
subscore1_new = abs(subtotal1_new-3*target);
subscore2_new = abs(subtotal2_new-3*target);
subscoret_new = subscore1_new + subscore2_new;
if (subscoret_new - subscoret_old) < (rand/6)
solution(index1r,index1c) = value2;
solution(index2r,index2c) = value1;
if (index1r~=index2r)
rowtotals(index1r)=rowtotals(index1r)+value2-value1;
rowtotals(index2r)=rowtotals(index2r)+value1-value2;
end
if (index1c~=index2c)
coltotals(index1c)=coltotals(index1c)+value2-value1;
coltotals(index2c)=coltotals(index2c)+value1-value2;
end
if (block1i~=block2i)
blocktotals(block1i)=blocktotals(block1i)+value2-value1;
blocktotals(block2i)=blocktotals(block2i)+value1-value2;
end
end;
end;
end;
sums = sum(solution(iX));
score = sum(abs(sum(sums)/27-sums));
if (score<bestscore2)
bestsol=solution;
bestscore2=score;
bestlist2=list;
if score<5
return
end
end
end;
%solution=bestsol;
%-----Further Improvement--------------
S_best=bestsol;
Result_best=checkresult(S_best);
solution1=rearrange(puzzle,S_best,1);
solution2=rearrange(puzzle,S_best,2);
Result1=checkresult(solution1);
Result2=checkresult(solution2);
[decideVal,decideIndex]=min([Result_best Result1 Result2]);
switch decideIndex
case 1
solution=S_best;
case 2
solution=solution1;
case 3
solution=solution2;
end
%---------------------------------------------
X=solution(iX);
list=bestlist2(bestlist2<1e6);
sX=sum(X);
mX=sum(sX)/27;
q=q';
[solution,X,sX,mX,list]=improver1(solution,X,sX,mX,list,iX,q,100);
[solution,X,sX,bv]=improver2(solution,X,sX,mX,iX,q,18);
if ~bv;return;end;
[solution,X,sX,mX,list,bv]=improver1(solution,X,sX,mX,list,iX,q,100);
if ~bv;return;end;
[solution,X,sX,bv]=improver2(solution,X,sX,mX,iX,q,20);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [puzzle,X,sX,mX,list,bv]=improver1(puzzle,X,sX,mX,list,iX,iList,n);
bv=0;
if numel(list)
for iRepeat=1:n
bVerb=0;
for i=iList
bV=0;
x=puzzle(i);
bbb=iX==i;
iii=sum(bbb)>0;
Xhoog=sX>mX;
nXh=sum(Xhoog);
nXhi=sum(Xhoog(iii));
dx=list-x;
dS=(nXhi-nXh/9)*dx;
j=find(dS<0);
if numel(j)
dx=dx(j);
Nh=dx;
Nhi=dx;
for k=1:length(j)
Nh(k)=sum(sX>mX+dx(k)/9);
Nhi(k)=sum(sX(iii)>mX-8/9*dx(k));
end
bNOK=(Nhi-Nh/9)*(nXhi-nXh/9)<0;
j(bNOK)=[];
if length(j)
dx(bNOK)=[];
if 1
if dx(1)>0
k=length(j);
else
k=1;
end
else
[mx,k]=max(abs(dx));
end
dx=dx(k);
j=j(k);
bV=1;
end
end
if bV
X(bbb)=list(j);
sX=sum(X);
mX1=mX+dx/9;
puzzle(i)=list(j);
bVerb=1;
list(j)=x;
mX=mX1;
bv=1;
end
end
if ~bVerb;break;end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [puzzle,X,sX,bv]=improver2(puzzle,X,sX,mX,iX,iList,n);
bv=1;
sX0=sum(abs(sX-mX));
for i=1:n
bVerb=0;
for iEl=iList
x=puzzle(iEl);
I=find(sum(iX==iEl));
s=sum(sX(I)-mX);
dx=puzzle(iList)-x;
S0=sX-mX;
if s>0
j=find(dx<0&dx+s*2>0);
else
j=find(dx>0&dx+s*2<0);
end
if numel(j)
T0=sum(abs(S0));
m=0;
for k=1:length(j)
S=S0;
S(I)=S(I)+dx(j(k));
J=find(sum(iX==iList(j(k))));
S(J)=S(J)-dx(j(k));
T=sum(abs(S));
if T<T0
m=k;
T0=T;
end
end
if m
j=iList(j(m));
puzzle(iEl)=puzzle(j);
puzzle(j)=x;
X=puzzle(iX);
sX=sum(X);
bVerb=1;
bv=1;
end
end
end
if ~bVerb;break;end
end
%--------------------------------------------------------
% Rearrange
function result=rearrange(a,ar,flag)
result=[];
R=checkresult(ar);
aar=ar;
if flag==1 %row wise
for i=1:9
pres=R;
[val,index]=setdiff(aar(i,:),a(i,:));
n=numel(val);
shift=ones(1,n);
% Only Circular shift
for j=1:n-1
val=circshift(val,shift);
aar(i,index)=val;
res=checkresult(aar);
if res<0.5*R
result=aar;
break
elseif res<R && res>=0.5*R && res<=pres
result=aar;
pres=res;
elseif res>R
aar=ar;
end
end
%if res>R, aar=ar; end
%Other shifts, swap 1st-last and shift
if n>2
[nval,index]=setdiff(aar(i,:),a(i,:));
temp=nval(1); nval(1)=nval(end); nval(end)=temp; %swap 1st and last
% Circular shift
for j=1:n-1
nval=circshift(nval,shift);
aar(i,index)=nval;
res=checkresult(aar);
if res<0.2*R
result=aar;
return
%break
elseif res<R && res>=0.2*R && res<=pres
result=aar;
pres=res;
elseif res>R
aar=ar;
end
end
end
%if res>R, aar=ar; end
end
else %columnwise
for i=1:9
pres=R;
[val,index]=setdiff(aar(:,i),a(:,i));
n=numel(val);
shift=ones(1,n);
% Only Circular shift
for j=1:n-1
val=circshift(val,shift);
aar(index,i)=val;
res=checkresult(aar);
if res<0.5*R
result=aar;
break
elseif res<R && res>=0.5*R && res<=pres
result=aar;
pres=res;
elseif res>R
aar=ar;
end
end
%if res>R, aar=ar; end
%Other shifts, swap 1st-last and shift
if n>2
[nval,index]=setdiff(aar(:,i),a(:,i));
temp=nval(1); nval(1)=nval(end); nval(end)=temp; %swap 1st and last
% Circular shift
for j=1:n-1
nval=circshift(nval,shift);
aar(index,i)=nval;
res=checkresult(aar);
if res<0.2*R
result=aar;
return
elseif res<R && res>=0.2*R && res<=pres
result=aar;
pres=res;
elseif res>R
aar=ar;
end
end
end
%if res>R, aar=ar; end
end
end
if result
result=result;
return;
else
c=setdiff(aar,ar,'rows');
if c
result=aar;
else
result=ar;
end
end
%--------------------------------------------------------
% Check Result
function result=checkresult(a)
%blocks
block=[1 2 3 10 11 12 19 20 21
4 5 6 13 14 15 22 23 24
7 8 9 16 17 18 25 26 27
28 29 30 37 38 39 46 47 48
31 32 33 40 41 42 49 50 51
34 35 36 43 44 45 52 53 54
55 56 57 64 65 66 73 74 75
58 59 60 67 68 69 76 77 78
61 62 63 70 71 72 79 80 81];
target=sum(a(:))/9;
regSum=sum(a(block'));
colSum=sum(a);
rowSum=sum(a');
result=sum(abs(colSum-target))+sum(abs(rowSum-target))+sum(abs(regSum-target));
|
