% inversio-ongelmien laskennallinen pk 2012

% Exercise 6

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);

% 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);



% We show here the solution for problem 3
%
%  m3_n = C1*x
%
% Other problems are solved analogiously

% Try first different parameter values
d = [0.001 0.01 0.1 1 10 100 1000 10000];

figure

for i=1:length(d),
    
    A=[C1; sqrt(d(i)) * eye(700)];
    m = [m3_n'; zeros(700,1)];
    
    regsol = A\m;
    
    discr(i) = norm(C1*regsol - m3_n');
    
    d(i)
    plot(s,'r')
    hold on
    plot(regsol)
    hold off
    
    pause
    
end


% Try another interval for which the discrepancies are near the error level e3

d = linspace(0.005,0.01,50);
for i=1:length(d),
    
    A=[C1; sqrt(d(i)) * eye(700)];
    m = [m3_n'; zeros(700,1)];
    
    regsol = A\m;
    
    discr(i) = norm(C1*regsol - m3_n');
    
    %d(i)
    %plot(s,'r')
    %hold on
    %plot(regsol)
    %hold off
    
    %pause
    
end

% Plot discrepancies and the error level

figure
plot(discr,'o')
hold on
plot([1 50],[e3 e3],'r')
hold off

% Plot the solution whose discrepancy  is nearest the error level

d = d(9);

 A=[C1; sqrt(d) * eye(700)];
    m = [m3_n'; zeros(700,1)];
    
    regsol = A\m;
    

figure
plot (s,'r')
hold on
plot(regsol)
hold off



% using the first order difference matrix
% We use diffmat1 fuction here


% Try first different parameter values
L = diffmat1(700);
d = [0.0001 0.001 0.01 0.1 1 10 100];


figure

for i=1:length(d),
    
    A=[C1; sqrt(d(i)) * L];
    m = [m3_n'; zeros(700,1)];
    
    regsol = A\m;
    
    discr(i) = norm(C1*regsol - m3_n');
    
    d(i)
    plot(s,'r')
    hold on
    plot(regsol)
    hold off
    
    pause
    
end

% Refine interval
d = linspace(1,10,50);
for i=1:length(d),
    
    A=[C1; sqrt(d(i)) * L];
    m = [m3_n'; zeros(700,1)];
    
    regsol = A\m;
    
    discr(i) = norm(C1*regsol - m3_n');
    
    %d(i)
    %plot(s,'r')
    %hold on
    %plot(regsol)
    %hold off
    
    %pause
    
end

figure
plot(discr,'o')
hold on
plot([1 50],[e3 e3],'r')
hold off


% Plot the solution whose discrepancy is neartest to the error level

dd=d(6);

A=[C1; sqrt(dd) * L];
m = [m3_n'; zeros(700,1)];
    
regsol = A\m;


figure
plot (s,'r')
hold on
plot(regsol)
hold off





