%%%%% Tikhonov.m %%%%%%%
%
%script for testing different codes for Tikhonov regularization

%construct ill conditioned A and b
m = 10;         %we assume m < n
n = 50;
for j=1:50
    [U,~] = qr(randn(m));
    [V,~] = qr(randn(n));
    d = sort(abs([5*randn(m-3,1)+0.1; (1e-9)*randn(3,1)]),'descend');
    D = [diag(d) zeros(m,n-m)];
    A = U*D*V';
    x0 = 5*randn(n,1)+0.1;
    b = A*x0;
    %b = 5*randn(m,1)+0.1;
    
    
    
    %my regularization
%     Q = V';
%     Lambda = D'*D;
%     u = Q'*A'*b;
%     f = @(x) maxcurve(x,Lambda,u);
%     options = optimoptions('fsolve','MaxFunctionEvaluations',500);
%     tau = fsolve(f,1);
%     
%     %my solution
%     AA = A'*A+0*eye(n);
%     RHS = A'*b;
%     xmy = AA\RHS;
    
    %backslash approach
    xback = A\b;
    
    %lsqr approach
    xlsqr = lsqr(A,b);
    
    %lsqr approach
    xlsmr = lsmr(A,b,[],0,0,0,50,[],true);
    
    %curvature with lsqr
    % tt = (1e-30:1e-31:1e-28);
    % normx = zeros(length(tt),1);
    % normout = zeros(length(tt),1);
    % for i=1:length(tt)
    %     Abig = [A; sqrt(tt(i))*eye(n)];
    %     bbig = [b; zeros(n,1)];
    %     xtemp = lsqr(Abig,bbig);
    %     normx(i) = norm(xtemp);
    %     normout(i) = norm(A*xtemp-b)^2;
    % end
    % figure(1)
    % clf
    % hold on
    % plot(tt,normout);
    % plot(tt,normx);
    % legend('Residual','Norm(x)','Location','best');
    % hold off
    
    %curvature with lsmr
    tt = (1e-5:1e-6:1e-3);
    normx = zeros(length(tt),1);
    normout = zeros(length(tt),1);
    for i=1:length(tt)
        xtemp = lsmr(A,b,tt(i));
        normx(i) = norm(xtemp);
        normout(i) = norm(A*xtemp-b)^2;
    end
    plot(normx,normout);
    
    
    fprintf('norm(xback): %g   norm(Axback-b): %g   \n',norm(xback),norm(A*xback-b));
    fprintf('norm(xlsqr): %g   norm(Axlsqr-b): %g   \n',norm(xlsqr),norm(A*xlsqr-b));
    fprintf('norm(xlsmr): %g   norm(Axlsmr-b): %g   \n',norm(xlsmr),norm(A*xlsmr-b));
    
    
end


