function [FLAG,varargout] = chksstate(THIS,varargin)
% chksstate  Check if equations hold for currently assigned steady0state values.
%
% Syntax
% =======
%
%     [FLAG,LIST] = chksstate(M,...)
%     [FLAG,DISCREP,LIST] = chksstate(M,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Model object.
%
% Output arguments
% =================
%
% * `FLAG` [ `true` | `false` ] - True if discrepancy between LHS and RHS is
% smaller than `'tolerance='` in each equation.
%
% * `DISCREP` [ numeric ] - Discrepancy between LHS and RHS in each
% equation evaluated at two consecutive times.
%
% * `LIST` [ cellstr ] - List of equations in which the discrepancy between
% LHS and RHS is greater than `'tolerance='`.
%
% Options
% ========
%
% * `'error='` [ `true` | *`false`* ] - Throw an error if one or more equations
% do not hold.
%
% * `'refresh='` [ *`true`* | `false` ] - Refresh dynamic links before
% evaluating the equations.
%
% * `'sort='` [ `true` | *`false`* ] - Undocumented option.
%
% * `'sstateEqtn='` [ `true` | *`false`* ] - If `false`, the dynamic model
% equations will be checked; if `true`, the steady-state versions of the
% equations (wherever available) will be checked.
%
% * `'tolerance='` [ numeric | `getrealsmall()` ] - Tolerance.
%
% * `'warning='` [ *`true`* | `false` ] - Display warnings produced by this
% function.
%
% Description
% ============
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright 2007-2012 Jaromir Benes.

opt = passvalopt('model.chksstate',varargin{:});

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

% Refresh dynamic links.
if opt.refresh && ~isempty(THIS.Refresh)
    THIS = refresh(THIS);
end

if opt.warning
    chk(THIS,Inf,'parameters','sstate','log');
end

t = THIS.tzero;
mint = 1 - t;
if issparse(THIS.occur)
    nt = size(THIS.occur,2) / length(THIS.name);
else
    nt = size(THIS.occur,3);
end
maxt = nt - t;
nalt = size(THIS.Assign,3);

FLAG = false([1,nalt]);
list = cell([1,nalt]);
mval = nan([sum(THIS.eqtntype <= 2),nalt]);

if ~opt.sstateeqtn
    discrep = nan([sum(THIS.eqtntype <= 2),2,nalt]);
    eqtn = THIS.eqtnF(THIS.eqtntype <= 2);
else
    discrep = nan([sum(THIS.eqtntype <= 2),nalt]);
    eqtn = THIS.eqtnS(THIS.eqtntype <= 2);
    % Remove the `exp` transformation applied to log variables. Otherwise,
    % we would have to log the sstate vector, and it would get de-logged
    % again in the equations.
    eqtn = strrep(eqtn,'exp?','');
    % Create anonymous funtions for sstate equations.
    for i = 1 : length(eqtn)
        eqtn{i} = str2func(['@(x,dx) ',eqtn{i}]);
    end
end

tvec = mint : maxt+1;
for ialt = 1 : nalt
    if ~opt.sstateeqtn
        % Check the full equations in two consecutive periods. This way we
        % can detect both incorrect levels and growth rates.
        x = mytrendarray(THIS,1:length(THIS.name),tvec,true,ialt);
        x = shiftdim(x,-1);
        % Levels and steady-state leveles at time t.
        x1 = x(1,:,1:end-1);
        L1 = x1(1,:,t);
        % Levels and steady-state levels at time t+1.
        x2 = x(1,:,2:end);
        L2 = x2(1,:,t);
        mval1 = (cellfun(@(fcn) fcn(x1,t,L1),eqtn)).';
        mval2 = (cellfun(@(fcn) fcn(x2,t,L2),eqtn)).';
        discrep(:,:,ialt) = [mval1,mval2];
        mval(:,ialt) = max(abs([mval1,mval2]),[],2);
    else
        x = real(THIS.Assign(1,:,ialt));
        dx = imag(THIS.Assign(1,:,ialt));
        dx(THIS.log & dx == 0) = 1;
        mval(:,ialt) = (cellfun(@(fcn) fcn(x,dx),eqtn)).';
        discrep(:,ialt) = mval(:,ialt);
    end
    index = abs(mval(:,ialt)) <= opt.tolerance;
    FLAG(ialt) = all(index == true);
    if nargout > 1 || opt.error
        list{ialt} = THIS.eqtn(~index);
    end
end

if opt.error && any(~FLAG)
    tmp = {};
    for i = find(~FLAG)
        for j = 1 : length(list{i})
            tmp{end+1} = i; %#ok<AGROW>
            tmp{end+1} = list{i}{j}; %#ok<AGROW>
        end
    end
    utils.error('model', ...
        'Steady-state error in this equation in%s: ''%s''.', ...
        tmp{:});
end

if opt.sort
    sortlist = cell([1,nalt]);
    for ialt = 1 : nalt
        [ans,index] = sort(mval(:,ialt),1,'descend'); %#ok<NOANS,ASGLU>
        discrep(:,:,ialt) = discrep(index,:,ialt);
        sortlist{ialt} = THIS.eqtn(index);
    end
end

if nalt == 1
    list = list{1};
    if opt.sort
        sortlist = sortlist{1};
    end
end

if nargout == 2
    varargout{1} = list;
elseif nargout > 2
    varargout{1} = discrep;
    varargout{2} = list;
    if opt.sort
        varargout{3} = sortlist;
    end
end

end