diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 767f4cb58aed85473b0e38b060c78f0bd9807503..5fe95c9432a88eea3c5c77a45e35bbfe195a927b 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 575c72aefa53a9616b6834297e5193d721ff6491..b7973aa4a728a988ec8c8e14973219801090b6f4 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 7f1e5d4a520c5bf757b12f4f3f0e5fdb5d7871d0..14f776a565fa435e3edf6efa3015dffc833ae22d 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 a0914aff03020fe2c33bed812f0b613498ce096c..e9b9d4a9055a5d0aa669260781f129e4f20691cf 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)
{