function [OBJ,V,DELTA,Y] = myfdlik(THIS,DATA,Y,OPT)
% MYFDLIK  [Not a public function] Approximate likelihood function in frequency domain.
%
% Backed IRIS function.
% No help provided.

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

% TODO: Allow for non-stationary measurement variables.

NARGOUT = nargout;
s = struct();
s.noutoflik = length(OPT.outoflik);

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

nalt = size(THIS.Assign,3);
ne = sum(THIS.nametype == 3);

% Number of original periods.
[~,nper,ndata] = size(DATA);
freq = 2*pi*(0 : nper-1)/nper;

% Number of fundemantal frequencies.
N = 1 + floor(nper/2);
freq = freq(1:N);

% Band of frequencies.
frqlo = 2*pi/max(OPT.band);
frqhi = 2*pi/min(OPT.band);
ixfrq = freq >= frqlo & freq <= frqhi;

% Drop zero frequency unless requested.
if ~OPT.zero
    ixfrq(freq == 0) = false;
end
ixfrq = find(ixfrq);

% Kronecker delta.
delta = ones(1,N);
if mod(nper,2) == 0
    delta(2:end-1) = 2;
else
    delta(2:end) = 2;
end

nloop = max(nalt,ndata);
OBJ = nan(1,nloop);
V = nan(1,nloop);

% Pre-allocate output data.
DELTA = zeros(1,s.noutoflik,nloop);
dorequestoutp();
if NARGOUT >= 4
    dopreallocoutp();
end

for iloop = 1 : nloop
    
    if iloop <= ndata
        y = DATA(:,:,iloop);
        excl = OPT.exclude(:) | any(isnan(y),2);
        nyincl = sum(~excl);
        diagix = logical(eye(nyincl));
    end
    
    if iloop <= nalt

        [T,R,K,Z,H,D,U,Omg] = mysspace(THIS,iloop,false); %#ok<ASGLU>
        [nx,nb] = size(T);
        nf = nx - nb;
        realsmall = getrealsmall();
        nunit = sum(abs(abs(THIS.eigval - 1)) <= realsmall);
        % Z(1:nunit,:) assumed to be zeros.
        if any(any(abs(Z(:,1:nunit)) > realsmall))
            utils.error('model', ...
                ['Cannot evalutate likelihood in frequency domain ', ...
                'with non-stationary measurement variables.']);
        end
        T = T(nf+nunit+1:end,nunit+1:end);
        R = R(nf+nunit+1:end,1:ne);
        Z = Z(~excl,nunit+1:end);
        H = H(~excl,:);
        Sa = R*Omg*transpose(R);
        Sy = H(~excl,:)*Omg*H(~excl,:).';
        
        % Fourier transform of steady state.
        issstate = false;
        if ~OPT.deviation
            S = mytrendarray(THIS,THIS.solutionid{1},1:nper,false,iloop);
            issstate = any(S(:) ~= 0);
            if issstate
                S = fft(S.').';
            end
        end
        
        % Fourier transform of deterministic trends.
        isdtrends = false;
        if OPT.dtrends
            [D,M] = mydtrends(THIS,OPT.ttrend,OPT.outoflik,iloop);
            isdtrends = any(D(:) ~= 0);
            if isdtrends
                D = fft(D.').';
            end
            isoutoflik = ~isempty(M) && any(M(:) ~= 0);
            if isoutoflik
                M = permute(M,[3,1,2]);
                M = fft(M);
                M = ipermute(M,[3,1,2]);
            end
        end
        
    end

    yy = y;
    
    % Subtract sstate trends from observations; note that fft(y-s)
    % equals fft(y) - fft(s).
    if ~OPT.deviation && issstate
        yy = yy - S;
    end
    
    % Subtract deterministic trends from observations.
    if OPT.dtrends && isdtrends
        yy = yy - D;
    end

    yy = yy(~excl,:);
    yy = yy / sqrt(nper);
    M = M(~excl,:,:);
    M = M / sqrt(nper);
    
    L0 = 0;
    L1 = 0;
    L2 = 0;
    L3 = 0;
    nobs = 0;
    
    for i = ixfrq
        dofrequency();
    end

    [OBJ(iloop),V(iloop),DELTA(1,:,iloop),Pdelta] = ...
        model.myoutoflik(L0,L1,L2,L3,nobs,OPT);
  
    if s.retpdelta
        Y.Pdelta(:,:,iloop) = Pdelta;
    end
    
end

% Nested functions.

%**************************************************************************
    function dorequestoutp()
        s.retpdelta = isfield(Y,'Pdelta');
    end
% dorequestoutp().

%**************************************************************************
    function dopreallocoutp()
        Y.Pdelta = zeros(s.noutoflik,s.noutoflik,nloop);
    end
% dopreallocoutp().

%**************************************************************************
    function dofrequency()
        nobs = nobs + delta(i)*nyincl;
        ZiW = Z / ((eye(size(T)) - T*exp(-1i*freq(i))));
        G = ZiW*Sa*ZiW' + Sy;
        G(diagix) = real(G(diagix));
        L0 = L0 + delta(i)*real(log(det(G)));
        L1 = L1 + delta(i)*real((yy(:,i)'/G)*yy(:,i));
        if isoutoflik
            MtGi = M(:,:,i)'/G;
            L2 = L2 + delta(i)*real(MtGi*M(:,:,i));
            L3 = L3 + delta(i)*real(MtGi*yy(:,i));
        end
    end
% dofrequency().

end
