%%
% this constructs 1D problem with anisotropic diffusion regularisaer
% and draws from prior
%
addpath('../','../Lin1D/');
n = 64;
sigma = 0.1; % std.dev. of added white noise
f = zeros(n,1);
f(12) = -1; f(20:28) = 1.5; f(32:36) = 2; % arbitrary function
[y,K] = linblur(f,0.04);
D1 = lindf(f);
df = D1*f;
kappa = exp(- (df).^2);
DKD = D1'*diag(kappa)*D1;

%%
yn = f + 0.1*randn(size(f));
figure(1);clf; hold on;
plot(f,'k');plot(yn,'r');

plot(kappa(1:n),'g');

%% draw from priors
CLapl = inv(D1'*D1);
CPM = inv(DKD);

L1 = chol(CLapl);
disp(num2str(norm(L1'*L1- CLapl)));
RL1 = L1'*randn(n,1000);
CC1 = RL1*RL1';
figure(3); clf;

L2 = chol(CPM);
disp(num2str(norm(L2'*L2- CPM)));
RL2 = L2'*randn(n,1000);
CC2 = RL2*RL2';
%%
figure(3); 
for j = 1:100
    subplot(1,2,1); plot(RL1(:,j));axis([0 n -9 9]);title('Isotropic Prior');
    subplot(1,2,2); plot(RL2(:,j));axis([0 n -9 9]);title('Anisotropic Prior');
    pause(0.5);
end

%%
figure(2);clf; 
subplot(2,3,1);imagesc((D1'*D1));('prior')
subplot(2,3,2);imagesc(CLapl);title('true covariance');
subplot(2,3,3);imagesc(CC1);title('estimated covariance');
subplot(2,3,4);imagesc((DKD));
subplot(2,3,5);imagesc(CPM);
subplot(2,3,6);imagesc(CC2);title('estimated covariance');