clear all
close all
%% Define input images
L=10; % L by L square image
X1=round(rand(L));%left eye - a random mix of bits
X2=round(rand(L));%right eye - a random mix of bits

%% Define priors/hyperparameters
b=ones(L)-1/2;% prior on S values when S in the range ...
tau=1/250 % kT=250
alpha=0 % extent to which zero values of pi are punished. Setting this to zero makes it as likely to be occluded as unoccluded. I set it to zero to make the switching behaviour more frequent.
beta=200 % degree of smoothness permitted in perceptions. tau times this is an inverse temperature.
gamma=200 % degree of smoothness permitted in occlusions. tau times this is an inverse temperature.
sigsqr1=1/200 % "trustworthiness" of right eye data (also a temperature modulator)
sigsqr2=1/200 % "trustworthiness" of left eye data (also a temperature modulator)
prior=0; % switches on/off prior terms. Since this complicates the dynamics (and is irrelvant when the images Xi are binary and b_n=1/2) we switch it off.

%% initializing the variables
W=ones(L)/2; % intitialization for burn-in. W is the left-right bias
S=X1.*W+(1-W).*X2; % initialization of S - S is derived from W but will be used for brevity
PI1=round(rand(L)); %intitialization for burn-in. Left eye occlusions initialized to zero
PI2=round(rand(L)); %intitialization for burn-in. Right eye occlusions initialized to zero

%% Sampling is asynchronous across all the variables.

%% Sampling parameters
burnin=200000 %L*L*3*500  % each site will be updated approximately 500 times. I call the variable burnin but it contains both genuine sampling and burn-in - we don't distinguish. How could we?

%% Variables we will use to log the simulation as it unfolds
avw=0.5 % this variable records the average support for model 1 at each time step
avp1=0.5 % this variable records the average level of left eye occlusion at each time step
avp2=0.5 % this variable records the average level of right eye occlusion at each time step
for i=1:burnin
    
    loc= randi(L^2,1); %picks a random pixel in image using linear indexing. Look up linear indexing for matlab if you're curious - it just avoids having to track to indices.
    neigbours=toroid(loc,L);% this function returns the 4 neighbouring linear indices of the linear index loc in toroidal co-ordinates
    p=rand; % decides if we have an update to perception or occlusion
    if p>1/3 % in this case we update the occlusion PI?
        q=rand; % decides which eye to update
        if q>0.5 % update left eye occlusion PI1
            
            cup=exp(-tau*(alpha*0+sum(gamma*(PI1(neigbours)-1).^2))); % A modified version of Hamiltonian in lectures. This fixes some bugs in that version.
            cdwn=exp(-tau*(alpha*1+sum(gamma*(PI1(neigbours)-0).^2)));% if the node has state "0" % in the original paper the term alpha*pi is a typo - if it punishes outliers it should be alpha*(1-pi)
            
            ctot=cup+cdwn; % finds normalization (easy! that's why Gibbs Sampling works - the marginals are sometimes easy to normalize)
            
            if rand>cup/ctot % assigns spin up and spin down appropriately
                PI1(loc)=0; % assigns new - spin down state to this point
            else
                PI1(loc)=1; % assigns new - spin up state to this point
            end
            
        else %update right eye occlusion PI2 - slightly more compact version of code above for PI1
            
            cup=exp(-tau*(alpha*0+sum(gamma*(PI2(neigbours)-1).^2))); % if the node has state "1"
            cdwn=exp(-tau*(alpha*1+sum(gamma*(PI2(neigbours)-0).^2))); % if the node has state "0"
            
            ctot=cup+cdwn; % finds normalization (easy! that's why Gibbs Sampling works - the marginals are sometimes easy to normalize)
            
            if rand>cup/ctot % assigns spin up and spin down appropriately
                PI2(loc)=0; % assigns new - spin down state to this point
            else
                PI2(loc)=1; % assigns new - spin up state to this point
            end
        end
        
    else % in this case we update the perception W (or S)
        %% in this version the perception is binary and forced to either favour the left eye or the right on a pixel-by-pixel basis.  
        wup=exp(-tau*( prior*(b(loc)-X1(loc))^2 + sum(beta*(W(neigbours)-1).^2)))*exp(-tau*((PI1(loc)*(X1(loc)-X1(loc))^2)/(2*sigsqr1) + (PI2(loc)*(X2(loc)-X1(loc))^2)/(2*sigsqr2) )); % left eye bias
        wdwn=exp(-tau*( prior*(b(loc)-X2(loc))^2 + sum(beta*(W(neigbours)-0).^2)))*exp(-tau*((PI1(loc)*(X1(loc)-X2(loc))^2)/(2*sigsqr1) + (PI2(loc)*(X2(loc)-X2(loc))^2)/(2*sigsqr2) )); % right eye bias
        
        wtot=wup+wdwn; % finds normalization
        
        if rand>wup/wtot % assigns spin up and spin down appropriately
            W(loc)=0; % assigns new - spin down state to this point
        else
            W(loc)=1; % assigns new - spin down state to this point
        end
        
        
    end
    avw=[avw sum(sum(W))]; % incrementing this order parameter of support for left eye image
    avp1=[avp1 sum(sum(PI1))]; % incrementing this order parameter for left eye occlusion
    avp2=[avp2 sum(sum(PI2))]; % incrementing this order parameter for right eye occlusion
end % ends burn-in loop

subplot(3,1,1); plot(avw)
title('%perceived support for left eye')
xlabel('Time')
ylabel('%perceived support for left eye')
subplot(3,1,2); plot(avp1)
title('%of left eye unoccluded')
xlabel('Time')
ylabel('%of left eye unoccluded')
subplot(3,1,3); plot(avp2)
title('%of right eye unoccluded')
xlabel('Time')
ylabel('%of right eye unoccluded')

