% 070331 gabriel@pads.ufrj.br

close all; clear all; rand('state',0);
xmin = ones(4,1);
x = floor(rand(4,1)*4); x0 = x; xc = x;
J = mean((x-xmin).^2); J0 = J; Jc = J;
edges = [0 1 3 5 7 9 11 13 15 17 19 22 25 32]/8;
C = [0 0.25 0.5 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 3.00 3.25 4.00];
N = 200000; BKJ = zeros(N,1); BKa = zeros(N,1); T = 0.5;

for k=1:N,    
    x = xc; a = floor(rand(1)*4)+1; b = ceil(rand(1)*3);
    x(a) = rem(x(a)+b,4);
    J = mean((x-xmin).^2);    
    if rand(1) < exp((Jc-J)/T),        
        xc = x; Jc = J;        
    end;    
    BKJ(k) = Jc; BKa(k) = a;   
end;
S = BKJ(ceil(0.5*N):N); H1 = histc(S,edges)/length(S); figure; subplot(2,1,1); stem(C,H1); grid on;

k=0; Z=0;
for a=0:3,
    for b=0:3,
        for c=0:3,
            for d=0:3,
                k=k+1; x = [a b c d]'; Jg(k) = mean((x-xmin).^2); Z = Z+(exp(-Jg(k)/T)); i(k) = 64*a+16*b+4*c+d+1;
            end;
        end;
    end;
end;
H2 = histc(Jg,edges)/length(Jg); hold on; subplot(2,1,1); stem(C,H2,'r-'); % w = sort(Jg); plot(w,exp(-w/T),'k.');
subplot(2,1,2); plot(C,H1./(H2'*256),'bo-'); hold on;
subplot(2,1,2); plot(C,(1/Z)*exp(-C/T),'ro:'); grid on;

% Invariant probability vector

p=(exp(-Jg/T)/Z);

% Markov Chain Transition Matrices

% M1

M1 = zeros(256,256);
for a=0:3,
    for b=0:3,
        for c=0:3,
            for d=0:3,
                j = 64*a+16*b+4*c+d+1;
                x = [a b c d]'; J = mean((x-xmin).^2);
                for n=1:3,
                    a2 = rem(a+n,4); i = 64*a2+16*b+4*c+d+1;
                    x2 = [a2 b c d]'; J2 = mean((x2-xmin).^2);
                    DeltaJ = J-J2;
                    if DeltaJ>0,
                        M1(i,j) = 1/3;
                    else
                        M1(i,j) = 1/3*exp(DeltaJ/T);
                        M1(j,j) = M1(j,j) + 1/3*(1-exp(DeltaJ/T));
                    end;
                end;                
            end;
        end;
    end;
end;
p1 = M1*p';
figure; imshow(M1);

% M2

M2 = zeros(256,256);
for a=0:3,
    for b=0:3,
        for c=0:3,
            for d=0:3,
                j = 64*a+16*b+4*c+d+1;
                x = [a b c d]'; J = mean((x-xmin).^2);
                for n=1:3,
                    b2 = rem(b+n,4); i = 64*a+16*b2+4*c+d+1;
                    x2 = [a b2 c d]'; J2 = mean((x2-xmin).^2);
                    DeltaJ = J-J2;
                    if DeltaJ>0,
                        M2(i,j) = 1/3;
                    else
                        M2(i,j) = 1/3*exp(DeltaJ/T);
                        M2(j,j) = M2(j,j) + 1/3*(1-exp(DeltaJ/T));
                    end;
                end;                
            end;
        end;
    end;
end;
p2 = M2*p';
figure; imshow(M2);

% M3

M3 = zeros(256,256);
for a=0:3,
    for b=0:3,
        for c=0:3,
            for d=0:3,
                j = 64*a+16*b+4*c+d+1;
                x = [a b c d]'; J = mean((x-xmin).^2);
                for n=1:3,
                    c2 = rem(c+n,4); i = 64*a+16*b+4*c2+d+1;
                    x2 = [a b c2 d]'; J2 = mean((x2-xmin).^2);
                    DeltaJ = J-J2;
                    if DeltaJ>0,
                        M3(i,j) = 1/3;
                    else
                        M3(i,j) = 1/3*exp(DeltaJ/T);
                        M3(j,j) = M3(j,j) + 1/3*(1-exp(DeltaJ/T));
                    end;
                end;                
            end;
        end;
    end;
end;
p3 = M3*p';
figure; imshow(M3);

% M4

M4 = zeros(256,256);
for a=0:3,
    for b=0:3,
        for c=0:3,
            for d=0:3,
                j = 64*a+16*b+4*c+d+1;
                x = [a b c d]'; J = mean((x-xmin).^2);
                for n=1:3,
                    d2 = rem(d+n,4); i = 64*a+16*b+4*c+d2+1;
                    x2 = [a b c d2]'; J2 = mean((x2-xmin).^2);
                    DeltaJ = J-J2;
                    if DeltaJ>0,
                        M4(i,j) = 1/3;
                    else
                        M4(i,j) = 1/3*exp(DeltaJ/T);
                        M4(j,j) = M4(j,j) + 1/3*(1-exp(DeltaJ/T));
                    end;
                end;                
            end;
        end;
    end;
end;
p4 = M4*p';
figure; imshow(M4);

% Note que p1, p2, p3, e p4 sao todos iguais a p.
