function [mloglik,Score,Info,se2] = ...
    mydiffloglik(THIS,DATA,RANGE,ASSIGNPOS,STDCORRPOS,OPT,KALMANOPT)
% MYDIFFLOGLIK  [Not a public function] Gradient and hessian of log-likelihood function.
%
% Backed IRIS function.
% No help provided.

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

if ~isfield(OPT,'progress')
    OPT.progress = false;
end

if ~isfield(OPT,'percent')
    OPT.percent = false;
end

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

OPT.fastsolve = true;
XRANGE = RANGE(1)-1:RANGE(end);

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

ny = sum(THIS.nametype == 1);

np = length(ASSIGNPOS);
[ans,nper,ndata] = size(DATA); %#ok<NOANS,ASGLU>

mloglik = zeros(1,ndata);
Score = zeros(1,np,ndata);
Info = zeros(np,np,ndata);
se2 = zeros(1,ndata);

p = nan(1,np);
assignNan = isnan(ASSIGNPOS);
stdcorrNan = isnan(STDCORRPOS);
p(~assignNan) = THIS.Assign(1,ASSIGNPOS(~assignNan));
p(~stdcorrNan) = THIS.stdcorr(1,STDCORRPOS(~stdcorrNan));

epsilon = eps()^(1/3);
step = max([abs(p);ones(size(p))],[],1)*epsilon;
pp = p + step;
mp = p - step;
twosteps = pp - mp;

% Create all parameterisations.
THIS(1:2*np+1) = THIS;
for i = 1 : np
    mindex = 1 + 2*(i-1) + 1;
    [THIS(mindex),npath,sstatesuccess] = myupdatemodel( ...
        THIS(mindex),pp(i),ASSIGNPOS(i),STDCORRPOS(i),OPT);
    if npath ~= 1 || ~sstatesuccess
        model.failed(THIS,npath,sstatesuccess,'diffloglik');
    end
    mindex = 1 + 2*(i-1) + 2;
    [THIS(mindex),npath,sstatesuccess] = myupdatemodel( ...
        THIS(mindex),mp(i),ASSIGNPOS(i),STDCORRPOS(i),OPT);
    if npath ~= 1 || ~sstatesuccess
        model.failed(THIS,npath,sstatesuccess,'diffloglik');
    end
end

% Horizontal vectorisation.
vechor = @(x) x(:)';

if OPT.progress
    % Create progress bar.
    progress = progressbar('IRIS model.diffloglik progress');
end

for idata = 1 : ndata
    domainloop();
end

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

% Nested functions.

%**************************************************************************
    function domainloop()
        
        dpe = cell(1,np);
        dpe(:) = {nan(ny,nper)};
        
        Fi_pe = zeros(ny,nper);
        X = zeros(ny);
        
        Fi_dpe = cell(1,np);
        Fi_dpe(1:np) = {nan(ny,nper)};
        
        dF = cell(1,np);
        dF(:) = {nan(ny,ny,nper)};
        
        dFvec = cell(1,np);
        dFvec(:) = {[]};
        
        Fi_dF = cell(1,np);
        Fi_dF(:) = {nan(ny,ny,nper)};

        OUTP = struct('F',[],'pe',[]);
        [mloglik(idata),se2(idata),~,OUTP] = ...
            mykalman(THIS(1),DATA(:,:,idata),OUTP,KALMANOPT);
        F = OUTP.F(:,:,2:end);
        pe = OUTP.pe(:,2:end);
        Fi = F;
        for ii = 1 : size(Fi,3)
            j = ~all(isnan(Fi(:,:,ii)),1);
            Fi(j,j,ii) = inv(Fi(j,j,ii));
        end
        
        for ii = 1 : np
            pm = THIS(1+2*(ii-1)+1);
            OUTP = struct('F',[],'pe',[]);
            [~,~,~,OUTP] = mykalman(pm,DATA(:,:,idata),OUTP,KALMANOPT);
            pF =  OUTP.F(:,:,2:end);
            ppe = OUTP.pe(:,2:end);
            
            mm = THIS(1+2*(ii-1)+2);
            OUTP = struct('F',[],'pe',[]);
            [~,~,~,OUTP] = mykalman(mm,DATA(:,:,idata),OUTP,KALMANOPT);
            mF =  OUTP.F(:,:,2:end);
            mpe = OUTP.pe(:,2:end);
            
            dF{ii}(:,:,:) = (pF - mF) / twosteps(ii);
            dpe{ii}(:,:) = (ppe - mpe) / twosteps(ii);
        end
        
        for t = 1 : nper
            o = ~isnan(pe(:,t));
            for ii = 1 : np
                Fi_dF{ii}(o,o,t) = Fi(o,o,t)*dF{ii}(o,o,t);
            end
        end
        
        for t = 1 : nper
            o = ~isnan(pe(:,t));
            for ii = 1 : np
                temp = dF{ii}(o,o,t);
                dFvec{t}(:,ii) = temp(:);
                for jj = 1 : ii
                    % Info(i,j,idata) =  ...
                    %     Info(i,j,idata) ...
                    %     + 0.5*trace(Fi_dF{i}(o,o,t)*Fi_dF{j}(o,o,t)) ...
                    %     + (transpose(dpe{i}(o,t))*Fi_dpe{j}(o,t));
                    % * the first term is data independent
                    % * trace A*B = vechor(A')*vec(B)
                    Xi = transpose(Fi_dF{ii}(o,o,t));
                    Xi = transpose(Xi(:));
                    Xj = Fi_dF{jj}(o,o,t);
                    Xj = Xj(:);
                    Info(ii,jj,idata) = Info(ii,jj,idata) + Xi*Xj/2;
                end
            end
        end
        
        % Score vector.
        for t = 1 : nper
            o = ~isnan(pe(:,t));
            Fi_pe(o,t) = Fi(o,o,t)*pe(o,t);
            X(o,o,t) = eye(sum(o)) - Fi_pe(o,t)*transpose(pe(o,t));
            dpevec = [];
            for ii = 1 : np
                dpevec = [dpevec,dpe{ii}(o,t)]; %#ok<AGROW>
                Fi_dpe{ii}(o,t) = Fi(o,o,t)*dpe{ii}(o,t);
            end
            Score(1,:,idata) = Score(1,:,idata) ...
                + vechor(Fi(o,o,t)*transpose(X(o,o,t)))*dFvec{t}/2 ...
                + transpose(Fi_pe(o,t))*dpevec;
        end
        
        % Information matrix.
        for t = 1 : nper
            o = ~isnan(pe(:,t));
            for ii = 1 : np
                for jj = 1 : ii
                    % Info(i,j,idata) =
                    %     Info(i,j,idata)
                    %     + 0.5*trace(Fi_dF{i}(o,o,t)*Fi_dF{j}(o,o,t))
                    %     + (transpose(dpe{i}(o,t))*Fi_dpe{j}(o,t));
                    % first term is data-independent and has been pre-computed.
                    Info(ii,jj,idata) = Info(ii,jj,idata) ...
                        + (transpose(dpe{ii}(o,t))*Fi_dpe{jj}(o,t));
                end
            end
        end
        
        Info(:,:,idata) = Info(:,:,idata) + transpose(tril(Info(:,:,idata),-1));
        
        if OPT.progress
            % Update progress bar.
            update(progress,idata/ndata);
        end
        
    end

end