function  [THIS,SUCCESS] = mysstatenonlin(THIS,USROPT)
% MYSSTATENONLIN [Not a public function] Steady-state solver for non-linear models.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

persistent ...
    NAMEBLKLEVEL NAMEBLKGROWTH EQTNBLK NBLK BLOCKFUNC OPT ISREFRESH LINK;

% Skip the following user option preprocessor if this is the second or
% further persistent call to `mysstatenonlin`.
if ~USROPT.persist || isempty(OPT)

    dooptions2fix();
    
    if ~isempty(USROPT.algorithm)
        USROPT.solver = USROPT.algorithm;
        utils.warning('obsolete', ...
            ['ALGORITHM is an obsolete option and will be removed in the future. ', ...
            'Use SOLVER instead.']);
    end

    if islogical(USROPT.display)
        if USROPT.display
            USROPT.display = 'iter';
        else
            USROPT.display = 'off';
        end
    end
    
    USROPT.optimset = optimset( ...
        'display',USROPT.display, ...
        'maxiter',USROPT.maxiter, ...
        'maxfunevals',USROPT.maxfunevals,  ...
        'tolx',USROPT.tolx, ...
        'tolfun',USROPT.tolfun, ...
        'algorithm','levenberg-marquardt', ...
        USROPT.optimset{:});
    % We must use Levenberg-Marquardt because it can handle underdetermined systems.
    
    OPT = USROPT;
    
end

%**************************************************************************

nalt = size(THIS.Assign,3);
SUCCESS = true(1,nalt);

% Set the level and growth of optimal policy multipliers to zero. We must
% do this before checking for NaNs in fixed variables.
if OPT.zeromultipliers
    THIS.Assign(1,THIS.multiplier,:) = 0;
end

% Check for levels and growth rate fixed to NaNs.
dochkfornans();

if ~OPT.persist ...
        || isempty(NAMEBLKLEVEL) || isempty(NAMEBLKGROWTH) ...
        || isempty(EQTNBLK) || isempty(NBLK) ...
        || isempty(BLOCKFUNC)
    
    % Initialise for each non-persistent call, or whenever the persistent
    % variables are empty.
    
    doinitialise();    
end

for ialt = 1 : nalt
    % Assign initial values.
    x = real(THIS.Assign(1,:,ialt));
    % Make sure steady state of shocks is always zero.
    x(THIS.nametype == 3) = 0;
    if ~isempty(OPT.resetinit)
        x(:) = real(OPT.resetinit);
    else
        % Assign NaN initial conditions.
        index = isnan(x) & THIS.nametype <= 2;
        x(index) = real(OPT.naninit);
        x(THIS.log) = log(x(THIS.log));
    end
    
    % Assign growth rates.
    dx = imag(THIS.Assign(1,:,ialt));
    % Makie sure steady-state growth of shocks is always zero.
    dx(THIS.nametype == 3) = 0;
    if ~isempty(OPT.resetinit)
        dx(:) = imag(OPT.resetinit);
    else
        % Assign NaN initial conditions.
        index = isnan(dx) & THIS.nametype <= 2;
        dx(index) = imag(OPT.naninit);
    end
    % Re-assign zero growth for log-variables to 1.
    dx(dx == 0 & THIS.log) = 1;
    dx(THIS.log) = log(dx(THIS.log));
        
    % Cycle over individual blocks.
    for iblk = 1 : NBLK
        if isempty(NAMEBLKLEVEL{iblk}) && isempty(NAMEBLKGROWTH{iblk})
            continue
        end
        xi = NAMEBLKLEVEL{iblk};
        dxi = NAMEBLKGROWTH{iblk};
        z0 = [x(xi),dx(dxi)];
        
        if OPT.warning
            % Call sstateeval_ to test for NaN and Inf.
            check = BLOCKFUNC{iblk}(x,dx);
            index = isnan(check) | isinf(check);
            if any(index)
                utils.warning('model', ...
                    'This equation evaluates to NaN or Inf: ''%s''.', ...
                    THIS.eqtn{EQTNBLK{iblk}(index)});
            end
        end
        
        % Number of levels; this variables is used also within `doobjfunc`.
        nxi = length(xi);
        
        % Function handles to equations in this block.
        f = BLOCKFUNC{iblk};
        
        % Call the solver.
        switch lower(char(OPT.solver))
            case 'lsqnonlin'
                [z,resnorm,residual,exitflag] = ...
                    lsqnonlin(@doobjfunc,z0,[],[],OPT.optimset); %#ok<ASGLU>
                if exitflag == -3
                    exitflag = 1;
                end
            case 'fsolve'
                [z,fval,exitflag] = ...
                    fsolve(@doobjfunc,z0,OPT.optimset); %#ok<ASGLU>
                if exitflag == -3
                    exitflag = 1;
                end
        end
        
        z(abs(z) <= OPT.optimset.TolX) = 0;
        x(xi) = z(1:nxi);
        dx(dxi) = z(nxi+1:end);
        thissuccess = ~any(isnan(z)) && double(exitflag) > 0;
        SUCCESS(ialt) = SUCCESS(ialt) && thissuccess;
    end

    % TODO: More thorough report on which equations and which variables failed.
    if OPT.warning && ~SUCCESS(ialt)
        utils.warning('model', ...
            'Steady state inaccurate or not returned for some variables.');
    end
    
    x(THIS.log) = exp(x(THIS.log));
    dx(THIS.log) = exp(dx(THIS.log));
    THIS.Assign(1,THIS.nametype <= 2,ialt) = ...
        x(THIS.nametype <= 2) + ...
        1i*dx(THIS.nametype <= 2);
    
end

if ~OPT.persist
    % Clean up after each non-persistent call.
    doreset();
end

% Nested functions.

%**********************************************************************
    function y = doobjfunc(p)
        % This is the objective function for the solver. Evaluate the
        % equations twice, at time t and t+10.
        x(xi) = p(1:nxi);
        dx(dxi) = p(nxi+1:end);
        % Refresh all dynamic links.
        if ISREFRESH
            dorefresh();
        end
        if any(dxi)
            % Evaluate the model equations at time t and t+10 if at least one growth
            % rate is needed.
            x1 = x + 10*dx;
            y = [f(x,dx);f(x1,dx)];
        else
            % Evaluate the model equations at time t.
            y = f(x,dx);
        end
        
        function dorefresh()
            x(THIS.log) = exp(x(THIS.log));
            dx(THIS.log) = exp(dx(THIS.log));
            THIS.Assign(1,:,ialt) = x + 1i*dx;
            THIS = refresh(THIS,ialt);
            x = real(THIS.Assign(1,:,ialt));
            dx = imag(THIS.Assign(1,:,ialt));
            dx(dx == 0 & THIS.log) = 1;
            x(THIS.log) = log(x(THIS.log));
            dx(THIS.log) = log(dx(THIS.log));
        end
        % dorefresh().
        
    end
% doobjfunc(),

%**********************************************************************
    function doreset()
        % Reset all persistent variables.
        NAMEBLKLEVEL = [];
        NAMEBLKGROWTH = [];
        EQTNBLK = [];
        NBLK = [];
        BLOCKFUNC = [];
        OPT = [];
    end
% doreset().

%**********************************************************************
    function dooptions2fix()
        % Process the fix, fixallbut, fixlevel, fixlevelallbut, fixgrowth,
        % and fixgrowthallbut options. All the user-supply information is
        % combined into fixlevel and fixgrowth.
        list = {'fix','fixlevel','fixgrowth'};        
        for i = 1 : length(list)
            fix = list{i};
            fixallbut = [fix,'allbut'];
            
            % Convert charlist to cellstr.
            if ischar(USROPT.(fix)) ...
                    && ~isempty(USROPT.(fix))
                USROPT.(fix) = regexp(USROPT.(fix),'\w+','match');
            end
            if ischar(USROPT.(fixallbut)) ...
                    && ~isempty(USROPT.(fixallbut))
                USROPT.(fixallbut) = ...
                    regexp(USROPT.(fixallbut),'\w+','match');
            end
            
            % Convert fixallbut to fix.
            if ~isempty(USROPT.(fixallbut))
                USROPT.(fix) = setdiff( ...
                    THIS.name(THIS.nametype <= 2 | THIS.nametype == 4), ...
                    USROPT.(fixallbut));
            end
            
            if ~isempty(USROPT.(fix))
                fixpos = mynameposition(THIS,USROPT.(fix));
                validate = ~isnan(fixpos);
                if all(validate)
                    validate = THIS.nametype(fixpos) <= 2 ...
                        | THIS.nametype(fixpos) == 4;
                end
                if any(~validate)
                    utils.error('model', ...
                        'Cannot fix this name: ''%s''.', ...
                        USROPT.(fix){~validate});
                end
                USROPT.(fix) = fixpos;
            else
                USROPT.(fix) = [];
            end
        end
        
        % Add the positions of optimal policy multipliers to the list of fied
        % variables. The level and growth of multipliers will be set to zero in the
        % main loop.
        if USROPT.zeromultipliers
            USROPT.fix = union(USROPT.fix,find(THIS.multiplier));
        end
        
        USROPT.fixlevel = union(USROPT.fix,USROPT.fixlevel);
        USROPT.fixgrowth = union(USROPT.fix,USROPT.fixgrowth);
        if ~USROPT.growth
            USROPT.fixgrowth = find(THIS.nametype <= 2);
        end
    end
% dooptions2fix().

%**********************************************************************
    function dochkfornans()
        % Check for levels fixed to NaN.
        fixlevelindex = false(1,length(THIS.name));
        fixlevelindex(OPT.fixlevel) = true;
        nansstate = any(isnan(real(THIS.Assign)),3) & fixlevelindex;
        if any(nansstate)
            utils.error('model', ...
                ['Cannot fix steady-state level for this variable ', ...
                'because it is NaN: ''%s''.'], ...
                THIS.name{nansstate});
        end
        % Check for growth rates fixed to NaN.
        fixgrowthindex = false(1,length(THIS.name));
        fixgrowthindex(OPT.fixgrowth) = true;
        nansstate = any(isnan(imag(THIS.Assign)),3) & fixgrowthindex;
        if any(nansstate)
            utils.error('model', ...
                ['Cannot fix steady-state growth for this variable ', ...
                'because it is NaN: ''%s''.'], ...
                THIS.name{nansstate});
        end
    end
% dochkfornans().

%**********************************************************************
    function doinitialise()
        % Initialise persistent variables.
        
        % Run BLAZER if it has not been run yet or if user requested
        % exogenise/endogenise.
        if OPT.blocks && ...
                (isempty(THIS.nameblk) || isempty(THIS.eqtnblk) ...
                || ~isempty(OPT.endogenise) || ~isempty(OPT.exogenise))
            % The `nameblk` and `eqtnblk` properties are meant to store the
            % blazer results without any swaps. If we run the blazer with a swap,
            % we need to set the `nameblk` and `eqntblk` back to what it was
            % before at the end of the file.
            nameblk = THIS.nameblk;
            eqtnblk = THIS.eqtnblk;
            THIS = myblazer(THIS);
        end
        
        % Prepare blocks of equations/names.
        if OPT.blocks
            NAMEBLKLEVEL = THIS.nameblk;
            NAMEBLKGROWTH = THIS.nameblk;
            EQTNBLK = THIS.eqtnblk;
        else
            % If 'blocks=' false, prepare two blocks: transition and measurement.
            NAMEBLKLEVEL = cell(1,2);
            NAMEBLKGROWTH = cell(1,2);
            % All transition equations and variables.
            EQTNBLK = cell(1,2);
            NAMEBLKLEVEL{1} = find(THIS.nametype == 2);
            NAMEBLKGROWTH{1} = find(THIS.nametype == 2);
            EQTNBLK{1} = find(THIS.eqtntype == 2);
            % All measurement equations and variables.
            NAMEBLKLEVEL{2} = find(THIS.nametype == 1);
            NAMEBLKGROWTH{2} = find(THIS.nametype == 1);
            EQTNBLK{2} = find(THIS.eqtntype == 1);
        end

        NBLK = length(NAMEBLKLEVEL);
        BLOCKFUNC = cell(1,NBLK);
        % Remove variables fixed by the user.
        % Prepare function handles to evaluate individual equation blocks.
        for ii = 1 : NBLK
            % Exclude fixed levels and growth rates from the list of optimised
            % names.
            NAMEBLKLEVEL{ii} = setdiff(NAMEBLKLEVEL{ii},OPT.fixlevel);
            NAMEBLKLEVEL{ii} = setdiff(NAMEBLKLEVEL{ii},THIS.Refresh);
            NAMEBLKGROWTH{ii} = setdiff(NAMEBLKGROWTH{ii},OPT.fixgrowth);
            NAMEBLKGROWTH{ii} = setdiff(NAMEBLKGROWTH{ii},THIS.Refresh);
            if isempty(NAMEBLKLEVEL{ii}) && isempty(NAMEBLKGROWTH{ii})
                continue
            end
            % Create an anonymous function handle for each block.
            eqtn = THIS.eqtnS(EQTNBLK{ii});
            eqtn = strrep(eqtn,'exp?','exp');
            % Replace log(exp(x(...))) with x(...). This helps a lot.
            eqtn = regexprep(eqtn,'log\(exp\(x\((\d+)\)\)\)','x($1)');
            % Create a function handle used to evaluate each block of
            % equations.
            BLOCKFUNC{ii} = eval(['@(x,dx) [',eqtn{:},']']);
        end

        % Set `nameblk` and `eqtnblk` to what it was in the input model if the
        % blazer was run with a swap.
        if OPT.blocks && ...
                (~isempty(OPT.endogenise) || ~isempty(OPT.exogenise))
            THIS.nameblk = nameblk;
            THIS.eqtnblk = eqtnblk;
        end
        
        ISREFRESH = ~isempty(THIS.Refresh);
        LINK = THIS.eqtnF(THIS.eqtntype == 4);

    end
% doinitialise().

end