% Selecting the data type

datatype = 'analytic';
%datatype = 'triangle';
%datatype = 'non-analytic';
%datatype = 'inverse crime';

% The Aristotelian prior


n = 150;    %  Length of the vector

% The second order finite difference matrix in interior points
col = zeros(n-2,1); col(1) = 1;
row = zeros(1,n); row(1:3) = [1,-2,1];
L = toeplitz(col,row);

% Augmenting to get the Dirichlet finite difference matrix
row1D = zeros(1,n); row1D(1:2) = [-2,1];
rownD = zeros(1,n); rownD(n-1:n) = [1,-2];
LD = [row1D;L;rownD];

% Calculating the maximal standard deviation at the center of the interval
ej = zeros(n,1); ej(n/2) = 1;
std = norm(LD'\ej);

% Augmenting the interior point matrix with the scaled white noise matrix
% at endpoints to get the Aristotelian difference matrix
row1A = zeros(1,n); row1A(1) = 1/std;
rownA = zeros(1,n); rownA(n) = 1/std;
LA = [row1A;L;rownA];

% The standard deviation of pixels using this LA is of the order std
% We scale it so that the target standard deviation is obtained

std_target = 1;
LA = 100*(std/std_target)*LA;

% Testing the result by  random draws

w = randn(n,5);
x = LA\w;
t = linspace(0,1,n);
figure(1)
plot(t,x,'-','LineWidth',1.5)

% Building the forward matrix

m = 40;
s = exp(log(10)*linspace(-4,0,m));
A = 1/n*exp(-s'*t);



% Calculating the data




if strcmp(datatype,'analytic')
   % Calculate the data analytically. The
   % true signal is a step
   a = 0.3;
   b0 = (1-exp(-a*s'))./s';
   xtrue = (t<a);
elseif strcmp(datatype,'non-analytic')
   % Numerivally computed data. A smooth
   % sigmoid is the true signal
   N = 500; % The dense grid used to calculate the data
   T = linspace(0,1,N);
   Adense = 1/N*exp(-s'*T);
   a = 0.3;
   w = 10;
   xdense = 1 ./(1 + exp(w*(T' - a)));
   xtrue = 1 ./(1 + exp(w*(t' - a)));
   b0 = Adense*xdense;
elseif strcmp(datatype,'inverse crime')
   a = 0.3;
   %xtrue = (t'<a);
   xtrue = 1 - t';
   b0 = A*xtrue;
   % estimating the discretization error
   b_analytic = (1-exp(-a*s'))./s';
   discr_err = b0 - b_analytic;
elseif strcmp(datatype,'triangle')
   b0 = -(1-exp(-s'))./(s'.*s') + 1 ./s';
   xtrue = 1 - t';
end  


% Adding low amplitude noise

noiselev = 1e-7*max(b0);
randn('state',1);
noise = noiselev*randn(m,1);
b = b0 + noise;
if strcmp(datatype,'inverse crime')
   % Plot the additive and the
   % discretixation error
   figure(5)
   loglog(s,abs(noise),'k-','LineWidth',1.2)
   hold on
   loglog(s,abs(discr_err),'r-','LineWidth',1.2)
   legend('Additive','Discretization')
   hold off
end


% CM estimate

x = (A'*A + (noiselev)^2*LA'*LA)\(A'*b);

% Calculating the diagonal of the covariance matrix

poststd = zeros(n,1);
aux = zeros(n,1);
GpostInv = (noiselev)^2* LA'*LA + A'*A;
for j = 1:n
    ej = aux; ej(j)=1;
    poststd(j)=noiselev*sqrt(ej'*(GpostInv\ej));
end

figure(3)
fill([t';t(n:-1:1)'],[x+2*poststd;x(n:-1:1)-2*poststd(n:-1:1)],'g')
hold on
plot(t,x,'r-','LineWidth',1.5)
plot(t,xtrue,'b-','LineWidth',1.5)
legend('2*STD envelope','CM','True')
hold off

% The same computation with LD 

LD = (std/std_target)*LD;
x = (A'*A + (noiselev)^2*LD'*LD)\(A'*b);

% Calculating the diagonal of the covariance matrix

poststd = zeros(n,1);
aux = zeros(n,1);
GpostInv = (noiselev)^2* LD'*LD + A'*A;
for j = 1:n
    ej = aux; ej(j)=1;
    poststd(j)=noiselev*sqrt(ej'*(GpostInv\ej));
end

figure(2)
fill([t';t(n:-1:1)'],[x+2*poststd;x(n:-1:1)-2*poststd(n:-1:1)],'g')
hold on
plot(t,x,'r-','LineWidth',1.5)
plot(t,xtrue,'b-','LineWidth',1.5)
legend('2*STD envelope','CM','True')
hold off

