diff --git a/src/CodeInterpreter.hh b/src/CodeInterpreter.hh
index 3331375ffc8ab7112d42fa3b3a28dcc19471de39..cbd5fc53e4dfb31c1aaf5b5de745bfe501f805dd 100644
--- a/src/CodeInterpreter.hh
+++ b/src/CodeInterpreter.hh
@@ -105,14 +105,6 @@ enum Tags
};
-enum class BlockType
- {
- simultans, //!< Simultaneous time separable block
- prologue, //!< Prologue block (one equation at the beginning, later merged)
- epilogue, //!< Epilogue block (one equation at the beginning, later merged)
- simultan //!< Simultaneous time unseparable block
- };
-
enum class EquationType
{
unknown, //!< Unknown equation type
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 6e1ed53cb69be10da0877efecdd6d0878e6b14a8..b7ceb6f1c7e9047f7c7b892c4048f92d5748815f 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -406,30 +406,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
else
output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl;
- BlockType block_type;
- if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
- || simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
- block_type = BlockType::simultan;
- else if (simulation_type == BlockSimulationType::solveForwardComplete
- || simulation_type == BlockSimulationType::solveBackwardComplete)
- block_type = BlockType::simultans;
- else if ((simulation_type == BlockSimulationType::solveForwardSimple
- || simulation_type == BlockSimulationType::solveBackwardSimple
- || simulation_type == BlockSimulationType::evaluateBackward
- || simulation_type == BlockSimulationType::evaluateForward)
- && blocks[block].first_equation < prologue)
- block_type = BlockType::prologue;
- else if ((simulation_type == BlockSimulationType::solveForwardSimple
- || simulation_type == BlockSimulationType::solveBackwardSimple
- || simulation_type == BlockSimulationType::evaluateBackward
- || simulation_type == BlockSimulationType::evaluateForward)
- && blocks[block].first_equation >= static_cast<int>(equations.size()) - epilogue)
- block_type = BlockType::epilogue;
- else
- block_type = BlockType::simultans;
+
output << " % ////////////////////////////////////////////////////////////////////////" << endl
- << " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
- << " //" << endl
+ << " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1
+ << " //" << endl
<< " % // Simulation type "
<< BlockSim(simulation_type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
@@ -4755,7 +4735,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
equationTypeDetermination(first_order_endo_derivatives, 0);
- reduceBlocksAndTypeDetermination();
+ reduceBlockDecomposition();
computeChainRuleJacobian();
@@ -4775,7 +4755,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
- computePrologueAndEpilogue();
+ auto [prologue, epilogue] = computePrologueAndEpilogue();
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
@@ -4783,9 +4763,9 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
cout << "Finding the optimal block decomposition of the model ..." << endl;
- computeBlockDecompositionAndFeedbackVariablesForEachBlock();
+ computeBlockDecomposition(prologue, epilogue);
- reduceBlocksAndTypeDetermination();
+ reduceBlockDecomposition();
printBlockDecomposition();
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 3a00a8ed3abd4a73b73d1a0b9586c7b61eae27d7..1b74c67ee9946a1409331f673ae9283576168fb4 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -161,8 +161,6 @@ ModelTree::ModelTree(const ModelTree &m) :
eq2block{m.eq2block},
is_equation_linear{m.is_equation_linear},
endo2eq{m.endo2eq},
- epilogue{m.epilogue},
- prologue{m.prologue},
cutoff{m.cutoff},
mfs{m.mfs}
{
@@ -210,8 +208,6 @@ ModelTree::operator=(const ModelTree &m)
eq2block = m.eq2block;
is_equation_linear = m.is_equation_linear;
endo2eq = m.endo2eq;
- epilogue = m.epilogue;
- prologue = m.prologue;
cutoff = m.cutoff;
mfs = m.mfs;
@@ -403,9 +399,6 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double
void
ModelTree::select_non_linear_equations_and_variables()
{
- prologue = 0;
- epilogue = 0;
-
endo2block.resize(endo2eq.size(), 1); // The 1 is a dummy value, distinct from 0
eq2block.resize(endo2eq.size(), 1);
int i = 0;
@@ -429,6 +422,7 @@ ModelTree::select_non_linear_equations_and_variables()
blocks[0].mfs_size = i;
blocks[0].first_equation = 0;
computeDynamicStructureOfBlock(0);
+ computeSimulationTypeOfBlock(0);
}
bool
@@ -459,7 +453,7 @@ ModelTree::computeNaturalNormalization()
return bool_result;
}
-void
+pair<int, int>
ModelTree::computePrologueAndEpilogue()
{
const int n = equations.size();
@@ -489,7 +483,7 @@ ModelTree::computePrologueAndEpilogue()
bool something_has_been_done;
// Find the prologue equations and place first the AR(1) shock equations first
- prologue = 0;
+ int prologue = 0;
do
{
something_has_been_done = false;
@@ -525,7 +519,7 @@ ModelTree::computePrologueAndEpilogue()
while (something_has_been_done);
// Find the epilogue equations
- epilogue = 0;
+ int epilogue = 0;
do
{
something_has_been_done = false;
@@ -559,6 +553,8 @@ ModelTree::computePrologueAndEpilogue()
while (something_has_been_done);
updateReverseVariableEquationOrderings();
+
+ return { prologue, epilogue };
}
void
@@ -671,6 +667,37 @@ ModelTree::computeDynamicStructureOfBlock(int blk)
}
}
+void
+ModelTree::computeSimulationTypeOfBlock(int blk)
+{
+ auto &type = blocks[blk].simulation_type;
+ if (blocks[blk].max_endo_lag > 0 && blocks[blk].max_endo_lead > 0)
+ {
+ if (blocks[blk].size == 1)
+ type = BlockSimulationType::solveTwoBoundariesSimple;
+ else
+ type = BlockSimulationType::solveTwoBoundariesComplete;
+ }
+ else if (blocks[blk].size > 1)
+ {
+ if (blocks[blk].max_endo_lead > 0)
+ type = BlockSimulationType::solveBackwardComplete;
+ else
+ type = BlockSimulationType::solveForwardComplete;
+ }
+ else
+ {
+ bool can_eval = (getBlockEquationType(blk, 0) == EquationType::evaluate
+ || getBlockEquationType(blk, 0) == EquationType::evaluate_s);
+ if (blocks[blk].max_endo_lead > 0)
+ type = can_eval ? BlockSimulationType::evaluateBackward :
+ BlockSimulationType::solveBackwardSimple;
+ else
+ type = can_eval ? BlockSimulationType::evaluateForward :
+ BlockSimulationType::solveForwardSimple;
+ }
+}
+
pair<lag_lead_vector_t, lag_lead_vector_t>
ModelTree::getVariableLeadLagByBlock() const
{
@@ -694,7 +721,7 @@ ModelTree::getVariableLeadLagByBlock() const
}
void
-ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock()
+ModelTree::computeBlockDecomposition(int prologue, int epilogue)
{
int nb_var = symbol_table.endo_nbr();
int nb_simvars = nb_var - prologue - epilogue;
@@ -812,7 +839,10 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock()
updateReverseVariableEquationOrderings();
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
- computeDynamicStructureOfBlock(blk);
+ {
+ computeDynamicStructureOfBlock(blk);
+ computeSimulationTypeOfBlock(blk);
+ }
}
void
@@ -843,92 +873,46 @@ ModelTree::printBlockDecomposition() const
}
void
-ModelTree::reduceBlocksAndTypeDetermination()
+ModelTree::reduceBlockDecomposition()
{
- for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
- {
- // Determine the block type
- BlockSimulationType Simulation_Type;
- if (blocks[blk].max_endo_lag > 0 && blocks[blk].max_endo_lead > 0)
- {
- if (blocks[blk].size == 1)
- Simulation_Type = BlockSimulationType::solveTwoBoundariesSimple;
- else
- Simulation_Type = BlockSimulationType::solveTwoBoundariesComplete;
- }
- else if (blocks[blk].size > 1)
- {
- if (blocks[blk].max_endo_lead > 0)
- Simulation_Type = BlockSimulationType::solveBackwardComplete;
- else
- Simulation_Type = BlockSimulationType::solveForwardComplete;
- }
- else
- {
- if (blocks[blk].max_endo_lead > 0)
- Simulation_Type = BlockSimulationType::solveBackwardSimple;
- else
- Simulation_Type = BlockSimulationType::solveForwardSimple;
- }
-
- if (blocks[blk].size == 1)
- {
- // Determine if the block can simply be evaluated
- if (getBlockEquationType(blk, 0) == EquationType::evaluate
- || getBlockEquationType(blk, 0) == EquationType::evaluate_s)
- {
- if (Simulation_Type == BlockSimulationType::solveBackwardSimple)
- Simulation_Type = BlockSimulationType::evaluateBackward;
- else if (Simulation_Type == BlockSimulationType::solveForwardSimple)
- Simulation_Type = BlockSimulationType::evaluateForward;
- }
-
- /* Try to merge this block with the previous one.
- 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 (blk > 0)
- {
- set<pair<int, int>> endos_and_lags;
- getBlockEquationExpr(blk, 0)->collectEndogenous(endos_and_lags);
- bool is_lead = false, is_lag = false;
- for (int var = 0; var < blocks[blk-1].size; var++)
- {
- 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[blk-1].simulation_type;
- if ((prev_Type == BlockSimulationType::evaluateForward
- && Simulation_Type == BlockSimulationType::evaluateForward
- && !is_lead)
- || (prev_Type == BlockSimulationType::evaluateBackward
- && Simulation_Type == BlockSimulationType::evaluateBackward
- && !is_lag))
- {
- // Merge the current block into the previous one
- 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). */
- computeDynamicStructureOfBlock(blk-1);
- blocks.erase(blocks.begin()+blk);
- for (auto &b : endo2block)
- if (b >= blk)
- b--;
- for (auto &b : eq2block)
- if (b >= blk)
- b--;
- blk--;
- continue;
- }
- }
- }
+ for (int blk = 1; blk < static_cast<int>(blocks.size()); blk++)
+ if (blocks[blk].size == 1)
+ {
+ /* Try to merge this block with the previous one.
+ 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. */
+ set<pair<int, int>> endos_and_lags;
+ getBlockEquationExpr(blk, 0)->collectEndogenous(endos_and_lags);
+ bool is_lead = false, is_lag = false;
+ for (int var = 0; var < blocks[blk-1].size; var++)
+ {
+ 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();
+ }
- blocks[blk].simulation_type = Simulation_Type;
- }
+ if ((blocks[blk-1].simulation_type == BlockSimulationType::evaluateForward
+ && blocks[blk].simulation_type == BlockSimulationType::evaluateForward
+ && !is_lead)
+ || (blocks[blk-1].simulation_type == BlockSimulationType::evaluateBackward
+ && blocks[blk].simulation_type == BlockSimulationType::evaluateBackward
+ && !is_lag))
+ {
+ // Merge the current block into the previous one
+ blocks[blk-1].size++;
+ blocks[blk-1].mfs_size = blocks[blk-1].size;
+ computeDynamicStructureOfBlock(blk-1);
+ blocks.erase(blocks.begin()+blk);
+ for (auto &b : endo2block)
+ if (b >= blk)
+ b--;
+ for (auto &b : eq2block)
+ if (b >= blk)
+ b--;
+ blk--;
+ continue;
+ }
+ }
}
void
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index c29dbd22e2349b70484de3c5f9e9fb943e58a68f..95764aaa2e2717e3f476d776908629afab7b33e6 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -263,9 +263,6 @@ protected:
/*! Maps endogenous type specific IDs to equation numbers */
vector<int> endo2eq;
- //! number of equation in the prologue and in the epilogue
- int epilogue, prologue;
-
/* Compute a pseudo-Jacobian whose all elements are either zero or one,
depending on whether the variable symbolically appears in the equation */
jacob_map_t computeSymbolicJacobian() const;
@@ -296,25 +293,27 @@ protected:
pair<jacob_map_t, jacob_map_t> evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose) const;
//! Select and reorder the non linear equations of the model
void select_non_linear_equations_and_variables();
- //! Search the equations and variables belonging to the prologue and the epilogue of the model
- void computePrologueAndEpilogue();
+ /* Search the equations and variables belonging to the prologue and the
+ epilogue of the model.
+ Initializes “eq_idx_block2orig” and “endo_idx_block2orig”.
+ Returns the sizes of the prologue and epilogue. */
+ pair<int, int> computePrologueAndEpilogue();
//! Determine the type of each equation of model and try to normalize the unnormalized equation
void equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs);
/* Fills the max lags/leads and n_{static,mixed,forward,backward} fields of a
given block.
Needs the fields size and first_equation. */
void computeDynamicStructureOfBlock(int blk);
+ /* Fills the simulation_type field of a given block.
+ Needs the fields size, max_endo_lag and max_endo_lead. */
+ void computeSimulationTypeOfBlock(int blk);
/* 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, first_equation,
- mfs_size, n_static, n_forward, n_backward, n_mixed, maximum lags/leads.
- Also initializes the endo2block and eq2block structures. */
- void computeBlockDecompositionAndFeedbackVariablesForEachBlock();
+ Initializes the “blocks”, “endo2block” and “eq2block” structures. */
+ void computeBlockDecomposition(int prologue, int epilogue);
/* Reduce the number of block by merging the same type of equations in the
- prologue and the epilogue, and determine the type of each block.
-
- Fills the “simulation_type” field of the “blocks” structure. */
- void reduceBlocksAndTypeDetermination();
+ prologue and the epilogue */
+ void reduceBlockDecomposition();
/* 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
@@ -486,24 +485,6 @@ public:
}
}
- inline static string
- BlockType0(BlockType type)
- {
- switch (type)
- {
- case BlockType::simultans:
- return "SIMULTANEOUS TIME SEPARABLE ";
- case BlockType::prologue:
- return "PROLOGUE ";
- case BlockType::epilogue:
- return "EPILOGUE ";
- case BlockType::simultan:
- return "SIMULTANEOUS TIME UNSEPARABLE";
- default:
- return "UNKNOWN ";
- }
- }
-
inline static string
BlockSim(BlockSimulationType type)
{
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 82ecd946079968a06a2d4be5e6ace8469ba1c42a..aafd62905d30abc1c00cb8f59adbc6c952ac5a83 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -298,27 +298,9 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
else
output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)" << endl;
- BlockType block_type;
- if (simulation_type == BlockSimulationType::solveForwardComplete
- || simulation_type == BlockSimulationType::solveBackwardComplete)
- block_type = BlockType::simultans;
- else if ((simulation_type == BlockSimulationType::solveForwardSimple
- || simulation_type == BlockSimulationType::solveBackwardSimple
- || simulation_type == BlockSimulationType::evaluateBackward
- || simulation_type == BlockSimulationType::evaluateForward)
- && blocks[block].first_equation < prologue)
- block_type = BlockType::prologue;
- else if ((simulation_type == BlockSimulationType::solveForwardSimple
- || simulation_type == BlockSimulationType::solveBackwardSimple
- || simulation_type == BlockSimulationType::evaluateBackward
- || simulation_type == BlockSimulationType::evaluateForward)
- && blocks[block].first_equation >= static_cast<int>(equations.size()) - epilogue)
- block_type = BlockType::epilogue;
- else
- block_type = BlockType::simultans;
output << " % ////////////////////////////////////////////////////////////////////////" << endl
- << " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
- << " //" << endl
+ << " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1
+ << " //" << endl
<< " % // Simulation type "
<< BlockSim(simulation_type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
@@ -1125,7 +1107,7 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
- computePrologueAndEpilogue();
+ auto [prologue, epilogue] = computePrologueAndEpilogue();
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
@@ -1133,9 +1115,9 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
cout << "Finding the optimal block decomposition of the model ..." << endl;
- computeBlockDecompositionAndFeedbackVariablesForEachBlock();
+ computeBlockDecomposition(prologue, epilogue);
- reduceBlocksAndTypeDetermination();
+ reduceBlockDecomposition();
printBlockDecomposition();