function [this,assign] = parse(this,p,assign,std)

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

this.fname = p.fname;
c = p.code;

s = sprintf( ...
    'parsing file(s) <a href="matlab: edit %s">%s</a>. ', ...
    strrep(this.fname,' & ',' '),this.fname);
errorParsing = ['Error ',s];
% warningParsing = ['Warning ',s];

blkKeyword = {
    '!endogenous_variables', ...
    '!exogenous_variables', ...
    '!residuals', ...
    '!parameters', ...
    '!equations', ...
    '!log_variables', ...
    };

nblk = length(blkKeyword);

endOfBlk = '(?=$';
for i = 1 : nblk
    endOfBlk = [endOfBlk,'|',blkKeyword{i}];
end
endOfBlk = [endOfBlk,')'];

blk = cell(1,nblk);
for i = 1 : nblk
    match = regexp(c,['(?<=',blkKeyword{i},').*?',endOfBlk],'match');
    blk{i} = [match{:}];
    if isempty(blk{i})
        blk{i} = '';
    end
end

%**************************************************************************
% Names.
name = cell(1,4);
namelabel = cell(1,4);
value = cell(1,4);
nametype = cell(1,4);
for i = 1 : 4
    if ~isempty(blk{i})
        [name{i},namelabel{i},value{i}] = theparser.parsenames(blk{i});
        nametype{i} = i*ones(1,length(name{i}));
    else
        name{i} = {};
        namelabel{i} = {};
        value{i} = [];
        nametype{i} = [];
    end
end
name = [name{:}];
namelabel = [namelabel{:}];
value = [value{:}];
nametype = [nametype{:}];

name0 = name;
[name,index] = unique(name);
if length(name) ~= length(name0)
    index = setdiff(1:length(name0),index);
    utils.error('bkwmodel', ...
        'This name is declared more than once: ''%s''.', ...
        name0{index});
end
name = name0;
nname = length(name);
ny = sum(nametype == 1);

%**************************************************************************
% Log variables.

[allbut,blk{6}] = strfun.findremove(blk{6},'!allbut');
islog = false(size(name));
if allbut
    islog(nametype == 1) = ~islog(nametype == 1);
end
loglist = regexp(blk{6},'\<[a-zA-Z]\w*\>','match');
loglist = unique(loglist);
valid = true(size(loglist));
for i = 1 : length(loglist)
    index = strcmp(loglist{i},name);
    if any(index) && nametype(index) == 1
        islog(index) = ~islog(index);
    else
        valid(i) = false;
    end
end
if any(~valid)
    utils.error('bkwmodel', ...
        ['This name found in the log-variables section ', ...
        'is not an endgenous variable: ''%s''.'], ...
        loglist{~valid});
end

%**************************************************************************
% Equations.
blk{5} = regexprep(blk{5},'\s+','');
% Do not include the semi-colon; it will be added when re-arranging LHS and
% RHS.
tokens = regexp(blk{5},'(#\(\d+\))?([^#;][^;]+);','tokens');
neqtn = length(tokens);
eqtn = cell(1,neqtn);
eqtnlabel = cell(1,neqtn);
for i = 1 : length(tokens)
    eqtn{i} = tokens{i}{2};
    if ~isempty(tokens{i}{1})
        eqtnlabel{i} = tokens{i}{1};
    else
        eqtnlabel{i} = '';
    end
end

if neqtn ~= sum(nametype == 1)
    utils.error('bkwmodel', ...
        [errorParsing, ...
        'The number of equations (%g) does not match the number ', ...
        'of endogenous variables (%g).'], ...
        neqtn,sum(nametype == 1));
end

% Remove white spaces.
eqtnEval = eqtn;
eqtnEval = regexprep(eqtnEval,'\s+','');

% Convert names to codes.
occurpos = [];
cannotParse = {};
leadsNotAllowed = {};
lagsNotAllowed = {};
replace = @name2code_; %#ok<NASGU>

for ieq = 1 : neqtn
    eqtnEval{ieq} = regexprep(eqtnEval{ieq}, ...
        '(\<[a-zA-Z]\w*\>)(\{[^\}]+\})?(?!\()', ...
        '${replace($0,$1,$2,true)}');
end

    % @ *******************************************************************
    function c = name2code_(thisexpression,thisname,thislag,capture)
        c = '';
        if isempty(thislag)
            thislag = 0;
        else
            thislag = -round(sscanf(thislag,'{%g}'));
        end
        thisid = find(strcmp(thisname,name));
        if isempty(thisid) || ~isnumericscalar(thislag)
            cannotParse{end+1} = thisexpression;
            return
        end
        if thislag < 0
            leadsNotAllowed{end+1} = thisexpression;
            return
        end
        if nametype(thisid) > 2 && thislag ~= 0
            lagsNotAllowed{end+1} = thisexpression;
            return
        end
        if thislag == 0
            c = ['x(',sprintf('%g',thisid),',t)'];
        else
            c = ['x(', ...
                sprintf('%g',thisid),',t-',sprintf('%g',thislag),')'];
        end
        if ~capture
            return
        end
        occurpos = [occurpos;ieq,thisid,thislag];
    end
    % @ name2code_().

if ~isempty(cannotParse)
    utils.error('bkwmodel', ...
        'Cannot parse this expression: ''%s''.',cannotParse{:});
end

if ~isempty(leadsNotAllowed)
    utils.error('bkwmodel', ...
        'Leads not allowed in bkwmodels: ''%s'',',leadsNotAllowed{:});
end

if ~isempty(lagsNotAllowed)
    utils.error('bkwmodel', ...
        'Lags only allowed for endogenous variables in bkwmodels: ''%s'',', ...
        lagsNotAllowed{:});
end

maxlag = max(occurpos(:,3));
if maxlag == 0
    maxlag = 1;
end
occur = false(neqtn,nname,maxlag);

for i = 1 : size(occurpos,1)
    occur(occurpos(i,1),occurpos(i,2),1+occurpos(i,3)) = true;
end

this.name = name;
this.nametype = nametype;
this.namelabel = labelsback(p,namelabel,'%s');
this.log = islog;
this.eqtnUsr = eqtn;
this.eqtnlabel = labelsback(p,eqtnlabel,'%s');
this.occur = occur;

% Set up the state vector.
do_metasystem();

% Solve for current dates of variables, and differentiate the equations wrt
% lagged endogenous variables.
do_findident();
do_diff();

this.eqtnEval = eqtnEval;
this.dimplicit = dimplicit;

this.Assign = nan(1,nname);
ne = sum(nametype == 3);
this.stdcorr = nan(1,ne + ne*(ne-1)/2);
this.stdcorr(1,1:ne) = std;
this.stdcorr(1,ne+1:end) = 0;

% Update the assign database with the values from the model code.
if isempty(assign)
    assign = struct();
end
for i = 1 : nname
    if isfield(assign,name{i})
        continue
    end
    thisvalue = []; %#ok<NASGU>
    try
        thisvalue = eval(value{i});
    catch %#ok<CTCH>
        thisvalue = NaN;
    end
    if ~isnumeric(thisvalue) || isempty(thisvalue)
        thisvalue = NaN;
    end
    if ~isnan(thisvalue)
        assign.(name{i}) = thisvalue;
    end
end

% Nested functions follow.
    
%{
    % @ *******************************************************************
    function findident_()
        xlist = name(nametype == 1);
        valid = true(1,neqtn);
        for ieq = 1 : neqtn
            pos = strfind(eqtnEval{ieq},':=');
            if isempty(pos)
                continue
            end
            tok = regexp(eqtnEval{ieq}, ...
                '^x\((\d+),t\):=(.*)$','tokens','once');
            if isempty(tok) || length(tok) ~= 2 ...
                    || isempty(tok{1}) || isempty(tok{2})
                valid(ieq) = false;
                continue
            end
            index = str2double(tok{1});
            rhs = strtrim(tok{2});
            eqtnEval{ieq} = strrep(eqtnEval{ieq},':=','=');
            if nametype(index) ~= 1                
                valid(ieq) = false;
                continue
            end
            identpoint(end+1) = index; %#ok<AGROW>
            eqtnident{end+1} = rhs; %#ok<AGROW>
        end
        if any(~valid)
            utils.error('rhsmodel', ...
                'Invalid LHS of this identity: ''%s''.', ...
                eqtn(~valid));
        end
    end
    % @ findident_().
%}
    
%**************************************************************************
    function do_findident()
        % Find identities because they need to be excluded when we solve
        % the system for residuals. Identities have := instead of =.
        n = cellfun(@length,eqtnEval);
        % Replace := with =.
        eqtnEval = strrep(eqtnEval,':=','=');
        this.isident = cellfun(@length,eqtnEval) ~= n;        
    end
    % do_findident


%**************************************************************************
    function do_diff()
    % Differentiate the implicit (user) equations wrt all occurences of
    % current and lagged variables, and mutliply the derivatives by the
    % respective variables if it is a log-variable. Create alpha vector.
        eqtnEval = regexprep(eqtnEval,'([^=]+)?=(.*)','$1-($2)');
        eqtnEval = do_code2symb(eqtnEval);
        symbeqtn = eqtnEval;
        for i = 1 : length(symbeqtn)
            symbeqtn{i} = sym(symbeqtn{i});
        end
        % Differentiate the implicit (user) equations wrt to current and
        % lagged endogenous variables and current residuals.
        [oname,olag] = findoccur(this);
        no = length(oname);
        dimplicit = cell(neqtn,no);
        for id = 1 : no
            if olag(id) == 0
                wrt = sprintf('x%g',oname(id));
            else
                wrt = sprintf('x%g_%g',oname(id),olag(id));
            end
            for ii = 1 : neqtn
                if ~this.occur(ii,oname,olag+1)
                    d = '0';
                else
                    sy = sym(symbeqtn{ii});
                    d = char(diff(sy,wrt));
                    if this.log(oname(id)) && ~isequal(d,'0')
                        d = ['(',d,')*',wrt]; %#ok<AGROW>
                    end
                end
                dimplicit{ii,id} = d;
            end
        end
        eqtnEval = do_symb2code(eqtnEval);
        dimplicit = do_symb2code(dimplicit);
        
        function x = do_code2symb(x)
            x = regexprep(x,'x\((\d+),t-(\d+)\)','x$1_$2');
            x = regexprep(x,'x\((\d+),t\)','x$1');
        end
        
        function x = do_symb2code(x)
            % Replace lags of endogenous variables with the corresponding
            % element of the system vector a. Otherwise, keep references to
            % exogenous variables and residuals.
            na = length(this.systemid);
            for ia = 1 : na
                realid = real(this.systemid(ia));
                imagid = imag(this.systemid(ia));
                pattern = sprintf('x%g_%g',realid,abs(imagid-1));
                replace = sprintf('a(%g,t-1)',ia);
                if this.systemlog(ia)
                    replace = ['exp(',replace,')']; %#ok<AGROW>
                end
                x = regexprep(x,['\<',pattern,'\>'],replace);
            end
            
            x = regexprep(x,'x(\d+)_(\d+)','x($1,t-$2)');
            
            % Current dated endogenous variables.
            for ia = 1 : ny
                pattern = sprintf('x%g',ia);
                replace = sprintf('a(%g,t)',ia);
                if this.systemlog(ia)
                    replace = ['exp(',replace,')']; %#ok<AGROW>
                end
                x = regexprep(x,['\<',pattern,'\>'],replace);
            end
            x = regexprep(x,'x(\d+)','x($1,t)');
        end
        
        
    end
    % do_diff

%**************************************************************************
    function do_metasystem()
        nt = size(this.occur,3);
        maxlag = nt - 1;
        
        this.systemid = [];
        for t = 0 : maxlag-1
            this.systemid = [this.systemid,find(this.nametype == 1) - 1i*t];
        end
        
        this.systemkeep = true(size(this.systemid));
        for ii = find(this.nametype == 1)
            thismaxlag = max([1,find(any(this.occur(:,ii,:),1),1,'last')-1]);
            index = real(this.systemid) == ii ...
                & abs(imag(this.systemid)) > thismaxlag-1;
            if any(index)
                this.systemkeep(index) = false;
            end
        end
        this.systemid = this.systemid(this.systemkeep);
        
        % Map the X vector to A vector. Find positions of A vector elements
        % in a window of X observations t-maxlag:t.
        this.systempos = nan(size(this.systemid));
        for ii = 1 : length(this.systemid)
            w = false(ny,maxlag+1);
            realid = real(this.systemid(ii));
            imagid = imag(this.systemid(ii));
            w(realid,end+imagid) = true;
            this.systempos(ii) = find(w);
        end
        
        %{
        Create logical vector systemlog: systemlog an index of log variables in
        the system vector a. It works as follows:

            a(this.systemlog,:) = log(a(this.systemlog,:));
        %}
        na = sum(this.systemkeep);
        this.systemlog = false(1,na);
        for ii = 1 : na
            realid = real(this.systemid(ii));
            this.systemlog(ii) = this.log(realid);
        end
        
        %{
        Create matrix J: J is a mapping from the lagged system vector a(t-1) to
        the lower part of the current-date system vector a(t) ; lower part
        means all lagged variables. It works as follows:

            a(ny+1:end,t) = J*a(:,t-1);
        %}
        nadd = ny*(maxlag-1);
        J = [nan(ny,nadd+ny);eye(nadd,nadd+ny)];
        J = J(this.systemkeep,this.systemkeep);
        this.systemj = J(ny+1:end,:);
        
    end
    % do_metasystem
    
end