function [t,f,s,g_all,ffit_all,g,ffit,lambda,c] = InvLap(t,f,Tmin,Tmax,Tn,method,iter,dim,varargin)

% [s,g,ffit] = InvLap(x,y,Tmin,Tmax,Tn,'method',iter,dim,OPTIONS)
%
% examples: InvLap(t,y,0.1,10000,100,'Press',100,'full');
%           v
% Input variables:
%    x:        vector array:    time
%    y:        vector array:    amplitude  
%    Tmin:     scalar:          lowest time constant in ms
%    Tmax:     scalar:          highest time constant in ms
%    Tn:       scalar:          number of time constants
%    method    string/scalar:   'SDCS' (1)(ref 1) or 'BRD' (2) (ref 4) 'SIMPLEX'   
%    iter:     scalar:          amount of iterations
%    dim:      string/scalar:   dimension of data vectors: 'full' or scalar 
%
% OPTIONS    
%
%   lambda:   string/scalar:        regularization parameter
%        scalar
%       'trace': weight of data and regularisation is comparable  (ref 1)
%       'lcur':  lambda is determined by L-curve method (ref 2)
%       'auto':  lambda is determined by S-curve method (ref 3)
%           'lambda_min': scalar
%           'lambda_max': scalar
%           'lambda_n':   scalar
%
%   plotmode:  scalar
%       0:      no plots  (default)
%       1:      intermediat plots
%
%   TOL: scalar:  (0<TOL<1) Tolerance parameter to determine lambda. Default = 0.005) (ref 3)
%
% Output variables
%   s:       vector array:   time constants
%   g:       vector array:   weigth factors of time constants
%   ffit:    vector array:   resulting decay
%
%
% Example: 
%
% References: 
% 1: W. H. Press, S. A. Teukolsky, W. T. Vetterling & B. P. Flannery (2002). Numerical Recipes
%    in C - The Art of Scientific Computing. Cambridge University Press, Cambridge, second
%    edition edition.
%
% 2: P. C. Hansen (1994). Regularization tools: A matlab package for analysis and solution
%    of discrete ill-posed problems. Numerical Algorithms, 6(1):1 – 35.
%
% 3: E. Fordham, S. A. & L. Hall (1995). Imaging multiexponential relaxation
%    in the (y, loget1) plane, with application to clay filtration in rock
%    cores. Journal of Magnetic Resonance, Series A, 113(2):139 – 150.
%
% 4: L. Venkataramanan, Y.-Q. Song & M. D. Hurlimann (2002). Solving fredhol integrals of the
%    firts kind with tensor product structure in 2 and 2.5 dimensions. IEEE Transactions on
%    signal processing, 50(5):1017 – 1026.

%%  set output grid s
s = logspace(log10(Tmin),log10(Tmax),Tn);

%% reducing dimensions if wanted

if strcmp(dim,'full')==1
    dim=length(t);
else
    if dim>length(t)
    error('info:dim','Dimension exceeds length data');
    end
    t=reddim(t,dim);
    f=reddim(f,dim);
end
 
%%  calculate K(t,s)=exp(-t/s): matrix a 
t=t(:);
s=s(:)';
A=exp(-t*(1./s));

%% calculate reg matrix H
m=length(s);
B=eye(m); %zeroth order regularisation
H=B'*B; 

%% Default values for options
c=1;
plotmode=0;
TOL=0.005;
lambda=trace(A'*A)/trace(H);
lambda_min=-3;
lambda_max=6;
lambda_n=10;

    
%% checking options
if nargin>8
for i=1:2:nargin-8
    pname = varargin{i};
    pval =  varargin{i+1};
    switch lower(pname)
        case('lambda')
            if ischar(pval)==0;
                lambda=pval;
            else           
                switch pval
                    case('trace')
                        lambda=trace(A'*A)/trace(H)
                    case('lcur')
                        [U1,s1,V1]=csvd(A);
                        [reg_corner,rho,eta,reg_param] = l_curve(U1,s1,f,'Tikh',B,V1);
                        lambda=reg_corner;
                    case('auto')
                        lambda=logspace(lambda_min,lambda_max,lambda_n);
                    otherwise
                        lambda;
                end
            end
             
        case('plotmode')
            switch pval
                case('on')
                    plotmode=1;
                case('off')
                    plotmode=0;
                otherwise
                    error('info:UnknownParameterName',...
                'Unknown plotmode: %s.',pname);
            end
        
        case('tol')
            TOL=pval;
            
        case('lambda_min')
            lambda_min=pval;
            lambda=logspace(lambda_min,lambda_max,lambda_n);
        case('lambda_max')
            lambda_max=pval;
            lambda=logspace(lambda_min,lambda_max,lambda_n);
        case('lambda_n')
            lambda_n=pval;
            lambda=logspace(lambda_min,lambda_max,lambda_n);
        
        otherwise
            error('info:UnknownParameterName',...
                'Unknown option: %s.',pname);
    end
end
end


chi=zeros(1,length(lambda));
g_all=zeros(length(s),length(lambda));
ffit_all=zeros(length(t),length(lambda));

for j=1:length(lambda)
    lambda_j=lambda(j);

switch method
    case('SDCS')    
        a=A'*A+lambda_j*H;
        y=A'*f;
        [L,U]=lu(a);
        g=U\(L\y);
        e=2/max(eig(a));

        a=A'*A;  %no regularisation now

        h=waitbar(0,['Calculation distribution with \alpha = ',num2str(lambda_j)]);
        if plotmode==1
            h2=figure;
        end

        for i=1:iter;
        waitbar((i+iter*(j-1))/(iter*length(lambda)));  
    
            g=(g>0).*g;   %map neg to zero
            g(1)=0;
            g(end)=0;
            g=(eye(m)-e*a)*g+e*y;   % Eq. 18.5.21   (b=0)
    
            if plotmode==1 && rem(i,iter/50)==0  %plotting
               title_string=['Alpha = ',num2str(lambda_j)];
               InvLap_plot(h2,title_string,s,g,t,f,A);
            end
        end
        close(h)
        g=(g>0).*g;    %map neg to zero again
        g(1)=0;
        g(end)=0;
        ffit=A*g;      
        
        
        case('Simplex')    
        g0=lognpdf(s,log(100),0.2);
        g0=g0(:);
        [g,tfit,ffit,eltime] = simplex(t,f,g0,s,lambda_j,1,iter,'test',@semilogx,'decay');
            
            if plotmode==1 && rem(i,iter/50)==0  %plotting
                   title_string=['Alpha = ',num2str(lambda_j)];
                   InvLap_plot(h2,title_string,s,g,t,f,A);
            end

            
    case('BRD')
        g0=lognpdf(s,log(100),0.2);
        g0=g0(:);
        c=(A*g0-f)/-lambda_j;
        h=waitbar(0,['calculation distribution with \alpha = ',num2str(lambda_j)]);

        if plotmode==1
           h2=figure;
        end

        for i=1:iter
            waitbar((i+iter*(j-1))/(iter*length(lambda)));  
            G=A*  diag((A'*c>0))  *A';
            Dchi=(G     +   lambda_j*eye(length(f)))*c-f;
            DDchi=G     +   lambda_j*eye(length(f));
            c=c-(DDchi)\Dchi;
            g=A'*c;
            g=(g>0).*g;
            ffit=A*g;
            if plotmode==1 %plotting
               title_string=['Alpha = ',num2str(lambda_j)];
               InvLap_plot(h2,title_string,s,g,t,f,A);
            end
        end
        close(h)
end
g_all(:,j)=g;
ffit_all(:,j)=ffit;
chi(j)=sum((A*g-f).^2)./dim;
end
s=s';

if j>1
figure
loglog(lambda,chi,'bx-')

[r c]=find((TOL-diff(log10(chi))./diff(log(lambda))).^2==min((TOL-diff(log10(chi))./diff(log(lambda))).^2));
lambda_j=lambda(c(end));
hold on
loglog(lambda_j,chi(c(end)),'o','Markerfacecolor','b');
xlabel('\alpha')
ylabel('\chi^2')
g=g_all(:,c);
ffit=ffit_all(:,c);
end

if j>1 || plotmode==0
h3=figure;
title_string=['Alpha = ',num2str(lambda_j)];
InvLap_plot(h3,title_string,s,g,t,f,A);
end

end

%% plotting function
function InvLap_plot(h,title_string,s,g,t,f,A)
h;
subplot(3,1,1)
semilogx(s,g,'xr','LineWidth',2)
title(title_string);
ylabel('Amplitude');xlabel('T_2 (ms)')
hold off

subplot(3,1,2)
plot(t,f,'.');hold on;plot(t,A*g,'r')
ylabel('Amplitude');xlabel('t (ms)')
hold off
            
subplot(3,1,3) 
plot(t,f-A*g,'k');hold on; line([0 max(t)],[0 0],'color','k')
ylabel('Residues');xlabel('t (ms)')
hold off    
        
drawnow
end

%% reducing dimensions
function Ar=reddim(A,l)
[r c]=size(A);

% reducing length must be dimension rows
if r<c
    A=A';
end
%reducing factor
rf=round(size(A,1)./l);
%new size A
Ar=zeros(round(size(A,1)/rf-1),size(A,2));
%reducing
for i1=1:round(size(A,1)/rf-1)
    for i2=1:rf
        Ar(i1,:)=Ar(i1,:)+A((i1-1)*rf+i2,:);
    end
end
Ar=Ar./rf;
if r<c
    Ar=Ar';
end
end