function [db, info] = simulatelongrun(modelname, db, maxperiods, info, cleanupflag, normalize, options)

% Computes the BGP of a model by simulating the model..
%
% INPUTS
% - modelname        [char]              1×n  Name of a mod file (with extension).
% - db               [dseries]           database used as an initial condition.
% - maxperiods       [integer]           scalar, number of periods in the simulation.
% - info             [cell,dseries]      m×3 array (returned by data.projectionsetup) describing the path for the exogenous
%                                        variables, or dseries objects containg paths for all the exogenous variables.
% - cleanupflag      [logical]           scalar, remove generated files if true.
% - options          [struct]            options for the nonlinear solver.
%
% OUTPUTS
% - db               [dseries]           updated database with simulation.
% - info             [cell]              updated info argument ()
global M_ options_

if nargin<3
    error('Three input arguments required (mod file name, database, ')
end

[filepath, basename, ext] = fileparts(modelname);

if ~isempty(ext) && ~ismember(ext, {'.mod','.dyn'})
    error('File %s is not a Dynare file.', modelname)
end

if ~isempty(ext)
    if ~exist(modelname, 'file')
        error('File %s not found.', modelname)
    end
else
    if ~exist(sprintf('%s.mod', modelname), 'file')
        error('Dynare file extension was not provided and %s.mod cannot be found.', modelname)
    end
end

if ~isscalar(maxperiods) || ~isint(maxperiods) || maxperiods<=0
    error('Third argument (the simulation max length) has to be a positive scalar integer.')
end

if nargin>3
    if ~(iscell(info) && isequal(columns(info), 4)) && ~isdseries(info)
        error('Fourth input argument has to be a cell array with four columns or a dseries object.')
    end
end

if nargin<3 || isempty(info)
    disp('Default path for the exogenous variables is used.')
    info = data.projectionsetup();
    cleanupflag = true;
end

if nargin<4
    cleanupflag = true;
end

normp = 1500;

if maxperiods>normp;
    % Dynare's simulate_backward_model will be called more than once (every normp periods).
    np = floor(maxperiods/normp);
    pp = normp*ones(np, 1);
    rr = rem(maxperiods, normp);
    if rr
        pp = [pp; rr];
    end
else
    pp = maxperiods;
end

% Set path
if ~isempty(filepath)
    origPath = pwd();
    cd(filepath)
end

% Call dynare.
dynare(basename, 'nopreprocessoroutput', 'notime', 'stochastic', 'json=transform')

if nargin>5
    options_.solve_algo = options.solve_algo;
    options_.steady.maxit = options.maxiter;
    options_.solve_tolx = options.tolx;
    options_.solve_tolf = options.tolf;
    if isfield(options, 'simulations_along_bgp')
        nbgp = options.simulations_along_bgp;
    else
        nbgp = 1;
    end
end

% Define periods for the initial conditions. Simulation will start
% after the last period in db, but we don't need to keep all the
% previous periods in the sample if they are not required for the
% initial conditions.
initialconditions = db(db.firstobservedperiod():db.lastobservedperiod());

% Project the exogenous variables if necessary or use the paths
% provided in the fourth argument.
if ~isdseries(info)
    exogenous = initialconditions{M_.exo_names{:}};         % Extract all the exogenous variables.
    exogenous.projection(info, maxperiods);                 % Call the projection method.
else
    if periods>info.nobs
        exogenous = info{M_.exo_names{:}}(info.dates(1:end));
    else
        exogenous = info{M_.exo_names{:}}(info.dates(1:periods));
    end
end

% Remove first observations in exogenous object (corresponding to
% the periods where the original sample is defined).
exogenous = exogenous(db.lastobservedperiod()+1:exogenous.lastdate());

% Read JSON output to identify the equations with a logged variable
% on the LHS.
if ismember(options_.solve_algo, [12,14])
    json = loadjson_(sprintf('%s/model/json/modfile.json', basename));
    lhs = cell(length(json.model),1);
    isauxdiffloggedrhs = false(length(json.model), 1);
    for i = 1:length(json.model)
        lhs{i} = json.model{i}.lhs;
        if i>M_.orig_endo_nbr &&  ~isempty(regexp(lhs{i}, '\<AUX_DIFF_(\d*)\>', 'once')) && ismember(lhs{i}, lhs(1:i-1)) && ...
                ~isempty(regexp(json.model{i}.rhs, 'log\(\w*\)-log\(\w*\(-1\)\)', 'once'))
            isauxdiffloggedrhs(i) = true;
        end
    end
    islog = @(x) ~isempty(regexp(x, 'log\(\w*\)', 'once'));
    M_.isloggedlhs = cellfun(islog, lhs);
    M_.lhs = lhs;
    M_.isauxdiffloggedrhs = isauxdiffloggedrhs;
end

% List of the variables to be updated in the database.
modelvariables = [M_.endo_names(1:M_.orig_endo_nbr); exogenous.name];
endogenousvariables = modelvariables(1:M_.orig_endo_nbr);

simdata = initialconditions{modelvariables{:}};

% Simulate the model
initp = exogenous.dates(1);
timep = 1;
condp = simdata;
sizep = condp.nobs;
exogp = exogenous(initp+(0:pp(1)-1));

convergence = false;
alongbgp = 0;

for i = 1:length(pp)
    if i>1 && normalize
        % Deflate trended variables if more than one iteration.
        condp = deflate(simdata(simdata.dates(end)-sizep+1:simdata.dates(end)), info, timep);
        exogp = deflate(exogenous(initp+(0:pp(i)-1)), info, timep);
    else
        exogp = exogenous(initp+(0:pp(1)-1));
    end
    names = condp.name;
    % Simulate model with deflated variables.
    tmpdata = simul_backward_model(condp, pp(i), exogp);
    % Set initial conditions for next iteration.
    if i>1 && normalize
        % Inflate variables if previously deflated.
        condp = inflate(tmpdata{names{:}}(tmpdata.dates(end-normp+1:end)), info, timep);
    else
        condp = tmpdata{names{:}}(tmpdata.dates(end-normp+1:end));
    end
    % Merge simulations in output argument.
    simdata = merge(simdata, condp{modelvariables{:}});
    % Check convergence to the BGP and exit loop if convergence is achieved
    if checkconvergence(tmpdata{endogenousvariables{:}}, info, M_.endo_names(1:M_.orig_endo_nbr), 1e-5, normalize)
        convergence = true;
        alongbgp = alongbgp+1;
        if isequal(alongbgp, nbgp)
            break
        end
    end
    % Update periods and number of simulations.
    initp = initp+pp(i);
    timep = sum(pp(1:i));
end

% Report results.
if convergence
    dprintf('Convergence achieved after %i periods.', timep);
end

% Update the database.
db = merge(db{modelvariables{:}}, simdata);

% Cleanup
if cleanupflag
    rmdir(sprintf('+%s', basename), 's');
    rmdir(sprintf('%s', basename), 's');
    delete(sprintf('%s.log', basename));
    delete(sprintf('%s_results.mat', basename));
end

% (re)Set path
if ~isempty(filepath)
    cd(origPath)
end

if nargout>1
    % Update NaN values in info and check that we find results consistent with non NaN values.
    idConstWithNaNs = find(strcmp(info(:,2), 'Constant') & cellfun(@isnan, info(:,3)));
    nameConstWithNaNs = info(idConstWithNaNs, 1);
    nameConstWithNaNs = intersect(nameConstWithNaNs, db.name);
    idConstWithValues = find(strcmp(info(:,2), 'Constant') & ~cellfun(@isnan, info(:,3)));
    nameConstWithValues = info(idConstWithValues, 1);
    nameConstWithValues = intersect(nameConstWithValues, db.name);
    for i=1:length(nameConstWithNaNs)
        info(strcmp(info(:,1), nameConstWithNaNs{i}),3) = {mean(db{nameConstWithNaNs{i}}.data(end-round(.1*normp):end))};
    end
    for i=1:length(nameConstWithValues)
        if abs(info{strcmp(info(:,1), nameConstWithValues{i}),3}-mean(db{nameConstWithValues{i}}.data(end-round(.1*normp):end)))>1e-5
            dprintf('Unexpected result for %s (%s vs %s)', nameConstWithValues{i}, num2str(info{strcmp(info(:,1), nameConstWithValues{i}),3}, 6), num2str(mean(db{nameConstWithValues{i}}.data(end-round(.1*normp):end)), 6))
            info(strcmp(info(:,1), nameConstWithValues{i}),3) = {mean(db{nameConstWithValues{i}}.data(end-round(.1*normp):end))};
        end
    end
end

end

function conv = checkconvergence(paths, info, endonames, precision, growth)
    p = 1.0;
    conv = true;
    idConst = find(strcmp(info(:,2), 'Constant'));
    nameConst = info(idConst, 1);
    nameConst = intersect(nameConst, paths.name);
    dataConst = paths{nameConst{:}}.data(end-round(p*paths.nobs)+1:end,:);
    diffConst = max(abs(dataConst(2:end,:)-dataConst(1:end-1,:)));
    [testConst, idTest] = max(diffConst);
    if testConst>precision
        conv = false;
        dprintf('Variable %s is still moving (%s)', nameConst{idTest}, testConst)
    end
    if growth
        conv = testTrends(conv, 'Real', paths, info, precision);
        conv = testTrends(conv, 'Price', paths, info, precision);
        conv = testTrends(conv, 'Nominal', paths, info, precision);
        conv = testTrends(conv, 'NominalPerCapita', paths, info, precision);
        conv = testTrends(conv, 'Tfp', paths, info, precision);
        conv = testTrends(conv, 'Population', paths, info, precision);
    end
end

function conv = testTrends(conv, type, paths, info, precision)
    if conv
        p = 0.2;
        idTrend = find(strcmp(info(:,2), 'Trend') & strcmp(info(:,4), type));
        nameTrend = info(idTrend, 1);
        nameTrend = intersect(nameTrend, paths.name);
        dataTrend = paths{nameTrend{:}}.data;
        diffTrend = max(dataTrend(end-round(p*paths.nobs)+2:end,:)./dataTrend(end-round(p*paths.nobs)+1:end-1,:)-info{idTrend(1),3});
        [testTrend, idTest] = max(diffTrend);
        if testTrend>precision
            conv = false;
            dprintf('Variable %s (of type %s) growth factor is still wrong (diff is %s)', nameTrend{idTest}, type, testTrend)
        end
    end
end