function [Output, Historico] = mmqnl(Input)
% MMQNL Non-Linear Least-Squares Method
%
% [Output] = mmqnl(Input)
% [Output, Historico] = mmqnl(Input)
%
% -------------------------------------------------------------------------
% About the 'Input' structure arguments.
% 'Input' is a structure with the following fields:
%
% 'x' - Vector or matrix with the experimental data 'xi'.
% 'y' - Vector with the experimental data 'yi'.
% 'model' - Adopted model (theorical function).
% 'beta0' - Vector with initial coefficient values.
% 's' - Vector with uncertainties of 'yi'.
% 'sx' - Vector with the uncertainties of 'xi'.
% 'cor' - [OPTIONAL] Correlation matrix.
% 'V' - [OPTIONAL] Covariance matrix. This field replace the data of 's'
%       and 'cor' fields.
% 'J' - [OPTIONAL] Functions to compute the derivates needed for the adjust 
%       planning matrix.
% 'MaxItera' - [OPTIONAL] If the model is non-linear, this field declares
%       the maximum number of iteractions of the gauss routine. The default
%       number of maximum iteractions is 10.
% 'USA_GM' - [OPTIONAL] If the model is non-linear, this field declares
%       whether the user wants or not to use the Gauss-Marquardt routine. 
%       The default configuration is USA_GM = 1 (TRUE), this way the mmqnl
%       function always use the Gauss-Marquardt in non-linear functions.
%       If the user doesn't want to use the Gauss-Maquardt routine, then 
%       he must declare this field as USA_GM = 0 (FALSE).
%
% -------------------------------------------------------------------------
% About the 'Output' structure arguments.
% 'Output' is a structure with the following fields:
%
% 'A' - Vector with the adjusted parameters.
% 'sA' - Vector with the uncertainties of the adjusted parameters.
% 'VA' - Covariance matrix of the adjusted parameters.
% 'corA' - Correlation matrix of the adjusted parameters.
% 'Q2' - Chi-squared.
% 'NGL' - Degrees of Freedom.
% 'Q2r' - Reduced Chi-squared (Q2/NGL).
% 'probQ2' - Probability of the chi-squared test (the complement of the 
%            chi-squared cumulative distribution function). 
% 'F' - Yi values calculated throgh the interpolated function.
% 'modelo' - Interpolated function calculated throgh the jacobian matrix
%           (or the given model) and the adjusted parameters.
% 'x' - Vector or matrix with the experimental data 'xi'.
% 'J' - Matrix of the jacobian J(i,j)=dF(xi)/dAj.
% 'y' - Vector with the experimental data 'yi'.
% 's' - Vector with the uncertainties of yi.
% 'Dif' - Vector with the absolute residual.
% 'Res' - Vector with the reduced residual.
%
% -------------------------------------------------------------------------
% About the 'Historico' structure arguments.
% 'Historico' is an optional structure with the following fields:
%
% 'n'
% 'Iter'
% 'Q2'
% 'A_atu'
% 'A'
% 'sA'
% 'Q2antes'
% 'Lambda'
% 'Aceitou'
%
% This structure is useful only to debug the mmqnl function while
%      studying the convergence of a non-linear case.
%
% -------------------------------------------------------------------------
% Prof. Dr. Zwinglio Guimaraes-Filho e Caike Crepaldi 
% Versao 4.3 - Ultima Modificacao: 08/08/2014 13:07 (Caike)
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% O que falta fazer:
% * Acrescentar mais modelos teoricos no "banco de dados".
% -------------------------------------------------------------------------

EH_LINEAR = 0;

TEM_J = 0;
TEM_modelo = 0;
TEM_geraJ = 0;

% Beginning of the verifications ------------------------------------------

if isfield(Input,'y') && ~isempty(Input.y)
    y = Input.y(:);
else
    error('It is necessary to provide the experimental data yi')
end

% -------------------------------------------------------------------------
TEM_INC_EM_X = 0;
if isfield(Input,'Vx') && ~isempty(Input.Vx)
    TEM_INC_EM_X = 1;
    Vx = Input.Vx;
    sx = sqrt( diag(Vx) ); %#ok<*NASGU>
else
    if isfield(Input,'sx') && ~isempty(Input.sx)
        TEM_INC_EM_X = 1;
        sx = Input.sx(:);
    else
        if isfield(Input,'corx') && ~isempty(Input.corx)
            sx = 1;
        else
            sx = 0;
        end
    end
    if length(sx)==1
        sx = repmat( sx, size(y) );
    end
    if isfield(Input,'corx') && ~isempty(Input.corx)
        TEM_INC_EM_X = 1;
        corx = Input.corx;
        Vx = (sx*sx').*corx;
    elseif TEM_INC_EM_X
        Vx = diag( sx.^2 );
    end
end
if (~TEM_INC_EM_X) || (max(max(abs(Vx-diag(diag(Vx)))))<10*eps)
    TEM_CORRELACAO_EM_X = 0;
else
    TEM_CORRELACAO_EM_X = 1;
end

% -------------------------------------------------------------------------
if isfield(Input,'V') && ~isempty(Input.V)
    Vy = Input.V;
    sy = sqrt( diag(Vy) );
else
    if isfield(Input,'s') && ~isempty(Input.s)
        sy = Input.s(:);
    elseif (TEM_INC_EM_X==0) || ( isfield(Input,'cor') && ~isempty(Input.cor) )
        sy = 1;
    else
        sy = 0;
    end
    if length(sy)==1
        sy = repmat( sy, size(y) );
    end
    if isfield(Input,'cor') && ~isempty(Input.cor)
        cory = Input.cor;
        Vy = (sy*sy').*cory;
    else
        Vy = diag( sy.^2 );
    end
end
if max(max(abs(Vy-diag(diag(Vy)))))<10*eps
    TEM_CORRELACAO_EM_Y = 0;
else
    TEM_CORRELACAO_EM_Y = 1;
end

% -------------------------------------------------------------------------

if isfield(Input,'x') && ~isempty(Input.x)
    x = Input.x;
    if length(size(x))==2 && min(size(x,1))==1
        x = x(:);
    end
else
    x = (1:length(y)).';
end

if isfield(Input,'model') && ~isempty(Input.model)
    model = Input.model;
    TEM_modelo = 1;
    if ~isa(model, 'function_handle') && model(1)=='#'
        if length(model)>4 && strcmpi( model(2:4), 'pol' )
            EH_LINEAR = 1;
            if model(5)~='['
                grau = str2double( model(5:end) );
                if grau<0
                    expoentes = 0:grau;
                else
                    expoentes = abs(grau):-1:0;
                end
            else
                expoentes = str2num( model(5:end) ); %#ok<ST2NM>
            end
            modelo = @(A,x)( repmat(x(:),[1 length(A)]).^repmat(expoentes(:).',[length(x) 1])*A(:) );
            NParam = length( expoentes );
            geraJ = @(qual,A,x)( x.^expoentes(qual) );
            TEM_geraJ = 1;
        end
    else
        modelo = model;
    end
end

if isfield(Input,'J') && ~isempty(Input.J)
    J = Input.J;
    EH_LINEAR = 1;
    TEM_J = 1;
    NParam = size(J,2);
    if ~TEM_modelo
        modelo = @(A,x)J*A;
    end
elseif ~TEM_modelo
    error( 'It is necessary to provide either the model or the jacobian (J)' )
end

if isfield(Input,'beta0') && ~isempty(Input.beta0)
    A_atu = Input.beta0(:);
    NParam = length( A_atu );
    if sum( isnan(A_atu) )==size(A_atu)
        EH_LINEAR = 1;        
    end
    A_atu( isnan(A_atu) ) = 1;
else
    EH_LINEAR = 1;
    if ~exist('NParam','var')
        if isfield(Input,'NParam') && ~isempty(Input.NParam)
            NParam = Input.NParam;
        else
            NParam = input( 'How many parameters?' );
        end
    end
    A_atu = zeros( NParam, 1 );
end

if isfield(Input,'MaxIter') && ~isempty(Input.MaxIter)
    MaxIter = Input.MaxIter;
else
    MaxIter = 100; % Default maximum number of iteractions
end

if isfield(Input,'MaxTime') && ~isempty(Input.MaxTime)
    MaxTime = Input.MaxTime;
else
    MaxTime = 30; % Default maximum time for iteractions
end

if isfield( Input,'USA_GM') && ~isempty(Input.USA_GM)
    USA_GM = Input.USA_GM;
else
    USA_GM = 1-EH_LINEAR;
end

if nargout>=2
    Historico.n = 0;
end

% End of the verifications ------------------------------------------------

NDados = length(y);
NGL = NDados - NParam;

% -------------------------------------------------------------------------
% -------------------------------------------------------------------------

GM_L0 = 1e-4;
GM_DL_ok = 0.01;
GM_DL_bad = sqrt(10);
GM_Lmax = 1e3;
GM_Lfinal = 1e-9;
GM_Lmin = 1e-9;

% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
TEM_CORRELACAO = TEM_CORRELACAO_EM_X + TEM_CORRELACAO_EM_Y;

CONVERGIU = 0;
Lambda = GM_L0;
Iter = 0;
TEMPO_INI = cputime;
while( (Iter<MaxIter) && ( (cputime-TEMPO_INI)<MaxTime ) && (CONVERGIU<=0) );
    Iter = Iter + 1;
    if (Iter==1) || (TEM_INC_EM_X)
        if TEM_INC_EM_X
            if (Iter>1) || (mean(A_atu==0)<1)
                dx = 1e-5;
                dFdx = (modelo( A_atu, x+dx )-modelo( A_atu, x ))/dx;
            else
                [ xSort, indSort ] = unique( x );
                ySort = y(indSort);
                dydxSort = diff(ySort)./diff(xSort);
                xdydxSort = (xSort(1:end-1)+xSort(2:end))/2;
                dFdx = interp1( xdydxSort, dydxSort, x, 'linear', 'extrap');
            end
            VxREB = (dFdx*dFdx.').*Vx;
        else
            VxREB = zeros(size(Vy));
        end
        V = VxREB + Vy;
        s = sqrt(diag(V));
        TEM_CORRELACAO = TEM_CORRELACAO_EM_X + TEM_CORRELACAO_EM_Y;
    end
    
    if EH_LINEAR
        Y = y;
    else
        F = modelo( A_atu, x );
        Y = y - F;
    end
    if ~EH_LINEAR || ~TEM_J
        if TEM_geraJ
            for qP = 1 : NParam
                J(:,qP) = geraJ( qP, A_atu, x ); 
            end
        else
            F0 = modelo( A_atu, x );
            for qP = 1 : NParam
                A_mais = A_atu;
                A_mais(qP) = A_mais(qP) + 1e-3*abs(A_mais(qP)) + 10*sqrt(eps);
                J(:,qP) = (modelo(A_mais,x) - F0)/(A_mais(qP)-A_atu(qP));
            end            
        end
        TEM_J = 1;
    end
    if  TEM_CORRELACAO == 0
        Js = J./(s*ones(1,NParam));
        JiVJ = ( Js.'*Js );
    else
        invV = V\eye(NDados);
%         JiVJ = ( J.' * invV * J );
        JiVJ = ( J.' * ( V \ J ) );
    end
    VA = JiVJ \ eye(NParam);

    % Checks if the correction of the parameters will be accepted.
    ACEITOU = 0;
    if USA_GM
        if Iter==1
            if TEM_CORRELACAO == 0
                Q2 = sum( ( (y-F)./s ).^2 );
            else
%                 Q2 = (y-F).' * invV * (y-F);
                Q2 = (y-F).' * ( V \ (y-F) );
            end
        end
        Q2antes = Q2;
    end
    while( ACEITOU==0 )
        ACEITOU=1;
        if USA_GM==0 || Lambda==0
            JiVJef = JiVJ;
            VAef = VA;
        else
            JiVJef = JiVJ + Lambda*diag( diag(JiVJ) );
            VAef = JiVJef \ eye(NParam);
        end
        
        if  TEM_CORRELACAO == 0
%             A = VAef * Js.' * (Y./s);
            A = JiVJef \ ( Js.' * (Y./s) );
        else
%             A = VAef * J.' * invV * Y;
            A = JiVJef \ ( J.' * ( V \ Y ) );
        end
        
        if EH_LINEAR
            F = modelo( A, x );
        else
            F = modelo( A_atu + A, x );
        end
        if  TEM_CORRELACAO == 0
            Q2 = sum( ( (y-F)./s ).^2 );
        else
%             Q2 = (y-F).' * invV * (y-F);
            Q2 = (y-F).' * ( V \ (y-F) );
        end
        if nargout>=2            
            Historico.n = Historico.n + 1;
            n = Historico.n;
            Historico.EH_LINEAR = EH_LINEAR;
            Historico.TEM_J = TEM_J;
            Historico.J{n} = J;            
            Historico.TEM_INC_EM_X = TEM_INC_EM_X;
            Historico.TEM_CORRELACAO = TEM_CORRELACAO;
            Historico.TEM_CORRELACAO_EM_Y = TEM_CORRELACAO_EM_Y;
            Historico.TEM_CORRELACAO_EM_X = TEM_CORRELACAO_EM_X;
            Historico.Iter(n) = Iter;
            Historico.Q2(n) = Q2;
            Historico.A_atu(:,n) = A_atu;
            Historico.A(:,n) = A;
            Historico.sA(:,n) = sqrt( diag(VA) );
        end
        if USA_GM
            if (Q2<=Q2antes) || ((abs(Q2-Q2antes)/(Q2antes+1e-15))<1e-3)
                if Lambda>GM_Lmin
                    Lambda = Lambda*GM_DL_ok;
                end
            else
                if Lambda<GM_Lmax
                    Lambda = Lambda*GM_DL_bad;
                    ACEITOU = 0;
                end
            end
            if nargout>=2
                Historico.Q2antes(n) = Q2antes;
                Historico.Lambda(n) = Lambda;
                Historico.ACEITOU(n) = ACEITOU;
            end
        end
    end
    
    % Checks if the adjustment has already converged.
    if (EH_LINEAR) && (TEM_INC_EM_X==0)
        CONVERGIU = 1;
    else
        if (Iter>1) && ( (USA_GM==0) || (Lambda<=GM_Lfinal) )
            if EH_LINEAR
                dA = A - A_atu;
            else
                dA = A;
            end
            if ( sum(abs(dA)./sqrt(diag(VA)) ) < (1e-2*sqrt(Q2/NGL)) ) || ...
                   ( sum(abs(dA)./sqrt(diag(VA)) ) < 1e-15 )
                CONVERGIU = 1;
            end
        end
        if EH_LINEAR
            A_atu = A;
        else
            A_atu = A_atu + A;
            A = A_atu;
        end
    end
    
end


sA = sqrt( diag(VA) );
corA = VA ./ ( sA *(sA.') );

% Creating the Output structure -------------------------------------------
Output.Input = Input;
Output.A = A;
Output.sA = sA;
Output.VA = VA;
Output.corA = corA;
Output.Q2 = Q2;
Output.NGL = NGL;
Output.Q2r = Q2 / NGL;
Output.probQ2 = 1 - chi2cdf( Q2, NGL );
Output.x = x;
Output.y = y;
Output.s = s;
Output.F = F;
Output.Dif = y - F; % residuo absoluto
Output.Res = Output.Dif ./ s; % residuo reduzido
Output.modelo = modelo;
Output.Iter = Iter;

Output.V = V;