function [F,Fi,delta,freq,G,step] = fisher(THIS,nper,plist,varargin)
% fisher  Approximate Fisher information matrix in frequency domain.
%
% Syntax
% =======
%
%     [F,Fi,delta,freq] = fisher(m,nper,list,...)
%
% Input arguments
% ================
%
% * `m` [ model ] - Solved model object.
%
% * `nper` [ numeric ] - Length of the hypothetical range for which the
% Fisher information will be computed.
%
% * `list` [ cellstr ] - List of parameters with respect to which the
% likelihood function will be differentiated.
%
% Output arguments
% =================
%
% * `F` [ numeric ] - Approximation of the Fisher information matrix.
%
% * `Fi` [ numeric ] - Contributions of individual frequencies to the total
% Fisher information matrix.
%
% * `delta` [ numeric ] - Kronecker delta by which the contributions in
% `Fi` need to be multiplied to sum up to `F`.
%
% * `freq` [ numeric ] - Vector of frequencies at which the Fisher
% information matrix is evaluated.
%
% Options
% ========
%
% * `'deviation='` [ *`true`* | `false` ] - Exclude the steady state effect at
% zero frequency.
%
% * `'exclude='` [ char | cellstr | *empty* ] - List of measurement
% variables that will be excluded from the likelihood function.
%
% * `'percent='` [ `true` | *`false`* ] - Report the overall Fisher matrix `F`
% as Hessian w.r.t. the log of variables; the interpretation is that the
% Fisher matrix then describes the changes in the log-likelihood function
% in reponse to percent, not absolute, changes in parameters.
%
% * `'progress='` [ `true` | *`false`* ] - Display progress bar in the command
% window.
%
% * `'refresh='` [ *`true`* | `false` ] - Refresh dynamic links in each
% differentiation step.
%
% * `'solve='` [ *`true`* | `false` ] - Re-solve model in each differentiation
% step.
%
% * `'sstate='` [ `true` | *`false`* | cell ] - Re-compute steady state in each
% differentiation step; if the model is non-linear, you can pass in a cell
% array with opt used in the `sstate` function.
%
% Description
% ============
%
% Example
% ========
%

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

% Validate required input arguments.
p = inputParser();
p.addRequired('this',@(x) isa(x,'model'));
p.addRequired('nper',@(x) isnumeric(x) && length(x) == 1);
p.addRequired('plist',@(x) iscellstr(x) || ischar(x));
p.parse(THIS,nper,plist);

% Read and validate optional input arguments.
opt = passvalopt('model.fisher',varargin{:});

% Old option 'sstateeffect' maintained for bkw compatibility.
if ~isempty(opt.sstateeffect)
    utils.warning('model',...
        ['The option ''sstateeffect'' in FISHER is deprecated. ', ...
        'Use ''deviation'' instead.']);
    opt.deviation = ~opt.sstateeffect;
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);

% Process the 'exclude' option.
exclude = false(ny,1);
if ~isempty(opt.exclude)
    if ischar(opt.exclude)
        opt.exclude = regexp(opt.exclude,'\w+','match');
    end
    ylist = THIS.name(THIS.nametype == 1);
    for i = 1 : length(opt.exclude)
        index = strcmp(ylist,opt.exclude{i});
        exclude(index) = true;
    end
end

% Get parameter cellstr list from a char list.
if ischar(plist)
    plist = regexp(plist,'\w+','match');
end

% Initialise persistent calls for non-lin steady-state solver.
opt.sstate = mysstateopt(THIS,opt.sstate);
clear('model/mysstatenonlin');

opt.fastsolve = true;

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

ny = ny - sum(exclude);
if ny == 0
    utils.warning('model', ...
        'No measurement variables included in computing Fisher matrix.');
end

ylog = THIS.log(THIS.nametype == 1);
ylog(exclude) = [];

[assignPos,stdcorrPos] = mynameposition(THIS,plist,'error');
assignNan = isnan(assignPos);
stdcorrNan = isnan(stdcorrPos);

nplist = length(plist);
nfreq = floor(nper/2) + 1;
freq = 2*pi*(0:nfreq-1)/nper;
realsmall = getrealsmall();

% Kronecker delta vector.
% Different for even or odd number of periods.
delta = ones(1,nfreq);
if mod(nper,2) == 0
    delta(2:end-1) = 2;
else
    delta(2:end) = 2;
end

Fi = nan(nplist,nplist,nfreq,nalt);
F = nan(nplist,nplist,nalt);

% Create a command-window progress bar.
if opt.progress
    progress = progressbar('IRIS model.fisher progress');
end

stationary = true(1,nalt);
for ialt = 1 : nalt    
    % Fetch the i-th parameterisation.
    m = THIS(ialt);
    
    % Minimum necessary state space.
    [T0,R0,Z0,H0,Omg0,nunit0,flag] = dogetsspace();
    if ~flag
        % Some measurement variables are non-stationary.
        stationary(ialt) = false;
        continue
    end
    
    % SGF and inverse SGF at p0.
    [G,Gi] = xxsgfy(T0,R0,Z0,H0,Omg0,nunit0,freq,opt);
    
    % Compute derivatives of SGF and steady state
    % wrt the selected parameters.
    dG = nan(ny,ny,nfreq,nplist);
    if ~opt.deviation
        dy = zeros(ny,nplist);
    end
    % Determine differentiation step.
    p0 = nan(1,nplist);
    p0(~assignNan) = m.Assign(1,assignPos(~assignNan));
    p0(~stdcorrNan) = m.stdcorr(1,stdcorrPos(~stdcorrNan));
    step = max([abs(p0);ones(1,nplist)],[],1)*eps()^opt.epspower;
    pp = p0 + step;
    pm = p0 - step;
    twosteps = pp - pm;
    for i = 1 : nplist
        dosstate = ~opt.deviation && ~isnan(assignPos(i));
        % Steady state, state space and SGF at p0(i) + step(i).
        [m,npath,sstatesuccess] = ...
            myupdatemodel(m,pp(i),assignPos(i),stdcorrPos(i),opt);
        % Force solve only in the first cycle.
        if npath ~= 1 || ~sstatesuccess
            model.failed(m,npath,sstatesuccess,'fisher');
        end
        if dosstate
            yp = dogetsstate();
        end
        [Tp,Rp,Zp,Hp,Omgp,nunitp] = dogetsspace();
        Gp = xxsgfy(Tp,Rp,Zp,Hp,Omgp,nunitp,freq,opt);
        % Steady state,state space and SGF at p0(i) - step(i).
        [m,npath,sstatesuccess] = myupdatemodel( ...
            m,pm(i),assignPos(i),stdcorrPos(i),opt);
        if npath ~= 1 || ~sstatesuccess
            model.failed(m,npath,sstatesuccess,'fisher');
        end
        if dosstate
            ym = dogetsstate();
        end
        [Tm,Rm,Zm,Hm,Omgm,nunitm] = dogetsspace();
        Gm = xxsgfy(Tm,Rm,Zm,Hm,Omgm,nunitm,freq,opt);
        % Differentiate SGF.
        dG(:,:,:,i) = (Gp - Gm) / twosteps(i);
        % Reset model parameters to `p0`.
        m.Assign(1,assignPos(~assignNan)) = p0(1,~assignNan);
        m.stdcorr(1,stdcorrPos(~stdcorrNan)) = p0(1,~stdcorrNan);
        if dosstate
            dy(:,i) = (yp(:) - ym(:)) / twosteps(i);
        end
        
        % Update the progress bar.
        if opt.progress
            update(progress,((ialt-1)*nplist+i)/(nalt*nplist));
        end
        
    end
    
    % Compute Fisher information matrix.
    % Steady-state-independent part.
    for i = 1 : nplist
        for j = i : nplist
            fi = zeros(1,nfreq);
            for k = 1 : nfreq
                fi(k) = ...
                    trace(real(Gi(:,:,k)*dG(:,:,k,i)*Gi(:,:,k)*dG(:,:,k,j)));
            end
            if ~opt.deviation
                % Add steady-state effect to zero frequency.
                % We don't divide the effect by 2*pi because
                % we skip dividing G by 2*pi, too.
                A = dy(:,i)*dy(:,j)';
                fi(1) = fi(1) + nper*trace(Gi(:,:,1)*(A + A'));
            end
            Fi(i,j,:,ialt) = fi;
            Fi(j,i,:,ialt) = fi;
            f = delta*fi';
            F(i,j,ialt) = f;
            F(j,i,ialt) = f;
        end
    end

    if opt.percent
        P0 = diag(p0);
        F(:,:,ialt) = P0*F(:,:,ialt)*P0;
    end
    
end
% End of main loop.

Fi = Fi / 2;
F = F / 2;

if any(~stationary)
    utils.warning('model', ...
        ['Cannot evaluate spectrum generating function because some ', ...
        'measurement variables are non-stationary:%s.'], ...
        sprintf(' #%g',find(~stationary)));
end

% Clean-up.
clear('model/mysstatenonlin');

% Nested functions.

%**************************************************************************
    function [T,R,Z,H,Omg,nunit,flag] = dogetsspace()
        T = m.solution{1};
        [nx,nb] = size(T);
        nf = nx - nb;
        nunit = sum(abs(abs(m.eigval(1,1:nb))-1) <= realsmall);
        Z = m.solution{4}(~exclude,:);
        % Check for (non)stationarity of observables
        stationaryindex = all(abs(Z(:,1:nunit)) <= realsmall,1);
        flag = all(stationaryindex);
        %{
        if ~flag
            T = [];
            R = [];
            Z = [];
            H = [];
            Omg = [];
            nunit = [];
            return
        end
        %}
        T = T(nf+nunit+1:end,nunit+1:end);
        Z = Z(:,nunit+1:end);
        % Cut off forward expansion.
        ne = sum(m.nametype == 3);
        R = m.solution{2}(nf+nunit+1:end,1:ne);
        H = m.solution{5}(~exclude,1:ne);
        Omg = omega(m);
    end
% dogetsspace().

%**************************************************************************
    function y = dogetsstate()
        % Get the steady-state levels for the measurement variables.
        y = m.Assign(THIS.nametype == 1);
        y = real(y);
        % Adjust for the excluded measurement variables.
        y(exclude) = [];
        % Take log of log-variables; `ylog` has been already adjusted for the
        % excluded measurement variables.
        y(ylog) = log(y(ylog));
    end
% dogetsstate().

end

% Subfunctions.

%**************************************************************************
function [G,Gi] = xxsgfy(T,R,Z,H,Omg,nunit,freq,opt) %#ok<INUSL>
% Spectrum generating function and its inverse.
% Computationally optimised for observables.
[ny,nb] = size(Z);
nfreq = length(freq(:));
Sgm1 = R*Omg*R';
Sgm2 = H*Omg*H';
G = nan(ny,ny,nfreq);
for i = 1 : nfreq
    X = Z/(eye(nb) - T*exp(-1i*freq(i)));
    G(:,:,i) = xxsymmetric(X*Sgm1*X' + Sgm2);
end
% Do not divide G by 2*pi.
% First, this cancels out in Gi*dG*Gi*dG
% and second, we do not divide the steady-state effect
% by 2*pi either.
if nargout > 1
    Gi = nan(ny,ny,nfreq);
    if opt.chksgf
        for i = 1 : nfreq
            Gi(:,:,i) = xxpinverse(G(:,:,i),opt.tolerance);
        end
    else
        for i = 1 : nfreq
            Gi(:,:,i) = inv(G(:,:,i));
        end
    end
end
end
% xxsgfy().

%**************************************************************************
function x = xxsymmetric(x)
% Minimise numerical inaccuracy between upper and lower parts
% of symmetric matrices.
index = eye(size(x)) == 1;
x = (x + x')/2;
x(index) = real(x(index));
end
% xxsymmetric().

%**************************************************************************
function X = xxpinverse(A,tol)
if isempty(A)
    X = zeros(size(A'),class(A));
    return
end
m = size(A,1);
s = svd(A);
r = sum(s/s(1) > tol);
if r == 0
    X = zeros(size(A'),class(A));
elseif r == m
    X = inv(A);
else
    [U,~,V] = svd(A,0);
    S = diag(1./s(1:r));
    X = V(:,1:r)*S*U(:,1:r)';
end
end
% xxpinverse().