Commit aedc60cb authored by sebastien's avatar sebastien

trunk preprocessor: added new option "block_mfs" to "steady"

* normalizes the static model
* computes its block decomposition, using topological order
* for each block, computes minimum feedback set of variables
* at this stage, only produces text output (no change in the computation of steady state)


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2798 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 762c5140
...@@ -32,6 +32,13 @@ SteadyStatement::SteadyStatement(const OptionsList &options_list_arg) : ...@@ -32,6 +32,13 @@ SteadyStatement::SteadyStatement(const OptionsList &options_list_arg) :
{ {
} }
void
SteadyStatement::checkPass(ModFileStructure &mod_file_struct)
{
if (options_list.num_options.find("block_mfs") != options_list.num_options.end())
mod_file_struct.steady_block_mfs_option = true;
}
void void
SteadyStatement::writeOutput(ostream &output, const string &basename) const SteadyStatement::writeOutput(ostream &output, const string &basename) const
{ {
...@@ -904,7 +911,7 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct) ...@@ -904,7 +911,7 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct)
void void
PlannerObjectiveStatement::computingPass() PlannerObjectiveStatement::computingPass()
{ {
model_tree->computingPass(true, false); model_tree->computingPass(false, true, false);
} }
void void
......
...@@ -34,6 +34,7 @@ private: ...@@ -34,6 +34,7 @@ private:
const OptionsList options_list; const OptionsList options_list;
public: public:
SteadyStatement(const OptionsList &options_list_arg); SteadyStatement(const OptionsList &options_list_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const; virtual void writeOutput(ostream &output, const string &basename) const;
}; };
......
...@@ -87,7 +87,7 @@ class ParsingDriver; ...@@ -87,7 +87,7 @@ class ParsingDriver;
%} %}
%token AR AUTOCORR %token AR AUTOCORR
%token BAYESIAN_IRF BETA_PDF BICGSTAB %token BAYESIAN_IRF BETA_PDF BICGSTAB BLOCK_MFS
%token BVAR_DENSITY BVAR_FORECAST %token BVAR_DENSITY BVAR_FORECAST
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA %token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN %token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
...@@ -631,6 +631,7 @@ steady_options_list : steady_options_list COMMA steady_options ...@@ -631,6 +631,7 @@ steady_options_list : steady_options_list COMMA steady_options
steady_options : o_solve_algo steady_options : o_solve_algo
| o_homotopy_mode | o_homotopy_mode
| o_homotopy_steps | o_homotopy_steps
| o_block_mfs
; ;
check : CHECK ';' check : CHECK ';'
...@@ -1434,6 +1435,7 @@ o_gsa_trans_ident : TRANS_IDENT EQUAL INT_NUMBER { driver.option_num("trans_iden ...@@ -1434,6 +1435,7 @@ o_gsa_trans_ident : TRANS_IDENT EQUAL INT_NUMBER { driver.option_num("trans_iden
o_homotopy_mode : HOMOTOPY_MODE EQUAL INT_NUMBER {driver.option_num("homotopy_mode",$3); }; o_homotopy_mode : HOMOTOPY_MODE EQUAL INT_NUMBER {driver.option_num("homotopy_mode",$3); };
o_homotopy_steps : HOMOTOPY_STEPS EQUAL INT_NUMBER {driver.option_num("homotopy_steps",$3); }; o_homotopy_steps : HOMOTOPY_STEPS EQUAL INT_NUMBER {driver.option_num("homotopy_steps",$3); };
o_block_mfs : BLOCK_MFS { driver.option_num("block_mfs", "1"); }
range : NAME ':' NAME range : NAME ':' NAME
{ {
......
...@@ -216,6 +216,7 @@ int sigma_e = 0; ...@@ -216,6 +216,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT>filename {return token::FILENAME;} <DYNARE_STATEMENT>filename {return token::FILENAME;}
<DYNARE_STATEMENT>diffuse_filter {return token::DIFFUSE_FILTER;} <DYNARE_STATEMENT>diffuse_filter {return token::DIFFUSE_FILTER;}
<DYNARE_STATEMENT>plot_priors {return token::PLOT_PRIORS;} <DYNARE_STATEMENT>plot_priors {return token::PLOT_PRIORS;}
<DYNARE_STATEMENT>block_mfs {return token::BLOCK_MFS;}
/* These four (var, varexo, varexo_det, parameters) are for change_type */ /* These four (var, varexo, varexo_det, parameters) are for change_type */
<DYNARE_STATEMENT>var { return token::VAR; } <DYNARE_STATEMENT>var { return token::VAR; }
......
...@@ -144,7 +144,7 @@ ModFile::computingPass(bool no_tmp_terms) ...@@ -144,7 +144,7 @@ ModFile::computingPass(bool no_tmp_terms)
{ {
// Compute static model and its derivatives // Compute static model and its derivatives
dynamic_model.toStatic(static_model); dynamic_model.toStatic(static_model);
static_model.computingPass(false, no_tmp_terms); static_model.computingPass(mod_file_struct.steady_block_mfs_option, false, no_tmp_terms);
// Set things to compute for dynamic model // Set things to compute for dynamic model
......
...@@ -30,7 +30,8 @@ ModFileStructure::ModFileStructure() : ...@@ -30,7 +30,8 @@ ModFileStructure::ModFileStructure() :
order_option(0), order_option(0),
bvar_density_present(false), bvar_density_present(false),
bvar_forecast_present(false), bvar_forecast_present(false),
identification_present(false) identification_present(false),
steady_block_mfs_option(false)
{ {
} }
......
...@@ -56,6 +56,8 @@ public: ...@@ -56,6 +56,8 @@ public:
bool bvar_forecast_present; bool bvar_forecast_present;
//! Whether an identification statement is present //! Whether an identification statement is present
bool identification_present; bool identification_present;
//! Whether the option "block_mfs" is used on steady statement
bool steady_block_mfs_option;
}; };
class Statement class Statement
......
...@@ -20,18 +20,22 @@ ...@@ -20,18 +20,22 @@
#include <cstdlib> #include <cstdlib>
#include <cassert> #include <cassert>
#include <deque>
#include <algorithm> #include <algorithm>
#include <functional> #include <functional>
/* #ifdef DEBUG
#include <ext/functional> # include <ext/functional>
using namespace __gnu_cxx; using namespace __gnu_cxx;
*/ #endif
#include <boost/graph/adjacency_list.hpp> #include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp> #include <boost/graph/max_cardinality_matching.hpp>
#include <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>
#include "StaticModel.hh" #include "StaticModel.hh"
#include "MinimumFeedbackSet.hh"
using namespace boost; using namespace boost;
...@@ -298,7 +302,7 @@ StaticModel::writeStaticFile(const string &basename) const ...@@ -298,7 +302,7 @@ StaticModel::writeStaticFile(const string &basename) const
} }
void void
StaticModel::computingPass(bool hessian, bool no_tmp_terms) StaticModel::computingPass(bool block_mfs, bool hessian, bool no_tmp_terms)
{ {
// Compute derivatives w.r. to all endogenous // Compute derivatives w.r. to all endogenous
set<int> vars; set<int> vars;
...@@ -319,24 +323,17 @@ StaticModel::computingPass(bool hessian, bool no_tmp_terms) ...@@ -319,24 +323,17 @@ StaticModel::computingPass(bool hessian, bool no_tmp_terms)
if (!no_tmp_terms) if (!no_tmp_terms)
computeTemporaryTerms(); computeTemporaryTerms();
/* if (block_mfs)
vector<int> endo2eq(equation_number());
computeNormalization(endo2eq);
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
for(int i = 0; i < symbol_table.endo_nbr(); i++)
{ {
if (natural_endo2eqs.count(i) == 0) vector<int> endo2eq(equation_number());
continue; computeNormalization(endo2eq);
pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i); vector<set<int> > blocks;
if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second) computeSortedBlockDecomposition(blocks, endo2eq);
cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
<< " not used." << endl; vector<set<int> > blocksMFS;
computeMFS(blocksMFS, blocks, endo2eq);
} }
*/
} }
int int
...@@ -358,9 +355,9 @@ StaticModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDExcepti ...@@ -358,9 +355,9 @@ StaticModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDExcepti
} }
void void
StaticModel::computeNormalization(vector<int> &endo_to_eq) const StaticModel::computeNormalization(vector<int> &endo2eq) const
{ {
int n = equation_number(); const int n = equation_number();
assert(n == symbol_table.endo_nbr()); assert(n == symbol_table.endo_nbr());
...@@ -384,15 +381,42 @@ StaticModel::computeNormalization(vector<int> &endo_to_eq) const ...@@ -384,15 +381,42 @@ StaticModel::computeNormalization(vector<int> &endo_to_eq) const
} }
// Compute maximum cardinality matching // Compute maximum cardinality matching
typedef vector<graph_traits<BipartiteGraph>::vertex_descriptor> mate_map_t; vector<int> mate_map(2*n);
mate_map_t mate_map(2*n);
#if 1
bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]); bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
#else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
fill(mate_map.begin(), mate_map.end(), graph_traits<BipartiteGraph>::null_vertex());
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
for(int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
int j = natural_endo2eqs.find(i)->second;
put(&mate_map[0], i, n+j);
put(&mate_map[0], n+j, i);
}
edmonds_augmenting_path_finder<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type> augmentor(g, &mate_map[0], get(vertex_index, g));
bool not_maximum_yet = true;
while(not_maximum_yet)
{
not_maximum_yet = augmentor.augment_matching();
}
augmentor.get_current_matching(&mate_map[0]);
bool check = maximum_cardinality_matching_verifier<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(vertex_index, g));
#endif
assert(check); assert(check);
// Check if all variables are normalized // Check if all variables are normalized
mate_map_t::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex()); vector<int>::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
if (it != mate_map.begin() + n) if (it != mate_map.begin() + n)
{ {
cerr << "ERROR: Could not normalize static model. Variable " cerr << "ERROR: Could not normalize static model. Variable "
...@@ -401,18 +425,44 @@ StaticModel::computeNormalization(vector<int> &endo_to_eq) const ...@@ -401,18 +425,44 @@ StaticModel::computeNormalization(vector<int> &endo_to_eq) const
exit(EXIT_FAILURE); exit(EXIT_FAILURE);
} }
#ifdef DEBUG
for(int i = 0; i < n; i++) for(int i = 0; i < n; i++)
cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i)) cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
<< " matched with equation " << (mate_map[i]-n+1) << endl; << " matched with equation " << (mate_map[i]-n+1) << endl;
#endif
assert((int) endo_to_eq.size() == n); assert((int) endo2eq.size() == n);
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to them // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
transform(mate_map.begin(), mate_map.begin() + n, endo_to_eq.begin(), bind2nd(minus<int>(), n)); transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));
#ifdef DEBUG
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
int n1 = 0, n2 = 0;
for(int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
n1++;
pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i);
if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second)
cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
<< " not used." << endl;
else
n2++;
}
cout << "Used " << n2 << " natural normalizations out of " << n1 << ", for a total of " << n << " equations." << endl;
#endif
} }
void void
StaticModel::computeNormalizedEquations(multimap<int, int> &endo_to_eqs) const StaticModel::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
{ {
for(int i = 0; i < equation_number(); i++) for(int i = 0; i < equation_number(); i++)
{ {
...@@ -429,7 +479,7 @@ StaticModel::computeNormalizedEquations(multimap<int, int> &endo_to_eqs) const ...@@ -429,7 +479,7 @@ StaticModel::computeNormalizedEquations(multimap<int, int> &endo_to_eqs) const
if (endo.find(make_pair(symb_id, 0)) != endo.end()) if (endo.find(make_pair(symb_id, 0)) != endo.end())
continue; continue;
endo_to_eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i)); endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl; cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
} }
} }
...@@ -439,3 +489,120 @@ StaticModel::writeLatexFile(const string &basename) const ...@@ -439,3 +489,120 @@ StaticModel::writeLatexFile(const string &basename) const
{ {
writeLatexModelFile(basename + "_static.tex", oLatexStaticModel); writeLatexModelFile(basename + "_static.tex", oLatexStaticModel);
} }
void
StaticModel::computeSortedBlockDecomposition(vector<set<int> > &blocks, const vector<int> &endo2eq) const
{
const int n = equation_number();
assert((int) endo2eq.size() == n);
// Compute graph representation of static model
typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
DirectedGraph g(n);
set<pair<int, int> > endo;
for(int i = 0; i < n; i++)
{
endo.clear();
equations[endo2eq[i]]->collectEndogenous(endo);
for(set<pair<int, int> >::const_iterator it = endo.begin();
it != endo.end(); it++)
add_edge(symbol_table.getTypeSpecificID(it->first), i, g);
}
// Compute strongly connected components
vector<int> endo2block(n);
int m = strong_components(g, &endo2block[0]);
// Create directed acyclic graph associated to the strongly connected components
DirectedGraph dag(m);
graph_traits<DirectedGraph>::edge_iterator ei, ei_end;
for(tie(ei, ei_end) = edges(g); ei != ei_end; ++ei)
{
int s = endo2block[source(*ei, g)];
int t = endo2block[target(*ei, g)];
if (s != t)
add_edge(s, t, dag);
}
// Compute topological sort of DAG (ordered list of unordered SCC)
deque<int> ordered2unordered;
topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order
// Construct mapping from unordered SCC to ordered SCC
vector<int> unordered2ordered(m);
for(int i = 0; i < m; i++)
unordered2ordered[ordered2unordered[i]] = i;
// Fill in data structure representing blocks
blocks.clear();
blocks.resize(m);
for(int i = 0; i < n; i++)
blocks[unordered2ordered[endo2block[i]]].insert(i);
#ifdef DEBUG
cout << "Found " << m << " blocks" << endl;
for(int i = 0; i < m; i++)
cout << " Block " << i << " of size " << blocks[i].size() << endl;
#endif
}
void
StaticModel::computeMFS(vector<set<int> > &blocksMFS, const vector<set<int> > &blocks, const vector<int> &endo2eq) const
{
const int n = equation_number();
assert((int) endo2eq.size() == n);
const int nblocks = blocks.size();
blocksMFS.clear();
blocksMFS.resize(nblocks);
// Iterate over blocks
for(int b = 0; b < nblocks; b++)
{
// Construct subgraph for MFS computation, where vertex number is position in the block
int p = blocks[b].size();
MFS::AdjacencyList_type g(p);
// Construct v_index and v_index1 properties, and a mapping between type specific IDs and vertex descriptors
property_map<MFS::AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, g);
property_map<MFS::AdjacencyList_type, vertex_index1_t>::type v_index1 = get(vertex_index1, g);
map<int, graph_traits<MFS::AdjacencyList_type>::vertex_descriptor> tsid2vertex;
int j = 0;
for(set<int>::const_iterator it = blocks[b].begin(); it != blocks[b].end(); ++it)
{
tsid2vertex[*it] = vertex(j, g);
put(v_index, vertex(j, g), *it);
put(v_index1, vertex(j, g), *it);
j++;
}
// Add edges, loop over endogenous in the block
set<pair<int, int> > endo;
for(set<int>::const_iterator it = blocks[b].begin(); it != blocks[b].end(); ++it)
{
endo.clear();
// Test if associated equation is in normalized form, and compute set of endogenous appearing in it
ExprNode *lhs = equations[endo2eq[*it]]->get_arg1();
VariableNode *lhs_var = dynamic_cast<VariableNode *>(lhs);
if (lhs_var == NULL || lhs_var->get_symb_id() != symbol_table.getID(eEndogenous, *it))
lhs->collectEndogenous(endo); // Only collect endogenous of LHS if not normalized form
ExprNode *rhs = equations[endo2eq[*it]]->get_arg2();
rhs->collectEndogenous(endo);
for(set<pair<int, int> >::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);
}
}
// Compute minimum feedback set
MFS::Minimal_set_of_feedback_vertex(blocksMFS[b], g);
cout << "Block " << b << ": " << blocksMFS[b].size() << "/" << blocks[b].size() << " in MFS" << endl;
}
}
...@@ -38,20 +38,36 @@ private: ...@@ -38,20 +38,36 @@ private:
virtual int computeDerivID(int symb_id, int lag); virtual int computeDerivID(int symb_id, int lag);
//! Computes normalization of the static model //! Computes normalization of the static model
/*! Maps each endogenous type specific ID to the equation which defines it */ /*! Maps each endogenous type specific ID to the equation to which it is associated */
void computeNormalization(vector<int> &endo_to_eq) const; void computeNormalization(vector<int> &endo2eq) const;
//! Computes blocks of the static model, sorted in topological order
/*!
\param blocks ordered list of blocks, each one containing a list of type specific IDs of endogenous variables
\param endo2eq chosen normalization (mapping from type specific IDs of endogenous to equations)
*/
void computeSortedBlockDecomposition(vector<set<int> > &blocks, const vector<int> &endo2eq) const;
//! For each block of the static model, computes minimum feedback set (MFS)
/*!
\param blocksMFS for each block, contains the subset of type specific IDs which are in the MFS
\param blocks ordered list of blocks, each one containing a list of type specific IDs of endogenous variables
\param endo2eq chosen normalization (mapping from type specific IDs of endogenous to equations)
*/
void computeMFS(vector<set<int> > &blocksMFS, const vector<set<int> > &blocks, const vector<int> &endo2eq) const;
//! Computes the list of equations which are already in normalized form //! Computes the list of equations which are already in normalized form
/*! Returns a multimap mapping endogenous which are normalized (represented by their type specific ID) to the equation(s) which define it */ /*! Returns a multimap mapping endogenous which are normalized (represented by their type specific ID) to the equation(s) which define it */
void computeNormalizedEquations(multimap<int, int> &endo_to_eqs) const; void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
public: public:
StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants); StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Execute computations (derivation) //! Execute computations (derivation)
/*! You must set computeStaticHessian before calling this function /*!
\param block_mfs whether block decomposition and minimum feedback set should be computed
\param hessian whether Hessian (w.r. to endogenous only) should be computed \param hessian whether Hessian (w.r. to endogenous only) should be computed
\param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */ \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
void computingPass(bool hessian, bool no_tmp_terms); void computingPass(bool block_mfs, bool hessian, bool no_tmp_terms);
//! Writes static model file //! Writes static model file
void writeStaticFile(const string &basename) const; void writeStaticFile(const string &basename) const;
......
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