Commit 81b6823b authored by sebastien's avatar sebastien

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
parent 8a6d750c
...@@ -155,14 +155,14 @@ public: ...@@ -155,14 +155,14 @@ public:
//! Computes the set of endogenous variables in the expression //! 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. They are added to the set given in argument.
*/ */
virtual void collectEndogenous(set<pair<int, int> > &result) const = 0; virtual void collectEndogenous(set<pair<int, int> > &result) const = 0;
//! Computes the set of exogenous variables in the expression //! 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. They are added to the set given in argument.
*/ */
virtual void collectExogenous(set<pair<int, int> > &result) const = 0; virtual void collectExogenous(set<pair<int, int> > &result) const = 0;
......
...@@ -201,6 +201,7 @@ StaticModel::computingPass(bool block_mfs_arg, bool hessian, bool no_tmp_terms) ...@@ -201,6 +201,7 @@ StaticModel::computingPass(bool block_mfs_arg, bool hessian, bool no_tmp_terms)
computeNormalization(); computeNormalization();
computeSortedBlockDecomposition(); computeSortedBlockDecomposition();
computeMFS(); computeMFS();
computeSortedRecursive();
computeBlockMFSJacobian(); computeBlockMFSJacobian();
} }
else if (!no_tmp_terms) else if (!no_tmp_terms)
...@@ -248,7 +249,7 @@ StaticModel::computeNormalization() ...@@ -248,7 +249,7 @@ StaticModel::computeNormalization()
equations[i]->collectEndogenous(endo); equations[i]->collectEndogenous(endo);
for(set<pair<int, int> >::const_iterator it = endo.begin(); for(set<pair<int, int> >::const_iterator it = endo.begin();
it != endo.end(); it++) 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 // Compute maximum cardinality matching
...@@ -346,7 +347,7 @@ StaticModel::computeNormalizedEquations(multimap<int, int> &endo2eqs) const ...@@ -346,7 +347,7 @@ StaticModel::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
set<pair<int, int> > endo; set<pair<int, int> > endo;
equations[i]->get_arg2()->collectEndogenous(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; continue;
endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i)); endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
...@@ -378,7 +379,7 @@ StaticModel::computeSortedBlockDecomposition() ...@@ -378,7 +379,7 @@ StaticModel::computeSortedBlockDecomposition()
equations[endo2eq[i]]->collectEndogenous(endo); equations[endo2eq[i]]->collectEndogenous(endo);
for(set<pair<int, int> >::const_iterator it = endo.begin(); for(set<pair<int, int> >::const_iterator it = endo.begin();
it != endo.end(); it++) it != endo.end(); it++)
add_edge(symbol_table.getTypeSpecificID(it->first), i, g); add_edge(it->first, i, g);
} }
// Compute strongly connected components // Compute strongly connected components
...@@ -463,37 +464,76 @@ StaticModel::computeMFS() ...@@ -463,37 +464,76 @@ StaticModel::computeMFS()
for(set<pair<int, int> >::const_iterator it2 = endo.begin(); for(set<pair<int, int> >::const_iterator it2 = endo.begin();
it2 != endo.end(); ++it2) it2 != endo.end(); ++it2)
{ if (blocks[b].find(it2->first) != blocks[b].end()) // Add edge only if vertex member of this block
const int tsid = symbol_table.getTypeSpecificID(it2->first); add_edge(tsid2vertex[it2->first], tsid2vertex[*it], g);
if (blocks[b].find(tsid) != blocks[b].end()) // Add edge only if vertex member of this block
add_edge(tsid2vertex[tsid], tsid2vertex[*it], g);
}
} }
// Compute minimum feedback set // Compute minimum feedback set
MFS::Minimal_set_of_feedback_vertex(blocksMFS[b], g); MFS::Minimal_set_of_feedback_vertex(blocksMFS[b], g);
#ifdef DEBUG
cout << "Block " << b << ": " << blocksMFS[b].size() << "/" << blocks[b].size() << " in MFS" << endl; cout << "Block " << b << ": " << blocksMFS[b].size() << "/" << blocks[b].size() << " in MFS" << endl;
#endif
} }
} }
void void
StaticModel::computeBlockMFSJacobian() StaticModel::computeSortedRecursive()
{ {
blocksMFSJacobian.clear(); const int nblocks = blocks.size();
for(int b = 0; b < (int) blocks.size(); b++) blocksRecursive.clear();
blocksRecursive.resize(nblocks);
for(int b = 0; b < nblocks; b++)
{ {
// Create the set of recursive variables (i.e. those not in the MFS) // Construct the set of recursive vars
set<int> recurs_vars; // The index in this vector will be the vertex descriptor in the graph
vector<int> recurs_vars;
set_difference(blocks[b].begin(), blocks[b].end(), set_difference(blocks[b].begin(), blocks[b].end(),
blocksMFS[b].begin(), blocksMFS[b].end(), blocksMFS[b].begin(), blocksMFS[b].end(),
inserter(recurs_vars, recurs_vars.begin())); back_inserter(recurs_vars));
// Create the map of recursive variables to their normalized equation // Construct graph representation of recursive vars
map<int, NodeID> recurs_vars_eqs; typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
for(set<int>::const_iterator it = recurs_vars.begin(); DirectedGraph dag(recurs_vars.size());
it != recurs_vars.end(); it++) set<pair<int, int> > endo;
recurs_vars_eqs[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]]; for(int i = 0; i < (int) recurs_vars.size(); i++)
{
endo.clear();
equations[endo2eq[recurs_vars[i]]]->get_arg2()->collectEndogenous(endo);
for(set<pair<int, int> >::const_iterator it = endo.begin();
it != endo.end(); it++)
{
vector<int>::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<int> 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<int>::const_iterator it = ordered_recurs_vertices.begin();
it != ordered_recurs_vertices.end(); it++)
blocksRecursive[b].push_back(recurs_vars[*it]);
}
}
void
StaticModel::computeBlockMFSJacobian()
{
blocksMFSJacobian.clear();
for(int b = 0; b < (int) blocks.size(); b++)
{
// Create the map of recursive vars to their normalized equation
map<int, NodeID> recursive2eq;
for(vector<int>::const_iterator it = blocksRecursive[b].begin();
it != blocksRecursive[b].end(); it++)
recursive2eq[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]];
for(set<int>::const_iterator it = blocksMFS[b].begin(); for(set<int>::const_iterator it = blocksMFS[b].begin();
it != blocksMFS[b].end(); it++) it != blocksMFS[b].end(); it++)
...@@ -503,7 +543,7 @@ StaticModel::computeBlockMFSJacobian() ...@@ -503,7 +543,7 @@ StaticModel::computeBlockMFSJacobian()
it2 != blocksMFS[b].end(); it2++) it2 != blocksMFS[b].end(); it2++)
{ {
int deriv_id = symbol_table.getID(eEndogenous, *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) if (d != Zero)
blocksMFSJacobian[make_pair(eq_no, *it2)] = d; blocksMFSJacobian[make_pair(eq_no, *it2)] = d;
} }
...@@ -536,12 +576,8 @@ StaticModel::writeStaticBlockMFSFile(ostream &output, const string &func_name) c ...@@ -536,12 +576,8 @@ StaticModel::writeStaticBlockMFSFile(ostream &output, const string &func_name) c
{ {
output << " case " << b+1 << endl output << " case " << b+1 << endl
<< " % Variables not in minimum feedback set" << endl; << " % Variables not in minimum feedback set" << endl;
set<int> recurs_vars; for(vector<int>::const_iterator it = blocksRecursive[b].begin();
set_difference(blocks[b].begin(), blocks[b].end(), it != blocksRecursive[b].end(); it++)
blocksMFS[b].begin(), blocksMFS[b].end(),
inserter(recurs_vars, recurs_vars.begin()));
for(set<int>::const_iterator it = recurs_vars.begin();
it != recurs_vars.end(); it++)
{ {
equations[endo2eq[*it]]->writeOutput(output, oMatlabStaticModel, temporary_terms_type()); equations[endo2eq[*it]]->writeOutput(output, oMatlabStaticModel, temporary_terms_type());
output << ";" << endl; output << ";" << endl;
......
...@@ -42,6 +42,10 @@ private: ...@@ -42,6 +42,10 @@ private:
/*! Elements of blocksMFS are subset of elements of blocks */ /*! Elements of blocksMFS are subset of elements of blocks */
vector<set<int> > blocksMFS; vector<set<int> > 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<vector<int> > blocksRecursive;
//! Jacobian for matrix restricted to MFS //! Jacobian for matrix restricted to MFS
/*! Maps a pair (equation ID, endogenous type specific ID) to the derivative expression. Stores only non-null derivatives. */ /*! Maps a pair (equation ID, endogenous type specific ID) to the derivative expression. Stores only non-null derivatives. */
map<pair<int, int>, NodeID> blocksMFSJacobian; map<pair<int, int>, NodeID> blocksMFSJacobian;
...@@ -65,8 +69,12 @@ private: ...@@ -65,8 +69,12 @@ private:
/*! Must be called after computeSortedBlockDecomposition() */ /*! Must be called after computeSortedBlockDecomposition() */
void computeMFS(); 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() */ /*! Must be called after computeMFS() */
void computeSortedRecursive();
//! Computes derivatives of each MFS
/*! Must be called after computeSortedRecursive() */
void computeBlockMFSJacobian(); void computeBlockMFSJacobian();
//! Computes the list of equations which are already in normalized form //! Computes the list of equations which are already in normalized form
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment