From 8eafd9ab4f1b46300eef09df47eccafdae012dae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Fri, 17 Apr 2020 19:21:37 +0200 Subject: [PATCH] Block decomposition: various refactorings --- src/DynamicModel.cc | 8 +-- src/ModelTree.cc | 123 +++++++++++++-------------------- src/ModelTree.hh | 22 +++--- src/StaticModel.cc | 2 +- src/VariableDependencyGraph.cc | 6 +- src/VariableDependencyGraph.hh | 6 +- 6 files changed, 70 insertions(+), 97 deletions(-) diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 34c7b21f..e0308aab 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -4809,9 +4809,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO if (!computeNaturalNormalization()) computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian); - lag_lead_vector_t equation_lag_lead, variable_lag_lead; - - tie(simblock_size, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) + tie(simblock_size, n_static, n_forward, n_backward, n_mixed) = select_non_linear_equations_and_variables(is_equation_linear); equationTypeDetermination(first_order_endo_derivatives, 0); @@ -4847,9 +4845,9 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO cout << "Finding the optimal block decomposition of the model ..." << endl; - lag_lead_vector_t equation_lag_lead, variable_lag_lead; + lag_lead_vector_t variable_lag_lead; - tie(simblock_size, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false, true); + tie(simblock_size, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false); reduceBlocksAndTypeDetermination(simblock_size, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition); diff --git a/src/ModelTree.cc b/src/ModelTree.cc index cb43b463..db2c35c4 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -456,8 +456,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double return { contemporaneous_jacobian, static_jacobian }; } -tuple>, lag_lead_vector_t, lag_lead_vector_t, - vector, vector, vector, vector> +tuple>, vector, vector, vector, vector> ModelTree::select_non_linear_equations_and_variables(const vector &is_equation_linear) { vector eq2endo(equations.size(), 0); @@ -476,7 +475,7 @@ ModelTree::select_non_linear_equations_and_variables(const vector &is_equa i++; j++; } - auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block, endo2block.size()); + auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block); vector n_static(endo2eq.size(), 0), n_forward(endo2eq.size(), 0), n_backward(endo2eq.size(), 0), n_mixed(endo2eq.size(), 0); for (int i = 0; i < static_cast(endo2eq.size()); i++) @@ -493,8 +492,7 @@ ModelTree::select_non_linear_equations_and_variables(const vector &is_equa cout.flush(); vector> simblock_size(1, {i, i}); updateReverseVariableEquationOrderings(); - return { simblock_size, equation_lag_lead, variable_lag_lead, - n_static, n_forward, n_backward, n_mixed }; + return { simblock_size, n_static, n_forward, n_backward, n_mixed }; } bool @@ -530,18 +528,14 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian) { const int n = equations.size(); - /* Compute reverse map (eq→endo) of normalization. Also initialize - “eq_idx_block2orig” and “endo_idx_block2orig” to the identity + /* Initialize “eq_idx_block2orig” and “endo_idx_block2orig” to the identity permutation. */ - vector eq2endo(n); eq_idx_block2orig.resize(n); endo_idx_block2orig.resize(n); for (int i = 0; i < n; i++) { - int it = endo2eq[i]; - eq2endo[it] = i; eq_idx_block2orig[i] = i; - endo_idx_block2orig[it] = i; + endo_idx_block2orig[endo2eq[i]] = i; } /* Compute incidence matrix, equations in rows, variables in columns. Row @@ -677,53 +671,42 @@ ModelTree::equationTypeDetermination(const map, expr_t> &fi } pair -ModelTree::getVariableLeadLagByBlock(const vector &endo2simblock, int num_simblocks) const +ModelTree::getVariableLeadLagByBlock(const vector &endo2simblock) const { int nb_endo = symbol_table.endo_nbr(); - lag_lead_vector_t variable_lead_lag(nb_endo, { 0, 0 }), equation_lead_lag(nb_endo, { 0, 0 }); - vector variable_blck(nb_endo), equation_blck(nb_endo); - for (int i = 0; i < nb_endo; i++) - { - if (i < prologue) - { - variable_blck[endo_idx_block2orig[i]] = i; - equation_blck[eq_idx_block2orig[i]] = i; - } - else if (i < static_cast(endo2simblock.size()) + prologue) - { - variable_blck[endo_idx_block2orig[i]] = endo2simblock[i-prologue] + prologue; - equation_blck[eq_idx_block2orig[i]] = endo2simblock[i-prologue] + prologue; - } - else - { - variable_blck[endo_idx_block2orig[i]] = i - (nb_endo - num_simblocks - prologue - epilogue); - equation_blck[eq_idx_block2orig[i]] = i - (nb_endo - num_simblocks - prologue - epilogue); - } - } + + auto belong_to_same_block = [&](int endo, int eq) + { + int endo2 = endo_idx_orig2block[endo]; + int eq2 = eq_idx_orig2block[eq]; + if (endo2 < prologue || endo2 >= nb_endo-epilogue + || eq2 < prologue || eq2 >= nb_endo-epilogue) + return endo2 == eq2; + else + return endo2simblock[endo2-prologue] == endo2simblock[eq2-prologue]; + }; + + lag_lead_vector_t variable_lag_lead(nb_endo, { 0, 0 }), equation_lag_lead(nb_endo, { 0, 0 }); for (const auto &[key, value] : dynamic_jacobian) { - auto [lag, j_1, i_1] = key; - if (variable_blck[i_1] == equation_blck[j_1]) + auto [lag, eq, endo] = key; + if (belong_to_same_block(endo, eq)) { - if (lag > variable_lead_lag[i_1].second) - variable_lead_lag[i_1] = { variable_lead_lag[i_1].first, lag }; - if (lag < -variable_lead_lag[i_1].first) - variable_lead_lag[i_1] = { -lag, variable_lead_lag[i_1].second }; - if (lag > equation_lead_lag[j_1].second) - equation_lead_lag[j_1] = { equation_lead_lag[j_1].first, lag }; - if (lag < -equation_lead_lag[j_1].first) - equation_lead_lag[j_1] = { -lag, equation_lead_lag[j_1].second }; + variable_lag_lead[endo].first = max(variable_lag_lead[endo].first, -lag); + variable_lag_lead[endo].second = max(variable_lag_lead[endo].second, lag); + equation_lag_lead[eq].first = max(equation_lag_lead[eq].first, -lag); + equation_lag_lead[eq].second = max(equation_lag_lead[eq].second, lag); } } - return { equation_lead_lag, variable_lead_lag }; + return { equation_lag_lead, variable_lag_lead }; } -tuple>, lag_lead_vector_t, lag_lead_vector_t, +tuple>, lag_lead_vector_t, vector, vector, vector, vector> -ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable) +ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_) { - int nb_var = endo_idx_block2orig.size(); - int n = nb_var - prologue - epilogue; + int nb_var = symbol_table.endo_nbr(); + int nb_simvars = nb_var - prologue - epilogue; /* Construct the graph representing the dependencies between all variables that do not belong to the prologue or the epilogue. @@ -731,7 +714,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob For detecting dependencies between variables, use the static jacobian, except when the cutoff is zero, in which case use the symbolic adjacency matrix */ - VariableDependencyGraph G(n); + VariableDependencyGraph G(nb_simvars); for (const auto &[key, value] : cutoff == 0 ? computeSymbolicJacobian() : static_jacobian) { auto [eq, endo] = key; @@ -758,25 +741,18 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob eqs_in_simblock[endo2simblock[i]].insert(i); } - auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2simblock, num_simblocks); + auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2simblock); /* Add a loop on vertices which could not be normalized or vertices related - to lead variables. This forces those vertices to belong to the feedback set */ - if (select_feedback_variable) - { - for (int i = 0; i < n; i++) - if (Equation_Type[eq_idx_block2orig[i+prologue]].first == EquationType::solve - || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0 - || variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0 - || equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0 - || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0 - || mfs == 0) - add_edge(vertex(i, G), vertex(i, G), G); - } - else - for (int i = 0; i < n; i++) - if (Equation_Type[eq_idx_block2orig[i+prologue]].first == EquationType::solve || mfs == 0) - add_edge(vertex(i, G), vertex(i, G), G); + to lead/lag variables. This forces those vertices to belong to the feedback set */ + for (int i = 0; i < nb_simvars; i++) + if (Equation_Type[eq_idx_block2orig[i+prologue]].first == EquationType::solve + || variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0 + || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0 + || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0 + || equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0 + || mfs == 0) + add_edge(vertex(i, G), vertex(i, G), G); int num_blocks = prologue+num_simblocks+epilogue; // Determines the dynamic structure of each equation @@ -786,8 +762,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob const vector old_eq_idx_block2orig(eq_idx_block2orig), old_endo_idx_block2orig(endo_idx_block2orig); for (int i = 0; i < prologue; i++) { - int max_lag = variable_lag_lead[old_endo_idx_block2orig[i]].first; - int max_lead = variable_lag_lead[old_endo_idx_block2orig[i]].second; + auto [max_lag, max_lead] = variable_lag_lead[old_endo_idx_block2orig[i]]; if (max_lag != 0 && max_lead != 0) n_mixed[i]++; else if (max_lag == 0 && max_lead != 0) @@ -806,7 +781,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob for (int i = 0; i < num_simblocks; i++) { auto subG = G.extractSubgraph(eqs_in_simblock[i]); - auto [G1, feed_back_vertices] = subG.minimalSetOfFeedbackVertices(); + auto feed_back_vertices = subG.minimalSetOfFeedbackVertices(); auto v_index1 = get(boost::vertex_index1, subG); simblock_size[i].second = feed_back_vertices.size(); auto reordered_vertices = subG.reorderRecursiveVariables(feed_back_vertices); @@ -815,8 +790,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob for (int j = 0; j < 4; j++) for (int its : reordered_vertices) { - int max_lag = variable_lag_lead[old_endo_idx_block2orig[its+prologue]].first; - int max_lead = variable_lag_lead[old_endo_idx_block2orig[its+prologue]].second; + auto [max_lag, max_lead] = variable_lag_lead[old_endo_idx_block2orig[its+prologue]]; auto reorder = [&]() { eq_idx_block2orig[ordidx] = old_eq_idx_block2orig[its+prologue]; @@ -850,8 +824,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob for (int fbvertex : feed_back_vertices) { int idx = v_index1[vertex(fbvertex, subG)]; - int max_lag = variable_lag_lead[old_endo_idx_block2orig[idx+prologue]].first; - int max_lead = variable_lag_lead[old_endo_idx_block2orig[idx+prologue]].second; + auto [max_lag, max_lead] = variable_lag_lead[old_endo_idx_block2orig[idx+prologue]]; auto reorder = [&]() { eq_idx_block2orig[ordidx] = old_eq_idx_block2orig[idx+prologue]; @@ -883,8 +856,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob for (int i = 0; i < epilogue; i++) { - int max_lag = variable_lag_lead[old_endo_idx_block2orig[prologue+n+i]].first; - int max_lead = variable_lag_lead[old_endo_idx_block2orig[prologue+n+i]].second; + auto [max_lag, max_lead] = variable_lag_lead[old_endo_idx_block2orig[prologue+nb_simvars+i]]; if (max_lag != 0 && max_lead != 0) n_mixed[prologue+num_simblocks+i]++; else if (max_lag == 0 && max_lead != 0) @@ -897,8 +869,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob updateReverseVariableEquationOrderings(); - return { simblock_size, equation_lag_lead, variable_lag_lead, - n_static, n_forward, n_backward, n_mixed }; + return { simblock_size, variable_lag_lead, n_static, n_forward, n_backward, n_mixed }; } void diff --git a/src/ModelTree.hh b/src/ModelTree.hh index 99127e3d..ea722a4c 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -302,22 +302,26 @@ protected: dynamic_jacobian. Elements below the cutoff are discarded. External functions are evaluated to 1. */ pair evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose); //! Select and reorder the non linear equations of the model - /*! Returns a tuple (blocks, equation_lag_lead, variable_lag_lead, n_static, - n_forward, n_backward, n_mixed) */ - tuple>, lag_lead_vector_t, lag_lead_vector_t, vector, vector, vector, vector> select_non_linear_equations_and_variables(const vector &is_equation_linear); + /*! Returns a tuple (blocks, n_static, n_forward, n_backward, n_mixed) */ + tuple>, vector, vector, vector, vector> select_non_linear_equations_and_variables(const vector &is_equation_linear); //! Search the equations and variables belonging to the prologue and the epilogue of the model void computePrologueAndEpilogue(const jacob_map_t &static_jacobian); //! Determine the type of each equation of model and try to normalize the unnormalized equation void equationTypeDetermination(const map, expr_t> &first_order_endo_derivatives, int mfs); //! Compute the block decomposition and for a non-recusive block find the minimum feedback set - /*! Returns a tuple (blocks, equation_lag_lead, variable_lag_lead, n_static, - n_forward, n_backward, n_mixed) */ - tuple>, lag_lead_vector_t, lag_lead_vector_t, vector, vector, vector, vector> computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable); + /*! Returns a tuple (blocks, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) */ + tuple>, lag_lead_vector_t, vector, vector, vector, vector> computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_); //! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block void reduceBlocksAndTypeDetermination(const vector> &simblock_size, const equation_type_and_normalized_equation_t &Equation_Type, const vector &n_static, const vector &n_forward, const vector &n_backward, const vector &n_mixed, bool linear_decomposition); - //! Determine the maximum number of lead and lag for the endogenous variable in a bloc - /*! Returns a pair { equation_lead,lag, variable_lead_lag } */ - pair getVariableLeadLagByBlock(const vector &endo2simblock, int num_simblocks) const; + /* The 1st output gives, for each equation (in original order) the (max_lag, + max_lead) across all endogenous that appear in the equation and that + belong to the same block (i.e. those endogenous are solved in the same + block). + + The 2nd output gives, for each type-specific endo IDs, its (max_lag, + max_lead) across all its occurences inside the equations of the block to + which it belongs. */ + pair getVariableLeadLagByBlock(const vector &endo2simblock) const; //! For each equation determine if it is linear or not vector equationLinear(const map, expr_t> &first_order_endo_derivatives) const; //! Print an abstract of the block structure of the model diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 616ec23f..eceaa37d 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -1142,7 +1142,7 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co cout << "Finding the optimal block decomposition of the model ..." << endl; - auto [blocks, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed] = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false, false); + auto [blocks, variable_lag_lead, n_static, n_forward, n_backward, n_mixed] = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false); reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, false); diff --git a/src/VariableDependencyGraph.cc b/src/VariableDependencyGraph.cc index c1de21f8..e0e951b6 100644 --- a/src/VariableDependencyGraph.cc +++ b/src/VariableDependencyGraph.cc @@ -320,7 +320,7 @@ VariableDependencyGraph::suppressionOfVerticesWithLoop(set &feed_back_verti return something_has_been_done; } -pair> +set VariableDependencyGraph::minimalSetOfFeedbackVertices() const { bool something_has_been_done = true; @@ -354,7 +354,7 @@ VariableDependencyGraph::minimalSetOfFeedbackVertices() const #ifdef verbose cout << "has_cycle=false\n"; #endif - return { G, feed_back_vertices }; + return feed_back_vertices; } if (num_vertices(G) > 0) { @@ -383,7 +383,7 @@ VariableDependencyGraph::minimalSetOfFeedbackVertices() const #ifdef verbose cout << "cut_=" << cut_ << "\n"; #endif - return { G, feed_back_vertices }; + return feed_back_vertices; } vector diff --git a/src/VariableDependencyGraph.hh b/src/VariableDependencyGraph.hh index 6974c020..668d2c18 100644 --- a/src/VariableDependencyGraph.hh +++ b/src/VariableDependencyGraph.hh @@ -56,7 +56,7 @@ public: */ VariableDependencyGraph extractSubgraph(const set &select_index) const; //! Return the feedback set - pair> minimalSetOfFeedbackVertices() const; + set minimalSetOfFeedbackVertices() const; //! Reorder the recursive variables /*! They appear first in a quasi triangular form and they are followed by the feedback variables */ vector reorderRecursiveVariables(const set &feedback_vertices) const; @@ -65,6 +65,8 @@ public: Returns the number of SCCs, and a mapping of vertex indices to sorted SCC indices. */ pair> sortedStronglyConnectedComponents() const; + // Print on stdout a description of the graph + void print() const; private: // Remove a vertex (including all edges to and from it); takes a vertex descriptor void suppress(vertex_descriptor vertex_to_eliminate); @@ -78,8 +80,6 @@ private: bool hasCycleDFS(vertex_descriptor u, color_t &color, vector &circuit_stack) const; // Determine whether the graph has a cycle bool hasCycle() const; - // Print on stdout a description of the graph - void print() const; bool vertexBelongsToAClique(vertex_descriptor vertex) const; bool eliminationOfVerticesWithOneOrLessIndegreeOrOutdegree(); bool eliminationOfVerticesBelongingToAClique();