diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 428592c9d..ce494e24a 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -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 SteadyStatement::writeOutput(ostream &output, const string &basename) const { @@ -904,7 +911,7 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct) void PlannerObjectiveStatement::computingPass() { - model_tree->computingPass(true, false); + model_tree->computingPass(false, true, false); } void diff --git a/preprocessor/ComputingTasks.hh b/preprocessor/ComputingTasks.hh index 7f52eaa87..bf0e31d86 100644 --- a/preprocessor/ComputingTasks.hh +++ b/preprocessor/ComputingTasks.hh @@ -34,6 +34,7 @@ private: const OptionsList options_list; public: SteadyStatement(const OptionsList &options_list_arg); + virtual void checkPass(ModFileStructure &mod_file_struct); virtual void writeOutput(ostream &output, const string &basename) const; }; diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 1f45870ac..b868694b2 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -87,7 +87,7 @@ class ParsingDriver; %} %token AR AUTOCORR -%token BAYESIAN_IRF BETA_PDF BICGSTAB +%token BAYESIAN_IRF BETA_PDF BICGSTAB BLOCK_MFS %token BVAR_DENSITY BVAR_FORECAST %token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA %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 steady_options : o_solve_algo | o_homotopy_mode | o_homotopy_steps + | o_block_mfs ; check : CHECK ';' @@ -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_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 { diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 5b0611b69..013e116ff 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -216,6 +216,7 @@ int sigma_e = 0; filename {return token::FILENAME;} diffuse_filter {return token::DIFFUSE_FILTER;} plot_priors {return token::PLOT_PRIORS;} +block_mfs {return token::BLOCK_MFS;} /* These four (var, varexo, varexo_det, parameters) are for change_type */ var { return token::VAR; } diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index d081a1e10..c66d53125 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -144,7 +144,7 @@ ModFile::computingPass(bool no_tmp_terms) { // Compute static model and its derivatives 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 diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 35df9360d..4350af575 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -30,7 +30,8 @@ ModFileStructure::ModFileStructure() : order_option(0), bvar_density_present(false), bvar_forecast_present(false), - identification_present(false) + identification_present(false), + steady_block_mfs_option(false) { } diff --git a/preprocessor/Statement.hh b/preprocessor/Statement.hh index 4bcba3c54..57c234e51 100644 --- a/preprocessor/Statement.hh +++ b/preprocessor/Statement.hh @@ -56,6 +56,8 @@ public: bool bvar_forecast_present; //! Whether an identification statement is present bool identification_present; + //! Whether the option "block_mfs" is used on steady statement + bool steady_block_mfs_option; }; class Statement diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index 1859b78d2..eb48ed312 100644 --- a/preprocessor/StaticModel.cc +++ b/preprocessor/StaticModel.cc @@ -20,18 +20,22 @@ #include #include +#include #include #include -/* -#include +#ifdef DEBUG +# include using namespace __gnu_cxx; -*/ +#endif #include #include +#include +#include #include "StaticModel.hh" +#include "MinimumFeedbackSet.hh" using namespace boost; @@ -298,7 +302,7 @@ StaticModel::writeStaticFile(const string &basename) const } 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 set vars; @@ -319,24 +323,17 @@ StaticModel::computingPass(bool hessian, bool no_tmp_terms) if (!no_tmp_terms) computeTemporaryTerms(); - /* - vector endo2eq(equation_number()); - computeNormalization(endo2eq); - - multimap natural_endo2eqs; - computeNormalizedEquations(natural_endo2eqs); - - for(int i = 0; i < symbol_table.endo_nbr(); i++) + if (block_mfs) { - if (natural_endo2eqs.count(i) == 0) - continue; + vector endo2eq(equation_number()); + computeNormalization(endo2eq); - pair::const_iterator, multimap::const_iterator> x = natural_endo2eqs.equal_range(i); - if (find_if(x.first, x.second, compose1(bind2nd(equal_to(), endo2eq[i]), select2nd::value_type>())) == x.second) - cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i)) - << " not used." << endl; + vector > blocks; + computeSortedBlockDecomposition(blocks, endo2eq); + + vector > blocksMFS; + computeMFS(blocksMFS, blocks, endo2eq); } - */ } int @@ -358,9 +355,9 @@ StaticModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDExcepti } void -StaticModel::computeNormalization(vector &endo_to_eq) const +StaticModel::computeNormalization(vector &endo2eq) const { - int n = equation_number(); + const int n = equation_number(); assert(n == symbol_table.endo_nbr()); @@ -384,15 +381,42 @@ StaticModel::computeNormalization(vector &endo_to_eq) const } // Compute maximum cardinality matching - typedef vector::vertex_descriptor> mate_map_t; - mate_map_t mate_map(2*n); + vector mate_map(2*n); +#if 1 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::null_vertex()); + + multimap 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::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::type>::verify_matching(g, &mate_map[0], get(vertex_index, g)); +#endif assert(check); // Check if all variables are normalized - mate_map_t::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits::null_vertex()); + vector::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits::null_vertex()); if (it != mate_map.begin() + n) { cerr << "ERROR: Could not normalize static model. Variable " @@ -401,18 +425,44 @@ StaticModel::computeNormalization(vector &endo_to_eq) const exit(EXIT_FAILURE); } +#ifdef DEBUG for(int i = 0; i < n; 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 - transform(mate_map.begin(), mate_map.begin() + n, endo_to_eq.begin(), bind2nd(minus(), n)); + transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus(), n)); + +#ifdef DEBUG + multimap 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::const_iterator, multimap::const_iterator> x = natural_endo2eqs.equal_range(i); + if (find_if(x.first, x.second, compose1(bind2nd(equal_to(), endo2eq[i]), select2nd::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 -StaticModel::computeNormalizedEquations(multimap &endo_to_eqs) const +StaticModel::computeNormalizedEquations(multimap &endo2eqs) const { for(int i = 0; i < equation_number(); i++) { @@ -429,7 +479,7 @@ StaticModel::computeNormalizedEquations(multimap &endo_to_eqs) const if (endo.find(make_pair(symb_id, 0)) != endo.end()) 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; } } @@ -439,3 +489,120 @@ StaticModel::writeLatexFile(const string &basename) const { writeLatexModelFile(basename + "_static.tex", oLatexStaticModel); } + +void +StaticModel::computeSortedBlockDecomposition(vector > &blocks, const vector &endo2eq) const +{ + const int n = equation_number(); + + assert((int) endo2eq.size() == n); + + // Compute graph representation of static model + typedef adjacency_list DirectedGraph; + DirectedGraph g(n); + + set > endo; + for(int i = 0; i < n; i++) + { + endo.clear(); + 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); + } + + // Compute strongly connected components + vector endo2block(n); + int m = strong_components(g, &endo2block[0]); + + // Create directed acyclic graph associated to the strongly connected components + DirectedGraph dag(m); + graph_traits::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 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 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 > &blocksMFS, const vector > &blocks, const vector &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::type v_index = get(vertex_index, g); + property_map::type v_index1 = get(vertex_index1, g); + map::vertex_descriptor> tsid2vertex; + int j = 0; + for(set::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 > endo; + for(set::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(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 >::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; + } +} diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh index 30302159b..1b2ce7ca8 100644 --- a/preprocessor/StaticModel.hh +++ b/preprocessor/StaticModel.hh @@ -38,20 +38,36 @@ private: virtual int computeDerivID(int symb_id, int lag); //! Computes normalization of the static model - /*! Maps each endogenous type specific ID to the equation which defines it */ - void computeNormalization(vector &endo_to_eq) const; + /*! Maps each endogenous type specific ID to the equation to which it is associated */ + void computeNormalization(vector &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 > &blocks, const vector &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 > &blocksMFS, const vector > &blocks, const vector &endo2eq) const; //! 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 */ - void computeNormalizedEquations(multimap &endo_to_eqs) const; + void computeNormalizedEquations(multimap &endo2eqs) const; public: StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants); //! 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 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 void writeStaticFile(const string &basename) const;