Bytecode w/ block decomposition: no longer output derivatives w.r.t. exogenous and endogenous outside the block

master
Sébastien Villemot 2023-01-17 13:47:02 +01:00
parent 1e377f061a
commit 2e09df90e7
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 16 additions and 262 deletions

View File

@ -88,18 +88,12 @@ operator<<(BytecodeWriter &code_file, const FBEGINBLOCK_ &instr)
}
code_file.write(reinterpret_cast<const char *>(&instr.nb_col_jacob), sizeof instr.nb_col_jacob);
code_file.write(reinterpret_cast<const char *>(&instr.det_exo_size), sizeof instr.det_exo_size);
code_file.write(reinterpret_cast<const char *>(&instr.nb_col_det_exo_jacob), sizeof instr.nb_col_det_exo_jacob);
code_file.write(reinterpret_cast<const char *>(&instr.exo_size), sizeof instr.exo_size);
code_file.write(reinterpret_cast<const char *>(&instr.nb_col_exo_jacob), sizeof instr.nb_col_exo_jacob);
code_file.write(reinterpret_cast<const char *>(&instr.other_endo_size), sizeof instr.other_endo_size);
code_file.write(reinterpret_cast<const char *>(&instr.nb_col_other_endo_jacob), sizeof instr.nb_col_other_endo_jacob);
for (int i{0}; i < instr.det_exo_size; i++)
code_file.write(reinterpret_cast<const char *>(&instr.det_exogenous[i]), sizeof instr.det_exogenous[i]);
for (int i{0}; i < instr.exo_size; i++)
code_file.write(reinterpret_cast<const char *>(&instr.exogenous[i]), sizeof instr.exogenous[i]);
for (int i{0}; i < instr.other_endo_size; i++)
code_file.write(reinterpret_cast<const char *>(&instr.other_endogenous[i]), sizeof instr.other_endogenous[i]);
return code_file;
}

View File

@ -94,7 +94,6 @@ enum class ExpressionType
TemporaryTerm,
ModelEquation,
FirstEndoDerivative,
FirstOtherEndoDerivative,
FirstExoDerivative,
FirstExodetDerivative,
};
@ -895,7 +894,6 @@ private:
BlockSimulationType type;
vector<int> variable;
vector<int> equation;
vector<int> other_endogenous;
vector<int> exogenous;
vector<int> det_exogenous;
bool is_linear{false};
@ -905,24 +903,25 @@ private:
int Max_Lead{0};
int u_count_int{0};
int nb_col_jacob{0};
int det_exo_size, exo_size, other_endo_size;
int nb_col_det_exo_jacob, nb_col_exo_jacob, nb_col_other_endo_jacob;
int det_exo_size, exo_size;
public:
FBEGINBLOCK_() : BytecodeInstruction{Tags::FBEGINBLOCK},
type{BlockSimulationType::unknown}
{
}
/* Constructor when derivatives w.r.t. exogenous are present (only makes
sense when there is no block-decomposition, since there is no provision for
derivatives w.r.t. endogenous not belonging to the block) */
FBEGINBLOCK_(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
const vector<int> &variable_arg, const vector<int> &equation_arg,
bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg, int nb_col_jacob_arg,
int det_exo_size_arg, int nb_col_det_exo_jacob_arg, int exo_size_arg, int nb_col_exo_jacob_arg, int other_endo_size_arg, int nb_col_other_endo_jacob_arg,
vector<int> det_exogenous_arg, vector<int> exogenous_arg, vector<int> other_endogenous_arg) :
int det_exo_size_arg, int exo_size_arg,
vector<int> det_exogenous_arg, vector<int> exogenous_arg) :
BytecodeInstruction{Tags::FBEGINBLOCK},
size{size_arg},
type{type_arg},
variable{variable_arg.begin()+first_element, variable_arg.begin()+(first_element+block_size)},
equation{equation_arg.begin()+first_element, equation_arg.begin()+(first_element+block_size)},
other_endogenous{move(other_endogenous_arg)},
exogenous{move(exogenous_arg)},
det_exogenous{move(det_exogenous_arg)},
is_linear{is_linear_arg},
@ -932,13 +931,10 @@ public:
u_count_int{u_count_int_arg},
nb_col_jacob{nb_col_jacob_arg},
det_exo_size{det_exo_size_arg},
exo_size{exo_size_arg},
other_endo_size{other_endo_size_arg},
nb_col_det_exo_jacob{nb_col_det_exo_jacob_arg},
nb_col_exo_jacob{nb_col_exo_jacob_arg},
nb_col_other_endo_jacob{nb_col_other_endo_jacob_arg}
exo_size{exo_size_arg}
{
}
// Constructor when derivatives w.r.t. exogenous are absent
FBEGINBLOCK_(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
const vector<int> &variable_arg, const vector<int> &equation_arg,
bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg, int nb_col_jacob_arg) :
@ -954,11 +950,7 @@ public:
u_count_int{u_count_int_arg},
nb_col_jacob{nb_col_jacob_arg},
det_exo_size{0},
exo_size{0},
other_endo_size{0},
nb_col_det_exo_jacob{0},
nb_col_exo_jacob{0},
nb_col_other_endo_jacob{0}
exo_size{0}
{
}
int
@ -1012,30 +1004,10 @@ public:
return exo_size;
};
int
get_nb_col_exo_jacob()
{
return nb_col_exo_jacob;
};
int
get_det_exo_size()
{
return det_exo_size;
};
int
get_nb_col_det_exo_jacob()
{
return nb_col_det_exo_jacob;
};
int
get_other_endo_size()
{
return other_endo_size;
};
int
get_nb_col_other_endo_jacob()
{
return nb_col_other_endo_jacob;
};
vector<int>
get_endogenous()
{
@ -1074,11 +1046,7 @@ public:
}
memcpy(&nb_col_jacob, code, sizeof(nb_col_jacob)); code += sizeof(nb_col_jacob);
memcpy(&det_exo_size, code, sizeof(det_exo_size)); code += sizeof(det_exo_size);
memcpy(&nb_col_det_exo_jacob, code, sizeof(nb_col_det_exo_jacob)); code += sizeof(nb_col_det_exo_jacob);
memcpy(&exo_size, code, sizeof(exo_size)); code += sizeof(exo_size);
memcpy(&nb_col_exo_jacob, code, sizeof(nb_col_exo_jacob)); code += sizeof(nb_col_exo_jacob);
memcpy(&other_endo_size, code, sizeof(other_endo_size)); code += sizeof(other_endo_size);
memcpy(&nb_col_other_endo_jacob, code, sizeof(nb_col_other_endo_jacob)); code += sizeof(nb_col_other_endo_jacob);
for (int i{0}; i < det_exo_size; i++)
{
@ -1092,12 +1060,6 @@ public:
memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
exogenous.push_back(tmp_i);
}
for (int i{0}; i < other_endo_size; i++)
{
int tmp_i;
memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
other_endogenous.push_back(tmp_i);
}
return code;
};
#endif

View File

@ -37,20 +37,6 @@ DynamicModel::copyHelper(const DynamicModel &m)
for (const auto &it : m.static_only_equations)
static_only_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
auto convert_block_derivative = [f](const map<tuple<int, int, int>, expr_t> &dt)
{
map<tuple<int, int, int>, expr_t> dt2;
for (const auto &it : dt)
dt2.emplace(it.first, f(it.second));
return dt2;
};
for (const auto &it : m.blocks_derivatives_other_endo)
blocks_derivatives_other_endo.emplace_back(convert_block_derivative(it));
for (const auto &it : m.blocks_derivatives_exo)
blocks_derivatives_exo.emplace_back(convert_block_derivative(it));
for (const auto &it : m.blocks_derivatives_exo_det)
blocks_derivatives_exo_det.emplace_back(convert_block_derivative(it));
}
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
@ -98,13 +84,7 @@ DynamicModel::DynamicModel(const DynamicModel &m) :
xref_exo_det{m.xref_exo_det},
nonzero_hessian_eqs{m.nonzero_hessian_eqs},
variableMapping{m.variableMapping},
blocks_other_endo{m.blocks_other_endo},
blocks_exo{m.blocks_exo},
blocks_exo_det{m.blocks_exo_det},
blocks_jacob_cols_endo{m.blocks_jacob_cols_endo},
blocks_jacob_cols_other_endo{m.blocks_jacob_cols_other_endo},
blocks_jacob_cols_exo{m.blocks_jacob_cols_exo},
blocks_jacob_cols_exo_det{m.blocks_jacob_cols_exo_det},
var_expectation_functions_to_write{m.var_expectation_functions_to_write}
{
copyHelper(m);
@ -148,16 +128,7 @@ DynamicModel::operator=(const DynamicModel &m)
xref_exo_det = m.xref_exo_det;
nonzero_hessian_eqs = m.nonzero_hessian_eqs;
variableMapping = m.variableMapping;
blocks_derivatives_other_endo.clear();
blocks_derivatives_exo.clear();
blocks_derivatives_exo_det.clear();
blocks_other_endo = m.blocks_other_endo;
blocks_exo = m.blocks_exo;
blocks_exo_det = m.blocks_exo_det;
blocks_jacob_cols_endo = m.blocks_jacob_cols_endo;
blocks_jacob_cols_other_endo = m.blocks_jacob_cols_other_endo;
blocks_jacob_cols_exo = m.blocks_jacob_cols_exo;
blocks_jacob_cols_exo_det = m.blocks_jacob_cols_exo_det;
var_expectation_functions_to_write = m.var_expectation_functions_to_write;
@ -166,52 +137,6 @@ DynamicModel::operator=(const DynamicModel &m)
return *this;
}
void
DynamicModel::additionalBlockTemporaryTerms(int blk,
vector<vector<temporary_terms_t>> &blocks_temporary_terms,
map<expr_t, tuple<int, int, int>> &reference_count) const
{
for (const auto &[ignore, d] : blocks_derivatives_exo[blk])
d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
for (const auto &[ignore, d] : blocks_derivatives_exo_det[blk])
d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
for (const auto &[ignore, d] : blocks_derivatives_other_endo[blk])
d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
}
void
DynamicModel::writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
const temporary_terms_t &temporary_terms_union,
const deriv_node_temp_terms_t &tef_terms) const
{
constexpr ExprNodeBytecodeOutputType output_type {ExprNodeBytecodeOutputType::dynamicModel};
/* FIXME: there is an inconsistency between endos and the following 3 other
variable types. For the latter, the index of equation within the block is
taken from FNUMEXPR, while it is taken from FSTPG3 for the former. */
for (const auto &[indices, d] : blocks_derivatives_exo[block])
{
const auto &[eq, var, lag] {indices};
code_file << FNUMEXPR_{ExpressionType::FirstExoDerivative, eq, 0, lag};
d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo[block].at({ var, lag })};
}
for (const auto &[indices, d] : blocks_derivatives_exo_det[block])
{
const auto &[eq, var, lag] {indices};
code_file << FNUMEXPR_{ExpressionType::FirstExodetDerivative, eq, 0, lag};
d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo_det[block].at({ var, lag })};
}
for (const auto &[indices, d] : blocks_derivatives_other_endo[block])
{
const auto &[eq, var, lag] {indices};
code_file << FNUMEXPR_{ExpressionType::FirstOtherEndoDerivative, eq, 0, lag};
d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag })};
}
}
void
DynamicModel::writeDynamicBytecode(const string &basename) const
{
@ -244,8 +169,6 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
[this](const auto &v)
{ return getTypeByDerivID(v.first) == SymbolType::endogenous; }))
};
int jacobian_ncols_exo {symbol_table.exo_nbr()};
int jacobian_ncols_exo_det {symbol_table.exo_det_nbr()};
vector<int> eq_idx(equations.size());
iota(eq_idx.begin(), eq_idx.end(), 0);
vector<int> endo_idx(symbol_table.endo_nbr());
@ -264,14 +187,9 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
u_count_int,
jacobian_ncols_endo,
symbol_table.exo_det_nbr(),
jacobian_ncols_exo_det,
symbol_table.exo_nbr(),
jacobian_ncols_exo,
0,
0,
exo_det,
exo,
{}};
exo};
writeBytecodeHelper<true>(code_file);
}
@ -315,16 +233,7 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
blocks[block].max_lag,
blocks[block].max_lead,
u_count,
static_cast<int>(blocks_jacob_cols_endo[block].size()),
static_cast<int>(blocks_exo_det[block].size()),
static_cast<int>(blocks_jacob_cols_exo_det[block].size()),
static_cast<int>(blocks_exo[block].size()),
static_cast<int>(blocks_jacob_cols_exo[block].size()),
static_cast<int>(blocks_other_endo[block].size()),
static_cast<int>(blocks_jacob_cols_other_endo[block].size()),
{ blocks_exo_det[block].begin(), blocks_exo_det[block].end() },
{ blocks_exo[block].begin(), blocks_exo[block].end() },
{ blocks_other_endo[block].begin(), blocks_other_endo[block].end() }};
static_cast<int>(blocks_jacob_cols_endo[block].size())};
writeBlockBytecodeHelper<true>(code_file, block);
}
@ -2506,43 +2415,6 @@ DynamicModel::computeChainRuleJacobian()
blocks_jacobian_sparse_colptr[blk] = computeCSCColPtr(blocks_jacobian_sparse_column_major_order[blk], (one_boundary ? 1 : 3)*mfs_size);
}
}
/* Also store information and derivatives w.r.t. other types of variables
(for the stochastic mode) */
blocks_derivatives_other_endo.resize(nb_blocks);
blocks_derivatives_exo.resize(nb_blocks);
blocks_derivatives_exo_det.resize(nb_blocks);
blocks_other_endo.resize(nb_blocks);
blocks_exo.resize(nb_blocks);
blocks_exo_det.resize(nb_blocks);
for (auto &[indices, d1] : derivatives[1])
{
auto [eq_orig, deriv_id] { vectorToTuple<2>(indices) };
int block_eq { eq2block[eq_orig] };
int eq { getBlockInitialEquationID(block_eq, eq_orig) };
int var { getTypeSpecificIDByDerivID(deriv_id) };
int lag { getLagByDerivID(deriv_id) };
switch (getTypeByDerivID(indices[1]))
{
case SymbolType::endogenous:
if (block_eq != endo2block[var])
{
blocks_derivatives_other_endo[block_eq][{ eq, var, lag }] = d1;
blocks_other_endo[block_eq].insert(var);
}
break;
case SymbolType::exogenous:
blocks_derivatives_exo[block_eq][{ eq, var, lag }] = d1;
blocks_exo[block_eq].insert(var);
break;
case SymbolType::exogenousDet:
blocks_derivatives_exo_det[block_eq][{ eq, var, lag }] = d1;
blocks_exo_det[block_eq].insert(var);
break;
default:
break;
}
}
}
void
@ -2580,27 +2452,10 @@ DynamicModel::computeBlockDynJacobianCols()
// Compute Jacobian column indices
blocks_jacob_cols_endo.resize(nb_blocks);
blocks_jacob_cols_other_endo.resize(nb_blocks);
blocks_jacob_cols_exo.resize(nb_blocks);
blocks_jacob_cols_exo_det.resize(nb_blocks);
for (size_t blk {0}; blk < nb_blocks; blk++)
{
for (int index{0};
auto [lag, var] : dynamic_endo[blk])
blocks_jacob_cols_endo[blk][{ var, lag }] = index++;
for (int index{0};
auto [lag, var] : dynamic_other_endo[blk])
blocks_jacob_cols_other_endo[blk][{ var, lag }] = index++;
for (int index{0};
auto [lag, var] : dynamic_exo[blk])
blocks_jacob_cols_exo[blk][{ var, lag }] = index++;
for (int index{0};
auto [lag, var] : dynamic_exo_det[blk])
blocks_jacob_cols_exo_det[blk][{ var, lag }] = index++;
}
for (int index{0};
auto [lag, var] : dynamic_endo[blk])
blocks_jacob_cols_endo[blk][{ var, lag }] = index++;
}
void

View File

@ -100,20 +100,9 @@ private:
//! Creates mapping for variables and equations they are present in
map<int, set<int>> variableMapping;
/* Derivatives of block equations with respect to: endogenous that do not
belong to the block, exogenous, deterministic exogenous.
Tuples are of the form (equation no. within the block, type-specific ID, lag) */
vector<map<tuple<int, int, int>, expr_t>> blocks_derivatives_other_endo,
blocks_derivatives_exo, blocks_derivatives_exo_det;
// For each block, gives type-specific other endos / exo / exo det that appear in it
vector<set<int>> blocks_other_endo, blocks_exo, blocks_exo_det;
/* For each block, and for each variable type, maps (variable ID, lag) to
Jacobian column.
For the endo version, the variable ID is the index within the block. For
the three others, its the type-specific ID */
vector<map<pair<int, int>, int>> blocks_jacob_cols_endo, blocks_jacob_cols_other_endo, blocks_jacob_cols_exo, blocks_jacob_cols_exo_det;
Jacobian column. The variable ID is the index within the block. */
vector<map<pair<int, int>, int>> blocks_jacob_cols_endo;
//! Used for var_expectation and var_model
map<string, set<int>> var_expectation_functions_to_write;
@ -122,10 +111,6 @@ private:
void writeDynamicMFile(const string &basename) const;
//! Writes the code of the block-decomposed model in virtual machine bytecode
void writeDynamicBlockBytecode(const string &basename) const;
// Writes derivatives w.r.t. exo, exo det and other endogenous
void writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
const temporary_terms_t &temporary_terms_union,
const deriv_node_temp_terms_t &tef_terms) const override;
//! Writes the code of the model in virtual machine bytecode
void writeDynamicBytecode(const string &basename) const;
@ -152,10 +137,6 @@ private:
string reform(const string &name) const;
void additionalBlockTemporaryTerms(int blk,
vector<vector<temporary_terms_t>> &blocks_temporary_terms,
map<expr_t, tuple<int, int, int>> &reference_count) const override;
SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
int getLagByDerivID(int deriv_id) const noexcept(false) override;
int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;

View File

@ -1014,8 +1014,6 @@ ModelTree::computeBlockTemporaryTerms(bool no_tmp_terms)
}
for (const auto &[ignore, d] : blocks_derivatives[blk])
d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
additionalBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
}
/* If the user has specified the notmpterms option, clear all temporary
@ -1035,13 +1033,6 @@ ModelTree::computeBlockTemporaryTerms(bool no_tmp_terms)
blocks_temporary_terms_idxs[tt] = idx++;
}
void
ModelTree::additionalBlockTemporaryTerms([[maybe_unused]] int blk,
[[maybe_unused]] vector<vector<temporary_terms_t>> &blocks_temporary_terms,
[[maybe_unused]] map<expr_t, tuple<int, int, int>> &reference_count) const
{
}
void
ModelTree::writeJsonTemporaryTerms(const temporary_terms_t &tt,
temporary_terms_t &temp_term_union,
@ -1903,14 +1894,6 @@ ModelTree::getRHSFromLHS(expr_t lhs) const
throw ExprNode::MatchFailureException{"Cannot find an equation with the requested LHS"};
}
void
ModelTree::writeBlockBytecodeAdditionalDerivatives([[maybe_unused]] BytecodeWriter &code_file,
[[maybe_unused]] int block,
[[maybe_unused]] const temporary_terms_t &temporary_terms_union,
[[maybe_unused]] const deriv_node_temp_terms_t &tef_terms) const
{
}
void
ModelTree::initializeMEXCompilationWorkers(int numworkers)
{

View File

@ -279,16 +279,6 @@ protected:
//! Computes temporary terms per block
void computeBlockTemporaryTerms(bool no_tmp_terms);
private:
/* Add additional temporary terms for a given block. This method is called by
computeBlockTemporaryTerms(). It does nothing by default, but is meant to
be overriden by subclasses (actually by DynamicModel, who needs extra
temporary terms for derivatives w.r.t. exogenous and other endogenous) */
virtual void additionalBlockTemporaryTerms(int blk,
vector<vector<temporary_terms_t>> &blocks_temporary_terms,
map<expr_t, tuple<int, int, int>> &reference_count) const;
protected:
//! Computes temporary terms for the file containing parameters derivatives
void computeParamsDerivativesTemporaryTerms();
//! Writes temporary terms
@ -349,13 +339,6 @@ protected:
template<bool dynamic>
void writeBlockBytecodeHelper(BytecodeWriter &code_file, int block) const;
/* Write additional derivatives w.r.t. to exogenous, exogenous det and other endo
in block+bytecode mode. Does nothing by default, but overriden by
DynamicModel which needs those. */
virtual void writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
const temporary_terms_t &temporary_terms_union,
const deriv_node_temp_terms_t &tef_terms) const;
// Helper for writing sparse derivatives indices in MATLAB/Octave driver file
template<bool dynamic>
void writeDriverSparseIndicesHelper(ostream &output) const;
@ -1851,10 +1834,6 @@ ModelTree::writeBlockBytecodeHelper(BytecodeWriter &code_file, int block) const
code_file << FSTPG2_{eq, getBlockJacobianEndoCol(block, var, lag)};
}
/* Write derivatives w.r.t. exo, exo det and other endogenous, but only in
dynamic mode */
writeBlockBytecodeAdditionalDerivatives(code_file, block, temporary_terms_union, tef_terms);
// Update jump offset for previous JMP
int pos_end_block {code_file.getInstructionCounter()};
code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1});