%% Problem Sheet 10
% Stochastic Simulation

close all; clear all; clc;

%% Exercise 5 a)

lambda = 5;
t_max = 10;

% 1.
S = [];
S_sum = 0;
while S_sum < t_max
    % simulate inter-arrival times
    S_i = -log(rand())/lambda;
    S = [S S_i];
    S_sum = S_sum + S_i;
end
% compute jump times
T = cumsum(S);
% for plot: add zero and truncate to t_max
T = [0 T];
T(end) = t_max;
N_t = 0:length(S)-1;

figure(1);
clf
for i = 1:length(T)-1
    line([T(i), T(i+1)], [N_t(i), N_t(i)]);
end
    
% 2.
% make a grid for the time
h = 0.05;
t = 0:h:t_max;

% simulate jump times (approximately)
T = (rand(1, length(t)-1) <= lambda*h);
% compute values of N_t on the grid
N_t = [0 cumsum(T)];

figure(2);
clf
for i = 0:max(N_t)
    line([min(t(N_t==i)), max(t(N_t==i))], [i i]);
end


%% Exercise 5 b)

N = 10^4;
N_1_1 = zeros(N, 1);
N_1_2 = zeros(N, 1);

% make a grid for the time (only once)
h = 0.05;
t = 0:h:t_max;

for i = 1:N
    
    % 1.
    S = [];
    S_sum = 0;
    while S_sum < t_max
        % simulate inter-arrival times
        S_i = -log(rand())/lambda;
        S = [S S_i];
        S_sum = S_sum + S_i;
    end
    % compute jump times
    T = cumsum(S);
    N_1_1(i) = sum(T <= 1);
    
    % 2.
    % simulate jump times (approximately)
    J = (rand(1, length(t)-1) <= lambda*h);
    % compute values of N_t on the grid
    N_t_2 = [0 cumsum(J)];
    N_1_2(i) = N_t_2(1/h + 1);
end

figure(1);
h = hist(N_1_1, 0:20);
bar(0:20, h/sum(h), 1);

figure(2);
h = hist(N_1_2, 0:20);
bar(0:20, h/sum(h), 1);

figure(3);
bar(0:20, poisspdf(0:20, lambda), 1);


%% Exercise 6

clear all; close all;

f = @(x) 0.5*x.^2-10*exp(-(x.^2)/20).*cos(3*x);
x = linspace(-10, 10, 500);

% figure(1);
% plot(x, f(x));

X_0 = 10;
T_0 = 10;
N = 10^3;
sigma = 0.1; % [0.1, 0.5, 1]
beta = 0.99;

X = zeros(N+1,1);
X(1) = X_0;
T = T_0;

for i=1:N
    Y = X(i) + sigma*randn();
    alpha = exp( (f(X(i))-f(Y)) / T );
    if rand() <= alpha
        X(i+1) = Y;
    else
        X(i+1) = X(i);
    end
    T = beta*T;
end
X(end)

figure(2);
[ax, p1, p2] = plotyy(x, f(x), X, 0:N);
title('Simulated Annealing')
xlabel('x');
ylabel(ax(1), 'f(x)');
ylabel(ax(2), 'iteration');
axis(ax(1), [-10 10 -20 60]);
axis(ax(2), [-10 10 -100 N+300])
axis(ax(2), 'ij');