function S = mysimulate(THIS,S,OPT,ILOOP)
% mysimulate  [Not a public function] Backend for simulate.
%
%
% Backed IRIS function.
% No help provided.

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

OPT.display = double(OPT.display);

if ischar(OPT.dtrends)
    OPT.dtrends = ~OPT.deviation;
end

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

ny = sum(THIS.nametype == 1);
nx = size(THIS.solution{1},1);
nb = size(THIS.solution{1},2);
nf = nx - nb;
ne = sum(THIS.nametype == 3);
nalt = size(THIS.Assign,3);
nalt = size(THIS.Assign,3);
neqtn = length(THIS.eqtn);
nn = sum(THIS.nonlin);
nper = size(S.e,2);

if OPT.anticipate
    anticipatedfunc = @real;
    unanticipatedfunc = @imag;
    complexfunc = @(ant,unant) complex(ant,unant);
else
    anticipatedfunc = @imag;
    unanticipatedfunc = @real;
    complexfunc = @(ant,unant) complex(unant,ant);
end

if ILOOP <= nalt
    S.Assign = THIS.Assign(1,:,ILOOP);
    S.T = THIS.solution{1}(:,:,ILOOP);
    S.R = THIS.solution{2}(:,:,ILOOP);
    S.K = THIS.solution{3}(:,:,ILOOP);
    S.Z = THIS.solution{4}(:,:,ILOOP);
    S.H = THIS.solution{5}(:,:,ILOOP);
    S.D = THIS.solution{6}(:,:,ILOOP);
    S.U = THIS.solution{7}(:,:,ILOOP);
    S.Y = []; % Effect of non-linear add-factors.
    if S.isnonlin
        S.Y = THIS.solution{8}(:,:,ILOOP);
    end
    % Expand solution forward up to t+k if needed.
    if S.tplusk > 0
        if S.isnonlin && (ne > 0 || nn > 0)
            % Expand solution forward to t+k for both shocks and non-linear
            % add-factors.
            [S.R,S.Y] = ...
                model.myexpand(S.R,S.Y,S.tplusk,THIS.Expand{1:6});
        elseif ne > 0
            % Expand solution forward to t+k for shocks only.
            S.R = ...
                model.myexpand(S.R,[],S.tplusk,THIS.Expand{1:5},[]);
        end
    end
end

% Solution not available, return immediately.
S.nansolution = any(isnan(S.T(:)));
if S.nansolution
    return
end

% Detect segmentation by uanticipated shocks if shocks are not ignored.
% Zero all shocks otherwise.
S.segment = 1;
if ~OPT.ignoreshocks
    % Positions of unanticipated shocks (segments).
    S.segment = find(any(unanticipatedfunc(S.e) ~= 0,1));
    if isempty(S.segment) || S.segment(1) ~= 1
        S.segment = [1,S.segment];
    end
end

% Subtract deterministic trends from measurement tunes.
if ny > 0 && S.istune && OPT.dtrends
    S.ytune = S.ytune - S.W;
end

S.y = [];
S.w = [];

if S.isnonlin
    dononlinear();
else
    S.count = 0;
    S.u = [];
    if OPT.contributions
        docontributions(nper);
    else
        dosimulate(nper);
    end
    if ~isempty(S.progress)
        update(S.progress,ILOOP/nloop);
    end
end

%**************************************************************************
    function dosimulate(nper)
        if S.lastexog == 0
            % Plain simulation.
            [S.y,S.w] = timedom.simulatemean(S.T,S.R,S.K,S.Z,S.H,S.D,S.U,...
                S.a0,S.e,nper,OPT.anticipate,OPT.deviation,S.Y,S.u);
        else
            % Exogenised simulation.
            % Plain simulation first.
            [S.y,S.w] = timedom.simulatemean(S.T,S.R,S.K,S.Z,S.H,S.D,S.U,...
                S.a0,S.e,S.lastexog,OPT.anticipate,OPT.deviation,S.Y,S.u);
            % Compute multiplier matrices in the first round only. No
            % need to re-calculate the matrices in the second and further
            % rounds of non-linear simulations.
            if S.count == 0
                S.M = [ ...
                    domultipliers(true), ...
                    domultipliers(false), ...
                    ];
            end
            % Back out residuals.
            doexogenise();
            % Add anticipated residuals.
            if OPT.anticipate
                S.addeu = 1i*S.addeu;
            else
                S.addea = 1i*S.addea;
            end
            S.e(:,1:S.lastendogu) = ...
                S.e(:,1:S.lastendogu) + S.addeu;
            S.e(:,1:S.lastendoga) = ...
                S.e(:,1:S.lastendoga) + S.addea;
            % Re-simulate with residuals added.
            [S.y,S.w] = timedom.simulatemean(S.T,S.R,S.K,S.Z,S.H,S.D,S.U,...
                S.a0,S.e,nper,OPT.anticipate,OPT.deviation,S.Y,S.u);
        end
    end
% dosimulate().

%**************************************************************************
    function docontributions(nper)
        % Compute contributions of shocks and
        % init.cond.+constant
        e0 = S.e;
        S.e = zeros(size(e0,1),size(e0,2),ne+1);
        S.y = nan(ny,nper,ne+1);
        S.w = nan(nx,nper,ne+1); % := [xf;a]
        for ii = 1 : ne
            S.e(ii,:,ii) = e0(ii,:);
            [S.y(:,:,ii),S.w(:,:,ii)] = timedom.simulatemean(S.T,S.R,S.K,S.Z,S.H,S.D,S.U,...
                zeros(nb,1),S.e(:,:,ii),nper,OPT.anticipate,true);
        end
        [S.y(:,:,ne+1),S.w(:,:,ne+1)] = ...
            timedom.simulatemean(S.T,S.R,S.K,S.Z,S.H,S.D,S.U,...
            S.a0,S.e(:,:,ne+1),nper,OPT.anticipate,OPT.deviation);
    end
% docontributions().

%**************************************************************************
    function M = domultipliers(ant)
        if ant
            lastendog = S.lastendoga;
        else
            lastendog = S.lastendogu;
        end
        % M := [My(1);Mx(1);My(1);Mx(1);...];
        nnzy = nnz(S.yanchors(:,1:S.lastexog));
        nnzx = nnz(S.xanchors(:,1:S.lastexog));
        if ant
            nnze = nnz(S.eanchorsA);
        else
            nnze = nnz(S.eanchorsU);
        end
        if S.lastexog == 0 || lastendog == 0
            M = zeros(nnzy+nnzx,nnze);
            return
        end
        ma = zeros(nb,ne*lastendog);
        Tf = S.T(1:nf,:);
        Ta = S.T(nf+1:end,:);
        M = zeros(0,nnze);
        if ant
            eanchors = S.eanchorsA(:,1:lastendog);
            eanchors = eanchors(:).';
            r = S.R(:,1:ne*lastendog);
        else
            eanchors = S.eanchorsU(:,1:lastendog);
            eanchors = eanchors(:).';
            r = S.R(:,1:ne);
        end
        for t = 1 : S.lastexog
            mf = Tf*ma;
            ma = Ta*ma;
            if ant
                mf(:,(t-1)*ne+1:end) = mf(:,(t-1)*ne+1:end) + r(1:nf,:);
                ma(:,(t-1)*ne+1:end) = ...
                    ma(:,(t-1)*ne+1:end) + r(nf+1:end,:);
                r = r(:,1:end-ne);
            elseif t <= lastendog
                mf(:,(t-1)*ne+(1:ne)) = mf(:,(t-1)*ne+(1:ne)) + r(1:nf,:);
                ma(:,(t-1)*ne+(1:ne)) = ...
                    ma(:,(t-1)*ne+(1:ne)) + r(nf+1:end,:);
            end
            my = S.Z*ma;
            if t <= lastendog
                my(:,(t-1)*ne+(1:ne)) = my(:,(t-1)*ne+(1:ne)) + S.H;
            end
            yanchors = S.yanchors(:,t);
            xfanchors = S.xanchors(1:nf,t);
            xbanchors = S.xanchors(nf+1:end,t);
            M = [M; ...
                my(yanchors,eanchors); ...
                mf(xfanchors,eanchors); ...
                S.U(xbanchors,:)*ma(:,eanchors); ...
                ]; %#ok<AGROW>
        end
    end
% domultipliers().

%**************************************************************************
    function doexogenise()
        % Convert w := [xf;a] vector to x := [xf;xb] vector.
        x = S.w;
        x(nf+1:end,:) = S.U*x(nf+1:end,:);
        % Compute prediction errors.
        % pe : = [ype(1);xpe(1);ype(2);xpe(2);...].
        pe = [];
        for t = 1 : S.lastexog
            pe = [pe; ...
                S.ytune(S.yanchors(:,t),t)-S.y(S.yanchors(:,t),t); ...
                S.xtune(S.xanchors(:,t),t)-x(S.xanchors(:,t),t); ...
                ]; %#ok<AGROW>
        end
        adde = S.M \ pe;
        nnzea = nnz(S.eanchorsA(:,1:S.lastendoga));
        eindexA = S.eanchorsA(:,1:S.lastendoga);
        eindexA = eindexA(:);
        eindexU = S.eanchorsU(:,1:S.lastendogu);
        eindexU = eindexU(:);
        S.addea = zeros(ne,S.lastendoga);
        S.addeu = zeros(ne,S.lastendogu);
        S.addea(eindexA) = adde(1:nnzea);
        S.addeu(eindexU) = adde(nnzea+1:end);
    end
% doexogenise().

%**************************************************************************
    function dononlinear()
        % dononlinear  Split non-linear simulation into segments of unanticipated
        % shocks, and call `dononlinearsegment()`.
        if OPT.deviation && OPT.addsstate
            S.nonlinxbar = mytrendarray(THIS, ...
                THIS.solutionid{2},0:S.lastnonlin,false,min(iloop,nalt));
        end
        npermax = max(nper,S.segment(end)+S.lastnonlin-1);
        ea = anticipatedfunc(S.e);
        eu = unanticipatedfunc(S.e);
        if size(ea,2) < npermax
            ea(:,end+1:npermax) = 0;
            eu(:,end+1:npermax) = 0;
        end
        % Temporary arrays to store results.
        y = zeros(ny,0);
        w = zeros(nx,0);
        e = zeros(ne,0);
        S.u = zeros(nn,0);
        nsegment = length(S.segment);
        stop = zeros(1,nsegment);
        addfact = nan(nn,npermax,nsegment);
        discrep = nan(nn,nper);
        for isegment = 1 : nsegment
            % The segment dates are defined by `first` to `last`, a total of `nper1`
            % periods. These are the dates that will be added to the output data.
            % However, the actual range to be simulated can be longer because
            % `lastnonlin` (the number of non-linearised periods) may go beyond `last`.
            % The number of periods simulated is therefore `nper1max`.
            first = S.segment(isegment);
            if isegment < nsegment
                last = S.segment(isegment+1) - 1;
            else
                last = nper;
            end
            S.nonlintime = first : last;
            nper1 = length(S.nonlintime);
            nper1max = max(nper1,S.lastnonlin);
            nper1min = min(nper1,S.lastnonlin);
            % Combine anticipated shocks on the whole segment with unanticipated shocks
            % in the initial period.
            S.e = complexfunc( ...
                ea(:,first:first+nper1max-1), ...
                [eu(:,first),zeros(ne,nper1max-1)]);
            % Reset.
            S.count = 0;
            S.stop = 0;
            S.lambda = OPT.lambda;
            % Re-use addfactors from the previous segment.
            S.u(:,end+1:S.lastnonlin) = 0;
            s = sprintf('%g:%g[%g]',...
                S.zerothsegment+first, ...
                S.zerothsegment+last, ...
                S.zerothsegment+first+S.lastnonlin-1);
            S.segmentstring = sprintf('%13s',s);
            dononlinearsegment();
            dosimulate(nper1max);
            % Store results in temporary arrays.
            y = [y,S.y(:,1:nper1)]; %#ok<AGROW>
            w = [w,S.w(:,1:nper1)]; %#ok<AGROW>
            e = [e,S.e(:,1:nper1)]; %#ok<AGROW>
            % Update initial condition for next segment.
            S.a0 = S.w(nf+1:end,nper1);
            discrep(:,first+(0:nper1min-1)) = S.discrep(:,1:nper1min);
            stop = [stop,S.stop]; %#ok<AGROW>
            addfact(:,first+(0:size(S.u,2)-1),isegment) = S.u;
            % Chop off add-factors within the current segment's reported range. The
            % add-factors beyond the reported range end will be used for initial
            % condition in the next segment. Note that `u` can be shorter than `nper1`.
            S.u(:,1:nper1min) = [];
            if S.progress
                update(progress, ...
                    ((iloop-1)*nsegment+isegment)/(nloop*nsegment));
            end
        end
        S.e = e;
        S.y = y;
        S.w = w;
        
        % Output arguments from `simulate` in non-linear simulations.
        S.EXITFLAG = stop;
        S.DISCREP = discrep;
        S.ADDFACT = addfact;
    end
% dononlinear().

%**************************************************************************
    function dononlinearsegment()
        % dononlinearsegment  Non-linear simulation of one segem
        
        S.histmindiscrep = Inf;
        S.histminu = [];
        S.histmincount = NaN;
        
        while true
            dosimulate(S.lastnonlin);
            dodiscrepancy();
            
            if S.maxdiscrep < S.histmindiscrep
                S.histmindiscrep = S.maxdiscrep;
                S.histminu = S.u;
                S.histmincount = S.count;
            end
            
            % Report discrepancies in this iteration if requested or if
            % this is the final iteration.
            if OPT.display > 0 && mod(S.count,OPT.display) == 0
                doreport();
            end
            
            if S.maxdiscrep <= OPT.tolerance
                S.stop = 1;
            elseif S.count >= OPT.maxiter;
                S.stop = -1;
            elseif ~isfinite(S.maxdiscrep)
                S.stop = -2;
            end
            
            if S.stop ~= 0
                if S.maxdiscrep > S.histmindiscrep
                    S.u = S.histminu;
                    dosimulate(S.lastnonlin);
                    if OPT.display > 0
                        doreportreverse();
                        dodiscrepancy();
                    end
                end
                %use.count = Inf;
                if OPT.display > 0
                    doreport();
                    fprintf('\n');
                end
                break
            end
            
            if S.maxdiscrep < 1.5*S.histmindiscrep
                addu = S.discrep;
                addu(abs(addu) <= OPT.tolerance) = 0;
                addu = S.lambda .* addu;
                S.u = S.u - addu;
            else
                % If the current discrepancy is twice the historical minimum (or more),
                % reverse the process to the historical minimum, and reduce `lambda`.
                S.u = S.histminu;
                S.lambda = S.lambda * OPT.reducelambda;
                doreportreverse();
                doreportlambdareduction();
            end
            S.count = S.count + 1;
        end
        
        % Failed to converge.
        if S.stop < 0
            if OPT.error
                messagefunc = @utils.error;
            else
                messagefunc = @utils.warning;
            end
            
            switch S.stop
                case -1
                    messagefunc('model', ...
                        ['Non-linear simulation #%g reached ', ...
                        'the max number of iterations ', ...
                        'without achieving convergence.'], ...
                        iloop);
                case -2
                    messagefunc('model', ...
                        ['Non-linear simulation #%g exploded ', ...
                        'to Inf or crashed at NaN.'], ...
                        iloop);
            end
        end
        
        function dodiscrepancy()
            tt = 1 : S.lastnonlin;
            LL = real(S.Assign);
            % Set up the vector of [xf;xb] including initial condition.
            xx = [[nan(nf,1);S.a0],S.w(:,tt)];
            xx(nf+1:end,:) = S.U*xx(nf+1:end,:);
            if OPT.deviation && OPT.addsstate
                xx = xx + S.nonlinxbar;
            end
            % Set up the vector of shocks including initial condition.
            ee = real(S.e(:,tt)) + imag(S.e(:,tt));
            ee = [zeros(ne,1),ee];
            % Get the current parameter values.
            pp = S.Assign(1,THIS.nametype == 4);
            d = zeros(neqtn,S.lastnonlin);
            for j = find(THIS.nonlin)
                tt = find(S.qanchors(j,:));
                if isempty(tt)
                    continue
                end
                try
                    d(j,tt) = THIS.eqtnN{j}(xx,ee,pp,1+tt,LL);
                catch Error
                    d(j,:) = NaN;
                    utils.warning('model', ...
                        ['Error evaluating non-linear equations.', ...
                        '\nMatlab says %s'], ...
                        Error.message);
                end
            end
            S.discrep = d(THIS.nonlin,:);
            % Maximum discrepancy and max addfactor.
            S.maxdiscrep2 = max(abs(S.discrep),[],2);
            S.maxdiscrep = max(S.maxdiscrep2);
            S.maxaddfactor2 = max(abs(S.u),[],2);
            S.maxaddfactor = max(S.maxaddfactor2);
        end
        % dodiscrepancy().
        
        function doreport()
            % Report nonlin simulation iteration.
            maxdiscrepeqtn = ...
                findnaninf(S.maxdiscrep2,S.maxdiscrep,1,'first');
            maxaddfactoreqtn = ....
                findnaninf(S.maxaddfactor2,S.maxaddfactor,1,'first');
            if S.count == 0 && S.stop == 0
                % This is the very first report line printed. Print the header first.
                fprintf(...
                    '%13s %6.6s %12.12s %-20.20s %7.7s %12.12s %-20.20s\n',...
                    'Segment','Iter','Max.discrep','Equation','Lambda', ...
                    'Max.addfact','Equation' ...
                    );
            end
            count = sprintf(' %5g',S.count);
            if S.stop ~= 0
                count = strrep(count,' ','=');
            end
            lambda = sprintf('%7g',S.lambda);
            maxdiscrep = sprintf('%12g',S.maxdiscrep);
            maxdiscreplabel = S.label{maxdiscrepeqtn};
            maxdiscreplabel = strfun.ellipsis(maxdiscreplabel,20);
            maxaddfactor = sprintf('%12g',S.maxaddfactor);
            maxaddfactorlabel = S.label{maxaddfactoreqtn};
            maxaddfactorlabel = strfun.ellipsis(maxaddfactorlabel,20);
            % Print current report line.
            fprintf(...
                '%s %s %s %s %s %s %s\n',...
                S.segmentstring,count, ...
                maxdiscrep,maxdiscreplabel,lambda, ...
                maxaddfactor,maxaddfactorlabel ...
                );
        end
        % doreport.
        
        function doreportreverse()
            fprintf('  Reversing to iteration %g.\n', ...
                S.histmincount);
        end
        % doreportreverse().
        
        function doreportlambdareduction()
            fprintf('  Lambda reduced to %g.\n', ...
                S.lambda);
        end
        % doreportlambdareduction().
        
    end
% dononlinear().

end