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-bf33cf982152time-shift
parent
1db8b0d277
commit
4ba6a505db
|
@ -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)
|
||||||
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()
|
||||||
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
|
||||||
|
|
||||||
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()
|
||||||
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,17 +464,62 @@ 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
|
||||||
|
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<int> 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<vecS, vecS, directedS> DirectedGraph;
|
||||||
|
DirectedGraph dag(recurs_vars.size());
|
||||||
|
set<pair<int, int> > endo;
|
||||||
|
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]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -483,17 +529,11 @@ StaticModel::computeBlockMFSJacobian()
|
||||||
blocksMFSJacobian.clear();
|
blocksMFSJacobian.clear();
|
||||||
for(int b = 0; b < (int) blocks.size(); b++)
|
for(int b = 0; b < (int) blocks.size(); b++)
|
||||||
{
|
{
|
||||||
// Create the set of recursive variables (i.e. those not in the MFS)
|
// Create the map of recursive vars to their normalized equation
|
||||||
set<int> recurs_vars;
|
map<int, NodeID> recursive2eq;
|
||||||
set_difference(blocks[b].begin(), blocks[b].end(),
|
for(vector<int>::const_iterator it = blocksRecursive[b].begin();
|
||||||
blocksMFS[b].begin(), blocksMFS[b].end(),
|
it != blocksRecursive[b].end(); it++)
|
||||||
inserter(recurs_vars, recurs_vars.begin()));
|
recursive2eq[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]];
|
||||||
|
|
||||||
// Create the map of recursive variables to their normalized equation
|
|
||||||
map<int, NodeID> recurs_vars_eqs;
|
|
||||||
for(set<int>::const_iterator it = recurs_vars.begin();
|
|
||||||
it != recurs_vars.end(); it++)
|
|
||||||
recurs_vars_eqs[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()
|
||||||
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
|
||||||
{
|
{
|
||||||
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:
|
||||||
/*! 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:
|
||||||
/*! 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
|
||||||
|
|
Loading…
Reference in New Issue