diff --git a/ExprNode.cc b/ExprNode.cc index eb4b6364..29099eea 100644 --- a/ExprNode.cc +++ b/ExprNode.cc @@ -276,7 +276,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type, { case oMatlabDynamicModel: case oCDynamicModel: - i = datatree.variable_table.getSortID(var_id) + OFFSET(output_type); + i = datatree.variable_table.getDynJacobianCol(var_id) + OFFSET(output_type); output << "y" << LPAR(output_type) << i << RPAR(output_type); break; case oMatlabStaticModel: diff --git a/ModelTree.cc b/ModelTree.cc index d475a87a..a50e0932 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -1694,14 +1694,10 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const << " g1 = NULL;" << endl << " if (nlhs >= 2)" << endl << " {" << endl - << " /* Set the output pointer to the output matrix g1. */" << endl; + << " /* Set the output pointer to the output matrix g1. */" << endl - if (computeJacobianExo) - mDynamicModelFile << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr() << ", mxREAL);" << endl; - else if (computeJacobian) - mDynamicModelFile << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.var_endo_nbr << ", mxREAL);" << endl; - - mDynamicModelFile << " /* Create a C pointer to a copy of the output matrix g1. */" << endl + << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.getDynJacobianColsNbr(computeJacobianExo) << ", mxREAL);" << endl + << " /* Create a C pointer to a copy of the output matrix g1. */" << endl << " g1 = mxGetPr(plhs[1]);" << endl << " }" << endl << endl @@ -1709,7 +1705,7 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const << " if (nlhs >= 3)" << endl << " {" << endl << " /* Set the output pointer to the output matrix g2. */" << endl - << " plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr()*variable_table.get_dyn_var_nbr() << ", mxREAL);" << endl + << " plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.getDynJacobianColsNbr(computeJacobianExo)*variable_table.getDynJacobianColsNbr(computeJacobianExo) << ", mxREAL);" << endl << " /* Create a C pointer to a copy of the output matrix g1. */" << endl << " g2 = mxGetPr(plhs[2]);" << endl << " }" << endl @@ -3295,11 +3291,7 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const writeModelEquations(model_output, output_type); int nrows = equations.size(); - int nvars; - if (computeJacobianExo) - nvars = variable_table.get_dyn_var_nbr(); - else - nvars = variable_table.var_endo_nbr; + int nvars = variable_table.getDynJacobianColsNbr(computeJacobianExo); int nvars_sq = nvars * nvars; // Writing Jacobian @@ -3315,7 +3307,7 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const { ostringstream g1; g1 << " g1"; - matrixHelper(g1, eq, variable_table.getSortID(var), output_type); + matrixHelper(g1, eq, variable_table.getDynJacobianCol(var), output_type); jacobian_output << g1.str() << "=" << g1.str() << "+"; d1->writeOutput(jacobian_output, output_type, temporary_terms); @@ -3333,8 +3325,8 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const int var2 = it->first.second.second; NodeID d2 = it->second; - int id1 = variable_table.getSortID(var1); - int id2 = variable_table.getSortID(var2); + int id1 = variable_table.getDynJacobianCol(var1); + int id2 = variable_table.getDynJacobianCol(var2); int col_nb = id1*nvars+id2; int col_nb_sym = id2*nvars+id1; @@ -3367,9 +3359,9 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const int var3 = it->first.second.second.second; NodeID d3 = it->second; - int id1 = variable_table.getSortID(var1); - int id2 = variable_table.getSortID(var2); - int id3 = variable_table.getSortID(var3); + int id1 = variable_table.getDynJacobianCol(var1); + int id2 = variable_table.getDynJacobianCol(var2); + int id3 = variable_table.getDynJacobianCol(var3); // Reference column number for the g3 matrix int ref_col = id1 * nvars_sq + id2 * nvars + id3; @@ -3507,7 +3499,7 @@ ModelTree::writeOutput(ostream &output) const try { int varID = variable_table.getID(eEndogenous, endoID, lag); - output << " " << variable_table.getSortID(varID) + 1; + output << " " << variable_table.getDynJacobianCol(varID) + 1; } catch(VariableTable::UnknownVariableKeyException &e) { @@ -3689,8 +3681,8 @@ ModelTree::computingPass(const eval_context_type &eval_context) { cout << equations.size() << " equation(s) found" << endl; - // Sorting variable table - variable_table.sort(); + // Computes dynamic jacobian columns + variable_table.computeDynJacobianCols(); // Determine derivation order int order = 1; diff --git a/VariableTable.cc b/VariableTable.cc index f70bce08..f994e9f6 100644 --- a/VariableTable.cc +++ b/VariableTable.cc @@ -17,6 +17,8 @@ * along with Dynare. If not, see . */ +#include + #include "VariableTable.hh" VariableTable::VariableTable(const SymbolTable &symbol_table_arg) : @@ -31,10 +33,10 @@ VariableTable::VariableTable(const SymbolTable &symbol_table_arg) : } int -VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySortedException) +VariableTable::addVariable(Type type, int symb_id, int lag) throw (DynJacobianColsAlreadyComputedException) { - if (sorted_ids_table.size() != 0) - throw AlreadySortedException(); + if (dyn_jacobian_cols_table.size() != 0) + throw DynJacobianColsAlreadyComputedException(); var_key_type key = make_pair(make_pair(type, lag), symb_id); @@ -48,13 +50,6 @@ VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySorted variable_table[key] = var_id; inv_variable_table[var_id] = key; - if (type == eEndogenous) - var_endo_nbr++; - if (type == eExogenous) - var_exo_nbr++; - if (type == eExogenousDet) - var_exo_det_nbr++; - // Setting maximum and minimum lags if (max_lead < lag) max_lead = lag; @@ -64,18 +59,21 @@ VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySorted switch(type) { case eEndogenous: + var_endo_nbr++; if (max_endo_lead < lag) max_endo_lead = lag; else if (-max_endo_lag > lag) max_endo_lag = -lag; break; case eExogenous: + var_exo_nbr++; if (max_exo_lead < lag) max_exo_lead = lag; else if (-max_exo_lag > lag) max_exo_lag = -lag; break; case eExogenousDet: + var_exo_det_nbr++; if (max_exo_det_lead < lag) max_exo_det_lead = lag; else if (-max_exo_det_lag > lag) @@ -88,21 +86,39 @@ VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySorted max_recur_lag = -lag; break; default: - ; + cerr << "VariableTable::addVariable(): forbidden variable type" << endl; + exit(-1); } return var_id; } void -VariableTable::sort() throw (AlreadySortedException) +VariableTable::computeDynJacobianCols() throw (DynJacobianColsAlreadyComputedException) { - if (sorted_ids_table.size() != 0) - throw AlreadySortedException(); + if (dyn_jacobian_cols_table.size() != 0) + throw DynJacobianColsAlreadyComputedException(); + dyn_jacobian_cols_table.resize(size()); + + variable_table_type::const_iterator it = variable_table.begin(); + + // Assign the first columns to endogenous, using the lexicographic order over (lag, symbol_id) implemented in variable_table map int sorted_id = 0; - sorted_ids_table.resize(size()); + while(it->first.first.first == eEndogenous) + { + dyn_jacobian_cols_table[it->second] = sorted_id++; + it++; + } - for(variable_table_type::const_iterator it = variable_table.begin(); - it != variable_table.end(); it++) - sorted_ids_table[it->second] = sorted_id++; + // Assign subsequent columns to exogenous and then exogenous deterministic, using an offset + symbol_id + while(it->first.first.first == eExogenous) + { + dyn_jacobian_cols_table[it->second] = var_endo_nbr + it->first.second; + it++; + } + while(it->first.first.first == eExogenousDet) + { + dyn_jacobian_cols_table[it->second] = var_endo_nbr + symbol_table.exo_nbr + it->first.second; + it++; + } } diff --git a/include/CodeInterpreter.hh b/include/CodeInterpreter.hh index 4c3df347..cfc78b2d 100644 --- a/include/CodeInterpreter.hh +++ b/include/CodeInterpreter.hh @@ -58,7 +58,7 @@ const int EVALUATE_FOREWARD_R=8; const int EVALUATE_BACKWARD_R=9; //! Enumeration of possible symbol types -/*! Warning: do not to change existing values: the order matters for VariableTable (at least ensure that eEndogenous is the first one), and the values matter for homotopy_setup command */ +/*! Warning: do not to change existing values: the order matters for VariableTable (at least for endogenous and exogenous types), and the values matter for homotopy_setup command */ enum Type { eEndogenous = 0, //!< Endogenous diff --git a/include/VariableTable.hh b/include/VariableTable.hh index 62bb44de..9ffb3990 100644 --- a/include/VariableTable.hh +++ b/include/VariableTable.hh @@ -35,7 +35,7 @@ private: //! A reference to the symbol table const SymbolTable &symbol_table; //! A variable is a tuple (type, lag, symbol_id) - /*! Warning: don't change the order of elements in the tuple, since this determines the lexicographic ordering in sort() */ + /*! Warning: don't change the order of elements in the tuple, since this determines the lexicographic ordering in computeDynJacobianCols() */ typedef pair, int> var_key_type; typedef map variable_table_type; @@ -45,16 +45,18 @@ private: typedef map inv_variable_table_type; //! Maps a variable ID to a tuple (type, lag, symbol_id) inv_variable_table_type inv_variable_table; - //! Contains the sorted IDs (indexed by variable IDs) - vector sorted_ids_table; -public: - VariableTable(const SymbolTable &symbol_table_arg); + //! Number of dynamic endogenous variables inside the model block int var_endo_nbr; //! Number of dynamic exogenous variables inside the model block int var_exo_nbr; //! Number of dynamic deterministic exogenous variables inside the model block int var_exo_det_nbr; + + //! Contains the columns indices for the dynamic jacobian (indexed by variable IDs) + vector dyn_jacobian_cols_table; +public: + VariableTable(const SymbolTable &symbol_table_arg); //! Maximum lag over all types of variables (positive value) int max_lag; //! Maximum lead over all types of variables @@ -91,17 +93,17 @@ public: int id; UnknownVariableIDException(int id_arg) : id(id_arg) {} }; - //! Thrown when getSortID() called before sort() - class NotYetSortedException + //! Thrown when getDynJacobianCol() called before computeDynJacobianCols() + class DynJacobianColsNotYetComputedException { }; - //! Thrown when sort() or addVariable() called after sort() - class AlreadySortedException + //! Thrown when computeDynJacobianCols() or addVariable() called after computeDynJacobianCols() + class DynJacobianColsAlreadyComputedException { }; //! Adds a variable in the table, and returns its (newly allocated) variable ID /*! Also works if the variable already exists */ - int addVariable(Type type, int symb_id, int lag) throw (AlreadySortedException); + int addVariable(Type type, int symb_id, int lag) throw (DynJacobianColsAlreadyComputedException); //! Return variable ID inline int getID(Type type, int symb_id, int lag) const throw (UnknownVariableKeyException); //! Return lag of variable @@ -112,26 +114,24 @@ public: inline Type getType(int var_id) const throw (UnknownVariableIDException); //! Get number of variables inline int size() const; - //! Get variable ID of sorted variable table - /*! In practice, only used for endogenous variables */ - inline int getSortID(int var_id) const throw (NotYetSortedException, UnknownVariableIDException); - //! Sorts variable table - /*! The order used is a lexicographic order over the tuple (type, lag, symbol_id) */ - void sort() throw (AlreadySortedException); - //! Get the number of dynamic variables - inline int get_dyn_var_nbr() const; + //! Get column index in dynamic jacobian + inline int getDynJacobianCol(int var_id) const throw (DynJacobianColsNotYetComputedException, UnknownVariableIDException); + //! Computes column indices in dynamic jacobian + void computeDynJacobianCols() throw (DynJacobianColsAlreadyComputedException); + //! Get the number of columns of dynamic jacobian (depending on whether we compute derivatives w.r. to exogenous) + inline int getDynJacobianColsNbr(bool computeJacobianExo) const; }; inline int -VariableTable::getSortID(int var_id) const throw (NotYetSortedException, UnknownVariableIDException) +VariableTable::getDynJacobianCol(int var_id) const throw (DynJacobianColsNotYetComputedException, UnknownVariableIDException) { - if (sorted_ids_table.size() == 0) - throw NotYetSortedException(); + if (dyn_jacobian_cols_table.size() == 0) + throw DynJacobianColsNotYetComputedException(); if (var_id < 0 || var_id >= size()) throw UnknownVariableIDException(var_id); - return sorted_ids_table[var_id]; + return dyn_jacobian_cols_table[var_id]; } inline int @@ -181,9 +181,12 @@ VariableTable::size() const } inline int -VariableTable::get_dyn_var_nbr() const +VariableTable::getDynJacobianColsNbr(bool computeJacobianExo) const { - return var_endo_nbr + symbol_table.exo_nbr + symbol_table.exo_det_nbr; + if (computeJacobianExo) + return var_endo_nbr + symbol_table.exo_nbr + symbol_table.exo_det_nbr; + else + return var_endo_nbr; } #endif