%% A bit of a TV on 1D inverse problems

addpath ../RegTools/  % use the subdirectory where you downloaded reg tools

% Total variation regularization for Ax=b amounts to minimizing
% ||Ax-b||_2^2 + \lambda ||Lx||_1
% Whre L is an approximate first derivative matrix and \lambda>0 a regularization factor. Unlike 2-norm (Tikhonov) regularization it allows the solution (considered as a function)  to jump if that fits the data better. The optimization problem is no longer quadratic so requires more sophisticated and costly solvers.
%
% First we set up the standard deriv2 problem but for an x with sharp corners and a jump.


n=100 %Size of problem
example=3;
[A,b,x]=deriv2(n,example);
xj = zeros(n,1); % an x with a step in it
xj((n/2):end)=1;
bj=A*xj;
xk= xj; xk(1:n/4)=(4/n)*((0:(-1+(n/4)))) ;  
xk(n/2:-1:(1+n/4))=(4/n)*((0:(-1+(n/4)))) ;   
plot(xk);title('An x with sharp corners)')

%
% We add 0.1% relative gaussian noise.

bk=A*xk;
bkn = bk + 0.001*norm(bk)*randn(n,1);% Add simulated noise
plot(bk); hold on; plot(bkn); title('Data with and without noise')
hold off
plot(A\bkn);title('Backslash solution.')

TikFac=0.001;
Aa = [A;TikFac*eye(n)]; %Augmented matrix for Tikhonov
bkna=[bkn;zeros(n,1)];

nA=norm(A,2) %Tik fac is relative to the norm of b
cond(A) % he condition number of A

%


TikFac=0.02*nA;
plot([A;TikFac*eye(n)]\bkna );title(['Tikhonov solution ',num2str(TikFac/nA)])

% Try different values of TikFac, by editing code above and run the section
% Now try some Total Variation regularization. See TVreg.m for details of the implementation. It uses standard optimization methods, note there are faster recent methods that are good for 2D and 3D inverse problems.
%
% The values are chosen to illustrate what TV typically does - captures jumps and flat bits and slopes well but flattens off corners.

 xTV = TVreg(A,bk,0.0001);plot(xTV);title('TV with no noise in data')
 xTV = TVreg(A,bkn,0.0001);plot(xTV);title('TV with noise in data and same fac as before')
 xTV = TVreg(A,bkn,0.001);plot(xTV);title('TV, noise fac 10 times larger')
 xTV = TVreg(A,bkn,0.005);plot(xTV);title('TV, noise fac 50 times larger')
 xTV = TVreg(A,bkn,0.05);plot(xTV);title('TV, noise fac 500 times larger')

%
% Now try different values and amounts of noise in the data yourself.
% A natural next step is to combine Tikhonov with TV.