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
  • 4.3
  • 4.4
  • 4.5
  • 4.6
  • 5.x
  • 6.x
  • asm
  • aux_func
  • clang+openmp
  • dates-and-dseries-improvements
  • declare_vars_in_model_block
  • dmm
  • dragonfly
  • dynare_minreal
  • eigen
  • error_msg_undeclared_model_vars
  • estim_params
  • exo_steady_state
  • gpm-optimal-policy
  • julia
  • madysson
  • master
  • mex-GetPowerDeriv
  • penalty
  • separateM_
  • slice
  • sphinx-doc-experimental
  • static_aux_vars
  • time-varying-information-set
  • various_fixes
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
  • 4.5.7
  • 4.6-beta1
  • 4.6.0
  • 4.6.0-rc1
  • 4.6.0-rc2
  • 4.6.1
  • 4.6.2
  • 4.6.3
  • 4.6.4
  • 4.7-beta1
  • 4.7-beta2
  • 4.7-beta3
  • 5.0
  • 5.0-rc1
  • 5.1
  • 5.2
  • 5.3
  • 5.4
  • 5.5
  • 6-beta1
  • 6-beta2
  • 6.0
  • 6.1
  • 6.2
  • 6.3
90 results

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
  • 4.3
  • 4.4
  • 4.5
  • DSMH
  • OneStep2
  • SMC
  • SMCsamplers
  • aux_func
  • dates-and-dseries-improvements
  • declare_vars_in_model_block
  • dmm
  • dynamic-striated
  • eigen
  • error_msg_undeclared_model_vars
  • estim_params
  • exceptions
  • exo_steady_state
  • filter_initial_state
  • gpm-optimal-policy
  • julia
  • master
  • merge-initvalfile-fix
  • mex-GetPowerDeriv
  • new_ep
  • nlf-fixes
  • nonlinear-filter-fixes
  • occbin
  • online-filter-as-a-sampler
  • penalty
  • remove-@dates
  • remove-@dseries
  • remove-utilities-tests
  • rmExtraExo
  • separateM_
  • slice
  • smc-sampler
  • sphinx-doc-experimental
  • static_aux_vars
  • temporary_terms
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
74 results
Show changes
Showing
with 0 additions and 4134 deletions
// Copyright (C) 2006, Ondra Kamenik
// $Id: matrix_parser.cpp 2269 2008-11-23 14:33:22Z michel $
#include "parser_exception.hh"
#include "matrix_parser.hh"
#include "location.hh"
#include "matrix_tab.hh"
#include <cstring>
using namespace ogp;
/** A global symbol for passing info to the MatrixParser from
* matrix_parse(). */
MatrixParser *mparser;
/** The declaration of functions defined in matrix_ll.cc and
* matrix_tab.cc generated from matrix.lex and matrix.y. */
void *matrix__scan_buffer(char *, size_t);
void matrix__destroy_buffer(void *);
int matrix_parse();
extern ogp::location_type matrix_lloc;
void
MatrixParser::parse(int length, const char *stream)
{
// reinitialize the object
data.clear();
row_lengths.clear();
nc = 0;
// allocate temporary buffer and parse
auto *buffer = new char[length+2];
strncpy(buffer, stream, length);
buffer[length] = '\0';
buffer[length+1] = '\0';
matrix_lloc.off = 0;
matrix_lloc.ll = 0;
void *p = matrix__scan_buffer(buffer, (unsigned int) length+2);
mparser = this;
matrix_parse();
delete [] buffer;
matrix__destroy_buffer(p);
}
void
MatrixParser::add_item(double v)
{
data.push_back(v);
if (row_lengths.size() == 0)
row_lengths.push_back(0);
(row_lengths.back())++;
if (row_lengths.back() > nc)
nc = row_lengths.back();
}
void
MatrixParser::start_row()
{
row_lengths.push_back(0);
}
void
MatrixParser::error(const char *mes) const
{
throw ParserException(mes, matrix_lloc.off);
}
int
MatrixParser::find_first_non_empty_row(int start) const
{
int r = start;
while (r < (int) row_lengths.size() && row_lengths[r] == 0)
r++;
return r;
}
MPIterator
MatrixParser::begin() const
{
MPIterator it(*this);
return it;
}
MPIterator
MatrixParser::end() const
{
MPIterator it(*this, "end");
return it;
}
MPIterator::MPIterator(const MatrixParser &mp)
: p(&mp), i(0), r(mp.find_first_non_empty_row())
{
}
MPIterator::MPIterator(const MatrixParser &mp, const char *dummy)
: p(&mp), i(mp.data.size()), r(mp.row_lengths.size())
{
}
MPIterator &
MPIterator::operator++()
{
i++;
c++;
if (p->row_lengths[r] <= c)
{
c = 0;
r = p->find_first_non_empty_row(r+1);
}
return *this;
}
// Copyright (C) 2006, Ondra Kamenik
// $Id: matrix_parser.h 762 2006-05-22 13:00:07Z kamenik $
#ifndef OGP_MATRIX_PARSER
#define OGP_MATRIX_PARSER
#include <cstdlib> // For NULL
#include <vector>
namespace ogp
{
using std::vector;
/** This class reads the given string and parses it as a
* matrix. The matrix is read row by row. The row delimiter is
* either a newline character or semicolon (first newline
* character after the semicolon is ignored), the column delimiter
* is either blank character or comma. A different number of items
* in the row is not reconciliated, we do not construct a matrix
* here. The class provides only an iterator to go through all
* read items, the iterator provides information on row number and
* column number of the item. */
class MPIterator;
class MatrixParser
{
friend class MPIterator;
protected:
/** Raw data as they were read. */
vector<double> data;
/** Number of items in each row. */
vector<int> row_lengths;
/** Maximum number of row lengths. */
int nc{0};
public:
MatrixParser()
= default;
MatrixParser(const MatrixParser &mp)
= default;
virtual ~MatrixParser()
= default;
/** Return a number of read rows. */
int
nrows() const
{
return (int) row_lengths.size();
}
/** Return a maximum number of items in the rows. */
int
ncols() const
{
return nc;
}
/** Parses a given data. This initializes the object data. */
void parse(int length, const char *stream);
/** Adds newly read item. This should be called from bison
* parser. */
void add_item(double v);
/** Starts a new row. This should be called from bison
* parser. */
void start_row();
/** Process a parse error from the parser. */
void error(const char *mes) const;
/** Return begin iterator. */
MPIterator begin() const;
/** Return end iterator. */
MPIterator end() const;
protected:
/** Returns an index of the first non-empty row starting at
* start. If the start row is non-empty, returns the start. If
* there is no other non-empty row, returns
* row_lengths.size(). */
int find_first_non_empty_row(int start = 0) const;
};
/** This is an iterator intended to iterate through a matrix parsed
* by MatrixParser. The iterator provides only read-only access. */
class MPIterator
{
friend class MatrixParser;
protected:
/** Reference to the matrix parser. */
const MatrixParser *p{nullptr};
/** The index of the pointed item in the matrix parser. */
unsigned int i{0};
/** The column number of the pointed item starting from zero. */
int c{0};
/** The row number of the pointed item starting from zero. */
int r{0};
public:
MPIterator()
= default;
/** Constructs an iterator pointing to the beginning of the
* parsed matrix. */
MPIterator(const MatrixParser &mp);
/** Constructs an iterator pointing to the past-the-end of the
* parsed matrix. */
MPIterator(const MatrixParser &mp, const char *dummy);
/** Return read-only reference to the pointed item. */
const double &
operator*() const
{
return p->data[i];
}
/** Return a row index of the pointed item. */
int
row() const
{
return r;
}
/** Return a column index of the pointed item. */
int
col() const
{
return c;
}
/** Assignment operator. */
MPIterator &
operator=(const MPIterator &it)
= default;
/** Return true if the iterators are the same, this is if they
* have the same underlying object and the same item index. */
bool
operator==(const MPIterator &it) const
{
return it.p == p && it.i == i;
}
/** Negative of the operator==. */
bool
operator!=(const MPIterator &it) const
{
return !(it == *this);
}
/** Increment operator. */
MPIterator &operator++();
};
};
#endif
// Copyright (C) 2006, Ondra Kamenik
// $Id: namelist.cpp 42 2007-01-22 21:53:24Z ondra $
#include "namelist.hh"
#include <cstring>
using namespace ogp;
/** A global symbol for passing info to NameListParser from its
* parser. */
NameListParser *name_list_parser;
void *namelist__scan_buffer(char *, unsigned int);
void namelist__destroy_buffer(void *);
void namelist_parse();
void
NameListParser::namelist_parse(int length, const char *stream)
{
auto *buffer = new char[length+2];
strncpy(buffer, stream, length);
buffer[length] = '\0';
buffer[length+1] = '\0';
void *p = namelist__scan_buffer(buffer, (unsigned int) length+2);
name_list_parser = this;
::namelist_parse();
delete [] buffer;
namelist__destroy_buffer(p);
}
// Copyright (C) 2007, Ondra Kamenik
// $Id: namelist.h 107 2007-05-10 22:35:04Z ondra $
#ifndef OGP_NAMELIST
#define OGP_NAMELIST
namespace ogp
{
/** Parent class of all parsers parsing a namelist. They must
* implement add_name() method and error() method, which is called
* when an parse error occurs.
*
* Parsing a name list is done as follows: implement
* NameListParser interface, create the object, and call
* NameListParser::namelist_parse(int lengt, const char*
* text). When implementing error(), one may consult global
* location_type namelist_lloc. */
class NameListParser
{
public:
virtual ~NameListParser()
= default;
virtual void add_name(const char *name) = 0;
virtual void namelist_error(const char *mes) = 0;
void namelist_parse(int length, const char *text);
};
};
#endif
/* -*- C++ -*- */
%{
#include "location.hh"
#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);
}
// -*- C++ -*-
// Copyright (C) 2007-2011, Ondra Kamenik
%{
#include "location.hh"
#include "namelist.hh"
#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);
}
// Copyright (C) 2006, Ondra Kamenik
// $Id: parser_exception.cpp 2269 2008-11-23 14:33:22Z michel $
#include "parser_exception.hh"
#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(nullptr),
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(nullptr),
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(nullptr),
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(nullptr),
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(nullptr),
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);
}
// 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
// Copyright (C) 2006, Ondra Kamenik
// $Id: static_atoms.cpp 1360 2007-07-10 11:44:20Z kamenik $
#include "static_atoms.hh"
#include "utils/cc/exception.hh"
using namespace ogp;
StaticAtoms::StaticAtoms(const StaticAtoms &a)
: Atoms(), Constants(a), varnames(a.varnames),
varorder(), vars(), indices()
{
// fill varorder
for (auto i : a.varorder)
{
const char *s = varnames.query(i);
varorder.push_back(s);
}
// fill vars
for (auto var : a.vars)
{
const char *s = varnames.query(var.first);
vars.insert(Tvarmap::value_type(s, var.second));
}
// fill indices
for (auto indice : a.indices)
{
const char *s = varnames.query(indice.second);
indices.insert(Tinvmap::value_type(indice.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);
if (da.is_referenced(name))
{
const DynamicAtoms::Tlagmap &lmap = da.lagmap(name);
for (auto it : lmap)
{
int told = it.second;
tmap.insert(Tintintmap::value_type(told, tnew));
}
}
}
}
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
{
auto it = vars.find(name);
if (it == vars.end())
return -1;
else
return (*it).second;
}
const char *
StaticAtoms::inv_index(int t) const
{
auto it = indices.find(t);
if (it == indices.end())
return nullptr;
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 (auto var : vars)
{
res.push_back(var.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 (auto var : vars)
printf("%s\t->\t%d\n", var.first, var.second);
}
// 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.hh"
namespace ogp
{
class StaticAtoms : public Atoms, public Constants
{
protected:
using Tvarmap = map<const char *, int, ltstr>;
using Tinvmap = map<int, const char *>;
/** 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. */
~StaticAtoms()
override = default;
/** 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 override;
/** 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) override;
int
nvar() const override
{
return varnames.num();
}
/** This returns a vector of all variables. */
vector<int> variables() const override;
/** 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 override;
/** 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
// Copyright (C) 2006, Ondra Kamenik
// $Id: static_fine_atoms.cpp 82 2007-04-19 11:33:30Z ondra $
#include "utils/cc/exception.hh"
#include "static_fine_atoms.hh"
#include "parser_exception.hh"
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 (auto fa_param : fa_params)
register_param(fa_param);
// endogenous
const vector<const char *> &fa_endovars = fa.get_endovars();
for (auto fa_endovar : fa_endovars)
register_endo(fa_endovar);
// exogenous
const vector<const char *> &fa_exovars = fa.get_exovars();
for (auto fa_exovar : fa_exovars)
register_exo(fa_exovar);
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 (auto fa_param : fa_params)
register_param(fa_param);
// 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 == nullptr)
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 (auto & endovar : endovars)
{
int t = index(endovar);
if (t != -1)
{
endo_atoms_map.push_back(der_atoms.size());
der_atoms.push_back(t);
}
}
for (auto & exovar : exovars)
{
int t = index(exovar);
if (t != -1)
{
exo_atoms_map.push_back(der_atoms.size());
der_atoms.push_back(t);
}
}
}
int
StaticFineAtoms::name2outer_param(const char *name) const
{
auto 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
{
auto 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
{
auto 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 == nullptr)
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 == nullptr)
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 == nullptr)
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));
}
// 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.hh"
#include "fine_atoms.hh"
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:
using Tintintmap = map<int, int>;
protected:
using Tvarintmap = map<const char *, int, ltstr>;
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()
= default;
/** 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);
}
~StaticFineAtoms()
override = default;
/** 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 override;
/** 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 override
{
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 override;
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
// Copyright (C) 2005-2011, Ondra Kamenik
#include "utils/cc/exception.hh"
#include "tree.hh"
#include <cstdlib>
#include <cmath>
#include <limits>
using namespace ogp;
/** 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);
auto 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);
auto 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
auto 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
auto 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 (auto & derivative : derivatives)
derivative.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 = erf(r1);
else if (op.getCode() == ERFC)
res = erfc(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;
}
}
// Copyright (C) 2005-2011, Ondra Kamenik
#ifndef OGP_TREE_H
#define OGP_TREE_H
#include <vector>
#include <set>
#include <map>
#include <unordered_map>
#include <unordered_set>
#include <cstdio>
namespace ogp
{
using std::unordered_set;
using std::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{NONE};
/** First operand. If none, then it is -1. */
int op1{-1};
/** Second operand. If none, then it is -1. */
int op2{-1};
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)
{
}
/** Constructs a nulary operation. */
Operation()
= default;
/** A copy constructor. */
Operation(const Operation &op)
= default;
/** Operator =. */
Operation &
operator=(const Operation &op)
= default;
/** 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()
= default;
};
/** 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. */
using _Topmap = unordered_map<Operation, int, ophash>;
using _Topval = _Topmap::value_type;
/** 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. */
using _Tintset = unordered_set<int>;
/** 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. */
using _Tderivmap = unordered_map<int, int>;
/** 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)
= default;
/** 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);
EvalTree(const EvalTree &) = delete;
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;
}
};
/** This is an interface describing how a given operation is
* formatted for output. */
class OperationFormatter
{
public:
/** Empty virtual destructor. */
virtual ~OperationFormatter()
= default;
/** 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) override;
/** 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()
= default;
/** 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()
= default;
/** 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
bin_PROGRAMS = dynare++
GENERATED_FILES = dynglob_ll.cc dynglob_tab.cc dynglob_tab.hh
dynare___SOURCES = \
main.cc \
dynare3.cc \
dynare_atoms.hh \
dynare_model.hh \
forw_subst_builder.hh \
planner_builder.cc \
dynare3.hh \
dynare_exception.hh \
dynare_params.cc \
planner_builder.hh \
dynare_atoms.cc \
dynare_model.cc \
dynare_params.hh \
forw_subst_builder.cc \
nlsolve.cc \
nlsolve.hh \
$(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 = $(AM_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)
dynare___CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
BUILT_SOURCES = $(GENERATED_FILES)
EXTRA_DIST = dynglob.ll dynglob.yy
dynglob_tab.cc dynglob_tab.hh: dynglob.yy
$(YACC) -d -odynglob_tab.cc dynglob.yy
dynglob_ll.cc: dynglob.ll
$(LEX) -i -odynglob_ll.cc dynglob.ll
#include <sstream>
#include <fstream>
#include "dynare3.hh"
#include "dynare_exception.hh"
#include "planner_builder.hh"
#include "forw_subst_builder.hh"
#include "utils/cc/exception.hh"
#include "parser/cc/parser_exception.hh"
#include "parser/cc/atom_substitutions.hh"
#include "../tl/cc/tl_exception.hh"
#include "../kord/kord_exception.hh"
#ifndef DYNVERSION
# define DYNVERSION "unknown"
#endif
/**************************************************************************************/
/* DynareNameList class */
/**************************************************************************************/
std::vector<int>
DynareNameList::selectIndices(const std::vector<const char *> &ns) const
{
std::vector<int> res;
for (auto n : ns)
{
int j = 0;
while (j < getNum() && strcmp(getName(j), n) != 0)
j++;
if (j == getNum())
throw DynareException(__FILE__, __LINE__,
std::string("Couldn't find name for ") + n
+" 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(nullptr), ysteady(nullptr), md(1), dnl(nullptr), denl(nullptr), dsnl(nullptr),
fe(nullptr), fde(nullptr), ss_tol(sstol)
{
std::ifstream f{modname};
if (f.fail())
throw DynareException(__FILE__, __LINE__, std::string{"Could not open model file "}+modname);
std::ostringstream buffer;
buffer << f.rdbuf();
std::string contents{buffer.str()};
try
{
model = new ogdyn::DynareParser(contents.c_str(), contents.length(), ord);
}
catch (const ogp::ParserException &pe)
{
// Compute line and column, given the offset in the file
int line = 1;
int col = 0;
size_t i = 0;
while (i < contents.length() && i < (size_t) pe.offset())
{
if (contents[i] == '\n')
{
line++;
col = 0;
}
i++;
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);
}
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(nullptr), ysteady(nullptr), md(1), dnl(nullptr), denl(nullptr), dsnl(nullptr),
fe(nullptr), fde(nullptr), 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(nullptr),
ysteady(nullptr), md(dynare.md),
dnl(nullptr), denl(nullptr), dsnl(nullptr), fe(nullptr), fde(nullptr),
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 ConstVector &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 ConstVector &yym, const ConstVector &yy,
const ConstVector &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++)
{
auto *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);
}
// $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.hh"
#include "../tl/cc/sparse_tensor.hh"
#include "../kord/decision_rule.hh"
#include "../kord/dynamic_model.hh"
#include "dynare_model.hh"
#include "nlsolve.hh"
#include <vector>
#include <matio.h>
class Dynare;
class DynareNameList : public NameList
{
std::vector<const char *> names;
public:
DynareNameList(const Dynare &dynare);
int
getNum() const override
{
return (int) names.size();
}
const char *
getName(int i) const override
{
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. */
std::vector<int> selectIndices(const std::vector<const char *> &ns) const;
};
class DynareExogNameList : public NameList
{
std::vector<const char *> names;
public:
DynareExogNameList(const Dynare &dynare);
int
getNum() const override
{
return (int) names.size();
}
const char *
getName(int i) const override
{
return names[i];
}
};
class DynareStateNameList : public NameList
{
std::vector<const char *> names;
public:
DynareStateNameList(const Dynare &dynare, const DynareNameList &dnl,
const DynareExogNameList &denl);
int
getNum() const override
{
return (int) names.size();
}
const char *
getName(int i) const override
{
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 override
{
return new Dynare(*this);
}
~Dynare() override;
int
nstat() const override
{
return model->getAtoms().nstat();
}
int
nboth() const override
{
return model->getAtoms().nboth();
}
int
npred() const override
{
return model->getAtoms().npred();
}
int
nforw() const override
{
return model->getAtoms().nforw();
}
int
nexog() const override
{
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 override
{
return model->getOrder();
}
const NameList &
getAllEndoNames() const override
{
return *dnl;
}
const NameList &
getStateNames() const override
{
return *dsnl;
}
const NameList &
getExogNames() const override
{
return *denl;
}
TwoDMatrix &
getVcov()
{
return model->getVcov();
}
const TwoDMatrix &
getVcov() const override
{
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 override
{
return md;
}
const Vector &
getSteady() const override
{
return *ysteady;
}
Vector &
getSteady() override
{
return *ysteady;
}
const ogdyn::DynareModel &
getModel() const
{
return *model;
}
// here is true public interface
void solveDeterministicSteady(Vector &steady);
void
solveDeterministicSteady() override
{
solveDeterministicSteady(*ysteady);
}
void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) override;
void evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy,
const ConstVector &yyp, const Vector &xx) override;
void calcDerivatives(const Vector &yy, const Vector &xx);
void calcDerivativesAtSteady() override;
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) override
{
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) override;
};
class DynareJacobian : public ogu::Jacobian, public ogp::FormulaDerEvalLoader
{
protected:
Dynare &d;
public:
DynareJacobian(Dynare &dyn);
~DynareJacobian()
override = default;
void load(int i, int iord, const int *vars, double res) override;
void eval(const Vector &in) override;
};
class DynareVectorFunction : public ogu::VectorFunction
{
protected:
Dynare &d;
public:
DynareVectorFunction(Dynare &dyn)
: d(dyn)
{
}
~DynareVectorFunction()
override = default;
int
inDim() const override
{
return d.ny();
}
int
outDim() const override
{
return d.ny();
}
void eval(const ConstVector &in, Vector &out) override;
};
#endif
// Copyright (C) 2006, Ondra Kamenik
// $Id: dynare_atoms.cpp 1765 2008-03-31 14:32:08Z kamenik $
#include "parser/cc/parser_exception.hh"
#include "utils/cc/exception.hh"
#include "dynare_atoms.hh"
#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 (nullptr == varnames.query(name))
throw ogp::ParserException(std::string("Unknown name <")+name+">", 0);
auto 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 (auto it : dda.atom_type)
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
{
auto 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 (auto it : atom_type)
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++)
{
if (atoms.is_referenced(atoms.get_params()[i]))
{
const ogp::DynamicAtoms::Tlagmap &lmap = atoms.lagmap(atoms.get_params()[i]);
for (auto it : lmap)
{
int t = it.second;
et.set_nulary(t, paramvals[i]);
}
}
}
// set endogenous
for (unsigned int outer_i = 0; outer_i < atoms.get_endovars().size(); outer_i++)
{
if (atoms.is_referenced(atoms.get_endovars()[outer_i]))
{
const ogp::DynamicAtoms::Tlagmap &lmap = atoms.lagmap(atoms.get_endovars()[outer_i]);
for (auto it : lmap)
{
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()]);
}
}
}
// set exogenous
for (unsigned int outer_i = 0; outer_i < atoms.get_exovars().size(); outer_i++)
{
if (atoms.is_referenced(atoms.get_exovars()[outer_i]))
{
const ogp::DynamicAtoms::Tlagmap &lmap = atoms.lagmap(atoms.get_exovars()[outer_i]);
for (auto it : lmap)
{
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]);
}
}
}
}
}
void
DynareStaticSteadyAtomValues::setValues(ogp::EvalTree &et) const
{
// set constants
atoms_static.setValues(et);
// set parameters
for (auto name : atoms_static.get_params())
{
int t = atoms_static.index(name);
if (t != -1)
{
int idyn = atoms.name2outer_param(name);
et.set_nulary(t, paramvals[idyn]);
}
}
// set endogenous
for (auto name : atoms_static.get_endovars())
{
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 (auto name : atoms_static.get_exovars())
{
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 (auto it : subst)
{
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 (auto it : subst)
{
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;
}
// 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.hh"
#include "parser/cc/static_atoms.hh"
#include "parser/cc/static_fine_atoms.hh"
#include "parser/cc/atom_substitutions.hh"
#include "parser/cc/tree.hh"
#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. */
using Tsubstmap = map<const char *, int, ogp::ltstr>;
class DynareStaticAtoms : public ogp::StaticAtoms
{
public:
DynareStaticAtoms()
: StaticAtoms()
{
}
DynareStaticAtoms(const DynareStaticAtoms &a)
= default;
~DynareStaticAtoms()
override = default;
/** 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) override;
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 override;
};
class DynareDynamicAtoms : public ogp::SAtoms, public ogp::NularyStringConvertor
{
public:
enum atype {endovar, exovar, param};
protected:
using Tatypemap = map<const char *, atype, ogp::ltstr>;
/** The map assigining a type to each name. */
Tatypemap atom_type;
public:
DynareDynamicAtoms()
: ogp::SAtoms()
{
}
DynareDynamicAtoms(const DynareDynamicAtoms &dda);
~DynareDynamicAtoms()
override = default;
/** This parses a variable of the forms: varname(+3),
* varname(3), varname, varname(-3), varname(0), varname(+0),
* varname(-0). */
void parse_variable(const char *in, std::string &out, int &ll) const override;
/** Registers unique name of endogenous variable. */
void register_uniq_endo(const char *name) override;
/** Registers unique name of exogenous variable. */
void register_uniq_exo(const char *name) override;
/** Registers unique name of parameter. */
void register_uniq_param(const char *name) override;
/** Return true if the name is a given type. */
bool is_type(const char *name, atype tp) const;
/** Debug print. */
void print() const override;
/** Implement NularyStringConvertor::convert. */
std::string convert(int t) const override;
};
/** 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 ConstVector &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 override;
};
/** 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 override
{
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 override;
};
/** 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) override;
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) override;
protected:
Vector &y;
vector<const char *> left_hand_sides;
vector<int> right_hand_sides;
};
};
#endif
// 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