Skip to content
Snippets Groups Projects

Compare revisions

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

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Show changes
Showing
with 0 additions and 3725 deletions
// Copyright (C) 2006, Ondra Kamenik
// $Id: atom_substitutions.h 42 2007-01-22 21:53:24Z ondra $
#ifndef OGP_ATOM_SUBSTITUTIONS_H
#define OGP_ATOM_SUBSTITUTIONS_H
#include "fine_atoms.h"
#include <string>
namespace ogp {
using std::string;
using std::map;
using std::pair;
/** This class tracts an information about the performed
* substitutions. In fact, there is only one number to keep track
* about, this is a number of substitutions. */
struct SubstInfo {
int num_substs;
SubstInfo() : num_substs(0) {}
};
/** This class tracks all atom substitutions during the job and
* then builds structures when all substitutions are finished. */
class AtomSubstitutions {
public:
typedef pair<const char*, int> Tshiftname;
typedef map<const char*, Tshiftname, ltstr> Tshiftmap;
typedef set<Tshiftname> Tshiftnameset;
typedef map<const char*, Tshiftnameset, ltstr> Toldnamemap;
protected:
/** This maps a new name to a shifted old name. This is, one
* entry looks as "a_m3 ==> a(-3)", saying that a variable
* "a_m3" corresponds to a variable "a" lagged by 3. */
Tshiftmap new2old;
/** This is inverse to new2old, which is not unique. For old
* name, say "a", it says what new names are derived with what
* shifts from the "a". For example, it can map "a" to a two
* element set {["a_m3", +3], ["a_p2", -2]}. This says that
* leading "a_m3" by 3 one gets old "a" and lagging "a_p2" by
* 2 one gets also old "a". */
Toldnamemap old2new;
/** This is a reference to old atoms with multiple leads and
* lags. They are supposed to be used with parsing finished
* being had called, so that the external ordering is
* available. */
const FineAtoms& old_atoms;
/** This is a reference to new atoms. All name pointers point
* to storage of these atoms. */
FineAtoms& new_atoms;
/** Substitutions information. */
SubstInfo info;
public:
/** Create the object with reference to the old and new
* atoms. In the beginning, old atoms are supposed to be with
* parsing_finished() called, and new atoms a simple copy of
* old atoms. The new atoms will be an instance of SAtoms. All
* substitution job is done by a substitution method of the
* new atoms. */
AtomSubstitutions(const FineAtoms& oa, FineAtoms& na)
: old_atoms(oa), new_atoms(na) {}
/** Construct a copy of the object using a different instances
* of old atoms and new atoms, which are supposed to be
* semantically same as the atoms from as. */
AtomSubstitutions(const AtomSubstitutions& as, const FineAtoms& oa, FineAtoms& na);
virtual ~AtomSubstitutions() {}
/** This is called during the substitution job from the
* substitution method of the new atoms. This says that the
* new name, say "a_m3" is a substitution of old name "a"
* shifted by -3. */
void add_substitution(const char* newname, const char* oldname, int tshift);
/** This is called when all substitutions are finished. This
* forms the new external ordering of the new atoms and calls
* parsing_finished() for the new atoms with the given ordering type. */
void substitutions_finished(VarOrdering::ord_type ot);
/** Returns a new name for old name and given tshift. For "a"
* and tshift=-3, it returns "a_m3". If there is no such
* substitution, it returns NULL. */
const char* get_new4old(const char* oldname, int tshift) const;
/** Return new2old. */
const Tshiftmap& get_new2old() const
{return new2old;}
/** Return old2new. */
const Toldnamemap& get_old2new() const
{return old2new;}
/** Return substitution info. */
const SubstInfo& get_info() const
{return info;}
/** Return old atoms. */
const FineAtoms& get_old_atoms() const
{return old_atoms;}
/** Return new atoms. */
const FineAtoms& get_new_atoms() const
{return new_atoms;}
/** Debug print. */
void print() const;
};
class SAtoms : public FineAtoms {
public:
SAtoms()
: FineAtoms() {}
SAtoms(const SAtoms& sa)
: FineAtoms(sa) {}
virtual ~SAtoms() {}
/** This substitutes all lags and leads for all exogenous and
* all lags and leads greater than 1 for all endogenous
* variables. This is useful for perfect foresight problems
* where we can do that. */
void substituteAllLagsAndLeads(FormulaParser& fp, AtomSubstitutions& as);
/** This substitutes all lags of all endo and exo and one step
* leads of all exo variables. This is useful for stochastic
* models where we cannot solve leads more than 1. */
void substituteAllLagsAndExo1Leads(FormulaParser& fp, AtomSubstitutions& as);
protected:
/** This finds an endogenous variable name which occurs between
* ll1 and ll2 included. */
const char* findEndoWithLeadInInterval(int ll1, int ll2) const
{return findNameWithLeadInInterval(get_endovars(), ll1, ll2);}
/** This finds an exogenous variable name which occurs between
* ll1 and ll2 included. */
const char* findExoWithLeadInInterval(int ll1, int ll2) const
{return findNameWithLeadInInterval(get_exovars(), ll1, ll2);}
/** This attempts to find a non registered name of the form
* <str>_m<abs(ll)> or <str>_p<abs(ll)>. A letter 'p' is
* chosen if ll is positive, 'm' if negative. If a name of
* such form is already registered, one more character (either
* 'p' or 'm') is added and the test is performed again. The
* resulting name is returned in a string out. */
void attemptAuxName(const char* str, int ll, string& out) const;
/** This makes auxiliary variables to eliminate all leads/lags
* greater/less than or equal to start up to the limit_lead
* for a variable with the given name. If the limit_lead is
* greater/less than the maxlead/minlag of the variable, than
* maxlead/minlag is used. This process is recorded in
* AtomSubstitutions. The new auxiliary variables and their
* atoms are created in this object. The auxiliary equations
* are created in the given FormulaParser. The value of step
* is allowed to be either -1 (lags) or +1 (leads). */
void makeAuxVariables(const char* name, int step, int start, int limit_lead,
FormulaParser& fp, AtomSubstitutions& as);
private:
/** This is a worker routine for findEndoWithLeadInInterval
* and findExoWithLeadInInterval. */
const char* findNameWithLeadInInterval(const vector<const char*>& names,
int ll1, int ll2) const;
};
};
#endif
// Local Variables:
// mode:C++
// End:
%{
// Copyright (C) 2007-2011, Ondra Kamenik
#include "location.h"
#include "csv_tab.hh"
extern YYLTYPE csv_lloc;
#define YY_USER_ACTION SET_LLOC(csv_);
%}
%option nounput
%option noyy_top_state
%option yylineno
%option prefix="csv_"
%option never-interactive
%%
, {return COMMA;}
\n {return NEWLINE;}
\r\n {return NEWLINE;}
[^,\n\r]+ {return ITEM;}
%%
int csv_wrap()
{
return 1;
}
void csv__destroy_buffer(void* p)
{
csv__delete_buffer((YY_BUFFER_STATE)p);
}
%{
#include "location.h"
#include "csv_parser.h"
#include "csv_tab.hh"
void csv_error(const char*);
int csv_lex(void);
extern int csv_lineno;
extern ogp::CSVParser* csv_parser;
extern YYLTYPE csv_lloc;
%}
%union {
char* string;
int integer;
}
%token COMMA NEWLINE BOGUS
%token <string> ITEM
%name-prefix="csv_";
%locations
%error-verbose
%%
csv_file : line_list | line_list line;
line_list : line_list line newline | line newline | line_list newline | newline;
line : line comma | line item | item | comma;
comma : COMMA {csv_parser->nextcol();};
newline : NEWLINE {csv_parser->nextrow();};
item : ITEM {csv_parser->item(@1.off, @1.ll);};
%%
void csv_error(const char* mes)
{
csv_parser->csv_error(mes);
}
#include "csv_parser.h"
#include "parser_exception.h"
#include "location.h"
#include "csv_tab.hh"
#include <cstring>
using namespace ogp;
/** A global symbol for passing info to the CSVParser from
* csv_parse(). */
CSVParser* csv_parser;
/** The declaration of functions defined in csv_ll.cc and
* csv_tab.cc generated from csv.lex and csv.y. */
void* csv__scan_buffer(char*, unsigned int);
void csv__destroy_buffer(void*);
int csv_parse();
extern ogp::location_type csv_lloc;
void CSVParser::csv_error(const char* mes)
{
throw ParserException(mes, csv_lloc.off);
}
void CSVParser::csv_parse(int length, const char* str)
{
// allocate temporary buffer and parse
char* buffer = new char[length+2];
strncpy(buffer, str, length);
buffer[length] = '\0';
buffer[length+1] = '\0';
csv_lloc.off = 0;
csv_lloc.ll = 0;
parsed_string = buffer;
void* p = csv__scan_buffer(buffer, (unsigned int)length+2);
csv_parser = this;
::csv_parse();
delete [] buffer;
csv__destroy_buffer(p);
parsed_string = NULL;
}
// Copyright (C) 2007, Ondra Kamenik
// $Id$
#ifndef OGP_CSV_PARSER
#define OGP_CSV_PARSER
namespace ogp {
class CSVParserPeer {
public:
virtual ~CSVParserPeer() {}
virtual void item(int irow, int icol, const char* str, int length) = 0;
};
class CSVParser {
private:
CSVParserPeer& peer;
int row;
int col;
const char* parsed_string;
public:
CSVParser(CSVParserPeer& p)
: peer(p), row(0), col(0), parsed_string(0) {}
CSVParser(const CSVParser& csvp)
: peer(csvp.peer), row(csvp.row),
col(csvp.col), parsed_string(csvp.parsed_string) {}
virtual ~CSVParser() {}
void csv_error(const char* mes);
void csv_parse(int length, const char* str);
void nextrow()
{row++; col = 0;}
void nextcol()
{col++;}
void item(int off, int length)
{peer.item(row, col, parsed_string+off, length);}
};
};
#endif
// Local Variables:
// mode:C++
// End:
// Copyright (C) 2005, Ondra Kamenik
// $Id: dynamic_atoms.cpp 1362 2007-07-10 11:50:18Z kamenik $
#include "utils/cc/exception.h"
#include "dynamic_atoms.h"
#include <cstring>
#include <climits>
using namespace ogp;
NameStorage::NameStorage(const NameStorage& stor)
{
for (unsigned int i = 0; i < stor.name_store.size(); i++) {
char* str = new char[strlen(stor.name_store[i])+1];
strcpy(str, stor.name_store[i]);
name_store.push_back(str);
name_set.insert(str);
}
}
NameStorage::~NameStorage()
{
while (name_store.size() > 0) {
delete [] name_store.back();
name_store.pop_back();
}
}
const char* NameStorage::query(const char* name) const
{
set<const char*, ltstr>::const_iterator it = name_set.find(name);
if (it == name_set.end())
return NULL;
else
return (*it);
}
const char* NameStorage::insert(const char* name)
{
set<const char*, ltstr>::const_iterator it = name_set.find(name);
if (it == name_set.end()) {
char* str = new char[strlen(name)+1];
strcpy(str, name);
name_store.push_back(str);
name_set.insert(str);
return str;
} else {
return (*it);
}
}
void NameStorage::print() const
{
for (unsigned int i = 0; i < name_store.size(); i++)
printf("%s\n", name_store[i]);
}
void Constants::import_constants(const Constants& c, OperationTree& otree, Tintintmap& tmap)
{
for (Tconstantmap::const_iterator it = c.cmap.begin();
it != c.cmap.end(); ++it) {
int told = (*it).first;
int tnew = otree.add_nulary();
tmap.insert(Tintintmap::value_type(told, tnew));
add_constant(tnew, (*it).second);
}
}
void Constants::setValues(EvalTree& et) const
{
Tconstantmap::const_iterator it;
for (it = cmap.begin(); it != cmap.end(); ++it)
et.set_nulary((*it).first, (*it).second);
}
void Constants::add_constant(int t, double val)
{
cmap.insert(Tconstantmap::value_type(t, val));
cinvmap.insert(Tconstantinvmap::value_type(val, t));
}
bool Constants::is_constant(int t) const
{
if (t < OperationTree::num_constants)
return true;
Tconstantmap::const_iterator it = cmap.find(t);
return (it != cmap.end());
}
double Constants::get_constant_value(int t) const
{
Tconstantmap::const_iterator it = cmap.find(t);
if (it != cmap.end())
return (*it).second;
else {
throw ogu::Exception(__FILE__,__LINE__,
"Tree index is not constant in Constants::get_constant_value");
return 0;
}
}
int Constants::check(const char* str) const
{
double d;
sscanf(str, "%lf", &d);
Tconstantinvmap::const_iterator it = cinvmap.find(d);
if (it != cinvmap.end())
return (*it).second;
else
return -1;
}
void Constants::print() const
{
Tconstantmap::const_iterator it;
for (it = cmap.begin(); it != cmap.end(); ++it)
printf("$%d: %8.4g\n", (*it).first, (*it).second);
}
DynamicAtoms::DynamicAtoms()
: nv(0), minlag(INT_MAX), maxlead(INT_MIN)
{
}
DynamicAtoms::DynamicAtoms(const DynamicAtoms& da)
: Constants(da),
varnames(da.varnames), vars(), indices(),
nv(da.nv), minlag(da.minlag), maxlead(da.maxlead)
{
// copy vars
for (Tvarmap::const_iterator it = da.vars.begin();
it != da.vars.end(); ++it)
vars.insert(Tvarmap::value_type(varnames.query((*it).first),
(*it).second));
// copy indices
for (Tindexmap::const_iterator it = da.indices.begin();
it != da.indices.end(); ++it)
indices.insert(Tindexmap::value_type((*it).first,
varnames.query((*it).second)));
}
int DynamicAtoms::check(const char* name) const
{
if (is_string_constant(name))
return Constants::check(name);
return check_variable(name);
}
int DynamicAtoms::check_variable(const char* name) const
{
string str;
int ll;
parse_variable(name, str, ll);
Tvarmap::const_iterator it = vars.find(str.c_str());
if (it != vars.end()) {
const Tlagmap& lmap = (*it).second;
Tlagmap::const_iterator itt = lmap.find(ll);
if (itt != lmap.end())
return (*itt).second;
}
return -1;
}
void DynamicAtoms::assign(const char* name, int t)
{
if (is_string_constant(name))
assign_constant(name, t);
else
assign_variable(name, t);
}
void DynamicAtoms::assign_constant(const char* name, int t)
{
double val;
sscanf(name, "%lf", &val);
add_constant(t, val);
}
// parse the name and then call assing_variable(varname, ll, t)
void DynamicAtoms::assign_variable(const char* name, int t)
{
int ll;
string str;
parse_variable(name, str, ll);
// here str is just name without lead/lag
const char* ss = varnames.insert(str.c_str());
assign_variable(ss, ll, t);
}
void DynamicAtoms::assign_variable(const char* varname, int ll, int t)
{
if (indices.end() != indices.find(t))
throw ogu::Exception(__FILE__,__LINE__,
"Attempt to assign already allocated tree index");
Tvarmap::iterator it = vars.find(varname);
if (it != vars.end()) {
Tlagmap& lmap = (*it).second;
if (lmap.end() != lmap.find(ll))
throw ogu::Exception(__FILE__,__LINE__,
"Attempt to assign already allocated variable");
lmap.insert(Tlagmap::value_type(ll, t));
} else {
Tlagmap lmap;
lmap.insert(Tlagmap::value_type(ll, t));
vars.insert(Tvarmap::value_type(varname, lmap));
}
indices.insert(Tindexmap::value_type(t, varname));
nv++;
if (ll < minlag)
minlag = ll;
if (ll > maxlead)
maxlead = ll;
}
void DynamicAtoms::unassign_variable(const char* varname, int ll, int t)
{
Tvarmap::iterator it = vars.find(varname);
if (it != vars.end()) {
Tlagmap& lmap = (*it).second;
Tlagmap::iterator itt = lmap.find(ll);
if (itt != lmap.end()) {
if ((*itt).second == t) {
// erase it from the lagmap; if it becomes empty,
// erase the lagmap from varmap
lmap.erase(itt);
if (lmap.size() == 0)
vars.erase(it);
// erase it from the indices
Tindexmap::iterator ittt = indices.find(t);
if (ittt != indices.end())
indices.erase(ittt);
nv--;
if (ll == minlag || ll == maxlead)
update_minmaxll();
} else
throw ogu::Exception(__FILE__,__LINE__,
"Tree index inconsistent in DynamicAtoms::unassign_variable");
} else
throw ogu::Exception(__FILE__,__LINE__,
"Lead/lag of the variable not found in DynamicAtoms::unassign_variable");
} else
throw ogu::Exception(__FILE__,__LINE__,
"Variable not found in DynamicAtoms::unassign_variable");
}
void DynamicAtoms::update_minmaxll()
{
minlag = INT_MAX;
maxlead =INT_MIN;
for (Tvarmap::const_iterator it = vars.begin(); it != vars.end(); ++it) {
const Tlagmap& lmap = (*it).second;
for (Tlagmap::const_iterator itt = lmap.begin(); itt != lmap.end(); ++itt) {
int ll = (*itt).first;
if (ll < minlag)
minlag = ll;
if (ll > maxlead)
maxlead = ll;
}
}
}
vector<int> DynamicAtoms::variables() const
{
vector<int> res;
for (Tvarmap::const_iterator it = vars.begin();
it != vars.end(); ++it) {
const Tlagmap& lmap = (*it).second;
for (Tlagmap::const_iterator itt = lmap.begin();
itt != lmap.end(); ++itt)
res.push_back((*itt).second);
}
return res;
}
void DynamicAtoms::varspan(int t, int& mlead, int& mlag) const
{
Tindexmap::const_iterator it = indices.find(t);
if (indices.end() == it) {
mlead = INT_MIN;
mlag = INT_MAX;
return;
}
varspan((*it).second, mlead, mlag);
}
void DynamicAtoms::varspan(const char* name, int& mlead, int& mlag) const
{
Tvarmap::const_iterator it = vars.find(name);
if (vars.end() == it) {
mlead = INT_MIN;
mlag = INT_MAX;
return;
}
const Tlagmap& lmap = (*it).second;
Tlagmap::const_iterator beg = lmap.begin();
Tlagmap::const_reverse_iterator end = lmap.rbegin();
mlag = (*beg).first;
mlead = (*end).first;
}
void DynamicAtoms::varspan(const vector<const char*>& names, int& mlead, int& mlag) const
{
mlead = INT_MIN;
mlag = INT_MAX;
for (unsigned int i = 0; i < names.size(); i++) {
int lag, lead;
varspan(names[i], lead, lag);
if (lead > mlead)
mlead = lead;
if (lag < mlag)
mlag = lag;
}
}
bool DynamicAtoms::is_named_atom(int t) const
{
return (indices.end() != indices.find(t));
}
int DynamicAtoms::index(const char* name, int ll) const
{
Tvarmap::const_iterator it = vars.find(name);
if (vars.end() != it) {
const Tlagmap& lmap = (*it).second;
Tlagmap::const_iterator itt = lmap.find(ll);
if (lmap.end() != itt)
return (*itt).second;
}
return -1;
}
const DynamicAtoms::Tlagmap& DynamicAtoms::lagmap(const char* name) const
{
Tvarmap::const_iterator it = vars.find(name);
if (vars.end() == it)
throw ogu::Exception(__FILE__,__LINE__,
std::string("Couldn't find the name ")
+ name + " in DynamicAtoms::lagmap");
return (*it).second;
}
const char* DynamicAtoms::name(int t) const
{
Tindexmap::const_iterator it = indices.find(t);
if (indices.end() == it)
throw ogu::Exception(__FILE__,__LINE__,
"Couldn't find tree index in DynamicAtoms::name");
return (*it).second;
}
int DynamicAtoms::lead(int t) const
{
const char* nam = name(t);
const Tlagmap& lmap = lagmap(nam);
Tlagmap::const_iterator it = lmap.begin();
while (it != lmap.end() && (*it).second != t)
++it;
if (lmap.end() == it)
throw ogu::Exception(__FILE__,__LINE__,
"Couldn't find the three index in DynamicAtoms::lead");
return (*it).first;
}
void DynamicAtoms::print() const
{
printf("names:\n");
varnames.print();
printf("constants:\n");
Constants::print();
printf("variables:\n");
for (Tvarmap::const_iterator it = vars.begin();
it != vars.end(); ++it) {
const Tlagmap& lmap = (*it).second;
for (Tlagmap::const_iterator itt = lmap.begin();
itt != lmap.end(); ++itt)
printf("$%d: %s(%d)\n", (*itt).second, (*it).first, (*itt).first);
}
printf("indices:\n");
for (Tindexmap::const_iterator it = indices.begin();
it != indices.end(); ++it)
printf("t=%d ==> %s\n", (*it).first, (*it).second);
}
/** Note that the str has been parsed by the lexicographic
* analyzer. It can be either a variable or a double. So it is easy to
* recognize it by the first character. */
bool DynamicAtoms::is_string_constant(const char* str)
{
return str[0] == '.' || str[0] == '-' || (str[0] >= '0' && str[0] <= '9');
}
VarOrdering::VarOrdering(const VarOrdering& vo, const vector<const char*>& vnames,
const DynamicAtoms& a)
: n_stat(vo.n_stat), n_pred(vo.n_pred), n_both(vo.n_both), n_forw(vo.n_forw),
der_atoms(vo.der_atoms), positions(vo.positions),
outer2y(vo.outer2y), y2outer(vo.y2outer), varnames(vnames), atoms(a)
{
}
bool VarOrdering::check(int t) const
{
map<int,int>::const_iterator it = positions.find(t);
return it != positions.end();
}
int VarOrdering::get_pos_of(int t) const
{
map<int,int>::const_iterator it = positions.find(t);
if (it != positions.end()) {
return (*it).second;
} else {
throw ogu::Exception(__FILE__,__LINE__,
"Couldn't find the tree index in VarOrdering::get_pos_of");
return -1;
}
}
void VarOrdering::do_general(ord_type ordering)
{
// auxiliary vectors for setting der_atoms and map
vector<int> pred_minus;
vector<int> both_minus;
vector<int> stat;
vector<int> pred_pad;
vector<int> both_pad;
vector<int> forw_pad;
vector<int> both_plus;
vector<int> forw_plus;
// auxiliary vectors for setting y2outer and outer2y
vector<int> y2o_stat;
vector<int> y2o_pred;
vector<int> y2o_both;
vector<int> y2o_forw;
for (unsigned int i = 0; i < varnames.size(); i++) {
const char* ss = varnames[i];
int lead;
int lag;
atoms.varspan(ss, lead, lag);
if (lag == 0 && lead == 0) {
stat.push_back(atoms.index(ss, 0));
y2o_stat.push_back(i);
} else if (lag == -1 && lead < 1) {
pred_minus.push_back(atoms.index(ss, -1));
pred_pad.push_back(atoms.index(ss, 0));
y2o_pred.push_back(i);
} else if (lag > -1 && lead == 1) {
forw_pad.push_back(atoms.index(ss, 0));
forw_plus.push_back(atoms.index(ss, 1));
y2o_forw.push_back(i);
} else if (lag == -1 && lead == 1) {
both_minus.push_back(atoms.index(ss, -1));
both_pad.push_back(atoms.index(ss, 0));
both_plus.push_back(atoms.index(ss, 1));
y2o_both.push_back(i);
} else {
throw ogu::Exception(__FILE__,__LINE__,
"A wrong lag/lead of a variable in VarOrdering::do_pbspbfbf");
}
}
// here we fill ords according to ordering
vector<int>* ords[8];
if (ordering == pbspbfbf) {
ords[0] = &pred_minus;
ords[1] = &both_minus;
ords[2] = &stat;
ords[3] = &pred_pad;
ords[4] = &both_pad;
ords[5] = &forw_pad;
ords[6] = &both_plus;
ords[7] = &forw_plus;
} else if (ordering == bfspbfpb) {
ords[0] = &both_plus;
ords[1] = &forw_plus;
ords[2] = &stat;
ords[3] = &pred_pad;
ords[4] = &both_pad;
ords[5] = &forw_pad;
ords[6] = &pred_minus;
ords[7] = &both_minus;
} else { // BEWARE: when implementing a new ordering, check also a
// code below setting y2outer
throw ogu::Exception(__FILE__,__LINE__,
"Ordering not implemented in VarOrdering::do_general");
}
// make der_atoms and positions
int off = 0;
for (unsigned int i = 0; i < 8; i++)
for (unsigned int j = 0; j < (ords[i])->size(); j++, off++)
if ((*(ords[i]))[j] != -1) {
der_atoms.push_back((*(ords[i]))[j]);
positions.insert(std::pair<int,int>((*(ords[i]))[j], off));
}
// set integer constants
n_stat = stat.size();
n_pred = pred_pad.size();
n_both = both_pad.size();
n_forw = forw_pad.size();
// make y2outer mapping
y2outer.insert(y2outer.end(), y2o_stat.begin(), y2o_stat.end());
y2outer.insert(y2outer.end(), y2o_pred.begin(), y2o_pred.end());
y2outer.insert(y2outer.end(), y2o_both.begin(), y2o_both.end());
y2outer.insert(y2outer.end(), y2o_forw.begin(), y2o_forw.end());
// make outer2y mapping
outer2y.resize(y2outer.size(), -1);
for (unsigned int i = 0; i < y2outer.size(); i++)
outer2y[y2outer[i]] = i;
}
void VarOrdering::do_increasing_time()
{
// get maxlead and minlag of the variables
int mlag, mlead;
atoms.varspan(varnames, mlead, mlag);
// setup the matrix of tree indices, if there is no occurrence,
// the index is set to -1
vector<int> ll_init(varnames.size(), -1);
vector<vector<int> > tree_ind(mlead-mlag+1, ll_init);
for (unsigned int iv = 0; iv < varnames.size(); iv++) {
try {
const DynamicAtoms::Tlagmap& lmap = atoms.lagmap(varnames[iv]);
for (DynamicAtoms::Tlagmap::const_iterator it = lmap.begin();
it != lmap.end(); ++it) {
int ll = (*it).first;
int t = (*it).second;
tree_ind[ll-mlag][iv] = t;
}
} catch (const ogu::Exception& e) {
// ignore the error of not found variable in the tree
}
}
// setup der_atoms and positions
for (int ll = mlag; ll <= mlead; ll++)
for (unsigned int iv = 0; iv < varnames.size(); iv++) {
int t = tree_ind[ll-mlag][iv];
if (t != -1) {
der_atoms.push_back(t);
int pos = (ll-mlag)*varnames.size() + iv;
positions.insert(map<int,int>::value_type(t, pos));
}
}
// set outer2y and y2outer to identities
for (unsigned int iv = 0; iv < varnames.size(); iv++) {
outer2y.push_back(iv);
y2outer.push_back(iv);
}
// set n_stat, n_pred, n_both, and n_forw
for (unsigned int iv = 0; iv < varnames.size(); iv++) {
int mmlag, mmlead;
atoms.varspan(varnames[iv], mmlead, mmlag);
if (mmlead == 0 && mmlag == 0) {
n_stat++;
} else if (mmlead <= 0 && mmlag < 0) {
n_pred++;
} else if (mmlead > 0 && mmlag >=0) {
n_forw++;
} else if (mmlead > 0 && mmlag < 0) {
n_both++;
} else if (mmlead < mmlag) {
// variable does not occur in the tree, cound as static
n_stat++;
} else {
throw ogu::Exception(__FILE__,__LINE__,
"A wrong lag/lead of a variable in VarOrdering::do_increasing_time");
}
}
}
void VarOrdering::print() const
{
printf("nstat=%d, npred=%d, nboth=%d, nforw=%d\n", n_stat, n_pred, n_both, n_forw);
printf("der_atoms:\n");
for (unsigned int i = 0; i < der_atoms.size(); i++)
printf(" %d", der_atoms[i]);
printf("\nmap:\n");
for (map<int,int>::const_iterator it = positions.begin(); it != positions.end(); ++it)
printf(" [%d->%d]", (*it).first, (*it).second);
printf("\ny2outer:\n");
for (unsigned int i = 0; i < y2outer.size(); i++)
printf(" %d", y2outer[i]);
printf("\nouter2y:\n");
for (unsigned int i = 0; i < outer2y.size(); i++)
printf(" %d", outer2y[i]);
printf("\n");
}
// Local Variables:
// mode:C++
// End:
// Copyright (C) 2005, Ondra Kamenik
// $Id: dynamic_atoms.h 2269 2008-11-23 14:33:22Z michel $
#ifndef OGP_DYNAMIC_ATOMS_H
#define OGP_DYNAMIC_ATOMS_H
#include "formula_parser.h"
#include <vector>
#include <map>
#include <set>
#include <string>
#include <cstring>
namespace ogp {
using std::vector;
using std::map;
using std::set;
using std::string;
struct ltstr {
bool operator()(const char* a1, const char* a2) const
{ return strcmp(a1, a2) < 0; }
};
/** Class storing names. We will keep names of variables in
* various places, and all these pointers will point to one
* storage, which will be responsible for allocation and
* deallocation. The main function of the class is to allocate
* space for names, and return a pointer of the stored name if
* required. */
class NameStorage {
protected:
/** Vector of names allocated, this is the storage. */
vector<char*> name_store;
/** Map useful to quickly decide if the name is already
* allocated or not. */
set<const char*, ltstr> name_set;
public:
NameStorage() {}
NameStorage(const NameStorage& stor);
virtual ~NameStorage();
/** Query for the name. If the name has been stored, it
* returns its address, otherwise 0. */
const char* query(const char* name) const;
/** Insert the name if it has not been inserted yet, and
* return its new or old allocation. */
const char* insert(const char* name);
int num() const
{return (int)name_store.size();}
const char* get_name(int i) const
{return name_store[i];}
/** Debug print. */
void print() const;
};
class Constants : public AtomValues {
public:
/** Type for a map mapping tree indices to double values. */
typedef map<int,double> Tconstantmap;
typedef map<int,int> Tintintmap;
protected:
/** Map mapping a tree index of a constant to its double value. */
Tconstantmap cmap;
public:
Constants() {}
/** Copy constructor. */
Constants(const Constants& c)
: cmap(c.cmap), cinvmap(c.cinvmap) {}
/** Copy constructor registering the constants in the given
* tree. The mapping from old tree indices to new ones is
* traced in tmap. */
Constants(const Constants& c, OperationTree& otree, Tintintmap& tmap)
{import_constants(c, otree, tmap);}
/** Import constants registering their tree indices in the
* given tree. The mapping form old tree indices to new ones
* is traced in tmap. */
void import_constants(const Constants& c, OperationTree& otree, Tintintmap& tmap);
/** Implements AtomValues interface. This sets the values to
* the evaluation tree EvalTree. */
void setValues(EvalTree& et) const;
/** This adds a constant with the given tree index. The
* constant must be checked previously and asserted that it
* does not exist. */
void add_constant(int t, double val);
/** Returns true if the tree index is either an hardwired
* constant (initial number OperationTree:num_constants in
* OperationTree) or the tree index is a registered constant
* by add_constant method. */
bool is_constant(int t) const;
double get_constant_value(int t) const;
/** Return -1 if the given string representation of a constant
* is not among the constants (double represenations). If it
* is, its tree index is returned. */
int check(const char* str) const;
/** Debug print. */
void print() const;
const Tconstantmap& get_constantmap() const
{return cmap;}
private:
/** Inverse map to Tconstantmap. */
typedef map<double,int> Tconstantinvmap;
/** This is an inverse map to cmap. This is only used for fast
* queries for the existing double constants in check
* method and add_constant. */
Tconstantinvmap cinvmap;
};
/** This class is a parent to Atoms classes which distinguish between
* constants (numerical literals), and variables with lags and
* leads. This abstraction does not distinguish between a parameter
* and a variable without lag or lead. In this sense, everything is a
* variable.*/
class DynamicAtoms : public Atoms, public Constants {
public:
/** Definition of a type mapping lags to the indices of the variables. */
typedef map<int,int> Tlagmap;
protected:
/** Definition of a type mapping names of the atoms to Tlagmap. */
typedef map<const char*, Tlagmap, ltstr> Tvarmap;
/** Definition of a type mapping indices of variables to the variable names. */
typedef map<int, const char*> Tindexmap;
/** This is just a storage for variable names, since all other
* instances of a variable name just point to the memory
* allocated by this object. */
NameStorage varnames;
/** This is the map for variables. Each variable name is
* mapped to the Tlagmap, which maps lags/leads to the nulary
* term indices in the tree. */
Tvarmap vars;
/** This is almost inverse map to the vars. It maps variable
* indices to the names. A returned name can be in turn used
* as a key in vars. */
Tindexmap indices;
/** Number of variables. */
int nv;
/** Minimum lag, if there is at least one lag, than this is a negative number. */
int minlag;
/** Maximum lead, if there is at least one lead, than this is a positive number. */
int maxlead;
public:
/** Construct empty DynamicAtoms. */
DynamicAtoms();
DynamicAtoms(const DynamicAtoms& da);
virtual ~DynamicAtoms() {}
/** Check the nulary term identified by its string
* representation. The nulary term can be either a constant or
* a variable. If constant, -1 is returned so that it could be
* assigned regardless if the same constant has already
* appeared or not. If variable, then -1 is returned only if
* the variable has not been assigned an index, otherwise the
* assigned index is returned. */
int check(const char* name) const;
/** Assign the nulary term identified by its string
* representation. This method should be called when check()
* returns -1. */
void assign(const char* name, int t);
/** Return a number of all variables. */
int nvar() const
{ return nv; }
/** Return the vector of variable indices. */
vector<int> variables() const;
/** Return max lead and min lag for a variable given by the
* index. If a variable cannot be found, the method retursn
* the smallest integer as maxlead and the largest integer as
* minlag. */
void varspan(int t, int& mlead, int& mlag) const;
/** Return max lead and min lag for a variable given by the
* name (without lead, lag). The same is valid if the variable
* name cannot be found. */
void varspan(const char* name, int& mlead, int& mlag) const;
/** Return max lead and min lag for a vector of variables given by the names. */
void varspan(const vector<const char*>& names, int& mlead, int& mlag) const;
/** Return true for all tree indices corresponding to a
* variable in the sense of this class. (This is parameters,
* exo and endo). Since the semantics of 'variable' will be
* changed in subclasses, we use name 'named atom'. These are
* all atoms but constants. */
bool is_named_atom(int t) const;
/** Return index of the variable described by the variable
* name and lag/lead. If it doesn't exist, return -1. */
int index(const char* name, int ll) const;
/** Return the lag map for the variable name. */
const Tlagmap& lagmap(const char* name) const;
/** Return the variable name for the tree index. It throws an
* exception if the tree index t is not a named atom. */
const char* name(int t) const;
/** Return the lead/lag for the tree index. It throws an
* exception if the tree index t is not a named atom. */
int lead(int t) const;
/** Return maximum lead. */
int get_maxlead() const
{return maxlead;}
/** Return minimum lag. */
int get_minlag() const
{return minlag;}
/** Return the name storage to allow querying to other
* classes. */
const NameStorage& get_name_storage() const
{return varnames;}
/** Assign the variable with a given lead. The varname must be
* from the varnames storage. The method checks if the
* variable iwht the given lead/lag is not assigned. If so, an
* exception is thrown. */
void assign_variable(const char* varname, int ll, int t);
/** Unassign the variable with a given lead and given tree
* index. The tree index is only provided as a check. An
* exception is thrown if the name, ll, and the tree index t
* are not consistent. The method also updates nv, indices,
* maxlead and minlag. The varname must be from the varnames
* storage. */
void unassign_variable(const char* varname, int ll, int t);
/** Debug print. */
void print() const;
protected:
/** Do the check for the variable. A subclass may need to
* reimplement this so that it could raise an error if the
* variable is not among a given list. */
virtual int check_variable(const char* name) const;
/** Assign the constant. */
void assign_constant(const char* name, int t);
/** Assign the variable. */
void assign_variable(const char* name, int t);
/** The method just updates minlag or/and maxlead. Note that
* when assigning variables, the update is done when inserting
* to the maps, however, if removing a variable, we need to
* call this method. */
void update_minmaxll();
/** The method parses the string to recover a variable name
* and lag/lead ll. The variable name doesn't contain a lead/lag. */
virtual void parse_variable(const char* in, string& out, int& ll) const = 0;
public:
/** Return true if the str represents a double.*/
static bool is_string_constant(const char* str);
};
/** This class is a parent of all orderings of the dynamic atoms
* of variables which can appear before t, at t, or after t. It
* encapsulates the ordering, and the information about the number
* of static (appearing only at time t) predetermined (appearing
* before t and possibly at t), both (appearing before t and after
* t and possibly at t) and forward looking (appearing after t and
* possibly at t).
*
* The constructor takes a list of variable names. The class also
* provides mapping from the ordering of the variables in the list
* (outer) to the new ordering (at time t) and back.
*
* The user of the subclass must call do_ordering() after
* initialization.
*
* The class contains a few preimplemented methods for
* ordering. The class is used in this way: Make a subclass, and
* implement pure virtual do_ordering() by just plugging a
* preimplemented method, or plugging your own implementation. The
* method do_ordering() is called by the user after the constructor.
*/
class VarOrdering {
protected:
/** Number of static variables. */
int n_stat;
/** Number of predetermined variables. */
int n_pred;
/** Number of both variables. */
int n_both;
/** Number of forward looking variables. */
int n_forw;
/** This is a set of tree indices corresponding to the
* variables at all times as they occur in the formulas. In
* fact, since this is used only for derivatives, the ordering
* of this vector is only important for ordering of the
* derivatives, in other contexts the ordering is not
* important, so it is rather a set of indices.*/
vector<int> der_atoms;
/** This maps tree index of the variable to the position in
* the row of the ordering. One should be careful with making
* space in the positions for variables not appearing at time
* t. For instance in the pred(t-1), both(t-1), stat(t),
* pred(t), both(t), forw(t), both(t+1), forw(t+1) ordering,
* the variables x(t-1), y(t-1), x(t+1), z(t-1), z(t), and
* z(t+1) having tree indices 6,5,4,3,2,1 will be ordered as
* follows: y(t-1), x(t-1), z(t-1), [y(t)], [x(t)], z(t),
* x(t+1), where a bracketed expresion means non-existent by
* occupying a space. The map thus will look as follows:
* {5->0, 6->1, 3->2, 2->5, 3->6}. Note that nothing is mapped
* to positions 3 and 4. */
map<int,int> positions;
/** This maps an ordering of the list of variables in
* constructor to the new ordering (at time t). The length is
* the number of variables. */
vector<int> outer2y;
/** This maps a new ordering to the ordering of the list of
* variables in constructor (at time t). The length is the
* number of variables. */
vector<int> y2outer;
/** This is just a reference for variable names to keep it
* from constructor to do_ordering() implementations. */
const vector<const char*>& varnames;
/** This is just a reference to atoms to keep it from
* constructor to do_ordering() implementations. */
const DynamicAtoms& atoms;
public:
/** This is an enum type for an ordering type implemented by
* do_general. */
enum ord_type {pbspbfbf, bfspbfpb};
/** Construct the ordering of the variables given by the names
* with their dynamic occurrences defined by the atoms. It
* calls the virtual method do_ordering which can be
* reimplemented. */
VarOrdering(const vector<const char*>& vnames, const DynamicAtoms& a)
: n_stat(0), n_pred(0), n_both(0), n_forw(0), varnames(vnames), atoms(a)
{}
VarOrdering(const VarOrdering& vo, const vector<const char*>& vnames,
const DynamicAtoms& a);
virtual VarOrdering* clone(const vector<const char*>& vnames,
const DynamicAtoms& a) const = 0;
/** Destructor does nothing here. */
virtual ~VarOrdering() {}
/** This is the method setting the ordering and the map. A
* subclass must reimplement it, possibly using a
* preimplemented ordering. This method must be called by the
* user after the class has been created. */
virtual void do_ordering() = 0;
/** Return number of static. */
int nstat() const
{return n_stat;}
/** Return number of predetermined. */
int npred() const
{return n_pred;}
/** Return number of both. */
int nboth() const
{return n_both;}
/** Return number of forward looking. */
int nforw() const
{return n_forw;}
/** Return the set of tree indices for derivatives. */
const vector<int>& get_der_atoms() const
{return der_atoms;}
/** Return the y2outer. */
const vector<int>& get_y2outer() const
{return y2outer;}
/** Return the outer2y. */
const vector<int>& get_outer2y() const
{return outer2y;}
/** Query the atom given by the tree index. True is returned
* if the atom is one of the variables in the object. */
bool check(int t) const;
/** Return the position of the atom (nulary term) given by a
* tree index. It is a lookup to the map. If the atom cannot
* be found, the exception is raised. */
int get_pos_of(int t) const;
/** This returns a length of ordered row of atoms. In all
* cases so far, it does not depend on the ordering and it is
* as follows. */
int length() const
{return n_stat+2*n_pred+3*n_both+2*n_forw;}
/** Debug print. */
void print() const;
protected:
/** This is a general ordering method which orders the
* variables by the given ordering ord_type. See documentation
* for respective do_ methods. */
void do_general(ord_type ordering);
/** This is a preimplemented ordering for do_ordering()
* method. It assumes that the variables appear only at time
* t-1, t, t+1. It orders the atoms as pred(t-1), both(t-1),
* stat(t), pred(t), both(t), forw(t), both(t+1),
* forw(t+1). It builds the der_atoms, the map of positions,
* as well as y2outer and outer2y. */
void do_pbspbfbf()
{do_general(pbspbfbf);}
/** This is a preimplemented ordering for do_ordering()
* method. It assumes that the variables appear only at time
* t-1, t, t+1. It orders the atoms as both(t+1), forw(t+1),
* stat(t), pred(t), both(t), forw(t), pred(t-1),
* both(t-1). It builds the der_atoms, the map of positions,
* as well as y2outer and outer2y. */
void do_bfspbfpb()
{do_general(bfspbfpb);}
/** This is a preimplemented ordering for do_ordering()
* method. It makes no assumptions about occurences of
* variables at different times. It orders the atoms with
* increasing time keeping the given ordering within one
* time. This implies that y2outer and outer2y will be
* identities. The der_atoms will be just a sequence of atoms
* from the least to the most time preserving the order of atoms
* within one time. */
void do_increasing_time();
private:
/** Declare this copy constructor as private to hide it. */
VarOrdering(const VarOrdering& vo);
};
};
#endif
// Local Variables:
// mode:C++
// End:
// Copyright (C) 2005, Ondra Kamenik
// $Id: fine_atoms.cpp 1759 2008-03-31 14:25:20Z kamenik $
#include "utils/cc/exception.h"
#include "parser_exception.h"
#include "fine_atoms.h"
using namespace ogp;
AllvarOuterOrdering::AllvarOuterOrdering(const vector<const char*>& allvar_outer,
const FineAtoms& a)
: atoms(a), allvar(),
endo2all(a.get_endovars().size(), -1),
exo2all(a.get_exovars().size(), -1)
{
// fill in the allvar from allvar_outer
for (unsigned int i = 0; i < allvar_outer.size(); i++) {
const char* s = atoms.varnames.query(allvar_outer[i]);
if (s)
allvar.push_back(s);
else
throw ogu::Exception(__FILE__, __LINE__,
string("Variable ") + allvar_outer[i] + " is not a declared symbol in AllvarOuterOrdering constructor");
}
// fill in endo2all and exo2all
for (unsigned int i = 0; i < allvar.size(); i++) {
Tvarintmap::const_iterator it = atoms.endo_outer_map.find(allvar[i]);
if (it != atoms.endo_outer_map.end())
endo2all[(*it).second] = i;
else {
it = atoms.exo_outer_map.find(allvar[i]);
if (it != atoms.exo_outer_map.end())
exo2all[(*it).second] = i;
else
throw ogu::Exception(__FILE__, __LINE__,
string("Name ") + allvar[i] + " is neither endogenous nor exogenous variable in AllvarOuterOrdering constructor");
}
}
// check whether everything has been filled
unsigned int iendo = 0;
while (iendo < endo2all.size() && endo2all[iendo] != -1) iendo++;
unsigned int iexo = 0;
while (iexo < exo2all.size() && exo2all[iexo] != -1) iexo++;
if (iendo < endo2all.size())
throw ogu::Exception(__FILE__, __LINE__,
string("Endogenous variable ") + atoms.get_endovars()[iendo] +
" not found in outer all ordering in AllvarOuterOrdering constructor");
if (iexo < exo2all.size())
throw ogu::Exception(__FILE__, __LINE__,
string("Exogenous variable ") + atoms.get_exovars()[iexo] +
" not found in outer all ordering in AllvarOuterOrdering constructor");
}
AllvarOuterOrdering::AllvarOuterOrdering(const AllvarOuterOrdering& avo,
const FineAtoms& a)
: atoms(a), allvar(),
endo2all(avo.endo2all),
exo2all(avo.exo2all)
{
// fill in the allvar from avo.allvar
for (unsigned int i = 0; i < avo.allvar.size(); i++) {
const char* s = atoms.varnames.query(avo.allvar[i]);
allvar.push_back(s);
}
}
FineAtoms::FineAtoms(const FineAtoms& fa)
: DynamicAtoms(fa), params(), endovars(), exovars(),
endo_order(NULL), exo_order(NULL), allvar_order(NULL),
der_atoms(fa.der_atoms),
endo_atoms_map(fa.endo_atoms_map),
exo_atoms_map(fa.exo_atoms_map)
{
// fill in params
for (unsigned int i = 0; i < fa.params.size(); i++) {
const char* s = varnames.query(fa.params[i]);
if (! s)
throw ogu::Exception(__FILE__, __LINE__,
string("Parameter ") + fa.params[i] + " does not exist in FineAtoms copy cosntructor");
params.push_back(s);
param_outer_map.insert(Tvarintmap::value_type(s, params.size()-1));
}
// fill in endovars
for (unsigned int i = 0; i < fa.endovars.size(); i++) {
const char* s = varnames.query(fa.endovars[i]);
if (! s)
throw ogu::Exception(__FILE__, __LINE__,
string("Endo variable ") + fa.endovars[i] + " does not exist in FineAtoms copy constructor");
endovars.push_back(s);
endo_outer_map.insert(Tvarintmap::value_type(s, endovars.size()-1));
}
// fill in exovars
for (unsigned int i = 0; i < fa.exovars.size(); i++) {
const char* s = varnames.query(fa.exovars[i]);
if (! s)
throw ogu::Exception(__FILE__, __LINE__,
string("Exo variable ") + fa.exovars[i] + " does not exist in FineAtoms copy cosntructor");
exovars.push_back(s);
exo_outer_map.insert(Tvarintmap::value_type(s, exovars.size()-1));
}
if (fa.endo_order)
endo_order = fa.endo_order->clone(endovars, *this);
if (fa.exo_order)
exo_order = fa.exo_order->clone(exovars, *this);
if (fa.allvar_order)
allvar_order = new AllvarOuterOrdering(*(fa.allvar_order), *this);
}
int FineAtoms::check_variable(const char* name) const
{
string str;
int ll;
parse_variable(name, str, ll);
if (varnames.query(str.c_str()))
return DynamicAtoms::check_variable(name);
else {
throw ParserException(string("Variable <")+str+"> not declared.",0);
return -1;
}
}
int FineAtoms::num_exo_periods() const
{
int mlead, mlag;
exovarspan(mlead, mlag);
return mlead-mlag+1;
}
void FineAtoms::parsing_finished(VarOrdering::ord_type ot)
{
make_internal_orderings(ot);
// by default, concatenate outer endo and outer exo and make it as
// allvar outer:
vector<const char*> allvar_tmp;
allvar_tmp.insert(allvar_tmp.end(), endovars.begin(), endovars.end());
allvar_tmp.insert(allvar_tmp.end(), exovars.begin(), exovars.end());
if (allvar_order)
delete allvar_order;
allvar_order = new AllvarOuterOrdering(allvar_tmp, *this);
}
void FineAtoms::parsing_finished(VarOrdering::ord_type ot,
const vector<const char*> allvar)
{
make_internal_orderings(ot);
if (allvar_order)
delete allvar_order;
allvar_order = new AllvarOuterOrdering(allvar, *this);
}
const vector<const char*>& FineAtoms::get_allvar() const
{
if (! allvar_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::get_allvars called before parsing_finished");
return allvar_order->get_allvar();
}
const vector<int>& FineAtoms::outer_endo2all() const
{
if (! allvar_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::outer_endo2all called before parsing_finished");
return allvar_order->get_endo2all();
}
const vector<int>& FineAtoms::outer_exo2all() const
{
if (! allvar_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::outer_exo2all called before parsing_finished");
return allvar_order->get_exo2all();
}
vector<int> FineAtoms::variables() const
{
if (endo_order) {
return der_atoms;
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::variables called before parsing_finished");
return vector<int>();
}
}
int FineAtoms::nstat() const
{
if (endo_order) {
return endo_order->nstat();
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::nstat called before parsing_finished");
return -1;
}
}
int FineAtoms::npred() const
{
if (endo_order) {
return endo_order->npred();
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::npred called before parsing_finished");
return -1;
}
}
int FineAtoms::nboth() const
{
if (endo_order) {
return endo_order->nboth();
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::nboth called before parsing_finished");
return -1;
}
}
int FineAtoms::nforw() const
{
if (endo_order) {
return endo_order->nforw();
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::nforw called before parsing_finished");
return -1;
}
}
int FineAtoms::get_pos_of_endo(int t) const
{
if (endo_order) {
return endo_order->get_pos_of(t);
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::get_pos_of_endo called before parsing_finished");
return -1;
}
}
int FineAtoms::get_pos_of_exo(int t) const
{
if (exo_order) {
return exo_order->get_pos_of(t);
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::get_pos_of_exo called before parsing_finished");
return -1;
}
}
int FineAtoms::get_pos_of_all(int t) const
{
if (endo_order && exo_order) {
if (endo_order->check(t))
return endo_order->get_pos_of(t);
else if (exo_order->check(t))
return endo_order->length() + exo_order->get_pos_of(t);
else {
throw ogu::Exception(__FILE__,__LINE__,
"Atom is not endo nor exo in FineAtoms::get_pos_of_all");
return -1;
}
} else {
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::get_pos_of_exo called before parsing_finished");
return -1;
}
}
const vector<int>& FineAtoms::y2outer_endo() const
{
if (! endo_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::y2outer_endo called before parsing_finished");
return endo_order->get_y2outer();
}
const vector<int>& FineAtoms::outer2y_endo() const
{
if (! endo_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::outer2y_endo called before parsing_finished");
return endo_order->get_outer2y();
}
const vector<int>& FineAtoms::y2outer_exo() const
{
if (! exo_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::y2outer_endo called before parsing_finished");
return exo_order->get_y2outer();
}
const vector<int>& FineAtoms::outer2y_exo() const
{
if (! exo_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::outer2y_exo called before parsing_finished");
return exo_order->get_outer2y();
}
const vector<int>& FineAtoms::get_endo_atoms_map() const
{
if (! endo_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::get_endo_atoms_map called before parsing_finished");
return endo_atoms_map;
}
const vector<int>& FineAtoms::get_exo_atoms_map() const
{
if (! exo_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::get_exo_atoms_map called before parsing_finished");
return exo_atoms_map;
}
int FineAtoms::name2outer_param(const char* name) const
{
Tvarintmap::const_iterator it = param_outer_map.find(name);
if (it == param_outer_map.end())
throw ogu::Exception(__FILE__,__LINE__,
"Name is not a parameter in FineAtoms::name2outer_param");
return (*it).second;
}
int FineAtoms::name2outer_endo(const char* name) const
{
Tvarintmap::const_iterator it = endo_outer_map.find(name);
if (it == endo_outer_map.end())
throw ogu::Exception(__FILE__,__LINE__,
"Name is not an endogenous variable in FineAtoms::name2outer_endo");
return (*it).second;
}
int FineAtoms::name2outer_exo(const char* name) const
{
Tvarintmap::const_iterator it = exo_outer_map.find(name);
if (it == exo_outer_map.end())
throw ogu::Exception(__FILE__,__LINE__,
"Name is not an exogenous variable in FineAtoms::name2outer_exo");
return (*it).second;
}
int FineAtoms::name2outer_allvar(const char* name) const
{
if (! allvar_order)
throw ogu::Exception(__FILE__,__LINE__,
"FineAtoms::name2outer_allvar called beore parsing_finished");
Tvarintmap::const_iterator it = endo_outer_map.find(name);
if (it != endo_outer_map.end())
return allvar_order->get_endo2all()[(*it).second];
else {
it = exo_outer_map.find(name);
if (it != exo_outer_map.end())
return allvar_order->get_exo2all()[(*it).second];
}
throw ogu::Exception(__FILE__,__LINE__,
string("Name ") + name + " is neither endo nor exo variable in FineAtoms::name2outer_allvar");
return -1;
}
void FineAtoms::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);
endovars.push_back(ss);
endo_outer_map.insert(Tvarintmap::value_type(ss, endovars.size()-1));
}
void FineAtoms::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);
exovars.push_back(ss);
exo_outer_map.insert(Tvarintmap::value_type(ss, exovars.size()-1));
}
void FineAtoms::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);
params.push_back(ss);
param_outer_map.insert(Tvarintmap::value_type(ss, params.size()-1));
}
void FineAtoms::make_internal_orderings(VarOrdering::ord_type ot)
{
bool endo_ordering_done = false;
bool exo_ordering_done = false;
order_type = ot;
int mlead, mlag;
endovarspan(mlead, mlag);
if (mlag >= -1 && mlead <= 1) {
// make endo ordering
if (endo_order)
delete endo_order;
if (ot == VarOrdering::pbspbfbf)
endo_order = new EndoVarOrdering1(endovars, *this);
else
endo_order = new EndoVarOrdering2(endovars, *this);
endo_order->do_ordering();
endo_ordering_done = true;
}
exovarspan(mlead, mlag);
if (mlag == 0 && mlead == 0) {
// make exo ordering
if (exo_order)
delete exo_order;
exo_order = new ExoVarOrdering(exovars, *this);
exo_order->do_ordering();
exo_ordering_done = true;
}
if (endo_ordering_done && exo_ordering_done) {
// concatenate der atoms from endo_order and exo_order
der_atoms.clear();
der_atoms.insert(der_atoms.end(),
endo_order->get_der_atoms().begin(),
endo_order->get_der_atoms().end());
der_atoms.insert(der_atoms.end(),
exo_order->get_der_atoms().begin(),
exo_order->get_der_atoms().end());
// create endo_atoms_map; der_atoms is a concatenation, so it is easy
int endo_atoms = endo_order->get_der_atoms().size();
endo_atoms_map.clear();
for (int i = 0; i < endo_atoms; i++)
endo_atoms_map.push_back(i);
// create exo_atoms_map
int exo_atoms = exo_order->get_der_atoms().size();
exo_atoms_map.clear();
for (int i = 0; i < exo_atoms; i++)
exo_atoms_map.push_back(endo_atoms + i);
}
}
void FineAtoms::print() const
{
DynamicAtoms::print();
if (endo_order) {
printf("Endo ordering:\n");
endo_order->print();
} else {
printf("Endo ordering not created.\n");
}
if (exo_order) {
printf("Exo ordering:\n");
exo_order->print();
} else {
printf("Exo ordering not created.\n");
}
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]);
}
// Copyright (C) 2005, Ondra Kamenik
// $Id: fine_atoms.h 1759 2008-03-31 14:25:20Z kamenik $
#ifndef OGP_FINE_ATOMS_H
#define OGP_FINE_ATOMS_H
#include "dynamic_atoms.h"
#include <vector>
#include <string>
namespace ogp {
using std::vector;
using std::string;
/** This is just ordering used for endogenous variables. It
* assumes that we have only time t-1, t, and t+1, orders them as
* pred(t-1), both(t-1), stat(t), pred(t), both(t), forw(t),
* both(t+1), forw(t+1). */
class EndoVarOrdering1 : public VarOrdering {
public:
EndoVarOrdering1(const vector<const char*>& vnames, const DynamicAtoms& a)
: VarOrdering(vnames, a) {}
EndoVarOrdering1(const EndoVarOrdering1& vo, const vector<const char*>& vnames,
const DynamicAtoms& a)
: VarOrdering(vo, vnames, a) {}
VarOrdering* clone(const vector<const char*>& vnames, const DynamicAtoms& a) const
{return new EndoVarOrdering1(*this, vnames, a);}
void do_ordering()
{do_pbspbfbf();}
};
/** This is just another ordering used for endogenous
* variables. It assumes that we have only time t-1, t, and t+1,
* orders them as both(t+1), forw(t+1), pred(t-1), both(t-1),
* stat(t), pred(t), both(t), forw(t). */
class EndoVarOrdering2 : public VarOrdering {
public:
EndoVarOrdering2(const vector<const char*>& vnames, const DynamicAtoms& a)
: VarOrdering(vnames, a) {}
EndoVarOrdering2(const EndoVarOrdering2& vo, const vector<const char*>& vnames,
const DynamicAtoms& a)
: VarOrdering(vo, vnames, a) {}
VarOrdering* clone(const vector<const char*>& vnames, const DynamicAtoms& a) const
{return new EndoVarOrdering2(*this, vnames, a);}
void do_ordering()
{do_bfspbfpb();}
};
/** This is just ordering used for exogenous variables. It makes
* no assumptions about their timing. It orders them from the
* least time to the latest time. */
class ExoVarOrdering : public VarOrdering {
public:
ExoVarOrdering(const vector<const char*>& vnames, const DynamicAtoms& a)
: VarOrdering(vnames, a) {}
ExoVarOrdering(const ExoVarOrdering& vo, const vector<const char*>& vnames,
const DynamicAtoms& a)
: VarOrdering(vo, vnames, a) {}
VarOrdering* clone(const vector<const char*>& vnames, const DynamicAtoms& a) const
{return new ExoVarOrdering(*this, vnames, a);}
void do_ordering()
{do_increasing_time();}
};
class FineAtoms;
/** This class provides an outer ordering of all variables (endo
* and exo). It maps the ordering to the particular outer
* orderings of endo and exo. It works tightly with the FineAtoms
* class. */
class AllvarOuterOrdering {
protected:
/** Type for a map mapping a variable name to an integer. */
typedef map<const char*, int, ltstr> Tvarintmap;
/** Reference to atoms. */
const FineAtoms& atoms;
/** The vector of all endo and exo variables in outer
* ordering. The pointers point to storage in atoms. */
vector<const char*> allvar;
/** The mapping from outer endogenous to outer all. For
* example endo2all[0] is the order of the first outer
* endogenous variable in the allvar ordering. */
vector<int> endo2all;
/** The mapping from outer exogenous to outer all. For example
* exo2all[0] is the order of the first outer exogenous
* variables in the allvar ordering. */
vector<int> exo2all;
public:
/** Construct the allvar outer ordering from the provided
* sequence of endo and exo names. The names can have an
* arbitrary storage, the storage is transformed to the atoms
* storage. An exception is thrown if either the list is not
* exhaustive, or some string is not a variable. */
AllvarOuterOrdering(const vector<const char*>& allvar_outer, const FineAtoms& a);
/** Copy constructor using the storage of provided atoms. */
AllvarOuterOrdering(const AllvarOuterOrdering& allvar_outer, const FineAtoms& a);
/** Return endo2all mapping. */
const vector<int>& get_endo2all() const
{return endo2all;}
/** Return exo2all mapping. */
const vector<int>& get_exo2all() const
{return exo2all;}
/** Return the allvar ordering. */
const vector<const char*>& get_allvar() const
{return allvar;}
};
/** This class refines the DynamicAtoms by distinguishing among
* parameters (no lag and leads) and endogenous and exogenous
* variables (with lags and leads). For parameters, endogenous and
* exogenous, it defines outer orderings and internal
* orderings. The internal orderings are created by
* parsing_finished() method when it is sure that no new variables
* would be registered. The outer orderings are given by the order
* of calls of registering methods.
*
* In addition, the class also defines outer ordering of
* endogenous and exogenous variables. This is input as a
* parameter to parsing_finished(). By default, this whole outer
* ordering is just a concatenation of outer ordering of
* endogenous and exogenous variables.
*
* The internal ordering of all endo and exo variables is just a
* concatenation of endo and exo variables in their internal
* orderings. This is the ordering with respect to which all
* derivatives are taken. */
class FineAtoms : public DynamicAtoms {
friend class AllvarOuterOrdering;
protected:
typedef map<const char*, int, ltstr> Tvarintmap;
private:
/** The vector of parameters names. The order gives the order
* the data is communicated with outside world. */
vector<const char*> params;
/** A map mapping a name of a parameter to an index in the outer
* 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 outer 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;
protected:
/** This is the internal ordering of all atoms corresponding
* to endogenous variables. It is constructed by
* parsing_finished() method, which should be called after all
* parsing jobs have been finished. */
VarOrdering* endo_order;
/** This is the internal ordering of all atoms corresponding
* to exogenous variables. It has the same handling as
* endo_order. */
VarOrdering* exo_order;
/** This is the all variables outer ordering. It is
* constructed by parsing finished. */
AllvarOuterOrdering* allvar_order;
/** This vector defines a set of atoms as tree indices used
* for differentiation. The order of the atoms in this vector
* defines ordering of the derivative tensors. The ordering is
* a concatenation of atoms from endo_order and then
* exo_order. 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:
FineAtoms()
: endo_order(NULL), exo_order(NULL), allvar_order(NULL) {}
FineAtoms(const FineAtoms& fa);
/** Deletes endo_order and exo_order. */
virtual ~FineAtoms()
{
if (endo_order) delete endo_order;
if (exo_order) delete exo_order;
if (allvar_order) delete allvar_order;
}
/** Overrides DynamicAtoms::check_variable so that the error
* would be raised if the variable name is not declared. A
* variable is declared by inserting it to
* DynamicAtoms::varnames. This is a responsibility of a
* subclass. */
int check_variable(const char* name) const;
/** This calculates min lag and max lead of endogenous variables. */
void endovarspan(int& mlead, int& mlag) const
{varspan(endovars, mlead, mlag);}
/** This calculates mim lag and max lead of exogenous variables. */
void exovarspan(int& mlead, int& mlag) const
{varspan(exovars, mlead, mlag);}
/** This calculates the number of periods in which at least
* one exogenous variable occurs. */
int num_exo_periods() const;
/** Return an (external) ordering of parameters. */
const vector<const char*>& get_params() const
{return params;}
/** Return an external ordering of endogenous variables. */
const vector<const char*>& get_endovars() const
{return endovars;}
/** Return an external ordering of exogenous variables. */
const vector<const char*>& get_exovars() const
{return exovars;}
/** This constructs internal orderings and makes the indices
* returned by variables method available. Further it
* constructs outer ordering of all variables by a simple
* concatenation of outer endogenous and outer exogenous. In
* addition, it makes nstat, npred, nboth, nforw available. */
void parsing_finished(VarOrdering::ord_type ot);
/** This does the same thing as
* parsing_finished(VarOrdering::ord_type) plus it allows for
* inputing a different outer ordering of all variables. The
* ordering is input as a list of strings, their storage can
* be arbitrary. */
void parsing_finished(VarOrdering::ord_type ot, const vector<const char*> avo);
/** Return the external ordering of all variables (endo and
* exo). This is either the second argument to
* parsing_finished or the default external ordering. This
* must be called only after parsing_finished. */
const vector<const char*>& get_allvar() const;
/** Return the map from outer ordering of endo variables to
* the allvar ordering. This must be called only after
* parsing_finished. */
const vector<int>& outer_endo2all() const;
/** Return the map from outer ordering of exo variables to
* the allvar ordering. This must be called only after
* parsing_finished. */
const vector<int>& outer_exo2all() const;
/** Return the atoms with respect to which we are going to
* differentiate. This must be called after
* parsing_finished. */
vector<int> variables() const;
/** Return the number of static. */
int nstat() const;
/** Return the number of predetermined. */
int npred() const;
/** Return the number of both. */
int nboth() const;
/** Return the number of forward looking. */
int nforw() const;
/** Return the index of an endogenous atom given by tree index in
* the endo ordering. This must be also called only after
* parsing_finished(). */
int get_pos_of_endo(int t) const;
/** Return the index of an exogenous atom given by tree index in
* the exo ordering. This must be also called only after
* parsing_finished(). */
int get_pos_of_exo(int t) const;
/** Return the index of either endogenous or exogenous atom
* given by tree index in the concatenated ordering of
* endogenous and exogenous atoms. This must be also called
* only after parsing_finished(). */
int get_pos_of_all(int t) const;
/** Return the mapping from endogenous at time t to outer
* ordering of endogenous. */
const vector<int>& y2outer_endo() const;
/** Return the mapping from the outer ordering of endogenous to endogenous
* at time t. */
const vector<int>& outer2y_endo() const;
/** Return the mapping from exogenous at time t to outer
* ordering of exogenous. */
const vector<int>& y2outer_exo() const;
/** Return the mapping from the outer ordering of exogenous to exogenous
* at time t. */
const vector<int>& outer2y_exo() const;
/** Return the endo_atoms_map. */
const vector<int>& get_endo_atoms_map() const;
/** Return the exo_atoms_map. */
const vector<int>& get_exo_atoms_map() const;
/** 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 an index in the outer ordering of all variables
* (endo and exo) for a given name. An exception is thrown if
* the name is not a variable. This must be called only after
* parsing_finished(). */
int name2outer_allvar(const char* name) const;
/** Return the number of endogenous variables at time t-1, these are state
* variables. */
int nys() const
{return npred()+nboth();}
/** Return the number of endogenous variables at time t+1. */
int nyss() const
{return nboth()+nforw();}
/** Return the number of endogenous variables. */
int ny() const
{return endovars.size();}
/** Return the number of exogenous variables. */
int nexo() const
{return (int)exovars.size();}
/** Return the number of parameters. */
int np() const
{return (int)(params.size());}
/** Register unique endogenous variable name. The order of
* calls defines the endo outer ordering. The method is
* virtual, since a superclass may want to do some additional
* action. */
virtual void register_uniq_endo(const char* name);
/** Register unique exogenous variable name. The order of
* calls defines the exo outer ordering. The method is
* virtual, since a superclass may want to do somem additional
* action. */
virtual void register_uniq_exo(const char* name);
/** Register unique parameter name. The order of calls defines
* the param outer ordering. The method is
* virtual, since a superclass may want to do somem additional
* action. */
virtual void register_uniq_param(const char* name);
/** Debug print. */
void print() const;
private:
/** This performs the common part of parsing_finished(), which
* is a construction of internal orderings. */
void make_internal_orderings(VarOrdering::ord_type ot);
protected:
/** This remembers the ordering type of the last call make_internal_ordering. */
VarOrdering::ord_type order_type;
};
};
#endif
// Local Variables:
// mode:C++
// End:
%{
#include "location.h"
#include "formula_tab.hh"
extern YYLTYPE fmla_lloc;
#define YY_USER_ACTION SET_LLOC(fmla_);
%}
%option nounput
%option noyy_top_state
%option stack
%option yylineno
%option prefix="fmla_"
%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]
[+] {return YPLUS;}
[-] {return YMINUS;}
[*] {return YTIMES;}
[/] {return YDIVIDE;}
[\^] {return YPOWER;}
exp {return YEXP;}
log {return YLOG;}
sin {return YSIN;}
cos {return YCOS;}
tan {return YTAN;}
sqrt {return YSQRT;}
erf {return YERF;}
erfc {return YERFC;}
diff {return YDIFF;}
/* names: parameters, variables (lagged/leaded) */
[A-Za-z_][A-Za-z0-9_]*([\(\{][+-]?[0-9]+[\)\}])? {
fmla_lval.string=fmla_text;
return NAME;
}
/* floating point numbers */
(([0-9]*\.?[0-9]+)|([0-9]+\.))([edED][-+]?[0-9]+)? {
fmla_lval.string=fmla_text;
return DNUMBER;
}
= {return EQUAL_SIGN;}
. {return fmla_text[0];}
%%
int fmla_wrap()
{
return 1;
}
void fmla__destroy_buffer(void* p)
{
fmla__delete_buffer((YY_BUFFER_STATE)p);
}
%{
/* Copyright (C) 2006-2011, Ondra Kamenik */
#include <cstdio>
#include "location.h"
#include "formula_parser.h"
#include "formula_tab.hh"
void fmla_error(const char*);
int fmla_lex(void);
extern int fmla_lineno;
extern ogp::FormulaParser* fparser;
extern YYLTYPE fmla_lloc;
// static void print_token_value (FILE *, int, YYSTYPE);
// #define YYPRINT(file, type, value) print_token_value (file, type, value)
%}
%union {
char* string;
double dvalue;
int integer;
}
%token EQUAL_SIGN
%left YPLUS YMINUS
%left YTIMES YDIVIDE
%left YUMINUS YUPLUS
%right YPOWER
%token YEXP YLOG YSIN YCOS YTAN YSQRT YERF YERFC YDIFF
%token <string> DNUMBER NAME
%type <integer> expression
%name-prefix="fmla_"
%locations
%error-verbose
%%
root : equation_list
| expression
{fparser->add_formula($1);}
;
equation_list : equation_list equation | equation ;
equation : expression EQUAL_SIGN expression ';'
{fparser->add_formula(fparser->add_binary(ogp::MINUS,$1,$3));}
| expression ';'
{fparser->add_formula($1);}
;
expression : '(' expression ')' { $$ = $2;}
| expression YPLUS expression {$$=fparser->add_binary(ogp::PLUS,$1,$3);}
| expression YMINUS expression {$$=fparser->add_binary(ogp::MINUS,$1,$3);}
| expression YTIMES expression {$$=fparser->add_binary(ogp::TIMES,$1,$3);}
| expression YDIVIDE expression {$$=fparser->add_binary(ogp::DIVIDE,$1,$3);}
| expression YPOWER expression {$$=fparser->add_binary(ogp::POWER,$1,$3);}
| YMINUS expression %prec YUMINUS {$$=fparser->add_unary(ogp::UMINUS,$2);}
| YPLUS expression %prec YUPLUS {$$ = $2;}
| YSIN '(' expression ')' {$$=fparser->add_unary(ogp::SIN,$3);}
| YCOS '(' expression ')' {$$=fparser->add_unary(ogp::COS,$3);}
| YTAN '(' expression ')' {$$=fparser->add_unary(ogp::TAN,$3);}
| YEXP '(' expression ')' {$$=fparser->add_unary(ogp::EXP,$3);}
| YLOG '(' expression ')' {$$=fparser->add_unary(ogp::LOG,$3);}
| YSQRT '(' expression ')' {$$=fparser->add_unary(ogp::SQRT,$3);}
| YERF '(' expression ')' {$$=fparser->add_unary(ogp::ERF,$3);}
| YERFC '(' expression ')' {$$=fparser->add_unary(ogp::ERFC,$3);}
| YDIFF '(' expression ',' NAME ')' {$$=fparser->add_derivative($3, fparser->add_nulary($5));}
| NAME {$$=fparser->add_nulary($1);}
| DNUMBER {$$=fparser->add_nulary($1);}
;
%%
void fmla_error(const char* s)
{
fparser->error(s);
}
/*
static void print_token_value(FILE* file, int type, YYSTYPE value)
{
if (type == NAME)
fprintf(file, "%s", value.string);
}
*/
// Copyright (C) 2005, Ondra Kamenik
// $Id: formula_parser.cpp 2268 2008-11-22 10:38:03Z michel $
#include "utils/cc/pascal_triangle.h"
#include "utils/cc/exception.h"
#include "parser_exception.h"
#include "location.h"
#include "formula_parser.h"
#include "formula_tab.hh"
#include <cmath>
using namespace ogp;
extern location_type fmla_lloc;
FormulaParser::FormulaParser(const FormulaParser& fp, Atoms& a)
: otree(fp.otree), atoms(a), formulas(fp.formulas), ders()
{
// create derivatives
for (unsigned int i = 0; i < fp.ders.size(); i++)
ders.push_back(new FormulaDerivatives(*(fp.ders[i])));
}
FormulaParser::~FormulaParser()
{
destroy_derivatives();
}
void FormulaParser::differentiate(int max_order)
{
destroy_derivatives();
vector<int> vars;
vars = atoms.variables();
for (unsigned int i = 0; i < formulas.size(); i++)
ders.push_back(new FormulaDerivatives(otree, vars, formulas[i], max_order));
}
const FormulaDerivatives& FormulaParser::derivatives(int i) const
{
if (i < (int)ders.size())
return *(ders[i]);
else
throw ogu::Exception(__FILE__,__LINE__,
"Wrong formula index in FormulaParser::derivatives");
return *(ders[0]); // just because of compiler
}
void FormulaParser::add_formula(int t)
{
formulas.push_back(t);
}
int FormulaParser::add_binary(code_t code, int t1, int t2)
{
return otree.add_binary(code, t1, t2);
}
int FormulaParser::add_unary(code_t code, int t)
{
return otree.add_unary(code, t);
}
int FormulaParser::add_nulary(const char* str)
{
int t = -1;
try {
t = atoms.check(str);
} catch (const ParserException& e) {
throw ParserException(e, fmla_lloc.off);
}
if (t == -1) {
t = otree.add_nulary();
atoms.assign(str, t);
}
return t;
}
void FormulaParser::add_subst_formulas(const map<int,int>& subst, const FormulaParser& fp)
{
for (int i = 0; i < fp.nformulas(); i++) {
int f = add_substitution(fp.formula(i), subst, fp);
add_formula(f);
}
}
void FormulaParser::substitute_formulas(const map<int,int>& smap)
{
for (int i = 0; i < nformulas(); i++) {
// make substitution and replace the formula for it
int f = add_substitution(formulas[i], smap);
formulas[i] = f;
// update the derivatives if any
if (i < (int)ders.size() && ders[i]) {
int order = ders[i]->get_order();
delete ders[i];
ders[i] = new FormulaDerivatives(otree, atoms.variables(), formulas[i], order);
}
}
}
/** Global symbols for passing info to parser. */
FormulaParser* fparser;
/** The declarations of functions defined in formula_ll.cc and
* formula_tab.cc generated from formula.lex and formula.y */
void* fmla__scan_buffer(char*, size_t);
void fmla__destroy_buffer(void*);
int fmla_parse();
extern location_type fmla_lloc;
/** This makes own copy of provided data, sets the buffer for the
* parser with fmla_scan_buffer, and launches fmla_parse(). Note that
* the pointer returned from fmla_scan_buffer must be freed at the
* end. */
void FormulaParser::parse(int length, const char* stream)
{
char* buffer = new char[length+2];
strncpy(buffer, stream, length);
buffer[length] = '\0';
buffer[length+1] = '\0';
fmla_lloc.off = 0;
fmla_lloc.ll = 0;
void* p = fmla__scan_buffer(buffer, (unsigned int)length+2);
fparser = this;
fmla_parse();
delete [] buffer;
fmla__destroy_buffer(p);
}
void FormulaParser::error(const char* mes) const
{
throw ParserException(mes, fmla_lloc.off);
}
int FormulaParser::last_formula() const
{
int res = -1;
for (unsigned int i = 0; i < formulas.size(); i++)
if (res < formulas[i])
res = formulas[i];
return std::max(res, otree.get_last_nulary());
}
int FormulaParser::pop_last_formula()
{
if (formulas.size() == 0)
return -1;
int t = formulas.back();
if (formulas.size() == ders.size()) {
delete ders.back();
ders.pop_back();
}
formulas.pop_back();
return t;
}
void FormulaParser::print() const
{
atoms.print();
for (unsigned int i = 0; i < formulas.size(); i++) {
printf("formula %d:\n", formulas[i]);
otree.print_operation(formulas[i]);
}
for (unsigned int i = 0; i < ders.size(); i++) {
printf("derivatives for the formula %d:\n", formulas[i]);
ders[i]->print(otree);
}
}
void FormulaParser::destroy_derivatives()
{
while (ders.size() > 0) {
delete ders.back();
ders.pop_back();
}
}
/** This constructor makes a vector of indices for formulas
* corresponding to derivatives of the given formula. The formula is
* supposed to belong to the provided tree, the created derivatives
* are added to the tree.
*
* The algorithm is as follows. todo: update description of the
* algorithm
*/
FormulaDerivatives::FormulaDerivatives(OperationTree& otree,
const vector<int>& vars, int f, int max_order)
: nvar(vars.size()), order(max_order)
{
FoldMultiIndex fmi_zero(nvar);
tder.push_back(f);
indices.push_back(fmi_zero);
unsigned int last_order_beg = 0;
unsigned int last_order_end = tder.size();
for (int k = 1; k <= order; k++) {
// interval <last_order_beg,last_order_end) is guaranteed
// here to contain at least one item
for (unsigned int run = last_order_beg; run < last_order_end; run++) {
// shift one order from the run
FoldMultiIndex fmi(indices[run], 1);
// set starting variable from the run, note that if k=1,
// the shift order ctor of fmi will set it to zero
int ivar_start = fmi[k-1];
for (int ivar = ivar_start; ivar < nvar; ivar++, fmi.increment()) {
int der = otree.add_derivative(tder[run], vars[ivar]);
if (der != OperationTree::zero) {
tder.push_back(der);
indices.push_back(fmi);
}
}
}
// set new last_order_beg and last_order_end
last_order_beg = last_order_end;
last_order_end = tder.size();
// if there was no new derivative, break out from the loop
if (last_order_beg >= last_order_end)
break;
}
// build ind2der map
for (unsigned int i = 0; i < indices.size(); i++)
ind2der.insert(Tfmiintmap::value_type(indices[i], i));
}
FormulaDerivatives::FormulaDerivatives(const FormulaDerivatives& fd)
: tder(fd.tder), indices(fd.indices), ind2der(fd.ind2der),
nvar(fd.nvar), order(fd.order)
{
}
int FormulaDerivatives::derivative(const FoldMultiIndex& mi) const
{
if (mi.order() > order)
throw ogu::Exception(__FILE__,__LINE__,
"Wrong order of multi-index in FormulaDerivatives::derivative");
if (mi.nv() != nvar)
throw ogu::Exception(__FILE__,__LINE__,
"Wrong multi-index variables in FormulaDerivatives::derivative");
Tfmiintmap::const_iterator it = ind2der.find(mi);
if (it == ind2der.end())
return OperationTree::zero;
else
return tder[(*it).second];
}
void FormulaDerivatives::print(const OperationTree& otree) const
{
for (Tfmiintmap::const_iterator it = ind2der.begin();
it != ind2der.end(); ++it) {
printf("derivative ");
(*it).first.print();
printf(" is formula %d\n", tder[(*it).second]);
otree.print_operation(tder[(*it).second]);
}
}
void FormulaCustomEvaluator::eval(const AtomValues& av, FormulaEvalLoader& loader)
{
etree.reset_all();
av.setValues(etree);
for (unsigned int i = 0; i < terms.size(); i++) {
double res = etree.eval(terms[i]);
loader.load((int)i, res);
}
}
FoldMultiIndex::FoldMultiIndex(int nv)
: nvar(nv), ord(0), data(new int[ord])
{
}
FoldMultiIndex::FoldMultiIndex(int nv, int ordd, int ii)
: nvar(nv), ord(ordd), data(new int[ord])
{
for (int i = 0; i < ord; i++)
data[i] = ii;
}
/** Note that a monotone sequence mapped by monotone mapping yields a
* monotone sequence. */
FoldMultiIndex::FoldMultiIndex(int nv, const FoldMultiIndex& mi, const vector<int>& mp)
: nvar(nv), ord(mi.ord), data(new int[ord])
{
for (int i = 0; i < ord; i++) {
if (i < ord-1 && mp[i+1] < mp[i])
throw ogu::Exception(__FILE__,__LINE__,
"Mapping not monotone in FoldMultiIndex constructor");
if (mp[mi[i]] >= nv || mp[mi[i]] < 0)
throw ogu::Exception(__FILE__,__LINE__,
"Mapping out of bounds in FoldMultiIndex constructor");
data[i] = mp[mi[i]];
}
}
FoldMultiIndex::FoldMultiIndex(const FoldMultiIndex& fmi, int new_orders)
: nvar(fmi.nvar),
ord(fmi.ord+new_orders),
data(new int[ord])
{
memcpy(data, fmi.data, fmi.ord*sizeof(int));
int new_item = (fmi.ord > 0)? fmi.data[fmi.ord-1] : 0;
for (int i = fmi.ord; i < ord; i++) {
data[i] = new_item;
}
}
FoldMultiIndex::FoldMultiIndex(const FoldMultiIndex& fmi)
: nvar(fmi.nvar),
ord(fmi.ord),
data(new int[fmi.ord])
{
memcpy(data, fmi.data, ord*sizeof(int));
}
const FoldMultiIndex& FoldMultiIndex::operator=(const FoldMultiIndex& fmi)
{
if (ord != fmi.ord) {
delete [] data;
data = new int[fmi.ord];
}
ord = fmi.ord;
nvar = fmi.nvar;
memcpy(data, fmi.data, ord*sizeof(int));
return *this;
}
bool FoldMultiIndex::operator<(const FoldMultiIndex& fmi) const
{
if (nvar != fmi.nvar)
ogu::Exception(__FILE__,__LINE__,
"Different nvar in FoldMultiIndex::operator<");
if (ord < fmi.ord)
return true;
if (ord > fmi.ord)
return false;
int i = 0;
while (i < ord && data[i] == fmi.data[i])
i++;
if (i == ord)
return false;
else
return data[i] < fmi.data[i];
}
bool FoldMultiIndex::operator==(const FoldMultiIndex& fmi) const
{
bool res = true;
res = res && (nvar == fmi.nvar) && (ord == fmi.ord);
if (res)
for (int i = 0; i < ord; i++)
if (data[i] != fmi.data[i])
return false;
return res;
}
void FoldMultiIndex::increment()
{
if (ord == 0)
return;
int k = ord-1;
data[k]++;
while (k > 0 && data[k] == nvar) {
data[k] = 0;
data[--k]++;
}
for (int kk = 1; kk < ord; kk++)
if (data[kk-1] > data[kk])
data[kk] = data[kk-1];
}
// For description of an algorithm for calculation of folded offset,
// see Tensor Library Documentation, Ondra Kamenik, 2005, description
// of FTensor::getOffsetRecurse().
int FoldMultiIndex::offset() const
{
// make copy for the recursions
int* tmp = new int[ord];
for (int i = 0; i < ord; i++)
tmp[i] = data[i];
// call the recursive algorithm
int res = offset_recurse(tmp, ord, nvar);
delete [] tmp;
return res;
}
void FoldMultiIndex::print() const
{
printf("[");
for (int i = 0; i < ord; i++)
printf("%d ", data[i]);
printf("]");
}
int FoldMultiIndex::offset_recurse(int* data, int len, int nv)
{
if (len == 0)
return 0;
// calculate length of initial constant indices
int prefix = 1;
while (prefix < len && data[0] == data[prefix])
prefix++;
int m = data[0];
int s1 = ptriang.noverk(nv+len-1, len) - ptriang.noverk(nv-m+len-1,len);
// cancel m from the rest of the sequence
for (int i = prefix; i < len; i++)
data[i] -= m;
// calculate offset of the remaining sequence
int s2 = offset_recurse(data+prefix, len-prefix, nv-m);
// return the sum
return s1+s2;
}
bool ltfmi::operator()(const FoldMultiIndex& i1, const FoldMultiIndex& i2) const
{
return i1 < i2;
}
FormulaDerEvaluator::FormulaDerEvaluator(const FormulaParser& fp)
: etree(fp.otree, -1)
{
for (unsigned int i = 0; i < fp.ders.size(); i++)
ders.push_back((const FormulaDerivatives*)(fp.ders[i]));
der_atoms = fp.atoms.variables();
}
void FormulaDerEvaluator::eval(const AtomValues& av, FormulaDerEvalLoader& loader, int order)
{
if (ders.size() == 0)
return;
int maxorder = ders[0]->order;
if (order > maxorder)
throw ogu::Exception(__FILE__,__LINE__,
"Wrong order in FormulaDerEvaluator::eval");
etree.reset_all();
av.setValues(etree);
int* vars = new int[order];
for (unsigned int i = 0; i < ders.size(); i++) {
for (FormulaDerivatives::Tfmiintmap::const_iterator it = ders[i]->ind2der.begin();
it != ders[i]->ind2der.end(); ++it) {
const FoldMultiIndex& mi = (*it).first;
if (mi.order() == order) {
// set vars from multiindex mi and variables
for (int k = 0; k < order; k++)
vars[k] = der_atoms[mi[k]];
// evaluate
double res = etree.eval(ders[i]->tder[(*it).second]);
// load
loader.load(i, order, vars, res);
}
}
}
delete [] vars;
}
void FormulaDerEvaluator::eval(const vector<int>& mp, const AtomValues& av,
FormulaDerEvalLoader& loader, int order)
{
etree.reset_all();
av.setValues(etree);
int nvar_glob = der_atoms.size();
int nvar = mp.size();
int* vars = new int[order];
for (unsigned int i = 0; i < ders.size(); i++) {
FoldMultiIndex mi(nvar, order);
do {
// find index of the derivative in the tensor
FoldMultiIndex mi_glob(nvar_glob, mi, mp);
int der = ders[i]->derivative(mi_glob);
if (der != OperationTree::zero) {
// set vars from the global multiindex
for (int k = 0; k < order; k++)
vars[k] = der_atoms[mi_glob[k]];
// evaluate derivative
double res = etree.eval(der);
// load
loader.load(i, order, vars, res);
}
mi.increment();
} while (! mi.past_the_end());
}
delete [] vars;
}
// Local Variables:
// mode:C++
// End:
// Copyright (C) 2005-2011, Ondra Kamenik
#ifndef OGP_FORMULA_PARSER_H
#define OGP_FORMULA_PARSER_H
#include "tree.h"
namespace ogp {
using std::vector;
/** Pure virtual class defining a minimal interface for
* representation of nulary terms within FormulaParser. */
class Atoms {
public:
Atoms() {}
virtual ~Atoms() {}
/** This returns previously assigned internal index to the
* given atom, or returns -1 if the atom has not been assigned
* yet. The method can raise an exception, if the Atoms
* implementation is strict and the name is not among
* prescribed possible values. */
virtual int check(const char* name) const = 0;
/** This method assigns an internal index to the nulary term
* described by the name. The internal index is allocated by
* OperationTree class. */
virtual void assign(const char* name, int t) = 0;
/** Returns a number of variables which will be used for
* differentiations. */
virtual int nvar() const = 0;
/** Returns a vector of variable's internal indices which will
* be used for differentiations. */
virtual vector<int> variables() const = 0;
/** Debug print. */
virtual void print() const = 0;
};
/** Pure virtual class defining interface for all classes able to
* set nulary terms to evaluation tree EvalTree. The
* implementations of this class will have to be connected with
* Atoms to have knowledge about the atoms and their indices in
* the tree, and will call EvalTree::set_nulary. */
class AtomValues {
public:
virtual ~AtomValues() {}
virtual void setValues(EvalTree& et) const = 0;
};
class FormulaDerEvaluator;
class FoldMultiIndex;
/** For ordering FoldMultiIndex in the std::map. */
struct ltfmi {
bool operator()(const FoldMultiIndex& i1, const FoldMultiIndex& i2) const;
};
/** This class stores derivatives (tree indices) of one formula
* for all orders upto a given one. It stores the derivatives as a
* sequence (vector) of these tree indices and sequence of the
* multidimensional indices of variables wrt which the derivatives
* were taken. In order to speed up querying for a derivative
* given the variables, we have a map mapping the multidimensional
* index to the order of the derivative in the sequence.
*
* The only reason we do not have only this map is that the
* iterators of the map do not survive the insertions to the map,
* and implementation of the constructor has to be very difficult.
*/
class FormulaDerivatives {
friend class FormulaDerEvaluator;
protected:
/** Vector of derivatives. This is a list of derivatives (tree
* indices), the ordering is given by the algorithm used to
* create it. Currently, it starts with zero-th derivative,
* the formula itself and carries with first order, second,
* etc. */
vector<int> tder;
/** Vector of multiindices corresponding to the vector of
* derivatives. */
vector<FoldMultiIndex> indices;
/** For retrieving derivatives via a multiindex, we have a map
* mapping a multiindex to a derivative in the tder
* ordering. This means that indices[ind2der[index]] == index. */
typedef map<FoldMultiIndex, int, ltfmi> Tfmiintmap;
Tfmiintmap ind2der;
/** The number of variables. */
int nvar;
/** The maximum order of derivatives. */
int order;
public:
/** The constructor allocates and fills the sequence of the
* indices of derivatives for a formula.
* @param otree the OperationTree for which all work is done
* and to which the derivatives are added.
* @param vars the vector of nulary terms in the tree; the
* derivatives are taken with respect to these variables in
* the ordering given by the vector.
* @param f the index of the formula being differentiated. The
* zero derivative is set to f.
* @param max_order the maximum order of differentiation.
*/
FormulaDerivatives(OperationTree& otree, const vector<int>& vars, int f, int max_order);
/** Copy constructor. */
FormulaDerivatives(const FormulaDerivatives& fd);
virtual ~FormulaDerivatives(){}
/** Random access to the derivatives via multiindex. */
int derivative(const FoldMultiIndex& mi) const;
/** Return the order. */
int get_order() const
{return order;}
/** Debug print. */
void print(const OperationTree& otree) const;
};
class FormulaEvaluator;
/** This class is able to parse a number of formulas and
* differentiate them. The life cycle of the object is as follows:
* After it is created, a few calls to parse will add formulas
* (zero derivatives) to the object. Then a method differentiate()
* can be called and a vector of pointers to derivatives for each
* formula is created. After this, no one should call other
* parse() or differentiate(). A const reference of the object can
* be used in constructors of FormulaEvaluator and
* FormulaDerEvaluator in order to evaluate formulas (zero
* derivatives) and higher derivatives resp. */
class FormulaParser {
friend class FormulaCustomEvaluator;
friend class FormulaDerEvaluator;
protected:
/** The OperationTree of all formulas, including derivatives. */
OperationTree otree;
/** Reference to Atoms. The Atoms are filled with nulary terms
* during execution of parse(). */
Atoms& atoms;
/** Vector of formulas (zero derivatives) in the order as they
* have been parsed. */
vector<int> formulas;
/** The vector to derivatives, each vector corresponds to a
* formula in the vector formulas. */
vector<FormulaDerivatives*> ders;
public:
/** Construct an empty formula parser. */
FormulaParser(Atoms& a)
: atoms(a) {}
/** Copy constructor using a different instance of Atoms. */
FormulaParser(const FormulaParser& fp, Atoms& a);
virtual ~FormulaParser();
/** Requires an addition of the formula; called from the
* parser. */
void add_formula(int t);
/** Requires an addition of the binary operation; called from
* the parser. */
int add_binary(code_t code, int t1, int t2);
/** Requires an addition of the unary operation; called from
* the parser. */
int add_unary(code_t code, int t);
/** Requires an addition of the nulary operation given by the
* string. The Atoms are consulted for uniquness and are given
* an internal index generated by the OperationTree. This is
* the channel through which the Atoms are filled. */
int add_nulary(const char* str);
/** Adds a derivative to the tree. This just calls
* OperationTree::add_derivative. */
int add_derivative(int t, int v)
{return otree.add_derivative(t, v);}
/** Adds a substitution. This just calls
* OperationTree::add_substitution. */
int add_substitution(int t, const map<int,int>& subst)
{return otree.add_substitution(t, subst);}
/** Add the substitution given by the map where left sides of
* substitutions come from another parser. The right sides are
* from this object. The given t is from the given parser fp. */
int add_substitution(int t, const map<int,int>& subst,
const FormulaParser& fp)
{return otree.add_substitution(t, subst, fp.otree);}
/** This adds formulas from the given parser with (possibly)
* different atoms applying substitutions from the given map
* mapping atoms from fp to atoms of the object. */
void add_subst_formulas(const map<int,int>& subst, const FormulaParser& fp);
/** Substitute formulas. For each i from 1 through all
* formulas, it adds a substitution of the i-th formula and
* make it to be i-th formula.*/
void substitute_formulas(const std::map<int,int>& subst);
/** This method turns the given term to nulary operation. It
* should be used with caution, since this method does not
* anything do with atoms, but usually some action is also
* needed (at leat to assign the tree index t to some
* atom). */
void nularify(int t)
{otree.nularify(t);}
/** Returns a set of nulary terms of the given term. Just
* calls OperationTree::nulary_of_term. */
const unordered_set<int>& nulary_of_term(int t) const
{return otree.nulary_of_term(t);}
/** Parse a given string containing one or more formulas. The
* formulas are parsed and added to the OperationTree and to
* the formulas vector. */
void parse(int length, const char* stream);
/** Processes a syntax error from bison. */
void error(const char* mes) const;
/** Differentiate all the formulas up to the given order. The
* variables with respect to which the derivatives are taken
* are obtained by Atoms::variables(). If the derivates exist,
* they are destroyed and created again (with possibly
* different order). */
void differentiate(int max_order);
/** Return i-th formula derivatives. */
const FormulaDerivatives& derivatives(int i) const;
/** This returns a maximum index of zero derivative formulas
* including all nulary terms. This is a mimumum length of the
* tree for which it is safe to evaluate zero derivatives of
* the formulas. */
int last_formula() const;
/** This returns a tree index of the i-th formula in the
* vector. */
int formula(int i) const
{return formulas[i];}
/** This returns a tree index of the last formula and pops its
* item from the formulas vector. The number of formulas is
* then less by one. Returns -1 if there is no formula. If
* there are derivatives of the last formula, they are
* destroyed and the vector ders is popped from the back. */
int pop_last_formula();
/** This returns a number of formulas. */
int nformulas() const
{return (int)(formulas.size());}
/** This returns a reference to atoms. */
const Atoms& getAtoms() const
{return atoms;}
Atoms& getAtoms()
{return atoms;}
/** This returns the tree. */
const OperationTree& getTree() const
{return otree;}
OperationTree& getTree()
{return otree;}
/** Debug print. */
void print() const;
private:
/** Hide this copy constructor declaration by declaring it as
* private. */
FormulaParser(const FormulaParser& fp);
/** Destroy all derivatives. */
void destroy_derivatives();
};
/** This is a pure virtual class defining an interface for all
* classes which will load the results of formula (zero
* derivative) evaluations. A primitive implementation of this
* class can be a vector of doubles. */
class FormulaEvalLoader {
public:
virtual ~FormulaEvalLoader() {}
/** Set the value res for the given formula. The formula is
* identified by an index corresponding to the ordering in
* which the formulas have been parsed (starting from
* zero). */
virtual void load(int i, double res) = 0;
};
/** This class evaluates a selected subset of terms of the
* tree. In the protected constructor, one can constraint the
* initialization of the evaluation tree to a given number of
* terms in the beginning. Using this constructor, one has to make
* sure, that the terms in the beginning do not refer to terms
* behind the initial part. */
class FormulaCustomEvaluator {
protected:
/** The evaluation tree. */
EvalTree etree;
/** The custom tree indices to be evaluated. */
vector<int> terms;
public:
/** Construct from FormulaParser and given list of terms. */
FormulaCustomEvaluator(const FormulaParser& fp, const vector<int>& ts)
: etree(fp.otree), terms(ts)
{}
/** Construct from OperationTree and given list of terms. */
FormulaCustomEvaluator(const OperationTree& ot, const vector<int>& ts)
: etree(ot), terms(ts)
{}
/** Evaluate the terms using the given AtomValues and load the
* results using the given loader. The loader is called for
* each term in the order of the terms. */
void eval(const AtomValues& av, FormulaEvalLoader& loader);
protected:
FormulaCustomEvaluator(const FormulaParser& fp)
: etree(fp.otree, fp.last_formula()), terms(fp.formulas)
{}
};
/** This class evaluates zero derivatives of the FormulaParser. */
class FormulaEvaluator : public FormulaCustomEvaluator {
public:
/** Construct from FormulaParser. */
FormulaEvaluator(const FormulaParser& fp)
: FormulaCustomEvaluator(fp) {}
};
/** This is a pure virtual class defining an interface for all
* classes which will load the results of formula derivative
* evaluations. */
class FormulaDerEvalLoader {
public:
virtual ~FormulaDerEvalLoader() {}
/** This loads the result of the derivative of the given
* order. The semantics of i is the same as in
* FormulaEvalLoader::load. The indices of variables with
* respect to which the derivative was taken are stored in
* memory pointed by vars. These are the tree indices of the
* variables. */
virtual void load(int i, int order, const int* vars, double res) = 0;
};
/** This class is a utility class representing the tensor
* multindex. It can basically increment itself, and calculate
* its offset in the folded tensor. */
class FoldMultiIndex {
/** Number of variables. */
int nvar;
/** Dimension. */
int ord;
/** The multiindex. */
int* data;
public:
/** Initializes to the zero derivative. Order is 0, data is
* empty. */
FoldMultiIndex(int nv);
/** Initializes the multiindex to zeros or given i. */
FoldMultiIndex(int nv, int order, int i = 0);
/** Makes a new multiindex of the same order applying a given
* mapping to the indices. The mapping is supposed to be monotone. */
FoldMultiIndex(int nv, const FoldMultiIndex& mi, const vector<int>& mp);
/** Shifting constructor. This adds a given number of orders
* to the end, copying the last item to the newly added items,
* keeping the index ordered. If the index was empty (zero-th
* dimension), then zeros are added. */
FoldMultiIndex(const FoldMultiIndex& fmi, int new_orders);
/** Copy constructor. */
FoldMultiIndex(const FoldMultiIndex& fmi);
/** Desctructor. */
virtual ~FoldMultiIndex()
{delete [] data;}
/** Assignment operator. */
const FoldMultiIndex& operator=(const FoldMultiIndex& fmi);
/** Operator < implementing lexicographic ordering within one
* order, increasing order across orders. */
bool operator<(const FoldMultiIndex& fmi) const;
bool operator==(const FoldMultiIndex& fmi) const;
/** Increment the multiindex. */
void increment();
/** Return offset of the multiindex in the folded tensor. */
int offset() const;
const int& operator[](int i) const
{return data[i];}
/** Return order of the multiindex, i.e. dimension of the
* tensor. */
int order() const
{return ord;}
/** Return the number of variables. */
int nv() const
{return nvar;}
/** Return the data. */
const int* ind() const
{return data;}
/** Return true if the end of the tensor is reached. The
* result of a subsequent increment should be considered
* unpredictable. */
bool past_the_end() const
{return (ord == 0) || (data[0] == nvar);}
/** Prints the multiindex in the brackets. */
void print() const;
private:
static int offset_recurse(int* data, int len, int nv);
};
/** This class evaluates derivatives of the FormulaParser. */
class FormulaDerEvaluator {
/** Its own instance of EvalTree. */
EvalTree etree;
/** The indices of derivatives for each formula. This is a
* const copy FormulaParser::ders. We do not allocate nor
* deallocate anything here. */
vector<const FormulaDerivatives*> ders;
/** A copy of tree indices corresponding to atoms to with
* respect the derivatives were taken. */
vector<int> der_atoms;
public:
/** Construct the object from FormulaParser. */
FormulaDerEvaluator(const FormulaParser& fp);
/** Evaluate the derivatives from the FormulaParser wrt to all
* atoms in variables vector at the given AtomValues. The
* given loader is used for output. */
void eval(const AtomValues& av, FormulaDerEvalLoader& loader, int order);
/** Evaluate the derivatives from the FormulaParser wrt to a
* selection of atoms of the atoms in der_atoms vector at the
* given AtomValues. The selection is given by a monotone
* mapping to the indices (not values) of the der_atoms. */
void eval(const vector<int>& mp, const AtomValues& av, FormulaDerEvalLoader& loader,
int order);
};
};
#endif
// Local Variables:
// mode:C++
// End:
// Copyright (C) 2006, Ondra Kamenik
// $Id: location.h 762 2006-05-22 13:00:07Z kamenik $
// Purpose: This file defines macros for lex and bison so that the
// very primitive location tracking would be enabled. The location of
// a token is given by offset of its first character. The offset is
// relative to the number which is (and must be) initialized before
// parsing. This file is to be included to the top of bison and lex
// sources.
// How to use: in preamble of bison and flex, you must include this
// file and declare extern YYLTYPE prefix##lloc. In addition, in flex,
// you must define int prefix##ll =0; and use macro SET_LLOC(prefix)
// in EVERY action consuming material (this can be done with #define
// YY_USER_ACTION) and in bison you must use option %locations.
#ifndef OG_LOCATION_H
#define OG_LOCATION_H
namespace ogp {
struct location_type {
int off; // offset of the token
int ll; // length ot the token
location_type() : off(0), ll(0) {}
};
};
#define YYLTYPE ogp::location_type
// set current off to the first off and add all lengths
#define YYLLOC_DEFAULT(Current, Rhs, N) \
{(Current).off = (Rhs)[1].off; \
(Current).ll = 0; \
for (int i = 1; i <= N; i++) (Current).ll += (Rhs)[i].ll;}
#define SET_LLOC(prefix) (prefix##lloc.off += prefix##lloc.ll, prefix##lloc.ll = prefix##leng)
#endif
// Local Variables:
// mode:C++
// End:
%{
// Copyright (C) 2006-2011, Ondra Kamenik
#include "location.h"
#include "matrix_tab.hh"
extern YYLTYPE matrix_lloc;
extern void matrix_error(const char*);
#define YY_USER_ACTION SET_LLOC(matrix_);
%}
%option nounput
%option noyy_top_state
%option stack
%option yylineno
%option prefix="matrix_"
%option never-interactive
%x CMT
%%
/* comments */
<*>"/*" {yy_push_state(CMT);}
<CMT>[^*\n]*
<CMT>"*"+[^*/\n]*
<CMT>"*"+"/" {yy_pop_state();}
<CMT>[\n]
"//".*\n
/* ignore spaces and commas */
[ \t,]
/* new row */
\r\n {return NEW_ROW;}
\n {return NEW_ROW;}
;[ \t]*\n {return NEW_ROW;}
;[ \t]*\r\n {return NEW_ROW;}
; {return NEW_ROW;}
[+-]?(([0-9]*\.?[0-9]+)|([0-9]+\.))([edED][-+]?[0-9]+)? {
matrix_lval.val = strtod(matrix_text, NULL);
return DNUMBER;
}
. {
char mes[300];
sprintf(mes, "Unrecognized character %s", matrix_text);
matrix_error(mes);
}
%%
int matrix_wrap()
{
return 1;
}
void matrix__destroy_buffer(void* p)
{
matrix__delete_buffer((YY_BUFFER_STATE)p);
}
// Copyright (C) 2006-2011, Ondra Kamenik
%{
#include "location.h"
#include "matrix_parser.h"
#include "matrix_tab.hh"
void matrix_error(const char*);
int matrix_lex(void);
extern int matrix_lineno;
extern ogp::MatrixParser* mparser;
extern YYLTYPE matrix_lloc;
// static void print_token_value (FILE *, int, YYSTYPE);
//#define YYPRINT(file, type, value) print_token_value (file, type, value)
%}
%union {
double val;
int integer;
}
%token NEW_ROW
%token <val> DNUMBER
%name-prefix="matrix_";
%locations
%error-verbose
%%
matrix : first_row other_rows
| first_row other_rows empty_rows
| first_row empty_rows other_rows empty_rows
| first_row empty_rows other_rows
| empty_rows first_row other_rows
| empty_rows first_row other_rows empty_rows
| empty_rows first_row empty_rows other_rows empty_rows
| empty_rows first_row empty_rows
| first_row empty_rows
| empty_rows first_row
| first_row
| empty_rows
;
empty_rows : empty_rows NEW_ROW | NEW_ROW;
lod : DNUMBER {mparser->add_item($1);}
| lod DNUMBER {mparser->add_item($2);}
;
first_row : lod;
other_rows : other_rows one_row | other_rows empty_rows one_row |one_row ;
one_row : NEW_ROW {mparser->start_row();} lod;
%%
void matrix_error(const char* s)
{
mparser->error(s);
}
// Copyright (C) 2006, Ondra Kamenik
// $Id: matrix_parser.cpp 2269 2008-11-23 14:33:22Z michel $
#include "parser_exception.h"
#include "matrix_parser.h"
#include "location.h"
#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
char* 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), c(0), r(mp.find_first_non_empty_row())
{}
MPIterator::MPIterator(const MatrixParser& mp, const char* dummy)
: p(&mp), i(mp.data.size()), c(0), 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;
public:
MatrixParser()
: nc(0) {}
MatrixParser(const MatrixParser& mp)
: data(mp.data), row_lengths(mp.row_lengths), nc(mp.nc) {}
virtual ~MatrixParser() {}
/** 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;
/** The index of the pointed item in the matrix parser. */
unsigned int i;
/** The column number of the pointed item starting from zero. */
int c;
/** The row number of the pointed item starting from zero. */
int r;
public:
MPIterator() : p(NULL), i(0), c(0), r(0) {}
/** 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. */
const MPIterator& operator=(const MPIterator& it)
{p = it.p; i = it.i; c = it.c; r = it.r; return *this;}
/** 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
// Local Variables:
// mode:C++
// End:
// Copyright (C) 2006, Ondra Kamenik
// $Id: namelist.cpp 42 2007-01-22 21:53:24Z ondra $
#include "namelist.h"
#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)
{
char* 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() {}
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
// Local Variables:
// mode:C++
// End: