%% Problem Sheet 5
% Stochastic Simulation

close all; clear all; clc;

%% Exercise 5

N = 10^5; % sample size
n = 5; % dimension
lambda = 0.5; % (0.5, 1.5)

% energy function (for legal configurations)
energy = @(x) -sum(sum(x))*log(lambda);

X = zeros(n,n); % initial state
n_of_ones = zeros(N,1); % vector to save number of ones in each iteration
for i=1:N
    % generate a proposal
    index = ceil(n*rand(1,2));
    Y = X;
    Y(index(1),index(2)) = 1-Y(index(1),index(2));
    % check if it is legal
    legal = true;
    for j=1:n
        for k=1:(n-1)
            if Y(j,k)==1 && Y(j,k+1)==1
                legal = false;
                break;
            end
            if Y(k,j)==1 && Y(k+1,j)==1
                legal = false;
                break;
            end
        end
        if ~legal
            break;
        end
    end
    % if yes, check if we accept
    if legal
        alpha = exp(-energy(Y)+energy(X));
        if rand() <= alpha
            X = Y;
        end
    else
        % otherwise do nothing and stay in X
    end
    % remember number of ones
    n_of_ones(i) = sum(sum(X==1));
end
X

fprintf('The estimated mean number of ones is %f\n', mean(n_of_ones));


%% Exercise 6

p = 0.6;
n = 7;

N = 10^4;

% b) + c)

tau = zeros(N,1);
for i=1:N
    % create empty vectors
    X = [];
    Y = [];
    % set initial values
    X(1) = 1;
    Y(1) = n;
    k=1; % next time we will simulate
    % run the Markov chain up to the stopping time
    while X(k) ~= Y(k)
        Z = rand() <= p;
        if Z == 1
            X(k+1) = min(X(k)+1, n);
            Y(k+1) = min(Y(k)+1, n);
        else
            X(k+1) = max(X(k)-1, 1);
            Y(k+1) = max(Y(k)-1, 1);
        end
        k = k+1;
    end
    tau(i) = k-1;
end

% plot the last realization (for b))
figure(1);
clf
for i=1:k
    line([i-1, i], [X(i), X(i)], 'Color', 'blue');
    line([i-1, i], [Y(i), Y(i)], 'Color', 'red');
end
axis([0, k+1, 0, n+1]);

% estimate the expected coupling time (for c))
fprintf('The estimated mean coupling time is %f\n', mean(tau));