%% ising1.m
% from https://engineering.purdue.edu/~bouman/publications/tutorials/mrf_tutorial/matlab/Matlab/Ising/Sim.m

beta = 1.0; % value of transition cost
N = 2^7; % size of image = (NxN) 
M = 100; % number of iterations = M
M0 = 10; % number of displayed frames = M0
samp_rate = ceil(M/M0);

%color = [0,0,0; 1,1,1];
%colormap(color);
%M_data = moviein(M0);

%% show some random bernoulli images; the last also serves as initial condition
if ~isvar('xb')
	rand('seed', sum(100*clock));

	nbern = 20;
	xb = zeros(N, N, nbern, 'single');
	for ii=1:nbern
		x = single(rand(N) > 0.5);
		xb(:,:,ii) = x;
		im(x)
		title 'initial random distribution'
		drawnow
		pause(0.1)
	%	prompt
	end

%	ir_savefig ising1a
%	movie2(xb, 'file', 'ising1a.avi', 'fps', 10)
prompt
end

%M_data(:,1) = getframe;
%print ising1 -deps

%% metropolis sampler
if ~isvar('xi')
	xi = zeros(N,N,M, 'single');
	xi(:,:,1) = x;
	for n=2:M
	for i=1:N
	for j=1:N
		ip = rem(N+i,N)+1;
		im0 = rem(N+i-2,N)+1;
		jp = rem(N+j,N)+1;
		jm = rem(N+j-2,N)+1;

		v = (1~=x(ip,j)) + (1~=x(i,jp)) + (1~=x(im0,j)) +(1~=x(i,jm));

		prob = 1 / ( 1+exp(2*beta*(v-2)) );
		x(i,j) = single(rand(1) < prob);
	end
	end
	if 1 %rem(n,samp_rate)==0
		xi(:,:,n) = x;
		im(x)
		titlef('Ising model: beta = %g, iteration = %d', beta, n)
drawnow
		pause(0.1)
%	tmp =
%	ir_savefig
%	M_data(:,n) = getframe;
	end
	end

%	movie2(xi, 'file', 'ising1b.avi', 'fps', 10)
%	ir_savefig ising1b
end

%movie(M_data,10,6)

if 0
	clf, im('notick', x(:,1:end/2), 'Ising process')
%	ir_savefig ising1c
end
