function x = SA(y,mu,sigma,x,beta,M,invT)
%% Modify the original code of C.A. Bouman because the original code
%% doesn't have annealing procedure. And the results show that my code is
%% correct while the original one is wrong. After ten iterations, my code
%% generates better output. However, after 20 steps, the two algorithms
%% give almost the same performance.
% Modified by Arthur Tang 2010-03-03

[row_max,col_max] = size(y);

for row = 1:row_max
for col = 1:col_max
 s = [row,col];
 for k = 1:2
  E(k) = energy(k,s,y,mu,sigma,x,beta);
 end

 deltaE = E(2)-E(1);
 prob = 1/(1+exp(-deltaE*invT));  % original code doesn't have this invT, added by TKT

 r = random('Uniform',0,1,1,1);
 if r<prob                       % the following is a Gibbs Sampler. prob is the conditional distribution of Ising Model
   x(row,col) = 1;
 else
   x(row,col) = 2;
 end
end
end

