From 81b6823bb31cd6903a364ab792a036ad5a912a73 Mon Sep 17 00:00:00 2001 From: sebastien Date: Fri, 10 Jul 2009 14:22:40 +0000 Subject: [PATCH] Preprocessor: bugfixes to block_mfs option of steady * sort recursive variables in topological order in the _static.m file * collectEndogenous() returs type specific IDs, not symb_id! git-svn-id: https://www.dynare.org/svn/dynare/trunk@2833 ac1d8469-bf42-47a9-8791-bf33cf982152 --- ExprNode.hh | 4 +-- StaticModel.cc | 88 +++++++++++++++++++++++++++++++++++--------------- StaticModel.hh | 10 +++++- 3 files changed, 73 insertions(+), 29 deletions(-) diff --git a/ExprNode.hh b/ExprNode.hh index d0c43986..d9d48854 100644 --- a/ExprNode.hh +++ b/ExprNode.hh @@ -155,14 +155,14 @@ public: //! Computes the set of endogenous variables in the expression /*! - Endogenous are stored as integer pairs of the form (symb_id, lag). + Endogenous are stored as integer pairs of the form (type_specific_id, lag). They are added to the set given in argument. */ virtual void collectEndogenous(set > &result) const = 0; //! Computes the set of exogenous variables in the expression /*! - Exogenous are stored as integer pairs of the form (symb_id, lag). + Exogenous are stored as integer pairs of the form (type_specific_id, lag). They are added to the set given in argument. */ virtual void collectExogenous(set > &result) const = 0; diff --git a/StaticModel.cc b/StaticModel.cc index fd6f4d33..cd182a86 100644 --- a/StaticModel.cc +++ b/StaticModel.cc @@ -201,6 +201,7 @@ StaticModel::computingPass(bool block_mfs_arg, bool hessian, bool no_tmp_terms) computeNormalization(); computeSortedBlockDecomposition(); computeMFS(); + computeSortedRecursive(); computeBlockMFSJacobian(); } else if (!no_tmp_terms) @@ -248,7 +249,7 @@ StaticModel::computeNormalization() equations[i]->collectEndogenous(endo); for(set >::const_iterator it = endo.begin(); it != endo.end(); it++) - add_edge(i + n, symbol_table.getTypeSpecificID(it->first), g); + add_edge(i + n, it->first, g); } // Compute maximum cardinality matching @@ -346,7 +347,7 @@ StaticModel::computeNormalizedEquations(multimap &endo2eqs) const set > endo; equations[i]->get_arg2()->collectEndogenous(endo); - if (endo.find(make_pair(symb_id, 0)) != endo.end()) + if (endo.find(make_pair(symbol_table.getTypeSpecificID(symb_id), 0)) != endo.end()) continue; endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i)); @@ -378,7 +379,7 @@ StaticModel::computeSortedBlockDecomposition() equations[endo2eq[i]]->collectEndogenous(endo); for(set >::const_iterator it = endo.begin(); it != endo.end(); it++) - add_edge(symbol_table.getTypeSpecificID(it->first), i, g); + add_edge(it->first, i, g); } // Compute strongly connected components @@ -463,17 +464,62 @@ StaticModel::computeMFS() for(set >::const_iterator it2 = endo.begin(); it2 != endo.end(); ++it2) - { - const int tsid = symbol_table.getTypeSpecificID(it2->first); - if (blocks[b].find(tsid) != blocks[b].end()) // Add edge only if vertex member of this block - add_edge(tsid2vertex[tsid], tsid2vertex[*it], g); - } + if (blocks[b].find(it2->first) != blocks[b].end()) // Add edge only if vertex member of this block + add_edge(tsid2vertex[it2->first], tsid2vertex[*it], g); } // Compute minimum feedback set MFS::Minimal_set_of_feedback_vertex(blocksMFS[b], g); +#ifdef DEBUG cout << "Block " << b << ": " << blocksMFS[b].size() << "/" << blocks[b].size() << " in MFS" << endl; +#endif + } +} + +void +StaticModel::computeSortedRecursive() +{ + const int nblocks = blocks.size(); + blocksRecursive.clear(); + blocksRecursive.resize(nblocks); + + for(int b = 0; b < nblocks; b++) + { + // Construct the set of recursive vars + // The index in this vector will be the vertex descriptor in the graph + vector recurs_vars; + set_difference(blocks[b].begin(), blocks[b].end(), + blocksMFS[b].begin(), blocksMFS[b].end(), + back_inserter(recurs_vars)); + + // Construct graph representation of recursive vars + typedef adjacency_list DirectedGraph; + DirectedGraph dag(recurs_vars.size()); + set > endo; + for(int i = 0; i < (int) recurs_vars.size(); i++) + { + endo.clear(); + equations[endo2eq[recurs_vars[i]]]->get_arg2()->collectEndogenous(endo); + for(set >::const_iterator it = endo.begin(); + it != endo.end(); it++) + { + vector::const_iterator it2 = find(recurs_vars.begin(), recurs_vars.end(), it->first); + if (it2 != recurs_vars.end()) + { + int source_vertex = it2 - recurs_vars.begin(); + add_edge(source_vertex, i, dag); + } + } + } + // Compute topological sort + deque ordered_recurs_vertices; + topological_sort(dag, front_inserter(ordered_recurs_vertices)); // We use a front inserter because topological_sort returns the inverse order + + // Construct the final order + for(deque::const_iterator it = ordered_recurs_vertices.begin(); + it != ordered_recurs_vertices.end(); it++) + blocksRecursive[b].push_back(recurs_vars[*it]); } } @@ -483,17 +529,11 @@ StaticModel::computeBlockMFSJacobian() blocksMFSJacobian.clear(); for(int b = 0; b < (int) blocks.size(); b++) { - // Create the set of recursive variables (i.e. those not in the MFS) - set recurs_vars; - set_difference(blocks[b].begin(), blocks[b].end(), - blocksMFS[b].begin(), blocksMFS[b].end(), - inserter(recurs_vars, recurs_vars.begin())); - - // Create the map of recursive variables to their normalized equation - map recurs_vars_eqs; - for(set::const_iterator it = recurs_vars.begin(); - it != recurs_vars.end(); it++) - recurs_vars_eqs[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]]; + // Create the map of recursive vars to their normalized equation + map recursive2eq; + for(vector::const_iterator it = blocksRecursive[b].begin(); + it != blocksRecursive[b].end(); it++) + recursive2eq[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]]; for(set::const_iterator it = blocksMFS[b].begin(); it != blocksMFS[b].end(); it++) @@ -503,7 +543,7 @@ StaticModel::computeBlockMFSJacobian() it2 != blocksMFS[b].end(); it2++) { int deriv_id = symbol_table.getID(eEndogenous, *it2); - NodeID d = equations[eq_no]->getChainRuleDerivative(deriv_id, recurs_vars_eqs); + NodeID d = equations[eq_no]->getChainRuleDerivative(deriv_id, recursive2eq); if (d != Zero) blocksMFSJacobian[make_pair(eq_no, *it2)] = d; } @@ -536,12 +576,8 @@ StaticModel::writeStaticBlockMFSFile(ostream &output, const string &func_name) c { output << " case " << b+1 << endl << " % Variables not in minimum feedback set" << endl; - set recurs_vars; - set_difference(blocks[b].begin(), blocks[b].end(), - blocksMFS[b].begin(), blocksMFS[b].end(), - inserter(recurs_vars, recurs_vars.begin())); - for(set::const_iterator it = recurs_vars.begin(); - it != recurs_vars.end(); it++) + for(vector::const_iterator it = blocksRecursive[b].begin(); + it != blocksRecursive[b].end(); it++) { equations[endo2eq[*it]]->writeOutput(output, oMatlabStaticModel, temporary_terms_type()); output << ";" << endl; diff --git a/StaticModel.hh b/StaticModel.hh index 9f2934c6..104eb12d 100644 --- a/StaticModel.hh +++ b/StaticModel.hh @@ -42,6 +42,10 @@ private: /*! Elements of blocksMFS are subset of elements of blocks */ vector > blocksMFS; + //! Variables not in minimum feedback set for each block, sorted in topological order + /*! This is the set difference blocks - blocksMFS. The variables are sorted in topological order. */ + vector > blocksRecursive; + //! Jacobian for matrix restricted to MFS /*! Maps a pair (equation ID, endogenous type specific ID) to the derivative expression. Stores only non-null derivatives. */ map, NodeID> blocksMFSJacobian; @@ -65,8 +69,12 @@ private: /*! Must be called after computeSortedBlockDecomposition() */ void computeMFS(); - //! Computes derivatives of each MFS + //! For each block of the static model, computes resursive variables (those not in minimum feedback set), and sort them in topological order /*! Must be called after computeMFS() */ + void computeSortedRecursive(); + + //! Computes derivatives of each MFS + /*! Must be called after computeSortedRecursive() */ void computeBlockMFSJacobian(); //! Computes the list of equations which are already in normalized form