Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 0 additions and 4023 deletions
dynare++/src/dynare_model.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2006-2011, Ondra Kamenik
#include "parser/cc/parser_exception.hh"
#include "parser/cc/location.hh"
#include "utils/cc/exception.hh"
#include "dynare_model.hh"
#include "dynare_exception.hh"
#include "planner_builder.hh"
#include "forw_subst_builder.hh"
#include <cstdlib>
#include <string>
#include <cmath>
#include <climits>
#include <ostream>
using namespace ogdyn;
ParsedMatrix::ParsedMatrix(const ogp::MatrixParser &mp)
: TwoDMatrix(mp.nrows(), mp.ncols())
{
zeros();
for (ogp::MPIterator it = mp.begin(); it != mp.end(); ++it)
get(it.row(), it.col()) = *it;
}
DynareModel::DynareModel()
: atoms(), eqs(atoms)
{
}
DynareModel::DynareModel(const DynareModel &dm)
: atoms(dm.atoms), eqs(dm.eqs, atoms), order(dm.order),
param_vals(nullptr), init_vals(nullptr), vcov_mat(nullptr),
t_plobjective(dm.t_plobjective),
t_pldiscount(dm.t_pldiscount),
pbuilder(nullptr), fbuilder(nullptr),
atom_substs(nullptr), old_atoms(nullptr)
{
if (dm.param_vals)
param_vals = new Vector((const Vector &) *(dm.param_vals));
if (dm.init_vals)
init_vals = new Vector((const Vector &) *(dm.init_vals));
if (dm.vcov_mat)
vcov_mat = new TwoDMatrix((const TwoDMatrix &) *(dm.vcov_mat));
if (dm.old_atoms)
old_atoms = new DynareDynamicAtoms((const DynareDynamicAtoms &) *(dm.old_atoms));
if (dm.atom_substs)
atom_substs = new ogp::AtomSubstitutions((*dm.atom_substs), *old_atoms, atoms);
if (dm.pbuilder)
pbuilder = new PlannerBuilder(*(dm.pbuilder), *this);
if (dm.fbuilder)
fbuilder = new ForwSubstBuilder(*(dm.fbuilder), *this);
}
DynareModel::~DynareModel()
{
if (param_vals)
delete param_vals;
if (init_vals)
delete init_vals;
if (vcov_mat)
delete vcov_mat;
if (old_atoms)
delete old_atoms;
if (atom_substs)
delete atom_substs;
if (pbuilder)
delete pbuilder;
if (fbuilder)
delete fbuilder;
}
const PlannerInfo *
DynareModel::get_planner_info() const
{
if (pbuilder)
return &(pbuilder->get_info());
return nullptr;
}
const ForwSubstInfo *
DynareModel::get_forw_subst_info() const
{
if (fbuilder)
return &(fbuilder->get_info());
return nullptr;
}
const ogp::SubstInfo *
DynareModel::get_subst_info() const
{
if (atom_substs)
return &(atom_substs->get_info());
return nullptr;
}
void
DynareModel::setInitOuter(const Vector &x)
{
if (x.length() != atoms.ny())
throw DynareException(__FILE__, __LINE__,
"Wrong length of vector in DynareModel::setInitOuter");
for (int i = 0; i < atoms.ny(); i++)
(*init_vals)[i] = x[atoms.y2outer_endo()[i]];
}
void
DynareModel::print() const
{
printf("all atoms:\n");
atoms.print();
printf("formulas:\n");
DebugOperationFormatter dof(*this);
for (int i = 0; i < eqs.nformulas(); i++)
{
int tf = eqs.formula(i);
printf("formula %d:\n", tf);
eqs.getTree().print_operation_tree(tf, stdout, dof);
}
}
void
DynareModel::dump_model(std::ostream &os) const
{
// endogenous variable declaration
os << "var";
for (auto i : atoms.get_endovars())
os << " " << i;
os << ";\n\n";
// exogenous variables
os << "varexo";
for (auto i : atoms.get_exovars())
os << " " << i;
os << ";\n\n";
// parameters
os << "parameters";
for (auto i : atoms.get_params())
os << " " << i;
os << ";\n\n";
// parameter values
os.precision(16);
for (int i = 0; i < (int) atoms.get_params().size(); i++)
os << atoms.get_params()[i] << "=" << getParams()[i] << ";\n";
os << "\n\n";
// model section
ogp::OperationStringConvertor osc(atoms, getParser().getTree());
os << "model;\n";
for (int i = 0; i < getParser().nformulas(); i++)
{
os << "// Equation " << i << "\n0 = ";
int t = getParser().formula(i);
os << osc.convert(getParser().getTree().operation(t), t);
os << ";\n";
}
os << "end;\n";
// initval as steady state
os << "initval;\n";
for (int i = 0; i < (int) atoms.get_endovars().size(); i++)
os << atoms.get_endovars()[atoms.y2outer_endo()[i]] << "=" << getInit()[i] << ";\n";
os << "end;\n";
}
void
DynareModel::add_name(const char *name, int flag)
{
if (flag == 1)
{
// endogenous
atoms.register_uniq_endo(name);
}
else if (flag == 2)
{
// exogenous
atoms.register_uniq_exo(name);
}
else if (flag == 3)
{
// parameter
atoms.register_uniq_param(name);
}
else
{
throw DynareException(__FILE__, __LINE__,
"Unrecognized flag value.");
}
}
void
DynareModel::check_model() const
{
if (order == -1)
throw DynareException(__FILE__, __LINE__,
"Order of approximation not set in DynareModel::check_model");
if (atoms.ny() != eqs.nformulas())
{
char mes[1000];
sprintf(mes, "Model has %d equations for %d endogenous variables", eqs.nformulas(), atoms.ny());
throw DynareException(__FILE__, __LINE__, mes);
}
// check whether all nulary terms of all formulas in eqs are
// either constant or assigned to a name
for (int i = 0; i < eqs.nformulas(); i++)
{
int ft = eqs.formula(i);
const unordered_set<int> &nuls = eqs.nulary_of_term(ft);
for (int nul : nuls)
if (!atoms.is_constant(nul) && !atoms.is_named_atom(nul))
throw DynareException(__FILE__, __LINE__,
"Dangling nulary term found, internal error.");
}
int mlag, mlead;
atoms.exovarspan(mlead, mlag);
if (atoms.nexo() > 0 && (mlead != 0 || mlag != 0))
throw DynareException(__FILE__, __LINE__,
"The model contains occurrences of lagged/leaded exogenous variables");
atoms.endovarspan(mlead, mlag);
if (mlead > 1 || mlag < -1)
throw DynareException(__FILE__, __LINE__,
"The model contains occurrences of too lagged/leaded endogenous variables");
// check the dimension of vcov matrix
if (getAtoms().nexo() != getVcov().nrows())
throw DynareException(__FILE__, __LINE__,
"Dimension of VCOV matrix does not correspond to the shocks");
}
int
DynareModel::variable_shift(int t, int tshift)
{
const char *name = atoms.name(t);
if (atoms.is_type(name, DynareDynamicAtoms::param)
|| atoms.is_constant(t))
throw DynareException(__FILE__, __LINE__,
"The tree index is not a variable in DynareModel::variable_shift");
int ll = atoms.lead(t) + tshift;
int res = atoms.index(name, ll);
if (res == -1)
{
std::string str(name);
str += '(';
char tmp[50];
sprintf(tmp, "%d", ll);
str += tmp;
str += ')';
res = eqs.add_nulary(str.c_str());
}
return res;
}
void
DynareModel::variable_shift_map(const unordered_set<int> &a_set, int tshift,
map<int, int> &s_map)
{
s_map.clear();
for (int t : a_set)
{
// make shift map only for non-constants and non-parameters
if (!atoms.is_constant(t))
{
const char *name = atoms.name(t);
if (atoms.is_type(name, DynareDynamicAtoms::endovar)
|| atoms.is_type(name, DynareDynamicAtoms::exovar))
{
int tt = variable_shift(t, tshift);
s_map.insert(map<int, int>::value_type(t, tt));
}
}
}
}
void
DynareModel::termspan(int t, int &mlead, int &mlag) const
{
mlead = INT_MIN;
mlag = INT_MAX;
const unordered_set<int> &nul_terms = eqs.nulary_of_term(t);
for (int nul_term : nul_terms)
{
if (!atoms.is_constant(nul_term)
&& (atoms.is_type(atoms.name(nul_term), DynareDynamicAtoms::endovar)
|| atoms.is_type(atoms.name(nul_term), DynareDynamicAtoms::exovar)))
{
int ll = atoms.lead(nul_term);
if (ll < mlag)
mlag = ll;
if (ll > mlead)
mlead = ll;
}
}
}
bool
DynareModel::is_constant_term(int t) const
{
const unordered_set<int> &nul_terms = eqs.nulary_of_term(t);
for (int nul_term : nul_terms)
if (!atoms.is_constant(nul_term)
&& !atoms.is_type(atoms.name(nul_term), DynareDynamicAtoms::param))
return false;
return true;
}
unordered_set<int>
DynareModel::get_nonlinear_subterms(int t) const
{
NLSelector nls(*this);
return eqs.getTree().select_terms(t, nls);
}
void
DynareModel::substitute_atom_for_term(const char *name, int ll, int t)
{
// if the term t is itself a named atom (parameter, exo, endo),
// then we have to unassign it first
if (atoms.is_named_atom(t))
atoms.unassign_variable(atoms.name(t), atoms.lead(t), t);
// assign allocated tree index
// for the term now to name(ll)
atoms.assign_variable(name, ll, t);
// make operation t nulary in operation tree
eqs.nularify(t);
}
void
DynareModel::final_job()
{
if (t_plobjective != -1 && t_pldiscount != -1)
{
// at this moment include all equations and all variables; in
// future we will exclude purely exogenous processes; todo:
PlannerBuilder::Tvarset vset;
for (int i = 0; i < atoms.ny(); i++)
vset.insert(atoms.get_endovars()[i]);
PlannerBuilder::Teqset eset;
for (int i = 0; i < eqs.nformulas(); i++)
eset.push_back(i);
// construct the planner builder, this adds a lot of stuff to
// the model
if (pbuilder)
delete pbuilder;
pbuilder = new PlannerBuilder(*this, vset, eset);
}
// construct ForwSubstBuilder
if (fbuilder)
delete fbuilder;
fbuilder = new ForwSubstBuilder(*this);
// call parsing_finished (this will define an outer ordering of all variables)
atoms.parsing_finished(ogp::VarOrdering::bfspbfpb);
// make a copy of atoms and name it old_atoms
if (old_atoms)
delete old_atoms;
old_atoms = new DynareDynamicAtoms(atoms);
// construct empty substitutions from old_atoms to atoms
if (atom_substs)
delete atom_substs;
atom_substs = new ogp::AtomSubstitutions(*old_atoms, atoms);
// do the actual substitution, it will also call
// parsing_finished for atoms which creates internal orderings
atoms.substituteAllLagsAndExo1Leads(eqs, *atom_substs);
}
extern ogp::location_type dynglob_lloc;
DynareParser::DynareParser(const char *stream, int len, int ord)
: DynareModel(),
pa_atoms(), paramset(pa_atoms),
ia_atoms(), initval(ia_atoms), vcov(),
model_beg(0), model_end(-1),
paramset_beg(0), paramset_end(-1),
initval_beg(0), initval_end(-1),
vcov_beg(0), vcov_end(-1),
order_beg(0), order_end(-1),
plobjective_beg(0), plobjective_end(-1),
pldiscount_beg(0), pldiscount_end(-1)
{
// global parse
try
{
parse_glob(len, stream);
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, dynglob_lloc.off);
}
// setting parameters parse
try
{
if (paramset_end > paramset_beg)
paramset.parse(paramset_end-paramset_beg, stream+paramset_beg);
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, paramset_beg);
}
// model parse
try
{
if (model_end > model_beg)
eqs.parse(model_end-model_beg, stream+model_beg);
else
throw ogp::ParserException("Model section not found.", 0);
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, model_beg);
}
// initval setting parse
try
{
if (initval_end > initval_beg)
initval.parse(initval_end-initval_beg, stream+initval_beg);
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, initval_beg);
}
// vcov parse
try
{
if (vcov_end > vcov_beg)
{
vcov.parse(vcov_end-vcov_beg, stream+vcov_beg);
}
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, vcov_beg);
}
// planner objective parse
try
{
if (plobjective_end > plobjective_beg)
{
eqs.parse(plobjective_end-plobjective_beg, stream+plobjective_beg);
t_plobjective = eqs.pop_last_formula();
}
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, plobjective_beg);
}
// planner discount parse
try
{
if (pldiscount_end > pldiscount_beg)
{
t_pldiscount = parse_pldiscount(pldiscount_end - pldiscount_beg,
stream + pldiscount_beg);
}
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, pldiscount_beg);
}
// order parse
try
{
if (order_end > order_beg)
{
order = parse_order(order_end > order_beg, stream + order_beg);
}
}
catch (const ogp::ParserException &e)
{
throw ogp::ParserException(e, order_beg);
}
// check the overridden order
if (ord != -1)
order = ord;
// end parsing job, add planner's FOCs, make substitutions
DynareModel::final_job();
// calculate parameters
calc_params();
// calculate initial values
calc_init();
if (vcov_end > vcov_beg)
vcov_mat = new ParsedMatrix(vcov);
else
{
// vcov has not been asserted, set it to unit matrix
vcov_mat = new TwoDMatrix(atoms.nexo(), atoms.nexo());
vcov_mat->unit();
}
// check the model
check_model();
// differentiate
if (order >= 1)
eqs.differentiate(order);
}
DynareParser::DynareParser(const DynareParser &dp)
: DynareModel(dp),
pa_atoms(dp.pa_atoms), paramset(dp.paramset, pa_atoms),
ia_atoms(dp.ia_atoms), initval(dp.initval, ia_atoms), vcov(dp.vcov),
model_beg(dp.model_beg), model_end(dp.model_end),
paramset_beg(dp.paramset_beg), paramset_end(dp.paramset_end),
initval_beg(dp.initval_beg), initval_end(dp.initval_end),
vcov_beg(dp.vcov_beg), vcov_end(dp.vcov_end),
order_beg(dp.order_beg), order_end(dp.order_end),
plobjective_beg(dp.plobjective_beg), plobjective_end(dp.plobjective_end),
pldiscount_beg(dp.pldiscount_beg), pldiscount_end(dp.pldiscount_end)
{
}
DynareParser::~DynareParser()
= default;
void
DynareParser::add_name(const char *name, int flag)
{
DynareModel::add_name(name, flag);
// register with static atoms used for atom assignements
if (flag == 1)
{
// endogenous
ia_atoms.register_name(name);
}
else if (flag == 2)
{
// exogenous
ia_atoms.register_name(name);
}
else if (flag == 3)
{
// parameter
pa_atoms.register_name(name);
ia_atoms.register_name(name);
}
else
{
throw DynareException(__FILE__, __LINE__,
"Unrecognized flag value.");
}
}
void
DynareParser::error(const char *mes)
{
// throwing zero offset since this exception will be caugth at
// constructor
throw ogp::ParserException(mes, 0);
}
void
DynareParser::print() const
{
DynareModel::print();
printf("parameter atoms:\n");
paramset.print();
printf("initval atoms:\n");
initval.print();
printf("model position: %d %d\n", model_beg, model_end);
printf("paramset position: %d %d\n", paramset_beg, paramset_end);
printf("initval position: %d %d\n", initval_beg, initval_end);
}
/** A global symbol for passing info to the DynareParser from
* parser. */
DynareParser *dynare_parser;
/** The declarations of functions defined in dynglob_ll.cc and
* dynglob_tab.cc generated from dynglob.lex and dynglob.y */
void *dynglob__scan_buffer(char *, size_t);
void dynglob__destroy_buffer(void *);
void dynglob_parse();
extern ogp::location_type dynglob_lloc;
void
DynareParser::parse_glob(int length, const char *stream)
{
auto *buffer = new char[length+2];
strncpy(buffer, stream, length);
buffer[length] = '\0';
buffer[length+1] = '\0';
void *p = dynglob__scan_buffer(buffer, (unsigned int) length+2);
dynare_parser = this;
dynglob_parse();
delete [] buffer;
dynglob__destroy_buffer(p);
}
int
DynareParser::parse_order(int len, const char *str)
{
auto *buf = new char[len+1];
strncpy(buf, str, len);
buf[len] = '\0';
int res;
sscanf(buf, "%d", &res);
delete [] buf;
return res;
}
int
DynareParser::parse_pldiscount(int len, const char *str)
{
auto *buf = new char[len+1];
strncpy(buf, str, len);
buf[len] = '\0';
if (!atoms.is_type(buf, DynareDynamicAtoms::param))
throw ogp::ParserException(std::string("Name ") + buf + " is not a parameter", 0);
int t = atoms.index(buf, 0);
if (t == -1)
t = eqs.add_nulary(buf);
delete [] buf;
return t;
}
void
DynareParser::calc_params()
{
if (param_vals)
delete param_vals;
param_vals = new Vector(atoms.np());
ogp::AtomAsgnEvaluator aae(paramset);
aae.eval();
for (int i = 0; i < atoms.np(); i++)
(*param_vals)[i] = aae.get_value(atoms.get_params()[i]);
for (unsigned int i = 0; i < atoms.get_params().size(); i++)
if (!std::isfinite((*param_vals)[i]))
printf("dynare++: warning: value for parameter %s is not finite\n",
atoms.get_params()[i]);
}
void
DynareParser::calc_init()
{
// update initval atoms assignings according to substitutions
if (atom_substs)
initval.apply_subst(atom_substs->get_old2new());
// calculate the vector of initial values
if (init_vals)
delete init_vals;
init_vals = new Vector(atoms.ny());
ogp::AtomAsgnEvaluator aae(initval);
// set parameters
for (int ip = 0; ip < atoms.np(); ip++)
aae.set_user_value(atoms.get_params()[ip], (*param_vals)[ip]);
// set exogenous to zeros
for (int ie = 0; ie < atoms.nexo(); ie++)
aae.set_user_value(atoms.get_exovars()[ie], 0.0);
// evaluate
aae.eval();
// set results to internally ordered vector init_vals
for (int outer = 0; outer < atoms.ny(); outer++)
{
int i = atoms.outer2y_endo()[outer];
(*init_vals)[i] = aae.get_value(atoms.get_endovars()[outer]);
}
// if the planner's FOCs have been added, then add estimate of
// Lagrange multipliers to the vector
if (pbuilder)
{
MultInitSS mis(*pbuilder, *param_vals, *init_vals);
}
// if forward substitution builder has been created, we have to
// its substitutions and evaluate them
if (fbuilder)
ogdyn::DynareSteadySubstitutions dss(atoms, eqs.getTree(),
fbuilder->get_aux_map(), *param_vals, *init_vals);
for (unsigned int i = 0; i < atoms.get_endovars().size(); i++)
if (!std::isfinite((*init_vals)[i]))
printf("dynare++: warning: initval for <%s> is not finite\n",
atoms.get_endovars()[atoms.y2outer_endo()[i]]);
}
// this returns false for linear functions
bool
NLSelector::operator()(int t) const
{
const ogp::Operation &op = model.getParser().getTree().operation(t);
const DynareDynamicAtoms &atoms = model.getAtoms();
// if the term is constant, return false
if (model.is_constant_term(t))
return false;
int nary = op.nary();
if (nary == 0)
{
if (atoms.is_type(atoms.name(t), DynareDynamicAtoms::endovar)
|| atoms.is_type(atoms.name(t), DynareDynamicAtoms::exovar))
return true;
else
return false;
}
else if (nary == 1)
{
if (op.getCode() == ogp::UMINUS)
return false;
else
return true;
}
else
{
if (op.getCode() == ogp::TIMES)
// if at least one operand is constant, than the TIMES is linear
if (model.is_constant_term(op.getOp1())
|| model.is_constant_term(op.getOp2()))
return false;
else
return true;
// both PLUS and MINUS are linear
if (op.getCode() == ogp::PLUS
|| op.getCode() == ogp::MINUS)
return false;
// POWER is linear if exponent or base is 0 or one
if (op.getCode() == ogp::POWER
&& (op.getOp1() == ogp::OperationTree::zero
|| op.getOp1() == ogp::OperationTree::one
|| op.getOp2() == ogp::OperationTree::zero
|| op.getOp2() == ogp::OperationTree::one))
return false;
else
return true;
// DIVIDE is linear if the denominator is constant, or if
// the nominator is zero
if (op.getCode() == ogp::DIVIDE
&& (op.getOp1() == ogp::OperationTree::zero
|| model.is_constant_term(op.getOp2())))
return false;
else
return true;
}
throw DynareException(__FILE__, __LINE__,
"Wrong operation in operation tree");
return false;
}
DynareSPModel::DynareSPModel(const char **endo, int num_endo,
const char **exo, int num_exo,
const char **par, int num_par,
const char *equations, int len,
int ord)
: DynareModel()
{
// set the order
order = ord;
// add names
for (int i = 0; i < num_endo; i++)
add_name(endo[i], 1);
for (int i = 0; i < num_exo; i++)
add_name(exo[i], 2);
for (int i = 0; i < num_par; i++)
add_name(par[i], 3);
// parse the equations
eqs.parse(len, equations);
// parsing finished
atoms.parsing_finished(ogp::VarOrdering::bfspbfpb);
// create what has to be created from DynareModel
param_vals = new Vector(atoms.np());
init_vals = new Vector(atoms.ny());
vcov_mat = new TwoDMatrix(atoms.nexo(), atoms.nexo());
// check the model
check_model();
// differentiate
if (order >= 1)
eqs.differentiate(order);
}
void
ModelSSWriter::write_der0(FILE *fd)
{
write_der0_preamble(fd);
write_atom_assignment(fd);
stop_set.clear();
for (int fi = 0; fi < model.eqs.nformulas(); fi++)
otree.print_operation_tree(model.eqs.formula(fi), fd, *this);
write_der0_assignment(fd);
}
void
ModelSSWriter::write_der1(FILE *fd)
{
write_der1_preamble(fd);
write_atom_assignment(fd);
stop_set.clear();
const vector<int> &variables = model.getAtoms().variables();
const vector<int> &eam = model.getAtoms().get_endo_atoms_map();
for (int i = 0; i < model.getParser().nformulas(); i++)
{
const ogp::FormulaDerivatives &fder = model.getParser().derivatives(i);
for (int j : eam)
{
int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, j));
if (t > 0)
otree.print_operation_tree(t, fd, *this);
}
}
write_der1_assignment(fd);
}
MatlabSSWriter::MatlabSSWriter(const DynareModel &dm, const char *idd)
: ModelSSWriter(dm), id(new char[strlen(idd)+1])
{
strcpy(id, idd);
}
void
MatlabSSWriter::write_der0_preamble(FILE *fd) const
{
fprintf(fd,
"%% Usage:\n"
"%% out = %s_f(params, y)\n"
"%% where\n"
"%% out is a (%d,1) column vector of the residuals\n"
"%% of the static system\n",
id, model.getAtoms().ny());
write_common1_preamble(fd);
fprintf(fd,
"function out = %s_f(params, y)\n", id);
write_common2_preamble(fd);
}
void
MatlabSSWriter::write_der1_preamble(FILE *fd) const
{
fprintf(fd,
"%% Usage:\n"
"%% out = %s_ff(params, y)\n"
"%% where\n"
"%% out is a (%d,%d) matrix of the first order\n"
"%% derivatives of the static system residuals\n"
"%% columns correspond to endo variables in\n"
"%% the ordering as declared\n",
id, model.getAtoms().ny(), model.getAtoms().ny());
write_common1_preamble(fd);
fprintf(fd,
"function out = %s_ff(params, y)\n", id);
write_common2_preamble(fd);
}
void
MatlabSSWriter::write_common1_preamble(FILE *fd) const
{
fprintf(fd,
"%% params is a (%d,1) vector of parameter values\n"
"%% in the ordering as declared\n"
"%% y is a (%d,1) vector of endogenous variables\n"
"%% in the ordering as declared\n"
"%%\n"
"%% Created by Dynare++ v. %s\n", model.getAtoms().np(),
model.getAtoms().ny(), DYNVERSION);
// write ordering of parameters
fprintf(fd, "\n%% params ordering\n%% =====================\n");
for (auto parname : model.getAtoms().get_params())
{
fprintf(fd, "%% %s\n", parname);
}
// write endogenous variables
fprintf(fd, "%%\n%% y ordering\n%% =====================\n");
for (auto endoname : model.getAtoms().get_endovars())
{
fprintf(fd, "%% %s\n", endoname);
}
fprintf(fd, "\n");
}
void
MatlabSSWriter::write_common2_preamble(FILE *fd) const
{
fprintf(fd, "if size(y) ~= [%d,1]\n\terror('Wrong size of y, must be [%d,1]');\nend\n",
model.getAtoms().ny(), model.getAtoms().ny());
fprintf(fd, "if size(params) ~= [%d,1]\n\terror('Wrong size of params, must be [%d,1]');\nend\n\n",
model.getAtoms().np(), model.getAtoms().np());
}
void
MatlabSSWriter::write_atom_assignment(FILE *fd) const
{
// write OperationTree::num_constants
fprintf(fd, "%% hardwired constants\n");
ogp::EvalTree etree(model.getParser().getTree(), ogp::OperationTree::num_constants-1);
for (int i = 0; i < ogp::OperationTree::num_constants; i++)
{
format_nulary(i, fd);
double g = etree.eval(i);
if (std::isnan(g))
fprintf(fd, " = NaN;\n");
else
fprintf(fd, " = %12.8g;\n", etree.eval(i));
}
// write numerical constants
fprintf(fd, "%% numerical constants\n");
const ogp::Constants::Tconstantmap &cmap = model.getAtoms().get_constantmap();
for (auto it : cmap)
{
format_nulary(it.first, fd);
fprintf(fd, " = %12.8g;\n", it.second);
}
// write parameters
fprintf(fd, "%% parameter values\n");
for (unsigned int ip = 0; ip < model.getAtoms().get_params().size(); ip++)
{
const char *parname = model.getAtoms().get_params()[ip];
int t = model.getAtoms().index(parname, 0);
if (t == -1)
{
fprintf(fd, "%% %s not used in the model\n", parname);
}
else
{
format_nulary(t, fd);
fprintf(fd, " = params(%d); %% %s\n", ip+1, parname);
}
}
// write exogenous variables
fprintf(fd, "%% exogenous variables to zeros\n");
for (unsigned int ie = 0; ie < model.getAtoms().get_exovars().size(); ie++)
{
const char *exoname = model.getAtoms().get_exovars()[ie];
try
{
const ogp::DynamicAtoms::Tlagmap &lmap = model.getAtoms().lagmap(exoname);
for (auto it : lmap)
{
format_nulary(it.second, fd);
fprintf(fd, " = 0.0; %% %s\n", exoname);
}
}
catch (const ogu::Exception &e)
{
// ignore the error of not found variable in the tree
}
}
// write endogenous variables
fprintf(fd, "%% endogenous variables to y\n");
for (unsigned int ie = 0; ie < model.getAtoms().get_endovars().size(); ie++)
{
const char *endoname = model.getAtoms().get_endovars()[ie];
const ogp::DynamicAtoms::Tlagmap &lmap = model.getAtoms().lagmap(endoname);
for (auto it : lmap)
{
format_nulary(it.second, fd);
fprintf(fd, " = y(%d); %% %s\n", ie+1, endoname);
}
}
fprintf(fd, "\n");
}
void
MatlabSSWriter::write_der0_assignment(FILE *fd) const
{
// initialize out variable
fprintf(fd, "%% setting the output variable\n");
fprintf(fd, "out = zeros(%d, 1);\n", model.getParser().nformulas());
// fill out with the terms
for (int i = 0; i < model.getParser().nformulas(); i++)
{
fprintf(fd, "out(%d) = ", i+1);
format_term(model.getParser().formula(i), fd);
fprintf(fd, ";\n");
}
}
void
MatlabSSWriter::write_der1_assignment(FILE *fd) const
{
// initialize out variable
fprintf(fd, "%% setting the output variable\n");
fprintf(fd, "out = zeros(%d, %d);\n", model.getParser().nformulas(), model.getAtoms().ny());
// fill out with the terms
const vector<int> &variables = model.getAtoms().variables();
const vector<int> &eam = model.getAtoms().get_endo_atoms_map();
for (int i = 0; i < model.getParser().nformulas(); i++)
{
const ogp::FormulaDerivatives &fder = model.getParser().derivatives(i);
for (int j : eam)
{
int tvar = variables[j];
const char *name = model.getAtoms().name(tvar);
int yi = model.getAtoms().name2outer_endo(name);
int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, j));
if (t != ogp::OperationTree::zero)
{
fprintf(fd, "out(%d,%d) = out(%d,%d) + ", i+1, yi+1, i+1, yi+1);
format_term(t, fd);
fprintf(fd, "; %% %s(%d)\n", name, model.getAtoms().lead(tvar));
}
}
}
}
void
MatlabSSWriter::format_term(int t, FILE *fd) const
{
fprintf(fd, "t%d", t);
}
void
MatlabSSWriter::format_nulary(int t, FILE *fd) const
{
fprintf(fd, "a%d", t);
}
void
DebugOperationFormatter::format_nulary(int t, FILE *fd) const
{
const DynareDynamicAtoms &a = model.getAtoms();
if (t == ogp::OperationTree::zero)
fprintf(fd, "0");
else if (t == ogp::OperationTree::one)
fprintf(fd, "1");
else if (t == ogp::OperationTree::nan)
fprintf(fd, "NaN");
else if (t == ogp::OperationTree::two_over_pi)
fprintf(fd, "2/sqrt(PI)");
else if (a.is_constant(t))
fprintf(fd, "%g", a.get_constant_value(t));
else
{
int ll = a.lead(t);
const char *name = a.name(t);
if (ll == 0)
fprintf(fd, "%s", name);
else
fprintf(fd, "%s(%d)", name, ll);
}
}
dynare++/src/dynare_model.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2005-2011, Ondra Kamenik
#ifndef OGDYN_DYNARE_MODEL
#define OGDYN_DYNARE_MODEL
#include "parser/cc/matrix_parser.hh"
#include "parser/cc/atom_assignings.hh"
#include "dynare_atoms.hh"
#include "twod_matrix.hh"
#include "Vector.hh"
#include "GeneralMatrix.hh"
#include <map>
#include <unordered_set>
namespace ogdyn
{
using std::unordered_set;
using std::map;
/** This represents an interval in a string by the pair of
* positions (including the first, excluding the second). A
* position is given by the line and the column within the line
* (both starting from 1). */
struct PosInterval
{
int fl;
int fc;
int ll;
int lc;
PosInterval()
= default;
PosInterval(int ifl, int ifc, int ill, int ilc)
: fl(ifl), fc(ifc), ll(ill), lc(ilc)
{
}
PosInterval &
operator=(const PosInterval &pi)
= default;
/** This returns the interval beginning and interval length
* within the given string. */
void translate(const char *beg, int len, const char * &ibeg, int &ilen) const;
/** Debug print. */
void
print() const
{
printf("fl=%d fc=%d ll=%d lc=%d\n", fl, fc, ll, lc);
}
};
/** This class is basically a GeneralMatrix but is created from
* parsed matrix data. */
class ParsedMatrix : public TwoDMatrix
{
public:
/** Construct the object from the parsed data of ogp::MatrixParser. */
ParsedMatrix(const ogp::MatrixParser &mp);
};
class PlannerBuilder;
class PlannerInfo;
class ForwSubstBuilder;
class ForwSubstInfo;
class MultInitSS;
class ModelSSWriter;
/** A subclass is responsible for creating param_vals, init_vals,
* and vcov_mat. */
class DynareModel
{
friend class PlannerBuilder;
friend class ForwSubstBuilder;
friend class MultInitSS;
friend class ModelSSWriter;
protected:
/** All atoms for whole model. */
DynareDynamicAtoms atoms;
/** Parsed model equations. */
ogp::FormulaParser eqs;
/** Order of approximation. */
int order{-1};
/** A vector of parameters values created by a subclass. It
* is stored with natural ordering (outer) of the parameters
* given by atoms. */
Vector *param_vals{nullptr};
/** A vector of initial values created by a subclass. It is
* stored with internal ordering given by atoms. */
Vector *init_vals{nullptr};
/** A matrix for vcov. It is created by a subclass. */
TwoDMatrix *vcov_mat{nullptr};
/** Tree index of the planner objective. If there was no
* planner objective keyword, the value is set to -1. */
int t_plobjective{-1};
/** Tree index of the planner discount. If there was no
* planner discount keyword, the value is set to -1. */
int t_pldiscount{-1};
/** Pointer to PlannerBuilder, which is created only if the
* planner's FOC are added to the model. */
PlannerBuilder *pbuilder{nullptr};
/** Pointer to an object which builds auxiliary variables and
* equations to rewrite a model containing multiple leads to
* an equivalent model having only +1 leads. */
ForwSubstBuilder *fbuilder{nullptr};
/** Pointer to AtomSubstitutions which are created when the
* atoms are being substituted because of multiple lags
* etc. It uses also an old copy of atoms, which is
* created. */
ogp::AtomSubstitutions *atom_substs{nullptr};
/** Pointer to a copy of original atoms before substitutions
* took place. */
ogp::SAtoms *old_atoms{nullptr};
public:
/** Initializes the object to an empty state. */
DynareModel();
/** Construct a new deep copy. */
DynareModel(const DynareModel &dm);
virtual
~DynareModel();
virtual DynareModel *clone() const = 0;
const DynareDynamicAtoms &
getAtoms() const
{
return atoms;
}
const ogp::FormulaParser &
getParser() const
{
return eqs;
}
int
getOrder() const
{
return order;
}
/** Return the vector of parameter values. */
const Vector &
getParams() const
{
return *param_vals;
}
Vector &
getParams()
{
return *param_vals;
}
/** Return the vector of initial values of endo variables. */
const Vector &
getInit() const
{
return *init_vals;
}
Vector &
getInit()
{
return *init_vals;
}
/** Return the vcov matrix. */
const TwoDMatrix &
getVcov() const
{
return *vcov_mat;
}
TwoDMatrix &
getVcov()
{
return *vcov_mat;
}
/** Return planner info. */
const PlannerInfo *get_planner_info() const;
/** Return forward substitutions info. */
const ForwSubstInfo *get_forw_subst_info() const;
/** Return substitutions info. */
const ogp::SubstInfo *get_subst_info() const;
/** This sets initial values given in outer ordering. */
void setInitOuter(const Vector &x);
/** This returns true if the given term is a function of
* hardwired constants, numerical constants and parameters. */
bool is_constant_term(int t) const;
/** Debug print. */
void print() const;
/** Dump the model to the output stream. This includes
* variable declarations, parameter values, model code,
* initval, vcov and order. */
void dump_model(std::ostream &os) const;
protected:
/** Adds a name of endogenous, exogenous or a parameter. The
* sort is governed by the flag. See dynglob.y for values of
* the flag. This is used by a subclass when declaring the
* names. */
void add_name(const char *name, int flag);
/** This checks the model consistency. Thus includes: number
* of endo variables and number of equations, min and max lag
* of endogenous variables and occurrrences of exogenous
* variables. It throws an exception, if there is a problem. */
void check_model() const;
/** This shifts the given variable identified by the tree
* index in time. So if the given tree index represents a(+3)
* and the tshift is -4, the method returns tree index of the
* a(-1). If a(-1) doesn't exist, it is added to the tree. If
* it exists, its tree index is returned. If the tree index
* doesn't correspond to an endogenous nor exogenous variable,
* an exception is thrown. */
int variable_shift(int t, int tshift);
/** For the given set of atoms identified by tree indices and
* given time shift, this method returns a map mapping each
* variable in the given set to its time shifted variable. The
* map is passed through the reference and is cleared in the
* beginning. */
void variable_shift_map(const unordered_set<int> &a_set, int tshift,
map<int, int> &s_map);
/** This returns maximum lead and minimum lag of an endogenous
* or exogenous variable in the given term. If there are no
* endo or exo variables, than it returns the least integer as
* max lead and the greatest integer as min lag. */
void termspan(int t, int &mlead, int &mlag) const;
/** This function returns a set of non-linear subterms of the
* given term, these are terms whose linear combination
* constitutes the given term. */
unordered_set<int> get_nonlinear_subterms(int t) const;
/** This method assigns already used tree index of some term
* to the not-yet used atom name with the given lead/lag. In
* this way, all occurrences of term t are substituted with
* the atom name(ll). The method handles also rewriting
* operation tree including derivatives of the term t. */
void substitute_atom_for_term(const char *name, int ll, int t);
/** This performs a final job after the model is parsed. It
* creates the PlannerBuilder object if the planner's FOC are
* needed, then it creates ForwSubstBuilder handling multiple
* leads and finally it creates the substitution object saving
* old atoms and performs the substitutions. */
void final_job();
};
/** This class constructs DynareModel from dynare++ model file. It
* parses variable declarations, model equations, parameter
* assignments, initval assignments, vcov matrix and order of
* approximation. */
class DynareParser : public DynareModel
{
protected:
/** Static atoms for parameter assignments. */
DynareStaticAtoms pa_atoms;
/** Assignments for the parameters. */
ogp::AtomAssignings paramset;
/** Static atoms for initval assignments. */
DynareStaticAtoms ia_atoms;
/** Assignments for the initval. */
ogp::AtomAssignings initval;
/** Matrix parser for vcov. */
ogp::MatrixParser vcov;
public:
/** This, in fact, creates DynareModel from the given string
* of the given length corresponding to the Dynare++ model
* file. If the given ord is not -1, then it overrides setting
* in the model file. */
DynareParser(const char *str, int len, int ord);
DynareParser(const DynareParser &p);
~DynareParser() override;
DynareModel *
clone() const override
{
return new DynareParser(*this);
}
/** Adds a name of endogenous, exogenous or a parameter. This
* addss the name to the parent class DynareModel and also
* registers the name to either paramset, or initval. */
void add_name(const char *name, int flag);
/** Sets position of the model section. Called from
* dynglob.y. */
void
set_model_pos(int off1, int off2)
{
model_beg = off1; model_end = off2;
}
/** Sets position of the section setting parameters. Called
* from dynglob.y. */
void
set_paramset_pos(int off1, int off2)
{
paramset_beg = off1; paramset_end = off2;
}
/** Sets position of the initval section. Called from
* dynglob.y. */
void
set_initval_pos(int off1, int off2)
{
initval_beg = off1; initval_end = off2;
}
/** Sets position of the vcov section. Called from
* dynglob.y. */
void
set_vcov_pos(int off1, int off2)
{
vcov_beg = off1; vcov_end = off2;
}
/** Parser the given string as integer and set to as the
* order. */
void
set_order_pos(int off1, int off2)
{
order_beg = off1; order_end = off2;
}
/** Sets position of the planner_objective section. Called
* from dynglob.y. */
void
set_pl_objective_pos(int off1, int off2)
{
plobjective_beg = off1; plobjective_end = off2;
}
/** Sets position of the planner_discount section. Called from
* dynglob.y. */
void
set_pl_discount_pos(int off1, int off2)
{
pldiscount_beg = off1; pldiscount_end = off2;
}
/** Processes a syntax error from bison. */
void error(const char *mes);
/** Debug print. */
void print() const;
protected:
void parse_glob(int length, const char *stream);
int parse_order(int length, const char *stream);
int parse_pldiscount(int length, const char *stream);
/** Evaluate paramset assignings and set param_vals. */
void calc_params();
/** Evaluate initval assignings and set init_vals. */
void calc_init();
/** Do the final job. This includes building the planner
* problem (if any) and substituting for multiple lags, and
* one period leads of exogenous variables, and calculating
* initial guess of lagrange multipliers in the social planner
* problem. Precondtion: everything parsed and calculated
* parameters, postcondition: calculated initvals vector and
* parsing_finished for expanded vectors. */
void final_job();
private:
int model_beg, model_end;
int paramset_beg, paramset_end;
int initval_beg, initval_end;
int vcov_beg, vcov_end;
int order_beg, order_end;
int plobjective_beg, plobjective_end;
int pldiscount_beg, pldiscount_end;
};
/** Semiparsed model. The equations are given by a string,
* everything other by C/C++ objects. The initial values are set
* manually after the creation of this object. This implies that
* no automatic substitutions cannot be done here, which in turn
* implies that we cannot do here a social planner nor substitutions
* of multiple lags. */
class DynareSPModel : public DynareModel
{
public:
DynareSPModel(const char **endo, int num_endo,
const char **exo, int num_exo,
const char **par, int num_par,
const char *equations, int len, int ord);
DynareSPModel(const DynareSPModel &dm)
= default;
~DynareSPModel()
override = default;
DynareModel *
clone() const override
{
return new DynareSPModel(*this);
}
};
/** This class implements a selector of operations which correspond
* to non-linear functions. This inherits from ogp::opselector and
* is used to calculate non-linear subterms in
* DynareModel::get_nonlinear_subterms(). */
class NLSelector : public ogp::opselector
{
private:
const DynareModel &model;
public:
NLSelector(const DynareModel &m) : model(m)
{
}
bool operator()(int t) const override;
};
/** This class writes a mathematical code evaluating the system of
* equations and the first derivatives at zero shocks and at the
* given (static) state. Static means that lags and leads are
* ignored. */
class ModelSSWriter : public ogp::DefaultOperationFormatter
{
protected:
const DynareModel &model;
public:
ModelSSWriter(const DynareModel &m)
: DefaultOperationFormatter(m.eqs.getTree()),
model(m)
{
}
/** This writes the evaluation of the system. It calls pure
* virtual methods for writing a preamble, then assignment of
* atoms, and then assignment for resulting object. These are
* language dependent and are implemented in the subclass. */
void write_der0(FILE *fd);
/** This writes the evaluation of the first order derivative of
the system. It calls pure virtual methods for writing a
preamble, assignment, and assignemnt of the resulting
objects. */
void write_der1(FILE *fd);
protected:
virtual void write_der0_preamble(FILE *fd) const = 0;
virtual void write_der1_preamble(FILE *fd) const = 0;
virtual void write_atom_assignment(FILE *fd) const = 0;
virtual void write_der0_assignment(FILE *fd) const = 0;
virtual void write_der1_assignment(FILE *fd) const = 0;
};
class MatlabSSWriter : public ModelSSWriter
{
protected:
/** Identifier used in function names. */
char *id;
public:
MatlabSSWriter(const DynareModel &dm, const char *idd);
~MatlabSSWriter() override
{
delete [] id;
}
protected:
// from ModelSSWriter
void write_der0_preamble(FILE *fd) const override;
void write_der1_preamble(FILE *fd) const override;
/** This writes atom assignments. We have four kinds of atoms
* set here: endogenous vars coming from one parameter,
* parameter values given by the second parameter, constants,
* and the OperationTree::num_constants hardwired constants in
* ogp::OperationTree. */
void write_atom_assignment(FILE *fd) const override;
void write_der0_assignment(FILE *fd) const override;
void write_der1_assignment(FILE *fd) const override;
/** This prints t10 for t=10. */
void format_term(int t, FILE *fd) const override;
/** This prints a10 for t=10. The atoms a10 are supposed to be
* set by write_atom_assignments(). */
void format_nulary(int t, FILE *fd) const override;
private:
void write_common1_preamble(FILE *fd) const;
void write_common2_preamble(FILE *fd) const;
};
/** This class implements OperationFormatter for debugging
* purposes. It renders atoms in a more friendly way than the
* ogp::DefaulOperationFormatter. */
class DebugOperationFormatter : public ogp::DefaultOperationFormatter
{
protected:
const DynareModel &model;
public:
DebugOperationFormatter(const DynareModel &m)
: DefaultOperationFormatter(m.getParser().getTree()),
model(m)
{
}
void format_nulary(int t, FILE *fd) const override;
};
};
#endif
dynare++/src/dynare_params.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2004-2011, Ondra Kamenik
#include "dynare_params.hh"
#include <getopt.h>
#include <cstdio>
#include <cstring>
#include <thread>
const char *help_str
= "usage: dynare++ [--help] [--version] [options] <model file>\n"
"\n"
" --help print this message and return\n"
" --version print version and return\n"
"\n"
"options:\n"
" --per <num> number of periods simulated after burnt [100]\n"
" --burn <num> number of periods burnt [0]\n"
" --sim <num> number of simulations [80]\n"
" --rtper <num> number of RT periods simulated after burnt [0]\n"
" --rtsim <num> number of RT simulations [0]\n"
" --condper <num> number of periods in cond. simulations [0]\n"
" --condsim <num> number of conditional simulations [0]\n"
" --steps <num> steps towards stoch. SS [0=deter.]\n"
" --centralize centralize the rule [do centralize]\n"
" --no-centralize do not centralize the rule [do centralize]\n"
" --prefix <string> prefix of variables in Mat-4 file [\"dyn\"]\n"
" --seed <num> random number generator seed [934098]\n"
" --order <num> order of approximation [no default]\n"
" --threads <num> number of max parallel threads [nb. of logical CPUs]\n"
" --ss-tol <num> steady state calcs tolerance [1.e-13]\n"
" --check pesPES check model residuals [no checks]\n"
" lower/upper case switches off/on\n"
" pP checking along simulation path\n"
" eE checking on ellipse\n"
" sS checking along shocks\n"
" --check-evals <num> max number of evals per residual [1000]\n"
" --check-num <num> number of checked points [10]\n"
" --check-scale <num> scaling of checked points [2.0]\n"
" --no-irfs shuts down IRF simulations [do IRFs]\n"
" --irfs performs IRF simulations [do IRFs]\n"
" --qz-criterium <num> threshold for stable eigenvalues [1.000001]\n"
"\n\n";
// returns the pointer to the first character after the last slash or
// backslash in the string
const char *dyn_basename(const char *str);
DynareParams::DynareParams(int argc, char **argv)
: modname(nullptr), num_per(100), num_burn(0), num_sim(80),
num_rtper(0), num_rtsim(0),
num_condper(0), num_condsim(0),
num_threads(std::thread::hardware_concurrency()), num_steps(0),
prefix("dyn"), seed(934098), order(-1), ss_tol(1.e-13),
check_along_path(false), check_along_shocks(false),
check_on_ellipse(false), check_evals(1000), check_num(10), check_scale(2.0),
do_irfs_all(true), do_centralize(true), qz_criterium(1.0+1e-6),
help(false), version(false)
{
if (argc == 1 || !strcmp(argv[1], "--help"))
{
help = true;
return;
}
if (argc == 1 || !strcmp(argv[1], "--version"))
{
version = true;
return;
}
modname = argv[argc-1];
argc--;
struct option const opts [] = {
{"periods", required_argument, nullptr, opt_per},
{"per", required_argument, nullptr, opt_per},
{"burn", required_argument, nullptr, opt_burn},
{"simulations", required_argument, nullptr, opt_sim},
{"sim", required_argument, nullptr, opt_sim},
{"rtperiods", required_argument, nullptr, opt_rtper},
{"rtper", required_argument, nullptr, opt_rtper},
{"rtsimulations", required_argument, nullptr, opt_rtsim},
{"rtsim", required_argument, nullptr, opt_rtsim},
{"condperiods", required_argument, nullptr, opt_condper},
{"condper", required_argument, nullptr, opt_condper},
{"condsimulations", required_argument, nullptr, opt_condsim},
{"condsim", required_argument, nullptr, opt_condsim},
{"prefix", required_argument, nullptr, opt_prefix},
{"threads", required_argument, nullptr, opt_threads},
{"steps", required_argument, nullptr, opt_steps},
{"seed", required_argument, nullptr, opt_seed},
{"order", required_argument, nullptr, opt_order},
{"ss-tol", required_argument, nullptr, opt_ss_tol},
{"check", required_argument, nullptr, opt_check},
{"check-scale", required_argument, nullptr, opt_check_scale},
{"check-evals", required_argument, nullptr, opt_check_evals},
{"check-num", required_argument, nullptr, opt_check_num},
{"qz-criterium", required_argument, nullptr, opt_qz_criterium},
{"no-irfs", no_argument, nullptr, opt_noirfs},
{"irfs", no_argument, nullptr, opt_irfs},
{"centralize", no_argument, nullptr, opt_centralize},
{"no-centralize", no_argument, nullptr, opt_no_centralize},
{"help", no_argument, nullptr, opt_help},
{"version", no_argument, nullptr, opt_version},
{nullptr, 0, nullptr, 0}
};
int ret;
int index;
while (-1 != (ret = getopt_long(argc, argv, "", opts, &index)))
{
switch (ret)
{
case opt_per:
if (1 != sscanf(optarg, "%d", &num_per))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_burn:
if (1 != sscanf(optarg, "%d", &num_burn))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_sim:
if (1 != sscanf(optarg, "%d", &num_sim))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_rtper:
if (1 != sscanf(optarg, "%d", &num_rtper))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_rtsim:
if (1 != sscanf(optarg, "%d", &num_rtsim))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_condper:
if (1 != sscanf(optarg, "%d", &num_condper))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_condsim:
if (1 != sscanf(optarg, "%d", &num_condsim))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_prefix:
prefix = optarg;
break;
case opt_threads:
if (1 != sscanf(optarg, "%d", &num_threads))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_steps:
if (1 != sscanf(optarg, "%d", &num_steps))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_seed:
if (1 != sscanf(optarg, "%d", &seed))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_order:
if (1 != sscanf(optarg, "%d", &order))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_ss_tol:
if (1 != sscanf(optarg, "%lf", &ss_tol))
fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg);
break;
case opt_check:
processCheckFlags(optarg);
break;
case opt_check_scale:
if (1 != sscanf(optarg, "%lf", &check_scale))
fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg);
break;
case opt_check_evals:
if (1 != sscanf(optarg, "%d", &check_evals))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_check_num:
if (1 != sscanf(optarg, "%d", &check_num))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
case opt_noirfs:
irf_list.clear();
do_irfs_all = false;
break;
case opt_irfs:
processIRFList(argc, argv);
if (irf_list.empty())
do_irfs_all = true;
else
do_irfs_all = false;
break;
case opt_centralize:
do_centralize = true;
break;
case opt_no_centralize:
do_centralize = false;
break;
case opt_qz_criterium:
if (1 != sscanf(optarg, "%lf", &qz_criterium))
fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg);
break;
case opt_help:
help = true;
break;
case opt_version:
version = true;
break;
case '?':
fprintf(stderr, "Unknown option, ignored\n");
break;
}
}
// make basename (get rid of the extension)
basename = dyn_basename(modname);
std::string::size_type i = basename.rfind('.');
if (i != std::string::npos)
basename.erase(i);
}
void
DynareParams::printHelp() const
{
printf("%s", help_str);
}
void
DynareParams::processCheckFlags(const char *flags)
{
for (unsigned int i = 0; i < strlen(flags); i++)
{
switch (flags[i])
{
case 'p':
check_along_path = false;
break;
case 'P':
check_along_path = true;
break;
case 'e':
check_on_ellipse = false;
break;
case 'E':
check_on_ellipse = true;
break;
case 's':
check_along_shocks = false;
break;
case 'S':
check_along_shocks = true;
break;
default:
fprintf(stderr, "Unknown check type selection character <%c>, ignored.\n", flags[i]);
}
}
}
void
DynareParams::processIRFList(int argc, char **argv)
{
irf_list.clear();
while (optind < argc && *(argv[optind]) != '-')
{
irf_list.push_back(argv[optind]);
optind++;
}
}
const char *
dyn_basename(const char *str)
{
int i = strlen(str);
while (i > 0 && str[i-1] != '/' && str[i-1] != '\\')
i--;
return str+i;
}
dynare++/src/dynare_params.hh deleted 100644 → 0
View file @ 219d2bb3
// $Id: dynare_params.h 2347 2009-03-24 11:54:29Z kamenik $
// Copyright 2004, Ondra Kamenik
/*
along shocks: m mult max_evals
ellipse: m mult max_evals (10*m) (0.5*mult)
simul: m max_evals (10*m)
--check-scale 2.0 --check-evals 1000 --check-num 10 --check PES
*/
#include <vector>
#include <string>
struct DynareParams
{
const char *modname;
std::string basename;
int num_per;
int num_burn;
int num_sim;
int num_rtper;
int num_rtsim;
int num_condper;
int num_condsim;
int num_threads;
int num_steps;
const char *prefix;
int seed;
int order;
/** Tolerance used for steady state calcs. */
double ss_tol;
bool check_along_path;
bool check_along_shocks;
bool check_on_ellipse;
int check_evals;
int check_num;
double check_scale;
/** Flag for doing IRFs even if the irf_list is empty. */
bool do_irfs_all;
/** List of shocks for which IRF will be calculated. */
std::vector<const char *> irf_list;
bool do_centralize;
double qz_criterium;
bool help;
bool version;
DynareParams(int argc, char **argv);
void printHelp() const;
int
getCheckShockPoints() const
{
return check_num;
}
double
getCheckShockScale() const
{
return check_scale;
}
int
getCheckEllipsePoints() const
{
return 10*check_num;
}
double
getCheckEllipseScale() const
{
return 0.5*check_scale;
}
int
getCheckPathPoints() const
{
return 10*check_num;
}
private:
enum {opt_per, opt_burn, opt_sim, opt_rtper, opt_rtsim, opt_condper, opt_condsim,
opt_prefix, opt_threads,
opt_steps, opt_seed, opt_order, opt_ss_tol, opt_check,
opt_check_along_path, opt_check_along_shocks, opt_check_on_ellipse,
opt_check_evals, opt_check_scale, opt_check_num, opt_noirfs, opt_irfs,
opt_help, opt_version, opt_centralize, opt_no_centralize, opt_qz_criterium};
void processCheckFlags(const char *flags);
/** This gathers strings from argv[optind] and on not starting
* with '-' to the irf_list. It stops one item before the end,
* since this is the model file. */
void processIRFList(int argc, char **argv);
};
dynare++/src/dynglob.ll deleted 100644 → 0
View file @ 219d2bb3
/* -*- C++ -*- */
%{
#include "parser/cc/location.hh"
#include "dynglob_tab.hh"
extern YYLTYPE dynglob_lloc;
#define YY_USER_ACTION SET_LLOC(dynglob_);
%}
%option nounput
%option noyy_top_state
%option stack
%option yylineno
%option prefix="dynglob_"
%option never-interactive
%x CMT
%%
/* comments */
<*>"/*" {yy_push_state(CMT);}
<CMT>[^*\n]*
<CMT>"*"+[^*/\n]*
<CMT>"*"+"/" {yy_pop_state();}
<CMT>[\n]
"//".*\n
/* initial spaces or tabs are ignored */
[ \t\r\n\0]
var {return VAR;}
varexo {return VAREXO;}
parameters {return PARAMETERS;}
model {return MODEL;}
end {return END;}
initval {return INITVAL;}
order {return ORDER;}
vcov {return VCOV;}
planner_objective {return PLANNEROBJECTIVE;}
planner_discount {return PLANNERDISCOUNT;}
/* names */
[A-Za-z_][A-Za-z0-9_]* {
dynglob_lval.string = dynglob_text;
return NAME;
}
; {return SEMICOLON;}
, {return COMMA;}
= {return EQUAL_SIGN;}
\[ {return LEFT_BRACKET;}
\] {return RIGHT_BRACKET;}
. {
dynglob_lval.character = dynglob_text[0];
return CHARACTER;
}
%%
int dynglob_wrap()
{
return 1;
}
void dynglob__destroy_buffer(void* p)
{
dynglob__delete_buffer((YY_BUFFER_STATE)p);
}
dynare++/src/dynglob.yy deleted 100644 → 0
View file @ 219d2bb3
// -*- C++ -*-
%{
// Copyright (C) 2006-2011, Ondra Kamenik
#include "parser/cc/location.hh"
#include "dynare_model.hh"
#include "dynglob_tab.hh"
#include <stdio.h>
void dynglob_error(const char*);
int dynglob_lex(void);
extern int dynglob_lineno;
extern ogdyn::DynareParser* dynare_parser;
int symblist_flag;
// static void print_token_value1 (FILE *, int, YYSTYPE);
//#define YYPRINT(file, type, value) print_token_value1 (file, type, value)
%}
%union {
int integer;
char *string;
char character;
}
%token END INITVAL MODEL PARAMETERS VAR VAREXO SEMICOLON COMMA EQUAL_SIGN CHARACTER
%token VCOV LEFT_BRACKET RIGHT_BRACKET ORDER PLANNEROBJECTIVE PLANNERDISCOUNT
%token <string> NAME;
%name-prefix="dynglob_"
%locations
%error-verbose
%%
dynare_file : preamble paramset model rest {
dynare_parser->set_paramset_pos(@2.off, @3.off);}
| preamble model rest {
dynare_parser->set_paramset_pos(0, 0);}
| preamble paramset planner model rest {
dynare_parser->set_paramset_pos(@2.off, @3.off);}
;
preamble : preamble preamble_statement | preamble_statement;
preamble_statement : var | varexo | parameters;
var : VAR {symblist_flag=1;} symblist SEMICOLON;
varexo : VAREXO {symblist_flag=2;} symblist SEMICOLON;
parameters : PARAMETERS {symblist_flag=3;} symblist SEMICOLON;
symblist : symblist NAME {dynare_parser->add_name($2,symblist_flag);}
| symblist COMMA NAME {dynare_parser->add_name($3,symblist_flag);}
| NAME {dynare_parser->add_name($1,symblist_flag);}
;
paramset : recnameset;
recnameset : recnameset onenameset | onenameset;
onenameset : NAME EQUAL_SIGN material SEMICOLON;
material : material CHARACTER | material NAME | NAME | CHARACTER;
model : MODEL SEMICOLON equations END SEMICOLON {
dynare_parser->set_model_pos(@3.off, @4.off);
};
equations : equations equation | equation;
equation : material EQUAL_SIGN material SEMICOLON | material SEMICOLON;
rest : rest_statement | rest rest_statement;
rest_statement : initval | vcov | order | planner;
initval : INITVAL SEMICOLON recnameset END SEMICOLON {
dynare_parser->set_initval_pos(@3.off, @4.off);
};
vcov : VCOV EQUAL_SIGN LEFT_BRACKET m_material RIGHT_BRACKET SEMICOLON {
dynare_parser->set_vcov_pos(@4.off, @5.off);
};
m_material : m_material CHARACTER | m_material NAME | m_material SEMICOLON | m_material COMMA | CHARACTER | NAME | SEMICOLON | COMMA;
order : ORDER EQUAL_SIGN material SEMICOLON {
dynare_parser->set_order_pos(@3.off, @4.off);
};
planner : planner_objective planner_discount
| planner_discount planner_objective
;
planner_objective : PLANNEROBJECTIVE material SEMICOLON {
dynare_parser->set_pl_objective_pos(@2.off, @3.off);
};
planner_discount : PLANNERDISCOUNT NAME SEMICOLON {
dynare_parser->set_pl_discount_pos(@2.off, @3.off);
};
%%
void dynglob_error(const char* mes)
{
dynare_parser->error(mes);
}
/*
static void print_token_value1(FILE* file, int type, YYSTYPE value)
{
if (type == NAME)
fprintf(file, "%s", value.string);
if (type == CHARACTER)
fprintf(file, "%c", value.character);
}
*/
dynare++/src/forw_subst_builder.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2006-2011, Ondra Kamenik
#include "forw_subst_builder.hh"
using namespace ogdyn;
ForwSubstBuilder::ForwSubstBuilder(DynareModel &m)
: model(m)
{
info.num_new_terms -= model.getParser().getTree().get_num_op();
// go through all equations
int neq = model.eqs.nformulas();
for (int i = 0; i < neq; i++)
{
int ft = model.eqs.formula(i);
int mlead, mlag;
model.termspan(ft, mlead, mlag);
// if equation is too forward looking
if (mlead > 1)
{
info.num_affected_equations++;
// break it to non-linear terms
unordered_set<int> nlt = model.get_nonlinear_subterms(ft);
int j = 0; // indexes subterms
// and make substitutions for all these non-linear subterms
for (unordered_set<int>::const_iterator it = nlt.begin();
it != nlt.end(); ++it, ++j)
substitute_for_term(*it, i, j);
}
}
// unassign all variables with lead greater than 1
unassign_gt_1_leads();
// forget the derivatives in the tree because some variables could
// have been unassigned
model.eqs.getTree().forget_derivative_maps();
info.num_new_terms += model.getParser().getTree().get_num_op();
}
void
ForwSubstBuilder::substitute_for_term(int t, int i, int j)
{
int mlead, mlag;
model.termspan(t, mlead, mlag);
if (mlead > 1)
{
info.num_subst_terms++;
// Example for comments: let t = f(x(+4))
// first make lagsubst be substitution setting f(x(+4)) to f(x(+1))
// this is lag = -3 (1-mlead)
map<int, int> lagsubst;
unordered_set<int> nult = model.eqs.nulary_of_term(t);// make copy of nult!
model.variable_shift_map(nult, 1-mlead, lagsubst);
int lagt = model.eqs.add_substitution(t, lagsubst);
// now maxlead of lagt is +1
// add AUXLD_*_*_1 = f(x(+1)) to the model
char name[100];
sprintf(name, "AUXLD_%d_%d_%d", i, j, 1);
model.atoms.register_uniq_endo(name);
info.num_aux_variables++;
const char *ss = model.atoms.get_name_storage().query(name);
int auxt = model.eqs.add_nulary(name);
model.eqs.add_formula(model.eqs.add_binary(ogp::MINUS, auxt, lagt));
aux_map.insert(Tsubstmap::value_type(ss, lagt));
// now add variables and equations
// AUXLD_*_*_2 = AUXLD_*_*_1(+1) through
// AUXLD_*_*_{mlead-1} = AUXLD_*_*_{mlead-2}(+1)
for (int ll = 1; ll <= mlead-2; ll++)
{
// create AUXLD_*_*_{ll}(+1)
sprintf(name, "AUXLD_%d_%d_%d(+1)", i, j, ll);
int lastauxt_lead = model.eqs.add_nulary(name);
// create AUXLD_*_*{ll+1}
sprintf(name, "AUXLD_%d_%d_%d", i, j, ll+1);
model.atoms.register_uniq_endo(name);
info.num_aux_variables++;
ss = model.atoms.get_name_storage().query(name);
auxt = model.eqs.add_nulary(name);
// add AUXLD_*_*_{ll+1} = AUXLD_*_*_{ll}(+1)
model.eqs.add_formula(model.eqs.add_binary(ogp::MINUS, auxt, lastauxt_lead));
// add substitution to the map; todo: this
// works well because in the context where
// aux_map is used the timing doesn't matter,
// however, it is misleading, needs to be
// changed
aux_map.insert(Tsubstmap::value_type(ss, lagt));
}
// now we have to substitute AUXLD_*_*{mlead-1}(+1) for t
sprintf(name, "AUXLD_%d_%d_%d", i, j, mlead-1);
ss = model.atoms.get_name_storage().query(name);
model.substitute_atom_for_term(ss, +1, t);
}
}
void
ForwSubstBuilder::unassign_gt_1_leads(const char *name)
{
const char *ss = model.atoms.get_name_storage().query(name);
int mlead, mlag;
model.atoms.varspan(name, mlead, mlag);
for (int ll = 2; ll <= mlead; ll++)
{
int t = model.atoms.index(ss, ll);
if (t != -1)
model.atoms.unassign_variable(ss, ll, t);
}
}
void
ForwSubstBuilder::unassign_gt_1_leads()
{
const vector<const char *> &endovars = model.atoms.get_endovars();
for (auto endovar : endovars)
unassign_gt_1_leads(endovar);
const vector<const char *> &exovars = model.atoms.get_exovars();
for (auto exovar : exovars)
unassign_gt_1_leads(exovar);
}
ForwSubstBuilder::ForwSubstBuilder(const ForwSubstBuilder &b, DynareModel &m)
: model(m)
{
for (auto it : b.aux_map)
{
const char *ss = m.atoms.get_name_storage().query(it.first);
aux_map.insert(Tsubstmap::value_type(ss, it.second));
}
}
dynare++/src/forw_subst_builder.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2006, Ondra Kamenik
// $Id$
#ifndef FORW_SUBST_BUILDER_H
#define FORW_SUBST_BUILDER_H
#include "dynare_model.hh"
namespace ogdyn
{
/** This struct encapsulates information about the process of
* forward substitutions. */
struct ForwSubstInfo
{
int num_affected_equations{0};
int num_subst_terms{0};
int num_aux_variables{0};
int num_new_terms{0};
ForwSubstInfo()
= default;
};
class ForwSubstBuilder
{
using Ttermauxmap = map<int, const char *>;
protected:
/** Reference to the model, to which we will add equations and
* change some equations. */
DynareModel &model;
/** A map mapping new auxiliary variables to the terms in the
* tree in the DynareModel. */
Tsubstmap aux_map;
/** Information about the substitutions. */
ForwSubstInfo info;
public:
/** Do all the jobs needed. This scans all equations in the
* model, and for equations containing forward looking
* variables greater than 1 lead, it makes corresponding
* substitutions. Basically, it breaks each equation to its
* non-linear components and creates substitutions for these
* components, not for whole equation. This is because the
* expectation operator can go through the linear part of the
* function. This will save us many occurrences of other
* variables involved in the equation. */
ForwSubstBuilder(DynareModel &m);
ForwSubstBuilder(const ForwSubstBuilder &b) = delete;
/** Copy constructor with a new instance of the model. */
ForwSubstBuilder(const ForwSubstBuilder &b, DynareModel &m);
/** Return the auxiliary variable mapping. */
const Tsubstmap &
get_aux_map() const
{
return aux_map;
}
/** Return the information. */
const ForwSubstInfo &
get_info() const
{
return info;
}
private:
/** This method takes a nonlinear term t, and if it has leads
* of greater than 1, then it substitutes the term for the new
* variable (or string of variables). Note that the
* substitution is done by DynamicAtoms::assign_variable. This
* means that the substitution is made for all other
* ocurrences of t in the model. So there is no need of
* tracking already substituted terms. The other two
* parameters are just for identification of the new auxiliary
* variables. When called from the constructor, i is an
* equation number, j is an order of the non-linear term in
* the equation. */
void substitute_for_term(int t, int i, int j);
/** This is called just at the end of the job. It unassigns
* all nulary terms with a lead greater than 1. */
void unassign_gt_1_leads();
/** This unassigns all leads greater than 1 of the given name. */
void unassign_gt_1_leads(const char *name);
};
};
#endif
dynare++/src/main.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2004-2011, Ondra Kamenik
#include "dynare3.hh"
#include "dynare_exception.hh"
#include "dynare_params.hh"
#include "utils/cc/exception.hh"
#include "parser/cc/parser_exception.hh"
#include "../sylv/cc/SylvException.hh"
#include "../kord/random.hh"
#include "../kord/global_check.hh"
#include "../kord/approximation.hh"
int
main(int argc, char **argv)
{
DynareParams params(argc, argv);
if (params.help)
{
params.printHelp();
return 0;
}
if (params.version)
{
printf("Dynare++ v. %s. Copyright (C) 2004-2011, Ondra Kamenik\n",
DYNVERSION);
printf("Dynare++ comes with ABSOLUTELY NO WARRANTY and is distributed under\n");
printf("GPL: modules integ, tl, kord, sylv, src, extern and documentation\n");
printf("LGPL: modules parser, utils\n");
printf(" for GPL see http://www.gnu.org/licenses/gpl.html\n");
printf(" for LGPL see http://www.gnu.org/licenses/lgpl.html\n");
return 0;
}
sthread::detach_thread_group::max_parallel_threads = params.num_threads;
try
{
// make journal name and journal
std::string jname(params.basename);
jname += ".jnl";
Journal journal(jname.c_str());
// make dynare object
Dynare dynare(params.modname, params.order, params.ss_tol, journal);
// make list of shocks for which we will do IRFs
std::vector<int> irf_list_ind;
if (params.do_irfs_all)
for (int i = 0; i < dynare.nexog(); i++)
irf_list_ind.push_back(i);
else
irf_list_ind = ((const DynareNameList &) dynare.getExogNames()).selectIndices(params.irf_list);
// write matlab files
FILE *mfd;
std::string mfile1(params.basename);
mfile1 += "_f.m";
if (nullptr == (mfd = fopen(mfile1.c_str(), "w")))
{
fprintf(stderr, "Couldn't open %s for writing.\n", mfile1.c_str());
exit(1);
}
ogdyn::MatlabSSWriter writer0(dynare.getModel(), params.basename.c_str());
writer0.write_der0(mfd);
fclose(mfd);
std::string mfile2(params.basename);
mfile2 += "_ff.m";
if (nullptr == (mfd = fopen(mfile2.c_str(), "w")))
{
fprintf(stderr, "Couldn't open %s for writing.\n", mfile2.c_str());
exit(1);
}
ogdyn::MatlabSSWriter writer1(dynare.getModel(), params.basename.c_str());
writer1.write_der1(mfd);
fclose(mfd);
// open mat file
std::string matfile(params.basename);
matfile += ".mat";
mat_t *matfd = Mat_Create(matfile.c_str(), nullptr);
if (matfd == nullptr)
{
fprintf(stderr, "Couldn't open %s for writing.\n", matfile.c_str());
exit(1);
}
// write info about the model (dimensions and variables)
dynare.writeMat(matfd, params.prefix);
// write the dump file corresponding to the input
dynare.writeDump(params.basename);
system_random_generator.initSeed(params.seed);
tls.init(dynare.order(),
dynare.nstat()+2*dynare.npred()+3*dynare.nboth()
+2*dynare.nforw()+dynare.nexog());
Approximation app(dynare, journal, params.num_steps, params.do_centralize, params.qz_criterium);
try
{
app.walkStochSteady();
}
catch (const KordException &e)
{
// tell about the exception and continue
printf("Caught (not yet fatal) Kord exception: ");
e.print();
JournalRecord rec(journal);
rec << "Solution routine not finished (" << e.get_message()
<< "), see what happens" << endrec;
}
std::string ss_matrix_name(params.prefix);
ss_matrix_name += "_steady_states";
ConstTwoDMatrix(app.getSS()).writeMat(matfd, ss_matrix_name.c_str());
// check the approximation
if (params.check_along_path || params.check_along_shocks
|| params.check_on_ellipse)
{
GlobalChecker gcheck(app, sthread::detach_thread_group::max_parallel_threads, journal);
if (params.check_along_shocks)
gcheck.checkAlongShocksAndSave(matfd, params.prefix,
params.getCheckShockPoints(),
params.getCheckShockScale(),
params.check_evals);
if (params.check_on_ellipse)
gcheck.checkOnEllipseAndSave(matfd, params.prefix,
params.getCheckEllipsePoints(),
params.getCheckEllipseScale(),
params.check_evals);
if (params.check_along_path)
gcheck.checkAlongSimulationAndSave(matfd, params.prefix,
params.getCheckPathPoints(),
params.check_evals);
}
// write the folded decision rule to the Mat-4 file
app.getFoldDecisionRule().writeMat(matfd, params.prefix);
// simulate conditional
if (params.num_condper > 0 && params.num_condsim > 0)
{
SimResultsDynamicStats rescond(dynare.numeq(), params.num_condper, 0);
Vector det_ss{app.getSS().getCol(0)};
rescond.simulate(params.num_condsim, app.getFoldDecisionRule(), det_ss, dynare.getVcov(), journal);
rescond.writeMat(matfd, params.prefix);
}
// simulate unconditional
//const DecisionRule& dr = app.getUnfoldDecisionRule();
const DecisionRule &dr = app.getFoldDecisionRule();
if (params.num_per > 0 && params.num_sim > 0)
{
SimResultsStats res(dynare.numeq(), params.num_per, params.num_burn);
res.simulate(params.num_sim, dr, dynare.getSteady(), dynare.getVcov(), journal);
res.writeMat(matfd, params.prefix);
// impulse response functions
if (!irf_list_ind.empty())
{
IRFResults irf(dynare, dr, res, irf_list_ind, journal);
irf.writeMat(matfd, params.prefix);
}
}
// simulate with real-time statistics
if (params.num_rtper > 0 && params.num_rtsim > 0)
{
RTSimResultsStats rtres(dynare.numeq(), params.num_rtper, params.num_burn);
rtres.simulate(params.num_rtsim, dr, dynare.getSteady(), dynare.getVcov(), journal);
rtres.writeMat(matfd, params.prefix);
}
Mat_Close(matfd);
}
catch (const KordException &e)
{
printf("Caugth Kord exception: ");
e.print();
return e.code();
}
catch (const TLException &e)
{
printf("Caugth TL exception: ");
e.print();
return 255;
}
catch (SylvException &e)
{
printf("Caught Sylv exception: ");
e.printMessage();
return 255;
}
catch (const DynareException &e)
{
printf("Caught Dynare exception: %s\n", e.message());
return 255;
}
catch (const ogu::Exception &e)
{
printf("Caught ogu::Exception: ");
e.print();
return 255;
}
catch (const ogp::ParserException &e)
{
printf("Caught parser exception: %s\n", e.message());
return 255;
}
return 0;
}
dynare++/src/nlsolve.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2006, Ondra Kamenik
// $Id: nlsolve.cpp 762 2006-05-22 13:00:07Z kamenik $
#include "nlsolve.hh"
#include "dynare_exception.hh"
#include <cmath>
using namespace ogu;
/** This should not be greater than DBL_EPSILON^(1/2). */
double GoldenSectionSearch::tol = 1.e-4;
/** This is equal to the golden section ratio. */
double GoldenSectionSearch::golden = (3.-std::sqrt(5.))/2;
double
GoldenSectionSearch::search(OneDFunction &f, double x1, double x2)
{
double b;
if (init_bracket(f, x1, x2, b))
{
double fb = f.eval(b);
double f1 = f.eval(x1);
f.eval(x2);
double dx;
do
{
double w = (b-x1)/(x2-x1);
dx = std::abs((1-2*w)*(x2-x1));
double x;
if (b-x1 > x2-b)
x = b - dx;
else
x = b + dx;
double fx = f.eval(x);
if (!std::isfinite(fx))
return x1;
if (b-x1 > x2-b)
{
// x is on the left from b
if (f1 > fx && fx < fb)
{
// pickup bracket [f1,fx,fb]
x2 = b;
fb = fx;
b = x;
}
else
{
// pickup bracket [fx,fb,fx2]
f1 = fx;
x1 = x;
}
}
else
{
// x is on the right from b
if (f1 > fb && fb < fx)
{
// pickup bracket [f1,fb,fx]
x2 = x;
}
else
{
// pickup bracket [fb,fx,f2]
f1 = fb;
x1 = b;
fb = fx;
b = x;
}
}
}
while (dx > tol);
}
return b;
}
bool
GoldenSectionSearch::init_bracket(OneDFunction &f, double x1, double &x2, double &b)
{
double f1 = f.eval(x1);
if (!std::isfinite(f1))
throw DynareException(__FILE__, __LINE__,
"Safer point not finite in GoldenSectionSearch::init_bracket");
int cnt = 0;
bool bracket_found = false;
do
{
bool finite_found = search_for_finite(f, x1, x2, b);
if (!finite_found)
{
b = x1;
return false;
}
double f2 = f.eval(x2);
double fb = f.eval(b);
double bsym = 2*x2 - b;
double fbsym = f.eval(bsym);
// now we know that f1, f2, and fb are finite
if (std::isfinite(fbsym))
{
// we have four numbers f1, fb, f2, fbsym, we test for the
// following combinations to find the bracket:
// [f1,f2,fbsym], [f1,fb,fbsym] and [f1,fb,fbsym]
if (f1 > f2 && f2 < fbsym)
{
bracket_found = true;
b = x2;
x2 = bsym;
}
else if (f1 > fb && fb < fbsym)
{
bracket_found = true;
x2 = bsym;
}
else if (f1 > fb && fb < f2)
{
bracket_found = true;
}
else
{
double newx2 = b;
// choose the smallest value in case we end
if (f1 > fbsym)
{
// the smallest value is on the other end, we do
// not want to continue
b = bsym;
return false;
}
else
b = x1;
// move x2 to b in case we continue
x2 = newx2;
}
}
else
{
// we have only three numbers, we test for the bracket,
// and if not found, we set b as potential result and
// shorten x2 as potential init value for next cycle
if (f1 > fb && fb < f2)
bracket_found = true;
else
{
double newx2 = b;
// choose the smaller value in case we end
if (f1 > f2)
b = x2;
else
b = x1;
// move x2 to b in case we continue
x2 = newx2;
}
}
cnt++;
}
while (!bracket_found && cnt < 5);
return bracket_found;
}
/** This moves x2 toward to x1 until the function at x2 is finite and
* b as a golden section between x1 and x2 yields also finite f. */
bool
GoldenSectionSearch::search_for_finite(OneDFunction &f, double x1, double &x2, double &b)
{
int cnt = 0;
bool found = false;
do
{
double f2 = f.eval(x2);
b = (1-golden)*x1 + golden*x2;
double fb = f.eval(b);
found = std::isfinite(f2) && std::isfinite(fb);
if (!found)
x2 = b;
cnt++;
}
while (!found && cnt < 5);
return found;
}
void
VectorFunction::check_for_eval(const ConstVector &in, Vector &out) const
{
if (inDim() != in.length() || outDim() != out.length())
throw DynareException(__FILE__, __LINE__,
"Wrong dimensions in VectorFunction::check_for_eval");
}
double
NLSolver::eval(double lambda)
{
Vector xx((const Vector &)x);
xx.add(1-lambda, xcauchy);
xx.add(lambda, xnewton);
Vector ff(func.outDim());
func.eval(xx, ff);
return ff.dot(ff);
}
bool
NLSolver::solve(Vector &xx, int &iter)
{
JournalRecord rec(journal);
rec << "Iter lambda residual" << endrec;
JournalRecord rec1(journal);
rec1 << "---------------------------" << endrec;
char tmpbuf[14];
x = (const Vector &) xx;
iter = 0;
// setup fx
Vector fx(func.outDim());
func.eval(x, fx);
if (!fx.isFinite())
throw DynareException(__FILE__, __LINE__,
"Initial guess does not yield finite residual in NLSolver::solve");
bool converged = fx.getMax() < tol;
JournalRecord rec2(journal);
sprintf(tmpbuf, "%10.6g", fx.getMax());
rec2 << iter << " N/A " << tmpbuf << endrec;
while (!converged && iter < max_iter)
{
// setup Jacobian
jacob.eval(x);
// calculate cauchy step
Vector g(func.inDim());
g.zeros();
ConstTwoDMatrix(jacob).multaVecTrans(g, fx);
Vector Jg(func.inDim());
Jg.zeros();
ConstTwoDMatrix(jacob).multaVec(Jg, g);
double m = -g.dot(g)/Jg.dot(Jg);
xcauchy = (const Vector &) g;
xcauchy.mult(m);
// calculate newton step
xnewton = (const Vector &) fx;
ConstTwoDMatrix(jacob).multInvLeft(xnewton);
xnewton.mult(-1);
// line search
double lambda = GoldenSectionSearch::search(*this, 0, 1);
x.add(1-lambda, xcauchy);
x.add(lambda, xnewton);
// evaluate func
func.eval(x, fx);
converged = fx.getMax() < tol;
// iter
iter++;
JournalRecord rec3(journal);
sprintf(tmpbuf, "%10.6g", fx.getMax());
rec3 << iter << " " << lambda << " " << tmpbuf << endrec;
}
xx = (const Vector &) x;
return converged;
}
dynare++/src/nlsolve.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2006, Ondra Kamenik
// $Id: nlsolve.h 762 2006-05-22 13:00:07Z kamenik $
#ifndef OGU_NLSOLVE_H
#define OGU_NLSOLVE_H
#include "twod_matrix.hh"
#include "journal.hh"
namespace ogu
{
class OneDFunction
{
public:
virtual ~OneDFunction()
= default;
virtual double eval(double) = 0;
};
class GoldenSectionSearch
{
protected:
static double tol;
static double golden;
public:
static double search(OneDFunction &f, double x1, double x2);
protected:
/** This initializes a bracket by moving x2 and b (as a golden
* section of x1,x2) so that f(x1)>f(b) && f(b)<f(x2). The point
* x1 is not moved, since it is considered as reliable and f(x1)
* is supposed to be finite. If initialization of a bracket
* succeeded, then [x1, b, x2] is the bracket and true is
* returned. Otherwise, b is the minimum found and false is
* returned. */
static bool init_bracket(OneDFunction &f, double x1, double &x2, double &b);
/** This supposes that f(x1) is finite and it moves x2 toward x1
* until x2 and b (as a golden section of x1,x2) are finite. If
* succeeded, the routine returns true and x2, and b. Otherwise,
* it returns false. */
static bool search_for_finite(OneDFunction &f, double x1, double &x2, double &b);
};
class VectorFunction
{
public:
VectorFunction()
= default;
virtual ~VectorFunction()
= default;
virtual int inDim() const = 0;
virtual int outDim() const = 0;
/** Check dimensions of eval parameters. */
void check_for_eval(const ConstVector &in, Vector &out) const;
/** Evaluate the vector function. */
virtual void eval(const ConstVector &in, Vector &out) = 0;
};
class Jacobian : public TwoDMatrix
{
public:
Jacobian(int n)
: TwoDMatrix(n, n)
{
}
~Jacobian()
override = default;
virtual void eval(const Vector &in) = 0;
};
class NLSolver : public OneDFunction
{
protected:
Journal &journal;
VectorFunction &func;
Jacobian &jacob;
const int max_iter;
const double tol;
private:
Vector xnewton;
Vector xcauchy;
Vector x;
public:
NLSolver(VectorFunction &f, Jacobian &j, int maxit, double tl, Journal &jr)
: journal(jr), func(f), jacob(j), max_iter(maxit), tol(tl),
xnewton(f.inDim()), xcauchy(f.inDim()), x(f.inDim())
{
xnewton.zeros(); xcauchy.zeros(); x.zeros();
}
~NLSolver()
override = default;
/** Returns true if the problem has converged. xx as input is the
* starting value, as output it is a solution. */
bool solve(Vector &xx, int &iter);
/** To implement OneDFunction interface. It returns
* func(xx)^T*func(xx), where
* xx=x+lambda*xcauchy+(1-lambda)*xnewton. It is non-const only
* because it calls func, x, xnewton, xcauchy is not changed. */
double eval(double lambda) override;
};
};
#endif
dynare++/src/planner_builder.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2006, Ondra Kamenik
// $Id$
#include "planner_builder.hh"
#include "dynare_exception.hh"
#include <cmath>
#include <utility>
using namespace ogdyn;
const IntegerMatrix &
IntegerMatrix::operator=(const IntegerMatrix &im)
{
if (nr != im.nr || nc != im.nc)
throw DynareException(__FILE__, __LINE__,
"Matrices have different dimensions in IntegerMatrix::operator=");
memcpy(data, im.data, nr*nc*sizeof(int));
return *this;
}
const IntegerArray3 &
IntegerArray3::operator=(const IntegerArray3 &ia3)
{
if (n1 != ia3.n1 || n2 != ia3.n2 || n3 != ia3.n3)
throw DynareException(__FILE__, __LINE__,
"Arrays have different dimensions in IntegerArray3::operator=");
memcpy(data, ia3.data, n1*n2*n3*sizeof(int));
return *this;
}
PlannerBuilder::PlannerBuilder(DynareModel &m, const Tvarset &yyset,
Teqset ffset)
: yset(), fset(std::move(ffset)), model(m),
tb(model.t_plobjective), tbeta(model.t_pldiscount),
maxlead(model.atoms.get_maxlead()),
minlag(model.atoms.get_minlag()),
diff_b(yyset.size(), 1-minlag),
diff_f(yyset.size(), fset.size(), 1+maxlead-minlag),
static_atoms(),
static_tree(),
diff_b_static(yyset.size(), 1-minlag),
diff_f_static(yyset.size(), fset.size(), 1+maxlead-minlag)
{
info.num_new_terms -= model.getParser().getTree().get_num_op();
fill_yset(m.atoms.get_name_storage(), yyset);
add_derivatives_of_b();
add_derivatives_of_f();
shift_derivatives_of_b();
shift_derivatives_of_f();
beta_multiply_b();
beta_multiply_f();
make_static_version();
lagrange_mult_f();
form_equations();
info.num_new_terms += model.getParser().getTree().get_num_op();
}
PlannerBuilder::PlannerBuilder(const PlannerBuilder &pb, ogdyn::DynareModel &m)
: yset(), fset(pb.fset), model(m),
tb(pb.tb), tbeta(pb.tbeta),
maxlead(pb.maxlead), minlag(pb.minlag),
diff_b(pb.diff_b), diff_f(pb.diff_f),
static_atoms(pb.static_atoms),
static_tree(pb.static_tree),
diff_b_static(pb.diff_b_static),
diff_f_static(pb.diff_f_static),
aux_map(), static_aux_map()
{
fill_yset(m.atoms.get_name_storage(), pb.yset);
fill_aux_map(m.atoms.get_name_storage(), pb.aux_map, pb.static_aux_map);
}
void
PlannerBuilder::add_derivatives_of_b()
{
int yi = 0;
for (Tvarset::const_iterator yname = yset.begin();
yname != yset.end(); ++yname, yi++)
for (int ll = minlag; ll <= 0; ll++)
{
int yt = model.atoms.index(*yname, ll);
if (yt != -1)
diff_b(yi, ll-minlag) = model.eqs.add_derivative(tb, yt);
else
diff_b(yi, ll-minlag) = ogp::OperationTree::zero;
}
}
void
PlannerBuilder::add_derivatives_of_f()
{
int yi = 0;
for (Tvarset::const_iterator yname = yset.begin();
yname != yset.end(); ++yname, yi++)
for (unsigned int fi = 0; fi < fset.size(); fi++)
for (int ll = minlag; ll <= maxlead; ll++)
{
int yt = model.atoms.index(*yname, ll);
if (yt != -1)
diff_f(yi, fi, ll-minlag)
= model.eqs.add_derivative(model.eqs.formula(fset[fi]), yt);
else
diff_f(yi, fi, ll-minlag) = ogp::OperationTree::zero;
}
}
void
PlannerBuilder::shift_derivatives_of_b()
{
map<int, int> subst;
for (int yi = 0; yi < diff_b.nrows(); yi++)
for (int ll = minlag; ll < 0; ll++)
if (diff_b(yi, ll-minlag) != ogp::OperationTree::zero)
{
model.variable_shift_map(model.eqs.nulary_of_term(diff_b(yi, ll-minlag)),
-ll, subst);
diff_b(yi, ll-minlag) = model.eqs.add_substitution(diff_b(yi, ll-minlag), subst);
}
}
void
PlannerBuilder::shift_derivatives_of_f()
{
map<int, int> subst;
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
{
// first do it leads which are put under expectation before t: no problem
for (int ll = 0; ll <= maxlead; ll++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
{
model.variable_shift_map(model.eqs.nulary_of_term(diff_f(yi, fi, ll-minlag)),
-ll, subst);
diff_f(yi, fi, ll-minlag)
= model.eqs.add_substitution(diff_f(yi, fi, ll-minlag), subst);
}
// now do it for lags, these are put as leads under
// expectations after time t, so we have to introduce
// auxiliary variables at time t, and make leads of them here
for (int ll = minlag; ll < 0; ll++)
{
int ft = diff_f(yi, fi, ll-minlag);
if (ft != ogp::OperationTree::zero)
{
// if the ft term has a lead, than we need to
// introduce an auxiliary variable z_t, define it
// as E_t[ft] and put z_{t-ll} to the
// equation. Otherwise, we just put leaded ft to
// the equation directly.
int ft_maxlead, ft_minlag;
model.termspan(ft, ft_maxlead, ft_minlag);
if (ft_maxlead > 0)
{
// make an auxiliary variable
char name[100];
sprintf(name, "AUX_%d_%d_%d", yi, fset[fi], -ll);
model.atoms.register_uniq_endo(name);
info.num_aux_variables++;
int taux = model.eqs.add_nulary(name);
sprintf(name, "AUX_%d_%d_%d(%d)", yi, fset[fi], -ll, -ll);
int taux_leaded = model.eqs.add_nulary(name);
// put aux_leaded to the equation
diff_f(yi, fi, ll-minlag) = taux_leaded;
// save auxiliary variable and the term
aux_map.insert(Tsubstmap::value_type(model.atoms.name(taux), ft));
}
else
{
// no auxiliary variable is needed and the
// term ft can be leaded in place
model.variable_shift_map(model.eqs.nulary_of_term(ft), -ll, subst);
diff_f(yi, fi, ll-minlag)
= model.eqs.add_substitution(ft, subst);
}
}
}
}
}
void
PlannerBuilder::beta_multiply_b()
{
int beta_pow = ogp::OperationTree::one;
for (int ll = 0; ll >= minlag; ll--,
beta_pow = model.eqs.add_binary(ogp::TIMES, beta_pow, tbeta))
for (int yi = 0; yi < diff_b.nrows(); yi++)
if (diff_b(yi, ll-minlag) != ogp::OperationTree::zero)
diff_b(yi, ll-minlag)
= model.eqs.add_binary(ogp::TIMES, beta_pow, diff_b(yi, ll-minlag));
}
void
PlannerBuilder::beta_multiply_f()
{
int beta_pow = ogp::OperationTree::one;
for (int ll = 0; ll <= maxlead; ll++,
beta_pow = model.eqs.add_binary(ogp::DIVIDE, beta_pow, tbeta))
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
diff_f(yi, fi, ll-minlag)
= model.eqs.add_binary(ogp::TIMES, beta_pow, diff_f(yi, fi, ll-minlag));
beta_pow = ogp::OperationTree::one;
for (int ll = 0; ll >= minlag; ll--,
beta_pow = model.eqs.add_binary(ogp::TIMES, beta_pow, tbeta))
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
diff_f(yi, fi, ll-minlag)
= model.eqs.add_binary(ogp::TIMES, beta_pow, diff_f(yi, fi, ll-minlag));
}
void
PlannerBuilder::make_static_version()
{
// map holding substitutions from dynamic to static
ogp::StaticFineAtoms::Tintintmap tmap;
// fill static atoms with outer ordering
static_atoms.import_atoms(model.atoms, static_tree, tmap);
// go through diff_b and fill diff_b_static
for (int ll = minlag; ll <= 0; ll++)
for (int yi = 0; yi < diff_b.nrows(); yi++)
diff_b_static(yi, ll-minlag)
= static_tree.add_substitution(diff_b(yi, ll-minlag),
tmap, model.eqs.getTree());
// go through diff_f and fill diff_f_static
for (int ll = minlag; ll <= maxlead; ll++)
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
diff_f_static(yi, fi, ll-minlag)
= static_tree.add_substitution(diff_f(yi, fi, ll-minlag),
tmap, model.eqs.getTree());
// go through aux_map and fill static_aux_map
for (Tsubstmap::const_iterator it = aux_map.begin();
it != aux_map.end(); ++it)
{
int tstatic = static_tree.add_substitution((*it).second, tmap, model.eqs.getTree());
const char *name = static_atoms.get_name_storage().query((*it).first);
static_aux_map.insert(Tsubstmap::value_type(name, tstatic));
}
}
void
PlannerBuilder::lagrange_mult_f()
{
// register multipliers
char mult_name[100];
for (int fi = 0; fi < diff_f.dim2(); fi++)
{
sprintf(mult_name, "MULT%d", fset[fi]);
model.atoms.register_uniq_endo(mult_name);
info.num_lagrange_mults++;
}
// multiply with the multipliers
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
for (int ll = minlag; ll <= maxlead; ll++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
{
sprintf(mult_name, "MULT%d(%d)", fset[fi], -ll);
int tm = model.eqs.add_nulary(mult_name);
diff_f(yi, fi, ll-minlag)
= model.eqs.add_binary(ogp::TIMES, tm, diff_f(yi, fi, ll-minlag));
}
}
void
PlannerBuilder::form_equations()
{
// add planner's FOCs
for (int yi = 0; yi < diff_f.dim1(); yi++)
{
int eq = ogp::OperationTree::zero;
for (int ll = minlag; ll <= 0; ll++)
eq = model.eqs.add_binary(ogp::PLUS, eq, diff_b(yi, ll-minlag));
for (int fi = 0; fi < diff_f.dim2(); fi++)
for (int ll = minlag; ll <= maxlead; ll++)
eq = model.eqs.add_binary(ogp::PLUS, eq, diff_f(yi, fi, ll-minlag));
model.eqs.add_formula(eq);
}
// add equations for auxiliary variables
for (Tsubstmap::const_iterator it = aux_map.begin();
it != aux_map.end(); ++it)
{
int t = model.atoms.index((*it).first, 0);
model.eqs.add_formula(model.eqs.add_binary(ogp::MINUS, t, (*it).second));
}
}
void
PlannerBuilder::fill_yset(const ogp::NameStorage &ns,
const PlannerBuilder::Tvarset &yyset)
{
for (auto it : yyset)
yset.insert(ns.query(it));
}
void
PlannerBuilder::fill_aux_map(const ogp::NameStorage &ns, const Tsubstmap &aaux_map,
const Tsubstmap &astatic_aux_map)
{
// fill aux_map
for (auto it : aaux_map)
aux_map.insert(Tsubstmap::value_type(ns.query(it.first), it.second));
// fill static_aux_map
for (auto it : astatic_aux_map)
static_aux_map.insert(Tsubstmap::value_type(static_atoms.get_name_storage().query(it.first),
it.second));
}
MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy)
: builder(pb), b(builder.diff_b_static.nrows()),
F(builder.diff_f_static.dim1(), builder.diff_f_static.dim2())
{
b.zeros();
F.zeros();
// first evaluate substitutions (auxiliary variables) from the builder
ogdyn::DynareStaticSteadySubstitutions dss(builder.model.atoms, builder.static_atoms,
builder.static_tree,
builder.static_aux_map, pvals, yy);
// gather all the terms from builder.diff_b_static and
// builder.diff_f_static to the vector, the ordering is important,
// since the index of this vector will have to be decoded to the
// position in b and F.
vector<int> terms;
for (int yi = 0; yi < builder.diff_b_static.nrows(); yi++)
for (int l = 0; l < builder.diff_b_static.ncols(); l++)
terms.push_back(builder.diff_b_static(yi, l));
for (int yi = 0; yi < builder.diff_f_static.dim1(); yi++)
for (int fi = 0; fi < builder.diff_f_static.dim2(); fi++)
for (int l = 0; l < builder.diff_f_static.dim3(); l++)
terms.push_back(builder.diff_f_static(yi, fi, l));
// evaluate the terms, it will call a series of load(i,res), which
// sum the results through lags/leads to b and F
DynareStaticSteadyAtomValues dssav(builder.model.atoms, builder.static_atoms, pvals, yy);
ogp::FormulaCustomEvaluator fe(builder.static_tree, terms);
fe.eval(dssav, *this);
// solve overdetermined system b+F*lambda=0 using SVD decomposition
SVDDecomp decomp(F);
Vector lambda(builder.diff_f_static.dim2());
decomp.solve(b, lambda);
lambda.mult(-1);
// take values of lambda and put it to yy
for (int fi = 0; fi < builder.diff_f_static.dim2(); fi++)
{
char mult_name[100];
sprintf(mult_name, "MULT%d", builder.fset[fi]);
int iouter = builder.model.atoms.name2outer_endo(mult_name);
int iy = builder.model.atoms.outer2y_endo()[iouter];
if (!std::isfinite(yy[iy]))
yy[iy] = lambda[fi];
// go through all substitutions of the multiplier and set them
// as well
if (builder.model.atom_substs)
{
const ogp::AtomSubstitutions::Toldnamemap &old2new
= builder.model.atom_substs->get_old2new();
const ogp::AtomSubstitutions::Toldnamemap::const_iterator it
= old2new.find(mult_name);
if (it != old2new.end())
{
const ogp::AtomSubstitutions::Tshiftnameset &sset = (*it).second;
for (const auto & itt : sset)
{
const char *newname = itt.first;
int iouter = builder.model.atoms.name2outer_endo(newname);
int iy = builder.model.atoms.outer2y_endo()[iouter];
if (!std::isfinite(yy[iy]))
yy[iy] = lambda[fi];
}
}
}
}
}
void
MultInitSS::load(int i, double res)
{
// we can afford it, since the evaluator sets res to exact zero if
// the term is zero
if (res == 0)
return;
// decode i and add to either b or F
if (i < builder.diff_b_static.nrows()*builder.diff_b_static.ncols())
{
// add to b
b[i / builder.diff_b_static.ncols()] += res;
}
else
{
// add to F
i -= builder.diff_b_static.nrows()*builder.diff_b_static.ncols();
int yifi = i / builder.diff_f_static.dim3();
int yi = yifi / builder.diff_f_static.dim2();
int fi = yifi % builder.diff_f_static.dim2();
F.get(yi, fi) += res;
}
}
dynare++/src/planner_builder.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2006-2011, Ondra Kamenik
#ifndef PLANNER_BUILDER_H
#define PLANNER_BUILDER_H
#include "dynare_model.hh"
namespace ogdyn
{
using std::unordered_set;
using std::map;
using std::vector;
/** This is a two dimensional array of integers. Nothing
* difficult. */
class IntegerMatrix
{
protected:
/** Number of rows. */
int nr;
/** Number of columns. */
int nc;
/** The pointer to the data. */
int *data;
public:
/** Construct uninitialized array. */
IntegerMatrix(int nrr, int ncc)
: nr(nrr), nc(ncc), data(new int[nr*nc])
{
}
/** Copy constructor. */
IntegerMatrix(const IntegerMatrix &im)
: nr(im.nr), nc(im.nc), data(new int[nr*nc])
{
memcpy(data, im.data, nr*nc*sizeof(int));
}
virtual ~IntegerMatrix()
{
delete [] data;
}
/** Assignment operator. It can only assing array with the
* same dimensions. */
const IntegerMatrix &operator=(const IntegerMatrix &im);
int &
operator()(int i, int j)
{
return data[i+j*nr];
}
const int &
operator()(int i, int j) const
{
return data[i+j*nr];
}
int
nrows() const
{
return nr;
}
int
ncols() const
{
return nc;
}
};
/** The three dimensional array of integers. Nothing difficult. */
class IntegerArray3
{
protected:
/** First dimension. */
int n1;
/** Second dimension. */
int n2;
/** Third dimension. */
int n3;
/** The data. */
int *data;
public:
/** Constrcut unitialized array. */
IntegerArray3(int nn1, int nn2, int nn3)
: n1(nn1), n2(nn2), n3(nn3), data(new int[n1*n2*n3])
{
}
/** Copy constructor. */
IntegerArray3(const IntegerArray3 &ia3)
: n1(ia3.n1), n2(ia3.n2), n3(ia3.n3), data(new int[n1*n2*n3])
{
memcpy(data, ia3.data, n1*n2*n3*sizeof(int));
}
virtual ~IntegerArray3()
{
delete [] data;
}
/** Assignment operator assigning the arrays with the same dimensions. */
const IntegerArray3 &operator=(const IntegerArray3 &ia3);
int &
operator()(int i, int j, int k)
{
return data[i+j*n1+k*n1*n2];
}
const int &
operator()(int i, int j, int k) const
{
return data[i+j*n1+k*n1*n2];
}
int
dim1() const
{
return n1;
}
int
dim2() const
{
return n2;
}
int
dim3() const
{
return n3;
}
};
/** This struct encapsulates information about the building of a
* planner's problem. */
struct PlannerInfo
{
int num_lagrange_mults{0};
int num_aux_variables{0};
int num_new_terms{0};
PlannerInfo()
= default;
};
class MultInitSS;
/** This class builds the first order conditions of the social
* planner problem with constraints being the equations in the
* model. The model is non-const parameter to the constructor
* which adds appropriate FOCs to the system. It also allows for
* an estimation of the lagrange multipliers given all other
* endogenous variables of the static system. For this purpose we
* need to create static atoms and static versions of all the tree
* index matrices. The algorithm and algebra are documented in
* dynare++-ramsey.pdf. */
class PlannerBuilder
{
friend class MultInitSS;
public:
/** Type for a set of variable names. */
using Tvarset = unordered_set<const char *>;
/** Type for a set of equations. An equation is identified by
* an index to an equation in the equation vector given by
* DynareModel::eqs. The tree index of the i-th formula is
* retrieved as DynareModel::egs.formula(i). */
using Teqset = vector<int>;
protected:
/** This is a set of variables wrt which the planner
* optimizes. These could be all endogenous variables, but it
* is beneficial to exclude all variables which are
* deterministic transformations of past exogenous variables,
* since the planner cannot influence them. This could save a
* few equations. This is not changed after it is constructed,
* but it is constructed manually, so it cannot be declared as
* const. */
Tvarset yset;
/** These are the equation indices constituing the constraints
* for the planner. Again, it is beneficial to exclude all
* equations defining exogenous variables excluded from
* yset. */
const Teqset fset;
/** Reference to the model. */
ogdyn::DynareModel &model;
/** Tree index of the planner objective. */
int tb;
/** Tree index of the planner discount parameter. */
int tbeta;
/** The maximum lead in the model including the planner's
* objective before building the planner's FOCs. */
const int maxlead;
/** The minimum lag in the model including the planner's objective
* before building the planner's FOCs. */
const int minlag;
/** Tree indices of formulas in the planner FOCs involving
* derivatives of the planner's objective. Rows correspond to the
* endogenous variables, columns correspond to lags in the
* objective function. The contents of the matrix will evolve as
* the algorithm proceeds. */
IntegerMatrix diff_b;
/** Tree indices of formulas in the planner FOCs involving
* derivatives of the model equations (constraints). The first
* dimension corresponds to endogenous variables, the second to
* the constraints, the third to lags or leads of endogenous
* variables in the constraints. The contents of the array will
* evolve as the algorithm proceeds.*/
IntegerArray3 diff_f;
/** Static version of the model atoms. It is needed to build
* static version of diff_b and diff_f. */
ogp::StaticFineAtoms static_atoms;
/** Static version of all the trees of diff_b and diff_f build
* over static_atoms. */
ogp::OperationTree static_tree;
/** Tree indices of static version of diff_b over static_atoms and static_tree. */
IntegerMatrix diff_b_static;
/** Tree indices of static version of diff_f over static_atoms
* and static_tree. This member is created before calling
* lagrange_mult_f(), so it does not contain the
* multiplication with the lagrange multipliers. */
IntegerArray3 diff_f_static;
/** Auxiliary variables mapping. During the algorithm, some
* auxiliary variables for the terms might be created, so we
* remember their names and tree indices of the terms. This
* maps a name to the tree index of an expression equal to the
* auxiliary variable at time zero. The auxiliary variables
* names point to the dynamic atoms storage, tree inidices to
* the dynamic model tree. */
Tsubstmap aux_map;
/** Static version of aux_map. The names point to static_atoms
* storage, the tree indices to the static_tree. */
Tsubstmap static_aux_map;
/** Information about the number of various things. */
PlannerInfo info;
public:
/** Build the planner problem for the given model optimizing
* through the given endogenous variables with the given
* constraints. We allow for a selection of a subset of
* equations and variables in order to eliminate exogenous
* predetermined process which cannot be influenced by the
* social planner. */
PlannerBuilder(ogdyn::DynareModel &m, const Tvarset &yyset,
Teqset ffset);
/** Construct a copy of the builder with provided model, which
* is supposed to be the copy of the model in the builder. */
PlannerBuilder(const PlannerBuilder &pb, ogdyn::DynareModel &m);
/** Avoid copying from only PlannerBuilder. */
PlannerBuilder(const PlannerBuilder &pb) = delete;
/** Return the information. */
const PlannerInfo &
get_info() const
{
return info;
}
protected:
/** Differentiate the planner objective wrt endogenous
* variables with different lags. */
void add_derivatives_of_b();
/** Differentiate the constraints wrt endogenous variables
* with different lags and leads. */
void add_derivatives_of_f();
/** Shift derivatives of diff_b. */
void shift_derivatives_of_b();
/** Shift derivatives of diff_ff. */
void shift_derivatives_of_f();
/** Multiply with the discount factor terms in diff_b. */
void beta_multiply_b();
/** Multiply with the discount factor terms in diff_f. */
void beta_multiply_f();
/** Fill static_atoms and static_tree and build diff_b_static,
* diff_f_static and aux_map_static with static versions of diff_b,
* diff_f and aux_map. */
void make_static_version();
/** Multiply diff_f with Langrange multipliers. */
void lagrange_mult_f();
/** Add the equations to the mode, including equation for auxiliary variables. */
void form_equations();
private:
/** Fill yset for a given yyset and given name storage. */
void fill_yset(const ogp::NameStorage &ns, const Tvarset &yyset);
/** Fill aux_map and aux_map_static for a given aaux_map and
* aaux_map_static for a given storage of dynamic atoms (used
* for aux_map) and static atoms storage from this object for
* aux_map_static. */
void fill_aux_map(const ogp::NameStorage &ns, const Tsubstmap &aaux_map,
const Tsubstmap &astatic_aux_map);
};
/** This class only calculates for the given initial guess of
* endogenous variables, initial guess of the Langrange
* multipliers of the social planner problem yielding the least
* square error. It is used by just calling its constructor. The
* constructor takes non-const reference to the vector of
* endogenous variables, calculates lambdas and put the values of
* lambdas to the vector. The algbera is found in
* dynare++-ramsey.pdf.
*
* The code can be run only after the parsing has been finished in
* atoms. */
class MultInitSS : public ogp::FormulaEvalLoader
{
protected:
/** The constant reference to the builder. */
const PlannerBuilder &builder;
/** The constant term of the problem. Its length is the number
* of endogenous variable wrt the planner optimizes. */
Vector b;
/** The matrix of the overdetermined problem. The number of
* rows is equal to the number of endogenous variables wrt
* which the planner optimizes, the number of columns is equal
* to the number of Langrange multipliers which is equal to
* the number of constraints which is smaller than the number
* of endogenous variables. Hence the system b+F*lambda=0 is
* overdetermined. */
GeneralMatrix F;
public:
/** The constructor of the object which does everything. Its
* main goal is to update yy. Note that if an item of yy
* corresponding to a lagrange multiplier is already set, it
* is not reset. */
MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy);
/** This loads evaluated parts of b or F and decodes i and
* advances b or F depending on the decoded i. The decoding is
* dependent on the way how the terms of builder.diff_b and
* builder.diff_f_save have been put the the
* ogp::FormulaCustomEvaluator. This is documented in the code
* of the constructor. */
void load(int i, double res) override;
};
};
#endif
dynare++/sylv/Makefile.am deleted 100644 → 0
View file @ 219d2bb3
SUBDIRS = cc testing
EXTRA_DIST = sylvester.tex change_log.html matlab
if HAVE_PDFLATEX
pdf-local: sylvester.pdf
sylvester.pdf: sylvester.tex
$(PDFLATEX) sylvester
$(PDFLATEX) sylvester
endif
CLEANFILES = *.pdf *.log *.aux
dynare++/sylv/cc/BlockDiagonal.hh deleted 100644 → 0
View file @ 219d2bb3
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/BlockDiagonal.h,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef BLOCK_DIAGONAL_H
#define BLOCK_DIAGONAL_H
#include <memory>
#include <vector>
#include "QuasiTriangular.hh"
class BlockDiagonal : public QuasiTriangular
{
std::vector<int> row_len, col_len;
public:
BlockDiagonal(ConstVector d, int d_size);
BlockDiagonal(int p, const BlockDiagonal &b);
BlockDiagonal(const BlockDiagonal &b) = default;
BlockDiagonal(const QuasiTriangular &t);
BlockDiagonal &operator=(const QuasiTriangular &t)
{
GeneralMatrix::operator=(t);
return *this;
}
BlockDiagonal &operator=(const BlockDiagonal &b) = default;
~BlockDiagonal() override = default;
void setZeroBlockEdge(diag_iter edge);
int getNumZeros() const;
int getNumBlocks() const;
int getLargestBlock() const;
void printInfo() const;
void multKron(KronVector &x) const override;
void multKronTrans(KronVector &x) const override;
const_col_iter col_begin(const DiagonalBlock &b) const override;
col_iter col_begin(const DiagonalBlock &b) override;
const_row_iter row_end(const DiagonalBlock &b) const override;
row_iter row_end(const DiagonalBlock &b) override;
std::unique_ptr<QuasiTriangular>
clone() const override
{
return std::make_unique<BlockDiagonal>(*this);
}
private:
void setZerosToRU(diag_iter edge);
const_diag_iter findBlockStart(const_diag_iter from) const;
static void savePartOfX(int si, int ei, const KronVector &x, Vector &work);
void multKronBlock(const_diag_iter start, const_diag_iter end,
KronVector &x, Vector &work) const;
void multKronBlockTrans(const_diag_iter start, const_diag_iter end,
KronVector &x, Vector &work) const;
};
#endif /* BLOCK_DIAGONAL_H */
dynare++/sylv/cc/GeneralSylvester.hh deleted 100644 → 0
View file @ 219d2bb3
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralSylvester.h,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef GENERAL_SYLVESTER_H
#define GENERAL_SYLVESTER_H
#include "SylvMatrix.hh"
#include "SimilarityDecomp.hh"
#include "SylvesterSolver.hh"
#include <memory>
class GeneralSylvester
{
SylvParams pars;
int order;
const SqSylvMatrix a;
const SylvMatrix b;
const SqSylvMatrix c;
SylvMatrix d;
bool solved;
std::unique_ptr<SchurDecompZero> bdecomp;
std::unique_ptr<SimilarityDecomp> cdecomp;
std::unique_ptr<SylvesterSolver> sylv;
public:
/* construct with my copy of d*/
GeneralSylvester(int ord, int n, int m, int zero_cols,
const ConstVector &da, const ConstVector &db,
const ConstVector &dc, const ConstVector &dd,
const SylvParams &ps);
GeneralSylvester(int ord, int n, int m, int zero_cols,
const ConstVector &da, const ConstVector &db,
const ConstVector &dc, const ConstVector &dd,
bool alloc_for_check = false);
/* construct with provided storage for d */
GeneralSylvester(int ord, int n, int m, int zero_cols,
const ConstVector &da, const ConstVector &db,
const ConstVector &dc, Vector &dd,
bool alloc_for_check = false);
GeneralSylvester(int ord, int n, int m, int zero_cols,
const ConstVector &da, const ConstVector &db,
const ConstVector &dc, Vector &dd,
const SylvParams &ps);
virtual ~GeneralSylvester() = default;
int
getM() const
{
return c.numRows();
}
int
getN() const
{
return a.numRows();
}
const double *
getResult() const
{
return d.base();
}
const SylvParams &
getParams() const
{
return pars;
}
SylvParams &
getParams()
{
return pars;
}
void solve();
void check(const ConstVector &ds);
private:
void init();
};
#endif /* GENERAL_SYLVESTER_H */
dynare++/sylv/cc/IterativeSylvester.hh deleted 100644 → 0
View file @ 219d2bb3
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/IterativeSylvester.h,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef ITERATIVE_SYLVESTER_H
#define ITERATIVE_SYLVESTER_H
#include "SylvesterSolver.hh"
#include "KronVector.hh"
#include "QuasiTriangular.hh"
#include "SimilarityDecomp.hh"
class IterativeSylvester : public SylvesterSolver
{
public:
IterativeSylvester(const QuasiTriangular &k, const QuasiTriangular &f)
: SylvesterSolver(k, f)
{
}
IterativeSylvester(const SchurDecompZero &kdecomp, const SchurDecomp &fdecomp)
: SylvesterSolver(kdecomp, fdecomp)
{
}
IterativeSylvester(const SchurDecompZero &kdecomp, const SimilarityDecomp &fdecomp)
: SylvesterSolver(kdecomp, fdecomp)
{
}
void solve(SylvParams &pars, KronVector &x) const override;
private:
double performFirstStep(KronVector &x) const;
static double performStep(const QuasiTriangular &k, const QuasiTriangular &f,
KronVector &x);
};
#endif /* ITERATIVE_SYLVESTER_H */
dynare++/sylv/cc/KronUtils.hh deleted 100644 → 0
View file @ 219d2bb3
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/KronUtils.h,v 1.1.1.1 2004/06/04 13:00:31 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef KRON_UTILS_H
#define KRON_UTILS_H
#include "KronVector.hh"
#include "QuasiTriangular.hh"
class KronUtils
{
public:
/* multiplies I_m\otimes..\I_m\otimes T\otimes I_m...I_m\otimes I_n
with given b and returns x. T must be (m,m), number of
\otimes is b.getDepth(), level is a number of I_m's between T
and I_n plus 1. If level=0, then we multiply
\I_m\otimes ..\otimes I_m\otimes T, T is (n,n) */
static void multAtLevel(int level, const QuasiTriangular &t,
KronVector &x);
static void multAtLevelTrans(int level, const QuasiTriangular &t,
KronVector &x);
/* multiplies x=(F'\otimes F'\otimes..\otimes K)x */
static void multKron(const QuasiTriangular &f, const QuasiTriangular &k,
KronVector &x);
};
#endif /* KRON_UTILS_H */
dynare++/sylv/cc/KronVector.cc deleted 100644 → 0
View file @ 219d2bb3
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/KronVector.cpp,v 1.1.1.1 2004/06/04 13:00:31 kamenik Exp $ */
/* Tag $Name: $ */
#include "KronVector.hh"
#include "SylvException.hh"
#include <utility>
int
power(int m, int depth)
{
int p = 1;
for (int i = 0; i < depth; i++)
p *= m;
return p;
}
KronVector::KronVector(int mm, int nn, int dp)
: Vector(power(mm, dp)*nn), m(mm), n(nn), depth(dp)
{
}
KronVector::KronVector(Vector &v, int mm, int nn, int dp)
: Vector(v), m(mm), n(nn), depth(dp)
{
if (length() != power(m, depth)*n)
throw SYLV_MES_EXCEPTION("Bad conversion KronVector from Vector.");
}
KronVector::KronVector(KronVector &v, int i)
: Vector(v, i*power(v.m, v.depth-1)*v.n, power(v.m, v.depth-1)*v.n), m(v.m), n(v.n),
depth(v.depth-1)
{
if (depth < 0)
throw SYLV_MES_EXCEPTION("Bad KronVector pick, depth < 0.");
}
KronVector::KronVector(const ConstKronVector &v)
: Vector(v), m(v.m), n(v.n), depth(v.depth)
{
}
KronVector &
KronVector::operator=(const ConstKronVector &v)
{
Vector::operator=(v);
m = v.m;
n = v.n;
depth = v.depth;
return *this;
}
KronVector &
KronVector::operator=(const Vector &v)
{
if (length() != v.length())
throw SYLV_MES_EXCEPTION("Wrong lengths for vector operator =.");
Vector::operator=(v);
return *this;
}
ConstKronVector::ConstKronVector(const KronVector &v)
: ConstVector(v), m(v.m), n(v.n), depth(v.depth)
{
}
ConstKronVector::ConstKronVector(const Vector &v, int mm, int nn, int dp)
: ConstVector(v), m(mm), n(nn), depth(dp)
{
if (length() != power(m, depth)*n)
throw SYLV_MES_EXCEPTION("Bad conversion KronVector from Vector.");
}
ConstKronVector::ConstKronVector(ConstVector v, int mm, int nn, int dp)
: ConstVector(std::move(v)), m(mm), n(nn), depth(dp)
{
if (length() != power(m, depth)*n)
throw SYLV_MES_EXCEPTION("Bad conversion KronVector from Vector.");
}
ConstKronVector::ConstKronVector(const KronVector &v, int i)
: ConstVector(v, i*power(v.getM(), v.getDepth()-1)*v.getN(),
power(v.getM(), v.getDepth()-1)*v.getN()),
m(v.getM()), n(v.getN()), depth(v.getDepth()-1)
{
if (depth < 0)
throw SYLV_MES_EXCEPTION("Bad KronVector pick, depth < 0.");
}
ConstKronVector::ConstKronVector(const ConstKronVector &v, int i)
: ConstVector(v, i*power(v.m, v.depth-1)*v.n, power(v.m, v.depth-1)*v.n),
m(v.getM()), n(v.getN()), depth(v.getDepth()-1)
{
if (depth < 0)
throw SYLV_MES_EXCEPTION("Bad KronVector pick, depth < 0.");
}
dynare++/sylv/cc/KronVector.hh deleted 100644 → 0
View file @ 219d2bb3
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/KronVector.h,v 1.1.1.1 2004/06/04 13:00:31 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef KRON_VECTOR_H
#define KRON_VECTOR_H
#include "Vector.hh"
class ConstKronVector;
class KronVector : public Vector
{
friend class ConstKronVector;
protected:
int m{0};
int n{0};
int depth{0};
public:
KronVector() = default;
KronVector(const KronVector &v) = default;
KronVector(KronVector &&v) = default;
KronVector(int mm, int nn, int dp); // new instance
KronVector(Vector &v, int mm, int nn, int dp); // conversion
KronVector(KronVector &, int i); // picks i-th subvector
// We don't want implict conversion from ConstKronVector, since it's expensive
explicit KronVector(const ConstKronVector &v); // new instance and copy
KronVector &operator=(const KronVector &v) = default;
KronVector &operator=(KronVector &&v) = default;
KronVector &operator=(const ConstKronVector &v);
KronVector &operator=(const Vector &v);
int
getM() const
{
return m;
}
int
getN() const
{
return n;
}
int
getDepth() const
{
return depth;
}
};
class ConstKronVector : public ConstVector
{
friend class KronVector;
protected:
int m;
int n;
int depth;
public:
// Implicit conversion from KronVector is ok, since it's cheap
ConstKronVector(const KronVector &v);
ConstKronVector(const ConstKronVector &v) = default;
ConstKronVector(ConstKronVector &&v) = default;
ConstKronVector(const Vector &v, int mm, int nn, int dp);
ConstKronVector(ConstVector v, int mm, int nn, int dp);
ConstKronVector(const KronVector &v, int i);
ConstKronVector(const ConstKronVector &v, int i);
ConstKronVector &operator=(const ConstKronVector &v) = delete;
ConstKronVector &operator=(ConstKronVector &&v) = delete;
int
getM() const
{
return m;
}
int
getN() const
{
return n;
}
int
getDepth() const
{
return depth;
}
};
int power(int m, int depth);
#endif /* KRON_VECTOR */