% inversio-ongelmien laskennallinen pk 2012

% Exercise 5


% Problem 1
u = [zeros(1,100) ones(1,500) zeros(1,100)];

x = linspace(0,1,500);
s = [zeros(1,100) 100*x.*(x-1).*(x-0.2).*(x-0.7) zeros(1,100)];


% Boxcar kernel
k1 = ones(1,100)/100;

% Hann window kernel
k2 = hann(100);



% Note that the first argument of conmatrix must be
% COLUMN vector!
C1 = conmatrix(k1',700);
C2 = conmatrix(k2,700);



% Problem 2

% Noiseless measurements, calculated using conv function
m1 = conv(k1,u);
m2 = conv(k2,u);
m3 = conv(k1,s);
m4 = conv(k2,s);


% Adding noise to the measuremnts
m1_n = m1 + 0.01 * max(u) * randn(size(m1));
m2_n = m2 + 0.01 * max(u) * randn(size(m2));
m3_n = m3 + 0.01 * max(s) * randn(size(m3));
m4_n = m4 + 0.01 * max(s) * randn(size(m4));


% Error norms
e1 = norm(m1_n - m1);
e2 = norm(m2_n - m2);
e3 = norm(m3_n - m3);
e4 = norm(m4_n - m4);




% Problem 3

% Trying to solve the deconvolution problems using pseudoinverse
% fails

s1 = C1\m1_n';
s2 = C2\m2_n';
s3 = C1\m3_n';
s4 = C2\m4_n';

% plotting the results
figure
plot(s1)
hold on
plot(u,'r')

figure
plot(s2)
hold on
plot(u,'r')

figure
plot(s3)
hold on
plot(s,'r')

figure
plot(s4)
hold on
plot(s,'r')




% Problem 4


% We show here the solution to the problem
%
%   m2_n = C2 * x
%
% Other problems are solved in the same way
% Note that results may vary because the noise
% is random!!

% First calculate the svd of the matrix C2

[U,S,V] = svd(C2);

% Plot the singular values to see how they look
figure
plot(diag(S))

% Let us choose some  different k's
kk = [1 10 50 75 100 200];

% Solve and plot different truncated solutions
% Calculate also the discrepancy

discr = zeros(1,6);
figure
plot(u,'r')
hold on
for i=1:6,
    
    tS = S(1:kk(i), 1:kk(i));
    
    kk(i)
    
    tpC = V(:,1:kk(i)) * inv(tS) * U(:,1:kk(i))';
    
    tsol = tpC * m2_n';
    
    plot(tsol)
    
    discr(i) = norm(C2*tsol - m2_n');
    
    pause;
end
    

% The optimal k seems to be between 50 and 75, so let's find
% it!

kk = 50:75;
discr = zeros(1,length(kk));
for i=1:length(kk),
    
    tS = S(1:kk(i), 1:kk(i));
    
    kk(i)
    
    tpC = V(:,1:kk(i)) * inv(tS) * U(:,1:kk(i))';
    
    tsol = tpC * m2_n';
    
    discr(i) = norm(C2*tsol - m2_n');
    

end

% Plot the discrepancies and the error level
figure
plot(kk,discr,'o')
hold on
plot([kk(1) kk(end)],[e2 e2],'r')


% Look for the optimal k from the plot or from dicr vector
% and plot the optimal solution, say k=61

tS = S(1:61, 1:61);
tpC = V(:,1:61) * inv(tS) * U(:,1:61)';
tsol = tpC * m2_n';

figure
plot(u,'r')
hold on
plot(tsol)



% Problem 5

% 200 first singular values are probably enough
kk = 1:200;

% Vectors for discrepancy norm and solution norm
dv = zeros(1,200);
sn = zeros(1,200);

for i=1:200,
    tS = S(1:kk(i), 1:kk(i));   
    tpC = V(:,1:kk(i)) * inv(tS) * U(:,1:kk(i))';
    tsol = tpC * m2_n';
    
    % Discrepancy norm
    dn(i) = norm(C2*tsol - m2_n');
    sn(i) = norm(tsol);
end

% Plot L-curve
figure
loglog(sn,dn,'o-')

% Try to find optimal k by looking the plot
% You maybe have to zoom in and plot only some of the points



