function [THIS,assign] = myparse(THIS,P,OPT)
% parse  [Not a public function] Parse model code.
%
% Backend IRIS function.
% No help provided.

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

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

assign = P.assign;
THIS.fname = P.fname;

% LINEAR OR NON-LINEAR MODEL
%============================

% Linear or non-linear model. First, check for the presence of th keyword
% `!linear` in the model code. However, if the user have specified the
% `'linear='` option in the `model` function, use that.
[THIS.linear,P.code] = strfun.findremove(P.code,'!linear');
if ~isempty(OPT.linear)
    THIS.linear = OPT.linear;
end

% RUN THE THETA PARSER
%======================

% Run theparser on the model file.
q = theparser('model',P);
S = parse(q,'error');

% LIST OF FUNCTIONS WITH USER DERIVATIVES
%=========================================

THIS.userdifflist = regexp(S(14).blk,'[a-zA-Z][\w\.]+','match');

% VARIABLES, SHOCKS AND PARAMETERS
%==================================

% Read the individual names of variables, shocks, and parameters.
blkorder = [1,2,11,12,4];
name = [S(blkorder).name];
nametype = [S(blkorder).nametype];
namelabel = [S(blkorder).namelabel];
namevalue = [S(blkorder).namevalue];
shocktype = nan(size(nametype));
shocktype(nametype == 11) = 1;
shocktype(nametype == 12) = 2;
nametype(nametype == 11 | nametype == 12) = 3;

% Check the following naming rules:
%
% * Names must not start with 0-9 or _.
% * The name `ttrend` is a reserved name for time trend in `!dtrends`.
% * Shock names must not contain double scores because of the way
% cross-correlations are referenced.
%
invalid = ~cellfun(@isempty,regexp(name,'^[0-9_]','once')) ...
    | strcmp(name,'ttrend') ...
    | (~cellfun(@isempty,strfind(name,'__')) & nametype == 3);
if any(invalid)
    % Invalid variable or parameter names.
    utils.error('model',[utils.errorparsing(THIS), ....
        'This is not a valid variable, shock, or parameter name: ''%s''.'], ...
        name{invalid});
end

% Evaluate values assigned in the model code and/or in the `assign`
% database. Evaluate parameters first so that they are available for
% steady-state expressions.
namevalue = strtrim(namevalue);
namevalue = regexprep(namevalue,'(\<[A-Za-z]\w*\>)(?![\(\.])', ...
    '${utils.iff(any(strcmpi($1,{''Inf'',''NaN''})),$1,[''assign.'',$1])}');
if isstruct(assign) && ~isempty(assign) 
    donoteval = fieldnames(assign);
else
    donoteval = {};
end
for itype = 4 : -1 : 1
    % Assign a value from declaration only if not in the input database.
    for j = find(nametype == itype)
        if isempty(namevalue{j}) || any(strcmp(name{j},donoteval))
            continue
        end
        try
            temp = eval(namevalue{j});
            if isnumeric(temp) && length(temp) == 1
                assign.(name{j}) = temp;
            end
        catch %#ok<CTCH>
            assign.(name{j}) = NaN;
        end
    end
end

% Remove user-declared std_ and corr_ from the list of params.
% Find all names starting with std_ or corr_.
stdindex = strncmp(name,'std_',4);
corrindex = strncmp(name,'corr_',5);

% Allowed are only parameter names starting with std_ or corr_, not
% variable or shock names.
invalid = (stdindex | corrindex) & nametype ~= 4;
if any(invalid)
    utils.error('model',[utils.errorparsing(THIS), ...
        'This is not a valid variable or shock name: ''%s''.'], ...
        name{invalid});
end

% Remove the declared std_ and corr_ names from the list of names.
if any(stdindex) || any(corrindex)
    stdname = name(stdindex);
    corrname = name(corrindex);
    name(stdindex | corrindex) = [];
    namelabel(stdindex | corrindex) = [];
    nametype(stdindex | corrindex) = [];
end

% Check for multiple names unless `'multiple=' true`.
if ~OPT.multiple
    nonunique = strfun.nonunique(name);
    if ~isempty(nonunique)
        utils.error('model',[utils.errorparsing(THIS), ...
            'This name is declared more than once: ''%s''.'], ...
            nonunique{:});
    end
else
    % Take the last defined/assigned unique name.
    [name,index] = unique(name,'last');
    nametype = nametype(index);
    shocktype = shocktype(index);
    namelabel = namelabel(index);
end

% Sort variable, shock and parameter names by the nametype.
[THIS.nametype,index] = sort(nametype);
THIS.name = name(index);
THIS.namelabel = namelabel(index);
shocktype = shocktype(index);
shocktype = shocktype(THIS.nametype == 3);

% Check that std and corr names refer to valid shock names.
dochkstdcorrnames();

% LOG VARIABLES
%===============

THIS.log = false(size(THIS.name));
THIS.log(THIS.nametype == 1) = S(1).nameflag;
THIS.log(THIS.nametype == 2) = S(2).nameflag;

% REPORTING EQUATIONS
%=====================

% TODO: Use theparser object instead of preparser object.
p1 = P;
p1.code = S(10).blk;
THIS.outside = reporting(p1);

% READ INDIVIDUAL EQUATIONS
%===========================

% There are four types of equations: measurement equations, transition
% equations, deterministic trends, and dynamic links.

% Read measurement equations.
[eqtn,eqtnF,eqtnS,eqtnlabel] = xxreadeqtns(S(6));
n = length(eqtn);
THIS.eqtn(end+(1:n)) = eqtn;
THIS.eqtnF(end+(1:n)) = eqtnF;
if ~THIS.linear
    THIS.eqtnS(end+(1:n)) = eqtnS;
else
    THIS.eqtnS(end+(1:n)) = {''};
end
THIS.eqtnlabel(end+(1:n)) = eqtnlabel;
THIS.eqtntype(end+(1:n)) = 1;
THIS.nonlin(end+(1:n)) = false;

% Read transition equations; loss function is always moved to the end.

[eqtn,eqtnF,eqtnS,eqtnlabel,nonlin,lossdisc,multiple] = xxreadeqtns(S(7));

if ischar(lossdisc) && isempty(lossdisc)
    utils.error('model',[utils.errorparsing(THIS), ...
        'Loss function discount factor is empty.']);
end

if multiple
    utils.error('model',[utils.errorparsing(THIS), ...
        'Multiple loss functions found in the transition equations.']);
end

n = length(eqtn);
THIS.eqtn(end+(1:n)) = eqtn;
THIS.eqtnF(end+(1:n)) = eqtnF;
if ~THIS.linear
    THIS.eqtnS(end+(1:n)) = eqtnS;
else
    THIS.eqtnS(end+(1:n)) = {''};
end
THIS.eqtnlabel(end+(1:n)) = eqtnlabel;
THIS.eqtntype(end+(1:n)) = 2;
THIS.nonlin(end+(1:n)) = nonlin;

% Check for empty dynamic equations. This may occur if the user types a
% semicolon between the full equations and its steady state version.
dochkemptyeqtn();

THIS.multiplier = false(size(THIS.name));
isloss = ischar(lossdisc) && ~isempty(lossdisc);
if isloss
    % Create placeholders for new transition names (mutlipliers) and new
    % transition equations (derivatives of the loss function wrt existing
    % variables).
    dolossfuncplaceholders();
end

% Introduce `nname` and `neqtn` only after we are done with placeholders
% for optimal policy variables and equations.
nname = length(THIS.name);
neqtn = length(THIS.eqtn);

% Read deterministic trend equaitons.

[THIS,logmissing,invalid,multiple] = xxreaddtrends(THIS,S(9));

if ~isempty(logmissing)
    utils.error('model',[utils.errorparsing(THIS), ...
        'The LHS variable must be logarithmised in this dtrend equation: ''%s''.'], ...
        logmissing{:});
end

if ~isempty(invalid)
    utils.error('model',[utils.errorparsing(THIS), ...
        'Invalid LHS in this dtrend equation: ''%s''.'], ...
        invalid{:});
end

if ~isempty(multiple)
    utils.error('model',[utils.errorparsing(THIS), ...
        'Mutliple dtrend equations ', ...
        'for this measurement variable: ''%s''.'], ...
        multiple{:});
end

% Read dynamic links.

[THIS,invalid] = xxreadlinks(THIS,S(13));

if ~isempty(invalid)
    utils.error('model',[utils.errorparsing(THIS), ...
        'Invalid LHS in this dynamic link: ''%s''.'], ...
        invalid{:});
end

% PROCESS EQUATIONS
%===================

% Remove label marks #(xx) from equations.
THIS.eqtn = regexprep(THIS.eqtn,'^\s*#\(\d+\)\s*','');

% Remove ! from math functions.
% This is for bkw compatibility only.
THIS.eqtnF = strrep(THIS.eqtnF,'!','');
if ~THIS.linear
    THIS.eqtnS = strrep(THIS.eqtnS,'!','');
end

% Remove blank spaces.
THIS.eqtn = regexprep(THIS.eqtn,{'\s+','".*?"'},{'',''});
THIS.eqtnF = regexprep(THIS.eqtnF,'\s+','');
if ~THIS.linear
    THIS.eqtnS = regexprep(THIS.eqtnS,'\s+','');
end
THIS.eqtnlabel = strrep(THIS.eqtnlabel,'"','');

% Replace names with code characters.
nameoffset = 1999;
namecode = char(nameoffset + (1:nname));

% Prepare patterns and their code substitutions for all variables,
% shocks, and parameter names.
codepatt = cell(1,nname);
coderepl = cell(1,nname);
len = cellfun(@length,THIS.name);
[aux,index] = sort(len,2,'descend');
for i = index
    codepatt{i} = ['\<',THIS.name{i},'\>'];
    coderepl{i} = namecode(i);
end
THIS.eqtnF = regexprep(THIS.eqtnF,codepatt,coderepl);
if ~THIS.linear
    THIS.eqtnS = regexprep(THIS.eqtnS,codepatt,coderepl);
end

% Try to catch undeclared names in all equations except dynamic links at
% this point; all valid names have been substituted for by the name codes.
% Do not do it in dynamic links because the links can contain std and corr
% names which have not been substituted for.
dochkundeclared();

% Check for sstate references occuring in wrong places. Also replace
% the old syntax & with $.
dochksstateref();

% MAX LAG AND LEAD
%==================

maxt = max([S.maxt]);
mint = min([S.mint]);
if isloss
    % Anticipate that multipliers will have leads as far as the greatest lag.
    maxt = max([maxt,-mint]);
end
maxt = maxt + 1;
mint = mint - 1;
tzero = 1 - mint;
THIS.tzero = tzero;
nt = maxt - mint + 1;

% REPLACE NAME CODES WITH WITH X(...)
%=====================================

% Allocate ise the `occur` and `occurS` properties. These need to be
% allocated before we run `xxtransformeqtn` for the first time.
THIS.occur = sparse(false(length(THIS.eqtnF),length(THIS.name)*nt));
THIS.occurS = sparse(false(length(THIS.eqtnF),length(THIS.name)));

THIS = xxtransformeqtn(THIS,[],namecode,Inf);

% Check equation syntax before we compute optimal policy.
dochksyntax(Inf);

if isloss
    % Parse the discount factor first, as it can be a general expression.
    lossdisc = regexprep(lossdisc,codepatt,coderepl);
    [THIS,lossdisc] = xxtransformeqtn(THIS,lossdisc,namecode);
    
    % Create optimal policy equations by adding the derivatives
    % of the lagrangian wrt to the original transition variables. These
    % `naddeqtn` new equation will be put in place of the loss function
    % and the `naddeqtn-1` empty placeholders.
    [neweqtn,neweqtnF,newnnonlin] = myoptpolicy(THIS,losspos,lossdisc);
    
    % Add the new equations to the model object, and parse them.
    last = find(THIS.eqtntype == 2,1,'last');
    THIS.eqtn(losspos:last) = neweqtn(losspos:last);
    THIS.eqtnF(losspos:last) = neweqtnF(losspos:last);
    THIS.eqtnF(losspos:last) = ...
        regexprep(THIS.eqtnF(losspos:last),codepatt,coderepl);
    
    % Add sstate equations. Note that we must at least replace the old equation
    % in `losspos` position (which was the objective function) with the new
    % equation (which is a derivative wrt to the first variables).
    THIS.eqtnS(losspos:last) = THIS.eqtnF(losspos:last);
    
    THIS.nonlin(losspos:last) = newnnonlin(losspos:last);
    
    % Replace name codes with x(...) in the new F and S equations.
    THIS = xxtransformeqtn(THIS,[],namecode,losspos:last);
    
    % Check syntax of newly created optimal policy equations.
    dochksyntax(losspos:last);
end

% FINISHING TOUCHES
%===================

% Vectorise *, /, \, ^ operators.
THIS.eqtnF = strfun.vectorise(THIS.eqtnF);

% Check the model structure.
[errorMessage,errorList] = dochkstructure();
if ~isempty(errorMessage)
    utils.error('model',[utils.errorparsing(THIS),errorMessage],errorList{:});
end

% Create placeholders for non-linearised equations.
THIS.eqtnN = cell(size(THIS.eqtn));
THIS.eqtnN(:) = {''};

% Make sure all equations end with semicolons.
for ieq = 1 : length(THIS.eqtn)
    if ~isempty(THIS.eqtn{ieq}) && THIS.eqtn{ieq}(end) ~= ';'
        THIS.eqtn{ieq}(end+1) = ';';
    end
    if ~isempty(THIS.eqtnF{ieq}) && THIS.eqtnF{ieq}(end) ~= ';'
        THIS.eqtnF{ieq}(end+1) = ';';
    end
    if ~isempty(THIS.eqtnS{ieq}) && THIS.eqtnS{ieq}(end) ~= ';'
        THIS.eqtnS{ieq}(end+1) = ';';
    end
end

% Nested functions.

%**************************************************************************
    function dochkstdcorrnames()
        
        if ~any(stdindex) && ~any(corrindex)
            % No std or corr names declared.
            return
        end
        
        if ~isempty(stdname)
            % Check that all std names declared by the user refer to a valid shock
            % name.
            [ans,pos] = mynameposition(THIS,stdname); %#ok<NOANS,ASGLU>
            invalid = stdname(isnan(pos));
            if ~isempty(invalid)
                utils.error('model',[utils.errorparsing(THIS), ...
                    'This is not a valid std deviation name: ''%s''.'], ...
                    invalid{:});
            end
        end
        
        if ~isempty(corrname)
            % Check that all corr names declared by the user refer to valid shock
            % names.
            [ans,pos] = mynameposition(THIS,corrname); %#ok<NOANS,ASGLU>
            invalid = corrname(isnan(pos));
            if ~isempty(invalid)
                utils.error('model',[utils.errorparsing(THIS), ...
                    'This is not a valid cross-correlation name: ''%s''.'], ...
                    invalid{:});
            end
        end
        
    end
% dochkstdcorrnames().

%**************************************************************************
    function [errorMessage,errorList] = dochkstructure()
        
        occurF = reshape(full(THIS.occur),[neqtn,nname,nt]);
        
        errorMessage = '';
        errorList = {};
        
        % At least one transition variable.
        if ~any(THIS.nametype == 2)
            errorMessage = 'No transition variable.';
            return
        end
        
        % Current dates of all transition variables.
        aux = ~any(occurF(THIS.eqtntype == 2,THIS.nametype == 2,tzero),1);
        if any(aux)
            errorList = THIS.name(THIS.nametype == 2);
            errorList = errorList(aux);
            errorMessage = ...
                'No current date of this transition variable: ''%s''.';
            return
        end
        
        % Current dates of all measurement variables.
        aux = ~any(occurF(THIS.eqtntype == 2,THIS.nametype == 2,tzero),1);
        if any(aux)
            errorList = THIS.name(THIS.nametype == 1);
            errorList = errorList(aux);
            errorMessage = ...
                'No current date of this measurement variable: ''%s''.';
            return
        end
        
        % # measurement equations == # measurement variables.
        nme = sum(THIS.eqtntype == 1);
        nmv = sum(THIS.nametype == 1);
        if nme ~= nmv
            errorMessage = sprintf( ...
                '%g measurement equation(s) for %g measurement variable(s).', ...
                nme,nmv);
            return
        end
        
        % # transition equations == # transition variables.
        nte = sum(THIS.eqtntype == 2);
        ntv = sum(THIS.nametype == 2);
        if nte ~= ntv
            errorMessage = sprintf(['%g transition equation(s) ', ...
                'for %g transition variable(s).'],nte,ntv);
            return
        end
        
        % No lags/leads of measurement variables.
        tt = true([1,size(occurF,3)]);
        tt(tzero) = false;
        aux = any(any(occurF(:,THIS.nametype == 1,tt),3),1);
        if any(aux)
            errorList = THIS.name(THIS.nametype == 1);
            errorList = errorList(aux);
            errorMessage = ...
                'This measurement variable has a lag or a lead: ''%s''.';
            return
        end
        
        % No lags/leads of shocks.
        aux = any(any(occurF(:,THIS.nametype == 3,tt),3),1);
        if any(aux)
            errorList = THIS.name(THIS.nametype == 3);
            errorList = errorList(aux);
            errorMessage = 'This shock has a lag or a lead: ''%s''.';
            return
        end
        
        % No lags/leads of parameters.
        aux = any(any(occurF(:,THIS.nametype == 4,tt),3),1);
        if any(aux)
            errorList = THIS.name(THIS.nametype == 4);
            errorList = errorList(aux);
            errorMessage = 'This parameter has a lag or a lead: ''%s''.';
            return
        end
        
        % No measurement variables in transition equations.
        aux = any(any(occurF(THIS.eqtntype == 2,THIS.nametype == 1,:),3),2);
        if any(aux)
            errorList = THIS.eqtn(THIS.eqtntype == 2);
            errorList = errorList(aux);
            errorMessage = ['This transition equation refers to ', ...
                'measurement variable(s): ''%s''.'];
            return
        end
        
        % No leads of transition variables in measurement equations.
        tt = true([1,size(occurF,3)]);
        tt(1:tzero) = false;
        aux = any(any(occurF(THIS.eqtntype == 1,THIS.nametype == 2,tt),3),2);
        if any(aux)
            errorList = THIS.eqtn(THIS.eqtntype == 1);
            errorList = errorList(aux);
            errorMessage = ['Lead(s) of transition variable(s) in this ', ...
                'measurement equation: ''%s''.'];
            return
        end
        
        % Current date of any measurement variable in each measurement
        % equation.
        aux = ~any(occurF(THIS.eqtntype == 1,THIS.nametype == 1,tzero),2);
        if any(aux)
            errorList = THIS.eqtn(THIS.eqtntype == 1);
            errorList = errorList(aux);
            errorMessage = ['No current-dated measurement variables ', ...
                'in this measurement equation: ''%s''.'];
            return
        end
        
        if any(THIS.nametype == 3)
            % Find shocks in measurement equations.
            aux1 = any(occurF(THIS.eqtntype == 1,THIS.nametype == 3,tzero),1);
            % Find shocks in transition equations.
            aux2 = any(occurF(THIS.eqtntype == 2,THIS.nametype == 3,tzero),1);
            % No measurement shock in transition equations.
            aux = aux2 & shocktype == 1;
            if any(aux)
                errorList = THIS.name(THIS.nametype == 3);
                errorList = errorList(aux);
                errorMessage = ['This measurement shock occurs ', ...
                    'in transition equation(s): ''%s''.'];
                return
            end
            % No transition shock in measurement equations.
            aux = aux1 & shocktype == 2;
            if any(aux)
                errorList = THIS.name(THIS.nametype == 3);
                errorList = errorList(aux);
                errorMessage = ['This transition shock occurs ', ...
                    'in measurement equation(s): ''%s''.'];
                return
            end
        end
        
        % Only parameters can occurF in deterministic trend equations.
        aux = any(any(occurF(THIS.eqtntype == 3,THIS.nametype ~= 4,:),3),2);
        if any(aux)
            errorList = THIS.eqtn(THIS.eqtntype == 3);
            errorList = errorList(aux);
            errorMessage = ['This deterministic trend refers to name(s) ', ...
                'other than parameters: ''%s''.'];
            return
        end
        
    end
% dochkstructure().

%**************************************************************************
    function dochkundeclared()
        % Undeclared names have not been substituted for by the name codes, except
        % std and corr names in dynamic links (std and corr names cannot be used in
        % other types of equations). Undeclared names in dynamic links will be
        % caught in `dochksyntax`. Distinguish variable names from function names
        % (func names are immediately followed by an opening bracket).
        % Unfortunately, `regexp` interprets high char codes as \w, so we need to
        % explicitly type the ranges.
        
        undeclared = {};
        stdcorr = {};
        for iieq = find(THIS.eqtntype ~= 4)
            list = regexp(THIS.eqtnF{iieq}, ...
                '\<[a-zA-Z][a-zA-Z0-9_]*\>(?![\(\.])','match');
            list = setdiff(unique(list),{'ttrend'});
            if ~isempty(list)
                for ii = 1 : length(list)
                    if strncmp(list{ii},'std_',4) ...
                        || strncmp(list{ii},'corr_',5)
                        stdcorr{end+1} = list{ii}; %#ok<AGROW>
                        stdcorr{end+1} = THIS.eqtn{iieq}; %#ok<AGROW>
                    else
                        undeclared{end+1} = list{ii}; %#ok<AGROW>
                        undeclared{end+1} = THIS.eqtn{iieq}; %#ok<AGROW>
                    end
                end
            end
        end
        
        % Report std or corr names used in equations other than links.
        if ~isempty(stdcorr)
            utils.error('model',[utils.errorparsing(THIS), ...
                'Std or corr name ''%s'' cannot be used in ''%s''.'], ...
                stdcorr{:});
        end
        
        % Report non-function names that have not been declared.
        if ~isempty(undeclared)
            utils.error('model',[utils.errorparsing(THIS), ...
                'Undeclared or mistyped name ''%s'' in ''%s''.'], ...
                undeclared{:});
        end
    end
% dochkundeclared().

%**************************************************************************
    function dochksstateref()
        % Check for sstate references in wrong places.
        func = @(x)~cellfun(@isempty,strfind(x,'&'));
        index = func(THIS.eqtnF);
        % Not allowed in linear models.
        if THIS.linear
            if any(index)
                utils.error('model',[utils.errorparsing(THIS), ...
                    'Steady-state references not allowed ', ...
                    'in linear models: ''%s''.'], ...
                    THIS.eqtn{index});
            end
            return
        end
        index = index | func(THIS.eqtnS);
        % Not allowed in deterministic trends.
        temp = index & THIS.eqtntype == 3;
        if any(temp)
            utils.error('model',[utils.errorparsing(THIS), ...
                'Steady-state references not allowed ', ...
                'in dtrends equations: ''%s''.'], ...
                THIS.eqtn{temp});
        end
        % Not allowed in dynamic links.
        temp = index & THIS.eqtntype == 4;
        if any(temp)
            utils.error('model',[utils.errorparsing(THIS), ...
                'Steady-state references not allowed ', ...
                'in dynamic links: ''%s''.'], ...
                THIS.eqtn{temp});
        end
    end
% dochksstateref().

%**************************************************************************
    function dolossfuncplaceholders()
        % Add new variables, i.e. the Lagrange multipliers associated with
        % all of the existing transition equations except the loss
        % function. These new names will be ordered first -- the logic is
        % that the final equations will be ordered as derivatives of the
        % lagrangian wrt to the individual variables.
        naddeqtn = sum(THIS.nametype == 2) - 1;
        naddname = sum(THIS.eqtntype == 2) - 1;
        newname = cell(1,naddname-1);
        newnametype = 2*ones(1,naddname);
        newnamelabel = cell(1,naddname);
        newnamelabel(:) = {''};
        newlog = false(1,naddname);
        newmultiplier = true(1,naddname);
        % The default name is |'Mu_Eq%g'| but can be changed through the
        % option `'multiplierName='`.
        for ii = 1 : naddname
            newname{ii} = sprintf(OPT.multipliername,ii);
        end
        % Insert the new names between at the beginning of the blocks of existing
        % transition variables.
        THIS.name = [THIS.name(THIS.nametype < 2),newname,THIS.name(THIS.nametype >= 2)];
        THIS.namelabel = [THIS.namelabel(THIS.nametype < 2), ...
            newnamelabel,THIS.namelabel(THIS.nametype >= 2)];
        THIS.log = [THIS.log(THIS.nametype < 2),newlog,THIS.log(THIS.nametype >= 2)];
        THIS.multiplier = [THIS.multiplier(THIS.nametype < 2), ...
            newmultiplier,THIS.multiplier(THIS.nametype >= 2)];
        % nametype must come last.
        THIS.nametype = [THIS.nametype(THIS.nametype < 2), ...
            newnametype,THIS.nametype(THIS.nametype >= 2)];
        % Loss function is always ordered last among transition equations.
        losspos = length(THIS.eqtn);
        % We will add `naddeqtn` new transition equations, i.e. the
        % derivatives of the Lagrangiag wrt the existing transition
        % variables. At the same time, we will remove the loss function so
        % we need to create only `naddeqtn-1` placeholders.
        THIS.eqtn(end+(1:naddeqtn)) = {''};
        THIS.eqtnF(end+(1:naddeqtn)) = {''};
        THIS.eqtnS(end+(1:naddeqtn)) = {''};
        THIS.eqtnlabel(end+(1:naddeqtn)) = {''};
        THIS.nonlin(end+(1:naddeqtn)) = false;
        THIS.eqtntype(end+(1:naddeqtn)) = 2;
    end
% dolossfuncplaceholders().

%**************************************************************************
    function dochksyntax(eqtnlist)
        t = THIS.tzero;
        nname = length(THIS.name);
        neqtn = length(THIS.eqtn);
        ne = sum(THIS.nametype == 3);
        std = double(THIS.linear)*1 + double(~THIS.linear)*0.1;
        x = rand(1,nname,nt);
        % This is the test vector for dynamic links. In dynamic links, we allow std
        % and corr names to occurs, and append them to the assign vector.
        if any(THIS.eqtntype == 4)
            x1 = [rand(1,nname,1),std*ones(1,ne),zeros(1,ne*(ne-1)/2)];
        end
        L = x(1,:,t);
        dx = zeros(1,nname,nt);
        ttrend = 0;
        undeclared = {};
        syntax = {};
        if isequal(eqtnlist,Inf)
            eqtnlist = 1 : neqtn;
        end
        for eq = eqtnlist
            try
                eqtn = THIS.eqtnF{eq};
                if isempty(eqtn)
                    continue
                end
                eqtn = strfun.vectorise(eqtn);
                eqtn = str2func(['@(x,dx,L,t,ttrend)',eqtn]);
                if THIS.eqtntype(eq) < 4
                    eqtn(x,dx,L,t,ttrend);
                else
                    % Evaluate RHS of dynamic links. They can refer to std or corr names, so we
                    % have to use the `x1` vector.
                    eqtn(x1,[],[],1,[]);
                end
                if eq <= length(THIS.eqtnS) && THIS.eqtntype(eq) <= 2
                    eqtn = THIS.eqtnS{eq};
                    eqtn = strrep(eqtn,'exp?','exp');
                    eqtn = str2func(['@(x,dx,L,t,ttrend)',THIS.eqtnF{eq}]);
                    eqtn(x,dx,L,t,ttrend);
                end
            catch E
                % Undeclared names should have been already caught. But a few exceptions
                % may still exist.
                [match,tokens] = ...
                    regexp(E.message,'Undefined function or variable ''(\w*)''','match','tokens','once');
                if ~isempty(match)
                    undeclared{end+1} = tokens{1}; %#ok<AGROW>
                    undeclared{end+1} = THIS.eqtn{eq}; %#ok<AGROW>
                else
                    syntax{end+1} = THIS.eqtn{eq}; %#ok<AGROW>
                    if E.message ~= '.'
                        E.message(end+1) = '.';
                    end
                    syntax{end+1} = E.message; %#ok<AGROW>
                end
            end
        end
        if ~isempty(undeclared)
            utils.error('model',[utils.errorparsing(THIS), ...
                'Undeclared or mistyped name ''%s'' in ''%s''.'], ...
                undeclared{:});
        end
        if ~isempty(syntax)
            utils.error('model',[utils.errorparsing(THIS), ...
                'Syntax error in ''%s''.\n\tMatlab says: %s'], ...
                syntax{:});
        end
    end
% dochksyntax().

%**************************************************************************
    function dochkemptyeqtn()
        % dochkemptyeqtn  Check for empty full equations.
        
        emptyeqtn = cellfun(@isempty,THIS.eqtnF);
        if any(emptyeqtn)
            utils.error('model',[utils.errorparsing(THIS), ...
                'This equation is empty: ''%s''.'], ...
                THIS.eqtn{emptyeqtn});
        end
        
    end
% dochkemptyeqtn().

end

% Subfunctions.

%**************************************************************************
function ...
    [EQTN,EQTNF,EQTNS,EQTNLABEL,EQTNNONLIN,LOSSDISC,MULTIPLE] = ...
    xxreadeqtns(S)
% xxreadeqtns  Read measurement or transition equations.

if isempty(S.eqtn) 
    EQTN = cell(1,0);
    EQTNLABEL = cell(1,0);
    EQTNF = cell(1,0);
    EQTNS = cell(1,0);
    EQTNNONLIN = false(1,0);
    return
end

% Check for a loss function and its discount factor if requested by the
% caller (done for transition equations only).
if nargout >= 6
    LOSSDISC = NaN;
    MULTIPLE = false;
    dolossfunc();
end

EQTN = S.eqtn;
EQTNLABEL = S.eqtnlabel;
EQTNNONLIN = strcmp(S.eqtnsign,'=#');

neqtn = length(S.eqtn);
EQTNF = strfun.emptycellstr(1,neqtn);
EQTNS = strfun.emptycellstr(1,neqtn);
for ieq = 1 : neqtn
    if ~isempty(S.eqtnlhs{ieq})
        EQTNF{ieq} = [S.eqtnlhs{ieq},'-(',S.eqtnrhs{ieq},')'];
    else
        EQTNF{ieq} = S.eqtnrhs{ieq};
    end
    if ~isempty(S.sstaterhs{ieq})
        if ~isempty(S.sstatelhs{ieq})
            EQTNS{ieq} = [S.sstatelhs{ieq},'-(',S.sstaterhs{ieq},')'];
        else
            EQTNS{ieq} = S.sstaterhs{ieq};
        end
    end
end

    function dolossfunc()
        % dolossfunc  Find loss function amongst equations.
        start = regexp(S.eqtnrhs,'^min#?\(','once');
        lossindex = ~cellfun(@isempty,start);
        if sum(lossindex) == 1
            % Order the loss function last.
            list = {'eqtn','eqtnlabel','eqtnlhs','eqtnrhs','eqtnsign', ...
                'sstatelhs','sstaterhs','sstatesign'};
            for i = 1 : length(list)
                S.(list{i}) = [S.(list{i})(~lossindex), ...
                    S.(list{i})(lossindex)];
            end
            S.eqtnlhs{end} = '';
            S.eqtnrhs{end} = strrep(S.eqtnrhs{end},'#','');
            % Get the discount factor from inside of the min(...) brackets.
            [close,LOSSDISC] = strfun.matchbrk(S.eqtnrhs{end},4);
            % Remove the min operator.
            S.eqtnrhs{end} = S.eqtnrhs{end}(close+1:end);
        elseif sum(lossindex) > 1
            MULTIPLE = true;
        end
    end

end
% xxreadeqtns().

%**************************************************************************
function [THIS,LOGMISSING,INVALID,MULTIPLE] = xxreaddtrends(THIS,S)

n = sum(THIS.nametype == 1);
eqtn = strfun.emptycellstr(1,n);
eqtnF = strfun.emptycellstr(1,n);
eqtnlabel = strfun.emptycellstr(1,n);

% Create list of measurement variable names against which the LHS of
% dtrends equations will be matched. Add log(...) for log variables.
list = THIS.name(THIS.nametype == 1);
islog = THIS.log(THIS.nametype == 1);
loglist = list;
loglist(islog) = regexprep(loglist(islog),'(.*)','log($1)','once');

neqtn = length(S.eqtn);
logmissing = false(1,neqtn);
invalid = false(1,neqtn);
multiple = false(1,neqtn);
for ieqtn = 1 : length(S.eqtn)
    index = find(strcmp(loglist,S.eqtnlhs{ieqtn}),1);
    if isempty(index)
        if any(strcmp(list,S.eqtnlhs{ieqtn}))
            logmissing(ieqtn) = true;
        else
            invalid(ieqtn) = true;
        end
        continue
    end
    if ~isempty(eqtn{index})
        multiple(ieqtn) = true;
        continue
    end
    eqtn{index} = S.eqtn{ieqtn};
    eqtnF{index} = S.eqtnrhs{ieqtn};
    eqtnlabel{index} = S.eqtnlabel{ieqtn};
end

LOGMISSING = S.eqtn(logmissing);
INVALID = S.eqtn(invalid);
MULTIPLE = S.eqtnlhs(multiple);
if any(multiple)
    MULTIPLE = unique(MULTIPLE);
end

THIS.eqtn(end+(1:n)) = eqtn;
THIS.eqtnF(end+(1:n)) = eqtnF;
THIS.eqtnS(end+(1:n)) = {''};
THIS.eqtnlabel(end+(1:n)) = eqtnlabel;
THIS.eqtntype(end+(1:n)) = 3;
THIS.nonlin(end+(1:n)) = false;

end
% xxreaddtrends().

%**************************************************************************
function [THIS,INVALID] = xxreadlinks(THIS,S)

nname = length(THIS.name);
neqtn = length(S.eqtn);

valid = false(1,neqtn);
refresh = nan(1,neqtn);
for ieq = 1 : neqtn
    if isempty(S.eqtn{ieq})
        continue
    end
    index = strcmp(THIS.name,S.eqtnlhs{ieq});
    if any(index)
        % The LHS name is a variable, shock, or parameter.
        valid(ieq) = true;
        refresh(ieq) = find(index);
    else
        % The LHS name is a std or corr. The refresh pointer will be |nname+n|.
        index = mystdcorrindex(THIS,S.eqtnlhs{ieq});
        if any(index)
            valid(ieq) = true;
            refresh(ieq) = nname + find(index);
        end
    end
end

INVALID = S.eqtn(~valid);
THIS.eqtn(end+(1:neqtn)) = S.eqtn;
THIS.eqtnF(end+(1:neqtn)) = S.eqtnrhs;
THIS.eqtnS(end+(1:neqtn)) = {''};
THIS.eqtnlabel(end+(1:neqtn)) = S.eqtnlabel;
THIS.eqtntype(end+(1:neqtn)) = 4;
THIS.nonlin(end+(1:neqtn)) = false;
THIS.Refresh = refresh;

end
% xxreadlinks().

%**************************************************************************
function [m,eqtn] = xxtransformeqtn(m,eqtn,namecode,eqs)
% xxtransformeqtn  Replace numerical codes with x() and the names of stds
% and corrs with s().

ftransform = @doftransform; %#ok<NASGU>
stransform = @dostransform; %#ok<NASGU>
stdcorr = @dostdcorr; %#ok<NASGU>
pattern = ['([&])?([',namecode(1),'-',namecode(end),'])(\{[\+\-]\d+\})?'];
stdcorrpattern = '\<(std|corr)_[a-zA-Z]\w*\>';

nameoffset = double(namecode(1)) - 1;
tzero = m.tzero;

if ischar(eqtn) && ~isempty(eqtn)
    % Transform a single equation passed in as a text string.
    eqtn = regexprep(eqtn,pattern,'${ftransform($1,$2,$3,NaN)}');
    return
end

neqtn = length(m.eqtnF);
nname = length(m.name);
nt = size(m.occur,2) / nname;
occurF = reshape(full(m.occur),[neqtn,nname,nt]);
occurS = full(m.occurS);

if isequal(eqs,Inf)
    eqs = 1 : neqtn;
end

% We need to pass the equation number, `ieq`, into the nested functions. We
% therefore use a `for` loop, and not a single `regexprep` command.
for ieq = eqs
    
    if isempty(m.eqtnF{ieq})
        continue
    end
    
    % If no steady-state version exists, copy the dynamic equation.
    if ~m.linear && m.eqtntype(ieq) <= 2 && isempty(m.eqtnS{ieq})
        m.eqtnS{ieq} = m.eqtnF{ieq};
    end
    
    % Steady-state equations.
    if ~isempty(m.eqtnS{ieq})
        m.eqtnS{ieq} = regexprep(m.eqtnS{ieq},pattern, ...
            '${stransform($1,$2,$3,ieq)}');
    end
    
    % Full dynamic equations.
    m.eqtnF{ieq} = regexprep(m.eqtnF{ieq},pattern, ...
        '${ftransform($1,$2,$3,ieq)}');
        
    % Allow std_ and corr_ names only in dynamic links.
    if m.eqtntype(ieq) == 4 
        m.eqtnF{ieq} = regexprep(m.eqtnF{ieq}, ...
            stdcorrpattern,'${stdcorr($0)}');
    end
    
end

m.occur = sparse(occurF(:,:));
m.occurS = sparse(occurS);

    function c = doftransform(c0,c1,c2,ieq)
        % Replace name codes with x vector in dynamic equations.
        % c0 can be empty or '&'.
        % c1 is the name code, e.g. char(highChar+number).
        % c2 is empty or the time subscript, e.g. {-1}.
        % Variable or parameter number.
        realid = double(c1) - nameoffset;
        if realid < 1 || realid > length(m.nametype)
            % Undeclared name. Will be captured later.
            c = [c1,c2];
            return
        end
        if isempty(c2)
            t = 0;
        else
            c2 = c2(2:end-1);
            t = sscanf(c2,'%g');
        end
        if m.nametype(realid) <= 2
            % Measurement and transition variables.
            if t == 0
                time = 't';
            else
                time = sprintf('t%+g',t);
            end
            if isempty(c0)
                c = sprintf('x(:,%g,%s)',realid,time);
                if isfinite(ieq)
                    occurF(ieq,realid,tzero+t) = true;
                end
            else
                % Steady-state level reference.
                c = sprintf('L(%g)',realid);
            end
        elseif m.nametype(realid) == 3
            % Shocks.
            if isempty(c0)
                c = sprintf('x(:,%g,t)',realid);
                if isfinite(ieq)
                    occurF(ieq,realid,tzero+t) = true;
                end
            else
                c = '0';
            end
        else
            % Parameters.
            c = sprintf('x(:,%g,t)',realid);
            if isfinite(ieq)
                occurF(ieq,realid,tzero+t) = true;
            end
        end
    end
% doftransform().

    function c = dostransform(c0,c1,c2,ieq) %#ok<INUSL>
        % Replace name codes with x vector in sstate equations.
        % c0 is not used.
        % Variable or parameter number.
        realid = double(c1) - nameoffset;
        if realid < 1 || realid > length(m.nametype)
            % Undeclared name. Will be captured later.
            c = [c1,c2];
            return
        end
        if isempty(c2)
            t = 0;
        else
            c2 = c2(2:end-1);
            t = sscanf(c2,'%g');
        end
        switch m.nametype(realid)
            case 1
                % Measurement variables.
                c = sprintf('x(%g)',realid);
                if m.log(realid)
                    c = ['exp?(',c,')'];
                end
            case 2
                % Transition variables.
                if t == 0
                    if ~m.log(realid)
                        c = sprintf('x(%g)',realid);
                    else
                        c = sprintf('exp?(x(%g))',realid);
                    end
                else
                    if ~m.log(realid)
                        c = sprintf('(x(%g)+(%g)*dx(%g))',realid,t,realid);
                    else
                        c = sprintf('(exp?(x(%g))*exp?(dx(%g))^(%g))', ...
                            realid,realid,t);
                    end
                end
            case 3
                % Shocks.
                c = '0';
            case 4
                % Parameters.
                c = sprintf('x(%g)',realid);
        end
        if ~isinf(ieq)
            occurS(ieq,realid) = true;
        end
    end
% dostransform().

    function c = dostdcorr(c)
        index = mystdcorrindex(m,c);
        if any(index)
            n = find(index);
            c = sprintf('x(:,%g,t)',nname+n);
        end
    end
% dostdcorr().

end
% xxtransformeqtn().