From 7f57821401f84668e44f0ae17a4220b68b562a45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 29 Apr 2020 18:48:42 +0200 Subject: [PATCH] Now compute blocks[].first_equation from ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock Also include various cosmetic changes. --- src/ModelTree.cc | 111 ++++++++++++++++++++++------------------------- src/ModelTree.hh | 4 +- 2 files changed, 53 insertions(+), 62 deletions(-) diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 37ebf56c..01f4ddaf 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -673,7 +673,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock() /* Identify the simultaneous blocks. Each simultaneous block is given an index, starting from 0, in recursive order */ - auto [num_simblocks, endo2simblock] = G.sortedStronglyConnectedComponents(); + auto [num_simblocks, simvar2simblock] = G.sortedStronglyConnectedComponents(); int num_blocks = prologue+num_simblocks+epilogue; @@ -683,29 +683,23 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock() eq2block.resize(nb_var); // Initialize size and mfs_size for prologue and epilogue, plus eq/endo→block mappings - for (int i = 0; i < prologue; i++) - { - blocks[i].size = 1; - blocks[i].mfs_size = 1; - endo2block[endo_idx_block2orig[i]] = i; - eq2block[eq_idx_block2orig[i]] = i; - } - for (int i = 0; i < epilogue; i++) - { - int blk = prologue+num_simblocks+i; - int var_eq = prologue+nb_simvars+i; - blocks[blk].size = 1; - blocks[blk].mfs_size = 1; - endo2block[endo_idx_block2orig[var_eq]] = blk; - eq2block[eq_idx_block2orig[var_eq]] = blk; - } + for (int blk = 0; blk < num_blocks; blk++) + if (blk < prologue || blk >= num_blocks-epilogue) + { + int var_eq = (blk < prologue ? blk : blk-num_simblocks+nb_simvars); + blocks[blk].size = 1; + blocks[blk].mfs_size = 1; + blocks[blk].first_equation = var_eq; + endo2block[endo_idx_block2orig[var_eq]] = blk; + eq2block[eq_idx_block2orig[var_eq]] = blk; + } - // Compute size and list of equations for simultaneous blocks - vector> eqs_in_simblock(num_simblocks); - for (int i = 0; i < static_cast(endo2simblock.size()); i++) + // Initialize size for simultaneous blocks, plus eq/endo→block mappings + vector> simblock2simvars(num_simblocks); + for (int i = 0; i < static_cast(simvar2simblock.size()); i++) { - eqs_in_simblock[endo2simblock[i]].push_back(i); - int blk = prologue+endo2simblock[i]; + simblock2simvars[simvar2simblock[i]].push_back(i); + int blk = prologue+simvar2simblock[i]; blocks[blk].size++; endo2block[endo_idx_block2orig[prologue+i]] = blk; eq2block[eq_idx_block2orig[prologue+i]] = blk; @@ -749,11 +743,12 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock() const vector old_eq_idx_block2orig(eq_idx_block2orig), old_endo_idx_block2orig(endo_idx_block2orig); int ordidx = prologue; - for (int i = 0; i < num_simblocks; i++) + for (int blk = prologue; blk < prologue+num_simblocks; blk++) { - auto subG = G.extractSubgraph(eqs_in_simblock[i]); + blocks[blk].first_equation = (blk == 0 ? 0 : blocks[blk-1].first_equation + blocks[blk-1].size); + auto subG = G.extractSubgraph(simblock2simvars[blk-prologue]); auto feed_back_vertices = subG.minimalSetOfFeedbackVertices(); - blocks[prologue+i].mfs_size = feed_back_vertices.size(); + blocks[blk].mfs_size = feed_back_vertices.size(); auto recursive_vertices = subG.reorderRecursiveVariables(feed_back_vertices); auto v_index1 = get(boost::vertex_index1, subG); @@ -816,19 +811,17 @@ ModelTree::printBlockDecomposition() const void ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) { - for (int i = 0; i < static_cast(blocks.size()); i++) + for (int blk = 0; blk < static_cast(blocks.size()); blk++) { - int first_eq = (i == 0) ? 0 : blocks[i-1].first_equation+blocks[i-1].size; - /* Compute the maximum lead and lag across all endogenous that appear in this block and that belong to it */ int max_lag = 0, max_lead = 0; - for (int eq = first_eq; eq < first_eq+blocks[i].size; eq++) + for (int eq = 0; eq < blocks[blk].size; eq++) { set> endos_and_lags; - equations[eq_idx_block2orig[eq]]->collectEndogenous(endos_and_lags); + getBlockEquationExpr(blk, eq)->collectEndogenous(endos_and_lags); for (const auto &[endo, lag] : endos_and_lags) - if (linear_decomposition || endo2block[endo] == i) + if (linear_decomposition || endo2block[endo] == blk) { max_lead = max(lag, max_lead); max_lag = max(-lag, max_lag); @@ -839,12 +832,12 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) BlockSimulationType Simulation_Type; if (max_lag > 0 && max_lead > 0) { - if (blocks[i].size == 1) + if (blocks[blk].size == 1) Simulation_Type = BlockSimulationType::solveTwoBoundariesSimple; else Simulation_Type = BlockSimulationType::solveTwoBoundariesComplete; } - else if (blocks[i].size > 1) + else if (blocks[blk].size > 1) { if (max_lead > 0) Simulation_Type = BlockSimulationType::solveBackwardComplete; @@ -859,11 +852,11 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) Simulation_Type = BlockSimulationType::solveForwardSimple; } - if (blocks[i].size == 1) + if (blocks[blk].size == 1) { // Determine if the block can simply be evaluated - if (equation_type_and_normalized_equation[eq_idx_block2orig[first_eq]].first == EquationType::evaluate - || equation_type_and_normalized_equation[eq_idx_block2orig[first_eq]].first == EquationType::evaluate_s) + if (getBlockEquationType(blk, 0) == EquationType::evaluate + || getBlockEquationType(blk, 0) == EquationType::evaluate_s) { if (Simulation_Type == BlockSimulationType::solveBackwardSimple) Simulation_Type = BlockSimulationType::evaluateBackward; @@ -875,19 +868,18 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) This is only possible if the two blocks can simply be evaluated (in the same direction), and if the merge does not break the restrictions on leads/lags. */ - if (i > 0) + if (blk > 0) { + set> endos_and_lags; + getBlockEquationExpr(blk, 0)->collectEndogenous(endos_and_lags); bool is_lead = false, is_lag = false; - for (int j = blocks[i-1].first_equation; - j < blocks[i-1].first_equation+blocks[i-1].size; j++) + for (int var = 0; var < blocks[blk-1].size; var++) { - set> endos_and_lags; - equations[eq_idx_block2orig[first_eq]]->collectEndogenous(endos_and_lags); - is_lag = endos_and_lags.find({ endo_idx_block2orig[j], -1 }) != endos_and_lags.end(); - is_lead = endos_and_lags.find({ endo_idx_block2orig[j], 1 }) != endos_and_lags.end(); + is_lag = endos_and_lags.find({ getBlockVariableID(blk-1, var), -1 }) != endos_and_lags.end(); + is_lead = endos_and_lags.find({ getBlockVariableID(blk-1, var), 1 }) != endos_and_lags.end(); } - BlockSimulationType prev_Type = blocks[i-1].simulation_type; + BlockSimulationType prev_Type = blocks[blk-1].simulation_type; if ((prev_Type == BlockSimulationType::evaluateForward && Simulation_Type == BlockSimulationType::evaluateForward && !is_lead) @@ -896,35 +888,34 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition) && !is_lag)) { // Merge the current block into the previous one - blocks[i-1].size++; - blocks[i-1].mfs_size = blocks[i-1].size; + blocks[blk-1].size++; + blocks[blk-1].mfs_size = blocks[blk-1].size; /* For max lag/lead, the max of the two blocks is not enough. We need to consider the case where a variable of the previous block appears with a lag/lead in the current one (the reverse case is excluded, by construction). */ - blocks[i-1].max_lag = is_lag ? 1 : max(blocks[i-1].max_lag, max_lag); - blocks[i-1].max_lead = is_lead ? 1 : max(blocks[i-1].max_lead, max_lead); - blocks[i-1].n_static += blocks[i].n_static; - blocks[i-1].n_forward += blocks[i].n_forward; - blocks[i-1].n_backward += blocks[i].n_backward; - blocks[i-1].n_mixed += blocks[i].n_mixed; - blocks.erase(blocks.begin()+i); + blocks[blk-1].max_lag = is_lag ? 1 : max(blocks[blk-1].max_lag, max_lag); + blocks[blk-1].max_lead = is_lead ? 1 : max(blocks[blk-1].max_lead, max_lead); + blocks[blk-1].n_static += blocks[blk].n_static; + blocks[blk-1].n_forward += blocks[blk].n_forward; + blocks[blk-1].n_backward += blocks[blk].n_backward; + blocks[blk-1].n_mixed += blocks[blk].n_mixed; + blocks.erase(blocks.begin()+blk); for (auto &b : endo2block) - if (b >= i) + if (b >= blk) b--; for (auto &b : eq2block) - if (b >= i) + if (b >= blk) b--; - i--; + blk--; continue; } } } - blocks[i].simulation_type = Simulation_Type; - blocks[i].first_equation = first_eq; - blocks[i].max_lag = max_lag; - blocks[i].max_lead = max_lead; + blocks[blk].simulation_type = Simulation_Type; + blocks[blk].max_lag = max_lag; + blocks[blk].max_lead = max_lead; } } diff --git a/src/ModelTree.hh b/src/ModelTree.hh index 87acc746..e33b3aa7 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -301,7 +301,7 @@ protected: 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 - Initializes the “blocks” structure, and fills the following fields: size, + Initializes the “blocks” structure, and fills the following fields: size, first_equation, mfs_size, n_static, n_forward, n_backward, n_mixed. Also initializes the endo2block and eq2block structures. */ void computeBlockDecompositionAndFeedbackVariablesForEachBlock(); @@ -309,7 +309,7 @@ protected: prologue and the epilogue, and determine the type of each block. Fills the following fields of the “blocks” structure: simulation_type, - first_equation, max_lead, max_lag. */ + max_lead, max_lag. */ void reduceBlocksAndTypeDetermination(bool linear_decomposition); /* 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