Block decomposition: reorganize data structures storing block information

issue#70
Sébastien Villemot 2020-04-21 18:10:46 +02:00
parent 8eafd9ab4f
commit 118ceab85b
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
4 changed files with 221 additions and 254 deletions

View File

@ -1249,7 +1249,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
block_size,
endo_idx_block2orig,
eq_idx_block2orig,
blocks_linear[block],
blocks[block].linear,
symbol_table.endo_nbr(),
block_max_lag,
block_max_lead,
@ -2083,7 +2083,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
<< " end;" << endl
<< " y = solve_one_boundary('" << basename << ".block.dynamic_" << block + 1 << "'"
<< ", y, x, params, steady_state, y_index, " << nze
<< ", options_.periods, " << blocks_linear[block]
<< ", options_.periods, " << blocks[block].linear
<< ", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0, M_, options_, oo_);" << endl
<< " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
<< " if any(isnan(tmp) | isinf(tmp))" << endl
@ -2115,7 +2115,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
<< " end;" << endl
<< " y = solve_one_boundary('" << basename << ".block.dynamic_" << block + 1 << "'"
<<", y, x, params, steady_state, y_index, " << nze
<<", options_.periods, " << blocks_linear[block]
<<", options_.periods, " << blocks[block].linear
<<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0, M_, options_, oo_);" << endl
<< " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
<< " if any(isnan(tmp) | isinf(tmp))" << endl
@ -2147,7 +2147,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
<<", y, x, params, steady_state, y_index, " << nze
<<", options_.periods, " << max_leadlag_block[block].first
<<", " << max_leadlag_block[block].second
<<", " << blocks_linear[block]
<<", " << blocks[block].linear
<<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, options_, M_, oo_);" << endl
<< " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
<< " if any(isnan(tmp) | isinf(tmp))" << endl
@ -4793,30 +4793,22 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
computeParamsDerivatives(paramsDerivsOrder);
}
jacob_map_t contemporaneous_jacobian, static_jacobian;
map<tuple<int, int, int>, expr_t> first_order_endo_derivatives;
// For each simultaneous block contains pair<Size, Num_Feedback_variable>
vector<pair<int, int>> simblock_size;
vector<int> n_static, n_forward, n_backward, n_mixed;
if (linear_decomposition)
{
first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
is_equation_linear = equationLinear(first_order_endo_derivatives);
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
equationLinear(first_order_endo_derivatives);
tie(contemporaneous_jacobian, static_jacobian) = evaluateAndReduceJacobian(eval_context, cutoff, false);
auto [contemporaneous_jacobian, static_jacobian] = evaluateAndReduceJacobian(eval_context, cutoff, false);
if (!computeNaturalNormalization())
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
tie(simblock_size, n_static, n_forward, n_backward, n_mixed)
= select_non_linear_equations_and_variables(is_equation_linear);
select_non_linear_equations_and_variables();
equationTypeDetermination(first_order_endo_derivatives, 0);
prologue = 0;
epilogue = 0;
reduceBlocksAndTypeDetermination(simblock_size, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
reduceBlocksAndTypeDetermination(linear_decomposition);
computeChainRuleJacobian();
@ -4830,26 +4822,23 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
if (!no_tmp_terms)
computeTemporaryTermsOrdered();
}
if (block)
else if (block)
{
tie(contemporaneous_jacobian, static_jacobian) = evaluateAndReduceJacobian(eval_context, cutoff, false);
auto [contemporaneous_jacobian, static_jacobian] = evaluateAndReduceJacobian(eval_context, cutoff, false);
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
computePrologueAndEpilogue(static_jacobian);
first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
equationTypeDetermination(first_order_endo_derivatives, mfs);
cout << "Finding the optimal block decomposition of the model ..." << endl;
lag_lead_vector_t variable_lag_lead;
auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false);
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);
reduceBlocksAndTypeDetermination(linear_decomposition);
printBlockDecomposition();

View File

@ -162,19 +162,16 @@ ModelTree::ModelTree(const ModelTree &m) :
eq_idx_orig2block{m.eq_idx_orig2block},
endo_idx_orig2block{m.endo_idx_orig2block},
map_idx{m.map_idx},
block_type_firstequation_size_mfs{m.block_type_firstequation_size_mfs},
blocks_linear{m.blocks_linear},
block_col_type{m.block_col_type},
endo_max_leadlag_block{m.endo_max_leadlag_block},
other_endo_max_leadlag_block{m.other_endo_max_leadlag_block},
exo_max_leadlag_block{m.exo_max_leadlag_block},
exo_det_max_leadlag_block{m.exo_det_max_leadlag_block},
max_leadlag_block{m.max_leadlag_block},
blocks{m.blocks},
is_equation_linear{m.is_equation_linear},
endo2eq{m.endo2eq},
epilogue{m.epilogue},
prologue{m.prologue},
block_lag_lead{m.block_lag_lead},
cutoff{m.cutoff},
mfs{m.mfs}
{
@ -215,25 +212,21 @@ ModelTree::operator=(const ModelTree &m)
first_chain_rule_derivatives.clear();
map_idx = m.map_idx;
equation_type_and_normalized_equation.clear();
block_type_firstequation_size_mfs = m.block_type_firstequation_size_mfs;
blocks_derivatives.clear();
blocks_linear = m.blocks_linear;
derivative_endo.clear();
derivative_other_endo.clear();
derivative_exo.clear();
derivative_exo_det.clear();
block_col_type = m.block_col_type;
endo_max_leadlag_block = m.endo_max_leadlag_block;
other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
exo_max_leadlag_block = m.exo_max_leadlag_block;
exo_det_max_leadlag_block = m.exo_det_max_leadlag_block;
max_leadlag_block = m.max_leadlag_block;
blocks = m.blocks;
is_equation_linear = m.is_equation_linear;
endo2eq = m.endo2eq;
epilogue = m.epilogue;
prologue = m.prologue;
block_lag_lead = m.block_lag_lead;
cutoff = m.cutoff;
mfs = m.mfs;
@ -456,43 +449,46 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double
return { contemporaneous_jacobian, static_jacobian };
}
tuple<vector<pair<int, int>>, vector<int>, vector<int>, vector<int>, vector<int>>
ModelTree::select_non_linear_equations_and_variables(const vector<bool> &is_equation_linear)
void
ModelTree::select_non_linear_equations_and_variables()
{
vector<int> eq2endo(equations.size(), 0);
int num = 0;
for (auto it : endo2eq)
if (!is_equation_linear[it])
num++;
vector<int> endo2block(endo2eq.size(), 1);
int i = 0, j = 0;
for (auto it : endo2eq)
if (!is_equation_linear[it])
{
eq_idx_block2orig[i] = it;
endo_idx_block2orig[i] = j;
endo2block[j] = 0;
i++;
j++;
}
auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block);
vector<int> 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<int>(endo2eq.size()); i++)
prologue = 0;
epilogue = 0;
vector<int> endo2block(endo2eq.size(), 1); // The 1 is a dummy value, distinct from 0
int i = 0;
for (int endo = 0; endo < static_cast<int>(endo2eq.size()); endo++)
{
if (variable_lag_lead[endo_idx_block2orig[i]].first != 0 && variable_lag_lead[endo_idx_block2orig[i]].second != 0)
n_mixed[i]++;
else if (variable_lag_lead[endo_idx_block2orig[i]].first == 0 && variable_lag_lead[endo_idx_block2orig[i]].second != 0)
n_forward[i]++;
else if (variable_lag_lead[endo_idx_block2orig[i]].first != 0 && variable_lag_lead[endo_idx_block2orig[i]].second == 0)
n_backward[i]++;
else if (variable_lag_lead[endo_idx_block2orig[i]].first == 0 && variable_lag_lead[endo_idx_block2orig[i]].second == 0)
n_static[i]++;
int eq = endo2eq[endo];
if (!is_equation_linear[eq])
{
eq_idx_block2orig[i] = eq;
endo_idx_block2orig[i] = endo;
endo2block[i] = 0;
i++;
}
}
cout.flush();
vector<pair<int, int>> simblock_size(1, {i, i});
updateReverseVariableEquationOrderings();
return { simblock_size, n_static, n_forward, n_backward, n_mixed };
blocks.clear();
blocks.resize(1);
blocks[0].size = i;
blocks[0].mfs_size = i;
auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block);
for (int i = 0; i < blocks[0].size; i++)
{
auto [max_lag, max_lead] = variable_lag_lead[endo_idx_block2orig[i]];
if (max_lag != 0 && max_lead != 0)
blocks[0].n_mixed++;
else if (max_lag == 0 && max_lead != 0)
blocks[0].n_forward++;
else if (max_lag != 0 && max_lead == 0)
blocks[0].n_backward++;
else
blocks[0].n_static++;
}
}
bool
@ -701,9 +697,8 @@ ModelTree::getVariableLeadLagByBlock(const vector<int> &endo2simblock) const
return { equation_lag_lead, variable_lag_lead };
}
tuple<vector<pair<int, int>>, lag_lead_vector_t,
vector<int>, vector<int>, vector<int>, vector<int>>
ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_)
lag_lead_vector_t
ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_)
{
int nb_var = symbol_table.endo_nbr();
int nb_simvars = nb_var - prologue - epilogue;
@ -731,13 +726,28 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
index, starting from 0, in recursive order */
auto [num_simblocks, endo2simblock] = G.sortedStronglyConnectedComponents();
vector<pair<int, int>> simblock_size(num_simblocks, { 0, 0 });
int num_blocks = prologue+num_simblocks+epilogue;
// equations belonging to the block
blocks.clear();
blocks.resize(num_blocks);
// Initialize size and mfs_size for prologue and epilogue
for (int i = 0; i < prologue; i++)
{
blocks[i].size = 1;
blocks[i].mfs_size = 1;
}
for (int i = 0; i < epilogue; i++)
{
blocks[prologue+num_simblocks+i].size = 1;
blocks[prologue+num_simblocks+i].mfs_size = 1;
}
// Compute size and list of equations for simultaneous blocks
vector<set<int>> eqs_in_simblock(num_simblocks);
for (int i = 0; i < static_cast<int>(endo2simblock.size()); i++)
{
simblock_size[endo2simblock[i]].first++;
blocks[prologue+endo2simblock[i]].size++;
eqs_in_simblock[endo2simblock[i]].insert(i);
}
@ -746,7 +756,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
/* Add a loop on vertices which could not be normalized or vertices related
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
if (equation_type_and_normalized_equation[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
@ -754,23 +764,19 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
|| 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
vector<int> n_static(num_blocks, 0), n_forward(num_blocks, 0),
n_backward(num_blocks, 0), n_mixed(num_blocks, 0);
const vector<int> old_eq_idx_block2orig(eq_idx_block2orig), old_endo_idx_block2orig(endo_idx_block2orig);
for (int i = 0; i < prologue; i++)
{
auto [max_lag, max_lead] = variable_lag_lead[old_endo_idx_block2orig[i]];
if (max_lag != 0 && max_lead != 0)
n_mixed[i]++;
blocks[i].n_mixed++;
else if (max_lag == 0 && max_lead != 0)
n_forward[i]++;
blocks[i].n_forward++;
else if (max_lag != 0 && max_lead == 0)
n_backward[i]++;
blocks[i].n_backward++;
else
n_static[i]++;
blocks[i].n_static++;
}
/* For each block, the minimum set of feedback variable is computed and the
@ -783,7 +789,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
auto subG = G.extractSubgraph(eqs_in_simblock[i]);
auto feed_back_vertices = subG.minimalSetOfFeedbackVertices();
auto v_index1 = get(boost::vertex_index1, subG);
simblock_size[i].second = feed_back_vertices.size();
blocks[prologue+i].mfs_size = feed_back_vertices.size();
auto reordered_vertices = subG.reorderRecursiveVariables(feed_back_vertices);
// First we have the recursive equations conditional on feedback variables
@ -799,22 +805,22 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
};
if (j == 2 && max_lag != 0 && max_lead != 0)
{
n_mixed[prologue+i]++;
blocks[prologue+i].n_mixed++;
reorder();
}
else if (j == 3 && max_lag == 0 && max_lead != 0)
{
n_forward[prologue+i]++;
blocks[prologue+i].n_forward++;
reorder();
}
else if (j == 1 && max_lag != 0 && max_lead == 0)
{
n_backward[prologue+i]++;
blocks[prologue+i].n_backward++;
reorder();
}
else if (j == 0 && max_lag == 0 && max_lead == 0)
{
n_static[prologue+i]++;
blocks[prologue+i].n_static++;
reorder();
}
}
@ -833,22 +839,22 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
};
if (j == 2 && max_lag != 0 && max_lead != 0)
{
n_mixed[prologue+i]++;
blocks[prologue+i].n_mixed++;
reorder();
}
else if (j == 3 && max_lag == 0 && max_lead != 0)
{
n_forward[prologue+i]++;
blocks[prologue+i].n_forward++;
reorder();
}
else if (j == 1 && max_lag != 0 && max_lead == 0)
{
n_backward[prologue+i]++;
blocks[prologue+i].n_backward++;
reorder();
}
else if (j == 0 && max_lag == 0 && max_lead == 0)
{
n_static[prologue+i]++;
blocks[prologue+i].n_static++;
reorder();
}
}
@ -858,18 +864,18 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
{
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]++;
blocks[prologue+num_simblocks+i].n_mixed++;
else if (max_lag == 0 && max_lead != 0)
n_forward[prologue+num_simblocks+i]++;
blocks[prologue+num_simblocks+i].n_forward++;
else if (max_lag != 0 && max_lead == 0)
n_backward[prologue+num_simblocks+i]++;
blocks[prologue+num_simblocks+i].n_backward++;
else
n_static[prologue+num_simblocks+i]++;
blocks[prologue+num_simblocks+i].n_static++;
}
updateReverseVariableEquationOrderings();
return { simblock_size, variable_lag_lead, n_static, n_forward, n_backward, n_mixed };
return variable_lag_lead;
}
void
@ -904,185 +910,152 @@ ModelTree::printBlockDecomposition() const
}
void
ModelTree::reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &simblock_size, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &n_static, const vector<int> &n_forward, const vector<int> &n_backward, const vector<int> &n_mixed, bool linear_decomposition)
ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition)
{
int count_equ = 0, simblock_counter = 0;
block_type_firstequation_size_mfs.clear();
BlockSimulationType prev_Type = BlockSimulationType::unknown;
int eq = 0;
int num_simblocks = simblock_size.size();
for (int i = 0; i < prologue+num_simblocks+epilogue; i++)
for (int i = 0; i < static_cast<int>(blocks.size()); i++)
{
int Blck_Size, MFS_Size;
int first_count_equ = count_equ;
if (i < prologue)
{
Blck_Size = 1;
MFS_Size = 1;
}
else if (i < prologue+num_simblocks)
{
Blck_Size = simblock_size[simblock_counter].first;
MFS_Size = simblock_size[simblock_counter].second;
simblock_counter++;
}
else if (i < prologue+num_simblocks+epilogue)
{
Blck_Size = 1;
MFS_Size = 1;
}
int first_eq = (i == 0) ? 0 : blocks[i-1].first_equation+blocks[i-1].size;
int Lag = 0, Lead = 0;
for (count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
/* 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++)
{
set<pair<int, int>> endos_and_lags;
equations[eq_idx_block2orig[count_equ]]->collectEndogenous(endos_and_lags);
for (const auto &[curr_variable, curr_lag] : endos_and_lags)
equations[eq_idx_block2orig[eq]]->collectEndogenous(endos_and_lags);
for (const auto &[endo, lag] : endos_and_lags)
{
if (linear_decomposition)
{
if (dynamic_jacobian.find({ curr_lag, eq_idx_block2orig[count_equ], curr_variable }) != dynamic_jacobian.end())
if (dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo })
!= dynamic_jacobian.end())
{
if (curr_lag > Lead)
Lead = curr_lag;
else if (-curr_lag > Lag)
Lag = -curr_lag;
max_lead = max(lag, max_lead);
max_lag = max(-lag, max_lag);
}
}
else
{
if (find(endo_idx_block2orig.begin()+first_count_equ, endo_idx_block2orig.begin()+(first_count_equ+Blck_Size), curr_variable)
!= endo_idx_block2orig.begin()+(first_count_equ+Blck_Size)
&& dynamic_jacobian.find({ curr_lag, eq_idx_block2orig[count_equ], curr_variable }) != dynamic_jacobian.end())
if (find(endo_idx_block2orig.begin()+first_eq,
endo_idx_block2orig.begin()+first_eq+blocks[i].size, endo)
!= endo_idx_block2orig.begin()+first_eq+blocks[i].size
&& dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo })
!= dynamic_jacobian.end())
{
if (curr_lag > Lead)
Lead = curr_lag;
else if (-curr_lag > Lag)
Lag = -curr_lag;
max_lead = max(lag, max_lead);
max_lag = max(-lag, max_lag);
}
}
}
}
// Determine the block type
BlockSimulationType Simulation_Type;
if (Lag > 0 && Lead > 0)
if (max_lag > 0 && max_lead > 0)
{
if (Blck_Size == 1)
if (blocks[i].size == 1)
Simulation_Type = BlockSimulationType::solveTwoBoundariesSimple;
else
Simulation_Type = BlockSimulationType::solveTwoBoundariesComplete;
}
else if (Blck_Size > 1)
else if (blocks[i].size > 1)
{
if (Lead > 0)
if (max_lead > 0)
Simulation_Type = BlockSimulationType::solveBackwardComplete;
else
Simulation_Type = BlockSimulationType::solveForwardComplete;
}
else
{
if (Lead > 0)
if (max_lead > 0)
Simulation_Type = BlockSimulationType::solveBackwardSimple;
else
Simulation_Type = BlockSimulationType::solveForwardSimple;
}
int l_n_static = n_static[i];
int l_n_forward = n_forward[i];
int l_n_backward = n_backward[i];
int l_n_mixed = n_mixed[i];
if (Blck_Size == 1)
if (blocks[i].size == 1)
{
if (Equation_Type[eq_idx_block2orig[eq]].first == EquationType::evaluate
|| Equation_Type[eq_idx_block2orig[eq]].first == EquationType::evaluate_s)
// 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 (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 (i > 0)
{
bool is_lead = false, is_lag = false;
int c_Size = get<2>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
int first_equation = get<1>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
if (c_Size > 0
&& ((prev_Type == BlockSimulationType::evaluateForward && Simulation_Type == BlockSimulationType::evaluateForward && !is_lead)
|| (prev_Type == BlockSimulationType::evaluateBackward && Simulation_Type == BlockSimulationType::evaluateBackward && !is_lag)))
for (int j = blocks[i-1].first_equation;
j < blocks[i-1].first_equation+blocks[i-1].size; j++)
{
for (int j = first_equation; j < first_equation+c_Size; j++)
{
if (dynamic_jacobian.find({ -1, eq_idx_block2orig[eq], endo_idx_block2orig[j] })
!= dynamic_jacobian.end())
is_lag = true;
if (dynamic_jacobian.find({ +1, eq_idx_block2orig[eq], endo_idx_block2orig[j] })
!= dynamic_jacobian.end())
is_lead = true;
}
if (dynamic_jacobian.find({ -1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] })
!= dynamic_jacobian.end())
is_lag = true;
if (dynamic_jacobian.find({ +1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] })
!= dynamic_jacobian.end())
is_lead = true;
}
if ((prev_Type == BlockSimulationType::evaluateForward && Simulation_Type == BlockSimulationType::evaluateForward && !is_lead)
|| (prev_Type == BlockSimulationType::evaluateBackward && Simulation_Type == BlockSimulationType::evaluateBackward && !is_lag))
BlockSimulationType prev_Type = blocks[i-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 with the previous one
BlockSimulationType c_Type = get<0>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
c_Size++;
block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1] = { c_Type, first_equation, c_Size, c_Size };
if (block_lag_lead[block_type_firstequation_size_mfs.size()-1].first > Lag)
Lag = block_lag_lead[block_type_firstequation_size_mfs.size()-1].first;
if (block_lag_lead[block_type_firstequation_size_mfs.size()-1].second > Lead)
Lead = block_lag_lead[block_type_firstequation_size_mfs.size()-1].second;
block_lag_lead[block_type_firstequation_size_mfs.size()-1] = { Lag, Lead };
auto tmp = block_col_type[block_col_type.size()-1];
block_col_type[block_col_type.size()-1] = { get<0>(tmp)+l_n_static, get<1>(tmp)+l_n_forward, get<2>(tmp)+l_n_backward, get<3>(tmp)+l_n_mixed };
}
else
{
block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
// Merge the current block into the previous one
blocks[i-1].size++;
blocks[i-1].mfs_size = blocks[i-1].size;
blocks[i-1].max_lag = max(blocks[i-1].max_lag, max_lag);
blocks[i-1].max_lead = 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);
i--;
continue;
}
}
else
{
block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
}
else
{
block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
prev_Type = Simulation_Type;
eq += Blck_Size;
blocks[i].simulation_type = Simulation_Type;
blocks[i].first_equation = first_eq;
blocks[i].max_lag = max_lag;
blocks[i].max_lead = max_lead;
}
}
vector<bool>
ModelTree::equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives) const
void
ModelTree::equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives)
{
vector<bool> is_linear(symbol_table.endo_nbr(), true);
for (const auto &it : first_order_endo_derivatives)
is_equation_linear.clear();
is_equation_linear.resize(symbol_table.endo_nbr(), true);
for (const auto &[indices, expr] : first_order_endo_derivatives)
{
expr_t Id = it.second;
set<pair<int, int>> endogenous;
Id->collectEndogenous(endogenous);
expr->collectEndogenous(endogenous);
if (endogenous.size() > 0)
{
int eq = get<0>(it.first);
is_linear[eq] = false;
int eq = get<0>(indices);
is_equation_linear[eq] = false;
}
}
return is_linear;
}
void
ModelTree::determineLinearBlocks()
{
int nb_blocks = getNbBlocks();
blocks_linear.clear();
blocks_linear.resize(nb_blocks, true);
for (int block = 0; block < nb_blocks; block++)
// Note that field “linear” in class BlockInfo defaults to true
for (int block = 0; block < getNbBlocks(); block++)
{
BlockSimulationType simulation_type = getBlockSimulationType(block);
int block_size = getBlockSize(block);
@ -1100,7 +1073,7 @@ ModelTree::determineLinearBlocks()
for (int l = 0; l < block_size; l++)
if (endogenous.find({ endo_idx_block2orig[first_variable_position+l], 0 }) != endogenous.end())
{
blocks_linear[block] = false;
blocks[block].linear = false;
goto the_end;
}
}
@ -1115,7 +1088,7 @@ ModelTree::determineLinearBlocks()
for (int l = 0; l < block_size; l++)
if (endogenous.find({ endo_idx_block2orig[first_variable_position+l], lag }) != endogenous.end())
{
blocks_linear[block] = false;
blocks[block].linear = false;
goto the_end;
}
}

View File

@ -55,9 +55,6 @@ using equation_type_and_normalized_equation_t = vector<pair<EquationType, expr_t
//! Vector describing variables: max_lag in the block, max_lead in the block
using lag_lead_vector_t = vector<pair<int, int>>;
//! for each block contains tuple<Simulation_Type, first_equation, Block_Size, Recursive_part_Size>
using block_type_firstequation_size_mfs_t = vector<tuple<BlockSimulationType, int, int, int>>;
//! for a block contains derivatives tuple<block_equation_number, block_variable_number, lead_lag, expr_t>
using block_derivatives_equation_variable_laglead_nodeid_t = vector<tuple<int, int, int, expr_t>>;
@ -190,27 +187,32 @@ protected:
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
equation_type_and_normalized_equation_t equation_type_and_normalized_equation;
//! For each block contains tuple<Simulation_Type, first_equation, Block_Size, Recursive_part_Size>
block_type_firstequation_size_mfs_t block_type_firstequation_size_mfs;
//! for all blocks derivatives description
blocks_derivatives_t blocks_derivatives;
//! Vector indicating if the block is linear in endogenous variable (true) or not (false)
vector<bool> blocks_linear;
//! Map the derivatives for a block tuple<lag, eq, var>
using derivative_t = map<tuple<int, int, int>, expr_t>;
//! Vector of derivative for each blocks
vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
//! for each block described the number of static, forward, backward and mixed variables in the block
/*! tuple<static, forward, backward, mixed> */
vector<tuple<int, int, int, int>> block_col_type;
//!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
vector<pair<int, int>> endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
class BlockInfo
{
public:
BlockSimulationType simulation_type;
int first_equation; // Stores a block-ordered equation ID
int size{0};
int mfs_size{0}; // Size of the minimal feedback set
bool linear{true}; // Whether the block is linear in endogenous variable
int n_static{0}, n_forward{0}, n_backward{0}, n_mixed{0};
int max_lag{0}, max_lead{0};
};
// Stores various informations on the blocks
vector<BlockInfo> blocks;
//! the file containing the model and the derivatives code
ofstream code_file;
@ -270,9 +272,6 @@ protected:
//! number of equation in the prologue and in the epilogue
int epilogue, prologue;
//! for each block contains pair< max_lag, max_lead>
lag_lead_vector_t block_lag_lead;
/* 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;
@ -302,17 +301,22 @@ protected:
dynamic_jacobian. Elements below the cutoff are discarded. External functions are evaluated to 1. */
pair<jacob_map_t, jacob_map_t> 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, n_static, n_forward, n_backward, n_mixed) */
tuple<vector<pair<int, int>>, vector<int>, vector<int>, vector<int>, vector<int>> select_non_linear_equations_and_variables(const vector<bool> &is_equation_linear);
void select_non_linear_equations_and_variables();
//! 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<tuple<int, int, int>, 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, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) */
tuple<vector<pair<int, int>>, lag_lead_vector_t, vector<int>, vector<int>, vector<int>, vector<int>> 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<pair<int, int>> &simblock_size, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &n_static, const vector<int> &n_forward, const vector<int> &n_backward, const vector<int> &n_mixed, bool linear_decomposition);
/* 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,
mfs_size, n_static, n_forward, n_backward, n_mixed. */
lag_lead_vector_t computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_);
/* 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 following fields of the blocks structure: simulation_type,
first_equation, 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
belong to the same block (i.e. those endogenous are solved in the same
@ -323,7 +327,7 @@ protected:
which it belongs. */
pair<lag_lead_vector_t, lag_lead_vector_t> getVariableLeadLagByBlock(const vector<int> &endo2simblock) const;
//! For each equation determine if it is linear or not
vector<bool> equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives) const;
void equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives);
//! Print an abstract of the block structure of the model
void printBlockDecomposition() const;
//! Determine for each block if it is linear or not
@ -337,25 +341,25 @@ protected:
int
getNbBlocks() const
{
return block_type_firstequation_size_mfs.size();
return blocks.size();
};
//! Determine the simulation type of each block
BlockSimulationType
getBlockSimulationType(int block_number) const
{
return get<0>(block_type_firstequation_size_mfs[block_number]);
return blocks[block_number].simulation_type;
};
//! Return the first equation number of a block
int
getBlockFirstEquation(int block_number) const
{
return get<1>(block_type_firstequation_size_mfs[block_number]);
return blocks[block_number].first_equation;
};
//! Return the size of the block block_number
int
getBlockSize(int block_number) const
{
return get<2>(block_type_firstequation_size_mfs[block_number]);
return blocks[block_number].size;
};
//! Return the number of exogenous variable in the block block_number
virtual int getBlockExoSize(int block_number) const = 0;
@ -365,61 +369,62 @@ protected:
int
getBlockMfs(int block_number) const
{
return get<3>(block_type_firstequation_size_mfs[block_number]);
return blocks[block_number].mfs_size;
};
//! Return the maximum lag in a block
int
getBlockMaxLag(int block_number) const
{
return block_lag_lead[block_number].first;
return blocks[block_number].max_lag;
};
//! Return the maximum lead in a block
int
getBlockMaxLead(int block_number) const
{
return block_lag_lead[block_number].second;
return blocks[block_number].max_lead;
};
inline void
setBlockLeadLag(int block, int max_lag, int max_lead)
{
block_lag_lead[block] = { max_lag, max_lead };
blocks[block].max_lag = max_lag;
blocks[block].max_lead = max_lead;
};
//! Return the type of equation (equation_number) belonging to the block block_number
EquationType
getBlockEquationType(int block_number, int equation_number) const
{
return equation_type_and_normalized_equation[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first;
return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].first;
};
//! Return true if the equation has been normalized
bool
isBlockEquationRenormalized(int block_number, int equation_number) const
{
return equation_type_and_normalized_equation[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == EquationType::evaluate_s;
return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].first == EquationType::evaluate_s;
};
//! Return the expr_t of the equation equation_number belonging to the block block_number
expr_t
getBlockEquationExpr(int block_number, int equation_number) const
{
return equations[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]];
return equations[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]];
};
//! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
expr_t
getBlockEquationRenormalizedExpr(int block_number, int equation_number) const
{
return equation_type_and_normalized_equation[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].second;
return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].second;
};
//! Return the original number of equation equation_number belonging to the block block_number
int
getBlockEquationID(int block_number, int equation_number) const
{
return eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number];
return eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number];
};
//! Return the original number of variable variable_number belonging to the block block_number
int
getBlockVariableID(int block_number, int variable_number) const
{
return endo_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+variable_number];
return endo_idx_block2orig[getBlockFirstEquation(block_number)+variable_number];
};
//! Return the original number of the exogenous variable varexo_number belonging to the block block_number
virtual int getBlockVariableExoID(int block_number, int variable_number) const = 0;
@ -427,13 +432,13 @@ protected:
int
getBlockInitialEquationID(int block_number, int equation_number) const
{
return eq_idx_orig2block[equation_number] - get<1>(block_type_firstequation_size_mfs[block_number]);
return eq_idx_orig2block[equation_number] - getBlockFirstEquation(block_number);
};
//! Return the position of variable_number in the block number belonging to the block block_number
int
getBlockInitialVariableID(int block_number, int variable_number) const
{
return endo_idx_orig2block[variable_number] - get<1>(block_type_firstequation_size_mfs[block_number]);
return endo_idx_orig2block[variable_number] - getBlockFirstEquation(block_number);
};
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int getBlockInitialExogenousID(int block_number, int variable_number) const = 0;

View File

@ -703,7 +703,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
block_size,
endo_idx_block2orig,
eq_idx_block2orig,
blocks_linear[block],
blocks[block].linear,
symbol_table.endo_nbr(),
0,
0,
@ -1142,9 +1142,9 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
cout << "Finding the optimal block decomposition of the model ..." << endl;
auto [blocks, variable_lag_lead, n_static, n_forward, n_backward, n_mixed] = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false);
auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false);
reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, false);
reduceBlocksAndTypeDetermination(false);
printBlockDecomposition();