diff --git a/parser.src/DataTree.cc b/parser.src/DataTree.cc
index 5a465b1ea9d005787318c8ce1f38a08454f6a783..f6d3f396232b84e6e3be3346101c01cdf3776af6 100644
--- a/parser.src/DataTree.cc
+++ b/parser.src/DataTree.cc
@@ -1,71 +1,348 @@
-/*! \file
- \version 1.0
- \date 04/09/2004
- \par This file implements the DataTree class methodes.
-*/
-//------------------------------------------------------------------------------
#include <iostream>
-#include <string>
-#include <vector>
-#include <map>
-using namespace std;
-#include <time.h>
-
-//------------------------------------------------------------------------------
-#include "DynareBison.hh"
-#include "VariableTable.hh"
-#include "NumericalConstants.hh"
+
#include "DataTree.hh"
-DataTree::DataTree(SymbolTable &symbol_table_arg) :
+DataTree::DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg) :
symbol_table(symbol_table_arg),
+ num_constants(num_constants_arg),
variable_table(symbol_table_arg),
- NoOpCode(-1), NullID(NULL)
-{
- offset = 1;
- current_order = 0;
-
- Zero = new MetaToken(reinterpret_cast <NodeID> (0), eNumericalConstant, NULL, -1);
- Zero->op_name = "";
- Zero->reference_count.resize(current_order+1,2);
- Zero->idx = 0;
- mModelTree.push_back(Zero);
- mIndexOfTokens[*Zero]=Zero;
-
- One = new MetaToken(reinterpret_cast <NodeID> (1), eNumericalConstant, NULL, -1); One->op_name = "";
- One->op_name = "";
- One->reference_count.resize(current_order+1,1);
- One->idx = 1;
- mModelTree.push_back(One);
- mIndexOfTokens[*One]=One;
-
- MinusOne = new MetaToken(One, eTempResult, NULL, token::UMINUS);
- MinusOne->op_name = OperatorTable::str(token::UMINUS);
- MinusOne->reference_count.resize(current_order+1,1);
- MinusOne->idx = 2;
- mModelTree.push_back(MinusOne);
- mIndexOfTokens[*MinusOne]=MinusOne;
-
- // Pushing "0=0" into mModelTree
- ZeroEqZero = new MetaToken(Zero, eTempResult, Zero, token::EQUAL);
- ZeroEqZero->op_name = OperatorTable::str(token::EQUAL);
- ZeroEqZero->reference_count.resize(current_order+1,1);
- ZeroEqZero->idx = 3;
- mModelTree.push_back(ZeroEqZero);
- mIndexOfTokens[*ZeroEqZero]=ZeroEqZero;
-
- // Initialise global node counter
- nodeCounter = 4;
-
- BeginModel = mModelTree.end();
- BeginModel--;
-}
-
-//------------------------------------------------------------------------------
+ offset(1)
+{
+ Zero = AddNumConstant("0.0");
+ One = AddNumConstant("1.0");
+
+ MinusOne = AddUMinus(One);
+}
+
DataTree::~DataTree()
{
- for (TreeIterator it = mModelTree.begin(); it != mModelTree.end(); it++)
+ for (node_list_type::iterator it = node_list.begin(); it != node_list.end(); it++)
+ delete *it;
+}
+
+NodeID
+DataTree::AddNumConstant(const string &value)
+{
+ int id = num_constants.AddConstant(value);
+
+ num_const_node_map_type::iterator it = num_const_node_map.find(id);
+ if (it != num_const_node_map.end())
+ return it->second;
+ else
+ return new NumConstNode(*this, id);
+}
+
+NodeID
+DataTree::AddVariable(const string &name, int lag)
+{
+ if (!symbol_table.Exist(name))
+ {
+ cerr << "Unknown symbol: " << name << endl;
+ exit(-1);
+ }
+
+ symbol_table.SetReferenced(name);
+
+ Type type = symbol_table.getType(name);
+ int id;
+ if (type == eEndogenous
+ || type == eExogenousDet
+ || type == eExogenous
+ || type == eRecursiveVariable)
+ id = variable_table.AddVariable(name, lag);
+ else
+ id = symbol_table.getID(name);
+
+ variable_node_map_type::iterator it = variable_node_map.find(make_pair(id, type));
+ if (it != variable_node_map.end())
+ return it->second;
+ else
+ return new VariableNode(*this, id, type);
+}
+
+NodeID
+DataTree::AddPlus(NodeID iArg1, NodeID iArg2)
+{
+ if (iArg1 != Zero && iArg2 != Zero)
+ {
+ // Simplify x+(-y) in x-y
+ UnaryOpNode *uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
+ if (uarg2 != NULL && uarg2->op_code == oUminus)
+ return AddMinus(iArg1, uarg2->arg);
+
+ // To treat commutativity of "+"
+ // Nodes iArg1 and iArg2 are sorted by index
+ if (iArg1->idx > iArg2->idx)
+ {
+ NodeID tmp = iArg1;
+ iArg1 = iArg2;
+ iArg2 = tmp;
+ }
+ return AddBinaryOp(iArg1, oPlus, iArg2);
+ }
+ else if (iArg1 != Zero)
+ return iArg1;
+ else if (iArg2 != Zero)
+ return iArg2;
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddMinus(NodeID iArg1, NodeID iArg2)
+{
+ if (iArg1 != Zero && iArg2 != Zero)
+ return AddBinaryOp(iArg1, oMinus, iArg2);
+ else if (iArg1 != Zero)
+ return iArg1;
+ else if (iArg2 != Zero)
+ return AddUMinus(iArg2);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddUMinus(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ {
+ // Simplify -(-x) in x
+ UnaryOpNode *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
+ if (uarg != NULL && uarg->op_code == oUminus)
+ return uarg->arg;
+
+ return AddUnaryOp(oUminus, iArg1);
+ }
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddTimes(NodeID iArg1, NodeID iArg2)
+{
+ if (iArg1 == MinusOne)
+ return AddUMinus(iArg2);
+ else if (iArg2 == MinusOne)
+ return AddUMinus(iArg1);
+ else if (iArg1 != Zero && iArg1 != One && iArg2 != Zero && iArg2 != One)
+ {
+ // To treat commutativity of "*"
+ // Nodes iArg1 and iArg2 are sorted by index
+ if (iArg1->idx > iArg2->idx)
+ {
+ NodeID tmp = iArg1;
+ iArg1 = iArg2;
+ iArg2 = tmp;
+ }
+ return AddBinaryOp(iArg1, oTimes, iArg2);
+ }
+ else if (iArg1 != Zero && iArg1 != One && iArg2 == One)
+ return iArg1;
+ else if (iArg2 != Zero && iArg2 != One && iArg1 == One)
+ return iArg2;
+ else if (iArg2 == One && iArg1 == One)
+ return One;
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddDivide(NodeID iArg1, NodeID iArg2)
+{
+ if (iArg1 != Zero && iArg2 != Zero && iArg2 != One)
+ return AddBinaryOp(iArg1, oDivide, iArg2);
+ else if (iArg2 == One)
+ return iArg1;
+ else if (iArg1 == Zero && iArg2 != Zero)
+ return Zero;
+ else
{
- if (*it != NullID) delete *it;
+ cerr << "Division by zero!" << endl;
+ exit(-1);
}
}
+
+NodeID
+DataTree::AddPower(NodeID iArg1, NodeID iArg2)
+{
+ if (iArg1 != Zero && iArg2 != Zero && iArg2 != One)
+ return AddBinaryOp(iArg1, oPower, iArg2);
+ else if (iArg2 == One)
+ return iArg1;
+ else if (iArg2 == Zero)
+ return One;
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddExp(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oExp, iArg1);
+ else
+ return One;
+}
+
+NodeID
+DataTree::AddLog(NodeID iArg1)
+{
+ if (iArg1 != Zero && iArg1 != One)
+ return AddUnaryOp(oLog, iArg1);
+ else if (iArg1 == One)
+ return Zero;
+ else
+ {
+ cerr << "log(0) isn't available" << endl;
+ exit(-1);
+ }
+}
+
+NodeID DataTree::AddLog10(NodeID iArg1)
+{
+ if (iArg1 != Zero && iArg1 != One)
+ return AddUnaryOp(oLog10, iArg1);
+ else if (iArg1 == One)
+ return Zero;
+ else
+ {
+ cerr << "log10(0) isn't available" << endl;
+ exit(-1);
+ }
+}
+
+NodeID
+DataTree::AddCos(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oCos, iArg1);
+ else
+ return One;
+}
+
+NodeID
+DataTree::AddSin(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oSin, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddTan(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oTan, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddACos(NodeID iArg1)
+{
+ if (iArg1 != One)
+ return AddUnaryOp(oAcos, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddASin(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oAsin, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddATan(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oAtan, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddCosH(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oCosh, iArg1);
+ else
+ return One;
+}
+
+NodeID
+DataTree::AddSinH(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oSinh, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddTanH(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oTanh, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddACosH(NodeID iArg1)
+{
+ if (iArg1 != One)
+ return AddUnaryOp(oAcosh, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddASinH(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oAsinh, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddATanH(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oAtanh, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddSqRt(NodeID iArg1)
+{
+ if (iArg1 != Zero)
+ return AddUnaryOp(oSqrt, iArg1);
+ else
+ return Zero;
+}
+
+NodeID
+DataTree::AddEqual(NodeID iArg1, NodeID iArg2)
+{
+ return AddBinaryOp(iArg1, oEqual, iArg2);
+}
+
+void
+DataTree::AddLocalParameter(const string &name, NodeID value) throw (LocalParameterException)
+{
+ int id = symbol_table.getID(name);
+
+ // Throw an exception if symbol already declared
+ map<int, NodeID>::iterator it = local_parameters_table.find(id);
+ if (it != local_parameters_table.end())
+ throw LocalParameterException(name);
+
+ local_parameters_table[id] = value;
+}
diff --git a/parser.src/DynareBison.yy b/parser.src/DynareBison.yy
index 5f55985f81b976607bdb5e03e1440d62ee19193d..58d3fca0493076f2138e389d387aabce23c14295 100644
--- a/parser.src/DynareBison.yy
+++ b/parser.src/DynareBison.yy
@@ -8,7 +8,7 @@ using namespace std;
class ParsingDriver;
#include "SymbolTableTypes.hh"
-#include "ModelTypes.hh"
+#include "ExprNode.hh"
//! Type for semantic value of non-derivable expressions
typedef pair<int, Type> ExpObj;
@@ -74,9 +74,6 @@ typedef pair<int, Type> ExpObj;
%left UMINUS
%nonassoc POWER
%token EXP LOG LOG10 SIN COS TAN ASIN ACOS ATAN SINH COSH TANH ASINH ACOSH ATANH SQRT
-/* isn't parsed from the *.mod file, but used to distinguish EQUAL in equation and EQUAL in assignment in operation codes
-*/
-%token ASSIGN
%type <exp_val> expression comma_expression
%type <model_val> equation hand_side model_var
diff --git a/parser.src/ExprNode.cc b/parser.src/ExprNode.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1662179823567ff92c162dba195c49eb91715cab
--- /dev/null
+++ b/parser.src/ExprNode.cc
@@ -0,0 +1,776 @@
+#include <iostream>
+#include <iterator>
+#include <algorithm>
+
+#include "ExprNode.hh"
+#include "DataTree.hh"
+
+ExprNode::ExprNode(DataTree &datatree_arg) : datatree(datatree_arg)
+{
+ // Add myself to datatree
+ datatree.node_list.push_back(this);
+
+ // Set my index and increment counter
+ idx = datatree.node_counter++;
+}
+
+ExprNode::~ExprNode()
+{
+}
+
+NodeID
+ExprNode::getDerivative(int varID)
+{
+ // Return zero if derivative is necessarily null (using symbolic a priori)
+ set<int>::const_iterator it = non_null_derivatives.find(varID);
+ if (it == non_null_derivatives.end())
+ return datatree.Zero;
+
+ // If derivative is stored in cache, use the cached value, otherwise compute it (and cache it)
+ map<int, NodeID>::const_iterator it2 = derivatives.find(varID);
+ if (it2 != derivatives.end())
+ return it2->second;
+ else
+ {
+ NodeID d = computeDerivative(varID);
+ derivatives[varID] = d;
+ return d;
+ }
+}
+
+int
+ExprNode::precedence(const temporary_terms_type &temporary_terms) const
+{
+ // For a constant, a variable, or a unary op, the precedence is maximal
+ return 100;
+}
+
+int
+ExprNode::cost(temporary_terms_type &temporary_terms) const
+{
+ // For a terminal node, the cost is null
+ return 0;
+}
+
+void
+ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+ temporary_terms_type &temporary_terms) const
+{
+ // Nothing to do for a terminal node
+}
+
+NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
+ ExprNode(datatree_arg),
+ id(id_arg)
+{
+ // Add myself to the num const map
+ datatree.num_const_node_map[id] = this;
+
+ // All derivatives are null, so non_null_derivatives is left empty
+}
+
+NodeID
+NumConstNode::computeDerivative(int varID)
+{
+ return datatree.Zero;
+}
+
+void
+NumConstNode::writeOutput(ostream &output, bool is_dynamic,
+ const temporary_terms_type &temporary_terms) const
+{
+ temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this));
+ if (it != temporary_terms.end())
+ output << "T" << idx;
+ else
+ output << datatree.num_constants.get(id);
+}
+
+VariableNode::VariableNode(DataTree &datatree_arg, int id_arg, Type type_arg) :
+ ExprNode(datatree_arg),
+ id(id_arg),
+ type(type_arg)
+{
+ // Add myself to the variable map
+ datatree.variable_node_map[make_pair(id, type)] = this;
+
+ // Fill in non_null_derivatives
+ switch(type)
+ {
+ case eEndogenous:
+ case eExogenous:
+ case eExogenousDet:
+ case eRecursiveVariable:
+ // For a variable, the only non-null derivative is with respect to itself
+ non_null_derivatives.insert(id);
+ break;
+ case eParameter:
+ // All derivatives are null, do nothing
+ break;
+ case eLocalParameter:
+ // Non null derivatives are those of the value of the local parameter
+ non_null_derivatives = datatree.local_parameters_table[id]->non_null_derivatives;
+ break;
+ case eNumericalConstant:
+ case eUNDEF:
+ case eTempResult:
+ // Impossible cases
+ cerr << "Incorrect symbol type used in VariableNode" << endl;
+ exit(-1);
+ }
+}
+
+NodeID
+VariableNode::computeDerivative(int varID)
+{
+ switch(type)
+ {
+ case eEndogenous:
+ case eExogenous:
+ case eExogenousDet:
+ case eRecursiveVariable:
+ if (varID == id)
+ return datatree.One;
+ else
+ return datatree.Zero;
+ case eParameter:
+ return datatree.Zero;
+ case eLocalParameter:
+ return datatree.local_parameters_table[id]->getDerivative(varID);
+ case eNumericalConstant:
+ case eUNDEF:
+ case eTempResult:
+ // Impossible cases
+ cerr << "Incorrect symbol type used in VariableNode" << endl;
+ exit(-1);
+ }
+ cerr << "Impossible case!" << endl;
+ exit(-1);
+}
+
+void
+VariableNode::writeOutput(ostream &output, bool is_dynamic,
+ const temporary_terms_type &temporary_terms) const
+{
+ // If node is a temporary term
+ temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<VariableNode *>(this));
+ if (it != temporary_terms.end())
+ {
+ output << "T" << idx;
+ return;
+ }
+
+ char &lpar = datatree.lpar;
+ char &rpar = datatree.rpar;
+
+ int idx;
+ switch(type)
+ {
+ case eParameter:
+ output << "params" << lpar << id + datatree.offset << rpar;
+ break;
+ case eLocalParameter:
+ output << datatree.symbol_table.getNameByID(eLocalParameter, id);
+ break;
+ case eEndogenous:
+ if (is_dynamic)
+ {
+ idx = datatree.variable_table.getPrintIndex(id) + datatree.offset;
+ output << "y" << lpar << idx << rpar;
+ }
+ else
+ {
+ idx = datatree.variable_table.getSymbolID(id) + datatree.offset;
+ output << "y" << lpar << idx << rpar;
+ }
+ break;
+ case eExogenous:
+ idx = datatree.variable_table.getSymbolID(id) + datatree.offset;
+ if (is_dynamic)
+ {
+ int lag = datatree.variable_table.getLag(id);
+ if (datatree.offset == 1)
+ if (lag != 0)
+ output << "x" << lpar << "it_ + " << lag << ", " << idx << rpar;
+ else
+ output << "x" << lpar << "it_, " << idx << rpar;
+ else
+ if (lag != 0)
+ output << "x" << lpar << "it_+" << lag << "+" << idx << "*nb_row_x" << rpar;
+ else
+ output << "x" << lpar << "it_+" << idx << "*nb_row_x" << rpar;
+ }
+ else
+ output << "x" << lpar << idx << rpar;
+ break;
+ case eExogenousDet:
+ idx = datatree.variable_table.getSymbolID(id) + datatree.symbol_table.exo_nbr
+ + datatree.offset;
+ if (is_dynamic)
+ {
+ int lag = datatree.variable_table.getLag(id);
+ if (datatree.offset == 1)
+ if (lag != 0)
+ output << "x" << lpar << "it_ + " << lag << ", " << idx << rpar;
+ else
+ output << "x" << lpar << "it_, " << idx << rpar;
+ else
+ if (lag != 0)
+ output << "x" << lpar << "it_ + " << lag << "+" << idx << "*nb_row_xd" << rpar;
+ else
+ output << "x" << lpar << "it_+" << idx << "*nb_row_xd" << rpar;
+ }
+ else
+ output << "x" << lpar << idx << rpar;
+ break;
+ case eRecursiveVariable:
+ cerr << "Recursive variable not implemented" << endl;
+ exit(-1);
+ case eTempResult:
+ case eNumericalConstant:
+ case eUNDEF:
+ // Impossible cases
+ cerr << "Incorrect symbol type used in VariableNode" << endl;
+ exit(-1);
+ }
+}
+
+UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) :
+ ExprNode(datatree_arg),
+ arg(arg_arg),
+ op_code(op_code_arg)
+{
+ // Add myself to the unary op map
+ datatree.unary_op_node_map[make_pair(arg, op_code)] = this;
+
+ // Non-null derivatives are those of the argument
+ non_null_derivatives = arg->non_null_derivatives;
+}
+
+NodeID
+UnaryOpNode::computeDerivative(int varID)
+{
+ NodeID darg = arg->getDerivative(varID);
+
+ NodeID t11, t12, t13;
+
+ switch(op_code)
+ {
+ case oUminus:
+ return datatree.AddUMinus(darg);
+ case oExp:
+ return datatree.AddTimes(darg, this);
+ case oLog:
+ return datatree.AddDivide(darg, arg);
+ case oLog10:
+ t11 = datatree.AddExp(datatree.One);
+ t12 = datatree.AddLog10(t11);
+ t13 = datatree.AddDivide(darg, arg);
+ return datatree.AddTimes(t12, t13);
+ case oCos:
+ t11 = datatree.AddSin(arg);
+ t12 = datatree.AddUMinus(t11);
+ return datatree.AddTimes(darg, t12);
+ case oSin:
+ t11 = datatree.AddCos(arg);
+ return datatree.AddTimes(darg, t11);
+ case oTan:
+ t11 = datatree.AddTimes(this, this);
+ t12 = datatree.AddPlus(t11, datatree.One);
+ return datatree.AddTimes(darg, t12);
+ case oAcos:
+ t11 = datatree.AddSin(this);
+ t12 = datatree.AddDivide(darg, t11);
+ return datatree.AddUMinus(t12);
+ case oAsin:
+ t11 = datatree.AddCos(this);
+ return datatree.AddDivide(darg, t11);
+ case oAtan:
+ t11 = datatree.AddTimes(arg, arg);
+ t12 = datatree.AddPlus(datatree.One, t11);
+ return datatree.AddDivide(darg, t12);
+ case oCosh:
+ t11 = datatree.AddSinH(arg);
+ return datatree.AddTimes(darg, t11);
+ case oSinh:
+ t11 = datatree.AddCosH(arg);
+ return datatree.AddTimes(darg, t11);
+ case oTanh:
+ t11 = datatree.AddTimes(this, this);
+ t12 = datatree.AddMinus(datatree.One, t11);
+ return datatree.AddTimes(darg, t12);
+ case oAcosh:
+ t11 = datatree.AddSinH(this);
+ return datatree.AddDivide(darg, t11);
+ case oAsinh:
+ t11 = datatree.AddCosH(this);
+ return datatree.AddDivide(darg, t11);
+ case oAtanh:
+ t11 = datatree.AddTimes(arg, arg);
+ t12 = datatree.AddMinus(datatree.One, t11);
+ return datatree.AddTimes(darg, t12);
+ case oSqrt:
+ t11 = datatree.AddPlus(this, this);
+ return datatree.AddDivide(darg, t11);
+ }
+ cerr << "Impossible case!" << endl;
+ exit(-1);
+}
+
+int
+UnaryOpNode::cost(temporary_terms_type &temporary_terms) const
+{
+ // For a temporary term, the cost is null
+ temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
+ if (it != temporary_terms.end())
+ return 0;
+
+ int cost = arg->cost(temporary_terms);
+
+ if (datatree.offset)
+ // Cost for Matlab files
+ switch(op_code)
+ {
+ case oUminus:
+ return cost + 70;
+ case oExp:
+ return cost + 160;
+ case oLog:
+ return cost + 300;
+ case oLog10:
+ return cost + 16000;
+ case oCos:
+ case oSin:
+ case oCosh:
+ return cost + 210;
+ case oTan:
+ return cost + 230;
+ case oAcos:
+ return cost + 300;
+ case oAsin:
+ return cost + 310;
+ case oAtan:
+ return cost + 140;
+ case oSinh:
+ return cost + 240;
+ case oTanh:
+ return cost + 190;
+ case oAcosh:
+ return cost + 770;
+ case oAsinh:
+ return cost + 460;
+ case oAtanh:
+ return cost + 350;
+ case oSqrt:
+ return cost + 570;
+ }
+ else
+ // Cost for C files
+ switch(op_code)
+ {
+ case oUminus:
+ return cost + 3;
+ case oExp:
+ case oAcosh:
+ return cost + 210;
+ case oLog:
+ return cost + 137;
+ case oLog10:
+ return cost + 139;
+ case oCos:
+ case oSin:
+ return cost + 160;
+ case oTan:
+ return cost + 170;
+ case oAcos:
+ case oAtan:
+ return cost + 190;
+ case oAsin:
+ return cost + 180;
+ case oCosh:
+ case oSinh:
+ case oTanh:
+ return cost + 240;
+ case oAsinh:
+ return cost + 220;
+ case oAtanh:
+ return cost + 150;
+ case oSqrt:
+ return cost + 90;
+ }
+ cerr << "Impossible case!" << endl;
+ exit(-1);
+}
+
+void
+UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+ temporary_terms_type &temporary_terms) const
+{
+ NodeID this2 = const_cast<UnaryOpNode *>(this);
+
+ map<NodeID, int>::iterator it = reference_count.find(this2);
+ if (it == reference_count.end())
+ {
+ reference_count[this2] = 1;
+ arg->computeTemporaryTerms(reference_count, temporary_terms);
+ }
+ else
+ {
+ reference_count[this2]++;
+ if (reference_count[this2] * cost(temporary_terms) > datatree.min_cost)
+ temporary_terms.insert(this2);
+ }
+}
+
+void
+UnaryOpNode::writeOutput(ostream &output, bool is_dynamic,
+ const temporary_terms_type &temporary_terms) const
+{
+ // If node is a temporary term
+ temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
+ if (it != temporary_terms.end())
+ {
+ output << "T" << idx;
+ return;
+ }
+
+ // Always put parenthesis around uminus nodes
+ if (op_code == oUminus)
+ output << "(";
+
+ switch(op_code)
+ {
+ case oUminus:
+ output << "-";
+ break;
+ case oExp:
+ output << "exp";
+ break;
+ case oLog:
+ output << "log";
+ break;
+ case oLog10:
+ output << "log10";
+ break;
+ case oCos:
+ output << "cos";
+ break;
+ case oSin:
+ output << "sin";
+ break;
+ case oTan:
+ output << "tan";
+ break;
+ case oAcos:
+ output << "acos";
+ break;
+ case oAsin:
+ output << "asin";
+ break;
+ case oAtan:
+ output << "atan";
+ break;
+ case oCosh:
+ output << "cosh";
+ break;
+ case oSinh:
+ output << "sinh";
+ break;
+ case oTanh:
+ output << "tanh";
+ break;
+ case oAcosh:
+ output << "acosh";
+ break;
+ case oAsinh:
+ output << "asinh";
+ break;
+ case oAtanh:
+ output << "atanh";
+ break;
+ case oSqrt:
+ output << "sqrt";
+ break;
+ }
+
+ bool close_parenthesis = false;
+
+ /* Enclose argument with parentheses if:
+ - current opcode is not uminus, or
+ - current opcode is uminus and argument has lowest precedence
+ */
+ if (op_code != oUminus
+ || (op_code == oUminus && arg->precedence(temporary_terms) < precedence(temporary_terms)))
+ {
+ output << "(";
+ close_parenthesis = true;
+ }
+
+ // Write argument
+ arg->writeOutput(output, is_dynamic, temporary_terms);
+
+ if (close_parenthesis)
+ output << ")";
+
+ // Close parenthesis for uminus
+ if (op_code == oUminus)
+ output << ")";
+}
+
+BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
+ BinaryOpcode op_code_arg, const NodeID arg2_arg) :
+ ExprNode(datatree_arg),
+ arg1(arg1_arg),
+ arg2(arg2_arg),
+ op_code(op_code_arg)
+{
+ datatree.binary_op_node_map[make_pair(make_pair(arg1, arg2), op_code)] = this;
+
+ // Non-null derivatives are the union of those of the arguments
+ // Compute set union of arg1->non_null_derivatives and arg2->non_null_derivatives
+ set_union(arg1->non_null_derivatives.begin(),
+ arg1->non_null_derivatives.end(),
+ arg2->non_null_derivatives.begin(),
+ arg2->non_null_derivatives.end(),
+ inserter(non_null_derivatives, non_null_derivatives.begin()));
+}
+
+NodeID
+BinaryOpNode::computeDerivative(int varID)
+{
+ NodeID darg1 = arg1->getDerivative(varID);
+ NodeID darg2 = arg2->getDerivative(varID);
+
+ NodeID t11, t12, t13, t14, t15;
+
+ switch(op_code)
+ {
+ case oPlus:
+ return datatree.AddPlus(darg1, darg2);
+ case oMinus:
+ return datatree.AddMinus(darg1, darg2);
+ case oTimes:
+ t11 = datatree.AddTimes(darg1, arg2);
+ t12 = datatree.AddTimes(darg2, arg1);
+ return datatree.AddPlus(t11, t12);
+ case oDivide:
+ t11 = datatree.AddTimes(darg1, arg2);
+ t12 = datatree.AddTimes(darg2, arg1);
+ t13 = datatree.AddMinus(t11, t12);
+ t14 = datatree.AddTimes(arg2, arg2);
+ return datatree.AddDivide(t13, t14);
+ case oPower:
+ if (darg2 == datatree.Zero)
+ {
+ if (darg1 == datatree.Zero)
+ return datatree.Zero;
+ else
+ {
+ t11 = datatree.AddMinus(arg2, datatree.One);
+ t12 = datatree.AddPower(arg1, t11);
+ t13 = datatree.AddTimes(arg2, t12);
+ return datatree.AddTimes(darg1, t13);
+ }
+ }
+ else
+ {
+ t11 = datatree.AddLog(arg1);
+ t12 = datatree.AddTimes(darg2, t11);
+ t13 = datatree.AddTimes(darg1, arg2);
+ t14 = datatree.AddDivide(t13, arg1);
+ t15 = datatree.AddPlus(t12, t14);
+ return datatree.AddTimes(t15, this);
+ }
+ case oEqual:
+ return datatree.AddMinus(darg1, darg2);
+ }
+ cerr << "Impossible case!" << endl;
+ exit(-1);
+}
+
+int
+BinaryOpNode::precedence(const temporary_terms_type &temporary_terms) const
+{
+ temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
+ // A temporary term behaves as a variable
+ if (it != temporary_terms.end())
+ return 100;
+
+ switch(op_code)
+ {
+ case oEqual:
+ case oPlus:
+ case oMinus:
+ return 0;
+ case oTimes:
+ case oDivide:
+ return 1;
+ case oPower:
+ if (datatree.offset)
+ // In C, power operator is of the form pow(a, b)
+ return 100;
+ else
+ return 3;
+ }
+ cerr << "Impossible case!" << endl;
+ exit(-1);
+}
+
+int
+BinaryOpNode::cost(temporary_terms_type &temporary_terms) const
+{
+ temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
+ // For a temporary term, the cost is null
+ if (it != temporary_terms.end())
+ return 0;
+
+ int cost = arg1->cost(temporary_terms);
+ cost += arg2->cost(temporary_terms);
+
+ if (datatree.offset)
+ // Cost for Matlab files
+ switch(op_code)
+ {
+ case oPlus:
+ case oMinus:
+ case oTimes:
+ return cost + 90;
+ case oDivide:
+ return cost + 990;
+ case oPower:
+ return cost + 1160;
+ case oEqual:
+ return cost;
+ }
+ else
+ // Cost for C files
+ switch(op_code)
+ {
+ case oPlus:
+ case oMinus:
+ case oTimes:
+ return cost + 4;
+ case oDivide:
+ return cost + 15;
+ case oPower:
+ return cost + 520;
+ case oEqual:
+ return cost;
+ }
+ cerr << "Impossible case!" << endl;
+ exit(-1);
+}
+
+void
+BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+ temporary_terms_type &temporary_terms) const
+{
+ NodeID this2 = const_cast<BinaryOpNode *>(this);
+ map<NodeID, int>::iterator it = reference_count.find(this2);
+ if (it == reference_count.end())
+ {
+ // If this node has never been encountered, set its ref count to one,
+ // and travel through its children
+ reference_count[this2] = 1;
+ arg1->computeTemporaryTerms(reference_count, temporary_terms);
+ arg2->computeTemporaryTerms(reference_count, temporary_terms);
+ }
+ else
+ {
+ // If the node has already been encountered, increment its ref count
+ // and declare it as a temporary term if it is too costly
+ reference_count[this2]++;
+ if (reference_count[this2] * cost(temporary_terms) > datatree.min_cost)
+ temporary_terms.insert(this2);
+ }
+}
+
+void
+BinaryOpNode::writeOutput(ostream &output, bool is_dynamic,
+ const temporary_terms_type &temporary_terms) const
+{
+ // If current node is a temporary term
+ temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
+ if (it != temporary_terms.end())
+ {
+ output << "T" << idx;
+ return;
+ }
+
+ // Treat special case of power operator in C
+ if (op_code == oPower && datatree.offset == 0)
+ {
+ output << "pow(";
+ arg1->writeOutput(output, is_dynamic, temporary_terms);
+ output << ",";
+ arg2->writeOutput(output, is_dynamic, temporary_terms);
+ output << ")";
+ return;
+ }
+
+ int prec = precedence(temporary_terms);
+
+ bool close_parenthesis = false;
+
+ // If left argument has a lower precedence, or if current and left argument are both power operators, add parenthesis around left argument
+ BinaryOpNode *barg1 = dynamic_cast<BinaryOpNode *>(arg1);
+ if (arg1->precedence(temporary_terms) < prec
+ || (op_code == oPower && barg1 != NULL && barg1->op_code == oPower))
+ {
+ output << "(";
+ close_parenthesis = true;
+ }
+
+ // Write left argument
+ arg1->writeOutput(output, is_dynamic, temporary_terms);
+
+ if (close_parenthesis)
+ output << ")";
+
+ // Write current operator symbol
+ switch(op_code)
+ {
+ case oPlus:
+ output << "+";
+ break;
+ case oMinus:
+ output << "-";
+ break;
+ case oTimes:
+ output << "*";
+ break;
+ case oDivide:
+ output << "/";
+ break;
+ case oPower:
+ output << "^";
+ break;
+ case oEqual:
+ output << "=";
+ break;
+ }
+
+ close_parenthesis = false;
+
+ /* Add parenthesis around right argument if:
+ - its precedence is lower than those of the current node
+ - it is a power operator and current operator is also a power operator
+ - it is a minus operator with same precedence than current operator
+ - it is a divide operator with same precedence than current operator */
+ BinaryOpNode *barg2 = dynamic_cast<BinaryOpNode *>(arg2);
+ int arg2_prec = arg2->precedence(temporary_terms);
+ if (arg2_prec < prec
+ || (op_code == oPower && barg2 != NULL && barg2->op_code == oPower)
+ || (op_code == oMinus && arg2_prec == prec)
+ || (op_code == oDivide && arg2_prec == prec))
+ {
+ output << "(";
+ close_parenthesis = true;
+ }
+
+ // Write right argument
+ arg2->writeOutput(output, is_dynamic, temporary_terms);
+
+ if (close_parenthesis)
+ output << ")";
+}
diff --git a/parser.src/Makefile b/parser.src/Makefile
index 02f1e61ab7aa165b02ef4f882d7508471e674664..55106485069c6b5bcbb42ce41dc913c2d272be4e 100644
--- a/parser.src/Makefile
+++ b/parser.src/Makefile
@@ -50,7 +50,8 @@ COMMON_OBJ=\
ParsingDriver.o\
DataTree.o \
ModFile.o \
- Statement.o
+ Statement.o \
+ ExprNode.o
MATLAB_OBJ = InterfaceMatlab.o
diff --git a/parser.src/ModFile.cc b/parser.src/ModFile.cc
index 7f0991f1b46f6eac474a857a7b12812fe2cde1e4..cd62be38b057af91034b12e11d6db5e9f5a1528e 100644
--- a/parser.src/ModFile.cc
+++ b/parser.src/ModFile.cc
@@ -1,3 +1,6 @@
+#include <iostream>
+#include <fstream>
+
#include "ModFile.hh"
#include "Interface.hh"
@@ -56,7 +59,7 @@ ModFile::computingPass()
}
void
-ModFile::writeOutputFiles(const string &basename, bool clear_all)
+ModFile::writeOutputFiles(const string &basename, bool clear_all) const
{
ofstream mOutputFile;
@@ -135,7 +138,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all)
model_tree.writeDynamicFile(basename);
// Print statements
- for(vector<Statement *>::iterator it = statements.begin();
+ for(vector<Statement *>::const_iterator it = statements.begin();
it != statements.end(); it++)
(*it)->writeOutput(mOutputFile, basename);
diff --git a/parser.src/ModelTree.cc b/parser.src/ModelTree.cc
index 9b722c4f3c8b009022dabb9f23aa1dffe3aad398..4c0fc49c16cbc50250dc30d2b781d4c9401afd8d 100644
--- a/parser.src/ModelTree.cc
+++ b/parser.src/ModelTree.cc
@@ -1,43 +1,13 @@
-/*! \file
- \version 1.0
- \date 04/09/2004
- \par This file implements the ModelTree class methodes.
-*/
-//------------------------------------------------------------------------------
#include <iostream>
-#include <string>
-#include <vector>
-#include <stack>
-#include <numeric>
-#include <stdio.h>
-#include <map>
-#include <time.h>
-using namespace std;
-//------------------------------------------------------------------------------
-#include "VariableTable.hh"
-#include "DynareBison.hh"
-#include "NumericalConstants.hh"
+#include <fstream>
+#include <sstream>
+
#include "ModelTree.hh"
#include "Interface.hh"
-inline NodeID MetaToken::getDerivativeAddress(int iVarID, const ModelTree &model_tree) const
-{
- std::map<int, NodeID, std::less<int> >::const_iterator iter = d1.find(iVarID);
- if (iter == d1.end())
- // No entry in map, derivative is therefore null
- if (op_code == token::EQUAL)
- return model_tree.ZeroEqZero;
- else
- return model_tree.Zero;
- else
- return iter->second;
-}
-
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
- const NumericalConstants &num_constants_arg) :
- DataTree(symbol_table_arg),
- num_constants(num_constants_arg),
- eq_nbr(0),
+ NumericalConstants &num_constants_arg) :
+ DataTree(symbol_table_arg, num_constants_arg),
computeJacobian(false),
computeJacobianExo(false),
computeHessian(false),
@@ -46,7 +16,122 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
}
void
-ModelTree::writeStaticMFile(const string &static_basename)
+ModelTree::derive(int order)
+{
+ cout << "Processing derivation ..." << endl;
+
+ cout << " Processing Order 1... ";
+ for(int var = 0; var < variable_table.size(); var++)
+ for(int eq = 0; eq < (int) equations.size(); eq++)
+ {
+ NodeID d1 = equations[eq]->getDerivative(var);
+ if (d1 == Zero)
+ continue;
+ first_derivatives[make_pair(eq, var)] = d1;
+ }
+ cout << "done" << endl;
+
+ if (order == 2)
+ {
+ cout << " Processing Order 2... ";
+ for(first_derivatives_type::const_iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int var1 = it->first.second;
+ NodeID d1 = it->second;
+
+ // Store only second derivatives with var2 <= var1
+ for(int var2 = 0; var2 <= var1; var2++)
+ {
+ NodeID d2 = d1->getDerivative(var2);
+ if (d2 == Zero)
+ continue;
+ second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
+ }
+ }
+ cout << "done" << endl;
+ }
+}
+
+void
+ModelTree::computeTemporaryTerms(int order)
+{
+ map<NodeID, int> reference_count;
+ temporary_terms.clear();
+
+ for(vector<BinaryOpNode *>::iterator it = equations.begin();
+ it != equations.end(); it++)
+ (*it)->computeTemporaryTerms(reference_count, temporary_terms);
+
+ for(first_derivatives_type::iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms);
+
+ if (order == 2)
+ for(second_derivatives_type::iterator it = second_derivatives.begin();
+ it != second_derivatives.end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms);
+}
+
+void
+ModelTree::writeTemporaryTerms(ostream &output, bool is_dynamic) const
+{
+ // A copy of temporary terms
+ temporary_terms_type tt2;
+
+ for(temporary_terms_type::const_iterator it = temporary_terms.begin();
+ it != temporary_terms.end(); it++)
+ {
+ (*it)->writeOutput(output, is_dynamic, temporary_terms);
+ output << " = ";
+
+ (*it)->writeOutput(output, is_dynamic, tt2);
+
+ // Insert current node into tt2
+ tt2.insert(*it);
+
+ output << ";" << endl;
+ }
+}
+
+void
+ModelTree::writeLocalParameters(ostream &output, bool is_dynamic) const
+{
+ for(map<int, NodeID>::const_iterator it = local_parameters_table.begin();
+ it != local_parameters_table.end(); it++)
+ {
+ int id = it->first;
+ NodeID value = it->second;
+ output << symbol_table.getNameByID(eLocalParameter, id) << " = ";
+ value->writeOutput(output, is_dynamic, temporary_terms);
+ output << ";" << endl;
+ }
+}
+
+void
+ModelTree::writeModelEquations(ostream &output, bool is_dynamic) const
+{
+ for(int eq = 0; eq < (int) equations.size(); eq++)
+ {
+ BinaryOpNode *eq_node = equations[eq];
+
+ NodeID lhs = eq_node->arg1;
+ output << "lhs =";
+ lhs->writeOutput(output, is_dynamic, temporary_terms);
+ output << ";" << endl;
+
+ NodeID rhs = eq_node->arg2;
+ output << "rhs =";
+ rhs->writeOutput(output, is_dynamic, temporary_terms);
+ output << ";" << endl;
+
+ output << "residual" << lpar << eq + 1 << rpar << "= lhs-rhs;" << endl;
+ }
+}
+
+void
+ModelTree::writeStaticMFile(const string &static_basename) const
{
string filename = static_basename + interfaces::function_file_extension();
@@ -74,7 +159,7 @@ ModelTree::writeStaticMFile(const string &static_basename)
void
-ModelTree::writeDynamicMFile(const string &dynamic_basename)
+ModelTree::writeDynamicMFile(const string &dynamic_basename) const
{
string filename = dynamic_basename + interfaces::function_file_extension();
@@ -100,7 +185,7 @@ ModelTree::writeDynamicMFile(const string &dynamic_basename)
}
void
-ModelTree::writeStaticCFile(const string &static_basename)
+ModelTree::writeStaticCFile(const string &static_basename) const
{
string filename = static_basename + ".c";
@@ -143,7 +228,7 @@ ModelTree::writeStaticCFile(const string &static_basename)
mStaticModelFile << " if (nlhs >= 1)\n";
mStaticModelFile << " {\n";
mStaticModelFile << " /* Set the output pointer to the output matrix residual. */\n";
- mStaticModelFile << " plhs[0] = mxCreateDoubleMatrix(" << eq_nbr << ",1, mxREAL);\n";
+ mStaticModelFile << " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);\n";
mStaticModelFile << " /* Create a C pointer to a copy of the output matrix residual. */\n";
mStaticModelFile << " residual = mxGetPr(plhs[0]);\n";
mStaticModelFile << " }\n\n";
@@ -151,7 +236,7 @@ ModelTree::writeStaticCFile(const string &static_basename)
mStaticModelFile << " if (nlhs >= 2)\n";
mStaticModelFile << " {\n";
mStaticModelFile << " /* Set the output pointer to the output matrix g1. */\n";
- mStaticModelFile << " plhs[1] = mxCreateDoubleMatrix(" << eq_nbr << ", " << symbol_table.endo_nbr << ", mxREAL);\n";
+ mStaticModelFile << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr << ", mxREAL);\n";
mStaticModelFile << " /* Create a C pointer to a copy of the output matrix g1. */\n";
mStaticModelFile << " g1 = mxGetPr(plhs[1]);\n";
mStaticModelFile << " }\n\n";
@@ -170,7 +255,7 @@ ModelTree::writeStaticCFile(const string &static_basename)
void
-ModelTree::writeDynamicCFile(const string &dynamic_basename)
+ModelTree::writeDynamicCFile(const string &dynamic_basename) const
{
string filename = dynamic_basename + ".c";
@@ -218,7 +303,7 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename)
mDynamicModelFile << " if (nlhs >= 1)\n";
mDynamicModelFile << " {\n";
mDynamicModelFile << " /* Set the output pointer to the output matrix residual. */\n";
- mDynamicModelFile << " plhs[0] = mxCreateDoubleMatrix(" << eq_nbr << ",1, mxREAL);\n";
+ mDynamicModelFile << " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);\n";
mDynamicModelFile << " /* Create a C pointer to a copy of the output matrix residual. */\n";
mDynamicModelFile << " residual = mxGetPr(plhs[0]);\n";
mDynamicModelFile << " }\n\n";
@@ -227,9 +312,9 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename)
mDynamicModelFile << " {\n";
mDynamicModelFile << " /* Set the output pointer to the output matrix g1. */\n";
if (computeJacobianExo)
- mDynamicModelFile << " plhs[1] = mxCreateDoubleMatrix(" << eq_nbr << ", " << variable_table.get_dyn_var_nbr() << ", mxREAL);\n";
+ mDynamicModelFile << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr() << ", mxREAL);\n";
else if (computeJacobian)
- mDynamicModelFile << " plhs[1] = mxCreateDoubleMatrix(" << eq_nbr << ", " << variable_table.var_endo_nbr << ", mxREAL);\n";
+ mDynamicModelFile << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.var_endo_nbr << ", mxREAL);\n";
mDynamicModelFile << " /* Create a C pointer to a copy of the output matrix g1. */\n";
mDynamicModelFile << " g1 = mxGetPr(plhs[1]);\n";
mDynamicModelFile << " }\n\n";
@@ -237,7 +322,7 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename)
mDynamicModelFile << " if (nlhs >= 3)\n";
mDynamicModelFile << " {\n";
mDynamicModelFile << " /* Set the output pointer to the output matrix g2. */\n";
- mDynamicModelFile << " plhs[2] = mxCreateDoubleMatrix(" << eq_nbr << ", " << variable_table.get_dyn_var_nbr()*variable_table.get_dyn_var_nbr() << ", mxREAL);\n";
+ mDynamicModelFile << " plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr()*variable_table.get_dyn_var_nbr() << ", mxREAL);\n";
mDynamicModelFile << " /* Create a C pointer to a copy of the output matrix g1. */\n";
mDynamicModelFile << " g2 = mxGetPr(plhs[2]);\n";
mDynamicModelFile << " }\n\n";
@@ -257,502 +342,70 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename)
mDynamicModelFile.close();
}
-
-void ModelTree::derive(int iOrder)
-{
- NodeID lToken; // To store current working token
- NodeID lD1, lD2; // To store derivative arguments of current argument
- NodeID lArg1, lArg2; // To store current arguments
- Type lType1; // Type of first argument
- NodeID t1,t11,t12,t13,
- t14, t15; // To store temoporary result arguments
- TreeIterator BeginIT; // Iterator of the 1st token to derive
- TreeIterator EndIT; // Iterator of the last token to derive
- TreeIterator currentIT; // Iterator (counter) for model tree loop
-
- vector<NodeID> EqualTokenIDs; // IDs of "equal token" in model Tree
- // Capturing equation IDs
- for (currentIT = BeginModel; currentIT != mModelTree.end(); currentIT++)
- {
- if ((*currentIT)->op_code == token::EQUAL)
- {
- EqualTokenIDs.push_back(*currentIT);
- // Equation is forced to be in Model Tree as referenced
- // This is useful to remove symmetric elements
- (*currentIT)->reference_count[0]++;
- }
- }
- //std::cout << "size " << EqualTokenIDs.size() << "\n";
- mDerivativeIndex.resize(iOrder);
- // Uncomment this to print model tree data
- /*
- //cout << "ModelTree==================================\n";
- for (currentIT = mModelTree.begin(); currentIT != mModelTree.end(); currentIT++)
- {
- lToken = *currentIT;
- int ID = lToken->idx;
- cout << ID << ":" << lToken << "->" << lToken->id1 << " " << lToken->type1 << " " <<
- lToken->id2 << " " << lToken->op_code << "\n";
- }
- */
-
- EndIT = mModelTree.begin();
- EndIT--;
- cout << "Processing derivation ...\n";
- // loop on order of derivation
- for(int Order = 1; Order <= iOrder; Order++)
- {
-
- cout << "\tProcessing Order " << Order << "... ";
- current_order = Order;
- BeginIT = EndIT;
- BeginIT++;
- EndIT = mModelTree.end();
- EndIT--;
- // Adding a reference counter for current order to tokens in mModelTree
- // and updating them
-
- for (TreeIterator it = mModelTree.begin(); it != mModelTree.end(); it++)
- {
- int s = (*it)->reference_count.size();
- for (int i = s; i <= current_order; i++)
- {
- int rc = (*it)->reference_count[i-1];
- (*it)->reference_count.push_back(rc);
- }
- }
-
- // Loop on tokens
- for (currentIT = BeginIT;; currentIT++)
- {
- lToken = *currentIT;
-
- // Loop on variables of derivation which may give non null result for this token
- for (vector<int>::iterator varIt = lToken->non_null_derivatives.begin();
- varIt != lToken->non_null_derivatives.end(); varIt++)
- {
- int var = *varIt;
- //cout << "Token " << (*currentIT)->idx << endl;
- if (accumulate((*currentIT)->reference_count.begin(), (*currentIT)->reference_count.end(),0) > 0)
- {
- lArg1 = lToken->id1;
- lArg2 = lToken->id2;
- lType1 = lToken->type1;
- lD1 = DeriveArgument(lArg1, lType1, var);
- lD2 = Zero;
- if (lArg2 != NullID)
- lD2 = DeriveArgument(lArg2, eTempResult, var);
- // Case where token is a terminal
- if (lToken->op_code == NoOpCode)
- {
- (*currentIT)->setDerivativeAddress(lD1, var);
- }
- else if (lD1 == Zero && lD2 == Zero)
- {
- (*currentIT)->setDerivativeAddress(Zero, var);
- }
- else
- {
- switch (lToken->op_code)
- {
- case token::UMINUS:
- t1 = AddUMinus(lD1);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::PLUS:
- t1 = AddPlus(lD1, lD2);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::MINUS:
- t1 = AddMinus(lD1, lD2);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::TIMES:
- t11 = AddTimes(lD1, lArg2);
- t12 = AddTimes(lD2, lArg1);
- t1 = AddPlus(t11, t12);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::DIVIDE:
- t11 = AddTimes(lD1, lArg2);
- t12 = AddTimes(lD2, lArg1);
- t13 = AddMinus(t11, t12);
- t14 = AddTimes(lArg2, lArg2);
- t1 = AddDivide(t13, t14);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::SQRT:
- t11 = AddPlus(*currentIT, *currentIT);
- t1 = AddDivide(lD1, t11);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::POWER:
- if (lD2 == Zero)
- {
- if (lD1 == Zero)
- t1 = Zero;
- else
- {
- t11 = AddMinus(lArg2, One);
- t12 = AddPower(lArg1, t11);
- t13 = AddTimes(lArg2, t12);
- t1 = AddTimes(lD1, t13);
- }
- }
- else
- {
- t11 = AddLog(lArg1);
- t12 = AddTimes(lD2, t11);
- t13 = AddTimes(lD1, lArg2);
- t14 = AddDivide(t13, lArg1);
- t15 = AddPlus(t12, t14);
- t1 = AddTimes(t15, *currentIT);
- }
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::EXP:
- t1 = AddTimes(lD1, *currentIT);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::LOG:
- t1 = AddDivide(lD1, lArg1);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::LOG10:
- t11 = AddExp(One);
- t12 = AddLog10(t11);
- t13 = AddDivide(lD1, lArg1);
- t1 = AddTimes(t12, t13);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::COS:
- t11 = AddSin(lArg1);
- t12 = AddUMinus(t11);
- t1 = AddTimes( lD1, t12);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::SIN:
- t11 = AddCos(lArg1);
- t1 = AddTimes(lD1,t11);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::TAN:
- t11 = AddTimes(*currentIT, *currentIT);
- t12 = AddPlus(t11, One);
- t1 = AddTimes(lD1, t12);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::ACOS:
- t11 = AddSin(*currentIT);
- t12 = AddDivide(lD1, t11);
- t1 = AddUMinus(t12);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::ASIN:
- t11 = AddCos(*currentIT);
- t1 = AddDivide(lD1, t11);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::ATAN:
- t11 = AddTimes(lArg1, lArg1);
- t12 = AddPlus(One, t11);
- t1 = AddDivide(lD1, t12);
- (*currentIT)->setDerivativeAddress(t1,var);
- break;
- case token::COSH:
- t11 = AddSinH(lArg1);
- t1 = AddTimes( lD1,t11);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::SINH:
- t11 = AddCosH(lArg1);
- t1 = AddTimes( lD1, t11);
- (*currentIT)->setDerivativeAddress(t1, var);
- break;
- case token::TANH:
- t11 = AddTimes(*currentIT, *currentIT);
- t12 = AddMinus(One, t11);
- t1 = AddTimes(lD1, t12);
- (*currentIT)->setDerivativeAddress(t1,var);
- break;
- case token::ACOSH:
- t11 = AddSinH(*currentIT);
- t1 = AddDivide(lD1, t11);
- (*currentIT)->setDerivativeAddress(t1,var);
- break;
- case token::ASINH:
- t11 = AddCosH(*currentIT);
- t1 = AddDivide(lD1, t11);
- (*currentIT)->setDerivativeAddress(t1,var);
- break;
- case token::ATANH:
- t11 = AddTimes(lArg1, lArg1);
- t12 = AddMinus(One, t11);
- t1 = AddTimes(lD1, t12);
- (*currentIT)->setDerivativeAddress(t1,var);
- break;
- case token::EQUAL:
- // Force the derivative to have zero on right hand side
- // (required for setStaticModel and setDynamicModel)
- if (lD1 == Zero && lD2 != Zero)
- {
- t11 = AddUMinus(lD2);
- t1 = AddEqual(t11, Zero);
- }
- else if (lD1 != Zero && lD2 == Zero)
- t1 = AddEqual(lD1, Zero);
- else
- {
- t11 = AddMinus(lD1, lD2);
- t1 = AddEqual(t11, Zero);
- }
- // The derivative is forced to be in Model Tree as referenced
- // This is useful to remove symmetric elements
- IncrementReferenceCount(t1);
- lToken->setDerivativeAddress(t1, var);
- break;
- }
- }
- }
- }
- if (currentIT == EndIT)
- break;
- }
-
- // Filling mDerivativeIndex
- // Loop on variables of derivation, skipping symmetric elements
- for (int var = 0; var < variable_table.size(); var++)
- {
- int starti = var*Order*(Order-1)*eq_nbr/2;
- for (unsigned int i = starti; i < EqualTokenIDs.size() ; i++ )
- {
- t1 = EqualTokenIDs[i]->getDerivativeAddress(var, *this);
- if (Order == 1)
- mDerivativeIndex[0].push_back(DerivativeIndex(t1, i-starti, var));
- else if (Order == 2)
- {
- int var1 = variable_table.getSortID(i/eq_nbr);
- int var2 = variable_table.getSortID(var);
- mDerivativeIndex[1].push_back(DerivativeIndex(t1,
- i-eq_nbr*(i/eq_nbr),
- var1*variable_table.get_dyn_var_nbr()+var2));
- }
- }
- }
-
- // Uncomment to debug : prints unreferenced tokens
- /*
- cout << "Order : " << Order << "\n";
- for (TokenCount = BeginModel; TokenCount < mModelTree.size() ; TokenCount++ )
- {
- if (accumulate(mModelTree[TokenCount].reference_count.begin(),mModelTree[TokenCount].reference_count.end(),0) == 0)
- cout << "\tNot referenced : token ID :" << TokenCount << endl;
- }
- */
- // Uncomment this to debug : mDerivative(1and2)Index data
- // before removing unreferenced tokens
- /*
- cout << "Contenence of mDerivative1Index\n";
- for (int i=0; i< mDerivativeIndex[0].size();i++)
- //if (mDerivativeIndex[0][i].token_id != 3)
- cout << "\t" << mDerivativeIndex[0][i].token_id << endl;
- cout << "Contenence of mDerivative2Index\n";
- for (int i=0; i< mDerivativeIndex[1].size();i++)
- //if (mDerivativeIndex[1][i].token_id != 3)
- cout << "\t" << mDerivativeIndex[1][i].token_id << endl;
- */
- //cout << "Removing unreferenced tokens range ids :" << CurrentID << " - " << mModelTree.size()-1 << endl;
- // Removing unreferenced tokens in last derivative
- // RemoveUnref(CurrentID, mModelTree.size()-1, Order);
- // Decrementing reference couter of unreferenced tokens in last derivative
- //DecrementUnref(CurrentID, mModelTree.size()-1, Order);
- /*
- cout << "Order : " << Order << "\n";
- for (TokenCount = BeginModel; TokenCount < mModelTree.size() ; TokenCount++ )
- {
- if (accumulate(mModelTree[TokenCount].reference_count.begin(),mModelTree[TokenCount].reference_count.end(),0) == 0)
- cout << "\tNot referenced : token ID :" << TokenCount << endl;
- }
- */
- EqualTokenIDs.clear();
- // Updating EqualTokenIDs
- for (unsigned int i=0; i< mDerivativeIndex[Order-1].size();i++)
- {
- EqualTokenIDs.push_back(mDerivativeIndex[Order-1][i].token_id);
- }
-
- // Uncomment this to debug : mDerivative(1and2)Index data
- // after removing unreferenced tokens
- /*
- cout << "Contenence of mDerivative1Index after removing\n";
- for (int i=0; i< mDerivativeIndex[0].size();i++)
- if (mDerivativeIndex[0][i].token_id != 3)
- cout << "\t" << mDerivativeIndex[0][i].token_id << endl;
- cout << "Contenence of mDerivative2Index after removing\n";
- for (int i=0; i< mDerivativeIndex[1].size();i++)
- //if (mDerivativeIndex[1][i].token_id != 3)
- cout << "\t" << mDerivativeIndex[1][i].token_id << endl;
- */
- cout << "done" << endl;
- }
-}
-
-//------------------------------------------------------------------------------
-inline NodeID ModelTree::DeriveArgument(NodeID iArg, Type iType, int iVarID)
-{
- switch(iType)
- {
- case eTempResult :
- return iArg->getDerivativeAddress(iVarID, *this);
- case eExogenous :
- case eExogenousDet :
- case eEndogenous :
- case eRecursiveVariable :
- if ((long int) iArg == iVarID)
- return One;
- else
- {
- return Zero;
- }
- break;
- case eNumericalConstant :
- case eParameter :
- return Zero;
- case eLocalParameter :
- return Zero;
- case eUNDEF :
- return NullID;
- default :
- cout << "ModelTree::DeriveArgument : Error: Unknown Type\n";
- exit(-1);
- };
-
-}
-
void
-ModelTree::writeStaticModel(ostream &StaticOutput)
+ModelTree::writeStaticModel(ostream &StaticOutput) const
{
ostringstream model_output; // Used for storing model equations
ostringstream jacobian_output; // Used for storing jacobian equations
ostringstream hessian_output;
ostringstream lsymetric; // For symmetric elements in hessian
- int d = current_order; // Minimum number of times a temparary expression apears in equations
- // Reference count of token "0=0" is set to 0
- // Not to be printed as a temp expression
- fill(ZeroEqZero->reference_count.begin(),
- ZeroEqZero->reference_count.end(),0);
- // Writing model Equations
- current_order = 1;
- TreeIterator tree_it = BeginModel;
- int lEquationNBR = 0;
- for (; tree_it != mModelTree.end(); tree_it++)
- {
- if ((*tree_it)->op_code == token::EQUAL || (*tree_it)->op_code == token::ASSIGN )
- {
- if ((*tree_it)->id1->type1 == eLocalParameter)
- {
- model_output << getExpression((*tree_it)->id1, eStaticEquations);
- model_output << " = ";
- model_output << getExpression((*tree_it)->id2, eStaticEquations) << ";" << endl;
- }
- else if (lEquationNBR < eq_nbr)
- {
- model_output << "lhs =";
- model_output << getExpression((*tree_it)->id1, eStaticEquations) << ";" << endl;
- model_output << "rhs =";
- model_output << getExpression((*tree_it)->id2, eStaticEquations) << ";" << endl;
- model_output << "residual" << lpar << lEquationNBR+1 << rpar << "= lhs-rhs;" << endl;
- lEquationNBR++;
- }
- else
- break;
- }
- else if ((*tree_it)->op_code != NoOpCode)
- {
- if (optimize(*tree_it))
- {
- model_output << "T" << (*tree_it)->idx << "=" << getExpression(*tree_it, eStaticEquations) << ";" << endl;
- (*tree_it)->tmp_status = 1;
- }
- else
- {
- (*tree_it)->tmp_status = 0;
- }
- }
- }
+ writeTemporaryTerms(model_output, false);
- for(; tree_it != mModelTree.end(); tree_it++)
- {
- if ((*tree_it)->op_code != NoOpCode
- && (*tree_it)->op_code != token::EQUAL
- && (*tree_it)->op_code != token::ASSIGN)
- {
- if (optimize(*tree_it) == 1)
- {
- jacobian_output << "T" << (*tree_it)->idx << "=" << getExpression(*tree_it, eStaticEquations) << ";" << endl;
- (*tree_it)->tmp_status = 1;
- }
- else
- {
- (*tree_it)->tmp_status = 0;
- }
- }
- }
+ writeLocalParameters(model_output, false);
+
+ writeModelEquations(model_output, false);
// Write Jacobian w.r. to endogenous only
- for (unsigned int i = 0; i < mDerivativeIndex[0].size(); i++)
+ for(first_derivatives_type::const_iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
{
- if (variable_table.getType(mDerivativeIndex[0][i].derivators) == eEndogenous)
+ int eq = it->first.first;
+ int var = it->first.second;
+ NodeID d1 = it->second;
+
+ if (variable_table.getType(var) == eEndogenous)
{
- NodeID startJacobian = mDerivativeIndex[0][i].token_id;
- if (startJacobian != ZeroEqZero)
- {
- string exp = getExpression(startJacobian->id1, eStaticDerivatives);
- ostringstream g1;
- g1 << " g1" << lpar << mDerivativeIndex[0][i].equation_id+1 << ", " <<
- variable_table.getSymbolID(mDerivativeIndex[0][i].derivators)+1 << rpar;
- jacobian_output << g1.str() << "=" << g1.str() << "+" << exp << ";\n";
- }
+ ostringstream g1;
+ g1 << " g1" << lpar << eq + 1 << ", " << variable_table.getSymbolID(var) + 1 << rpar;
+
+ jacobian_output << g1.str() << "=" << g1.str() << "+";
+ d1->writeOutput(jacobian_output, false, temporary_terms);
+ jacobian_output << ";" << endl;
}
}
// Write Hessian w.r. to endogenous only
if (computeStaticHessian)
- {
- for (unsigned int i = 0; i < mDerivativeIndex[1].size(); i++)
- {
- NodeID startHessian = mDerivativeIndex[1][i].token_id;
- if (startHessian != ZeroEqZero)
- {
- string exp = getExpression(startHessian->id1, eStaticDerivatives);
-
- int varID1 = mDerivativeIndex[1][i].derivators / variable_table.get_dyn_var_nbr();
- int varID2 = mDerivativeIndex[1][i].derivators - varID1 * variable_table.get_dyn_var_nbr();
-
- // Keep only derivatives w.r. to endogenous variables
- if (variable_table.getType(varID1) != eEndogenous
- || variable_table.getType(varID2) != eEndogenous)
- continue;
+ for(second_derivatives_type::const_iterator it = second_derivatives.begin();
+ it != second_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int var1 = it->first.second.first;
+ int var2 = it->first.second.second;
+ NodeID d2 = it->second;
+
+ // Keep only derivatives w.r. to endogenous variables
+ if (variable_table.getType(var1) == eEndogenous
+ && variable_table.getType(var2) == eEndogenous)
+ {
+ int id1 = variable_table.getSymbolID(var1);
+ int id2 = variable_table.getSymbolID(var2);
- int id1 = variable_table.getSymbolID(varID1);
- int id2 = variable_table.getSymbolID(varID2);
+ int col_nb = id1*symbol_table.endo_nbr+id2+1;
+ int col_nb_sym = id2*symbol_table.endo_nbr+id1+1;
- int col_nb = id1*symbol_table.endo_nbr+id2+1;
- int col_nb_sym = id2*symbol_table.endo_nbr+id1+1;
+ hessian_output << " g2" << lpar << eq+1 << ", " << col_nb << rpar << " = ";
+ d2->writeOutput(hessian_output, false, temporary_terms);
+ hessian_output << ";" << endl;
- hessian_output << " g2" << lpar << mDerivativeIndex[1][i].equation_id+1 << ", " <<
- col_nb << rpar << " = " << exp << ";\n";
- // Treating symetric elements
- if (varID1 != varID2)
- lsymetric << " g2" << lpar << mDerivativeIndex[1][i].equation_id+1 << ", " <<
- col_nb_sym << rpar << " = " <<
- "g2" << lpar << mDerivativeIndex[1][i].equation_id+1 << ", " <<
- col_nb << rpar << ";\n";
- }
+ // Treating symetric elements
+ if (var1 != var2)
+ lsymetric << " g2" << lpar << eq+1 << ", " << col_nb_sym << rpar << " = "
+ << "g2" << lpar << eq+1 << ", " << col_nb << rpar << ";" << endl;
+ }
- }
- }
+ }
// Writing ouputs
if (offset == 1)
@@ -760,7 +413,7 @@ ModelTree::writeStaticModel(ostream &StaticOutput)
StaticOutput << "global M_ \n";
StaticOutput << "if M_.param_nbr > 0\n params = M_.params;\nend\n";
- StaticOutput << " residual = zeros( " << eq_nbr << ", 1);\n";
+ StaticOutput << " residual = zeros( " << equations.size() << ", 1);\n";
StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
StaticOutput << "Model equations\n\t";
StaticOutput << interfaces::comment() + "\n\n";
@@ -770,7 +423,7 @@ ModelTree::writeStaticModel(ostream &StaticOutput)
StaticOutput << "end\n";
StaticOutput << "if nargout >= 2,\n";
StaticOutput << " g1 = " <<
- "zeros(" << eq_nbr << ", " <<
+ "zeros(" << equations.size() << ", " <<
symbol_table.endo_nbr << ");\n" ;
StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
StaticOutput << "Jacobian matrix\n\t";
@@ -786,7 +439,7 @@ ModelTree::writeStaticModel(ostream &StaticOutput)
// Writing initialization instruction for matrix g2
int ncols = symbol_table.endo_nbr * symbol_table.endo_nbr;
StaticOutput << " g2 = " <<
- "sparse([],[],[]," << eq_nbr << ", " << ncols << ", " <<
+ "sparse([],[],[]," << equations.size() << ", " << ncols << ", " <<
5*ncols << ");\n";
StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
StaticOutput << "Hessian matrix\n\t";
@@ -814,155 +467,69 @@ ModelTree::writeStaticModel(ostream &StaticOutput)
StaticOutput << " }\n";
StaticOutput << "}\n\n";
}
- current_order = d;
}
void
-ModelTree::writeDynamicModel(ostream &DynamicOutput)
+ModelTree::writeDynamicModel(ostream &DynamicOutput) const
{
ostringstream lsymetric; // Used when writing symetric elements in Hessian
ostringstream model_output; // Used for storing model equations
ostringstream jacobian_output; // Used for storing jacobian equations
ostringstream hessian_output; // Used for storing Hessian equations
- int d = current_order;
-
- // Reference count of token "0=0" is set to 0
- // Not to be printed as a temp expression
- fill(ZeroEqZero->reference_count.begin(),
- ZeroEqZero->reference_count.end(),0);
-
- // Getting equations from model tree
- // Starting from the end of equation
- // Searching for the next '=' operator,
- current_order = 1;
- int lEquationNBR = 0;
- cout << "\tequations .. ";
- TreeIterator tree_it = BeginModel;
- for (; tree_it != mModelTree.end(); tree_it++)
- {
- if ((*tree_it)->op_code == token::EQUAL || (*tree_it)->op_code == token::ASSIGN)
- {
- if ((*tree_it)->id1->type1 == eLocalParameter)
- {
- model_output << getExpression((*tree_it)->id1, eStaticEquations);
- model_output << " = ";
- model_output << getExpression((*tree_it)->id2, eStaticEquations) << ";" << endl;
- }
- else if (lEquationNBR < eq_nbr)
- {
- model_output << "lhs =";
- model_output << getExpression(((*tree_it)->id1), eDynamicEquations) << ";" << endl;
- model_output << "rhs =";
- model_output << getExpression(((*tree_it)->id2), eDynamicEquations) << ";" << endl;
- model_output << "residual" << lpar << lEquationNBR+1 << rpar << "= lhs-rhs;" << endl;
- lEquationNBR++;
- }
- else
- break;
- }
- else if ((*tree_it)->op_code != NoOpCode)
- {
- if (optimize(*tree_it))
- {
- model_output << "T" << (*tree_it)->idx << "=" << getExpression(*tree_it, eDynamicEquations) << ";" << endl;
- (*tree_it)->tmp_status = 1;
- }
- else
- {
- (*tree_it)->tmp_status = 0;
- }
- }
- }
+ writeTemporaryTerms(model_output, true);
- for(; tree_it != mModelTree.end(); tree_it++)
- {
- if ((*tree_it)->op_code != NoOpCode
- && (*tree_it)->op_code != token::EQUAL
- && (*tree_it)->op_code != token::ASSIGN)
- {
- if (optimize(*tree_it) == 1)
- {
- jacobian_output << "T" << (*tree_it)->idx << "=" << getExpression(*tree_it, eDynamicEquations) << ";" << endl;
- (*tree_it)->tmp_status = 1;
- }
- else
- {
- (*tree_it)->tmp_status = 0;
- }
- }
- }
- cout << "done \n";
- // Getting Jacobian from model tree
+ writeLocalParameters(model_output, true);
+
+ writeModelEquations(model_output, true);
+
+ // Writing Jacobian
if (computeJacobian || computeJacobianExo)
- {
- cout << "\tJacobian .. ";
- for(; tree_it != mModelTree.end(); tree_it++)
- {
- if ((*tree_it)->op_code != NoOpCode
- && (*tree_it)->op_code != token::EQUAL
- && (*tree_it)->op_code != token::ASSIGN)
- {
- if (optimize(*tree_it) == 1)
- {
- jacobian_output << "T" << (*tree_it)->idx << "=" << getExpression(*tree_it, eDynamicEquations) << ";" << endl;
- (*tree_it)->tmp_status = 1;
- }
- else
- {
- (*tree_it)->tmp_status = 0;
- }
- }
- }
+ for(first_derivatives_type::const_iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int var = it->first.second;
+ NodeID d1 = it->second;
- for (unsigned int i = 0; i < mDerivativeIndex[0].size(); i++)
- {
- if (computeJacobianExo || variable_table.getType(mDerivativeIndex[0][i].derivators) == eEndogenous)
- {
- NodeID startJacobian = mDerivativeIndex[0][i].token_id;
- if (startJacobian != ZeroEqZero)
- {
- string exp = getExpression(startJacobian->id1, eDynamicDerivatives);
- ostringstream g1;
- g1 << " g1" << lpar << mDerivativeIndex[0][i].equation_id+1 << ", " <<
- variable_table.getSortID(mDerivativeIndex[0][i].derivators)+1 << rpar;
- jacobian_output << g1.str() << "=" << g1.str() << "+" << exp << ";\n";
- }
- }
- }
- cout << "done \n";
- }
+ if (computeJacobianExo || variable_table.getType(var) == eEndogenous)
+ {
+ ostringstream g1;
+ g1 << " g1" << lpar << eq + 1 << ", " << variable_table.getSortID(var) + 1 << rpar;
+
+ jacobian_output << g1.str() << "=" << g1.str() << "+";
+ d1->writeOutput(jacobian_output, true, temporary_terms);
+ jacobian_output << ";" << endl;
+ }
+ }
+ // Writing Hessian
if (computeHessian)
- {
- // Getting Hessian from model tree
- // Starting from the end of equation
- // Searching for the next '=' operator,
- cout << "\tHessian .. ";
- for (unsigned int i = 0; i < mDerivativeIndex[1].size(); i++)
- {
- NodeID startHessian = mDerivativeIndex[1][i].token_id;
- //cout << "ID = " << startHessian << " exp = " << exp << "\n";
- if (startHessian != ZeroEqZero)
- {
- string exp = getExpression(startHessian->id1, eDynamicDerivatives);
-
- int varID1 = mDerivativeIndex[1][i].derivators / variable_table.get_dyn_var_nbr();
- int varID2 = mDerivativeIndex[1][i].derivators - varID1 * variable_table.get_dyn_var_nbr();
- hessian_output << " g2" << lpar << mDerivativeIndex[1][i].equation_id+1 << ", " <<
- mDerivativeIndex[1][i].derivators+1 << rpar << " = " << exp << ";\n";
- // Treating symetric elements
- if (varID1 != varID2)
- lsymetric << " g2" << lpar << mDerivativeIndex[1][i].equation_id+1 << ", " <<
- varID2*variable_table.get_dyn_var_nbr()+varID1+1 << rpar << " = " <<
- "g2" << lpar << mDerivativeIndex[1][i].equation_id+1 << ", " <<
- mDerivativeIndex[1][i].derivators+1 << rpar << ";\n";
- }
+ for(second_derivatives_type::const_iterator it = second_derivatives.begin();
+ it != second_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int var1 = it->first.second.first;
+ int var2 = it->first.second.second;
+ NodeID d2 = it->second;
- }
- cout << "done \n";
- }
- int nrows = eq_nbr;
+ int id1 = variable_table.getSortID(var1);
+ int id2 = variable_table.getSortID(var2);
+
+ int col_nb = id1*variable_table.get_dyn_var_nbr()+id2+1;
+ int col_nb_sym = id2*variable_table.get_dyn_var_nbr()+id1+1;
+
+ hessian_output << " g2" << lpar << eq+1 << ", " << col_nb << rpar << " = ";
+ d2->writeOutput(hessian_output, true, temporary_terms);
+ hessian_output << ";" << endl;
+
+ // Treating symetric elements
+ if (id1 != id2)
+ lsymetric << " g2" << lpar << eq+1 << ", " << col_nb_sym << rpar << " = "
+ << "g2" << lpar << eq+1 << ", " << col_nb << rpar << ";" << endl;
+ }
+
+ int nrows = equations.size();
int nvars = variable_table.var_endo_nbr;
if (computeJacobianExo)
nvars += symbol_table.exo_nbr + symbol_table.exo_det_nbr;
@@ -1030,225 +597,6 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput)
}
DynamicOutput << "}\n\n";
}
- current_order = d;
-}
-
-string
-ModelTree::getExpression(NodeID StartID, EquationType iEquationType)
-{
-
- // Stack of tokens
- stack <int, vector<NodeID> > stack_token;
- ostringstream exp;
- NodeID current_token_ID;
-
- stack_token.push(StartID);
- int precedence_last_op = 0;
- int last_op_code = 0;
- int on_the_right_of_upper_node = 0;
-
- while (stack_token.size() > 0)
- {
- current_token_ID = stack_token.top();
- // defining short hand for current op code
- int current_op_code = current_token_ID->op_code;
- if ( current_token_ID->tmp_status == 1)
- {
- exp << "T" << current_token_ID->idx;
- // set precedence of terminal token to highest value
- precedence_last_op = 100;
- stack_token.pop();
- continue;
- }
- // else if current token is final
- else if ( current_op_code == NoOpCode)
- {
- exp << getArgument(current_token_ID->id1, current_token_ID->type1, iEquationType);
- // set precedence of terminal token to highest value
- precedence_last_op = 100;
- stack_token.pop();
- continue;
- }
-
- int precedence_current_op = OperatorTable::precedence(current_op_code);
- // deal with left argument first
- if (current_token_ID->left_done == 0)
- {
- // if current operator is not a temporary variable and
- // of lesser precedence than previous one, insert '('
- if (precedence_current_op < precedence_last_op
- || (on_the_right_of_upper_node == 1
- && (last_op_code == token::MINUS || last_op_code == token::DIVIDE)
- && (precedence_current_op == precedence_last_op))
- || (last_op_code == token::POWER && current_op_code == token::POWER)
- || current_op_code == token::UMINUS)
- {
- exp << "(";
- current_token_ID->close_parenthesis = 1;
- }
- // set flag: left argument has been explored
- current_token_ID->left_done = 1;
- precedence_last_op = precedence_current_op;
- last_op_code = current_op_code;
- if (offset == 0 && current_op_code == token::POWER)
- {
- exp << "pow(";
- precedence_last_op = 0;
- current_token_ID->close_parenthesis = 1;
- }
- else if (current_op_code == token::UMINUS)
- {
- exp << "-";
- current_token_ID->close_parenthesis = 1;
- }
- else if ( OperatorTable::isfunction(current_op_code) == true)
- {
- exp << current_token_ID->op_name << "(";
- precedence_last_op = 0;
- current_token_ID->close_parenthesis = 1;
- }
- on_the_right_of_upper_node = 0;
- stack_token.push(current_token_ID->id1);
- }
- // deal with right argument when left branch is entirely explored
- else if ( current_token_ID->right_done == 0 )
- {
- current_token_ID->right_done = 1;
- if (offset == 0 && current_op_code == token::POWER)
- {
- exp << ",";
- }
-
- if ( current_token_ID->id2 != NullID )
- {
- exp << current_token_ID->op_name;
- precedence_last_op = precedence_current_op;
- last_op_code = current_op_code;
- on_the_right_of_upper_node = 1;
- stack_token.push(current_token_ID->id2);
- }
- }
- else
- {
- if ( current_token_ID->close_parenthesis == 1)
- {
- exp << ")";
- }
- precedence_last_op = precedence_current_op;
- last_op_code = current_op_code;
- current_token_ID->left_done=0;
- current_token_ID->right_done=0;
- current_token_ID->close_parenthesis=0;
- stack_token.pop();
- }
- }
- return exp.str();
-}
-
-//------------------------------------------------------------------------------
-inline string ModelTree::getArgument(NodeID id, Type type, EquationType iEquationType)
-{
-
- stringstream argument;
-
- if (type == eParameter)
- {
- argument << "params" << lpar << (long int)id+offset << rpar;
- }
- else if (type == eLocalParameter)
- {
- argument << symbol_table.getNameByID(eLocalParameter, (long int) id);
- }
- else if (type == eNumericalConstant)
- {
- argument << num_constants.get((long int) id);
- }
- else if (type == eEndogenous || type == eExogenous || type == eExogenousDet)
- if (iEquationType == eStaticEquations || iEquationType == eStaticDerivatives)
- {
- int idx = variable_table.getSymbolID((long int) id)+offset;
- if (type == eEndogenous)
- {
- argument << "y" << lpar << idx << rpar;
- }
- else if (type == eExogenous)
- {
- argument << "x" << lpar << idx << rpar;
- }
- else if (type == eExogenousDet)
- {
- idx += symbol_table.exo_nbr;
- argument << "x" << lpar << idx << rpar;
- }
- }
- else
- {
- if (type == eEndogenous)
- {
- int idx = variable_table.getPrintIndex((long int) id) + offset;
- argument << "y" << lpar << idx << rpar;
- }
- else if (type == eExogenous)
- {
- int idx = variable_table.getSymbolID((long int) id) + offset;
- int lag = variable_table.getLag((long int) id);
- if (offset == 1)
- {
- if ( lag != 0)
- {
- argument << "x" << lpar << "it_ + " << lag
- << ", " << idx << rpar;
- }
- else
- {
- argument << "x" << lpar << "it_, " << idx << rpar;
- }
- }
- else
- {
- if ( lag != 0)
- {
- argument << "x" << lpar << "it_+" << lag
- << "+" << idx << "*nb_row_x" << rpar;
- }
- else
- {
- argument << "x" << lpar << "it_+" << idx << "*nb_row_x" << rpar;
- }
- }
- }
- else if (type == eExogenousDet)
- {
- int idx = variable_table.getSymbolID((long int) id) + symbol_table.exo_nbr + offset;
- int lag = variable_table.getLag((long int) id);
- if (offset == 1)
- {
- if (lag != 0)
- {
- argument << "x" << lpar << "it_ + " << lag
- << ", " << idx << rpar;
- }
- else
- {
- argument << "x" << lpar << "it_, " << idx << rpar;
- }
- }
- else
- {
- if (lag != 0)
- {
- argument << "x" << lpar << "it_ + " << lag
- << "+" << idx << "*nb_row_xd" << rpar;
- }
- else
- {
- argument << "x" << lpar << "it_+" << idx << "*nb_row_xd" << rpar;
- }
- }
-
- }
- }
- return argument.str();
}
void
@@ -1315,42 +663,25 @@ ModelTree::writeOutput(ostream &output) const
output << "M_.params = zeros(" << symbol_table.parameter_nbr << ", 1);\n";
}
-inline int ModelTree::optimize(NodeID node)
+void
+ModelTree::addEquation(NodeID eq)
{
- int cost;
- int tmp_status = 0;
- if (node->op_code != NoOpCode)
- {
- cost = OperatorTable::cost(node->op_code,offset);
- if (node->id1 != NullID && node->id1->op_code != NoOpCode)
- {
- cost += OperatorTable::cost(node->id1->op_code,offset);
- }
- if (node->id2 != NullID && node->id2->op_code != NoOpCode)
- {
- cost += OperatorTable::cost(node->id2->op_code,offset);
- }
- }
- cost *= node->reference_count[current_order];
- if (cost > min_cost)
- {
- tmp_status = 1;
- node->cost = 0;
- }
- else
+ BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
+
+ if (beq == NULL || beq->op_code != oEqual)
{
- tmp_status = 0;
- node->cost = cost;
+ cerr << "ModelTree::addEquation: you didn't provide an equal node!" << endl;
+ exit(-1);
}
- node->tmp_status = 0;
- return tmp_status;
+
+ equations.push_back(beq);
}
void
ModelTree::checkPass() const
{
// Exit if there is no equation in model file
- if (eq_nbr == 0)
+ if (equations.size() == 0)
{
cerr << "No equation found in model file" << endl;
exit(-1);
@@ -1360,34 +691,38 @@ ModelTree::checkPass() const
void
ModelTree::computingPass()
{
- cout << eq_nbr << " equation(s) found \n";
+ cout << equations.size() << " equation(s) found" << endl;
+
// Sorting variable table
variable_table.Sort();
- // Setting number of equations in ModelParameters class
- // Here no derivative are computed
- BeginModel++;
- min_cost = 40 * OperatorTable::cost(token::PLUS, offset);
- // Setting format of parentheses
if (offset == 1)
{
+ min_cost = 40 * 90;
lpar = '(';
rpar = ')';
}
else
{
+ min_cost = 40 * 4;
lpar = '[';
rpar = ']';
}
if (computeHessian || computeStaticHessian)
- derive(2);
+ {
+ derive(2);
+ computeTemporaryTerms(2);
+ }
else
- derive(1);
+ {
+ derive(1);
+ computeTemporaryTerms(1);
+ }
}
void
-ModelTree::writeStaticFile(const string &basename)
+ModelTree::writeStaticFile(const string &basename) const
{
if (offset)
writeStaticMFile(basename + "_static");
@@ -1396,7 +731,7 @@ ModelTree::writeStaticFile(const string &basename)
}
void
-ModelTree::writeDynamicFile(const string &basename)
+ModelTree::writeDynamicFile(const string &basename) const
{
if (offset)
writeDynamicMFile(basename + "_dynamic");
diff --git a/parser.src/ParsingDriver.cc b/parser.src/ParsingDriver.cc
index cb592d7c0b6ad173446e499feba92cdd9cdb83dd..d67f779d1cb12b9fa40ca04c48badf9adbc0dae3 100644
--- a/parser.src/ParsingDriver.cc
+++ b/parser.src/ParsingDriver.cc
@@ -42,7 +42,6 @@ ParsingDriver::parse(const string &f)
{
mod_file = new ModFile();
- mod_file->model_tree.error = error;
mod_file->symbol_table.error = error;
expression.setNumericalConstants(&mod_file->num_constants);
@@ -118,16 +117,16 @@ ParsingDriver::add_expression_constant(string *constant)
NodeID
ParsingDriver::add_model_constant(string *constant)
{
- int id = mod_file->num_constants.AddConstant(*constant);
+ NodeID id = model_tree->AddNumConstant(*constant);
delete constant;
- return model_tree->AddTerminal((NodeID) id, eNumericalConstant);
+ return id;
}
NodeID
ParsingDriver::add_model_variable(string *name)
{
check_symbol_existence(*name);
- NodeID id = model_tree->AddTerminal(*name);
+ NodeID id = model_tree->AddVariable(*name);
delete name;
return id;
}
@@ -145,7 +144,7 @@ ParsingDriver::add_model_variable(string *name, string *olag)
<< *name
<< " has lag " << lag << endl;
}
- NodeID id = model_tree->AddTerminal(*name, lag);
+ NodeID id = model_tree->AddVariable(*name, lag);
delete name;
delete olag;
return id;
@@ -920,8 +919,8 @@ void
ParsingDriver::end_planner_objective(NodeID expr)
{
// Add equation corresponding to expression
- model_tree->AddEqual(expr, model_tree->Zero);
- model_tree->eq_nbr++;
+ NodeID eq = model_tree->AddEqual(expr, model_tree->Zero);
+ model_tree->addEquation(eq);
mod_file->addStatement(new PlannerObjectiveStatement(model_tree));
}
@@ -938,7 +937,7 @@ NodeID
ParsingDriver::add_model_equal(NodeID arg1, NodeID arg2)
{
NodeID id = model_tree->AddEqual(arg1, arg2);
- model_tree->eq_nbr++;
+ model_tree->addEquation(id);
return id;
}
@@ -952,8 +951,14 @@ void
ParsingDriver::declare_and_init_local_parameter(string *name, NodeID rhs)
{
mod_file->symbol_table.AddSymbolDeclar(*name, eLocalParameter, *name);
- NodeID id = model_tree->AddTerminal(*name);
- model_tree->AddAssign(id, rhs);
+ try
+ {
+ model_tree->AddLocalParameter(*name, rhs);
+ }
+ catch(DataTree::LocalParameterException &e)
+ {
+ error("Local parameter " + e.name + " declared twice");
+ }
delete name;
}
diff --git a/parser.src/SymbolTable.cc b/parser.src/SymbolTable.cc
index c710b94d9ffcea96bd779765c6a7b0506bf7d055..cd4b186045a098a8c9daf6c004fa12e98e69ad64 100644
--- a/parser.src/SymbolTable.cc
+++ b/parser.src/SymbolTable.cc
@@ -100,7 +100,7 @@ Reference SymbolTable::isReferenced(const std::string &name) const
}
void
-SymbolTable::writeOutput(ostream &output)
+SymbolTable::writeOutput(ostream &output) const
{
if (exo_nbr > 0)
{
diff --git a/parser.src/include/DataTree.hh b/parser.src/include/DataTree.hh
index 0e2cf50bcf0b581eba14d1e3db6163f7a031b6fa..50d19e7861fccc9e2d2632e25435bc5475573b30 100644
--- a/parser.src/include/DataTree.hh
+++ b/parser.src/include/DataTree.hh
@@ -1,778 +1,144 @@
#ifndef _DATATREE_HH
#define _DATATREE_HH
-//------------------------------------------------------------------------------
-/*! \file
- \version 1.0
- \date 04/09/2004
- \par This file defines the DataTree class.
-*/
-//------------------------------------------------------------------------------
-#include <iostream>
+
+using namespace std;
+
#include <string>
-#include <list>
-#include <stack>
-#include <sstream>
-#include <fstream>
#include <map>
-#include <stdio.h>
-//------------------------------------------------------------------------------
+#include <list>
+
#include "SymbolTable.hh"
-#include "OperatorTable.hh"
#include "NumericalConstants.hh"
-#include "ModelTypes.hh"
#include "VariableTable.hh"
+#include "ExprNode.hh"
-typedef std::map<MToken, NodeID, MTokenLess> TreeMap;
-typedef std::list<NodeID> TreeList;
-typedef TreeList::iterator TreeIterator;
-/*!
- \class DataTree
- \brief Provides data storage for model tree
-*/
class DataTree
{
-protected :
+ friend class ExprNode;
+ friend class NumConstNode;
+ friend class VariableNode;
+ friend class UnaryOpNode;
+ friend class BinaryOpNode;
+protected:
//! A reference to the symbol table
SymbolTable &symbol_table;
+ //! Reference to numerical constants table
+ NumericalConstants &num_constants;
//! The variable table
VariableTable variable_table;
-
- /*! A list of structures "token" */
- TreeList mModelTree;
- /*! matches key with entry id in list of token */
- TreeMap mIndexOfTokens;
- /*! A counter for filling MetatToken's idx field */
- int nodeCounter;
- /*! ID of first token in model tree (first tokens are "0", "1" and "0=" */
- TreeIterator BeginModel;
- /*!
- Used with field reference_count to know if
- dealing with derivative 0, 1, 2 etc
- */
- int current_order;
- /*! Pushs token into model tree */
- inline NodeID PushToken(NodeID iArg1,int iOpCode, NodeID iArg2 = NULL, Type iType1 = eTempResult);
- inline NodeID getIDOfToken(const MToken &iToken);
-public :
- /*! Flag for empty operator (final toekn) */
- const int NoOpCode;
- const NodeID NullID;
- NodeID Zero;
- NodeID One;
- NodeID MinusOne;
- NodeID ZeroEqZero;
- /*! Type of output 0 for C and 1 for Matlab (default) , also used as matrix index offset*/
+
+ typedef list<NodeID> node_list_type;
+ //! The list of nodes
+ node_list_type node_list;
+ //! A counter for filling ExprNode's idx field
+ int node_counter;
+
+ //! Stores local parameters value
+ map<int, NodeID> local_parameters_table;
+
+ //! Computing cost above which a node can be declared a temporary term
+ int min_cost;
+
+ //! Left indexing parenthesis
+ char lpar;
+ //! Right indexing parenthesis
+ char rpar;
+
+ typedef map<int, NodeID> num_const_node_map_type;
+ num_const_node_map_type num_const_node_map;
+ typedef map<pair<int, Type>, NodeID> variable_node_map_type;
+ variable_node_map_type variable_node_map;
+ typedef map<pair<NodeID, int>, NodeID> unary_op_node_map_type;
+ unary_op_node_map_type unary_op_node_map;
+ typedef map<pair<pair<NodeID, NodeID>, int>, NodeID> binary_op_node_map_type;
+ binary_op_node_map_type binary_op_node_map;
+
+ inline NodeID AddUnaryOp(UnaryOpcode op_code, NodeID arg);
+ inline NodeID AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2);
+public:
+ DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg);
+ virtual ~DataTree();
+ NodeID Zero, One, MinusOne;
+ //! Type of output 0 for C and 1 for Matlab (default), also used as matrix index offset
int offset;
- /*! Pointer to error function of parser class */
- void (* error) (const char* m);
- /*! Increment reference count of given token */
- inline void IncrementReferenceCount(NodeID token);
- /*! Adds terminal token to model tree */
- inline NodeID AddTerminal(NodeID iArg, Type type);
- /*! Adds terminal token (a variable) to model tree */
- inline NodeID AddTerminal(std::string iArgName, int iLag = 0);
- /*! Adds "arg1+arg2" to model tree */
- inline NodeID AddPlus(NodeID iArg1, NodeID iArg2);
- /*! Adds "arg1-arg2" to model tree */
- inline NodeID AddMinus(NodeID iArg1, NodeID iArg2);
- /*! Adds "-arg" to model tree */
- inline NodeID AddUMinus(NodeID iArg1);
- /*! Adds "arg1*arg2" to model tree */
- inline NodeID AddTimes(NodeID iArg1, NodeID iArg2);
- /*! Adds "arg1/arg2" to model tree */
- inline NodeID AddDivide(NodeID iArg1, NodeID iArg2);
- /*! Adds "arg1^arg2" to model tree */
- inline NodeID AddPower(NodeID iArg1, NodeID iArg2);
- /*! Adds "exp(arg)" to model tree */
- inline NodeID AddExp(NodeID iArg1);
- /*! Adds "log(arg)" to model tree */
- inline NodeID AddLog(NodeID iArg1);
- /*! Adds "log10(arg)" to model tree */
- inline NodeID AddLog10(NodeID iArg1);
- /*! Adds "cos(arg)" to model tree */
- inline NodeID AddCos(NodeID iArg1);
- /*! Adds "sin(arg)" to model tree */
- inline NodeID AddSin(NodeID iArg1);
- /*! Adds "tan(arg)" to model tree */
- inline NodeID AddTan(NodeID iArg1);
- /*! Adds "acos(arg)" to model tree */
- inline NodeID AddACos(NodeID iArg1);
- /*! Adds "asin(arg)" to model tree */
- inline NodeID AddASin(NodeID iArg1);
- /*! Adds "atan(arg)" to model tree */
- inline NodeID AddATan(NodeID iArg1);
- /*! Adds "cosh(arg)" to model tree */
- inline NodeID AddCosH(NodeID iArg1);
- /*! Adds "sinh(arg)" to model tree */
- inline NodeID AddSinH(NodeID iArg1);
- /*! Adds "tanh(arg)" to model tree */
- inline NodeID AddTanH(NodeID iArg1);
- /*! Adds "acosh(arg)" to model tree */
- inline NodeID AddACosH(NodeID iArg1);
- /*! Adds "asinh(arg)" to model tree */
- inline NodeID AddASinH(NodeID iArg1);
- /*! Adds "atanh(args)" to model tree */
- inline NodeID AddATanH(NodeID iArg1);
- /*! Adds "sqrt(arg)" to model tree */
- inline NodeID AddSqRt(NodeID iArg1);
- /*! Adds "arg1=arg2" to model tree */
- inline NodeID AddEqual(NodeID iArg1, NodeID iArg2);
- /*! Adds "arg1=arg2" as assignment to model tree */
- inline NodeID AddAssign(NodeID iArg1, NodeID iArg2);
-public :
- /*! Constructor */
- DataTree(SymbolTable &symbol_table_arg);
- /*! Destructor */
- ~DataTree();
-};
-//------------------------------------------------------------------------------
-inline NodeID DataTree::PushToken(NodeID iArg1,int iOpCode, NodeID iArg2, Type iType1)
-{
- MetaToken* lToken = new MetaToken(iArg1, iType1, iArg2, iOpCode);
-
- if (iOpCode != NoOpCode)
- {
- lToken->op_name = OperatorTable::str(iOpCode);
- }
- else
- {
- lToken->op_name = "";
- }
- lToken->reference_count.resize(current_order+1,0);
- lToken->idx = nodeCounter++;
- mModelTree.push_back(lToken);
- //Updating reference counters and time costs
- if (iType1 == eTempResult )
- {
- lToken->cost = OperatorTable::cost(iOpCode,offset)+iArg1->cost;
- IncrementReferenceCount(iArg1);
- }
- if(iArg2 != NullID)
- {
- lToken->cost += iArg2->cost;
- IncrementReferenceCount(iArg2);
- }
- mIndexOfTokens[*lToken] = lToken;
-
- /*
- std::cout << "ID = " << ID << " / " << mIndexOfTokens.size()-1<< " - " << getIDOfToken(lToken2) << " " <<
- lToken.op_code << " " << lToken.id1 << " " <<
- lToken.id2 << " " << lToken.type1<< "\n";
- */
- /*
- std::cout << "\nmIndexOfTokens------------\n";
- for (iter=mIndexOfTokens.begin();iter != mIndexOfTokens.end(); iter++)
- std::cout << (*iter).second << " -> " << (*iter).first << std::endl;
- */
- return lToken;
-}
-
-//------------------------------------------------------------------------------
-inline void DataTree::IncrementReferenceCount(NodeID token)
-{
- int s = token->reference_count.size();
- for (int i = s; i <= current_order; i++)
- {
- int rc = token->reference_count[i-1];
- token->reference_count.push_back(rc);
- }
- token->reference_count[current_order]++;
-
-}
-
-//------------------------------------------------------------------------------
-inline NodeID DataTree::getIDOfToken(const MToken &iToken)
-{
- TreeMap::iterator iter = mIndexOfTokens.find(iToken);
- if (iter != mIndexOfTokens.end())
- return iter->second;
- else
- return NullID;
-}
-
-//------------------------------------------------------------------------------
-inline NodeID DataTree::AddTerminal(NodeID iArg1, Type iType1)
-{
-
- if (iType1 == eTempResult)
- {
- return iArg1;
- }
- else
- {
- MToken lToken(iArg1, iType1, NullID, NoOpCode);
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, NoOpCode, NullID, iType1);
- }
-
-}
-
-inline NodeID DataTree::AddTerminal(std::string iArgName, int iLag)
-{
- int id;
- Type type;
-
- if (!symbol_table.Exist(iArgName))
- {
- std::string msg = "unknown symbol: " + iArgName;
- (* error) (msg.c_str());
- exit(-1);
- }
- type = symbol_table.getType(iArgName);
- if (type != eUNDEF)
- {
- symbol_table.SetReferenced(iArgName);
- if (type == eEndogenous ||
- type == eExogenousDet ||
- type == eExogenous ||
- type == eRecursiveVariable)
- id = variable_table.AddVariable(iArgName,iLag);
- else
- id = symbol_table.getID(iArgName);
-
- }
- else
- {
- std::string msg = "unknown parameter: " + iArgName;
- (* error) (msg.c_str());
- exit(-1);
- }
- return AddTerminal((NodeID) id, type);
-}
-
-inline NodeID DataTree::AddPlus(NodeID iArg1, NodeID iArg2)
-{
-
- if (iArg1 != Zero && iArg2 != Zero)
- {
- MToken lToken;
-
- if (iArg1 <= iArg2)
- lToken = MToken(iArg1, eTempResult, iArg2, token::PLUS);
- else
- lToken = MToken(iArg2, eTempResult, iArg1, token::PLUS);
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
-
- // To treat commutativity of "+"
- // Token id (iArg1 and iArg2) are sorted
- if (iArg1 <= iArg2)
- return PushToken(iArg1, token::PLUS, iArg2);
- else
- return PushToken(iArg2, token::PLUS, iArg1);
- }
- else if (iArg1 != Zero)
- {
- return iArg1;
- }
- else if (iArg2 != Zero)
- {
- return iArg2;
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddMinus(NodeID iArg1, NodeID iArg2)
-{
-
- if (iArg1 != Zero && iArg2 != Zero)
- {
- MToken lToken(iArg1, eTempResult, iArg2, token::MINUS);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::MINUS, iArg2);
- }
- else if (iArg1 != Zero)
- {
- return iArg1;
- }
- else if (iArg2 != Zero)
- {
- return AddUMinus(iArg2);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddUMinus(NodeID iArg1)
-{
- if (iArg1 != Zero )
- {
- MToken lToken(iArg1, eTempResult, NullID, token::UMINUS);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- if (iArg1->type1 == eTempResult &&
- iArg1->op_code == token::UMINUS)
- {
- IncrementReferenceCount(iArg1->id1);
- return iArg1->id1;
-
- }
- return PushToken(iArg1, token::UMINUS);
- }
- else
- {
- return Zero;
- }
-
-}
-
-inline NodeID DataTree::AddTimes(NodeID iArg1, NodeID iArg2)
-{
-
- if (iArg1 != Zero && iArg1 != One && iArg2 != Zero && iArg2 != One)
- {
- MToken lToken;
-
- if (iArg1 <= iArg2)
- lToken = MToken(iArg1, eTempResult, iArg2, token::TIMES);
- else
- lToken = MToken(iArg2, eTempResult, iArg1, token::TIMES);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- // Testing if arg[1 or 2] is "-1", in this case add "-arg" instead of "-1*iarg"
- // Here 2 is ID of "-1"
- if (iArg1 == MinusOne)
- return AddUMinus(iArg2);
- else if (iArg2 == MinusOne)
- return AddUMinus(iArg1);
- else
- {
- // To treat commutativity of "*"
- // Token id (iArg1 and iArg2) are sorted
- if (iArg1 <= iArg2)
- return PushToken(iArg1, token::TIMES, iArg2);
- else
- return PushToken(iArg2, token::TIMES, iArg1);
- }
- }
- else if (iArg1 != Zero && iArg1 != One && iArg2 == One)
- {
- return iArg1;
- }
- else if (iArg2 != Zero && iArg2 != One && iArg1 == One)
- {
- return iArg2;
- }
- else if (iArg2 == One && iArg1 == One)
- {
- return One;
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddDivide(NodeID iArg1, NodeID iArg2)
-{
-
- if (iArg1 != Zero && iArg2 != Zero && iArg2 != One)
- {
- MToken lToken(iArg1, eTempResult, iArg2, token::DIVIDE);
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::DIVIDE, iArg2);
- }
- else if (iArg2 == One)
- {
- return iArg1;
- }
- else if (iArg1 == Zero && iArg2 != Zero)
- {
- return Zero;
- }
- else
- {
- std::cout << "DIVIDE 0/0 non available\n";
- exit(-1);
- }
-}
-
-inline NodeID DataTree::AddPower(NodeID iArg1, NodeID iArg2)
-{
- if (iArg1 != Zero && iArg2 != Zero && iArg2 != One)
- {
- MToken lToken(iArg1, eTempResult, iArg2, token::POWER);
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::POWER, iArg2);
- }
- else if (iArg2 == One)
- {
- return iArg1;
- }
- else if (iArg2 == Zero)
- {
- return One;
- }
- else
- {
- return Zero;
- }
-
-}
-
-inline NodeID DataTree::AddExp(NodeID iArg1)
-{
-
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::EXP);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
-
- return PushToken(iArg1, token::EXP);
- }
- else
- {
- return One;
- }
-}
-inline NodeID DataTree::AddLog(NodeID iArg1)
-{
- if (iArg1 != Zero && iArg1 != One)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::LOG);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::LOG);
- }
- else if (iArg1 == One)
- {
- return Zero;
- }
- else
- {
- std::cout << "log(0) isn't available\n";
- exit(-1);
- }
-}
-
-inline NodeID DataTree::AddLog10(NodeID iArg1)
-{
- if (iArg1 != Zero && iArg1 != One)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::LOG10);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
-
- return PushToken(iArg1, token::LOG10);
- }
- else if (iArg1 == One)
- {
- return Zero;
- }
- else
- {
- std::cout << "log10(0) isn't available\n";
- exit(-1);
- }
-}
-
-inline NodeID DataTree::AddCos(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::COS);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::COS);
- }
- else
- {
- return One;
- }
-}
-
-inline NodeID DataTree::AddSin(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::SIN);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::SIN);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddTan(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::TAN);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::TAN);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddACos(NodeID iArg1)
-{
- MToken lToken(iArg1, eTempResult, NullID, token::ACOS);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::ACOS);
-}
-
-inline NodeID DataTree::AddASin(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::ASIN);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::ASIN);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddATan(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::ATAN);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::ATAN);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddCosH(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::COSH);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::COSH);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddSinH(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::SINH);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::SINH);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddTanH(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::TANH);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::TANH);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddACosH(NodeID iArg1)
-{
- MToken lToken(iArg1, eTempResult, NullID, token::ACOSH);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::ACOSH);
-}
-
-inline NodeID DataTree::AddASinH(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::ASINH);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::ASINH);
- }
- else
- {
- return Zero;
- }
-}
-
-inline NodeID DataTree::AddATanH(NodeID iArg1)
-{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::ATANH);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::ATANH);
- }
- else
- {
- return Zero;
- }
-}
+ //! Raised when a local parameter is declared twice
+ class LocalParameterException
+ {
+ public:
+ string name;
+ LocalParameterException(const string &name_arg) : name(name_arg) {}
+ };
+
+ NodeID AddNumConstant(const string &value);
+ NodeID AddVariable(const string &name, int lag = 0);
+ //! Adds "arg1+arg2" to model tree
+ NodeID AddPlus(NodeID iArg1, NodeID iArg2);
+ //! Adds "arg1-arg2" to model tree
+ NodeID AddMinus(NodeID iArg1, NodeID iArg2);
+ //! Adds "-arg" to model tree
+ NodeID AddUMinus(NodeID iArg1);
+ //! Adds "arg1*arg2" to model tree
+ NodeID AddTimes(NodeID iArg1, NodeID iArg2);
+ //! Adds "arg1/arg2" to model tree
+ NodeID AddDivide(NodeID iArg1, NodeID iArg2);
+ //! Adds "arg1^arg2" to model tree
+ NodeID AddPower(NodeID iArg1, NodeID iArg2);
+ //! Adds "exp(arg)" to model tree
+ NodeID AddExp(NodeID iArg1);
+ //! Adds "log(arg)" to model tree
+ NodeID AddLog(NodeID iArg1);
+ //! Adds "log10(arg)" to model tree
+ NodeID AddLog10(NodeID iArg1);
+ //! Adds "cos(arg)" to model tree
+ NodeID AddCos(NodeID iArg1);
+ //! Adds "sin(arg)" to model tree
+ NodeID AddSin(NodeID iArg1);
+ //! Adds "tan(arg)" to model tree
+ NodeID AddTan(NodeID iArg1);
+ //! Adds "acos(arg)" to model tree
+ NodeID AddACos(NodeID iArg1);
+ //! Adds "asin(arg)" to model tree
+ NodeID AddASin(NodeID iArg1);
+ //! Adds "atan(arg)" to model tree
+ NodeID AddATan(NodeID iArg1);
+ //! Adds "cosh(arg)" to model tree
+ NodeID AddCosH(NodeID iArg1);
+ //! Adds "sinh(arg)" to model tree
+ NodeID AddSinH(NodeID iArg1);
+ //! Adds "tanh(arg)" to model tree
+ NodeID AddTanH(NodeID iArg1);
+ //! Adds "acosh(arg)" to model tree
+ NodeID AddACosH(NodeID iArg1);
+ //! Adds "asinh(arg)" to model tree
+ NodeID AddASinH(NodeID iArg1);
+ //! Adds "atanh(args)" to model tree
+ NodeID AddATanH(NodeID iArg1);
+ //! Adds "sqrt(arg)" to model tree
+ NodeID AddSqRt(NodeID iArg1);
+ //! Adds "arg1=arg2" to model tree
+ NodeID AddEqual(NodeID iArg1, NodeID iArg2);
+ void AddLocalParameter(const string &name, NodeID value) throw (LocalParameterException);
+};
-inline NodeID DataTree::AddSqRt(NodeID iArg1)
+inline NodeID
+DataTree::AddUnaryOp(UnaryOpcode op_code, NodeID arg)
{
- if (iArg1 != Zero)
- {
- MToken lToken(iArg1, eTempResult, NullID, token::SQRT);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::SQRT);
- }
+ unary_op_node_map_type::iterator it = unary_op_node_map.find(make_pair(arg, op_code));
+ if (it != unary_op_node_map.end())
+ return it->second;
else
- {
- return Zero;
- }
+ return new UnaryOpNode(*this, op_code, arg);
}
-inline NodeID DataTree::AddEqual(NodeID iArg1, NodeID iArg2)
+inline NodeID
+DataTree::AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2)
{
- if (iArg1 == Zero && iArg2 == Zero)
- return ZeroEqZero;
+ binary_op_node_map_type::iterator it = binary_op_node_map.find(make_pair(make_pair(arg1, arg2), op_code));
+ if (it != binary_op_node_map.end())
+ return it->second;
else
- {
- MToken lToken(iArg1, eTempResult, iArg2, token::EQUAL);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- return ID;
-
- return PushToken(iArg1, token::EQUAL, iArg2);
- }
-}
-
-inline NodeID DataTree::AddAssign(NodeID iArg1, NodeID iArg2)
-{
- MToken lToken(iArg1, eTempResult, iArg2, token::ASSIGN);
-
- NodeID ID = getIDOfToken(lToken);
- if (ID != NullID)
- {
- return ID;
- }
- return PushToken(iArg1, token::ASSIGN, iArg2);
+ return new BinaryOpNode(*this, arg1, op_code, arg2);
}
-//------------------------------------------------------------------------------
#endif
diff --git a/parser.src/include/ExprNode.hh b/parser.src/include/ExprNode.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dfd1e16eebf6ded514d9d9774374b4fb94e88eb4
--- /dev/null
+++ b/parser.src/include/ExprNode.hh
@@ -0,0 +1,174 @@
+#ifndef _EXPR_NODE_HH
+#define _EXPR_NODE_HH
+
+using namespace std;
+
+#include <set>
+#include <map>
+
+#include "SymbolTableTypes.hh"
+
+class DataTree;
+
+typedef class ExprNode *NodeID;
+
+struct ExprNodeLess;
+
+//! Type for set of temporary terms
+/*! They are ordered by index number thanks to ExprNodeLess */
+typedef set<NodeID, ExprNodeLess> temporary_terms_type;
+
+//! Base class for expression nodes
+class ExprNode
+{
+ friend class DataTree;
+ friend class ModelTree;
+ friend class ExprNodeLess;
+ friend class NumConstNode;
+ friend class VariableNode;
+ friend class UnaryOpNode;
+ friend class BinaryOpNode;
+
+private:
+ //! Computes derivative w.r. to variable varID (but doesn't store it in derivatives map)
+ /*! You shoud use getDerivative() to get the benefit of symbolic a priori and of caching */
+ virtual NodeID computeDerivative(int varID) = 0;
+
+protected:
+ //! Reference to the enclosing DataTree
+ DataTree &datatree;
+
+ //! Index number
+ int idx;
+
+ //! Set of variable IDs with respect to which the derivative is potentially non-null
+ set<int> non_null_derivatives;
+
+ //! Used for caching of first order derivatives (when non-null)
+ map<int, NodeID> derivatives;
+
+ //! Cost of computing current node
+ /*! Nodes included in temporary_terms are considered having a null cost */
+ virtual int cost(temporary_terms_type &temporary_terms) const;
+
+public:
+ ExprNode(DataTree &datatree_arg);
+ virtual ~ExprNode();
+
+ //! Returns derivative w.r. to variable varID
+ /*! Uses a symbolic a priori to pre-detect null derivatives, and caches the result for other derivatives (to avoid computing it several times)
+ For an equal node, returns the derivative of lhs minus rhs */
+ NodeID getDerivative(int varID);
+
+ //! Returns precedence of node
+ /*! Equals 100 for constants, variables, unary ops, and temporary terms */
+ virtual int precedence(const temporary_terms_type &temporary_terms) const;
+
+ //! Fills temporary_terms set, using reference counts
+ /*! A node will be marked as a temporary term if it is referenced at least two times (i.e. has at least two parents), and has a computing cost (multiplied by reference count) greater to datatree.min_cost */
+ virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms) const;
+
+ //! Writes output of node, using a Txxx notation for nodes in temporary_terms
+ virtual void writeOutput(ostream &output, bool is_dynamic, const temporary_terms_type &temporary_terms) const = 0;
+};
+
+//! Object used to compare two nodes (using their indexes)
+struct ExprNodeLess
+{
+ bool operator()(NodeID arg1, NodeID arg2) const
+ {
+ return arg1->idx < arg2->idx;
+ }
+};
+
+//! Numerical constant node
+class NumConstNode : public ExprNode
+{
+private:
+ //! Id from numerical constants table
+ const int id;
+ virtual NodeID computeDerivative(int varID);
+public:
+ NumConstNode(DataTree &datatree_arg, int id_arg);
+ virtual void writeOutput(ostream &output, bool is_dynamic, const temporary_terms_type &temporary_terms) const;
+};
+
+//! Symbol or variable node
+class VariableNode : public ExprNode
+{
+private:
+ //! Id of the symbol/variable
+ /*! For an endogenous, exogenous or recursive variable, the id is taken from the variable table (before sorting). For other types of symbols, it's the id from the symbol table. */
+ const int id;
+ const Type type;
+ virtual NodeID computeDerivative(int varID);
+public:
+ VariableNode(DataTree &datatree_arg, int id_arg, Type type_arg);
+ virtual void writeOutput(ostream &output, bool is_dynamic, const temporary_terms_type &temporary_terms) const;
+};
+
+enum UnaryOpcode
+ {
+ oUminus,
+ oExp,
+ oLog,
+ oLog10,
+ oCos,
+ oSin,
+ oTan,
+ oAcos,
+ oAsin,
+ oAtan,
+ oCosh,
+ oSinh,
+ oTanh,
+ oAcosh,
+ oAsinh,
+ oAtanh,
+ oSqrt
+ };
+
+//! Unary operator node
+class UnaryOpNode : public ExprNode
+{
+ friend class DataTree;
+private:
+ const NodeID arg;
+ const UnaryOpcode op_code;
+ virtual NodeID computeDerivative(int varID);
+
+ int cost(temporary_terms_type &temporary_terms) const;
+public:
+ UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg);
+ virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms) const;
+ virtual void writeOutput(ostream &output, bool is_dynamic, const temporary_terms_type &temporary_terms) const;
+};
+
+enum BinaryOpcode
+ {
+ oPlus,
+ oMinus,
+ oTimes,
+ oDivide,
+ oPower,
+ oEqual
+ };
+
+//! Binary operator node
+class BinaryOpNode : public ExprNode
+{
+ friend class ModelTree;
+private:
+ const NodeID arg1, arg2;
+ const BinaryOpcode op_code;
+ virtual NodeID computeDerivative(int varID);
+ int cost(temporary_terms_type &temporary_terms) const;
+public:
+ BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
+ BinaryOpcode op_code_arg, const NodeID arg2_arg);
+ virtual int precedence(const temporary_terms_type &temporary_terms) const;
+ virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms) const;
+ virtual void writeOutput(ostream &output, bool is_dynamic, const temporary_terms_type &temporary_terms) const;
+};
+
+#endif
diff --git a/parser.src/include/ModFile.hh b/parser.src/include/ModFile.hh
index 96060c86dccb1a7bd207ed258d4556d28d7babe1..edfdd8181cefc5b8637ebfc7fa4af3735c0aaeca 100644
--- a/parser.src/include/ModFile.hh
+++ b/parser.src/include/ModFile.hh
@@ -45,7 +45,7 @@ public:
\param clear_all Should a "clear all" instruction be written to output ?
\todo make this method const
*/
- void writeOutputFiles(const string &basename, bool clear_all);
+ void writeOutputFiles(const string &basename, bool clear_all) const;
};
#endif // ! MOD_FILE_HH
diff --git a/parser.src/include/ModelTree.hh b/parser.src/include/ModelTree.hh
index 378acc4b20745f29e32ae1f872ab0690f56b4c5e..c95157fef0c31a0f5244386d271128e93503b021 100644
--- a/parser.src/include/ModelTree.hh
+++ b/parser.src/include/ModelTree.hh
@@ -1,72 +1,73 @@
#ifndef _MODELTREE_HH
#define _MODELTREE_HH
-//------------------------------------------------------------------------------
-/*! \file
- \version 1.0
- \date 04/13/2003
- \par This file defines the ModelTree class.
-*/
-//------------------------------------------------------------------------------
+
+using namespace std;
+
#include <string>
#include <vector>
-#include <list>
-#include <stack>
-#include <sstream>
-#include <fstream>
+#include <map>
#include <ostream>
-//------------------------------------------------------------------------------
+
#include "SymbolTable.hh"
#include "NumericalConstants.hh"
-#include "ModelTypes.hh"
#include "DataTree.hh"
-//------------------------------------------------------------------------------
-/*!
- \class ModelTree
- \brief Stores a model equations and derivatives.
-*/
+
+//! Stores a model's equations and derivatives
class ModelTree : public DataTree
{
-private :
- /*! Stores ID of equations and their derivatives */
- std::vector<std::vector<DerivativeIndex> > mDerivativeIndex;
- /*!
- A token is writen as a temporary expression
- if its cost is greater than min_cost
+private:
+ //! Stores declared equations
+ vector<BinaryOpNode *> equations;
+
+ typedef map<pair<int, int>, NodeID> first_derivatives_type;
+ //! First order derivatives
+ /*! First index is equation number, second is variable w.r. to which is computed the derivative.
+ Only non-null derivatives are stored in the map.
+ Variable indexes used are those of the variable_table, before sorting.
*/
- int min_cost;
- /*! left and right parentheses can be (,[ or ),] */
- char lpar, rpar;
- //! Reference to numerical constants table
- const NumericalConstants &num_constants;
+ first_derivatives_type first_derivatives;
- /*! Computes argument derivative */
- inline NodeID DeriveArgument(NodeID iArg, Type iType, int iVarID);
- /*! Gets output argument of terminal token */
- inline std::string getArgument(NodeID id, Type type, EquationType iEquationType);
- /*! Gets expression of part of model tree */
- std::string getExpression(NodeID StartID, EquationType iEquationType);
- inline int optimize(NodeID id);
- /*! Computes derivatives of ModelTree */
- void derive(int iOrder);
+ typedef map<pair<int, pair<int, int> >, NodeID> second_derivatives_type;
+ //! Second order derivatives
+ /*! First index is equation number, second and third are variables w.r. to which is computed the derivative.
+ Only non-null derivatives are stored in the map.
+ Contains only second order derivatives where var1 >= var2 (for obvious symmetry reasons).
+ Variable indexes used are those of the variable_table, before sorting.
+ */
+ second_derivatives_type second_derivatives;
+
+ //! Temporary terms (those which will be noted Txxxx)
+ temporary_terms_type temporary_terms;
+
+ //! Computes derivatives of ModelTree
+ void derive(int order);
+ //! Computes temporary terms
+ void computeTemporaryTerms(int order);
+
+ //! Writes temporary terms
+ void writeTemporaryTerms(ostream &output, bool is_dynamic) const;
+ //! Writes local parameters
+ void writeLocalParameters(ostream &output, bool is_dynamic) const;
+ //! Writes model equations
+ void writeModelEquations(ostream &output, bool is_dynamic) const;
//! Writes the static model equations and its derivatives
/*! \todo handle hessian in C output */
- void writeStaticModel(std::ostream &StaticOutput);
+ void writeStaticModel(ostream &StaticOutput) const;
//! Writes the dynamic model equations and its derivatives
- void writeDynamicModel(std::ostream &DynamicOutput);
+ void writeDynamicModel(ostream &DynamicOutput) const;
//! Writes static model file (Matlab version)
- void writeStaticMFile(const std::string &static_basename);
+ void writeStaticMFile(const string &static_basename) const;
//! Writes static model file (C version)
- void writeStaticCFile(const std::string &static_basename);
+ void writeStaticCFile(const string &static_basename) const;
//! Writes dynamic model file (Matlab version)
- void writeDynamicMFile(const std::string &dynamic_basename);
+ void writeDynamicMFile(const string &dynamic_basename) const;
//! Writes dynamic model file (C version)
- void writeDynamicCFile(const std::string &dynamic_basename);
+ void writeDynamicCFile(const string &dynamic_basename) const;
public:
- //! Constructor
- ModelTree(SymbolTable &symbol_table_arg, const NumericalConstants &num_constants);
- //! Number of equations contained in this model tree
- int eq_nbr;
+ ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
+ //! Declare a node as an equation of the model
+ void addEquation(NodeID eq);
//! Do some checking
void checkPass() const;
//! Whether dynamic Jacobian (w.r. to endogenous) should be written
@@ -82,13 +83,11 @@ public:
calling this function */
void computingPass();
//! Writes model initialization and lead/lag incidence matrix to output
- void writeOutput(std::ostream &output) const;
+ void writeOutput(ostream &output) const;
//! Writes static model file
- /*! \todo make this method const */
- void writeStaticFile(const std::string &basename);
+ void writeStaticFile(const string &basename) const;
//! Writes dynamic model file
- /*! \todo make this method const */
- void writeDynamicFile(const std::string &basename);
+ void writeDynamicFile(const string &basename) const;
};
-//------------------------------------------------------------------------------
+
#endif
diff --git a/parser.src/include/ModelTypes.hh b/parser.src/include/ModelTypes.hh
deleted file mode 100644
index d7efbcd5ddded4babcbc365e26feb17c86f37489..0000000000000000000000000000000000000000
--- a/parser.src/include/ModelTypes.hh
+++ /dev/null
@@ -1,230 +0,0 @@
-#ifndef _MODELTYPES_HH
-#define _MODELTYPES_HH
-//------------------------------------------------------------------------------
-/*! \file
- \version 1.0
- \date 04/13/2004
- \par This file defines types for ModelTree class.
-*/
-//------------------------------------------------------------------------------
-#include <string>
-#include <map>
-
-typedef class MetaToken* NodeID;
-
-/*!
- \struct DerivativeIndex
- \brief Derivative index structure
-*/
-//------------------------------------------------------------------------------
-struct DerivativeIndex
-{
- /*!/ id of "equal token" in mModelTree */
- NodeID token_id;
- /*! id of equation being derived */
- int equation_id;
- /*! column number of Jacobian or Hessian matrix */
- int derivators;
- /*! Constructor */
- inline DerivativeIndex(NodeID iToken_id, int iEquationID, int iDerivator)
- {
- token_id = iToken_id;
- equation_id = iEquationID;
- derivators = iDerivator;
- }
-};
-/*!
- \struct MToken
- \brief Basic token structure, essentially for computing hash keys
-*/
-//------------------------------------------------------------------------------
-struct MToken
-{
-public :
- /*! ID of first operand */
- NodeID id1;
- /*! Type of first operand */
- Type type1;
- /*! ID of second operand */
- NodeID id2;
- /*! Operator code */
- int op_code;
- /*! Constructor */
- inline MToken(NodeID iId1 = NULL, Type iType1 = eUNDEF, NodeID iId2 = NULL, int iOpCode = -1)
- {
- id1 = iId1;
- type1 = iType1;
- id2 = iId2;
- op_code = iOpCode;
-
- }
- /*! Copy constructor */
- inline MToken(const MToken& t)
- {
- id1 = t.id1;
- id2 = t.id2;
- type1 = t.type1;
- op_code = t.op_code;
- }
- /*! Destructor */
- ~MToken() { }
-};
-
-/*!
- \struct MTokenLess
- \brief Class which defines a lexicographic order over MToken, used in std::map
-*/
-struct MTokenLess
-{
- bool operator()(const MToken &n1, const MToken &n2) const
- {
- if (n1.id1 != n2.id1)
- return(n1.id1 < n2.id1);
- else if (n1.id2 != n2.id2)
- return(n1.id2 < n2.id2);
- else if (n1.type1 != n2.type1)
- return(n1.type1 < n2.type1);
- else
- return(n1.op_code < n2.op_code);
- }
-};
-
-// For MetaToken::getDerivativeAddress()
-class ModelTree;
-
-/*!
- \struct MetaToken
- \brief Meta token structure
-*/
-struct MetaToken : public MToken
-{
-private:
- /*! Address of first order partial derivative */
- std::map<int, NodeID, std::less<int> > d1;
-public :
- /*! Index of token starting from zero */
- int idx;
- /*!
- Number of times a given token is referenced?
- Each element of vector refer to derivative 0, 1, 2, etc
- */
- std::vector<int> reference_count;
- /*! Commulative computation time cost */
- int cost;
- /*!
- If the token can be a temporary result, this flag is set to 1.
- If it can not be a temporary result (operation with constants etc)
- it is set to 0.
- After writing corresponding expression for the first time it is set to -1.
- */
- int tmp_status;
- /*! Name of unkown functions */
- std::string func_name;
- /*! Expression of token, exp[0] for static model
- exp[1] for dynamic model */
- std::string exp[2];
- /*! Flag : operand 1 is followed or not */
- bool followed1;
- /*! Flag : operand 2 is followed or not */
- bool followed2;
- /*! Operator name */
- std::string op_name;
- /*! set to one when left node has been treated */
- int left_done;
- /*! set to one when right node has been treated */
- int right_done;
- /*! set to one if closing parenthesis after token */
- int close_parenthesis;
-
- /*! Ordered list of variable IDs with respect to which the derivative is potentially non-null */
- std::vector<int> non_null_derivatives;
-
- /*! Constructor */
- inline MetaToken(NodeID iId1 = NULL, Type iType1 = eUNDEF,NodeID iId2= NULL, int iOpCode = -1)
- : MToken(iId1, iType1,iId2,iOpCode)
- {
- cost = 0;
- followed1 = followed2 = false;
- tmp_status = 0;
- left_done = 0;
- right_done = 0;
- close_parenthesis = 0;
-
- switch(type1)
- {
- case eEndogenous:
- case eExogenous:
- case eExogenousDet:
- case eRecursiveVariable:
- // For a variable, the only non-null derivative is with respect to itself
- non_null_derivatives.push_back((long int) id1);
- break;
- case eParameter:
- case eLocalParameter:
- case eNumericalConstant:
- case eUNDEF:
- // All derivatives are null
- break;
- case eTempResult:
- if (!id2)
- non_null_derivatives = id1->non_null_derivatives;
- else
- {
- // Compute set union of id1->non_null_derivates and id2->non_null_derivatives
- // The result is ordered by ascending order
- std::vector<int>::iterator it1, it2;
- it1 = id1->non_null_derivatives.begin();
- it2 = id2->non_null_derivatives.begin();
- while(it1 != id1->non_null_derivatives.end()
- && it2 != id2->non_null_derivatives.end())
- {
- if (*it1 == *it2)
- {
- non_null_derivatives.push_back(*it1++);
- it2++;
- }
- else if (*it1 < *it2)
- non_null_derivatives.push_back(*it1++);
- else
- non_null_derivatives.push_back(*it2++);
- }
- while(it1 != id1->non_null_derivatives.end())
- non_null_derivatives.push_back(*it1++);
- while(it2 != id2->non_null_derivatives.end())
- non_null_derivatives.push_back(*it2++);
- }
- break;
- }
- }
- /*! Destructor */
- ~MetaToken()
- {
- }
-
- /*! Sets derivative with respect to iVarID equal to iDerivative */
- inline void setDerivativeAddress(NodeID iDerivative, int iVarID)
- {
- d1[iVarID] = iDerivative;
- }
-
- //! Get derivative with respect to iVarID.
- /*
- Defined in ModelTree.cc because it needs a reference to the enclosing ModelTree,
- and is only used in that source file.
-
- \param iVarID variable with respect to get derivative
- \param model_tree enclosing ModelTree
- */
- inline NodeID getDerivativeAddress(int iVarID, const ModelTree &model_tree) const;
-};
-//------------------------------------------------------------------------------
-/*! Equation type enum */
-enum EquationType
- {
- eStaticEquations = 0, //!< Equation of static model
- eDynamicEquations = 1, //!< Equation of dynamic model
- eStaticDerivatives = 2, //!< Derivative of static model
- eDynamicDerivatives = 3 //!< Derivative of dynamic model
- };
-//------------------------------------------------------------------------------
-#endif
diff --git a/parser.src/include/SymbolTable.hh b/parser.src/include/SymbolTable.hh
index a6f0c576748e4c153a53564cc805aa7f7dbf04d4..bdd80b2a5581d372ee078f9ba4afec3e8b52bb41 100644
--- a/parser.src/include/SymbolTable.hh
+++ b/parser.src/include/SymbolTable.hh
@@ -76,15 +76,15 @@ public :
*/
inline bool Exist(const std::string &name) const;
/*! Gets name by type and ID */
- inline std::string getNameByID(Type type, int id);
+ inline std::string getNameByID(Type type, int id) const;
/*! Gets tex name by type and ID */
- inline std::string getTexNameByID(Type type, int id);
+ inline std::string getTexNameByID(Type type, int id) const;
/*! Gets type by name */
inline Type getType(const std::string &name) const;
/*! Gets ID by name */
inline int getID(const std::string &name) const;
//! Write output of this class
- void writeOutput(std::ostream &output);
+ void writeOutput(std::ostream &output) const;
};
inline bool SymbolTable::Exist(const std::string &name) const
@@ -94,7 +94,7 @@ inline bool SymbolTable::Exist(const std::string &name) const
}
//------------------------------------------------------------------------------
-inline std::string SymbolTable::getNameByID(Type type,int id)
+inline std::string SymbolTable::getNameByID(Type type,int id) const
{
if (id >= 0 && (int)name_table[type].size() > id)
return(name_table[type][id]);
@@ -102,7 +102,7 @@ inline std::string SymbolTable::getNameByID(Type type,int id)
}
//------------------------------------------------------------------------------
-inline std::string SymbolTable::getTexNameByID(Type type,int id)
+inline std::string SymbolTable::getTexNameByID(Type type,int id) const
{
if (id >= 0 && (int)tex_name_table[type].size() > id)
return(tex_name_table[type][id]);