Commit a0af8f7e authored by sebastien's avatar sebastien
Browse files

trunk preprocessor: restructuring the way we decide which derivatives to compute


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2610 ac1d8469-bf42-47a9-8791-bf33cf982152
parent d99229ca
......@@ -907,8 +907,7 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct)
void
PlannerObjectiveStatement::computingPass()
{
model_tree->computeStaticHessian = true;
model_tree->computingPass(false);
model_tree->computingPass(true, false);
}
void
......
This diff is collapsed.
......@@ -39,10 +39,6 @@ private:
/*! Contains only endogenous, exogenous and exogenous deterministic */
map<int, int> dyn_jacobian_cols_table;
//! Number of dynamic endogenous variables inside the model block
/*! Set by computeDerivID() */
int var_endo_nbr;
//! Maximum lag and lead over all types of variables (positive values)
/*! Set by computeDerivID() */
int max_lag, max_lead;
......@@ -56,6 +52,10 @@ private:
/*! Set by computeDerivID() */
int max_exo_det_lag, max_exo_det_lead;
//! Number of columns of dynamic jacobian
/*! Set by computeDerivID() and computeDynJacobianCols() */
int dynJacobianColsNbr;
//! Writes dynamic model file (Matlab version)
void writeDynamicMFile(const string &dynamic_basename) const;
//! Writes dynamic model file (C version)
......@@ -82,21 +82,19 @@ private:
//! Build The incidence matrix form the modeltree
void BuildIncidenceMatrix();
void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock);
void computeTemporaryTermsOrdered(Model_Block *ModelBlock);
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const;
virtual int computeDerivID(int symb_id, int lag);
//! Get the number of columns of dynamic jacobian
int getDynJacobianColsNbr() const;
//! Get the type corresponding to a derivation ID
int getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the lag corresponding to a derivation ID
int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the symbol ID corresponding to a derivation ID
int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the dynamic Jacobian
void computeDynJacobianCols();
void computeDynJacobianCols(bool jacobianExo);
public:
DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
......@@ -109,18 +107,16 @@ public:
double markowitz;
//! the file containing the model and the derivatives code
ofstream code_file;
//! Whether dynamic Jacobian (w.r. to endogenous) should be written
bool computeJacobian;
//! Whether dynamic Jacobian (w.r. to endogenous and exogenous) should be written
bool computeJacobianExo;
//! Whether dynamic Hessian (w.r. to endogenous and exogenous) should be written
bool computeHessian;
//! Whether dynamic third order derivatives (w.r. to endogenous and exogenous) should be written
bool computeThirdDerivatives;
//! Execute computations (variable sorting + derivation)
/*! You must set computeJacobian, computeJacobianExo, computeHessian and computeThirdDerivatives to correct values before calling this function
\param no_tmp_terms if true, no temporary terms will be computed in the dynamic files */
void computingPass(const eval_context_type &eval_context, bool no_tmp_terms);
/*!
\param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
\param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
\param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
\param eval_context evaluation context for normalization
\param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
*/
void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives,
const eval_context_type &eval_context, bool no_tmp_terms);
//! Writes model initialization and lead/lag incidence matrix to output
void writeOutput(ostream &output) const;
//! Complete set to block decompose the model
......
......@@ -144,11 +144,12 @@ ModFile::computingPass(bool no_tmp_terms)
{
// Compute static model and its derivatives
dynamic_model.toStatic(static_model);
static_model.computingPass(no_tmp_terms);
static_model.computingPass(false, no_tmp_terms);
// Set things to compute for dynamic model
if (mod_file_struct.simul_present)
dynamic_model.computeJacobian = true;
dynamic_model.computingPass(false, false, false, global_eval_context, no_tmp_terms);
else
{
if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3)
......@@ -156,14 +157,12 @@ ModFile::computingPass(bool no_tmp_terms)
cerr << "Incorrect order option..." << endl;
exit(EXIT_FAILURE);
}
dynamic_model.computeJacobianExo = true;
if (mod_file_struct.order_option >= 2)
dynamic_model.computeHessian = true;
if (mod_file_struct.order_option == 3)
dynamic_model.computeThirdDerivatives = true;
bool hessian = mod_file_struct.order_option >= 2;
bool thirdDerivatives = mod_file_struct.order_option == 3;
dynamic_model.computingPass(true, hessian, thirdDerivatives, global_eval_context, no_tmp_terms);
}
dynamic_model.computingPass(global_eval_context, no_tmp_terms);
}
for(vector<Statement *>::iterator it = statements.begin();
it != statements.end(); it++)
(*it)->computingPass();
......
......@@ -48,65 +48,77 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
}
void
ModelTree::derive(int order)
ModelTree::computeJacobian(const set<int> &vars)
{
for (int var = 0; var < getDerivIDNbr(); var++)
for(set<int>::const_iterator it = vars.begin();
it != vars.end(); it++)
for (int eq = 0; eq < (int) equations.size(); eq++)
{
NodeID d1 = equations[eq]->getDerivative(var);
NodeID d1 = equations[eq]->getDerivative(*it);
if (d1 == Zero)
continue;
first_derivatives[make_pair(eq, var)] = d1;
first_derivatives[make_pair(eq, *it)] = d1;
}
}
if (order >= 2)
void
ModelTree::computeHessian(const set<int> &vars)
{
for (first_derivatives_type::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
for (first_derivatives_type::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
int eq = it->first.first;
int var1 = it->first.second;
NodeID d1 = it->second;
// Store only second derivatives with var2 <= var1
for(set<int>::const_iterator it2 = vars.begin();
it2 != vars.end(); it2++)
{
int eq = it->first.first;
int var1 = it->first.second;
NodeID d1 = it->second;
// Store only second derivatives with var2 <= var1
for (int var2 = 0; var2 <= var1; var2++)
{
NodeID d2 = d1->getDerivative(var2);
if (d2 == Zero)
continue;
second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
}
int var2 = *it2;
if (var2 > var1)
continue;
NodeID d2 = d1->getDerivative(var2);
if (d2 == Zero)
continue;
second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
}
}
}
if (order >= 3)
void
ModelTree::computeThirdDerivatives(const set<int> &vars)
{
for (second_derivatives_type::const_iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
{
for (second_derivatives_type::const_iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
int eq = it->first.first;
int var1 = it->first.second.first;
int var2 = it->first.second.second;
// By construction, var2 <= var1
NodeID d2 = it->second;
// Store only third derivatives such that var3 <= var2 <= var1
for(set<int>::const_iterator it2 = vars.begin();
it2 != vars.end(); it2++)
{
int eq = it->first.first;
int var1 = it->first.second.first;
int var2 = it->first.second.second;
// By construction, var2 <= var1
NodeID d2 = it->second;
// Store only third derivatives such that var3 <= var2 <= var1
for (int var3 = 0; var3 <= var2; var3++)
{
NodeID d3 = d2->getDerivative(var3);
if (d3 == Zero)
continue;
third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
}
int var3 = *it2;
if (var3 > var2)
continue;
NodeID d3 = d2->getDerivative(var3);
if (d3 == Zero)
continue;
third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
}
}
}
void
ModelTree::computeTemporaryTerms(int order)
ModelTree::computeTemporaryTerms()
{
map<NodeID, int> reference_count;
temporary_terms.clear();
......@@ -121,15 +133,13 @@ ModelTree::computeTemporaryTerms(int order)
it != first_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
if (order >= 2)
for (second_derivatives_type::iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
for (second_derivatives_type::iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
if (order >= 3)
for (third_derivatives_type::iterator it = third_derivatives.begin();
it != third_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
for (third_derivatives_type::iterator it = third_derivatives.begin();
it != third_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
}
void
......
......@@ -74,12 +74,20 @@ protected:
//! Temporary terms (those which will be noted Txxxx)
temporary_terms_type temporary_terms;
//! Computes derivatives of ModelTree
void derive(int order);
//! Computes 1st derivatives
/*! \param vars the derivation IDs w.r. to which compute the derivatives */
void computeJacobian(const set<int> &vars);
//! Computes 2nd derivatives
/*! \param vars the derivation IDs w.r. to which derive the 1st derivatives */
void computeHessian(const set<int> &vars);
//! Computes 3rd derivatives
/*! \param vars the derivation IDs w.r. to which derive the 2nd derivatives */
void computeThirdDerivatives(const set<int> &vars);
//! Write derivative of an equation w.r. to a variable
void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
//! Computes temporary terms
void computeTemporaryTerms(int order);
//! Computes temporary terms (for all equations and derivatives)
void computeTemporaryTerms();
//! Writes temporary terms
void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const;
//! Writes model local variables
......
......@@ -23,8 +23,7 @@
StaticModel::StaticModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :
ModelTree(symbol_table_arg, num_constants_arg),
computeStaticHessian(false)
ModelTree(symbol_table_arg, num_constants_arg)
{
}
......@@ -151,38 +150,37 @@ StaticModel::writeStaticModel(ostream &StaticOutput) const
jacobian_output << ";" << endl;
}
// Write Hessian w.r. to endogenous only
if (computeStaticHessian)
for (second_derivatives_type::const_iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
{
int eq = it->first.first;
int symb_id1 = inv_deriv_id_table[it->first.second.first];
int symb_id2 = inv_deriv_id_table[it->first.second.second];
NodeID d2 = it->second;
int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
hessian_output << " g2";
matrixHelper(hessian_output, eq, col_nb, output_type);
hessian_output << " = ";
d2->writeOutput(hessian_output, output_type, temporary_terms);
hessian_output << ";" << endl;
// Treating symetric elements
if (symb_id1 != symb_id2)
{
lsymetric << " g2";
matrixHelper(lsymetric, eq, col_nb_sym, output_type);
lsymetric << " = " << "g2";
matrixHelper(lsymetric, eq, col_nb, output_type);
lsymetric << ";" << endl;
}
}
// Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
for (second_derivatives_type::const_iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
{
int eq = it->first.first;
int symb_id1 = inv_deriv_id_table[it->first.second.first];
int symb_id2 = inv_deriv_id_table[it->first.second.second];
NodeID d2 = it->second;
int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
hessian_output << " g2";
matrixHelper(hessian_output, eq, col_nb, output_type);
hessian_output << " = ";
d2->writeOutput(hessian_output, output_type, temporary_terms);
hessian_output << ";" << endl;
// Treating symetric elements
if (symb_id1 != symb_id2)
{
lsymetric << " g2";
matrixHelper(lsymetric, eq, col_nb_sym, output_type);
lsymetric << " = " << "g2";
matrixHelper(lsymetric, eq, col_nb, output_type);
lsymetric << ";" << endl;
}
}
// Writing ouputs
if (mode != eDLLMode)
......@@ -208,9 +206,11 @@ StaticModel::writeStaticModel(ostream &StaticOutput) const
<< " g1 = real(g1)+2*imag(g1);" << endl
<< " end" << endl
<< "end" << endl;
if (computeStaticHessian)
// If 2nd order derivatives have been computed
if (second_derivatives.size())
{
StaticOutput << "if nargout >= 3,\n";
StaticOutput << "if nargout >= 3," << endl;
// Writing initialization instruction for matrix g2
int ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
StaticOutput << " g2 = sparse([],[],[], " << equations.size() << ", " << ncols << ", " << 5*ncols << ");" << endl
......@@ -266,20 +266,26 @@ StaticModel::writeStaticFile(const string &basename) const
}
void
StaticModel::computingPass(bool no_tmp_terms)
StaticModel::computingPass(bool hessian, bool no_tmp_terms)
{
// Determine derivation order
int order = 1;
if (computeStaticHessian)
order = 2;
// Compute derivatives w.r. to all derivation IDs (i.e. all endogenous)
set<int> vars;
for(int i = 0; i < getDerivIDNbr(); i++)
vars.insert(i);
// Launch computations
cout << "Computing static model derivatives at order " << order << "...";
derive(order);
cout << " done" << endl;
cout << "Computing static model derivatives:" << endl
<< " - order 1" << endl;
computeJacobian(vars);
if (hessian)
{
cout << " - order 2" << endl;
computeHessian(vars);
}
if (!no_tmp_terms)
computeTemporaryTerms(order);
computeTemporaryTerms();
}
int
......
......@@ -46,12 +46,11 @@ private:
public:
StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Whether static Hessian (w.r. to endogenous only) should be written
bool computeStaticHessian;
//! Execute computations (derivation)
/*! You must set computeStaticHessian before calling this function
\param hessian whether Hessian (w.r. to endogenous only) should be computed
\param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
void computingPass(bool no_tmp_terms);
void computingPass(bool hessian, bool no_tmp_terms);
//! Writes static model file
void writeStaticFile(const string &basename) const;
......
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