% We simulate 1-dimensional deconvolution data that does not involve 
% inverse crime. This is done by computing data using a very fine
% discretization and a perturbed point spread function. That data is then
% sampled at the desired computational grid, and random noise is added.
%
% Samuli Siltanen, February 2011


% Choose number of discretization points
n  = 512;

% Choose two noise levels. We compute three kinds of data: no added noise,
% and noisy data with noise amplitudes given by these sigma1 and sigma2.
sigma1 = 0.01;
sigma2 = 0.03;


%%%%%%% STEP ONE: COMPUTE HIGH RESOLUTION DATA WITH PERTURBED PARAMETERS %%

% Parameter a specifies the width of the PSF, must satisfy 0<a<1/2. 
% Value a=0.04 is used in the matrix model, so here we use a slightly 
% different a.
a = 0.041;

% Choose the points at which to compute the convolved function
Nxx = 2000;
xx  = linspace(0,1,Nxx);

% Create integration points
Nxxp = 500;
xxp  = linspace(-a,a,Nxxp);
Dxxp = xxp(2)-xxp(1);

% Evaluate normalized PSF at integration points
psf = zeros(size(xxp));
ind = abs(xxp)<a;
psf(ind) = DC_PSF(xxp(ind),a);
Ca  = 1/(Dxxp*trapz(psf)); % Normalization constant
psf = Ca*psf;

% Initialize result
Af = zeros(size(xx));

% Compute convolution by integration using trapezoidal rule
for iii = 1:Nxx
   targ = DC_target(xx(iii)-xxp);
   Af(iii) = Dxxp*trapz(psf.*targ);
   % Monitor the run
   if mod(iii,100)==0
       disp([iii Nxx])
   end
end


%%%%%%% STEP TWO: CONSTRUCT MATRIX MODEL AND DATA WITH INVERSE CRIME %%%%%%

% Construct sample points
x  = [0:(n-1)]/n;
Dx = x(2)-x(1);

% Parameter a specifies the width of the PSF, must satisfy 0<a<1/2. 
a = 0.04;

% Construct normalized discrete point spread function
nPSF = ceil(a/Dx); 
xPSF = [-nPSF:nPSF]*Dx;
PSF  = zeros(size(xPSF));
ind   = abs(xPSF)<a;
PSF(ind) = DC_PSF(xPSF(ind),a);
Ca   = 1/(Dx*trapz(PSF));
PSF  = Ca*PSF;

% Construct convolution matrix
A = Dx*DC_convmtx(PSF,n);

% Compute ideal data WITH INVERSE CRIME
f   = DC_target(x);
mIC = A*f(:);



%%%%%%% STEP THREE: CONSTRUCT NOISY DISCRETE DATA WITH NO INVERSE CRIME %%%

% Create sample points
y = [0:(n-1)]/n;
y = y(:);

% Interpolate sample values from precomputed convolution. Now the
% measurement m contains only simulated modelling error but no added
% simulated measurement noise. There is no inverse crime committed.
m = interp1(xx,Af,y,'spline');

% Reset random number generator to make it possible to repeat experiments
% with exactly same noise. If you want to experiment with different noise
% components, comment this line.
reset(RandStream.getDefaultStream);

% Create data with small measurement noise and no inverse crime
noise1 = sigma1*max(abs(m))*randn(n,1);
mn1    = m + noise1;

% Create data with large measurement noise and no inverse crime
noise2 = sigma2*max(abs(m))*randn(n,1);
mn2    = m + noise2;

% Compute the amount of simulated modelling error in m
moderr = max(abs(m-mIC))/max(abs(m));
moderr2 = norm(m-mIC)/norm(m);
disp(['Relative sup norm error is ', num2str(moderr)])
disp(['Relative square norm error is ', num2str(moderr2)])
disp(' ')

% Compute the amount of simulated measurement noise in mn1
relerr = max(abs(m-mn1))/max(abs(m));
relerr2 = norm(m-mn1)/norm(m);
disp(['Relative sup norm error in mn1 is ', num2str(relerr)])
disp(['Relative square norm error in mn1 is ', num2str(relerr2)])
disp(' ')

% Compute the amount of simulated measurement noise in mn1
relerrB = max(abs(m-mn2))/max(abs(m));
relerrB2 = norm(m-mn2)/norm(m);
disp(['Relative sup norm error in mn2 is ', num2str(relerrB)])
disp(['Relative square norm error in mn2 is ', num2str(relerrB2)])
disp(' ')

% Save results to file
save data/DC_nocrimedata xx Af A y n m mn1 mn2 sigma1 sigma2 

% Plot results
DC3_nocrimedata_plot



