%% SVisTSP demo
% define pattern

N = 20;

disp('TSP problem size')
disp((factorial(N))^2);
disp('Press a key');
pause

y = (0:1:(N-1))/N;
x = y';

cscheme = autumn;
numcols = size(cscheme,1);

figure(1) 
clf
cor = cos(3*pi*x-1)*sin(pi*y+2);
m = min(min(cor));
M = max(max(cor));
cor = numcols*(cor-m)/(M-m)


I = randperm(N);
J = randperm(N);
pcor = cor(I,J);

subplot(3,1,1)
imshow(cor,cscheme);
title('theoretical desired pattern')

subplot(3,1,2)
imshow(pcor,cscheme);
title('initial, experimental pattern')


R = [20,1:19];
Ri = [2:20,1];
M = pcor;
M1 = M-M(:,R);
M2 = M-M(:,Ri);
M3 = M-M(R,:);
M4 = M-M(Ri,:);

sd = sum(sum(M1.^2))+sum(sum(M2.^2))+sum(sum(M2.^2))+sum(sum(M2.^2));

subplot(3,1,3)
imshow(M,cscheme);
title('best current pattern')
disp('Press a key');
pause

for n=1:100000 
I = randperm(20);
J = randperm(20);
M = M(I,J);
M1 = M-M(:,R);
M2 = M-M(:,Ri);
M3 = M-M(R,:);
M4 = M-M(Ri,:);
sdn = sum(sum(M1.^2))+sum(sum(M2.^2))+sum(sum(M2.^2))+sum(sum(M2.^2));
if sdn < sd 
   subplot(3,1,3)
   imshow(M,cscheme);
   title('best current pattern')
   sd = sdn;
   disp('Press a key');
   pause
   n
   sdn
end
end