Commit bd6eee93 authored by Sébastien Villemot's avatar Sébastien Villemot

Block decomposition: refactor the prototypes of various functions

— return output arguments on the left-hand side
— do not pass class members as input/output arguments

By the way, fix a (benign) vector allocation bug in
{Static,Dynamic}Model::computeChainRuleJacobian().
parent 76c2e87c
Pipeline #3388 passed with stages
in 4 minutes and 6 seconds
......@@ -4822,27 +4822,25 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
is_equation_linear = equationLinear(first_order_endo_derivatives);
evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
tie(contemporaneous_jacobian, static_jacobian) = evaluateAndReduceJacobian(eval_context, cutoff, false);
if (!computeNaturalNormalization())
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
lag_lead_vector_t equation_lag_lead, variable_lag_lead;
blocks = select_non_linear_equations_and_variables(is_equation_linear, dynamic_jacobian, equation_reordered, variable_reordered,
inv_equation_reordered, inv_variable_reordered,
equation_lag_lead, variable_lag_lead,
n_static, n_forward, n_backward, n_mixed);
tie(blocks, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed)
= select_non_linear_equations_and_variables(is_equation_linear);
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, 0);
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, 0);
prologue = 0;
epilogue = 0;
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, linear_decomposition);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
computeChainRuleJacobian(blocks_derivatives);
computeChainRuleJacobian();
blocks_linear = BlockLinear(blocks_derivatives, variable_reordered);
blocks_linear = BlockLinear();
collect_block_first_order_derivatives();
......@@ -4855,29 +4853,29 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
if (block)
{
evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
tie(contemporaneous_jacobian, static_jacobian) = evaluateAndReduceJacobian(eval_context, cutoff, false);
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
computePrologueAndEpilogue(static_jacobian);
first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, mfs);
cout << "Finding the optimal block decomposition of the model ..." << endl;
lag_lead_vector_t equation_lag_lead, variable_lag_lead;
computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
tie(blocks, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false, true);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, linear_decomposition);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
printBlockDecomposition(blocks);
computeChainRuleJacobian(blocks_derivatives);
computeChainRuleJacobian();
blocks_linear = BlockLinear(blocks_derivatives, variable_reordered);
blocks_linear = BlockLinear();
collect_block_first_order_derivatives();
......@@ -4888,7 +4886,8 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
computeTemporaryTermsOrdered();
int k = 0;
equation_block.resize(equations.size());
variable_block_lead_lag = vector<tuple<int, int, int>>(equations.size());
variable_block_lead_lag.clear();
variable_block_lead_lag.resize(equations.size());
for (unsigned int i = 0; i < getNbBlocks(); i++)
{
for (unsigned int j = 0; j < getBlockSize(i); j++)
......@@ -5072,11 +5071,11 @@ DynamicModel::get_Derivatives(int block)
}
void
DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_endo_derivatives)
DynamicModel::computeChainRuleJacobian()
{
map<int, expr_t> recursive_variables;
unsigned int nb_blocks = getNbBlocks();
blocks_endo_derivatives = blocks_derivatives_t(nb_blocks);
blocks_derivatives.resize(nb_blocks);
for (unsigned int block = 0; block < nb_blocks; block++)
{
block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
......@@ -5084,7 +5083,6 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_endo_derivat
int block_size = getBlockSize(block);
int block_nb_mfs = getBlockMfs(block);
int block_nb_recursives = block_size - block_nb_mfs;
blocks_endo_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
for (int i = 0; i < block_nb_recursives; i++)
{
if (getBlockEquationType(block, i) == E_EVALUATE_S)
......@@ -5110,7 +5108,7 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_endo_derivat
}
tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]);
}
blocks_endo_derivatives[block] = tmp_derivatives;
blocks_derivatives[block] = tmp_derivatives;
}
}
......
......@@ -138,7 +138,7 @@ private:
//! return a map on the block jacobian
map<tuple<int, int, int, int, int>, int> get_Derivatives(int block);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
void computeChainRuleJacobian();
string reform(const string &name) const;
......
This diff is collapsed.
......@@ -53,7 +53,7 @@ vectorToTuple(const vector<T> &v)
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>>;
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>>;
......@@ -279,34 +279,38 @@ protected:
If no matching is found using a strictly positive cutoff, then a zero cutoff is applied (i.e. use a symbolic normalization); in that case, the method adds zeros in the jacobian matrices to reflect all the edges in the symbolic incidence matrix.
If no matching is found with a zero cutoff close to zero an error message is printout.
*/
void computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian);
void computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian);
//! Try to find a natural normalization if all equations are matched to an endogenous variable on the LHS
bool computeNaturalNormalization();
//! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
multimap<int, int> computeNormalizedEquations() const;
//! Evaluate the jacobian and suppress all the elements below the cutoff
void evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose);
/*! Returns a pair (contemporaneous_jacobian, static_jacobian). Also fills
dynamic_jacobian. */
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
vector<pair<int, int>> select_non_linear_equations_and_variables(vector<bool> is_equation_linear, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered,
vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered,
lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead,
vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed);
/*! Returns a tuple (blocks, equation_lag_lead, variable_lag_lead, n_static,
n_forward, n_backward, n_mixed) */
tuple<vector<pair<int, int>>, lag_lead_vector_t, lag_lead_vector_t, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>> select_non_linear_equations_and_variables(const vector<bool> &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, vector<int> &equation_reordered, vector<int> &variable_reordered);
void computePrologueAndEpilogue(const jacob_map_t &static_jacobian);
//! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
equation_type_and_normalized_equation_t equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
equation_type_and_normalized_equation_t equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs) const;
//! Compute the block decomposition and for a non-recusive block find the minimum feedback set
void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead_t, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const;
/*! Returns a tuple (blocks, equation_lag_lead, variable_lag_lead, n_static,
n_forward, n_backward, n_mixed) */
tuple<vector<pair<int, int>>, lag_lead_vector_t, lag_lead_vector_t, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>> computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable);
//! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<tuple<int, int, int, int>> &block_col_type, bool linear_decomposition);
block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<unsigned int> &n_static, const vector<unsigned int> &n_forward, const vector<unsigned int> &n_backward, const vector<unsigned int> &n_mixed, bool linear_decomposition);
//! Determine the maximum number of lead and lag for the endogenous variable in a bloc
void getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const;
/*! Returns a pair { equation_lead,lag, variable_lead_lag } */
pair<lag_lead_vector_t, lag_lead_vector_t> getVariableLeadLagByBlock(const vector<int> &components_set, int nb_blck_sim) const;
//! For each equation determine if it is linear or not
vector<bool> equationLinear(map<tuple<int, int, int>, expr_t> first_order_endo_derivatives) const;
vector<bool> equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives) const;
//! Print an abstract of the block structure of the model
void printBlockDecomposition(const vector<pair<int, int>> &blocks) const;
//! Determine for each block if it is linear or not
vector<bool> BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector<int> &variable_reordered) const;
vector<bool> BlockLinear() const;
//! Remove equations specified by exclude_eqs
vector<int> includeExcludeEquations(set<pair<string, string>> &eqs, bool exclude_eqs,
vector<BinaryOpNode *> &equations, vector<int> &equations_lineno,
......
......@@ -1144,35 +1144,27 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
if (block)
{
jacob_map_t contemporaneous_jacobian, static_jacobian;
vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
auto [contemporaneous_jacobian, static_jacobian] = evaluateAndReduceJacobian(eval_context, cutoff, false);
// for each block contains pair<Size, Feddback_variable>
vector<pair<int, int>> blocks;
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
computePrologueAndEpilogue(static_jacobian);
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
auto first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
map<tuple<int, int, int>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, mfs);
cout << "Finding the optimal block decomposition of the model ..." << endl;
lag_lead_vector_t equation_lag_lead, variable_lag_lead;
computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
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);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, false);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, false);
printBlockDecomposition(blocks);
computeChainRuleJacobian(blocks_derivatives);
computeChainRuleJacobian();
blocks_linear = BlockLinear(blocks_derivatives, variable_reordered);
blocks_linear = BlockLinear();
collect_block_first_order_derivatives();
......@@ -2176,11 +2168,11 @@ StaticModel::get_Derivatives(int block)
}
void
StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
StaticModel::computeChainRuleJacobian()
{
map<int, expr_t> recursive_variables;
unsigned int nb_blocks = getNbBlocks();
blocks_derivatives = blocks_derivatives_t(nb_blocks);
blocks_derivatives.resize(nb_blocks);
for (unsigned int block = 0; block < nb_blocks; block++)
{
block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
......@@ -2191,7 +2183,6 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
int block_nb_recursives = block_size - block_nb_mfs;
if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
{
blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
for (int i = 0; i < block_nb_recursives; i++)
{
if (getBlockEquationType(block, i) == E_EVALUATE_S)
......@@ -2220,7 +2211,6 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
}
else
{
blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
for (int i = 0; i < block_nb_recursives; i++)
{
if (getBlockEquationType(block, i) == E_EVALUATE_S)
......
......@@ -90,7 +90,7 @@ private:
//! return a map on the block jacobian
map<tuple<int, int, int, int, int>, int> get_Derivatives(int block);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
void computeChainRuleJacobian();
//! Collect only the first derivatives
map<tuple<int, int, int>, expr_t> collect_first_order_derivatives_endogenous();
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment