function GMRFtest;


n=100;

t=[1:n]';

% Distances
D=sqrt((t*ones(1,n)-ones(n,1)*t').^2);
% Covariance matrix
Sigma=exp(-0.1*D);

figure(1);
clf;
subplot(1,2,1), imagesc(Sigma); axis image; colorbar; title('Covariance matrix');

% Cholesky factor
L=chol(Sigma)';
subplot(1,2,2), imagesc(L); axis image; colorbar; title('Cholesky L of covariance matrix');

x=L*randn(n,1); 
figure(2); 
clf; 
plot(x); title('Random GP sample'); 


% Precision matrix
Q=inv(Sigma);

figure(3);
clf;
subplot(1,2,1), imagesc(Q); axis image; colorbar; title('Precision matrix');

% Cholesky factor
LQ=chol(Q);
subplot(1,2,2), imagesc(LQ'); axis image; colorbar; title('Cholesky L of precision matrix');

xQ=LQ\randn(n,1); 

figure(4);
clf; 
plot(xQ); title('Random GP sample'); 


% Change order of indeces
ii=randsample(n,n);
Sigmanew=Sigma(ii,ii); 
Lnew=chol(Sigmanew)';

Qnew=inv(Sigmanew);
LQnew=chol(Qnew)';

figure(5);
clf;
subplot(1,2,1), imagesc(Sigmanew); axis image; colorbar; title('Covariance for changed order');
subplot(1,2,2), imagesc(Lnew); axis image; colorbar; title('Cholesky for changed order');

figure(6);
clf;
subplot(1,2,1), imagesc(Qnew); axis image; colorbar; title('Precision matrix for changed order');
subplot(1,2,2), imagesc(LQnew); axis image; colorbar; title('Cholesky Q for changed order');



