stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Block decomposition: simplify computation of block derivatives 

— use a std::map for storing block derivatives
— remove redundant ModelTree::first_chain_rule_derivatives structure
— remove unused codepaths in StaticModel
— DynamicModel: simplify code that determines the type of derivatives in a
  block. We now use a slightly different categorization.
— by the way, fix the max lead/lag information for blocks that are obtained via
  merging. A workaround was previously implemented in
  DynamicModel::get_Derivative(), but it is no longer needed with this fix.
issue#70
Sébastien Villemot 2020-04-24 12:29:02 +02:00
parent 5bd41f8599
commit 2e136cae56
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 162 additions and 282 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

203
src/DynamicModel.cc
View File

@ -165,9 +165,10 @@ DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_n
}
void
DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const
DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const map_idx_t &map_idx) const
{
if (auto it = first_chain_rule_derivatives.find({ eqr, varr, lag }); it != first_chain_rule_derivatives.end())
if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
it != blocks_derivatives[blk].end())
it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
else
{
@ -182,7 +183,6 @@ DynamicModel::computeTemporaryTermsOrdered()
map<expr_t, pair<int, int>> first_occurence;
map<expr_t, int> reference_count;
BinaryOpNode *eq_node;
first_chain_rule_derivatives_t::const_iterator it_chr;
ostringstream tmp_s;
v_temporary_terms.clear();
map_idx.clear();
@ -211,11 +211,8 @@ DynamicModel::computeTemporaryTermsOrdered()
eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
}
}
for (const auto &it : blocks_derivatives[block])
{
expr_t id = get<3>(it);
id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
}
for (const auto &[ignore, id] : blocks_derivatives[block])
id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
for (const auto &it : derivative_endo[block])
it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
for (const auto &it : derivative_other_endo[block])
@ -241,11 +238,8 @@ DynamicModel::computeTemporaryTermsOrdered()
eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
}
}
for (const auto &it : blocks_derivatives[block])
{
expr_t id = get<3>(it);
id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
}
for (const auto &[ignore, id] : blocks_derivatives[block])
id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
for (const auto &it : derivative_endo[block])
it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
for (const auto &it : derivative_other_endo[block])
@ -267,11 +261,8 @@ DynamicModel::computeTemporaryTermsOrdered()
eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
}
}
for (const auto &it : blocks_derivatives[block])
{
expr_t id = get<3>(it);
id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
}
for (const auto &[ignore, id] : blocks_derivatives[block])
id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
for (const auto &it : derivative_endo[block])
it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
for (const auto &it : derivative_other_endo[block])
@ -339,8 +330,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
int prev_lag;
int prev_var, count_col, count_col_endo, count_col_exo, count_col_exo_det, count_col_other_endo;
map<tuple<int, int, int>, expr_t> tmp_block_endo_derivative;
for (const auto &it : blocks_derivatives[block])
tmp_block_endo_derivative[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it);
for (const auto &[idx, d] : blocks_derivatives[block])
tmp_block_endo_derivative[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d;
prev_var = 999999999;
prev_lag = -9999999;
count_col_endo = 0;
@ -762,8 +753,9 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
output << " else" << endl;
for (const auto &[eq, var, lag, id] : blocks_derivatives[block])
for (const auto &[indices, id] : blocks_derivatives[block])
{
auto [eq, var, lag] = indices;
int eqr = getBlockEquationID(block, eq);
int varr = getBlockVariableID(block, var);
if (lag == 0)
@ -782,8 +774,9 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
case BlockSimulationType::solveTwoBoundariesSimple:
case BlockSimulationType::solveTwoBoundariesComplete:
output << " else" << endl;
for (const auto &[eq, var, lag, id] : blocks_derivatives[block])
for (const auto &[indices, id] : blocks_derivatives[block])
{
auto [eq, var, lag] = indices;
int eqr = getBlockEquationID(block, eq);
int varr = getBlockVariableID(block, var);
ostringstream tmp_output;
@ -1185,8 +1178,8 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
file_open = true;
}
map<tuple<int, int, int>, expr_t> tmp_block_endo_derivative;
for (const auto &it : blocks_derivatives[block])
tmp_block_endo_derivative[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it);
for (const auto &[idx, d] : blocks_derivatives[block])
tmp_block_endo_derivative[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d;
map<tuple<int, int, int>, expr_t> tmp_exo_derivative;
for (const auto &it : derivative_exo[block])
tmp_exo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
@ -1393,8 +1386,9 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
case BlockSimulationType::solveTwoBoundariesComplete:
case BlockSimulationType::solveTwoBoundariesSimple:
count_u = feedback_variables.size();
for (const auto &[eq, var, lag, ignore] : blocks_derivatives[block])
for (const auto &[indices, ignore] : blocks_derivatives[block])
{
auto [eq, var, lag] = indices;
int eqr = getBlockEquationID(block, eq);
int varr = getBlockVariableID(block, var);
if (eq >= block_recursive and var >= block_recursive)
@ -1419,7 +1413,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
Uf[eqr].Ufl->lag = lag;
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag);
fnumexpr.write(code_file, instruction_number);
compileChainRuleDerivative(code_file, instruction_number, eqr, varr, lag, map_idx);
compileChainRuleDerivative(code_file, instruction_number, block, eq, var, lag, map_idx);
FSTPU_ fstpu(count_u);
fstpu.write(code_file, instruction_number);
count_u++;
@ -1808,8 +1802,9 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, int num,
int block_size = blocks[num].size;
int block_mfs = blocks[num].mfs_size;
int block_recursive = blocks[num].getRecursiveSize();
for (const auto &[eq, var, lag, ignore] : blocks_derivatives[num])
for (const auto &[indices, ignore] : blocks_derivatives[num])
{
auto [eq, var, lag] = indices;
if (lag != 0 && !is_two_boundaries)
continue;
if (eq >= block_recursive && var >= block_recursive)
@ -3240,8 +3235,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
tmp_s.str("");
count_lead_lag_incidence = 0;
dynamic_jacob_map_t reordered_dynamic_jacobian;
for (const auto &it : blocks_derivatives[block])
reordered_dynamic_jacobian[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it);
for (const auto &[idx, d] : blocks_derivatives[block])
reordered_dynamic_jacobian[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d;
output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];" << endl;
int last_var = -1;
vector<int> local_state_var;
@ -3415,10 +3410,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
if (block == 0)
{
set<pair<int, int>> row_state_var_incidence;
for (const auto &it : blocks_derivatives[block])
if (auto it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(block, get<1>(it))+1);
for (const auto &[idx, ignore] : blocks_derivatives[block])
if (auto it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(block, get<1>(idx))+1);
it_state_var != state_var.end())
if (auto it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(block, get<0>(it))+1);
if (auto it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(block, get<0>(idx))+1);
it_state_equ != state_equ.end())
row_state_var_incidence.emplace(it_state_equ - state_equ.begin(), it_state_var - state_var.begin());
auto row_state_var_incidence_it = row_state_var_incidence.begin();
@ -4973,108 +4968,80 @@ DynamicModel::writeRevXrefs(ostream &output, const map<pair<int, int>, set<int>>
}
}
map<tuple<int, int, int, int, int>, int>
DynamicModel::get_Derivatives(int block)
{
int max_lag, max_lead;
map<tuple<int, int, int, int, int>, int> Derivatives;
BlockSimulationType simulation_type = blocks[block].simulation_type;
if (simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
{
max_lag = 1;
max_lead = 1;
blocks[block].max_lag = max_lag;
blocks[block].max_lead = max_lead;
}
else
{
max_lag = blocks[block].max_lag;
max_lead = blocks[block].max_lead;
}
int block_size = blocks[block].size;
int block_nb_recursive = blocks[block].getRecursiveSize();
for (int lag = -max_lag; lag <= max_lead; lag++)
{
for (int eq = 0; eq < block_size; eq++)
{
int eqr = getBlockEquationID(block, eq);
for (int var = 0; var < block_size; var++)
{
int varr = getBlockVariableID(block, var);
if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end())
{
bool OK = true;
if (auto its = Derivatives.find({ lag, eq, var, eqr, varr });
its != Derivatives.end() && its->second == 2)
OK = false;
if (OK)
{
if (getBlockEquationType(block, eq) == EquationType::evaluate_s
&& eq < block_nb_recursive)
//It's a normalized equation, we have to recompute the derivative using chain rule derivative function
Derivatives[{ lag, eq, var, eqr, varr }] = 1;
else
//It's a feedback equation we can use the derivatives
Derivatives[{ lag, eq, var, eqr, varr }] = 0;
}
if (var < block_nb_recursive)
{
int eqs = getBlockEquationID(block, var);
for (int vars = block_nb_recursive; vars < block_size; vars++)
{
int varrs = getBlockVariableID(block, vars);
//A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation)
if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end())
Derivatives[{ lag, eq, vars, eqr, varrs }] = 2;
}
}
}
}
}
}
return Derivatives;
map<tuple<int, int, int>, DynamicModel::BlockDerivativeType>
DynamicModel::determineBlockDerivativesType(int blk)
{
map<tuple<int, int, int>, BlockDerivativeType> derivType;
int size = blocks[blk].size;
int nb_recursive = blocks[blk].getRecursiveSize();
for (int lag = -blocks[blk].max_lag; lag <= blocks[blk].max_lead; lag++)
for (int eq = 0; eq < size; eq++)
for (int var = 0; var < size; var++)
if (int eq_orig = getBlockEquationID(blk, eq), var_orig = getBlockVariableID(blk, var);
dynamic_jacobian.find({ lag, eq_orig, var_orig }) != dynamic_jacobian.end())
{
if (getBlockEquationType(blk, eq) == EquationType::evaluate_s
&& eq < nb_recursive)
/* It’s a normalized recursive equation, we have to recompute
the derivative using the chain rule */
derivType[{ lag, eq, var }] = BlockDerivativeType::normalizedChainRule;
else if (derivType.find({ lag, eq, var }) == derivType.end())
derivType[{ lag, eq, var }] = BlockDerivativeType::standard;
if (var < nb_recursive)
for (int feedback_var = nb_recursive; feedback_var < size; feedback_var++)
if (derivType.find({ lag, var, feedback_var }) != derivType.end())
/* A new derivative needs to be computed using the chain rule
(a feedback variable appears in the recursive equation
defining the current variable) */
derivType[{ lag, eq, feedback_var }] = BlockDerivativeType::chainRule;
}
return derivType;
}
void
DynamicModel::computeChainRuleJacobian()
{
map<int, expr_t> recursive_variables;
int nb_blocks = blocks.size();
blocks_derivatives.resize(nb_blocks);
for (int block = 0; block < nb_blocks; block++)
for (int blk = 0; blk < nb_blocks; blk++)
{
block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
recursive_variables.clear();
int block_nb_recursives = blocks[block].getRecursiveSize();
for (int i = 0; i < block_nb_recursives; i++)
int nb_recursives = blocks[blk].getRecursiveSize();
// Create a map from recursive vars to their defining (normalized) equation
map<int, expr_t> recursive_vars;
for (int i = 0; i < nb_recursives; i++)
{
if (getBlockEquationType(block, i) == EquationType::evaluate_s)
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(blk, i)), 0);
if (getBlockEquationType(blk, i) == EquationType::evaluate_s)
recursive_vars[deriv_id] = getBlockEquationRenormalizedExpr(blk, i);
else
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
recursive_vars[deriv_id] = getBlockEquationExpr(blk, i);
}
auto Derivatives = get_Derivatives(block);
for (const auto &it : Derivatives)
// Compute the block derivatives
for (const auto &[indices, derivType] : determineBlockDerivativesType(blk))
{
int Deriv_type = it.second;
auto [lag, eq, var, eqr, varr] = it.first;
if (Deriv_type == 0)
first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
else if (Deriv_type == 1)
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
auto [lag, eq, var] = indices;
int eq_orig = getBlockEquationID(blk, eq), var_orig = getBlockVariableID(blk, var);
int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, var_orig), lag);
expr_t d{nullptr};
switch (derivType)
{
if (getBlockEquationType(block, eq) == EquationType::evaluate_s
&& eq < block_nb_recursives)
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else
first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
case BlockDerivativeType::standard:
d = derivatives[1][{ eq_orig, deriv_id }];
break;
case BlockDerivativeType::chainRule:
d = equations[eq_orig]->getChainRuleDerivative(deriv_id, recursive_vars);
break;
case BlockDerivativeType::normalizedChainRule:
d = equation_type_and_normalized_equation[eq_orig].second->getChainRuleDerivative(deriv_id, recursive_vars);
break;
}
tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]);
blocks_derivatives[blk][{ eq, var, lag }] = d;
}
blocks_derivatives[block] = tmp_derivatives;
}
}

22
src/DynamicModel.hh
View File

@ -128,15 +128,19 @@ private:
void writeSetAuxiliaryVariables(const string &basename, bool julia) const;
void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
//! Computes jacobian and prepares for equation normalization
/*! Using values from initval/endval blocks and parameter initializations:
- computes the jacobian for the model w.r. to contemporaneous variables
- removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
*/
//void evaluateJacobian(const eval_context_t &eval_context, jacob_map *j_m, bool dynamic);
// Used by determineBlockDerivativesType()
enum class BlockDerivativeType
{
standard,
chainRule,
normalizedChainRule
};
/* For each tuple (lag, eq, var) within the given block, determine the type
of the derivative to be computed. Indices are within the block (i.e.
between 0 and blocks[blk].size-1). */
map<tuple<int, int, int>, BlockDerivativeType> determineBlockDerivativesType(int blk);
//! 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();
@ -150,7 +154,7 @@ private:
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, const map_idx_t &map_idx) const;
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const map_idx_t &map_idx) const;
//! Get the type corresponding to a derivation ID
SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;

23
src/ModelTree.cc
View File

@ -103,17 +103,15 @@ ModelTree::copyHelper(const ModelTree &m)
nonstationary_symbols_map[it.first] = {it.second.first, f(it.second.second)};
for (const auto &it : m.dynamic_jacobian)
dynamic_jacobian[it.first] = f(it.second);
for (const auto &it : m.first_chain_rule_derivatives)
first_chain_rule_derivatives[it.first] = f(it.second);
for (const auto &it : m.equation_type_and_normalized_equation)
equation_type_and_normalized_equation.emplace_back(it.first, f(it.second));
for (const auto &it : m.blocks_derivatives)
{
block_derivatives_equation_variable_laglead_nodeid_t v;
map<tuple<int, int, int>, expr_t> v;
for (const auto &it2 : it)
v.emplace_back(get<0>(it2), get<1>(it2), get<2>(it2), f(get<3>(it2)));
v[it2.first] = f(it2.second);
blocks_derivatives.push_back(v);
}
@ -209,7 +207,6 @@ ModelTree::operator=(const ModelTree &m)
endo_idx_block2orig = m.endo_idx_block2orig;
eq_idx_orig2block = m.eq_idx_orig2block;
endo_idx_orig2block = m.endo_idx_orig2block;
first_chain_rule_derivatives.clear();
map_idx = m.map_idx;
equation_type_and_normalized_equation.clear();
blocks_derivatives.clear();
@ -965,8 +962,12 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition)
// 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);
/* 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;
@ -1010,12 +1011,13 @@ ModelTree::determineLinearBlocks()
{
BlockSimulationType simulation_type = blocks[block].simulation_type;
int block_size = blocks[block].size;
block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives[block];
auto derivatives_block = blocks_derivatives[block];
int first_variable_position = blocks[block].first_equation;
if (simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete)
for (const auto &[ignore, ignore2, lag, d1] : derivatives_block)
for (const auto &[indices, d1] : derivatives_block)
{
int lag = get<2>(indices);
if (lag == 0)
{
set<pair<int, int>> endogenous;
@ -1031,8 +1033,9 @@ ModelTree::determineLinearBlocks()
}
else if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
for (const auto &[ignore, ignore2, lag, d1] : derivatives_block)
for (const auto &[indices, d1] : derivatives_block)
{
int lag = get<2>(indices);
set<pair<int, int>> endogenous;
d1->collectEndogenous(endogenous);
if (endogenous.size() > 0)

16
src/ModelTree.hh
View File

@ -55,12 +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 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>>;
//! for all blocks derivatives description
using blocks_derivatives_t = vector<block_derivatives_equation_variable_laglead_nodeid_t>;
//! Shared code for static and dynamic models
class ModelTree : public DataTree
{
@ -178,17 +172,15 @@ protected:
Set by updateReverseVariableEquationOrderings() */
vector<int> eq_idx_orig2block, endo_idx_orig2block;
//! Store the derivatives or the chainrule derivatives:map<tuple<equation, variable, lead_lag>, expr_t>
using first_chain_rule_derivatives_t = map<tuple<int, int, int>, expr_t>;
first_chain_rule_derivatives_t first_chain_rule_derivatives;
map_idx_t map_idx;
//! 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 all blocks derivatives description
blocks_derivatives_t blocks_derivatives;
/* Stores derivatives of each block.
The tuple is: equation number (inside the block), variable number (inside
the block), lead/lag */
vector<map<tuple<int, int, int>, expr_t>> blocks_derivatives;
//! Map the derivatives for a block tuple<lag, eq, var>
using derivative_t = map<tuple<int, int, int>, expr_t>;

176
src/StaticModel.cc
View File

@ -139,10 +139,10 @@ StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_nu
}
void
StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
{
if (auto it = first_chain_rule_derivatives.find({ eqr, varr, lag });
it != first_chain_rule_derivatives.end())
if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
it != blocks_derivatives[blk].end())
it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
else
{
@ -157,7 +157,6 @@ StaticModel::computeTemporaryTermsOrdered()
map<expr_t, pair<int, int>> first_occurence;
map<expr_t, int> reference_count;
BinaryOpNode *eq_node;
first_chain_rule_derivatives_t::const_iterator it_chr;
ostringstream tmp_s;
map_idx.clear();
@ -192,11 +191,8 @@ StaticModel::computeTemporaryTermsOrdered()
eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, i);
}
}
for (const auto &it : blocks_derivatives[block])
{
expr_t id = get<3>(it);
id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, block_size-1);
}
for (const auto &[ignore, id] : blocks_derivatives[block])
id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, block_size-1);
v_temporary_terms_inuse[block] = {};
computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
}
@ -218,11 +214,8 @@ StaticModel::computeTemporaryTermsOrdered()
eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
}
}
for (const auto &it : blocks_derivatives[block])
{
expr_t id = get<3>(it);
id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
}
for (const auto &[ignore, id] : blocks_derivatives[block])
id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
}
for (int block = 0; block < nb_blocks; block++)
@ -241,11 +234,8 @@ StaticModel::computeTemporaryTermsOrdered()
eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
}
}
for (const auto &it : blocks_derivatives[block])
{
expr_t id = get<3>(it);
id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
}
for (const auto &[ignore, id] : blocks_derivatives[block])
id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
for (int i = 0; i < blocks[block].size; i++)
for (const auto &it : v_temporary_terms[block][i])
it->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
@ -453,8 +443,9 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
case BlockSimulationType::solveForwardSimple:
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
for (const auto &[eq, var, ignore, id] : blocks_derivatives[block])
for (const auto &[indices, id] : blocks_derivatives[block])
{
auto &[eq, var, ignore] = indices;
int eqr = getBlockEquationID(block, eq);
int varr = getBlockVariableID(block, var);
output << " g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
@ -817,8 +808,9 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
count_u = feedback_variables.size();
for (const auto &[eq, var, ignore, ignore2] : blocks_derivatives[block])
for (const auto &[indices, ignore2] : blocks_derivatives[block])
{
auto [eq, var, ignore] = indices;
int eqr = getBlockEquationID(block, eq);
int varr = getBlockVariableID(block, var);
if (eq >= block_recursive && var >= block_recursive)
@ -838,7 +830,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
Uf[eqr].Ufl->var = varr;
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
fnumexpr.write(code_file, instruction_number);
compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx, temporary_terms);
compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, map_idx, temporary_terms);
FSTPSU_ fstpsu(count_u);
fstpsu.write(code_file, instruction_number);
count_u++;
@ -1005,14 +997,15 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
count_u = feedback_variables.size();
for (const auto &[eq, var, ignore, ignore2] : blocks_derivatives[block])
for (const auto &[indices, ignore2] : blocks_derivatives[block])
{
auto &[eq, var, ignore] = indices;
int eqr = getBlockEquationID(block, eq);
int varr = getBlockVariableID(block, var);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
fnumexpr.write(code_file, instruction_number);
compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx2[block], tt2 /*temporary_terms*/);
compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, map_idx2[block], tt2 /*temporary_terms*/);
FSTPG2_ fstpg2(eq, var);
fstpg2.write(code_file, instruction_number);
@ -1055,8 +1048,9 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, int num,
int block_size = blocks[num].size;
int block_mfs = blocks[num].mfs_size;
int block_recursive = blocks[num].getRecursiveSize();
for (const auto &[eq, var, ignore, ignore2] : blocks_derivatives[num])
for (const auto &[indices, ignore2] : blocks_derivatives[num])
{
auto [eq, var, ignore] = indices;
int lag = 0;
if (eq >= block_recursive && var >= block_recursive)
{
@ -2104,121 +2098,43 @@ StaticModel::addAllParamDerivId(set<int> &deriv_id_set)
deriv_id_set.insert(i + symbol_table.endo_nbr());
}
map<tuple<int, int, int, int, int>, int>
StaticModel::get_Derivatives(int block)
{
map<tuple<int, int, int, int, int>, int> Derivatives;
int block_size = blocks[block].size;
int block_nb_recursive = blocks[block].getRecursiveSize();
int lag = 0;
for (int eq = 0; eq < block_size; eq++)
{
int eqr = getBlockEquationID(block, eq);
for (int var = 0; var < block_size; var++)
{
int varr = getBlockVariableID(block, var);
if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end())
{
bool OK = true;
if (auto its = Derivatives.find({ lag, eq, var, eqr, varr });
its != Derivatives.end() && its->second == 2)
OK = false;
if (OK)
{
if (getBlockEquationType(block, eq) == EquationType::evaluate_s
&& eq < block_nb_recursive)
//It's a normalized equation, we have to recompute the derivative using chain rule derivative function
Derivatives[{ lag, eq, var, eqr, varr }] = 1;
else
//It's a feedback equation we can use the derivatives
Derivatives[{ lag, eq, var, eqr, varr }] = 0;
}
if (var < block_nb_recursive)
{
int eqs = getBlockEquationID(block, var);
for (int vars = block_nb_recursive; vars < block_size; vars++)
{
int varrs = getBlockVariableID(block, vars);
//A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation)
if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end())
Derivatives[{ lag, eq, vars, eqr, varrs }] = 2;
}
}
}
}
}
return Derivatives;
}
void
StaticModel::computeChainRuleJacobian()
{
map<int, expr_t> recursive_variables;
int nb_blocks = blocks.size();
blocks_derivatives.resize(nb_blocks);
for (int block = 0; block < nb_blocks; block++)
for (int blk = 0; blk < nb_blocks; blk++)
{
block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
recursive_variables.clear();
BlockSimulationType simulation_type = blocks[block].simulation_type;
int block_size = blocks[block].size;
int block_nb_recursives = blocks[block].getRecursiveSize();
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
int nb_recursives = blocks[blk].getRecursiveSize();
map<int, expr_t> recursive_vars;
for (int i = 0; i < nb_recursives; i++)
{
for (int i = 0; i < block_nb_recursives; i++)
int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(blk, i)), 0);
if (getBlockEquationType(blk, i) == EquationType::evaluate_s)
recursive_vars[deriv_id] = getBlockEquationRenormalizedExpr(blk, i);
else
recursive_vars[deriv_id] = getBlockEquationExpr(blk, i);
}
assert(blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesSimple
&& blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesComplete);
int size = blocks[blk].size;
for (int eq = nb_recursives; eq < size; eq++)
{
int eq_orig = getBlockEquationID(blk, eq);
for (int var = nb_recursives; var < size; var++)
{
if (getBlockEquationType(block, i) == EquationType::evaluate_s)
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
else
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
}
auto Derivatives = get_Derivatives(block);
for (const auto &it : Derivatives)
{
int Deriv_type = it.second;
auto [lag, eq, var, eqr, varr] = it.first;
if (Deriv_type == 0)
first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
else if (Deriv_type == 1)
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
{
if (getBlockEquationType(block, eq) == EquationType::evaluate_s
&& eq < block_nb_recursives)
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else
first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
}
tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]);
int var_orig = getBlockVariableID(blk, var);
expr_t d1 = equations[eq_orig]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, var_orig), 0), recursive_vars);
if (d1 == Zero)
continue;
blocks_derivatives[blk][{ eq, var, 0 }] = d1;
}
}
else
{
for (int i = 0; i < block_nb_recursives; i++)
{
if (getBlockEquationType(block, i) == EquationType::evaluate_s)
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
else
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
}
for (int eq = block_nb_recursives; eq < block_size; eq++)
{
int eqr = getBlockEquationID(block, eq);
for (int var = block_nb_recursives; var < block_size; var++)
{
int varr = getBlockVariableID(block, var);
expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), 0), recursive_variables);
if (d1 == Zero)
continue;
first_chain_rule_derivatives[{ eqr, varr, 0 }] = d1;
tmp_derivatives.emplace_back(eq, var, 0, first_chain_rule_derivatives[{ eqr, varr, 0 }]);
}
}
}
blocks_derivatives[block] = tmp_derivatives;
}
}

4
src/StaticModel.hh
View File

@ -77,7 +77,7 @@ private:
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
//! Get the type corresponding to a derivation ID
SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
@ -87,8 +87,6 @@ private:
int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
//! Compute the column indices of the static Jacobian
void computeStatJacobianCols();
//! 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();