function [THIS,DATA,INDEX,CRIT] = sort(THIS,DATA,SORTBY,varargin)
% sort  Sort SVAR parameterisations by squared distance of shock reponses to median.
%
% Syntax
% =======
%
%     [B,DATA,INDEX,CRIT] = sort(A,DATA,SORTBY,...)
%
% Input arguments
% ================
%
% * `A` [ SVAR ] - SVAR object with multiple parameterisations that will
% be sorted.
%
% * `DATA` [ tseries | struct | empty ] - SVAR data (endogenous variables
% and structural shocks); the structural shocks will be re-ordered
% according to the SVAR parameterisations. You can leave the input argument
% `DATA` empty.
%
% * `SORTBY` [ char ] - Text string that will be evaluated to compute the
% criterion by which the parameterisations will be sorted; see Description
% for how to write `SORTBY`.
%
% Output arguments
% =================
%
% * `B` [ SVAR ] - SVAR object with parameterisations sorted by the
% specified criterion.
%
% * `DATA` [ tseries | struct | empty ] - SVAR data with the structural
% shocks re-ordered to correspond to the order of parameterisations.
%
% * `INDEX` [ numeric ] - Vector of indices so that `B = A(INDEX)`.
%
% * `CRIT` [ numeric ] - The value of the criterion based on the string
% `SORTBY` for each parameterisation.
%
% Options
% ========
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% Description
% ============
%
% The individual parameterisations within the SVAR object `A` are sorted by
% the sum squared distances of selected shock responses to the respective
% median reponses. Formally, the following criterion is evaluated for each
% parameterisation
%
% $$ \sum_{i\in I,j\in J,k\in K} \left[ S_{i,j}(k) - M_{i,j}(k) \right]^2 $$
%
% where $S_{i,j}(k)$ denotes the response of the i-th variable to the j-th
% shock in period k, and $M_{i,j}(k)$ is the median responses. The sets of
% variables, shocks and periods, i.e. `I`, `J`, `K`, respectively, over
% which the summation runs are determined by the user in the `SORTBY`
% string.
%
% How do you select the shock responses that enter the criterion in
% `SORTBY`? The input argument `SORTBY` is a text string that refers to
% array `S`, whose element `S(i,j,k)` is the response of the i-th
% variable to the j-th shock in period k.
%
% Example
% ========
%
% Sort the parameterisations by squared distance to median of shock
% responses of all variables to the first shock in the first four periods.
% The parameterisation that is closest to the median responses
%
%     S2 = sort(S1,[],'S(:,1,1:4)')
%

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

opt = passvalopt('SVAR.sort',varargin{:});

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

[ny,p,nalt] = sizeof(THIS); %#ok<ASGLU>

% Handle residuals.
isdata = nargin > 1 && nargout > 1 && ~isempty(DATA);
if isdata
    % Get data.
    [outputformat,range,y,e] = VAR.datarequest(THIS,DATA,Inf,opt);
end

% Look for the simulation horizon and the presence of asymptotic responses
% in the `SORTBY` string.
[h,isy] = mywoonvav(THIS,SORTBY);

if opt.progress
    pbar = progressbar('IRIS SVAR.sort progress');
end

XX = [];
for ialt = 1 : nalt
    [S,Y] = dosimulate(); %#ok<NASGU,ASGLU>
    doevalsort();
    if opt.progress
        update(pbar,ialt/nalt);
    end
end

INDEX = dosort();
THIS = mysubsalt(THIS,INDEX);

if isdata
    e = e(:,:,INDEX);
    DATA = VAR.outputdata(THIS,outputformat,range,[y;e],[], ...
        [THIS.ynames,THIS.enames]);
end

% Nested functions.

%**************************************************************************
    function [S,Y] = dosimulate()
        % Simulate the test statistics.
        S = zeros(ny,ny,0);
        Y = nan(ny,ny,1);
        % Impulse responses.
        if h > 0
            S = timedom.var2vma(THIS.A(:,:,ialt),THIS.B(:,:,ialt),h);
        end
        % Asymptotic impulse responses.
        if isy
            A = poly.var2poly(THIS.A(:,:,ialt));
            C = sum(A,3);
            Y = C\THIS.B(:,:,ialt);
        end
    end
% dosimulate().

%**************************************************************************
    function doevalsort()
        % Evalutate the sort criterion.
        try
            X = eval(SORTBY);
            XX = [XX,X(:)];
        catch err
            utils.error('VAR', ...
                ['Error evaluating the sort string ''%s''.', ...
                '\nMatlab says ''%s'''], ...
                SORTBY,err.message);
        end
    end
% doevalsort().

%**************************************************************************
    function index = dosort()
        % Sort by the distance from median.
        n = size(XX,2);
        if n > 0
            MM = median(XX,2);
            CRIT = nan(1,n);
            for ii = 1 : n
                CRIT(ii) = sum((XX(:,ii) - MM).^2 / n);
            end
            [CRIT,index] = sort(CRIT,'ascend');
        end
    end
% dosort().

end