function [THIS,DERIV,NANDERIV] = myderiv(THIS,EQSELECT,IALT,SYMBOLIC)
% MYDERIV  [Not a public function] Compute first-order expansion of equations around current steady state.
%
% Backed IRIS function.
% No help provided.

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

isnanderiv = nargout > 2;

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

% Copy last computed derivatives.
DERIV = struct();
DERIV.c = THIS.deriv0.c;
DERIV.f = THIS.deriv0.f;
DERIV.n = THIS.deriv0.n;

assign = THIS.Assign(1,:,IALT);
nname = length(THIS.name);
neqtn = length(THIS.eqtn);
EQSELECT(THIS.eqtntype >= 3) = false;

NANDERIV = false(1,neqtn);

% Prepare 3D occur array limited to occurences of variables and shocks in
% measurement and transition equations.
occur = full(THIS.occur);
occur = reshape(occur,[neqtn,nname,size(THIS.occur,2)/nname]);
occur = occur(THIS.eqtntype <= 2,THIS.nametype <= 3,:);
occur = permute(occur,[3,2,1]);

if any(EQSELECT)
    nt = size(THIS.occur,2) / nname;
    nvar = sum(THIS.nametype <= 3);
    t = THIS.tzero;
    issymb = ~cellfun(@isempty,THIS.deqtnF);
    if ~SYMBOLIC
        issymb(:) = false;
    end
    symbselect = issymb & EQSELECT;
    numselect = ~issymb & EQSELECT;
    if any(symbselect)
        % Symbolic derivatives.
        dosymbderiv();
    end
    if any(numselect)
        % Numerical derivatives.
        donumderiv();
    end
    % Reset the add-factors in non-linear equations to 1.
    tempeye = -eye(sum(THIS.eqtntype <= 2));
    DERIV.n(EQSELECT,:) = tempeye(EQSELECT,THIS.nonlin);
    % Normalise derivatives by largest number in non-linear models.
    if ~THIS.linear
        for ieq = find(EQSELECT)
            index = DERIV.f(ieq,:) ~= 0;
            if any(index)
                norm = max(abs(DERIV.f(ieq,index)));
                DERIV.f(ieq,index) = DERIV.f(ieq,index) / norm;
                DERIV.n(ieq,:) = DERIV.n(ieq,:) / norm;
            end
        end
    end
end

if IALT == 1
    THIS.Assign0(:) = THIS.Assign(1,:,IALT);
    THIS.deriv0.c(:) = DERIV.c;
    THIS.deriv0.f(:) = DERIV.f;
    THIS.deriv0.n(:) = DERIV.n;
end

% Nested functions.

%**************************************************************************
    function donumderiv()
        
        mint = 1 - t;
        maxt = nt - t;
        tvec = mint : maxt;
        
        if THIS.linear
            init = zeros(1,nname);
            init(1,THIS.nametype == 4) = real(assign(THIS.nametype == 4));
            init = init(1,:,ones(1,nt));
            h = ones(1,nname,nt);
        else
            init = mytrendarray(THIS,1:nname,tvec,false,IALT);
            init = shiftdim(init,-1);
            h = abs(THIS.epsilon(1))*max([init;ones([1,nname,nt])],[],1);
        end
        
        plus = init + h;
        minus = init - h;
        step = plus - minus;
        
        if any(THIS.log)
            init(1,THIS.log,:) = exp(init(1,THIS.log,:));
            plus(1,THIS.log,:) = exp(plus(1,THIS.log,:));
            minus(1,THIS.log,:) = exp(minus(1,THIS.log,:));
        end
        
        % References to steady-state levels and growth rates.
        if ~THIS.linear
            L = init(:,:,t);
        else
            L = [];
        end
        
        fgrid = cell(1,3);
        fevaluate = cell(1,3);
        
        for iieq = find(numselect)
            % Get occurences of variables in this equation.
            [tmocc,nmocc] = find(occur(:,:,iieq));
            
            % Total number of derivatives to be computed in this equation.
            n = length(nmocc);
            fgrid{2} = init;
            for ii = [1,3]
                fgrid{ii} = init(ones(1,n),:,:);
            end
            for ii = 1 : n
                fgrid{1}(ii,nmocc(ii),tmocc(ii)) = ...
                    minus(1,nmocc(ii),tmocc(ii));
                fgrid{3}(ii,nmocc(ii),tmocc(ii)) = ...
                    plus(1,nmocc(ii),tmocc(ii));
            end
            
            for ii = [1,3]
                x = fgrid{ii};
                fevaluate{ii} = THIS.eqtnF{iieq}(x,t,L);
            end
            
            % Constant in linear models.
            if THIS.linear
                x = fgrid{2};
                DERIV.c(iieq) = THIS.eqtnF{iieq}(x,t,L);
            end
            
            value = zeros(1,n);
            for ii = 1 : n
                value(ii) = (fevaluate{3}(ii)-fevaluate{1}(ii)) ...
                    / step(1,nmocc(ii),tmocc(ii));
            end
            
            % Round numbers close to integers.
            % roundIndex = abs(value - round(value)) <= realsmall;
            % value(roundIndex) = round(value(roundIndex));
            
            % Assign values to the array of derivatives.
            index = (tmocc-1)*nvar + nmocc;
            DERIV.f(iieq,index) = value;
            
            % Check for NaN derivatives.
            if isnanderiv && any(~isfinite(value))
                NANDERIV(iieq) = true;
            end

        end
        
    end
% donumderiv().

%**************************************************************************
    function dosymbderiv()

        if THIS.linear
            x = zeros(1,nname);
            if any(THIS.log)
                x(1,THIS.log) = 1;
            end
            x(1,THIS.nametype == 4) = real(assign(THIS.nametype == 4));
            x = x(1,:,ones(1,nt));
            % References to steady-state levels and growth rates.
            L = [];
        else
            mint = 1 - t;
            maxt = nt - t;
            tvec = mint : maxt;
            x = mytrendarray(THIS,1:nname,tvec,true,IALT);
            x = shiftdim(x,-1);
            % References to steady-state levels and growth rates.
            L = x;
        end
   
        for iieq = find(symbselect)
            % Get occurences of variables in this equation.
            [tmocc,nmocc] = find(occur(:,:,iieq));
            
            % Constant in linear models. Becuase all variables are set to
            % zero, evaluating the equations gives the constant.
            if THIS.linear
                if isnumeric(THIS.ceqtnF{iieq})
                    c = THIS.ceqtnF{iieq};
                else
                    c = THIS.ceqtnF{iieq}(x,t,L);
                end
                DERIV.c(iieq) = c;
            end
            % Evaluate all derivatives of the equation at once.
            value = THIS.deqtnF{iieq}(x,t,L);
            
            % Assign values to the array of derivatives.
            index = (tmocc-1)*nvar + nmocc;
            DERIV.f(iieq,index) = value;
            
            % Check for NaN derivatives.
            if isnanderiv && any(~isfinite(value))
                NANDERIV(iieq) = true;
            end

        end
        
    end
% dosymbderiv().

end
