function [NEWEQTN,NEWEQTNF,NEWNONLIN] = myoptpolicy(THIS,LOSSPOS,LOSSDISC)
% MYOPTPOLICY  [Not a public function] Calcualte equations for discretionary optimal policy.
%
% Backend IRIS function.
% No help provided.

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

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

% Make the model names visible inside dynamic regexps.
name = THIS.name; %#ok<NASGU>

eqtn = cell(size(THIS.eqtnF));
eqtn(:) = {''};
eqtn(THIS.nametype == 2) = THIS.eqtnF(THIS.nametype == 2);

% Replace x(:,n,t+k) with xN, xNpK, or xNmK, and &x(n) with Ln.
eqtn = sydney.myeqtn2symb(eqtn);
LOSSDISC = sydney.myeqtn2symb(LOSSDISC);

% First transition equation.
first = find(THIS.eqtntype == 2,1);

% The loss function is always the last equation. After the loss function,
% there are empty place holders. The new equations will be put in place of
% the loss function and the placeholders.
NEWEQTN = cell(size(eqtn));
NEWEQTN(:) = {''};
NEWEQTNF = NEWEQTN;
NEWNONLIN = false(size(eqtn));
zdisc = sydney(LOSSDISC);

% The lagrangian is
%
%     Mu_Eq1*eqtn1 + Mu_Eq2*eqtn2 + ... + lossfunc.
%
% The new k-th (k runs only for the original transition variables) equation
% is the derivative of the lagrangian wrt to the k-th variable, and is
% given by
%
%     Mu_Eq1*diff(eqtn1,namek) + Mu_Eq2*diff(eqtn2,namek) + ...
%     + diff(lossfunc,namek) = 0.
%
% We loop over the equations, and build gradually the new equations from
% the individual derivatives.

for eq = first : LOSSPOS
    
    % Get the list of all variables and shocks (current dates, lags, leads) in
    % this equation.
    [tmocc,nmocc] = myfindoccur(THIS,eq,'variables_shocks');
    tmocc = tmocc(:).';
    nmocc = nmocc(:).';
    
    % This is a discretionary policy. We only differentiate wrt to current
    % dates or lags of transition variables. Remove leads from the list of
    % variables wrt which we will differentiate.
    index = THIS.nametype(nmocc) == 2 & tmocc <= THIS.tzero;
    tmocc = tmocc(index);
    nmocc = nmocc(index);
    nocc = length(tmocc);
    
    % Write a cellstr with the symbolic names of variables wrt which we will
    % differentiate.
    unknown = cell(1,nocc);
    for j = 1 : nocc
        if tmocc(j) == THIS.tzero
            % Time index == 0: replace x(1,23,t) with x23.
            unknown{j} = sprintf('x%g',nmocc(j));
        else
            % Time index < 0: replace x(1,23,t-1) with x23m1.
            unknown{j} = sprintf('x%gm%g', ...
                nmocc(j),round(THIS.tzero-tmocc(j)));
        end
    end
    
    z = sydney(eqtn{eq});
    
    % Differentiate this equation wrt to all variables at once in an Inf
    % mode, that means, return a cell array of individual sydney objects.
    dz = diff(z,unknown,Inf);

    for j = 1 : nocc
        
        shift = -(tmocc(j) - THIS.tzero);
        neweq = nmocc(j);
        
        % Multiply derivatives wrt lagged variables by the discount factor.
        if shift == 1
            dz{j} = zdisc * dz{j};
        elseif shift > 1
            dz{j} = power(zdisc,shift) * dz{j};
        end
        
        % If this is not the loss function, multiply the derivative by
        % the multiplier.
        if eq < LOSSPOS
            if shift == 0
                dz{j} = sydney(sprintf('x%g',eq)) * dz{j};
            elseif shift == 1
                dz{j} = sydney(sprintf('x%gp1',eq)) * dz{j};
            else
                dz{j} = sydney(sprintf('x%gp%g',eq,shift)) * dz{j};
            end
        end
        
        deqtn = char(dz{j},'human');
        
        deqtn = regexprep(deqtn,'x(\d+)p(\d+)', ...
            '${[name{sscanf($1,''%g'')},''{+'',$2,''}'']}');
        deqtn = regexprep(deqtn,'x(\d+)m(\d+)', ...
            '${[name{sscanf($1,''%g'')},''{-'',$2,''}'']}');
        deqtn = regexprep(deqtn,'x(\d+)', ...
            '${name{sscanf($1,''%g'')}}');
        deqtn = regexprep(deqtn,'L(\d+)','${[''&'',name{sscanf($1,''%g'')}]}');
        
        % Put together the derivative of the Lagrangian wrt to variable
        % #neweq.
        if isempty(NEWEQTN{neweq})
            NEWEQTN{neweq} = [deqtn,'=0;'];
            NEWEQTNF{neweq} = [deqtn,';'];
        else
            NEWEQTN{neweq} = [deqtn,'+',NEWEQTN{neweq}];
            NEWEQTNF{neweq} = [deqtn,'+',NEWEQTNF{neweq}];
        end
        
        % Earmark the derivative for non-linear simulation if at least one equation
        % in it is non-linear and the derivative is non-zero. The derivative of the
        % loss function is supposed to be treated as non-linear if the loss
        % function itself has been introduced by min#() and not min().
        thisnonlin = THIS.nonlin(eq) && ~isequal(deqtn,'0');
        NEWNONLIN(neweq) = NEWNONLIN(neweq) || thisnonlin;

    end
end

% Replace = with #= in non-linear human equations.
NEWEQTN(NEWNONLIN) = strrep(NEWEQTN(NEWNONLIN),'=0;','=#0;');

end