function [OUTP,EXITFLAG,ADDFACT,DISCREP] = simulate(THIS,INP,RANGE,varargin)
% simulate  Simulate model.
%
% Syntax
% =======
%
%     S = simulate(M,D,RANGE,...)
%     [S,FLAG,ADDF,DISCREP] = simulate(M,D,RANGE,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Solved model object.
%
% * `D` [ struct | cell ] - Input database or datapack from which the
% initial conditions and shocks from within the simulation range will be
% read.
%
% * `RANGE` [ numeric ] - Simulation range.
%
% Output arguments
% =================
%
% * `S` [ struct | cell ] - Database with simulation results.
%
% Output arguments in simulations with a non-linear plan
% =======================================================
%
% * `FLAG` [ cell ] - Cell array with exit flags for non-linearised
% simulations.
%
% * `ADDF` [ cell ] - Cell array of tseries with final add-factors added to
% first-order approximate equations to make non-linear equations hold.
%
% * `DISCREP` [ cell ] - Cell array of tseries with final discrepancies
% between LHS and RHS in equations earmarked for non-linear simulations by
% a double-equal sign.
%
% Options
% ========
%
% * `'anticipate='` [ *`true`* | `false` ] - If `true`, real future shocks are
% anticipated, imaginary are unanticipated; vice versa if `false`.
%
% * `'contributions='` [ `true` | *`false`* ] - Decompose the simulated paths
% into contributions of individual shocks.
%
% * `'deviation='` [ `true` | *`false`* ] - Treat input and output data as
% deviations from balanced-growth path.
%
% * `'dbExtend='` [ `true` | *`false`* | struct ] - Use the function
% `dbextend` to combine the simulated output data with the input database,
% or with another database, at the end.
%
% * `'dTrends='` [ *'auto'* | `true` | `false` ] - Add deterministic trends to
% measurement variables.
%
% * `'ignoreShocks='` [ `true` | *`false`* ] - Read only initial conditions from
% input data, and ignore any shocks within the simulation range.
%
% * `'plan='` [ plan ] - Specify a simulation plan to swap endogeneity
% and exogeneity of some variables and shocks temporarily, and/or to
% simulate some of the non-linear equations accurately.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% Options for models with non-linearised equations
% =================================================
%
% * `'addSstate='` [ *`true`* | `false` ] - Add steady state levels to simulated
% paths before evaluating non-linear equations; this option is used only if
% `'deviation=' true`.
%
% * `'display='` [ *`true`* | `false` | numeric ] - Display one-line report on each
% iteration; if `'display='` is a number `n`, display a report after every
% `n` iterations and a final report.
%
% * `'error='` [ `true` | *`false`* ] - Throw an error whenever a non-linear
% simulation fails converge; if `false`, only an warning will display.
%
% * `'lambda='` [ numeric | *`1`* ] - Step size (between `0` and `1`)
% for add factors added to non-linearised equations in every iteration.
%
% * `'reduceLambda='` [ numeric | *`0.5`* ] - Factor (between `0` and
% `1`) by which `lambda` will be multiplied if the non-linear simulation
% gets on an divergence path.
%
% * `'maxIter='` [ numeric | *`100`* ] - Maximum number of iterations.
%
% * `'tolerance='` [ numeric | *`1e-5`* ] - Convergence tolerance.
%
% Description
% ============
%
% The time series in the output database, `S`, are are defined on the
% simulation range, `RANGE`, plus include any necessary initial conditions
% (for those variables that occur with lags in the model).
%
% Example
% ========
%

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

% Parse required inputs.
p = inputParser();
p.addRequired('m',@ismodel);
p.addRequired('data',@(x) isstruct(x) || iscell(x));
p.addRequired('range',@isnumeric);
p.parse(THIS,INP,RANGE);

% Parse options.
opt = passvalopt('model.simulate',varargin{:});

use = struct();

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);
RANGE = RANGE(1) : RANGE(end);
nper = length(RANGE);
neqtn = length(THIS.eqtn);

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

% Simulation plan.
use.isplan = isa(opt.plan,'plan');
use.istune = use.isplan && nnzendog(opt.plan) > 0 && nnzexog(opt.plan) > 0;
use.isnonlinplan = use.isplan && nnznonlin(opt.plan) > 0;
use.isnonlinopt = ~isempty(opt.nonlinearise) && opt.nonlinearise > 0;
use.isnonlin = use.isnonlinplan || use.isnonlinopt;

% Check for option conflicts.
dochkconflicts();

% Bkw compatibility: obsolete option names.
if ~isempty(opt.ignoreresiduals)
    opt.ignoreshocks = opt.ignoreresiduals;
end
if ~isempty(opt.lambdafactor)
    opt.reducelambda = opt.lambdafactor;
end

% Find position of last logical true, or return zero.
findlastfunc = @(x) max([0,find(any(any(x,3),1),1,'last')]);

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

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

% Get initial condition for alpha.
% alpha is always expanded to match nalt within `datarequest`.
[ainit,xinit,naninit] = datarequest('init',THIS,INP,RANGE);
if ~isempty(naninit)
    if isnan(opt.missing)
        naninit = unique(naninit);
        utils.error('model', ...
            'Initial condition not available for this variable: ''%s''.', ...
            naninit{:});
    else
        ainit(isnan(ainit)) = opt.missing;
    end
end
ninit = size(ainit,3);

% Get shocks; both reals and imags are checked for NaNs within
% `datarequest`.
if ~opt.ignoreshocks
    Ee = datarequest('e',THIS,INP,RANGE);
    use.lastea = findlastfunc(anticipatedfunc(Ee) ~= 0);
    use.lasteu = findlastfunc(unanticipatedfunc(Ee) ~= 0);
    nshock = size(Ee,3);    
else
    use.lastea = 0;
    use.lasteu = 0;
    nshock = 0;
end

% Simulation range and plan range must be identical.
if use.isplan
    [yanchors,xanchors,eareal,eaimag,~,~,use.qanchors] = ...
        myanchors(THIS,opt.plan,RANGE);
    if use.isnonlinplan
        use.label = opt.plan.qlist;
    end
end

% Nonlinearised simulation through the option `'nonlinearise='`.
if use.isnonlinopt
    if isnumericscalar(opt.nonlinearise) && isround(opt.nonlinearise)
        qstart = 1;
        qend = opt.nonlinearise;
    else
        qstart = round(opt.nonlinearise(1) - RANGE(1) + 1);
        qend = round(opt.nonlinearise(end) - RANGE(1) + 1);
    end
    use.qanchors = false(neqtn,max(nper,qend));
    use.qanchors(THIS.nonlin,qstart:qend) = true;
    use.label = myget(THIS,'canbenonlinearised');
end

if use.istune
    use.yanchors = yanchors;
    use.xanchors = xanchors;
    if opt.anticipate
        use.eanchorsA = eareal;
        use.eanchorsU = eaimag;
    else
        use.eanchorsA = eaimag;
        use.eanchorsU = eareal;
    end
    use.lastendoga = findlastfunc(use.eanchorsA);
    use.lastendogu = findlastfunc(use.eanchorsU);
    use.lastexog = findlastfunc([use.yanchors;use.xanchors]);
    % Get actual values for exogenised data points.
    Yy = datarequest('y',THIS,INP,RANGE);
    Xx = datarequest('x',THIS,INP,RANGE);
    % Check for NaNs in exogenised variables.
    dochknanexog();
    % Check the number of exogenised and endogenised data points
    % (exogenising must always be an exactly determined system).
    dochkdetermined();
    ntune = max(size(Yy,3),size(Xx,3));
else
    ntune = 0;
    use.lastendoga = 0;
    use.lastendogu = 0;
    use.lastexog = 0;
end

% Total number of cycles.
nloop = max([1,nalt,ninit,nshock,ntune]);

if use.isnonlin
    use.lastnonlin = findlastfunc(use.qanchors);
    % The field `zerothsegment` is used by the Kalman filter to report the
    % correct period.
    use.zerothsegment = 0;
    % Prepare output arguments for non-linear simulations.
    EXITFLAG = cell(1,nloop);
    ADDFACT = cell(1,nloop);
    DISCREP = cell(1,nloop);
    dochknonlinconflicts();
else
    % Output arguments for non-linear simulations.
    use.lastnonlin = 0;
    use.display = 0;
    EXITFLAG = {};
    ADDFACT = {};
    DISCREP = {};
end

if opt.contributions
    noutput = ne + 1;
else
    noutput = nloop;
end
xrange = RANGE(1)-1 : RANGE(end);
nxper = length(xrange);
hdata = myhdatainit(THIS,nxper,noutput);

% Maximum expansion needed.
use.tplusk = max([1,use.lastea,use.lastendoga,use.lastnonlin]) - 1;

%**************************************************************************
% Main loop.

nansolution = false(1,nloop);
use.progress = opt.progress && use.display == 0;

if use.progress
    use.progress = progressbar('IRIS model.simulate progress');
else
    use.progress = [];
end

for iloop = 1 : nloop
    
    % Get current initial condition for the transformed state vector, and
    % current shocks.
    use.a0 = ainit(:,1,min(iloop,end));
    if ~opt.ignoreshocks
        use.e = Ee(:,:,min(iloop,end));
    else
        use.e = zeros(ne,nper);
    end
    
    % Current tunes on measurement and transition variables.
    use.ytune = [];
    use.xtune = [];
    if use.istune
        use.ytune = Yy(:,:,min(iloop,end));
        use.xtune = Xx(:,:,min(iloop,end));
    end
    
    % Compute deterministic trends if requested. We don't compute the dtrends
    % inside `mysimulate` because they are dealt with differently when called
    % from within the Kalman filter.
    use.W = [];
    if ny > 0 && opt.dtrends
        use.W = mydtrendsrequest(THIS,'range',RANGE,iloop);
    end

    % CALL BACKEND SIMULATION FUNCTION
    %==================================
    use = mysimulate(THIS,use,opt,iloop);
    
    % Simulation didn't run because solution not available, return immediately.
    if use.nansolution
        nansolution(iloop) = true;
        continue
    end

    % Extra output arguments in non-linear simulations.
    if use.isnonlin
        EXITFLAG{iloop} = use.EXITFLAG;
        DISCREP{iloop} = use.DISCREP;
        ADDFACT{iloop} = use.ADDFACT;
    end
    
    % Add measurement detereministic trends.
    if ny > 0 && opt.dtrends
        % Add to simulation.
        % When contributions == true, add to last simulation.
        use.y(:,:,end) = use.y(:,:,end) + use.W;
    end
    
    % Initial condition for the original state vector.
    use.x0 = xinit(:,1,min(iloop,end));
    
    % Assign output data.
    doassignoutput();

    % Add equation labels to add-factor and discrepancy series.
    if use.isnonlin && nargout > 2
        label = use.label;
        nsegment = length(use.segment);
        ADDFACT{iloop} = tseries(RANGE(1), ...
            permute(ADDFACT{iloop},[2,1,3]),label(1,:,ones(1,nsegment)));
        DISCREP{iloop} = tseries(RANGE(1), ...
            permute(DISCREP{iloop},[2,1,3]),label);
    end
    
end
% End of main loop.

%**************************************************************************
% Post mortem.

% Report solutions not available.
if any(nansolution)
    utils.warning('model', ...
        '#Solution_not_available', ...
        sprintf(' #%g',find(nansolution)));
end

% Convert hdataobj to struct. The comments assigned to the output series
% depend on whether this is a contributions simulation or not.
OUTP = myhdata2tseries(THIS,hdata,xrange);
if opt.contributions
    % Change comments.
end

% Add parameters to output database.
for i = find(THIS.nametype == 4)
    OUTP.(THIS.name{i}) = permute(THIS.Assign(1,i,:),[1,3,2]);
end

if isequal(opt.dbextend,true)
    OUTP = dbextend(INP,OUTP);
elseif isstruct(opt.dbextend)
    OUTP = dbextend(opt.dbextend,OUTP);
end

% Nested functions.


%**************************************************************************
    function dochknanexog()
        % Check for NaNs in exogenised variables.
        index1 = [use.yanchors;use.xanchors];
        index2 = [any(~isfinite(Yy),3);any(~isfinite(Xx),3)];
        index3 = [any(imag(Yy) ~= 0,3);any(imag(Xx) ~= 0,3)];
        index = any(index1 & (index2 | index3),2);
        if any(index)
            list = [THIS.solutionvector{1:2}];
            utils.error('model', ...
                'This variable is exogenised to NaN, Inf or complex number: ''%s''.', ...
                list{index});
        end
    end
% dochknanexog().

%**************************************************************************
    function dochkdetermined()
        if nnzexog(opt.plan) ~= nnzendog(opt.plan)
            utils.warning('model', ...
                ['Number of exogenised data points (%g) does not match ', ...
                'number of endogenised data points (%g).'], ...
                nnzexog(opt.plan),nnzendog(opt.plan));
        end
    end
% dochkdetermined().

%**************************************************************************
    function doassignoutput()
        n = size(use.w,3);
        xf = [nan(nf,1,n),use.w(1:nf,:,:)];
        xb = use.w(nf+1:end,:,:);
        for ii = 1 : n
            xb(:,:,ii) = use.U*xb(:,:,ii);
        end
        tmpinit = zeros(nb,1,n);
        tmpinit(:,1,end) = use.x0;
        xb = [tmpinit,xb];    
        % Columns to place results in output data.
        if opt.contributions
            cols = 1 : ne+1;
        else
            cols = iloop;
        end
        % Add current results to output data.
        myhdataassign(THIS,hdata,cols, ...
            [nan(ny,1,n),use.y], ...
            [xf;xb], ...
            [nan(ne,1,n),use.e]);
    end
% doassignoutput().

%**************************************************************************
    function dochkconflicts()
        % The option `'contributions='` option cannot be used with the `'plan='`
        % option or with multiple parameterisations.
        if opt.contributions
            if use.istune || use.isnonlin
                utils.error('model', ...
                    ['Cannot run SIMULATE with ''CONTRIBUTIONS='' true ', ...
                    'and ''PLAN='' non-empty.']);
            end
            if nalt > 1
                utils.error('model', ...
                    ['Cannot run SIMULATE with ''CONTRIBUTIONS='' true ', ...
                    'on multiple parameterisations or multiple data sets.']);
            end
        end
    end
% dochkconflicts().

%**************************************************************************
    function dochknonlinconflicts()
        if use.lastendogu > 0
            utils.error('model', ...
                ['Non-linearised simulations cannot have ', ...
                'unanticipated endogenised shocks.']);
        end
        if (use.lastendogu > 0 || use.lastendoga > 0) ...
                && use.lasteu > 0
            utils.error('model', ...
                ['Non-linearised simulations cannot combine ', ...
                'unanticipated shocks and endogenised shocks.']);
        end
    end
% dochknonlinconflicts().

end