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 5074 deletions
dynare++/parser/cc/namelist.lex deleted 100644 → 0
View file @ 6ca649c4
%{
#include "location.h"
#include "namelist_tab.hh"
extern YYLTYPE namelist_lloc;
#define YY_USER_ACTION SET_LLOC(namelist_);
%}
%option nounput
%option noyy_top_state
%option stack
%option prefix="namelist_"
%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]
/* names */
[A-Za-z_][A-Za-z0-9_]* {
namelist_lval.string = namelist_text;
return NAME;
}
, {return COMMA;}
. {
namelist_lval.character = namelist_text[0];
return CHARACTER;
}
%%
int namelist_wrap()
{
return 1;
}
void namelist__destroy_buffer(void* p)
{
namelist__delete_buffer((YY_BUFFER_STATE)p);
}
dynare++/parser/cc/namelist.y deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2007-2011, Ondra Kamenik
%{
#include "location.h"
#include "namelist.h"
#include "namelist_tab.hh"
void namelist_error(const char*);
int namelist_lex(void);
extern ogp::NameListParser* name_list_parser;
%}
%union {
int integer;
char *string;
char character;
}
%token COMMA CHARACTER
%token <string> NAME;
%name-prefix="namelist_"
%locations
%error-verbose
%%
namelist : namelist NAME {name_list_parser->add_name($2);}
| namelist COMMA NAME {name_list_parser->add_name($3);}
| NAME {name_list_parser->add_name($1);}
;
%%
void namelist_error(const char* mes)
{
name_list_parser->namelist_error(mes);
}
dynare++/parser/cc/parser_exception.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: parser_exception.cpp 2269 2008-11-23 14:33:22Z michel $
#include "parser_exception.h"
#include <cstring>
#include <cstdio>
using namespace ogp;
ParserException::ParserException(const char* m, int offset)
: mes(new char[strlen(m)+1]), off(offset),
aux_i1(-1), aux_i2(-1), aux_i3(-1)
{
strcpy(mes, m);
}
ParserException::ParserException(const string& m, int offset)
: mes(new char[m.size()+1]), off(offset),
aux_i1(-1), aux_i2(-1), aux_i3(-1)
{
strncpy(mes, m.c_str(), m.size());
mes[m.size()] = '\0';
}
ParserException::ParserException(const string& m, const char* dum, int i1)
: mes(new char[m.size()+1]), off(0),
aux_i1(i1), aux_i2(-1), aux_i3(-1)
{
strncpy(mes, m.c_str(), m.size());
mes[m.size()] = '\0';
}
ParserException::ParserException(const string& m, const char* dum, int i1, int i2)
: mes(new char[m.size()+1]), off(0),
aux_i1(i1), aux_i2(i2), aux_i3(-1)
{
strncpy(mes, m.c_str(), m.size());
mes[m.size()] = '\0';
}
ParserException::ParserException(const string& m, const char* dum, int i1, int i2, int i3)
: mes(new char[m.size()+1]), off(0),
aux_i1(i1), aux_i2(i2), aux_i3(i3)
{
strncpy(mes, m.c_str(), m.size());
mes[m.size()] = '\0';
}
ParserException::ParserException(const ParserException& m, int plus_offset)
: mes(NULL),
aux_i1(-1), aux_i2(-1), aux_i3(-1)
{
copy(m);
off += plus_offset;
}
ParserException::ParserException(const ParserException& m, const char* dum, int i)
: mes(NULL),
aux_i1(-1), aux_i2(-1), aux_i3(-1)
{
copy(m);
aux_i3 = m.aux_i2;
aux_i2 = m.aux_i1;
aux_i1 = i;
}
ParserException::ParserException(const ParserException& m, const char* dum, int i1, int i2)
: mes(NULL),
aux_i1(-1), aux_i2(-1), aux_i3(-1)
{
copy(m);
aux_i3 = m.aux_i1;
aux_i2 = i2;
aux_i1 = i1;
}
ParserException::ParserException(const ParserException& m, const char* dum, int i1, int i2, int i3)
: mes(NULL),
aux_i1(-1), aux_i2(-1), aux_i3(-1)
{
copy(m);
aux_i3 = i3;
aux_i2 = i2;
aux_i1 = i1;
}
ParserException::ParserException(const ParserException& e)
: mes(NULL),
aux_i1(-1), aux_i2(-1), aux_i3(-1)
{
copy(e);
}
ParserException::~ParserException()
{
delete [] mes;
}
void ParserException::copy(const ParserException& e)
{
if (mes)
delete [] mes;
mes = new char[strlen(e.mes)+1];
strcpy(mes, e.mes);
off = e.off;
aux_i1 = e.aux_i1;
aux_i2 = e.aux_i2;
aux_i3 = e.aux_i3;
}
void ParserException::print(FILE* fd) const
{
// todo: to be refined
fprintf(fd, "%s: offset %d\n", mes, off);
}
dynare++/parser/cc/parser_exception.h deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: parser_exception.h 1761 2008-03-31 14:27:13Z kamenik $
#ifndef OG_FORMULA_PARSER_H
#define OG_FORMULA_PARSER_H
#include <string>
namespace ogp {
using std::string;
/** This is an easy exception, which, besides the message, stores
* also an offset of the parse error. Since we might need to track
* the argument number and for example the filed in the argument
* which caused the error, we add three integers, which have no
* semantics here. They should be documented in the function which
* throws an exception and sets them. Their default value is -1,
* which means they have not been set. */
class ParserException {
protected:
char* mes;
int off;
int aux_i1;
int aux_i2;
int aux_i3;
public:
ParserException(const char* m, int offset);
ParserException(const string& m, int offset);
ParserException(const string& m, const char* dum, int i1);
ParserException(const string& m, const char* dum, int i1, int i2);
ParserException(const string& m, const char* dum, int i1, int i2, int i3);
ParserException(const ParserException& e, int plus_offset);
/** Makes a copy and pushes given integer to aux_i1 shuffling
* others and forgetting the last. */
ParserException(const ParserException& e, const char* dum, int i);
/** Makes a copy and pushes given two integers to aux_i1 and aux_i2 shuffling
* others and forgetting the last two. */
ParserException(const ParserException& e, const char* dum, int i1, int i2);
/** Makes a copy and pushes given three integers to aux_i1, aux_i2, aus_i3 shuffling
* others and forgetting the last three. */
ParserException(const ParserException& e, const char* dum, int i1, int i2, int i3);
ParserException(const ParserException& e);
virtual ~ParserException();
void print(FILE* fd) const;
const char* message() const
{return mes;}
int offset() const
{return off;}
const int& i1() const
{return aux_i1;}
int& i1()
{return aux_i1;}
const int& i2() const
{return aux_i2;}
int& i2()
{return aux_i2;}
const int& i3() const
{return aux_i3;}
int& i3()
{return aux_i3;}
protected:
void copy(const ParserException& e);
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/parser/cc/static_atoms.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: static_atoms.cpp 1360 2007-07-10 11:44:20Z kamenik $
#include "static_atoms.h"
#include "utils/cc/exception.h"
using namespace ogp;
StaticAtoms::StaticAtoms(const StaticAtoms& a)
: Atoms(), Constants(a), varnames(a.varnames),
varorder(), vars(), indices()
{
// fill varorder
for (unsigned int i = 0; i < a.varorder.size(); i++) {
const char* s = varnames.query(a.varorder[i]);
varorder.push_back(s);
}
// fill vars
for (Tvarmap::const_iterator it = a.vars.begin();
it != a.vars.end(); ++it) {
const char* s = varnames.query((*it).first);
vars.insert(Tvarmap::value_type(s, (*it).second));
}
// fill indices
for (Tinvmap::const_iterator it = a.indices.begin();
it != a.indices.end(); ++it) {
const char* s = varnames.query((*it).second);
indices.insert(Tinvmap::value_type((*it).first, s));
}
}
void StaticAtoms::import_atoms(const DynamicAtoms& da, OperationTree& otree, Tintintmap& tmap)
{
Constants::import_constants(da, otree, tmap);
for (int i = 0; i < da.get_name_storage().num(); i++) {
const char* name = da.get_name_storage().get_name(i);
register_name(name);
int tnew = otree.add_nulary();
assign(name, tnew);
try {
const DynamicAtoms::Tlagmap& lmap = da.lagmap(name);
for (DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
it != lmap.end(); ++it) {
int told = (*it).second;
tmap.insert(Tintintmap::value_type(told, tnew));
}
} catch (const ogu::Exception& e) {
}
}
}
int StaticAtoms::check(const char* name) const
{
if (DynamicAtoms::is_string_constant(name)) {
return Constants::check(name);
} else {
return check_variable(name);
}
}
int StaticAtoms::index(const char* name) const
{
Tvarmap::const_iterator it = vars.find(name);
if (it == vars.end())
return -1;
else
return (*it).second;
}
const char* StaticAtoms::inv_index(int t) const
{
Tinvmap::const_iterator it = indices.find(t);
if (it == indices.end())
return NULL;
else
return (*it).second;
}
void StaticAtoms::assign(const char* name, int t)
{
if (DynamicAtoms::is_string_constant(name)) {
double val;
sscanf(name, "%lf", &val);
add_constant(t, val);
} else {
const char* ss = varnames.insert(name);
vars.insert(Tvarmap::value_type(ss, t));
indices.insert(Tinvmap::value_type(t, ss));
}
}
vector<int> StaticAtoms::variables() const
{
vector<int> res;
for (Tvarmap::const_iterator it = vars.begin();
it != vars.end(); ++it) {
res.push_back((*it).second);
}
return res;
}
void StaticAtoms::register_name(const char* name)
{
const char* ss = varnames.insert(name);
varorder.push_back(ss);
}
void StaticAtoms::print() const
{
printf("constants:\n");
Constants::print();
printf("variable names:\n");
varnames.print();
printf("map to tree indices:\n");
for (Tvarmap::const_iterator it = vars.begin(); it != vars.end(); ++it)
printf("%s\t->\t%d\n", (*it).first, (*it).second);
}
dynare++/parser/cc/static_atoms.h deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: static_atoms.h 1218 2007-03-19 21:52:49Z kamenik $
#ifndef OGP_STATIC_ATOMS
#define OGP_STATIC_ATOMS
#include "dynamic_atoms.h"
namespace ogp {
class StaticAtoms : public Atoms, public Constants {
protected:
typedef map<const char*, int, ltstr> Tvarmap;
typedef map<int, const char*> Tinvmap;
/** Storage for names. */
NameStorage varnames;
/** Outer order of variables. */
vector<const char*> varorder;
/** This is the map mapping a variable name to the tree
* index. */
Tvarmap vars;
/** This is the inverse mapping. It maps a tree index to the
* variable name. */
Tinvmap indices;
public:
StaticAtoms() : Atoms(), Constants(), varnames(), varorder(), vars()
{}
/* Copy constructor. */
StaticAtoms(const StaticAtoms& a);
/** Conversion from DynamicAtoms. This takes all atoms from
* the DynamicAtoms and adds its static version. The new tree
* indices are allocated in the passed OperationTree. Whole
* the process is traced in the map mapping old tree indices
* to new tree indices. */
StaticAtoms(const DynamicAtoms& da, OperationTree& otree, Tintintmap& tmap)
: Atoms(), Constants(), varnames(), varorder(), vars()
{import_atoms(da, otree, tmap);}
/* Destructor. */
virtual ~StaticAtoms() {}
/** This imports atoms from dynamic atoms inserting the new
* tree indices to the given tree (including constants). The
* mapping from old atoms to new atoms is traced in tmap. */
void import_atoms(const DynamicAtoms& da, OperationTree& otree,
Tintintmap& tmap);
/** If the name is constant, it returns its tree index if the
* constant is registered in Constants, it returns -1
* otherwise. If the name is not constant, it returns result
* from check_variable, which is implemented by a subclass. */
int check(const char* name) const;
/** This assigns a given tree index to the variable name. The
* name should have been checked before the call. */
void assign(const char* name, int t);
int nvar() const
{return varnames.num();}
/** This returns a vector of all variables. */
vector<int> variables() const;
/** This returns a tree index of the given variable. */
int index(const char* name) const;
/** This returns a name from the given tree index. NULL is
* returned if the tree index doesn't exist. */
const char* inv_index(int t) const;
/** This returns a name in a outer ordering. (There is no other ordering.) */
const char* name(int i) const
{return varorder[i];}
/** Debug print. */
void print() const;
/** This registers a variable. A subclass can reimplement
* this, for example, to ensure uniqueness of the
* name. However, this method should be always called in
* overriding methods to do the registering job. */
virtual void register_name(const char* name);
/** Return the name storage to allow querying to other
* classes. */
const NameStorage& get_name_storage() const
{return varnames;}
protected:
/** This checks the variable. The implementing subclass might
* want to throw an exception if the variable has not been
* registered. */
virtual int check_variable(const char* name) const = 0;
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/parser/cc/static_fine_atoms.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: static_fine_atoms.cpp 82 2007-04-19 11:33:30Z ondra $
#include "utils/cc/exception.h"
#include "static_fine_atoms.h"
#include "parser_exception.h"
using namespace ogp;
StaticFineAtoms::StaticFineAtoms(const StaticFineAtoms& sfa)
: StaticAtoms(sfa),
params(), param_outer_map(),
endovars(), endo_outer_map(),
exovars(), exo_outer_map(),
der_atoms(sfa.der_atoms),
endo_atoms_map(sfa.endo_atoms_map),
exo_atoms_map(sfa.exo_atoms_map)
{
for (unsigned int i = 0; i < sfa.params.size(); i++) {
const char* name = varnames.query(sfa.params[i]);
params.push_back(name);
param_outer_map.insert(Tvarintmap::value_type(name, i));
}
for (unsigned int i = 0; i < sfa.endovars.size(); i++) {
const char* name = varnames.query(sfa.endovars[i]);
endovars.push_back(name);
endo_outer_map.insert(Tvarintmap::value_type(name, i));
}
for (unsigned int i = 0; i < sfa.exovars.size(); i++) {
const char* name = varnames.query(sfa.exovars[i]);
exovars.push_back(name);
exo_outer_map.insert(Tvarintmap::value_type(name, i));
}
}
void StaticFineAtoms::import_atoms(const FineAtoms& fa, OperationTree& otree, Tintintmap& tmap)
{
StaticAtoms::import_atoms(fa, otree, tmap);
// we just need to put parameters, endovars, and exovars to
// respective vectors, the names are already in the storage
// parameters
const vector<const char*>& fa_params = fa.get_params();
for (unsigned int i = 0; i < fa_params.size(); i++)
register_param(fa_params[i]);
// endogenous
const vector<const char*>& fa_endovars = fa.get_endovars();
for (unsigned int i = 0; i < fa_endovars.size(); i++)
register_endo(fa_endovars[i]);
// exogenous
const vector<const char*>& fa_exovars = fa.get_exovars();
for (unsigned int i = 0; i < fa_exovars.size(); i++)
register_exo(fa_exovars[i]);
parsing_finished();
}
void StaticFineAtoms::import_atoms(const FineAtoms& fa, OperationTree& otree, Tintintmap& tmap,
const char* dummy)
{
StaticAtoms::import_atoms(fa, otree, tmap);
// we just need to put parameters, endovars, and exovars to
// respective vectors, the names are already in the storage
// parameters
const vector<const char*>& fa_params = fa.get_params();
for (unsigned int i = 0; i < fa_params.size(); i++)
register_param(fa_params[i]);
// endogenous
const vector<const char*>& fa_endovars = fa.get_endovars();
for (unsigned int i = 0; i < fa_endovars.size(); i++)
register_endo(fa_endovars[fa.y2outer_endo()[i]]);
// exogenous
const vector<const char*>& fa_exovars = fa.get_exovars();
for (unsigned int i = 0; i < fa_exovars.size(); i++)
register_exo(fa_exovars[fa.y2outer_exo()[i]]);
parsing_finished();
}
int StaticFineAtoms::check_variable(const char* name) const
{
const char* ss = varnames.query(name);
if (ss == NULL)
throw ParserException(string("Variable <")+name+"> not declared.",0);
return index(name);
}
void StaticFineAtoms::parsing_finished()
{
// build der_atoms, and endo_atoms_map and exo_atoms_map
der_atoms.clear();
endo_atoms_map.clear();
exo_atoms_map.clear();
// go through all endo and exo insert tree indices, ignore names
// whose tree index is -1 (those which are not referenced)
for (unsigned int i = 0; i < endovars.size(); i++) {
int t = index(endovars[i]);
if (t != -1) {
endo_atoms_map.push_back(der_atoms.size());
der_atoms.push_back(t);
}
}
for (unsigned int i = 0; i < exovars.size(); i++) {
int t = index(exovars[i]);
if (t != -1) {
exo_atoms_map.push_back(der_atoms.size());
der_atoms.push_back(t);
}
}
}
int StaticFineAtoms::name2outer_param(const char* name) const
{
Tvarintmap::const_iterator it = param_outer_map.find(name);
if (it == param_outer_map.end())
throw ogu::Exception(__FILE__,__LINE__,
"Name is not a parameter in StaticFineAtoms::name2outer_param");
return (*it).second;
}
int StaticFineAtoms::name2outer_endo(const char* name) const
{
Tvarintmap::const_iterator it = endo_outer_map.find(name);
if (it == endo_outer_map.end())
throw ogu::Exception(__FILE__,__LINE__,
"Name is not an endogenous variable in StaticFineAtoms::name2outer_endo");
return (*it).second;
}
int StaticFineAtoms::name2outer_exo(const char* name) const
{
Tvarintmap::const_iterator it = exo_outer_map.find(name);
if (it == exo_outer_map.end())
throw ogu::Exception(__FILE__,__LINE__,
"Name is not an exogenous variable in StaticFineAtoms::name2outer_exo");
return (*it).second;
}
void StaticFineAtoms::register_uniq_endo(const char* name)
{
if (varnames.query(name))
throw ogp::ParserException(string("Endogenous variable <")+name+"> is not unique.",0);
const char* ss = varnames.insert(name);
register_endo(ss);
}
void StaticFineAtoms::register_uniq_exo(const char* name)
{
if (varnames.query(name))
throw ogp::ParserException(string("Exogenous variable <")+name+"> is not unique.",0);
const char* ss = varnames.insert(name);
register_exo(ss);
}
void StaticFineAtoms::register_uniq_param(const char* name)
{
if (varnames.query(name))
throw ogp::ParserException(string("Parameter <")+name+"> is not unique.",0);
const char* ss = varnames.insert(name);
register_param(ss);
}
void StaticFineAtoms::print() const
{
StaticAtoms::print();
printf("endo atoms map:\n");
for (unsigned int i = 0; i < endo_atoms_map.size(); i++)
printf("%d --> %d\n", i, endo_atoms_map[i]);
printf("exo atoms map:\n");
for (unsigned int i = 0; i < exo_atoms_map.size(); i++)
printf("%d --> %d\n", i, exo_atoms_map[i]);
printf("der atoms:\n");
for (unsigned int i = 0; i < der_atoms.size(); i++)
printf("%d\t%d\n",i, der_atoms[i]);
}
void StaticFineAtoms::register_endo(const char* name)
{
const char* ss = varnames.query(name);
if (ss == NULL)
throw ogp::ParserException(string("Endogenous variable <")
+name+"> not found in storage.",0);
endovars.push_back(ss);
endo_outer_map.insert(Tvarintmap::value_type(ss, endovars.size()-1));
}
void StaticFineAtoms::register_exo(const char* name)
{
const char* ss = varnames.query(name);
if (ss == NULL)
throw ogp::ParserException(string("Exogenous variable <")
+name+"> not found in storage.",0);
exovars.push_back(ss);
exo_outer_map.insert(Tvarintmap::value_type(ss, exovars.size()-1));
}
void StaticFineAtoms::register_param(const char* name)
{
const char* ss = varnames.query(name);
if (ss == NULL)
throw ogp::ParserException(string("Parameter <")+name+"> not found in storage.",0);
params.push_back(ss);
param_outer_map.insert(Tvarintmap::value_type(ss, params.size()-1));
}
dynare++/parser/cc/static_fine_atoms.h deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: static_fine_atoms.h 42 2007-01-22 21:53:24Z ondra $
#ifndef OGP_STATIC_FINE_ATOMS_H
#define OGP_STATIC_FINE_ATOMS_H
#include "static_atoms.h"
#include "fine_atoms.h"
namespace ogp {
/** This class represents static atoms distinguishing between
* parameters, endogenous and exogenous variables. The class
* maintains also ordering of all three categories (referenced as
* outer or inner, since there is only one ordering). It can be
* constructed either from scratch, or from fine dynamic atoms. In
* the latter case, one can decide if the ordering of this static
* atoms should be internal or external ordering of the original
* dynamic fine atoms. */
class StaticFineAtoms : public StaticAtoms {
public:
typedef map<int,int> Tintintmap;
protected:
typedef map<const char*, int, ltstr> Tvarintmap;
private:
/** The vector of parameter names, gives the parameter
* ordering. */
vector<const char*> params;
/** A map mappping a parameter name to an index in the ordering. */
Tvarintmap param_outer_map;
/** The vector of endogenous variables. This defines the order
* like parameters. */
vector<const char*> endovars;
/** A map mapping a name of an endogenous variable to an index
* in the ordering. */
Tvarintmap endo_outer_map;
/** The vector of exogenous variables. Also defines the order
* like parameters and endovars. */
vector<const char*> exovars;
/** A map mapping a name of an exogenous variable to an index
* in the outer ordering. */
Tvarintmap exo_outer_map;
/** This vector defines a set of atoms as tree indices used
* for differentiation. The order of the atoms in is the
* concatenation of the outer ordering of endogenous and
* exogenous. This vector is setup by parsing_finished() and
* is returned by variables(). */
vector<int> der_atoms;
/** This is a mapping from endogenous atoms to all atoms in
* der_atoms member. The mapping maps index in endogenous atom
* ordering to index (not value) in der_atoms. It is useful if
* one wants to evaluate derivatives wrt only endogenous
* variables. It is set by parsing_finished(). By definition,
* it is monotone. */
vector<int> endo_atoms_map;
/** This is a mapping from exogenous atoms to all atoms in
* der_atoms member. It is the same as endo_atoms_map for
* atoms of exogenous variables. */
vector<int> exo_atoms_map;
public:
StaticFineAtoms() {}
/** Copy constructor making a new storage for atom names. */
StaticFineAtoms(const StaticFineAtoms& sfa);
/** Conversion from dynamic FineAtoms taking its outer
* ordering as ordering of parameters, endogenous and
* exogenous. A biproduct is an integer to integer map mapping
* tree indices of the dynamic atoms to tree indices of the
* static atoms. */
StaticFineAtoms(const FineAtoms& fa, OperationTree& otree, Tintintmap& tmap)
{StaticFineAtoms::import_atoms(fa, otree, tmap);}
/** Conversion from dynamic FineAtoms taking its internal
* ordering as ordering of parameters, endogenous and
* exogenous. A biproduct is an integer to integer map mapping
* tree indices of the dynamic atoms to tree indices of the
* static atoms. */
StaticFineAtoms(const FineAtoms& fa, OperationTree& otree, Tintintmap& tmap,
const char* dummy)
{StaticFineAtoms::import_atoms(fa, otree, tmap, dummy);}
virtual ~StaticFineAtoms() {}
/** This adds atoms from dynamic atoms inserting new tree
* indices to the given tree and tracing the mapping from old
* atoms to new atoms in tmap. The ordering of the static
* atoms is the same as outer ordering of dynamic atoms. */
void import_atoms(const FineAtoms& fa, OperationTree& otree, Tintintmap& tmap);
/** This adds atoms from dynamic atoms inserting new tree
* indices to the given tree and tracing the mapping from old
* atoms to new atoms in tmap. The ordering of the static
* atoms is the same as internal ordering of dynamic atoms. */
void import_atoms(const FineAtoms& fa, OperationTree& otree, Tintintmap& tmap,
const char* dummy);
/** Overrides StaticAtoms::check_variable so that the error
* would be raised if the variable name is not declared. A
* variable is declared by inserting it to
* StaticAtoms::varnames, which is done with registering
* methods. This a responsibility of a subclass. */
int check_variable(const char* name) const;
/** Return an (external) ordering of parameters. */
const vector<const char*>& get_params() const
{return params;}
/** Return an external ordering of endogenous variables. */
const vector<const char*>& get_endovars() const
{return endovars;}
/** Return an external ordering of exogenous variables. */
const vector<const char*>& get_exovars() const
{return exovars;}
/** This constructs der_atoms, and the endo_endoms_map and
* exo_atoms_map, which can be created only after the parsing
* is finished. */
void parsing_finished();
/** Return the atoms with respect to which we are going to
* differentiate. */
vector<int> variables() const
{return der_atoms;}
/** Return the endo_atoms_map. */
const vector<int>& get_endo_atoms_map() const
{return endo_atoms_map;}
/** Return the exo_atoms_map. */
const vector<int>& get_exo_atoms_map() const
{return endo_atoms_map;}
/** Return an index in the outer ordering of a given
* parameter. An exception is thrown if the name is not a
* parameter. */
int name2outer_param(const char* name) const;
/** Return an index in the outer ordering of a given
* endogenous variable. An exception is thrown if the name is not a
* and endogenous variable. */
int name2outer_endo(const char* name) const;
/** Return an index in the outer ordering of a given
* exogenous variable. An exception is thrown if the name is not a
* and exogenous variable. */
int name2outer_exo(const char* name) const;
/** Return the number of endogenous variables. */
int ny() const
{return endovars.size();}
/** Return the number of exogenous variables. */
int nexo() const
{return (int)exovars.size();}
/** Return the number of parameters. */
int np() const
{return (int)(params.size());}
/** Register unique endogenous variable name. The order of
* calls defines the endo outer ordering. The method is
* virtual, since a superclass may want to do some additional
* action. */
virtual void register_uniq_endo(const char* name);
/** Register unique exogenous variable name. The order of
* calls defines the exo outer ordering. The method is
* virtual, since a superclass may want to do somem additional
* action. */
virtual void register_uniq_exo(const char* name);
/** Register unique parameter name. The order of calls defines
* the param outer ordering. The method is
* virtual, since a superclass may want to do somem additional
* action. */
virtual void register_uniq_param(const char* name);
/** Debug print. */
void print() const;
private:
/** Add endogenous variable name, which is already in the name
* storage. */
void register_endo(const char* name);
/** Add exogenous variable name, which is already in the name
* storage. */
void register_exo(const char* name);
/** Add parameter name, which is already in the name
* storage. */
void register_param(const char* name);
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/parser/cc/tree.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2005-2011, Ondra Kamenik
#include "utils/cc/exception.h"
#include "tree.h"
#include <cstdlib>
#include <cmath>
#include <limits>
using namespace ogp;
/** Here we just implement complementary error function without
* declaring it for uses from outside this unit. The implementation is taken from "Numerical Recipes in C" 2nd ed. 1992 p. 221, */
double erffc(double x)
{
double z = std::abs(x);
double t = 1/(1+0.5*z);
double r = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))));
return x >= 0 ? r : 2-r;
}
/** Here we initialize OperationTree to contain only zero, one, nan
* and two_over_pi terms. */
OperationTree::OperationTree()
{
last_nulary = -1;
// allocate space for the constants
for (int i = 0; i < num_constants; i++)
add_nulary();
}
int OperationTree::add_nulary()
{
int op = terms.size();
Operation nulary;
terms.push_back(nulary);
_Tintset s;
s.insert(op);
nul_incidence.push_back(s);
_Tderivmap empty;
derivatives.push_back(empty);
last_nulary = op;
return op;
}
int OperationTree::add_unary(code_t code, int op)
{
if (op == zero &&
(code == UMINUS ||
code == SIN ||
code == TAN ||
code == SQRT ||
code == ERF))
return zero;
if (op == zero && code == LOG || op == nan)
return nan;
if (op == zero && (code == EXP ||
code == COS ||
code == ERFC))
return one;
Operation unary(code, op);
_Topmap::const_iterator i = ((const _Topmap&)opmap).find(unary);
if (i == opmap.end()) {
int newop = terms.size();
// add to the terms
terms.push_back(unary);
// copy incidence of the operand
nul_incidence.push_back(nul_incidence[op]);
// insert it to opmap
opmap.insert(_Topval(unary, newop));
// add empty map of derivatives
_Tderivmap empty;
derivatives.push_back(empty);
return newop;
}
return (*i).second;
}
int OperationTree::add_binary(code_t code, int op1, int op2)
{
// quick exits for special values
if (op1 == nan || op2 == nan)
return nan;
// for plus
if (code == PLUS)
if (op1 == zero && op2 == zero)
return zero;
else if (op1 == zero)
return op2;
else if (op2 == zero)
return op1;
// for minus
if (code == MINUS)
if (op1 == zero && op2 == zero)
return zero;
else if (op1 == zero)
return add_unary(UMINUS, op2);
else if (op2 == zero)
return op1;
// for times
if (code == TIMES)
if (op1 == zero || op2 == zero)
return zero;
else if (op1 == one)
return op2;
else if (op2 == one)
return op1;
// for divide
if (code == DIVIDE)
if (op1 == op2)
return one;
else if (op1 == zero)
return zero;
else if (op2 == zero)
return nan;
// for power
if (code == POWER)
if (op1 == zero && op2 == zero)
return nan;
else if (op1 == zero)
return zero;
else if (op2 == zero)
return one;
else if (op1 == one)
return one;
else if (op2 == one)
return op1;
// order operands of commutative operations
if (code == TIMES || code == PLUS)
if (op1 > op2) {
int tmp = op1;
op1 = op2;
op2 = tmp;
}
// construct operation and check/add it
Operation binary(code, op1, op2);
_Topmap::const_iterator i = ((const _Topmap&)opmap).find(binary);
if (i == opmap.end()) {
int newop = terms.size();
terms.push_back(binary);
// sum both sets of incidenting nulary operations
nul_incidence.push_back(nul_incidence[op1]);
nul_incidence.back().insert(nul_incidence[op2].begin(), nul_incidence[op2].end());
// add to opmap
opmap.insert(_Topval(binary, newop));
// add empty map of derivatives
_Tderivmap empty;
derivatives.push_back(empty);
return newop;
}
return (*i).second;
}
int OperationTree::add_derivative(int t, int v)
{
if (t < 0 || t >= (int) terms.size())
throw ogu::Exception(__FILE__,__LINE__,
"Wrong value for tree index in OperationTree::add_derivative");
// quick returns for nulary terms or empty incidence
if (terms[t].nary() == 0 && t != v) {
return zero;
}
if (terms[t].nary() == 0 && t == v) {
return one;
}
if (nul_incidence[t].end() == nul_incidence[t].find(v)) {
return zero;
}
// quick return if the derivative has been registered
_Tderivmap::const_iterator i = derivatives[t].find(v);
if (i != derivatives[t].end())
return (*i).second;
int res = -1;
switch (terms[t].getCode()) {
case UMINUS:
{
int tmp = add_derivative(terms[t].getOp1(), v);
res = add_unary(UMINUS, tmp);
break;
}
case LOG:
{
int tmp = add_derivative(terms[t].getOp1(), v);
res = add_binary(DIVIDE, tmp, terms[t].getOp1());
break;
}
case EXP:
{
int tmp = add_derivative(terms[t].getOp1(), v);
res = add_binary(TIMES, t, tmp);
break;
}
case SIN:
{
int tmp = add_derivative(terms[t].getOp1(), v);
res = add_binary(TIMES, add_unary(COS, terms[t].getOp1()), tmp);
break;
}
case COS:
{
int tmp = add_derivative(terms[t].getOp1(), v);
res = add_unary(UMINUS, add_binary(TIMES, add_unary(SIN, terms[t].getOp1()), tmp));
break;
}
case TAN:
{
int tmp = add_derivative(terms[t].getOp1(), v);
int tmp2 = add_unary(COS, terms[t].getOp1());
res = add_binary(DIVIDE, tmp, add_binary(TIMES, tmp2, tmp2));
break;
}
case SQRT:
{
int tmp = add_derivative(terms[t].getOp1(), v);
res = add_binary(DIVIDE, tmp,
add_binary(PLUS, t, t));
break;
}
case ERF:
{
int tmp = add_binary(TIMES, terms[t].getOp1(), terms[t].getOp1());
tmp = add_unary(UMINUS, tmp);
tmp = add_unary(EXP, tmp);
int der = add_derivative(terms[t].getOp1(), v);
tmp = add_binary(TIMES, tmp, der);
res = add_binary(TIMES, two_over_pi, tmp);
break;
}
case ERFC:
{
int tmp = add_binary(TIMES, terms[t].getOp1(), terms[t].getOp1());
tmp = add_unary(UMINUS, tmp);
tmp = add_unary(EXP, tmp);
int der = add_derivative(terms[t].getOp1(), v);
tmp = add_binary(TIMES, tmp, der);
tmp = add_binary(TIMES, two_over_pi, tmp);
res = add_unary(UMINUS, tmp);
break;
}
case PLUS:
{
int tmp1 = add_derivative(terms[t].getOp1(), v);
int tmp2 = add_derivative(terms[t].getOp2(), v);
res = add_binary(PLUS, tmp1, tmp2);
break;
}
case MINUS:
{
int tmp1 = add_derivative(terms[t].getOp1(), v);
int tmp2 = add_derivative(terms[t].getOp2(), v);
res = add_binary(MINUS, tmp1, tmp2);
break;
}
case TIMES:
{
int tmp1 = add_derivative(terms[t].getOp1(), v);
int tmp2 = add_derivative(terms[t].getOp2(), v);
int res1 = add_binary(TIMES, terms[t].getOp1(), tmp2);
int res2 = add_binary(TIMES, tmp1, terms[t].getOp2());
res = add_binary(PLUS, res1, res2);
break;
}
case DIVIDE:
{
int tmp1 = add_derivative(terms[t].getOp1(), v);
int tmp2 = add_derivative(terms[t].getOp2(), v);
if (tmp2 == zero)
res = add_binary(DIVIDE, tmp1, terms[t].getOp2());
else {
int nom = add_binary(MINUS,
add_binary(TIMES, tmp1, terms[t].getOp2()),
add_binary(TIMES, tmp2, terms[t].getOp1()));
int den = add_binary(TIMES, terms[t].getOp2(), terms[t].getOp2());
res = add_binary(DIVIDE, nom, den);
}
break;
}
case POWER:
{
int tmp1 = add_derivative(terms[t].getOp1(), v);
int tmp2 = add_derivative(terms[t].getOp2(), v);
int s1 = add_binary(TIMES, tmp2,
add_binary(TIMES, t,
add_unary(LOG, terms[t].getOp1())));
int s2 = add_binary(TIMES, tmp1,
add_binary(TIMES, terms[t].getOp2(),
add_binary(POWER, terms[t].getOp1(),
add_binary(MINUS, terms[t].getOp2(), one))));
res = add_binary(PLUS, s1, s2);
break;
}
case NONE:
break;
}
if (res == -1)
throw ogu::Exception(__FILE__,__LINE__,
"Unknown operation code.");
register_derivative(t, v, res);
return res;
}
int OperationTree::add_substitution(int t, const map<int,int>& subst)
{
return add_substitution(t, subst, *this);
}
int OperationTree::add_substitution(int t, const map<int,int>& subst,
const OperationTree& otree)
{
// return substitution of t if it is in the map
map<int,int>::const_iterator it = subst.find(t);
if (subst.end() != it)
return (*it).second;
int nary = otree.terms[t].nary();
if (nary == 2) {
// return the binary operation of the substituted terms
int t1 = add_substitution(otree.terms[t].getOp1(), subst, otree);
int t2 = add_substitution(otree.terms[t].getOp2(), subst, otree);
return add_binary(otree.terms[t].getCode(), t1, t2);
} else if (nary == 1) {
// return the unary operation of the substituted term
int t1 = add_substitution(otree.terms[t].getOp1(), subst, otree);
return add_unary(otree.terms[t].getCode(), t1);
} else {
// if t is not the first num_constants, and otree is not this
// tree, then raise and exception. Otherwise return t, since
// it is either a special term (having the same semantics in
// both trees), or the trees are the same, hence t has the
// same semantics
if (t < num_constants || this == &otree)
return t;
else {
throw ogu::Exception(__FILE__,__LINE__,
"Incomplete substitution map in OperationTree::add_substitution");
return -1;
}
}
}
void OperationTree::nularify(int t)
{
// remove the original operation from opmap
_Topmap::iterator it = opmap.find(terms[t]);
if (it != opmap.end())
opmap.erase(it);
// turn the operation to nulary
Operation nulary_op;
terms[t] = nulary_op;
// update last nulary
if (last_nulary < t)
last_nulary = t;
// update nul_incidence information for all terms including t
update_nul_incidence_after_nularify(t);
}
void OperationTree::register_derivative(int t, int v, int tder)
{
// todo: might check that the insert inserts a new pair
derivatives[t].insert(_Tderivmap::value_type(v, tder));
}
unordered_set<int> OperationTree::select_terms(int t, const opselector& sel) const
{
unordered_set<int> subterms;
select_terms(t, sel, subterms);
return subterms;
}
void OperationTree::select_terms(int t, const opselector& sel, unordered_set<int>& subterms) const
{
const Operation& op = terms[t];
if (sel(t))
subterms.insert(t);
else
if (op.nary() == 2) {
select_terms(op.getOp1(), sel, subterms);
select_terms(op.getOp2(), sel, subterms);
} else if (op.nary() == 1) {
select_terms(op.getOp1(), sel, subterms);
}
}
unordered_set<int> OperationTree::select_terms_inv(int t, const opselector& sel) const
{
unordered_set<int> subterms;
select_terms_inv(t, sel, subterms);
return subterms;
}
bool OperationTree::select_terms_inv(int t, const opselector& sel, unordered_set<int>& subterms) const
{
const Operation& op = terms[t];
if (op.nary() == 2) {
bool a1 = select_terms_inv(op.getOp1(), sel, subterms);
bool a2 = select_terms_inv(op.getOp2(), sel, subterms);
if (a1 && a2 && sel(t)) {
subterms.insert(t);
return true;
}
} else if (op.nary() == 1) {
bool a1 = select_terms_inv(op.getOp1(), sel, subterms);
if (a1 && sel(t)) {
subterms.insert(t);
return true;
}
} else {
if (sel(t)) {
subterms.insert(t);
return true;
}
}
return false;
}
void OperationTree::forget_derivative_maps()
{
for (unsigned int i = 0; i < derivatives.size(); i++)
derivatives[i].clear();
}
void OperationTree::print_operation_tree(int t, FILE* fd, OperationFormatter& f) const
{
f.format(terms[t], t, fd);
}
void OperationTree::print_operation(int t) const
{
DefaultOperationFormatter dof(*this);
print_operation_tree(t, stdout, dof);
}
void OperationTree::update_nul_incidence_after_nularify(int t)
{
unordered_set<int> updated;
for (int tnode = num_constants; tnode < (int)terms.size(); tnode++) {
const Operation& op = terms[tnode];
if (op.nary() == 2) {
int op1 = op.getOp1();
int op2 = op.getOp2();
if (op1 >= tnode || op2 >= tnode)
throw ogu::Exception(__FILE__,__LINE__,
"Tree disorder asserted");
bool updated1 = (updated.end() != updated.find(op1));
bool updated2 = (updated.end() != updated.find(op2));
if (updated1 || updated2) {
nul_incidence[tnode] = nul_incidence[op1];
nul_incidence[tnode].insert(nul_incidence[op2].begin(), nul_incidence[op2].end());
updated.insert(tnode);
}
} else if (op.nary() == 1) {
int op1 = op.getOp1();
if (op1 >= tnode)
throw ogu::Exception(__FILE__,__LINE__,
"Tree disorder asserted");
bool updated1 = (updated.end() != updated.find(op1));
if (updated1) {
nul_incidence[tnode] = nul_incidence[op1];
updated.insert(tnode);
}
} else if (op.nary() == 0) {
if (tnode == t) {
nul_incidence[tnode].clear();
nul_incidence[tnode].insert(tnode);
updated.insert(tnode);
}
}
}
}
EvalTree::EvalTree(const OperationTree& ot, int last)
: otree(ot),
values(new double[(last==-1)? ot.terms.size() : last+1]),
flags(new bool[(last==-1)? ot.terms.size() : last+1]),
last_operation((last==-1)? ot.terms.size()-1 : last)
{
if (last_operation < OperationTree::num_constants-1 ||
last_operation > (int)ot.terms.size()-1)
throw ogu::Exception(__FILE__,__LINE__,
"Wrong last in EvalTree constructor.");
values[0] = 0.0;
flags[0] = true;
values[1] = 1.0;
flags[1] = true;
values[2] = std::numeric_limits<double>::quiet_NaN();
flags[2] = true;
values[3] = 2.0/sqrt(M_PI);
flags[3] = true;
// this sets from num_constants on
reset_all();
}
void EvalTree::reset_all()
{
for (int i = OperationTree::num_constants; i <= last_operation; i++)
flags[i] = false;
}
void EvalTree::set_nulary(int t, double val)
{
if (t < 0 || t > last_operation)
throw ogu::Exception(__FILE__,__LINE__,
"The tree index out of bounds in EvalTree::set_nulary");
if (t < OperationTree::num_constants || otree.terms[t].nary() != 0)
throw ogu::Exception(__FILE__,__LINE__,
"The term is not nulary assignable in EvalTree::set_nulary");
values[t] = val;
flags[t] = true;
}
double EvalTree::eval(int t)
{
if (t < 0 || t > last_operation)
throw ogu::Exception(__FILE__,__LINE__,
"The tree index out of bounds in EvalTree::eval");
if (otree.terms[t].nary() == 0 && flags[t] == false)
throw ogu::Exception(__FILE__,__LINE__,
"Nulary term has not been assigned a value in EvalTree::eval");
if (! flags[t]) {
const Operation& op = otree.terms[t];
if (op.nary() == 1) {
double r1 = eval(op.getOp1());
double res;
if (op.getCode() == UMINUS)
res = -r1;
else if (op.getCode() == LOG)
res = log(r1);
else if (op.getCode() == EXP)
res = exp(r1);
else if (op.getCode() == SIN)
res = sin(r1);
else if (op.getCode() == COS)
res = cos(r1);
else if (op.getCode() == TAN)
res = tan(r1);
else if (op.getCode() == SQRT)
res = sqrt(r1);
else if (op.getCode() == ERF)
res = 1-erffc(r1);
else if (op.getCode() == ERFC)
res = erffc(r1);
else {
throw ogu::Exception(__FILE__,__LINE__,
"Unknown unary operation code in EvalTree::eval");
res = 0.0;
}
values[t] = res;
flags[t] = true;
} else if (op.nary() == 2) {
double res;
if (op.getCode() == PLUS) {
double r1 = eval(op.getOp1());
double r2 = eval(op.getOp2());
res = r1 + r2;
} else if (op.getCode() == MINUS) {
double r1 = eval(op.getOp1());
double r2 = eval(op.getOp2());
res = r1 - r2;
} else if (op.getCode() == TIMES) {
// pickup less complex formula first
unsigned int nul1 = otree.nulary_of_term(op.getOp1()).size();
unsigned int nul2 = otree.nulary_of_term(op.getOp2()).size();
if (nul1 < nul2) {
double r1 = eval(op.getOp1());
if (r1 == 0.0)
res = 0.0;
else {
double r2 = eval(op.getOp2());
res = r1 * r2;
}
} else {
double r2 = eval(op.getOp2());
if (r2 == 0)
res = 0.0;
else {
double r1 = eval(op.getOp1());
res = r1*r2;
}
}
} else if (op.getCode() == DIVIDE) {
double r1 = eval(op.getOp1());
if (r1 == 0)
res = 0.0;
else {
double r2 = eval(op.getOp2());
res = r1 / r2;
}
} else if (op.getCode() == POWER) {
// suppose that more complex is the first op in average
double r2 = eval(op.getOp2());
if (r2 == 0.0)
res = 1.0;
else {
double r1 = eval(op.getOp1());
res = pow(r1, r2);
}
} else {
throw ogu::Exception(__FILE__,__LINE__,
"Unknown binary operation code in EvalTree::eval");
res = 0.0;
}
values[t] = res;
flags[t] = true;
}
return values[t];
}
// if (! std::isfinite(values[t]))
// printf("Tree value t=%d is not finite = %f\n", t, values[t]);
return values[t];
}
void EvalTree::print() const
{
printf("last_op=%d\n", last_operation);
printf(" 0 1 2 3 4 5 6 7 8 9\n");
printf("----------------------------------------------------------------\n");
for (int i = 0; i <= (last_operation+1)/10; i++) {
printf("%-3d|", i);
int j = 0;
while (j < 10 && 10*i+j < last_operation+1) {
int k = 10*i+j;
if (flags[k])
printf(" %5.1g", values[k]);
else
printf(" -----");
j++;
}
printf("\n");
}
}
void DefaultOperationFormatter::format(const Operation& op, int t, FILE* fd)
{
// add to the stop_set
if (stop_set.end() == stop_set.find(t))
stop_set.insert(t);
else
return;
// call recursively non-nulary terms of the operation
if (op.nary() == 2) {
int t1 = op.getOp1();
const Operation& op1 = otree.terms[t1];
int t2 = op.getOp2();
const Operation& op2 = otree.terms[t2];
if (op1.nary() > 0)
format(op1, t1, fd);
if (op2.nary() > 0)
format(op2, t2, fd);
}
if (op.nary() == 1) {
int t1 = op.getOp1();
const Operation& op1 = otree.terms[t1];
if (op1.nary() > 0)
format(op1, t1, fd);
}
// print 'term ='
format_term(t, fd);
fprintf(fd, " = ");
if (op.nary() == 0) {
format_nulary(t, fd);
} else if (op.nary() == 1) {
int t1 = op.getOp1();
const Operation& op1 = otree.terms[t1];
const char* opname = "unknown";
switch (op.getCode()) {
case UMINUS:
opname = "-";
break;
case LOG:
opname = "log";
break;
case EXP:
opname = "exp";
break;
case SIN:
opname = "sin";
break;
case COS:
opname = "cos";
break;
case TAN:
opname = "tan";
break;
case SQRT:
opname = "sqrt";
break;
case ERF:
opname = "erf";
break;
case ERFC:
opname = "erfc";
break;
default:
break;
}
fprintf(fd, "%s(", opname);
if (op1.nary() == 0)
format_nulary(t1, fd);
else
format_term(t1, fd);
fprintf(fd, ")");
} else {
int t1 = op.getOp1();
const Operation& op1 = otree.terms[t1];
int t2 = op.getOp2();
const Operation& op2 = otree.terms[t2];
const char* opname = "unknown";
switch (op.getCode()) {
case PLUS:
opname = "+";
break;
case MINUS:
opname = "-";
break;
case TIMES:
opname = "*";
break;
case DIVIDE:
opname = "/";
break;
case POWER:
opname = "^";
break;
default:
break;
}
if (op1.nary() == 0)
format_nulary(t1, fd);
else
format_term(t1, fd);
fprintf(fd, " %s ", opname);
if (op2.nary() == 0)
format_nulary(t2, fd);
else
format_term(t2, fd);
}
print_delim(fd);
}
void DefaultOperationFormatter::format_term(int t, FILE* fd) const
{
fprintf(fd, "$%d", t);
}
void DefaultOperationFormatter::format_nulary(int t, FILE* fd) const
{
if (t == OperationTree::zero)
fprintf(fd, "0");
else if (t == OperationTree::one)
fprintf(fd, "1");
else if (t == OperationTree::nan)
fprintf(fd, "NaN");
else
fprintf(fd, "$%d", t);
}
void DefaultOperationFormatter::print_delim(FILE* fd) const
{
fprintf(fd, ";\n");
}
std::string OperationStringConvertor::convert(const Operation& op, int t) const
{
if (op.nary() == 0) {
if (t < OperationTree::num_constants)
if (t == OperationTree::zero)
return std::string("0");
else if (t == OperationTree::one)
return std::string("1");
else if (t == OperationTree::nan)
return std::string("NaN");
else if (t == OperationTree::two_over_pi) {
char buf[100];
sprintf(buf, "%20.16g", 2.0/std::sqrt(M_PI));
return std::string(buf);
} else {
return std::string("error!error");
}
else
return nulsc.convert(t);
} else if (op.nary() == 1) {
int t1 = op.getOp1();
const Operation& op1 = otree.operation(t1);
const char* opname = "unknown";
switch (op.getCode()) {
case UMINUS:
opname = "-";
break;
case LOG:
opname = "log";
break;
case EXP:
opname = "exp";
break;
case SIN:
opname = "sin";
break;
case COS:
opname = "cos";
break;
case TAN:
opname = "tan";
break;
case SQRT:
opname = "sqrt";
break;
case ERF:
opname = "erf";
break;
case ERFC:
opname = "erfc";
break;
default:
break;
}
std::string s1 = convert(op1, t1);
return std::string(opname) + "(" + s1 + ")";
} else {
int t1 = op.getOp1();
const Operation& op1 = otree.operation(t1);
int t2 = op.getOp2();
const Operation& op2 = otree.operation(t2);
const char* opname = "unknown";
switch (op.getCode()) {
case PLUS:
opname = "+";
break;
case MINUS:
opname = "-";
break;
case TIMES:
opname = "*";
break;
case DIVIDE:
opname = "/";
break;
case POWER:
opname = "^";
break;
default:
break;
}
// decide about parenthesis
bool op1_par = true;
bool op2_par = true;
if (op.getCode() == PLUS) {
op1_par = false;
op2_par = false;
} else if (op.getCode() == MINUS) {
op1_par = false;
if (op2.getCode() != MINUS && op2.getCode() != PLUS)
op2_par = false;
} else {
if (op1.nary() < 2)
op1_par = false;
if (op2.nary() < 2)
op2_par = false;
}
std::string res;
if (op1_par)
res += "(";
res += convert(op1, t1);
if (op1_par)
res += ")";
res += " ";
res += opname;
res += " ";
if (op2_par)
res += "(";
res += convert(op2, t2);
if (op2_par)
res += ")";
return res;
}
}
// Local Variables:
// mode:C++
// End:
dynare++/parser/cc/tree.h deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2005-2011, Ondra Kamenik
#ifndef OGP_TREE_H
#define OGP_TREE_H
#include <vector>
#include <set>
#include <map>
#include <boost/unordered_map.hpp>
#include <boost/unordered_set.hpp>
#include <cstdio>
namespace ogp {
using boost::unordered_set;
using boost::unordered_map;
using std::vector;
using std::set;
using std::map;
/** Enumerator representing nulary, unary and binary operation
* codes. For nulary, 'none' is used. When one is adding a new
* codes, he should update the code of #OperationTree::add_unary,
* #OperationTree::add_binary, and of course
* #OperationTree::add_derivative. */
enum code_t {NONE, UMINUS, LOG, EXP, SIN, COS, TAN, SQRT, ERF,
ERFC, PLUS, MINUS, TIMES, DIVIDE, POWER};
/** Class representing a nulary, unary, or binary operation. */
class Operation {
protected:
/** Code of the operation. */
code_t code;
/** First operand. If none, then it is -1. */
int op1;
/** Second operand. If none, then it is -1. */
int op2;
public:
/** Constructs a binary operation. */
Operation(code_t cd, int oper1, int oper2)
: code(cd), op1(oper1), op2(oper2) {}
/** Constructs a unary operation. */
Operation(code_t cd, int oper1)
: code(cd), op1(oper1), op2(-1) {}
/** Constructs a nulary operation. */
Operation()
: code(NONE), op1(-1), op2(-1) {}
/** A copy constructor. */
Operation(const Operation& op)
: code(op.code), op1(op.op1), op2(op.op2) {}
/** Operator =. */
const Operation& operator=(const Operation& op)
{
code = op.code;
op1 = op.op1;
op2 = op.op2;
return *this;
}
/** Operator ==. */
bool operator==(const Operation& op) const
{
return code == op.code && op1 == op.op1 && op2 == op.op2;
}
/** Operator < implementing lexicographic ordering. */
bool operator<(const Operation& op) const
{
return (code < op.code ||
code == op.code &&
(op1 < op.op1 || op1 == op.op1 && op2 < op.op2));
}
/** Returns a number of operands. */
int nary() const
{
return (op2 == -1)? ((op1 == -1) ? 0 : 1) : 2;
}
/** Returns a hash value of the operation. */
size_t hashval() const
{
return op2+1 + (op1+1)^15 + code^30;
}
code_t getCode() const
{ return code; }
int getOp1() const
{ return op1; }
int getOp2() const
{ return op2; }
};
/** This struct is a predicate for ordering of the operations in
* OperationTree class. now obsolete */
struct ltoper {
bool operator()(const Operation& oper1, const Operation& oper2) const
{return oper1 < oper2;}
};
/** Hash function object for Operation. */
struct ophash {
size_t operator()(const Operation& op) const
{ return op.hashval(); }
};
/** This struct is a function object selecting some
* operations. The operation is given by a tree index. */
struct opselector {
virtual bool operator()(int t) const = 0;
virtual ~opselector() {}
};
/** Forward declaration of OperationFormatter. */
class OperationFormatter;
class DefaultOperationFormatter;
/** Forward declaration of EvalTree to make it friend of OperationTree. */
class EvalTree;
/** Class representing a set of trees for terms. Each term is
* given a unique non-negative integer. The terms are basically
* operations whose (integer) operands point to another terms in
* the tree. The terms are stored in the vector. Equivalent unary
* and binary terms are stored only once. This class guarantees
* the uniqueness. The uniqueness of nulary terms is guaranteed by
* the caller, since at this level of Operation abstraction, one
* cannot discriminate between different nulary operations
* (constants, variables). The uniqueness is enforced by the
* unordered_map whose keys are operations and values are integers
* (indices of the terms).
* This class can also make derivatives of a given term with
* respect to a given nulary term. I order to be able to quickly
* recognize zero derivativates, we maintain a list of nulary
* terms contained in the term. A possible zero derivative is then quickly
* recognized by looking at the list. The list is implemented as a
* unordered_set of integers.
*
* In addition, many term can be differentiated multiple times wrt
* one variable since they can be referenced multiple times. To
* avoid this, for each term we maintain a map mapping variables
* to the derivatives of the term. As the caller will
* differentiate wrt more and more variables, these maps will
* become richer and richer.
*/
class OperationTree {
friend class EvalTree;
friend class DefaultOperationFormatter;
protected:
/** This is the vector of the terms. An index to this vector
* uniquelly determines the term. */
vector<Operation> terms;
/** This defines a type for a map mapping the unary and binary
* operations to their indices. */
typedef unordered_map<Operation, int, ophash> _Topmap;
typedef _Topmap::value_type _Topval;
/** This is the map mapping the unary and binary operations to
* the indices of the terms.*/
_Topmap opmap;
/** This is a type for a set of integers. */
typedef unordered_set<int> _Tintset;
/** This is a vector of integer sets corresponding to the
* nulary terms contained in the term. */
vector<_Tintset> nul_incidence;
/** This is a type of the map from variables (nulary terms) to
* the terms. */
typedef unordered_map<int, int> _Tderivmap;
/** This is a vector of derivative mappings. For each term, it
* maps variables to the derivatives of the term with respect
* to the variables. */
vector<_Tderivmap> derivatives;
/** The tree index of the last nulary term. */
int last_nulary;
public:
/** This is a number of constants set in the following
* enum. This number reserves space in a vector of terms for
* the constants. */
static const int num_constants = 4;
/** Enumeration for special terms. We need zero, one, nan and
* 2/pi. These will be always first four terms having indices
* zero, one and two, three. If adding anything to this
* enumeration, make sure you have updated num_constants above.*/
enum {zero=0, one=1, nan=2, two_over_pi=3};
/** The unique constructor which initializes the object to
* contain only zero, one and nan and two_over_pi.*/
OperationTree();
/** Copy constructor. */
OperationTree(const OperationTree& ot)
: terms(ot.terms), opmap(ot.opmap), nul_incidence(ot.nul_incidence),
derivatives(ot.derivatives),
last_nulary(ot.last_nulary)
{}
/** Add a nulary operation. The caller is responsible for not
* inserting two semantically equivalent nulary operations.
* @return newly allocated index
*/
int add_nulary();
/** Add a unary operation. The uniqness is checked, if it
* already exists, then it is not added.
* @param code the code of the unary operation
* @param op the index of the operand
* @return the index of the operation
*/
int add_unary(code_t code, int op);
/** Add a binary operation. The uniqueness is checked, if it
* already exists, then it is not added.
* @param code the code of the binary operation
* @param op1 the index of the first operand
* @param op2 the index of the second operand
* @return the index of the operation
*/
int add_binary(code_t code, int op1, int op2);
/** Add the derivative of the given term with respect to the
* given nulary operation.
* @param t the index of the operation being differentiated
* @param v the index of the nulary operation
* @return the index of the derivative
*/
int add_derivative(int t, int v);
/** Add the substitution given by the map. This adds a new
* term which is equal to the given term with applied
* substitutions given by the map replacing each term on the
* left by a term on the right. We do not check that the terms
* on the left are not subterms of the terms on the right. If
* so, the substituted terms are not subject of further
* substitution. */
int add_substitution(int t, const map<int,int>& subst);
/** Add the substitution given by the map where left sides of
* substitutions come from another tree. The right sides are
* from this tree. The given t is from the given otree. */
int add_substitution(int t, const map<int,int>& subst,
const OperationTree& otree);
/** This method turns the given term to a nulary
* operation. This is an only method, which changes already
* existing term (all other methods add something new). User
* should use this with caution and must make sure that
* something similar has happened for atoms. In addition, it
* does not do anything with derivatives, so it should not be
* used after some derivatives were created, and derivatives
* already created and saved in derivatives mappings should be
* forgotten with forget_derivative_maps. */
void nularify(int t);
/** Return the set of nulary terms of the given term. */
const unordered_set<int>& nulary_of_term(int t) const
{return nul_incidence[t];}
/** Select subterms of the given term according a given
* operation selector and return the set of terms that
* correspond to the compounded operations. The given term is
* a compound function of the returned subterms and the
* function consists only from operations which yield false in
* the selector. */
unordered_set<int> select_terms(int t, const opselector& sel) const;
/** Select subterms of the given term according a given
* operation selector and return the set of terms that
* correspond to the compounded operations. The given term is
* a compound function of the returned subterms and the
* subterms are maximal subterms consisting from operations
* yielding true in the selector. */
unordered_set<int> select_terms_inv(int t, const opselector& sel) const;
/** This forgets all the derivative mappings. It is used after
* a term has been nularified, and then the derivative
* mappings carry wrong information. Note that the derivatives
* mappings serve only as a tool for quick returns in
* add_derivative. Resseting the mappings is harmless, all the
* information is rebuilt in add_derivative without any
* additional nodes (trees). */
void forget_derivative_maps();
/** This returns an operation of a given term. */
const Operation& operation(int t) const
{return terms[t];}
/** This outputs the operation to the given file descriptor
* using the given OperationFormatter. */
void print_operation_tree(int t, FILE* fd, OperationFormatter& f) const;
/** Debug print of a given operation: */
void print_operation(int t) const;
/** Return the last tree index of a nulary term. */
int get_last_nulary() const
{return last_nulary;}
/** Get the number of all operations. */
int get_num_op() const
{return (int)(terms.size());}
private:
/** This registers a calculated derivative of the term in the
* #derivatives vector.
* @param t the index of the term for which we register the derivative
* @param v the index of the nulary term (variable) to which
* respect the derivative was taken
* @param tder the index of the resulting derivative
*/
void register_derivative(int t, int v, int tder);
/** This does the same job as select_terms with the only
* difference, that it adds the terms to the given set and
* hence can be used recursivelly. */
void select_terms(int t, const opselector& sel, unordered_set<int>& subterms) const;
/** This does the same job as select_terms_inv with the only
* difference, that it adds the terms to the given set and
* hence can be used recursivelly and returns true if the term
* was selected. */
bool select_terms_inv(int t, const opselector& sel, unordered_set<int>& subterms) const;
/** This updates nul_incidence information after the term t
* was turned to a nulary term in all terms. It goes through
* the tree from simplest terms to teh more complex ones and
* changes the nul_incidence information where necesary. It
* maintains a set where the changes have been made.*/
void update_nul_incidence_after_nularify(int t);
};
/** EvalTree class allows for an evaluation of the given tree for
* a given values of nulary terms. For each term in the
* OperationTree the class maintains a resulting value and a flag
* if the value has been calculated or set. The life cycle of the
* class is the following: After it is initialized, the user must
* set values for necessary nulary terms. Then the object can be
* requested to evaluate particular terms. During this process,
* the number of evaluated terms is increasing. Then the user can
* request overall reset of evaluation flags, set the nulary terms
* to new values and evaluate a number of terms.
*
* Note that currently the user cannot request a reset of
* evaluation flags only for those terms depending on a given
* nulary term. This might be added in future and handeled by a
* subclasses of OperationTree and EvalTree, since we need a
* support for this in OperationTree.
*/
class EvalTree {
protected:
/** Reference to the OperationTree over which all evaluations
* are done. */
const OperationTree& otree;
/** The array of values. */
double* const values;
/** The array of evaluation flags. */
bool* const flags;
/** The index of last operation in the EvalTree. Length of
* values and flags will be then last_operation+1. */
int last_operation;
public:
/** Initializes the evaluation tree for the given operation
* tree. If last is greater than -1, that the evaluation tree
* will contain only formulas up to the given last index
* (included). */
EvalTree(const OperationTree& otree, int last = -1);
virtual ~EvalTree()
{ delete [] values; delete [] flags; }
/** Set evaluation flag to all terms (besides the first
* special terms) to false. */
void reset_all();
/** Set value for a given nulary term. */
void set_nulary(int t, double val);
/** Evaluate the given term with nulary terms set so far. */
double eval(int t);
/** Debug print. */
void print() const;
/* Return the operation tree. */
const OperationTree& getOperationTree() const
{return otree;}
private:
EvalTree(const EvalTree&);
};
/** This is an interface describing how a given operation is
* formatted for output. */
class OperationFormatter {
public:
/** Empty virtual destructor. */
virtual ~OperationFormatter() {}
/** Print the formatted operation op with a given tree index t
* to a given descriptor. (See class OperationTree to know
* what is a tree index.) This prints all the tree. This
* always writes equation, left hand side is a string
* represenation (a variable, temporary, whatever) of the
* term, the right hand side is a string representation of the
* operation (which will refer to other string representation
* of subterms). */
virtual void format(const Operation& op, int t, FILE* fd)=0;
};
/** The default formatter formats the formulas with a usual syntax
* (for example Matlab). A formatting of atoms and terms might be
* reimplemented by a subclass. In addition, during its life, the
* object maintains a set of tree indices which have been output
* and they are not output any more. */
class DefaultOperationFormatter : public OperationFormatter {
protected:
const OperationTree& otree;
set<int> stop_set;
public:
DefaultOperationFormatter(const OperationTree& ot)
: otree(ot) {}
/** Format the operation with the default syntax. */
void format(const Operation& op, int t, FILE* fd);
/** This prints a string represenation of the given term, for
* example 'tmp10' for term 10. In this implementation it
* prints $10. */
virtual void format_term(int t, FILE* fd) const;
/** Print a string representation of the nulary term. */
virtual void format_nulary(int t, FILE* fd) const;
/** Print a delimiter between two statements. By default it is
* "\n". */
virtual void print_delim(FILE* fd) const;
};
class NularyStringConvertor {
public:
virtual ~NularyStringConvertor() {}
/** Return the string representation of the atom with the tree
* index t. */
virtual std::string convert(int t) const = 0;
};
/** This class converts the given term to its mathematical string representation. */
class OperationStringConvertor {
protected:
const NularyStringConvertor& nulsc;
const OperationTree& otree;
public:
OperationStringConvertor(const NularyStringConvertor& nsc, const OperationTree& ot)
: nulsc(nsc), otree(ot) {}
/** Empty virtual destructor. */
virtual ~OperationStringConvertor() {}
/** Convert the operation to the string mathematical
* representation. This does not write any equation, just
* returns a string representation of the formula. */
std::string convert(const Operation& op, int t) const;
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/src/Makefile.am deleted 100644 → 0
View file @ 6ca649c4
bin_PROGRAMS = dynare++
GENERATED_FILES = dynglob_ll.cc dynglob_tab.cc dynglob_tab.hh
dynare___SOURCES = \
main.cpp \
dynare3.cpp \
dynare_atoms.h \
dynare_model.h \
forw_subst_builder.h \
planner_builder.cpp \
dynare3.h \
dynare_exception.h \
dynare_params.cpp \
planner_builder.h \
dynare_atoms.cpp \
dynare_model.cpp \
dynare_params.h \
forw_subst_builder.cpp \
nlsolve.cpp \
nlsolve.h \
$(GENERATED_FILES)
dynare___CPPFLAGS = -I../sylv/cc -I../tl/cc -I../kord -I../integ/cc -I.. -I$(top_srcdir)/mex/sources -DDYNVERSION=\"$(PACKAGE_VERSION)\" $(BOOST_CPPFLAGS) $(CPPFLAGS_MATIO)
dynare___LDFLAGS = $(LDFLAGS_MATIO) $(BOOST_LDFLAGS)
dynare___LDADD = ../kord/libkord.a ../integ/cc/libinteg.a ../tl/cc/libtl.a ../parser/cc/libparser.a ../utils/cc/libutils.a ../sylv/cc/libsylv.a $(LIBADD_MATIO) $(noinst_LIBRARIES) $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(PTHREAD_LIBS)
dynare___CXXFLAGS = $(PTHREAD_CFLAGS)
BUILT_SOURCES = $(GENERATED_FILES)
EXTRA_DIST = dynglob.lex dynglob.y
dynglob_tab.cc dynglob_tab.hh: dynglob.y
$(YACC) -d -odynglob_tab.cc dynglob.y
dynglob_ll.cc: dynglob.lex
$(LEX) -i -odynglob_ll.cc dynglob.lex
dynare++/src/dynare3.cpp deleted 100644 → 0
View file @ 6ca649c4
#include "dynare3.h"
#include "dynare_exception.h"
#include "planner_builder.h"
#include "forw_subst_builder.h"
#include "utils/cc/memory_file.h"
#include "utils/cc/exception.h"
#include "parser/cc/parser_exception.h"
#include "parser/cc/atom_substitutions.h"
#include "../tl/cc/tl_exception.h"
#include "../kord/kord_exception.h"
#ifndef DYNVERSION
#define DYNVERSION "unknown"
#endif
/**************************************************************************************/
/* DynareNameList class */
/**************************************************************************************/
vector<int> DynareNameList::selectIndices(const vector<const char*>& ns) const
{
vector<int> res;
for (unsigned int i = 0; i < ns.size(); i++) {
int j = 0;
while (j < getNum() && strcmp(getName(j), ns[i]) != 0)
j++;
if (j == getNum())
throw DynareException(__FILE__, __LINE__,
string("Couldn't find name for ") + ns[i] +
" in DynareNameList::selectIndices");
res.push_back(j);
}
return res;
}
/**************************************************************************************/
/* Dynare class */
/**************************************************************************************/
Dynare::Dynare(const char* modname, int ord, double sstol, Journal& jr)
: journal(jr), model(NULL), ysteady(NULL), md(1), dnl(NULL), denl(NULL), dsnl(NULL),
fe(NULL), fde(NULL), ss_tol(sstol)
{
// make memory file
ogu::MemoryFile mf(modname);
if (mf.exists()) {
try {
model = new ogdyn::DynareParser(mf.base(), mf.length(), ord);
} catch (const ogp::ParserException& pe) {
int line;
int col;
mf.line_and_col(pe.offset(), line, col);
throw DynareException(pe.message(), modname, line, col);
}
ysteady = new Vector(model->getAtoms().ny());
dnl = new DynareNameList(*this);
denl = new DynareExogNameList(*this);
dsnl = new DynareStateNameList(*this, *dnl, *denl);
fe = new ogp::FormulaEvaluator(model->getParser());
fde = new ogp::FormulaDerEvaluator(model->getParser());
writeModelInfo(journal);
} else {
throw DynareException(__FILE__, __LINE__, string("Could not open model file ")+modname);
}
}
Dynare::Dynare(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,
double sstol, Journal& jr)
: journal(jr), model(NULL), ysteady(NULL), md(1), dnl(NULL), denl(NULL), dsnl(NULL),
fe(NULL), fde(NULL), ss_tol(sstol)
{
try {
model = new ogdyn::DynareSPModel(endo, num_endo, exo, num_exo, par, num_par,
equations, len, ord);
} catch (const ogp::ParserException& pe) {
throw DynareException(pe.message(), pe.offset());
}
ysteady = new Vector(model->getAtoms().ny());
dnl = new DynareNameList(*this);
denl = new DynareExogNameList(*this);
dsnl = new DynareStateNameList(*this, *dnl, *denl);
fe = new ogp::FormulaEvaluator(model->getParser());
fde = new ogp::FormulaDerEvaluator(model->getParser());
writeModelInfo(journal);
}
Dynare::Dynare(const Dynare& dynare)
: journal(dynare.journal), model(NULL),
ysteady(NULL), md(dynare.md),
dnl(NULL), denl(NULL), dsnl(NULL), fe(NULL), fde(NULL),
ss_tol(dynare.ss_tol)
{
model = dynare.model->clone();
ysteady = new Vector(*(dynare.ysteady));
dnl = new DynareNameList(*this);
denl = new DynareExogNameList(*this);
dsnl = new DynareStateNameList(*this, *dnl, *denl);
fe = new ogp::FormulaEvaluator(model->getParser());
fde = new ogp::FormulaDerEvaluator(model->getParser());
}
Dynare::~Dynare()
{
if (model)
delete model;
if (ysteady)
delete ysteady;
if (dnl)
delete dnl;
if (dsnl)
delete dsnl;
if (denl)
delete denl;
if (fe)
delete fe;
if (fde)
delete fde;
}
void Dynare::writeMat(mat_t* fd, const char* prefix) const
{
char tmp[100];
sprintf(tmp, "%s_vars", prefix);
getAllEndoNames().writeMat(fd, tmp);
getAllEndoNames().writeMatIndices(fd, prefix);
sprintf(tmp, "%s_state_vars", prefix);
getStateNames().writeMat(fd, tmp);
sprintf(tmp, "%s_shocks", prefix);
getExogNames().writeMat(fd, tmp);
getExogNames().writeMatIndices(fd, prefix);
sprintf(tmp, "%s_vcov_exo", prefix);
model->getVcov().writeMat(fd, tmp);
TwoDMatrix aux(1,1);
sprintf(tmp, "%s_nstat", prefix);
aux.get(0,0) = nstat();
aux.writeMat(fd, tmp);
sprintf(tmp, "%s_npred", prefix);
aux.get(0,0) = npred();
aux.writeMat(fd, tmp);
sprintf(tmp, "%s_nboth", prefix);
aux.get(0,0) = nboth();
aux.writeMat(fd, tmp);
sprintf(tmp, "%s_nforw", prefix);
aux.get(0,0) = nforw();
aux.writeMat(fd, tmp);
}
void Dynare::writeDump(const std::string& basename) const
{
std::string fname(basename);
fname += ".dump";
std::ofstream out(fname.c_str());
model->dump_model(out);
out.close();
}
void Dynare::solveDeterministicSteady(Vector& steady)
{
JournalRecordPair pa(journal);
pa << "Non-linear solver for deterministic steady state" << endrec;
steady = (const Vector&) model->getInit();
DynareVectorFunction dvf(*this);
DynareJacobian dj(*this);
ogu::NLSolver nls(dvf, dj, 500, ss_tol, journal);
int iter;
if (! nls.solve(steady, iter))
throw DynareException(__FILE__, __LINE__,
"Could not obtain convergence in non-linear solver");
}
// evaluate system at given y_t=y_{t+1}=y_{t-1}, and given shocks x_t
void Dynare::evaluateSystem(Vector& out, const Vector& yy, const Vector& xx)
{
ConstVector yym(yy, nstat(), nys());
ConstVector yyp(yy, nstat()+npred(), nyss());
evaluateSystem(out, yym, yy, yyp, xx);
}
// evaluate system at given y^*_{t-1}, y_t, y^{**}_{t+1} and at
// exogenous x_t, all three vectors yym, yy, and yyp have the
// respective lengths of y^*_{t-1}, y_t, y^{**}_{t+1}
void Dynare::evaluateSystem(Vector& out, const Vector& yym, const Vector& yy,
const Vector& yyp, const Vector& xx)
{
ogdyn::DynareAtomValues dav(model->getAtoms(), model->getParams(), yym, yy, yyp, xx);
DynareEvalLoader del(model->getAtoms(), out);
fe->eval(dav, del);
}
void Dynare::calcDerivatives(const Vector& yy, const Vector& xx)
{
ConstVector yym(yy, nstat(), nys());
ConstVector yyp(yy, nstat()+npred(), nyss());
ogdyn::DynareAtomValues dav(model->getAtoms(), model->getParams(), yym, yy, yyp, xx);
DynareDerEvalLoader ddel(model->getAtoms(), md, model->getOrder());
for (int iord = 1; iord <= model->getOrder(); iord++)
fde->eval(dav, ddel, iord);
}
void Dynare::calcDerivativesAtSteady()
{
Vector xx(nexog());
xx.zeros();
calcDerivatives(*ysteady, xx);
}
void Dynare::writeModelInfo(Journal& jr) const
{
// write info on variables
{
JournalRecordPair rp(journal);
rp << "Information on variables" << endrec;
JournalRecord rec1(journal);
rec1 << "Number of endogenous: " << ny() << endrec;
JournalRecord rec2(journal);
rec2 << "Number of exogenous: " << nexog() << endrec;
JournalRecord rec3(journal);
rec3 << "Number of static: " << nstat() << endrec;
JournalRecord rec4(journal);
rec4 << "Number of predetermined: " << npred()+nboth() << endrec;
JournalRecord rec5(journal);
rec5 << "Number of forward looking: " << nforw()+nboth() << endrec;
JournalRecord rec6(journal);
rec6 << "Number of both: " << nboth() << endrec;
}
// write info on planner variables
const ogdyn::PlannerInfo* pinfo = model->get_planner_info();
if (pinfo) {
JournalRecordPair rp(journal);
rp << "Information on planner variables" << endrec;
JournalRecord rec1(journal);
rec1 << "Number of Lagrange multipliers: " << pinfo->num_lagrange_mults << endrec;
JournalRecord rec2(journal);
rec2 << "Number of auxiliary variables: " << pinfo->num_aux_variables << endrec;
JournalRecord rec3(journal);
rec3 << "Number of new terms in the tree: " << pinfo->num_new_terms << endrec;
}
// write info on forward substitutions
const ogdyn::ForwSubstInfo* finfo = model->get_forw_subst_info();
if (finfo) {
JournalRecordPair rp(journal);
rp << "Information on forward substitutions" << endrec;
JournalRecord rec1(journal);
rec1 << "Number of affected equations: " << finfo->num_affected_equations << endrec;
JournalRecord rec2(journal);
rec2 << "Number of substituted terms: " << finfo->num_subst_terms << endrec;
JournalRecord rec3(journal);
rec3 << "Number of auxiliary variables: " << finfo->num_aux_variables << endrec;
JournalRecord rec4(journal);
rec4 << "Number of new terms in the tree: " << finfo->num_new_terms << endrec;
}
// write info on substitutions
const ogp::SubstInfo* sinfo = model->get_subst_info();
if (sinfo) {
JournalRecordPair rp(journal);
rp << "Information on substitutions" << endrec;
JournalRecord rec1(journal);
rec1 << "Number of substitutions: " << sinfo->num_substs << endrec;
}
}
DynareNameList::DynareNameList(const Dynare& dynare)
{
for (int i = 0; i < dynare.ny(); i++) {
int j = dynare.model->getAtoms().y2outer_endo()[i];
const char* name = dynare.model->getAtoms().get_endovars()[j];
names.push_back(name);
}
}
DynareStateNameList::DynareStateNameList(const Dynare& dynare, const DynareNameList& dnl,
const DynareExogNameList& denl)
{
for (int i = 0; i < dynare.nys(); i++)
names.push_back(dnl.getName(i+dynare.nstat()));
for (int i = 0; i < dynare.nexog(); i++)
names.push_back(denl.getName(i));
}
DynareExogNameList::DynareExogNameList(const Dynare& dynare)
{
for (int i = 0; i < dynare.nexog(); i++) {
int j = dynare.model->getAtoms().y2outer_exo()[i];
const char* name = dynare.model->getAtoms().get_exovars()[j];
names.push_back(name);
}
}
DynareEvalLoader::DynareEvalLoader(const ogp::FineAtoms& a, Vector& out)
: Vector(out)
{
if (a.ny() != out.length())
throw DynareException(__FILE__, __LINE__, "Wrong length of out vector in DynareEvalLoader constructor");
}
/** This clears the container of model derivatives and initializes it
* inserting empty sparse tensors up to the given order. */
DynareDerEvalLoader::DynareDerEvalLoader(const ogp::FineAtoms& a,
TensorContainer<FSSparseTensor>& mod_ders,
int order)
: atoms(a), md(mod_ders)
{
md.clear();
for (int iord = 1; iord <= order; iord++) {
FSSparseTensor* t = new FSSparseTensor(iord, atoms.ny()+atoms.nys()+atoms.nyss()+atoms.nexo(), atoms.ny());
md.insert(t);
}
}
void DynareDerEvalLoader::load(int i, int iord, const int* vars, double res)
{
FSSparseTensor* t = md.get(Symmetry(iord));
IntSequence s(iord, 0);
for (int j = 0; j < iord; j++)
s[j] = atoms.get_pos_of_all(vars[j]);
t->insert(s, i, res);
}
DynareJacobian::DynareJacobian(Dynare& dyn)
: Jacobian(dyn.ny()), d(dyn)
{
zeros();
}
void DynareJacobian::eval(const Vector& yy)
{
ogdyn::DynareSteadyAtomValues
dav(d.getModel().getAtoms(), d.getModel().getParams(), yy);
zeros();
d.fde->eval(dav, *this, 1);
}
void DynareJacobian::load(int i, int iord, const int* vars, double res)
{
if (iord != 1)
throw DynareException(__FILE__, __LINE__,
"Derivative order different from order=1 in DynareJacobian::load");
int t = vars[0];
int j = d.getModel().getAtoms().get_pos_of_all(t);
if (j < d.nyss())
get(i, j+d.nstat()+d.npred()) += res;
else if (j < d.nyss()+d.ny())
get(i, j-d.nyss()) += res;
else if (j < d.nyss()+d.ny()+d.nys())
get(i, j-d.nyss()-d.ny()+d.nstat()) += res;
}
void DynareVectorFunction::eval(const ConstVector& in, Vector& out)
{
check_for_eval(in, out);
Vector xx(d.nexog());
xx.zeros();
d.evaluateSystem(out, in, xx);
}
dynare++/src/dynare3.h deleted 100644 → 0
View file @ 6ca649c4
// $Id: dynare3.h 1764 2008-03-31 14:30:55Z kamenik $
// Copyright 2005, Ondra Kamenik
#ifndef DYNARE3_H
#define DYNARE3_H
#include "../tl/cc/t_container.h"
#include "../tl/cc/sparse_tensor.h"
#include "../kord/decision_rule.h"
#include "../kord/dynamic_model.h"
#include "dynare_model.h"
#include "nlsolve.h"
#include <vector>
#include <matio.h>
class Dynare;
class DynareNameList : public NameList {
vector<const char*> names;
public:
DynareNameList(const Dynare& dynare);
int getNum() const
{return (int)names.size();}
const char* getName(int i) const
{return names[i];}
/** This for each string of the input vector calculates its index
* in the names. And returns the resulting vector of indices. If
* the name cannot be found, then an exception is raised. */
vector<int> selectIndices(const vector<const char*>& ns) const;
};
class DynareExogNameList : public NameList {
vector<const char*> names;
public:
DynareExogNameList(const Dynare& dynare);
int getNum() const
{return (int)names.size();}
const char* getName(int i) const
{return names[i];}
};
class DynareStateNameList : public NameList {
vector<const char*> names;
public:
DynareStateNameList(const Dynare& dynare, const DynareNameList& dnl,
const DynareExogNameList& denl);
int getNum() const
{return (int)names.size();}
const char* getName(int i) const
{return names[i];}
};
// The following only implements DynamicModel with help of ogdyn::DynareModel
class DynareJacobian;
class Dynare : public DynamicModel {
friend class DynareNameList;
friend class DynareExogNameList;
friend class DynareStateNameList;
friend class DynareJacobian;
Journal& journal;
ogdyn::DynareModel* model;
Vector* ysteady;
TensorContainer<FSSparseTensor> md;
DynareNameList* dnl;
DynareExogNameList* denl;
DynareStateNameList* dsnl;
ogp::FormulaEvaluator* fe;
ogp::FormulaDerEvaluator* fde;
const double ss_tol;
public:
/** Parses the given model file and uses the given order to
* override order from the model file (if it is != -1). */
Dynare(const char* modname, int ord, double sstol, Journal& jr);
/** Parses the given equations with explicitly given names. */
Dynare(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,
double sstol, Journal& jr);
/** Makes a deep copy of the object. */
Dynare(const Dynare& dyn);
DynamicModel* clone() const
{return new Dynare(*this);}
virtual ~Dynare();
int nstat() const
{return model->getAtoms().nstat();}
int nboth() const
{return model->getAtoms().nboth();}
int npred() const
{return model->getAtoms().npred();}
int nforw() const
{return model->getAtoms().nforw();}
int nexog() const
{return model->getAtoms().nexo();}
int nys() const
{return model->getAtoms().nys();}
int nyss() const
{return model->getAtoms().nyss();}
int ny() const
{return model->getAtoms().ny();}
int order() const
{return model->getOrder();}
const NameList& getAllEndoNames() const
{return *dnl;}
const NameList& getStateNames() const
{return *dsnl;}
const NameList& getExogNames() const
{return *denl;}
TwoDMatrix& getVcov()
{return model->getVcov();}
const TwoDMatrix& getVcov() const
{return model->getVcov();}
Vector& getParams()
{return model->getParams();}
const Vector& getParams() const
{return model->getParams();}
void setInitOuter(const Vector& x)
{model->setInitOuter(x);}
const TensorContainer<FSSparseTensor>& getModelDerivatives() const
{return md;}
const Vector& getSteady() const
{return *ysteady;}
Vector& getSteady()
{return *ysteady;}
const ogdyn::DynareModel& getModel() const
{return *model;}
// here is true public interface
void solveDeterministicSteady(Vector& steady);
void solveDeterministicSteady()
{solveDeterministicSteady(*ysteady);}
void evaluateSystem(Vector& out, const Vector& yy, const Vector& xx);
void evaluateSystem(Vector& out, const Vector& yym, const Vector& yy,
const Vector& yyp, const Vector& xx);
void calcDerivatives(const Vector& yy, const Vector& xx);
void calcDerivativesAtSteady();
void writeMat(mat_t *fd, const char* prefix) const;
void writeDump(const std::string& basename) const;
private:
void writeModelInfo(Journal& jr) const;
};
class DynareEvalLoader : public ogp::FormulaEvalLoader, public Vector {
public:
DynareEvalLoader(const ogp::FineAtoms& a, Vector& out);
void load(int i, double res)
{operator[](i) = res;}
};
class DynareDerEvalLoader : public ogp::FormulaDerEvalLoader {
protected:
const ogp::FineAtoms& atoms;
TensorContainer<FSSparseTensor>& md;
public:
DynareDerEvalLoader(const ogp::FineAtoms& a, TensorContainer<FSSparseTensor>& mod_ders,
int order);
void load(int i, int iord, const int* vars, double res);
};
class DynareJacobian : public ogu::Jacobian, public ogp::FormulaDerEvalLoader {
protected:
Dynare& d;
public:
DynareJacobian(Dynare& dyn);
virtual ~DynareJacobian() {}
void load(int i, int iord, const int* vars, double res);
void eval(const Vector& in);
};
class DynareVectorFunction : public ogu::VectorFunction {
protected:
Dynare& d;
public:
DynareVectorFunction(Dynare& dyn)
: d(dyn) {}
virtual ~DynareVectorFunction() {}
int inDim() const
{return d.ny();}
int outDim() const
{return d.ny();}
void eval(const ConstVector& in, Vector& out);
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/src/dynare_atoms.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: dynare_atoms.cpp 1765 2008-03-31 14:32:08Z kamenik $
#include "parser/cc/parser_exception.h"
#include "utils/cc/exception.h"
#include "dynare_atoms.h"
#include <string>
#include <cmath>
using namespace ogdyn;
using std::string;
void DynareStaticAtoms::register_name(const char* name)
{
if (varnames.query(name))
throw ogp::ParserException(string("The name ")+name+" is not unique.", 0);
StaticAtoms::register_name(name);
}
int DynareStaticAtoms::check_variable(const char* name) const
{
if (0 == varnames.query(name))
throw ogp::ParserException(std::string("Unknown name <")+name+">", 0);
Tvarmap::const_iterator it = vars.find(name);
if (it == vars.end())
return -1;
else
return (*it).second;
}
DynareDynamicAtoms::DynareDynamicAtoms(const DynareDynamicAtoms& dda)
: SAtoms(dda)
{
// fill atom_type
for (Tatypemap::const_iterator it = dda.atom_type.begin();
it != dda.atom_type.end(); ++it)
atom_type.insert(Tatypemap::value_type(varnames.query((*it).first), (*it).second));
}
void DynareDynamicAtoms::parse_variable(const char* in, std::string& out, int& ll) const
{
ll = 0;
std::string str = in;
int left = str.find_first_of("({");
if (left != -1) {
out = str.substr(0, left);
left++;
int right = str.find_first_of(")}", left);
if ((int)string::npos == right)
throw ogp::ParserException(
string("Syntax error when parsing Dynare atom <")+in+">.", 0);
std::string tmp(str, left, right-left);
sscanf(tmp.c_str(), "%d", &ll);
} else {
out = in;
}
}
void DynareDynamicAtoms::register_uniq_endo(const char* name)
{
FineAtoms::register_uniq_endo(name);
atom_type.insert(Tatypemap::value_type(varnames.query(name), endovar));
}
void DynareDynamicAtoms::register_uniq_exo(const char* name)
{
FineAtoms::register_uniq_exo(name);
atom_type.insert(Tatypemap::value_type(varnames.query(name), exovar));
}
void DynareDynamicAtoms::register_uniq_param(const char* name)
{
FineAtoms::register_uniq_param(name);
atom_type.insert(Tatypemap::value_type(varnames.query(name), param));
}
bool DynareDynamicAtoms::is_type(const char* name, atype tp) const
{
Tatypemap::const_iterator it = atom_type.find(name);
if (it != atom_type.end() && (*it).second == tp)
return true;
else
return false;
}
void DynareDynamicAtoms::print() const
{
SAtoms::print();
printf("Name types:\n");
for (Tatypemap::const_iterator it = atom_type.begin();
it != atom_type.end(); ++it)
printf("name=%s type=%s\n", (*it).first,
((*it).second == endovar) ? "endovar" : (((*it).second == exovar)? "exovar" : "param"));
}
std::string DynareDynamicAtoms::convert(int t) const
{
if (t < ogp::OperationTree::num_constants) {
throw ogu::Exception(__FILE__,__LINE__,
"Tree index is a built-in constant in DynareDynamicAtoms::convert");
return std::string();
}
if (is_constant(t)) {
double v = get_constant_value(t);
char buf[100];
sprintf(buf, "%20.16g", v);
const char* s = buf;
while (*s == ' ')
++s;
return std::string(s);
}
const char* s = name(t);
if (is_type(s, endovar)) {
int ll = lead(t);
char buf[100];
if (ll)
sprintf(buf, "%s(%d)", s, ll);
else
sprintf(buf, "%s", s);
return std::string(buf);
}
return std::string(s);
}
void DynareAtomValues::setValues(ogp::EvalTree& et) const
{
// set constants
atoms.setValues(et);
// set parameteres
for (unsigned int i = 0; i < atoms.get_params().size(); i++) {
try {
const ogp::DynamicAtoms::Tlagmap& lmap = atoms.lagmap(atoms.get_params()[i]);
for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
it != lmap.end(); ++it) {
int t = (*it).second;
et.set_nulary(t, paramvals[i]);
}
} catch (const ogu::Exception& e) {
// ignore non-referenced parameters; there is no
// lagmap for them
}
}
// set endogenous
for (unsigned int outer_i = 0; outer_i < atoms.get_endovars().size(); outer_i++) {
try {
const ogp::DynamicAtoms::Tlagmap& lmap = atoms.lagmap(atoms.get_endovars()[outer_i]);
for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
it != lmap.end(); ++it) {
int ll = (*it).first;
int t = (*it).second;
int i = atoms.outer2y_endo()[outer_i];
if (ll == -1) {
et.set_nulary(t, yym[i-atoms.nstat()]);
}
else if (ll == 0)
et.set_nulary(t, yy[i]);
else
et.set_nulary(t, yyp[i-atoms.nstat()-atoms.npred()]);
}
} catch (const ogu::Exception& e) {
// ignore non-referenced endogenous variables; there is no
// lagmap for them
}
}
// set exogenous
for (unsigned int outer_i = 0; outer_i < atoms.get_exovars().size(); outer_i++) {
try {
const ogp::DynamicAtoms::Tlagmap& lmap = atoms.lagmap(atoms.get_exovars()[outer_i]);
for (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
it != lmap.end(); ++it) {
int ll = (*it).first;
if (ll == 0) { // this is always true because of checks
int t = (*it).second;
int i = atoms.outer2y_exo()[outer_i];
et.set_nulary(t, xx[i]);
}
}
} catch (const ogu::Exception& e) {
// ignore non-referenced variables
}
}
}
void DynareStaticSteadyAtomValues::setValues(ogp::EvalTree& et) const
{
// set constants
atoms_static.setValues(et);
// set parameters
for (unsigned int i = 0; i < atoms_static.get_params().size(); i++) {
const char* name = atoms_static.get_params()[i];
int t = atoms_static.index(name);
if (t != -1) {
int idyn = atoms.name2outer_param(name);
et.set_nulary(t, paramvals[idyn]);
}
}
// set endogenous
for (unsigned int i = 0; i < atoms_static.get_endovars().size(); i++) {
const char* name = atoms_static.get_endovars()[i];
int t = atoms_static.index(name);
if (t != -1) {
int idyn = atoms.outer2y_endo()[atoms.name2outer_endo(name)];
et.set_nulary(t, yy[idyn]);
}
}
// set exogenous
for (unsigned int i = 0; i < atoms_static.get_exovars().size(); i++) {
const char* name = atoms_static.get_exovars()[i];
int t = atoms_static.index(name);
if (t != -1)
et.set_nulary(t, 0.0);
}
}
DynareSteadySubstitutions::DynareSteadySubstitutions(const ogp::FineAtoms& a,
const ogp::OperationTree& tree,
const Tsubstmap& subst,
const Vector& pvals, Vector& yy)
: atoms(a), y(yy)
{
// fill the vector of left and right hand sides
for (Tsubstmap::const_iterator it = subst.begin();
it != subst.end(); ++it) {
left_hand_sides.push_back((*it).first);
right_hand_sides.push_back((*it).second);
}
// evaluate right hand sides
DynareSteadyAtomValues dsav(atoms, pvals, y);
ogp::FormulaCustomEvaluator fe(tree, right_hand_sides);
fe.eval(dsav, *this);
}
void DynareSteadySubstitutions::load(int i, double res)
{
const char* name = left_hand_sides[i];
int iouter = atoms.name2outer_endo(name);
int iy = atoms.outer2y_endo()[iouter];
if (! std::isfinite(y[iy]))
y[iy] = res;
}
DynareStaticSteadySubstitutions::
DynareStaticSteadySubstitutions(const ogp::FineAtoms& a, const ogp::StaticFineAtoms& sa,
const ogp::OperationTree& tree,
const Tsubstmap& subst,
const Vector& pvals, Vector& yy)
: atoms(a), atoms_static(sa), y(yy)
{
// fill the vector of left and right hand sides
for (Tsubstmap::const_iterator it = subst.begin();
it != subst.end(); ++it) {
left_hand_sides.push_back((*it).first);
right_hand_sides.push_back((*it).second);
}
// evaluate right hand sides
DynareStaticSteadyAtomValues dsav(atoms, atoms_static, pvals, y);
ogp::FormulaCustomEvaluator fe(tree, right_hand_sides);
fe.eval(dsav, *this);
}
void DynareStaticSteadySubstitutions::load(int i, double res)
{
const char* name = left_hand_sides[i];
int iouter = atoms.name2outer_endo(name);
int iy = atoms.outer2y_endo()[iouter];
if (! std::isfinite(y[iy]))
y[iy] = res;
}
dynare++/src/dynare_atoms.h deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: dynare_atoms.h 1765 2008-03-31 14:32:08Z kamenik $
#ifndef OGDYN_DYNARE_ATOMS_H
#define OGDYN_DYNARE_ATOMS_H
#include "sylv/cc/Vector.h"
#include "parser/cc/static_atoms.h"
#include "parser/cc/static_fine_atoms.h"
#include "parser/cc/atom_substitutions.h"
#include "parser/cc/tree.h"
#include <map>
#include <vector>
namespace ogdyn {
using std::map;
using std::vector;
/** A definition of a type mapping a string to an integer. Used as
* a substitution map, saying what names are substituted for what
* expressions represented by tree indices. */
typedef map<const char*, int, ogp::ltstr> Tsubstmap;
class DynareStaticAtoms : public ogp::StaticAtoms {
public:
DynareStaticAtoms()
: StaticAtoms() {}
DynareStaticAtoms(const DynareStaticAtoms& a)
: StaticAtoms(a) {}
virtual ~DynareStaticAtoms() {}
/** This registers a unique varname identifier. It throws an
* exception if the variable name is duplicate. It checks the
* uniqueness and then it calls StaticAtoms::register_name. */
void register_name(const char* name);
protected:
/** This returns a tree index of the given variable, and if
* the variable has not been registered, it throws an
* exception. */
int check_variable(const char* name) const;
};
class DynareDynamicAtoms : public ogp::SAtoms, public ogp::NularyStringConvertor {
public:
enum atype {endovar, exovar, param};
protected:
typedef map<const char*, atype, ogp::ltstr> Tatypemap;
/** The map assigining a type to each name. */
Tatypemap atom_type;
public:
DynareDynamicAtoms()
: ogp::SAtoms() {}
DynareDynamicAtoms(const DynareDynamicAtoms& dda);
virtual ~DynareDynamicAtoms() {}
/** This parses a variable of the forms: varname(+3),
* varname(3), varname, varname(-3), varname(0), varname(+0),
* varname(-0). */
virtual void parse_variable(const char* in, std::string& out, int& ll) const;
/** Registers unique name of endogenous variable. */
void register_uniq_endo(const char* name);
/** Registers unique name of exogenous variable. */
void register_uniq_exo(const char* name);
/** Registers unique name of parameter. */
void register_uniq_param(const char* name);
/** Return true if the name is a given type. */
bool is_type(const char* name, atype tp) const;
/** Debug print. */
void print() const;
/** Implement NularyStringConvertor::convert. */
std::string convert(int t) const;
};
/** This class represents the atom values for dynare, where
* exogenous variables can occur only at time t, and endogenous at
* times t-1, t, and t+1. */
class DynareAtomValues : public ogp::AtomValues {
protected:
/** Reference to the atoms (we suppose that they are only at
* t-1,t,t+1. */
const ogp::FineAtoms& atoms;
/** De facto reference to the values of parameters. */
const ConstVector paramvals;
/** De facto reference to the values of endogenous at time t-1. Only
* predetermined and both part. */
const ConstVector yym;
/** De facto reference to the values of endogenous at time t. Ordering
* given by the atoms. */
const ConstVector yy;
/** De facto reference to the values of endogenous at time t+1. Only
* both and forward looking part. */
const ConstVector yyp;
/** De facto reference to the values of exogenous at time t. */
const ConstVector xx;
public:
DynareAtomValues(const ogp::FineAtoms& a, const Vector& pvals, const Vector& ym,
const Vector& y, const Vector& yp, const Vector& x)
: atoms(a), paramvals(pvals), yym(ym), yy(y), yyp(yp), xx(x) {}
DynareAtomValues(const ogp::FineAtoms& a, const Vector& pvals, const ConstVector& ym,
const Vector& y, const ConstVector& yp, const Vector& x)
: atoms(a), paramvals(pvals), yym(ym), yy(y), yyp(yp), xx(x) {}
void setValues(ogp::EvalTree& et) const;
};
/** This class represents the atom values at the steady state. It
* makes only appropriate subvector yym and yyp of the y vector,
* makes a vector of zero exogenous variables and uses
* DynareAtomValues with more general interface. */
class DynareSteadyAtomValues : public ogp::AtomValues {
protected:
/** Subvector of yy. */
const ConstVector yym;
/** Subvector of yy. */
const ConstVector yyp;
/** Vector of zeros for exogenous variables. */
Vector xx;
/** Atom values using this yym, yyp and xx. */
DynareAtomValues av;
public:
DynareSteadyAtomValues(const ogp::FineAtoms& a, const Vector& pvals, const Vector& y)
: yym(y, a.nstat(), a.nys()),
yyp(y, a.nstat()+a.npred(), a.nyss()),
xx(a.nexo()),
av(a, pvals, yym, y, yyp, xx)
{xx.zeros();}
void setValues(ogp::EvalTree& et) const
{av.setValues(et);}
};
class DynareStaticSteadyAtomValues : public ogp::AtomValues {
protected:
/** Reference to static atoms over which the tree, where the
* values go, is defined. */
const ogp::StaticFineAtoms& atoms_static;
/** Reference to dynamic atoms for which the class gets input
* data. */
const ogp::FineAtoms& atoms;
/** De facto reference to input data, this is a vector of
* endogenous variables in internal ordering of the dynamic
* atoms. */
ConstVector yy;
/** De facto reference to input parameters corresponding to
* ordering defined by the dynamic atoms. */
ConstVector paramvals;
public:
/** Construct the object. */
DynareStaticSteadyAtomValues(const ogp::FineAtoms& a, const ogp::StaticFineAtoms& sa,
const Vector& pvals, const Vector& yyy)
: atoms_static(sa),
atoms(a),
yy(yyy),
paramvals(pvals) {}
/** Set the values to the tree defined over the static atoms. */
void setValues(ogp::EvalTree& et) const;
};
/** This class takes a vector of endogenous variables and a
* substitution map. It supposes that variables at the right hand
* sides of the substitutions are set in the endogenous vector. It
* evaluates the substitutions and if the variables corresponding
* to left hand sides are not set in the endogenous vector it sets
* them to calculated values. If a variable is already set, it
* does not override its value. It has no methods, everything is
* done in the constructor. */
class DynareSteadySubstitutions : public ogp::FormulaEvalLoader {
protected:
const ogp::FineAtoms& atoms;
public:
DynareSteadySubstitutions(const ogp::FineAtoms& a, const ogp::OperationTree& tree,
const Tsubstmap& subst,
const Vector& pvals, Vector& yy);
void load(int i, double res);
protected:
Vector& y;
vector<const char*> left_hand_sides;
vector<int> right_hand_sides;
};
/** This class is a static version of DynareSteadySustitutions. It
* works for static atoms and static tree and substitution map
* over the static tree. It also needs dynamic version of the
* atoms, since it defines ordering of the vectors pvals, and
* yy. */
class DynareStaticSteadySubstitutions : public ogp::FormulaEvalLoader {
protected:
const ogp::FineAtoms& atoms;
const ogp::StaticFineAtoms& atoms_static;
public:
DynareStaticSteadySubstitutions(const ogp::FineAtoms& a,
const ogp::StaticFineAtoms& sa,
const ogp::OperationTree& tree,
const Tsubstmap& subst,
const Vector& pvals, Vector& yy);
void load(int i, double res);
protected:
Vector& y;
vector<const char*> left_hand_sides;
vector<int> right_hand_sides;
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/src/dynare_exception.h deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006, Ondra Kamenik
// $Id: dynare_exception.h 853 2006-08-01 08:42:42Z kamenik $
#ifndef DYNARE_EXCEPTION_H
#define DYNARE_EXCEPTION_H
#include <string>
class DynareException {
char* mes;
public:
DynareException(const char* m, const char* fname, int line, int col)
{
mes = new char[strlen(m) + strlen(fname) + 100];
sprintf(mes, "Parse error at %s, line %d, column %d: %s", fname, line, col, m);
}
DynareException(const char* fname, int line, const std::string& m)
{
mes = new char[m.size() + strlen(fname) + 50];
sprintf(mes, "%s:%d: %s", fname, line, m.c_str());
}
DynareException(const char* m, int offset)
{
mes = new char[strlen(m) + 100];
sprintf(mes, "Parse error in provided string at offset %d: %s", offset, m);
}
DynareException(const DynareException& e)
: mes(new char[strlen(e.mes)+1])
{strcpy(mes, e.mes);}
virtual ~DynareException()
{delete [] mes;}
const char* message() const
{return mes;}
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/src/dynare_model.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2006-2011, Ondra Kamenik
#include "parser/cc/parser_exception.h"
#include "parser/cc/location.h"
#include "utils/cc/exception.h"
#include "dynare_model.h"
#include "dynare_exception.h"
#include "planner_builder.h"
#include "forw_subst_builder.h"
#include <cstdlib>
#include <string>
#include <cmath>
#include <climits>
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), order(-1),
param_vals(0), init_vals(0), vcov_mat(0),
t_plobjective(-1), t_pldiscount(-1),
pbuilder(NULL), fbuilder(NULL),
atom_substs(NULL), old_atoms(NULL)
{}
DynareModel::DynareModel(const DynareModel& dm)
: atoms(dm.atoms), eqs(dm.eqs, atoms), order(dm.order),
param_vals(0), init_vals(0), vcov_mat(0),
t_plobjective(dm.t_plobjective),
t_pldiscount(dm.t_pldiscount),
pbuilder(NULL), fbuilder(NULL),
atom_substs(NULL), old_atoms(NULL)
{
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 NULL;
}
const ForwSubstInfo* DynareModel::get_forw_subst_info() const
{
if (fbuilder)
return &(fbuilder->get_info());
return NULL;
}
const ogp::SubstInfo* DynareModel::get_subst_info() const
{
if (atom_substs)
return &(atom_substs->get_info());
return NULL;
}
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 (int i = 0; i < (int)atoms.get_endovars().size(); i++)
os << " " << atoms.get_endovars()[i];
os << ";\n\n";
// exogenous variables
os << "varexo";
for (int i = 0; i < (int)atoms.get_exovars().size(); i++)
os << " " << atoms.get_exovars()[i];
os << ";\n\n";
// parameters
os << "parameters";
for (int i = 0; i < (int)atoms.get_params().size(); i++)
os << " " << atoms.get_params()[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 (unordered_set<int>::const_iterator it = nuls.begin();
it != nuls.end(); ++it)
if (! atoms.is_constant(*it) && ! atoms.is_named_atom(*it))
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 (unordered_set<int>::const_iterator it = a_set.begin();
it != a_set.end(); ++it) {
int t = *it;
// 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 (unordered_set<int>::const_iterator ni = nul_terms.begin();
ni != nul_terms.end(); ++ni) {
if (!atoms.is_constant(*ni) &&
(atoms.is_type(atoms.name(*ni), DynareDynamicAtoms::endovar) ||
atoms.is_type(atoms.name(*ni), DynareDynamicAtoms::exovar))) {
int ll = atoms.lead(*ni);
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 (unordered_set<int>::const_iterator ni = nul_terms.begin();
ni != nul_terms.end(); ++ni)
if (! atoms.is_constant(*ni) &&
! atoms.is_type(atoms.name(*ni), 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()
{
}
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)
{
char* 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)
{
char* 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)
{
char* 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 (unsigned int j = 0; j < eam.size(); j++) {
int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, eam[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 (unsigned int ip = 0; ip < model.getAtoms().get_params().size(); ip++) {
const char* parname = model.getAtoms().get_params()[ip];
fprintf(fd, "%% %s\n", parname);
}
// write endogenous variables
fprintf(fd, "%%\n%% y ordering\n%% =====================\n");
for (unsigned int ie = 0; ie < model.getAtoms().get_endovars().size(); ie++) {
const char* endoname = model.getAtoms().get_endovars()[ie];
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 (ogp::Constants::Tconstantmap::const_iterator it = cmap.begin();
it != cmap.end(); ++it) {
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 (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
it != lmap.end(); ++it) {
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 (ogp::DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
it != lmap.end(); ++it) {
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 (unsigned int j = 0; j < eam.size(); j++) {
int tvar = variables[eam[j]];
const char* name = model.getAtoms().name(tvar);
int yi = model.getAtoms().name2outer_endo(name);
int t = fder.derivative(ogp::FoldMultiIndex(variables.size(), 1, eam[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.h deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2005-2011, Ondra Kamenik
#ifndef OGDYN_DYNARE_MODEL
#define OGDYN_DYNARE_MODEL
#include "parser/cc/matrix_parser.h"
#include "parser/cc/atom_assignings.h"
#include "dynare_atoms.h"
#include "twod_matrix.h"
#include "Vector.h"
#include "GeneralMatrix.h"
#include <map>
#include <boost/unordered_set.hpp>
namespace ogdyn {
using boost::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() {}
PosInterval(int ifl, int ifc, int ill, int ilc)
: fl(ifl), fc(ifc), ll(ill), lc(ilc) {}
const PosInterval& operator=(const PosInterval& pi)
{fl = pi.fl; fc = pi.fc; ll = pi.ll; lc = pi.lc; return *this;}
/** 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;
/** 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;
/** A vector of initial values created by a subclass. It is
* stored with internal ordering given by atoms. */
Vector* init_vals;
/** A matrix for vcov. It is created by a subclass. */
TwoDMatrix* vcov_mat;
/** Tree index of the planner objective. If there was no
* planner objective keyword, the value is set to -1. */
int t_plobjective;
/** Tree index of the planner discount. If there was no
* planner discount keyword, the value is set to -1. */
int t_pldiscount;
/** Pointer to PlannerBuilder, which is created only if the
* planner's FOC are added to the model. */
PlannerBuilder* pbuilder;
/** 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;
/** 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;
/** Pointer to a copy of original atoms before substitutions
* took place. */
ogp::SAtoms* old_atoms;
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);
virtual ~DynareParser();
DynareModel* clone() const
{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)
: DynareModel(dm) {}
~DynareSPModel() {}
virtual DynareModel* clone() const
{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;
};
/** 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);
virtual ~MatlabSSWriter()
{delete [] id;}
protected:
// from ModelSSWriter
void write_der0_preamble(FILE* fd) const;
void write_der1_preamble(FILE* fd) const;
/** 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;
void write_der0_assignment(FILE* fd) const;
void write_der1_assignment(FILE* fd) const;
/** This prints t10 for t=10. */
void format_term(int t, FILE* fd) const;
/** 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;
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;
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/src/dynare_params.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2004-2011, Ondra Kamenik
#include "dynare_params.h"
#include <getopt.h>
#include <cstdio>
#include <cstring>
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 [2]\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(NULL), num_per(100), num_burn(0), num_sim(80),
num_rtper(0), num_rtsim(0),
num_condper(0), num_condsim(0),
num_threads(2), 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, NULL, opt_per},
{"per", required_argument, NULL, opt_per},
{"burn", required_argument, NULL, opt_burn},
{"simulations", required_argument, NULL, opt_sim},
{"sim", required_argument, NULL, opt_sim},
{"rtperiods", required_argument, NULL, opt_rtper},
{"rtper", required_argument, NULL, opt_rtper},
{"rtsimulations", required_argument, NULL, opt_rtsim},
{"rtsim", required_argument, NULL, opt_rtsim},
{"condperiods", required_argument, NULL, opt_condper},
{"condper", required_argument, NULL, opt_condper},
{"condsimulations", required_argument, NULL, opt_condsim},
{"condsim", required_argument, NULL, opt_condsim},
{"prefix", required_argument, NULL, opt_prefix},
{"threads", required_argument, NULL, opt_threads},
{"steps", required_argument, NULL, opt_steps},
{"seed", required_argument, NULL, opt_seed},
{"order", required_argument, NULL, opt_order},
{"ss-tol", required_argument, NULL, opt_ss_tol},
{"check", required_argument, NULL, opt_check},
{"check-scale", required_argument, NULL, opt_check_scale},
{"check-evals", required_argument, NULL, opt_check_evals},
{"check-num", required_argument, NULL, opt_check_num},
{"qz-criterium",required_argument, NULL, opt_qz_criterium},
{"no-irfs", no_argument, NULL, opt_noirfs},
{"irfs", no_argument, NULL, opt_irfs},
{"centralize", no_argument, NULL, opt_centralize},
{"no-centralize", no_argument, NULL, opt_no_centralize},
{"help", no_argument, NULL, opt_help},
{"version", no_argument, NULL, opt_version},
{NULL, 0, NULL, 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;
}
// Local Variables:
// mode:C++
// End:
dynare++/src/dynare_params.h deleted 100644 → 0
View file @ 6ca649c4
// $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);
};
// Local Variables:
// mode:C++
// End: