From 2af26ee2c2815f2ec9d75a48fd562b8e90d37e2c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 13 Dec 2019 17:15:03 +0100
Subject: [PATCH] Ramsey: use information from transformed model for filling
M_.nonzero_hessian_eqs
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Since commit 9c9e8f816f056fb4f9f2b17234ba7c100f3a996a, it’s the information
from the original model which was in this field, which is not what is expected.
By the way, do not output this field (and the related M_.hessian_eq_zero) when
the Hessian is not computed by the preprocessor (i.e. in practice for perfect
foresight), since they would otherwise contain incorrect information.
Ref. dynare#1681
---
src/DynamicModel.cc | 30 +++++++-------------------
src/DynamicModel.hh | 18 +++++++++++++---
src/ModFile.cc | 51 ++++++++++++++++++++++++++-------------------
src/ModelTree.hh | 11 ++++++++++
4 files changed, 63 insertions(+), 47 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 767f4cb5..5fe95c94 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1875,33 +1875,12 @@ DynamicModel::printNonZeroHessianEquations(ostream &output) const
{
if (it != nonzero_hessian_eqs.begin())
output << " ";
- output << it->first;
+ output << *it;
}
if (nonzero_hessian_eqs.size() != 1)
output << "]";
}
-void
-DynamicModel::setNonZeroHessianEquations(map<int, string> &eqs)
-{
- for (const auto &it : derivatives[2])
- {
- int eq = it.first[0];
- if (nonzero_hessian_eqs.find(eq) == nonzero_hessian_eqs.end())
- {
- nonzero_hessian_eqs[eq] = "";
- for (auto & equation_tag : equation_tags)
- if (equation_tag.first == eq)
- if (equation_tag.second.first == "name")
- {
- nonzero_hessian_eqs[eq] = equation_tag.second.second;
- break;
- }
- }
- }
- eqs = nonzero_hessian_eqs;
-}
-
void
DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num,
int &u_count_int, bool &file_open, bool is_two_boundaries, const bool linear_decomposition) const
@@ -5003,6 +4982,13 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
computeDerivatives(derivsOrder, vars);
+ if (derivsOrder > 1)
+ {
+ hessian_computed = true;
+ for (const auto &it : derivatives[2])
+ nonzero_hessian_eqs.insert(it.first[0]);
+ }
+
if (paramsDerivsOrder > 0)
{
cout << "Computing dynamic model derivatives w.r.t. parameters (order " << paramsDerivsOrder << ")." << endl;
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 575c72ae..b7973aa4 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -98,7 +98,9 @@ private:
map<pair<int, int>, set<int>> xref_exo_det;
//! Nonzero equations in the Hessian
- map<int, string> nonzero_hessian_eqs;
+ set<int> nonzero_hessian_eqs;
+ //! Whether the hessian has actually been computed
+ bool hessian_computed{false};
//! Number of columns of dynamic jacobian
/*! Set by computeDerivID()s and computeDynJacobianCols() */
@@ -354,8 +356,18 @@ public:
//! Print equations that have non-zero second derivatives
void printNonZeroHessianEquations(ostream &output) const;
- //! Set the equations that have non-zero second derivatives
- void setNonZeroHessianEquations(map<int, string> &eqs);
+ //! Tells whether Hessian has been computed
+ /*! This is needed to know whether no non-zero equation in Hessian means a
+ zero Hessian or Hessian not computed */
+ inline bool isHessianComputed() const
+ {
+ return hessian_computed;
+ }
+ //! Returns equations that have non-zero second derivatives
+ inline set<int> getNonZeroHessianEquations() const
+ {
+ return nonzero_hessian_eqs;
+ }
//! Fill Autoregressive Matrix for var_model
map<string, map<tuple<int, int, int>, expr_t>> fillAutoregressiveMatrix(bool is_var) const;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 7f1e5d4a..14f776a5 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -802,24 +802,31 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
else // No computing task requested, compute derivatives up to 2nd order by default
dynamic_model.computingPass(true, 2, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, linear_decomposition);
- map<int, string> eqs;
- if (mod_file_struct.ramsey_model_present)
- orig_ramsey_dynamic_model.setNonZeroHessianEquations(eqs);
- else
- dynamic_model.setNonZeroHessianEquations(eqs);
-
- if (linear && !eqs.empty())
+ /* Check that the model is linear.
+ FIXME: this check always passes if derivsOrder = 1, i.e. for a perfect
+ foresight model, because the Hessian is not computed in that case. */
+ if (linear)
{
- cerr << "ERROR: If the model is declared linear the second derivatives must be equal to zero." << endl
- << " The following equations had non-zero second derivatives:" << endl;
- for (map<int, string >::const_iterator it = eqs.begin(); it != eqs.end(); it++)
+ set<int> eqs;
+ if (mod_file_struct.ramsey_model_present)
+ eqs = orig_ramsey_dynamic_model.getNonZeroHessianEquations();
+ else
+ eqs = dynamic_model.getNonZeroHessianEquations();
+
+ if (!eqs.empty())
{
- cerr << " * Eq # " << it->first+1;
- if (!it->second.empty())
- cerr << " [" << it->second << "]";
- cerr << endl;
+ cerr << "ERROR: If the model is declared linear the second derivatives must be equal to zero." << endl
+ << " The following equations have non-zero second derivatives:" << endl;
+ for (const auto &it : eqs)
+ {
+ cerr << " * Eq # " << it+1;
+ auto tags = dynamic_model.getEquationTags(it);
+ if (auto it2 = tags.find("name"); it2 != tags.end())
+ cerr << " [" << it2->second << "]";
+ cerr << endl;
+ }
+ exit(EXIT_FAILURE);
}
- exit(EXIT_FAILURE);
}
}
@@ -992,13 +999,13 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
mOutputFile << "};" << endl;
}
- mOutputFile << "M_.nonzero_hessian_eqs = ";
- if (mod_file_struct.ramsey_model_present)
- orig_ramsey_dynamic_model.printNonZeroHessianEquations(mOutputFile);
- else
- dynamic_model.printNonZeroHessianEquations(mOutputFile);
- mOutputFile << ";" << endl
- << "M_.hessian_eq_zero = isempty(M_.nonzero_hessian_eqs);" << endl;
+ if (dynamic_model.isHessianComputed())
+ {
+ mOutputFile << "M_.nonzero_hessian_eqs = ";
+ dynamic_model.printNonZeroHessianEquations(mOutputFile);
+ mOutputFile << ";" << endl
+ << "M_.hessian_eq_zero = isempty(M_.nonzero_hessian_eqs);" << endl;
+ }
if (!onlymodel)
config_file.writeCluster(mOutputFile);
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index a0914aff..e9b9d4a9 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -382,6 +382,17 @@ public:
/*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[2]]
If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[3]] */
void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
+
+ //! Returns all the equation tags associated to an equation
+ inline map<string, string> getEquationTags(int eq) const
+ {
+ map<string, string> r;
+ for (auto &[eq2, tagpair]: equation_tags)
+ if (eq2 == eq)
+ r[tagpair.first] = tagpair.second;
+ return r;
+ }
+
inline static std::string
c_Equation_Type(int type)
{
--
GitLab