clear all; clc;

c1 = @(X,Y,R) exp(-abs(X-Y)/R);
c2 = @(X,Y,R) (1 + sqrt(3) * abs(X-Y)/R) .* exp(-sqrt(3)*abs(X-Y)/R);
c3 = @(X,Y,R) exp(- (X-Y).^2/R);

R = 1;
R = 0.1;
R = 0.01;
%% (a)

% Diskretisierung
h = 0.001;
x = 0:h:1;
N = length(x);

% Matrix der Diskretisierungspunkt in 2D
e = ones(1,N);
X = x'*e;
Y = e'*x;

% Punktweise Kovarianzmatrizen
C1 = c1(X,Y,R);
C2 = c2(X,Y,R);
C3 = c3(X,Y,R);

% Cholesky-Faktoren der Kovarianzmatrizen
L1 = chol(C1,'lower');
L2 = chol(C2,'lower');
L3 = chol(C3+1e3*eps*eye(N),'lower'); 
% Bemerkung: Korrektur von C3 um Maschinengenauigkeit eps, da sonst 
% numerisch negative Eigenwerte auftreten

% Erzeugen von 3 Realisierungen mittels Cholesky-Faktor
figure();
subplot(3,1,1);
plot(x, L1 * randn(N,3),'LineWidth',2);
title('Realisierungen fuer c_1','FontSize',20); xlabel('x')
set(gca,'FontSize',20)

subplot(3,1,2);
plot(x, L2 * randn(N,3),'LineWidth',2);
title('Realisierungen fuer c_2','FontSize',20); xlabel('x')
set(gca,'FontSize',20)

subplot(3,1,3);
plot(x, L3 * randn(N,3),'LineWidth',2);
title('Realisierungen fuer c_3','FontSize',20); xlabel('x')
set(gca,'FontSize',20)

%% (b)

% KL berechnen
c = @(X,Y) c1(X,Y,R);
[Lambda1, KL1] = ComputeKL1D(c, 1000);

c = @(X,Y) c2(X,Y,R);
[Lambda2, KL2] = ComputeKL1D(c, 1000);

c = @(X,Y) c3(X,Y,R);
[Lambda3, KL3] = ComputeKL1D(c, 1000);

% Plotten der berechneten Eigenwerte
figure();
loglog(1:N, Lambda1,'-or', 1:N, Lambda2,'-+b', ...
    1:N, Lambda3,'-sg','LineWidth',2);
legend('c_1','c_2','c_3')
title('Eigenwerte der drei Kovarianzfunktionen','FontSize',20)

% Abfallraten
r1 = polyfit(log(1:100), log(Lambda1(1:100)'),1)
r2 = polyfit(log(1:100), log(Lambda2(1:100)'),1)

hold on;
loglog(10:100, exp(r1(2)+1 + r1(1)*log(10:100)),'--r',...
    10:100, exp(r2(2)+1 + r2(1)*log(10:100)),'--b');
hold off;
legend('c_1','c_2','c_3',sprintf('Rate: %.2f',r1(1)), sprintf('Rate: %.2f',r2(1)))
set(gca,'FontSize',20)
