function json = cherrypick(infile, outfold, eqtags, noresids, json, substitutetarget, donotpicktarget)

% Extract some equations in infile (mod file used for estimation)
% and write them in outfile (mod file used for simulation).
%
% INPUTS
% - infile            [string]    Name of the mod file where all the equations used for estimation are available.
% - outfold           [string]    Name of the folder where the generated files are saveda subset of the equations is to be printed.
% - eqtags            [cell]      Equation tags of the selected equations.
% - noresids          [logical]   Removes estimation residuals (not to be used in simulation) if true.
% - json              [char]      Content of a JSON file.
% - substitutetarget  [logical]   Substitute expression of a composite target in the error correction term of PAC equation if true. Default is false.
% - donotpicktarget   [logical]   Do not cherrypick auxiliary equation associated to the composite target. Default is false (equation is automagically cherrypicked).
%
% OUTPUTS
% - json              [char]      Content of a JSON file.
%
% SPECIAL REQUIREMENTS
% It is expected that the file infile.mod has already been run, and
% that the associated JSON output is available.

% Copyright © 2019-2025 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

global M_

% Set default value
if nargin<4 || isempty(noresids)
    noresids = true;
end

% Delete outfold subdirectory if it already exists
if isfolder(outfold)
    rmdir(outfold, 's');
end

% Create the subdirectoty where the generated files will be saved.
mkdir(outfold);

% Check that infile.mod and the related JSON output exist.
if ~isfile(sprintf('%s.mod', infile))
    error('Cannot find %s.mod.', infile)
end
if ~isfolder(sprintf('%s/model/json', infile))
    error('Cannot find %s/model/json folder. Did you run %s.mod with the json option?', infile, infile);
end

% Check if some variables have to be renamed.
rename = M_.equations_tags(strcmp('rename',M_.equations_tags(:,2)),[1,3]);
isrename = ~isempty(rename);

if nargin<5 || isempty(json)
    % Load json file (original mod file)
    json = loadjson_(sprintf('%s/model/json/modfile-original.json', M_.dname));
end

if nargin<6 || isempty(substitutetarget)
    substitutetarget = false;
    % Otherwise substitute target expression (in case of composite target) in the error correction term
end

if nargin<7 || isempty(donotpicktarget)
    donotpicktarget = false;
    % By default the auxiliary equation for the target in PAC equation is cherrypicked
end


% Create a new file.
fid = fopen(sprintf('%s/model.inc', outfold), 'w');

plist = {};
elist = {};
xlist = {};

try
    for i=1:length(eqtags)
        rhs = [];
        lhs = [];
        % Get equation number.
        eqnum = get_equation_number_by_tag(eqtags{i}, M_);
        % Get the original equation.
        [LHS, RHS] = get_lhs_and_rhs(eqtags{i}, M_, true, json);
        % Get the parameters, endogenous and exogenous variables in the current equation.
        [pnames, enames, xnames] = get_variables_and_parameters_in_equation(LHS, RHS, M_);
        lhs_expression = LHS;
        LHS = get_variables_and_parameters_in_expression(LHS);
        enames = union(enames, LHS);
        if length(LHS)>1
            error('Expressions with more than one variable on the LHS are not allowed.')
        end
        LHS = LHS{1};
        if isrename
            [variable_has_to_be_renamed, id] = ismember(eqnum, [rename{:,1}]);
            if variable_has_to_be_renamed
                TMP = strsplit(rename{id,2}, ',');
                for j=1:length(TMP)
                    tmp = strsplit(TMP{j}, '->');
                    lhs_expression = exactstrrep(lhs_expression, tmp{1}, tmp{2});
                    RHS = exactstrrep(RHS, tmp{1}, tmp{2});
                    rep = strcmp(tmp{1}, enames);
                    if any(rep)
                        enames(rep) = tmp(2);
                    end
                    rep = strcmp(tmp{1}, xnames);
                    if any(rep)
                        xnames(rep) = tmp(2);
                    end
                end
            end
        end
        % Remove residual from equation if required.
        if noresids
            if isfield(M_, 'simulation_exo_names')
                exogenous_variables_to_be_removed = ~ismember(xnames, M_.simulation_exo_names);
            else
                exogenous_variables_to_be_removed = false(size(xnames));
            end
            if any(exogenous_variables_to_be_removed)
                switch sum(exogenous_variables_to_be_removed)
                  case 1
                    RHS = regexprep(RHS, sprintf('\\ *\\+\\ *%s', xnames{exogenous_variables_to_be_removed}), '');
                    RHS = regexprep(RHS, sprintf('%s', xnames{exogenous_variables_to_be_removed}), '');
                  case 0
                    % Nothing to do.
                  otherwise
                    error('Cannot remove more than one exogenous variable in an equation (%s).', eqtags{i})
                end
                xnames = setdiff(xnames, xnames{exogenous_variables_to_be_removed});
            end
        end
        % Unroll expectation terms if any.
        isvar = regexp(RHS, 'var_expectation\(model_name = (?<name>\w+)\)', 'names');
        ispac = regexp(RHS, 'pac_expectation\(model_name = (?<name>\w+)\)', 'names');
        istar = regexp(RHS, 'pac_target_nonstationary\(model_name = (?<name>\w+)\)', 'names');
        if ~isempty(isvar)
            rhs = write_expectations(isvar.name, 'var');
            auxlhs = sprintf('%s_VE', eqtags{i});
            RHS = strrep(RHS, sprintf('var_expectation(model_name = %s)', isvar.name), auxlhs);
        end
        if ~isempty(ispac)
            if isfield(M_.pac.(ispac.name), 'components')
                [rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac', false, false);
            else
                [rhs, growthneutralitycorrection] = write_expectations(ispac.name, 'pac');
            end
            if ~isempty(rhs)
                if iscell(rhs) % PAC expectations are decomposed.
                    auxlhs = cell(size(rhs));
                    for k=1:length(M_.pac.(ispac.name).components)
                        if isequal(k, 1)
                            lhs = M_.lhs{M_.pac.(ispac.name).components(k).aux_id};
                            auxlhs{k} = lhs;
                            if ~isempty(growthneutralitycorrection{k})
                                lhs = sprintf('%s+%s', lhs, growthneutralitycorrection{k});
                            end
                            if  ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
                                if isempty(growthneutralitycorrection{k})
                                    lhs = sprintf('%s*%s', M_.pac.(ispac.name).components(k).coeff_str, lhs);
                                else
                                    lhs = sprintf('%s*(%s)', M_.pac.(ispac.name).components(k).coeff_str, lhs);
                                end
                            end
                        else
                            lhs_ = M_.lhs{M_.pac.(ispac.name).components(k).aux_id};
                            auxlhs{k} = lhs_;
                            if ~isempty(growthneutralitycorrection{k})
                                lhs_ = sprintf('%s+%s', lhs_, growthneutralitycorrection{k});
                            end
                            if  ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
                                if isempty(growthneutralitycorrection{k})
                                    lhs_ = sprintf('%s*%s', M_.pac.(ispac.name).components(k).coeff_str, lhs_);
                                else
                                    lhs_ = sprintf('%s*(%s)', M_.pac.(ispac.name).components(k).coeff_str, lhs_);
                                end
                            end
                            lhs = sprintf('%s+%s', lhs, lhs_);
                        end
                    end
                    RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), lhs);
                else
                    auxlhs = M_.lhs{M_.pac.(ispac.name).aux_id};
                    if isempty(growthneutralitycorrection)
                        RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), auxlhs);
                    else
                        RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', auxlhs, growthneutralitycorrection));
                    end
                end
            else
                % MCE version of the PAC equation.
                if ~isfield(M_.pac.(ispac.name), 'components')
                    auxlhs = M_.endo_names{M_.pac.(ispac.name).mce.z};
                    [rhs, growthneutralitycorrection] = write_pac_mce_expectations(eqtags{i}, ispac.name, auxlhs);
                    if isempty(growthneutralitycorrection)
                        RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), auxlhs);
                    else
                        RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), sprintf('%s+%s', auxlhs, growthneutralitycorrection));
                    end
                else
                    n = length(M_.pac.(ispac.name).components);
                    auxlhs = cell(n, 1);
                    rhs = cell(n, 1);
                    growthneutralitycorrection = cell(n, 1);
                    for c=1:n
                        auxlhs{c} = M_.endo_names{M_.pac.(ispac.name).mce.z(c)};
                        [rhs{c}, growthneutralitycorrection{c}] = write_pac_mce_expectations(eqtags{i}, ispac.name, auxlhs{c}, c);
                    end
                    if isempty(growthneutralitycorrection{1})
                        if strcmp(M_.pac.(ispac.name).components(1).coeff_str, '1')
                            RHS__ =  auxlhs{1};
                        else
                            RHS__ =  sprintf('%s*%s', M_.pac.(ispac.name).components(1).coeff_str, auxlhs{1});
                        end
                    else
                        if strcmp(M_.pac.(ispac.name).components(1).coeff_str, '1')
                            RHS__ = sprintf('%s+%s', auxlhs{1}, growthneutralitycorrection{1});
                        else
                            RHS__ = sprintf('%s*(%s+%s)', M_.pac.(ispac.name).components(1).coeff_str, auxlhs{1}, growthneutralitycorrection{1});
                        end
                    end
                    for c=2:n
                        if isempty(growthneutralitycorrection{c})
                            if strcmp(M_.pac.(ispac.name).components(c).coeff_str, '1')
                                RHS__ =  sprintf('%s+%s', RHS__, auxlhs{c});
                            else
                                RHS__ =  sprintf('%s+%s*%s', RHS__, M_.pac.(ispac.name).components(c).coeff_str, auxlhs{c});
                            end
                        else
                            if strcmp(M_.pac.(ispac.name).components(c).coeff_str, '1')
                                RHS__ = sprintf('%s+%s+%s', RHS__,  auxlhs{c}, growthneutralitycorrection{c});
                            else
                                RHS__ = sprintf('%s+%s*(%s+%s)', RHS__, M_.pac.(ispac.name).components(c).coeff_str, auxlhs{c}, growthneutralitycorrection{c});
                            end
                        end
                    end
                    RHS = strrep(RHS, sprintf('pac_expectation(model_name = %s)', ispac.name), RHS__);
                end
            end
        end
        if ~isempty(istar)
            if substitutetarget
                RHS = strrep(RHS, sprintf('pac_target_nonstationary(model_name = %s)', ispac.name), targetexpr());
            else
                if ~donotpicktarget
                    auxlhs{end+1} = M_.endo_names{M_.pac.(ispac.name).ec.vars(M_.pac.(ispac.name).ec.istarget)};
                    rhs{end+1} = M_.aux_vars(strmatch(auxlhs{end}, M_.endo_names, 'exact')==[M_.aux_vars(:).endo_index]).orig_expr;
                    rhs{end} = remove_aux_variables_from_expression(rhs{end}, M_);
                    RHS = strrep(RHS, sprintf('pac_target_nonstationary(model_name = %s)', ispac.name), sprintf('%s(-1)', auxlhs{end}));
                else
                    RHS = strrep(RHS, sprintf('pac_target_nonstationary(model_name = %s)', ispac.name), sprintf('%s(-1)', M_.endo_names{M_.pac.(ispac.name).ec.vars(M_.pac.(ispac.name).ec.istarget)}));
                end
            end
        end
        % Print equation for unrolled PAC/VAR-expectation and update
        % list of parameters and endogenous variables (if any).
        if ~isempty(rhs)
            % Note that the call to get_variables_and_parameters_in_equation()
            % will not return the lhs variable in expectation_enames since
            % the name is created on the fly and is not a  member of M_.endo_names.
            expectation_pnames = get_variables_and_parameters_in_equation('', rhs, M_);
            expectation_enames = get_variables_and_parameters_in_expression(auxlhs);
            expectation_xnames = get_variables_and_parameters_in_expression(rhs);
            pnames = union(pnames, expectation_pnames);
            xnames = union(xnames, setdiff(expectation_xnames, expectation_pnames));
            enames = union(enames, expectation_enames);
            if ischar(rhs)
                fprintf(fid, '[name=''%s'']\n', auxlhs);
                fprintf(fid, '%s = %s;\n\n', auxlhs, rhs);
            else
                for k=1:length(rhs)
                    fprintf(fid, '[name=''%s'']\n', auxlhs{k});
                    fprintf(fid, '%s = %s;\n\n', auxlhs{k}, rhs{k});
                end
            end
            if ~isempty(ispac) && isfield(M_.pac.(ispac.name), 'components')
                for k=1:length(M_.pac.(ispac.name).components)
                    if ~isequal(M_.pac.(ispac.name).components(k).coeff_str, '1')
                        params = get_variables_and_parameters_in_expression(M_.pac.(ispac.name).components(k).coeff_str);
                        pnames = union(pnames, params);
                    end
                end
            end
        else
            pRHS = get_variables_and_parameters_in_equation('', RHS, M_);
            xRHS = get_variables_and_parameters_in_expression(RHS);
            xnames = union(xnames, setdiff(xRHS, pRHS));
            pnames = union(pnames, pRHS);
        end
        % Update pnames, enames and xnames if PAC with growth neutrality correction.
        if ~isempty(ispac) && ~isempty(growthneutralitycorrection)
            if iscell(growthneutralitycorrection)
                growthneutralitycorrection_pnames = [];
                growthneutralitycorrection_enames = [];
                growthneutralitycorrection_xnames = [];
                for component=1:length(growthneutralitycorrection)
                    if ~isempty(growthneutralitycorrection{component})
                        [tmp_pnames, tmp_enames, tmp_xnames] = get_variables_and_parameters_in_equation('', growthneutralitycorrection{component}, M_);
                        growthneutralitycorrection_pnames = union(growthneutralitycorrection_pnames, tmp_pnames);
                        growthneutralitycorrection_xnames = union(growthneutralitycorrection_xnames, tmp_xnames);
                        growthneutralitycorrection_enames = union(growthneutralitycorrection_enames, tmp_enames);
                    end
                end
            else
                [growthneutralitycorrection_pnames, ...
                 growthneutralitycorrection_enames, ...
                 growthneutralitycorrection_xnames] = get_variables_and_parameters_in_equation('', growthneutralitycorrection, M_);
            end
            if ~isempty(growthneutralitycorrection_pnames)
                pnames = union(pnames, growthneutralitycorrection_pnames);
            end
            if ~isempty(growthneutralitycorrection_enames)
                xnames = union(xnames, growthneutralitycorrection_enames);
            end
            if ~isempty(growthneutralitycorrection_xnames)
                xnames = union(xnames, growthneutralitycorrection_xnames);
            end
        end
        % Print tags
        if iscell(json.model)
            tfields = fieldnames(json.model{eqnum}.tags);
            tags = sprintf('%s=''%s''', tfields{1}, json.model{eqnum}.tags.(tfields{1}));
            for j=2:length(tfields)
                if ~isempty(json.model{eqnum}.tags.(tfields{j}))
                    tags = sprintf('%s, %s=''%s''', tags, tfields{j}, json.model{eqnum}.tags.(tfields{j}));
                end
            end
        else
            tfields = fieldnames(json.model.tags);
            tags = sprintf('%s=''%s''', tfields{1}, json.model.tags.(tfields{1}));
            for j=2:length(tfields)
                if ~isempty(json.model.tags.(tfields{j}))
                    tags = sprintf('%s, %s=''%s''', tags, tfields{j}, json.model.tags.(tfields{j}));
                end
            end
        end
        fprintf(fid, '[%s]\n', tags);
        % Print equation.
        fprintf(fid, '%s = %s;\n\n', lhs_expression, RHS);
        % Update lists of parameters, endogenous variables and exogenous variables.
        plist = union(plist, pnames);
        elist = union(elist, enames);
            xlist = union(xlist, xnames);
        end
catch e
    fclose(fid);
    fprintf(2, '%s\n', e.message)
    fprintf(2, 'Delete ''%s'' file.\n', sprintf('%s/model.inc', outfold))
    delete(sprintf('%s/model.inc', outfold))
    return
end

fclose(fid);

% Export parameters
if ~isempty(plist)
    fid = fopen(sprintf('%s/parameters.inc', outfold), 'w');
    fprintf(fid, 'parameters %s;', sprintf('%s ', plist{:}));
    fclose(fid);
end

% Export endogegnous variables
fid = fopen(sprintf('%s/endogenous.inc', outfold), 'w');
printlistofvariables(fid, 'endo', elist, M_, elist);
fclose(fid);

% Export exogenous variables
if ~isempty(xlist)
    fid = fopen(sprintf('%s/exogenous.inc', outfold), 'w');
    printlistofvariables(fid, 'exo', xlist, M_, xlist);
    fclose(fid);
end

% Export parameter values
if ~isempty(plist)
    fid = fopen(sprintf('%s/parameter-values.inc', outfold), 'w');
    for i=1:length(plist)
        id = strcmp(plist{i}, M_.param_names);
        if any(id)
            if isnan(M_.params(id))
                warning('Parameter %s has no value.', plist{i})
            else
                fprintf(fid, '%s = %s;\n', plist{i}, num2str(M_.params(id), 16));
            end
        end
    end
    fclose(fid);
end


function printlistofvariables(fid, kind, list, DynareModel, vappend)
    if isfield(DynareModel, sprintf('%s_partitions', kind))
        % Some endogenous variables are tagged.
        switch kind
          case 'exo'
            tfields = fieldnames(DynareModel.exo_partitions);
            vlist = 'varexo';
            vnames = DynareModel.exo_names;
            partitions = DynareModel.exo_partitions;
          case 'endo'
            tfields = fieldnames(DynareModel.endo_partitions);
            vlist = 'var';
            vnames = DynareModel.endo_names(1:DynareModel.orig_endo_nbr);
            partitions = DynareModel.endo_partitions;
          otherwise
            error('Illegal value for second input argument.')
        end
        for i = 1:length(list)
            id = strmatch(list{i}, vnames, 'exact');
            if ~isempty(id)
                tags = '';
                for j=1:length(tfields)
                    if ~isempty(partitions.(tfields{j}){id})
                        tags = sprintf('%s, %s=''%s''', tags, tfields{j}, partitions.(tfields{j}){id});
                    end
                end
                if ~isempty(tags)
                    tags = sprintf('(%s)', tags(3:end));
                end
            elseif ~isempty(strmatch(list{i}, vappend, 'exact'))
                % Nothing to do, this variable was renamed by cherrypick
                tags = '';
            else
                if isequal(kind, 'endo') &&  (isequal(list{i}(end-2:end), '_PE') || isequal(list{i}(end-2:end), '_VE'))
                    if isequal(list{i}(end-2:end), '_PE')
                        tags = sprintf('(expectation_kind=''%s'')', 'pac');
                    else
                        tags = sprintf('(expectation_kind=''%s'')', 'var');
                    end
                else
                    error('Unknown variable.')
                end
            end
            if isempty(tags)
                vlist = sprintf('%s\n\t%s', vlist, list{i});
            else
                vlist = sprintf('%s\n\t%s %s', vlist, list{i}, tags);
            end
        end
        fprintf(fid, '%s;', vlist);
    else
        switch kind
          case 'exo'
            vlist = 'varexo';
          case 'endo'
            vlist = 'var';
          otherwise
            error('Illegal value for second input argument.')
        end
        for i=1:length(list)
            vlist = sprintf('%s\n\t%s', vlist, list{i});
        end
        fprintf(fid, '%s;', vlist);
    end
end

function expr = targetexpr()
    if isfield(M_.pac.(ispac.name), 'components')
        ns = ~strcmp({M_.pac.(ispac.name).components.kind}, 'll'); % non stationary components
        expr = '';
        for i=1:length(M_.pac.(ispac.name).components)
            if ns(i)
                variable = rmauxiliary(M_.endo_names{M_.pac.(ispac.name).components(i).endo_var}, 1);
                if isempty(expr)
                    if strcmp(M_.pac.(ispac.name).components(i).coeff_str, '1')
                        expr = variable;
                    else
                        expr = sprintf('%s*%s', M_.pac.(ispac.name).components(i).coeff_str, variable);
                    end
                else
                    if strcmp(M_.pac.(ispac.name).components(i).coeff_str, '1')
                        expr = sprintf('%s+%s', expr, variable);
                    else
                        expr = sprintf('%s+%s*%s', expr, M_.pac.(ispac.name).components(i).coeff_str, variable);
                    end
                end
            end
        end
    else
        expr = rmauxiliary(M_.endo_names{M_.pac.(ispac.name).ec.vars(M_.pac.(ispac.name).ec.istarget)}, 1);
    end
end


function variable = rmauxiliary(variable, lag)
    transformations = {};
    if isauxiliary(variable)
        ida = get_aux_variable_id(variable);
        op = 0;
        while ida
            op = op+1;
            if isequal(M_.aux_vars(ida).type, 8)
                transformations(op) = {'diff'};
                variable = M_.endo_names{M_.aux_vars(ida).orig_index};
                ida = get_aux_variable_id(variable);
            elseif isequal(M_.aux_vars(ida).type, 10)
                transformations(op) = {M_.aux_vars(ida).unary_op};
                variable = M_.endo_names{M_.aux_vars(ida).orig_index};
                ida = get_aux_variable_id(variable);
            else
                error('This case is not implemented.')
            end
        end
    end
    if nargin>1
        variable = sprintf('%s(-%u)', variable, lag);
    end
    for k=length(transformations):-1:1
        variable = sprintf('%s(%s)', transformations{k}, variable);
    end
end

end