function varargout = isstationary(this,varargin)
% isstationary  True if model or specified combination of variables is stationary.
%
% Syntax
% =======
%
%     FLAG = isstationary(M)
%     FLAG = isstationary(M,EXPRESSION)
%
% Input arguments
% ================
%
% * `M` [ model ] - Model object.
%
% * `EXPRESSION` [ char ] - Combination of transition variables to be tested.
%
% Output arguments
% =================
%
% * `'flag='` [ `true` | `false` ] - True if the model (if called without a
% second input argument) or the specified combination of variables (if
% called with a second input argument) is stationary.
%
% Description
% ============
%
% Example
% ========
%

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

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

if isempty(varargin)
    % Called flag = isstationary(model).
    if isempty(this.solution{1})
        varargout{1} = NaN;
    else
        nb = size(this.solution{1},2);
        realsmall = getrealsmall();
        varargout{1} = ...
            permute(all(abs(this.eigval(1,1:nb,:)) < 1-realsmall,2), ...
            [1,3,2]);
    end
else
    % Called [flag,...] = isstationary(model,expression,...).
    [varargout{1:length(varargin)}] = ...
        xxiscointegrated(this,varargin{:});
end

end

% Subfunctions.

%**************************************************************************
function varargout = xxiscointegrated(this,varargin)

realsmall = getrealsmall();
xvector = get(this,'xvector');
xvector = strrep(xvector,'(','\(');
xvector = strrep(xvector,')','\)');
nxvector = length(xvector);

varargout = cell(1,length(varargin));
for iarg = 1 : length(varargin)
    flag = dothisargin(varargin{iarg});
    varargout{iarg} = dothisargin(varargin{iarg});
end

%**************************************************************************
    function flag = dothisargin(comb)
        present = false(size(xvector));
        comb = regexprep(comb,'\s+','');
        for j = 1 : nxvector
            pattern = ['\<',xvector{j},'\>'];
            present(j) = ~isempty(regexp(comb,pattern,'once'));
            if present(j)
                replace = sprintf('x(%g)',j);
                comb = regexprep(comb,pattern,replace);
            end
        end
        w = xxevalcomb(comb,present);
        [nx,nb,nalt] = size(this.solution{1});
        nf = nx - nb;
        flag = false(1,nalt);
        for iloop = 1 : nalt
            Tf = this.solution{1}(1:nf,:,iloop);
            U = this.solution{7}(:,:,iloop);
            nunit = sum(abs(abs(this.eigval(1,1:nb,iloop)) - 1) <= realsmall);
            test = w*[Tf(:,1:nunit);U(:,1:nunit)];
            flag(iloop) = all(abs(test) <= realsmall);
        end
    end
% xxthisargin().

end
% xxiscointegrated().

%**************************************************************************
function coeff = xxevalcomb(expr,present)

nxvector = length(present);
x = zeros(1,nxvector);
const = eval(expr);
coeff = zeros(1,nxvector);
for i = find(present)
    x(i) = 1;
    tmp = eval(expr);
    coeff(i) = tmp - const;
    x(i) = 0;
end

end
% xxevalcomb().