From 22000da5eac20dec7c61511422a24dc8e6199ac8 Mon Sep 17 00:00:00 2001
From: sebastien <sebastien@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Mon, 21 Dec 2009 10:29:21 +0000
Subject: [PATCH] Changes by Ferhat: * fix options stack_solve_algo={2,3,4}
(closes #68) * fix crashes for singular normalizations (closes #44) and
implement decreasing cutoff * fail for stack_solve_algo=2 under Octave
(because there is no gmres function in Octave)
git-svn-id: https://www.dynare.org/svn/dynare/trunk@3279 ac1d8469-bf42-47a9-8791-bf33cf982152
---
DynamicModel.cc | 2 +-
ModelTree.cc | 109 +++++++++++++++++++++++++++++++-----------------
ModelTree.hh | 29 ++++++-------
StaticModel.cc | 6 +--
StaticModel.hh | 2 +-
5 files changed, 89 insertions(+), 59 deletions(-)
diff --git a/DynamicModel.cc b/DynamicModel.cc
index 26b7c2c9..1bfb09d6 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -1995,7 +1995,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
- computePossiblySingularNormalization(contemporaneous_jacobian, cutoff == 0);
+ computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue);
diff --git a/ModelTree.cc b/ModelTree.cc
index 0b285315..d7d6f1e2 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -32,8 +32,8 @@
using namespace boost;
using namespace MFS;
-void
-ModelTree::computeNormalization(const set<pair<int, int> > &endo_eqs_incidence) throw (NormalizationException)
+bool
+ModelTree::computeNormalization(const jacob_map &contemporaneous_jacobian, bool verbose)
{
const int n = equation_number();
@@ -50,8 +50,8 @@ ModelTree::computeNormalization(const set<pair<int, int> > &endo_eqs_incidence)
// Fill in the graph
set<pair<int, int> > endo;
- for (set<pair<int, int> >::const_iterator it = endo_eqs_incidence.begin(); it != endo_eqs_incidence.end(); it++)
- add_edge(it->first + n, it->second, g);
+ for (jacob_map::const_iterator it = contemporaneous_jacobian.begin(); it != contemporaneous_jacobian.end(); it++)
+ add_edge(it->first.first + n, it->first.second, g);
// Compute maximum cardinality matching
vector<int> mate_map(2*n);
@@ -125,63 +125,96 @@ ModelTree::computeNormalization(const set<pair<int, int> > &endo_eqs_incidence)
// Check if all variables are normalized
vector<int>::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
if (it != mate_map.begin() + n)
- throw NormalizationException(symbol_table.getID(eEndogenous, it - mate_map.begin()));
+ {
+ if (verbose)
+ cerr << "ERROR: Could not normalize the model. Variable "
+ << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
+ << " is not in the maximum cardinality matching." << endl;
+ check = false;
+ }
+ return check;
}
void
-ModelTree::computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic)
+ModelTree::computeNonSingularNormalization(jacob_map &contemporaneous_jacobian, double cutoff, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian)
{
+ bool check = false;
+
cout << "Normalizing the model..." << endl;
- set<pair<int, int> > endo_eqs_incidence;
+ int n = equation_number();
- for (jacob_map::const_iterator it = contemporaneous_jacobian.begin();
- it != contemporaneous_jacobian.end(); it++)
- endo_eqs_incidence.insert(make_pair(it->first.first, it->first.second));
+ // compute the maximum value of each row of the contemporaneous Jacobian matrix
+ //jacob_map normalized_contemporaneous_jacobian;
+ jacob_map normalized_contemporaneous_jacobian(contemporaneous_jacobian);
+ vector<double> max_val(n, 0.0);
+ for (jacob_map::const_iterator iter = contemporaneous_jacobian.begin(); iter != contemporaneous_jacobian.end(); iter++)
+ if (fabs(iter->second) > max_val[iter->first.first])
+ max_val[iter->first.first] = fabs(iter->second);
- try
- {
- computeNormalization(endo_eqs_incidence);
- return;
- }
- catch (NormalizationException &e)
+ for (jacob_map::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
+ iter->second /= max_val[iter->first.first];
+
+ //We start with the highest value of the cutoff and try to normalize the model
+ double current_cutoff = 0.99999999;
+
+ int suppressed = 0;
+ while (!check && current_cutoff > 1e-19)
{
- if (try_symbolic)
- cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
- else
+ jacob_map tmp_normalized_contemporaneous_jacobian;
+ int suppress = 0;
+ for (jacob_map::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
+ if (fabs(iter->second) > max(current_cutoff, cutoff))
+ tmp_normalized_contemporaneous_jacobian[make_pair(iter->first.first, iter->first.second)] = iter->second;
+ else
+ suppress++;
+
+ if (suppress != suppressed)
+ check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
+ suppressed = suppress;
+ if (!check)
{
- cerr << "ERROR: Could not normalize the model. Variable "
- << symbol_table.getName(e.symb_id)
- << " is not in the maximum cardinality matching. Try to decrease the cutoff." << endl;
- exit(EXIT_FAILURE);
+ current_cutoff /= 2;
+ // In this last case try to normalize with the complete jacobian
+ if (current_cutoff <= 1e-19)
+ check = computeNormalization(normalized_contemporaneous_jacobian, false);
}
}
- // If no non-singular normalization can be found, try to find a normalization even with a potential singularity
- if (try_symbolic)
+ if (!check)
{
- endo_eqs_incidence.clear();
+ cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
+ //if no non-singular normalization can be found, try to find a normalization even with a potential singularity
+ jacob_map tmp_normalized_contemporaneous_jacobian;
set<pair<int, int> > endo;
- for (int i = 0; i < equation_number(); i++)
+ for (int i = 0; i < n; i++)
{
endo.clear();
equations[i]->collectEndogenous(endo);
for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
- endo_eqs_incidence.insert(make_pair(i, it->first));
+ tmp_normalized_contemporaneous_jacobian[make_pair(i, it->first)] = 1;
}
-
- try
+ check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
+ if (check)
{
- computeNormalization(endo_eqs_incidence);
- }
- catch (NormalizationException &e)
- {
- cerr << "ERROR: Could not normalize the model even with zero cutoff. Variable "
- << symbol_table.getName(e.symb_id)
- << " is not in the maximum cardinality matching." << endl;
- exit(EXIT_FAILURE);
+ // Update the jacobian matrix
+ for (jacob_map::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
+ {
+ if (static_jacobian.find(make_pair(it->first.first, it->first.second)) == static_jacobian.end())
+ static_jacobian[make_pair(it->first.first, it->first.second)] = 0;
+ if (dynamic_jacobian.find(make_pair(0, make_pair(it->first.first, it->first.second))) == dynamic_jacobian.end())
+ dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0;
+ if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end())
+ contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0;
+ }
}
}
+
+ if (!check)
+ {
+ cerr << "No normalization could be computed. Aborting." << endl;
+ exit(EXIT_FAILURE);
+ }
}
void
diff --git a/ModelTree.hh b/ModelTree.hh
index 7b1d7036..4614a12c 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -136,24 +136,21 @@ protected:
//! for each block contains pair< max_lag, max_lead>
t_lag_lead_vector block_lag_lead;
- //! Exception thrown when normalization fails
- class NormalizationException
- {
- public:
- //! A variable missing from the maximum cardinal matching
- int symb_id;
- NormalizationException(int symb_id_arg) : symb_id(symb_id_arg)
- {
- }
- };
-
//! Compute the matching between endogenous and variable using the jacobian contemporaneous_jacobian
- /*! \param endo_eqs_incidence A set indicating which endogenous appear in which equation. First element of pairs is equation number, second is type specific endo ID */
- void computeNormalization(const set<pair<int, int> > &endo_eqs_incidence) throw (NormalizationException);
+ /*!
+ \param contemporaneous_jacobian Jacobian used as an incidence matrix: all elements declared in the map (even if they are zero), are used as vertices of the incidence matrix
+ \return True if a complete normalization has been achieved
+ */
+ bool computeNormalization(const jacob_map &contemporaneous_jacobian, bool verbose);
+
//! Try to compute the matching between endogenous and variable using a decreasing cutoff
- /*! applied to the jacobian contemporaneous_jacobian and stop when a matching is found.*/
- /*! if no matching is found with a cutoff close to zero an error message is printout */
- void computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic);
+ /*!
+ Applied to the jacobian contemporaneous_jacobian and stop when a matching is found.
+ If no matching is found using a strictly positive cutoff, then a zero cutoff is applied (i.e. use a symbolic normalization); in that case, the method adds zeros in the jacobian matrices to reflect all the edges in the symbolic incidence matrix.
+ If no matching is found with a zero cutoff close to zero an error message is printout.
+ */
+ void computeNonSingularNormalization(jacob_map &contemporaneous_jacobian, double cutoff, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian);
+
//! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
//! Evaluate the jacobian and suppress all the elements below the cutoff
diff --git a/StaticModel.cc b/StaticModel.cc
index 784403e7..6bce4e2b 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -45,7 +45,7 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg,
}
void
-StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const
+StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_type &map_idx) const
{
first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id)));
if (it != first_derivatives.end())
@@ -539,7 +539,7 @@ StaticModel::writeModelEquationsCodeOrdered(const string file_name, const string
{
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
- compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
+ compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
{
FSTPG_ fstpg(0);
fstpg.write(code_file);
@@ -727,7 +727,7 @@ StaticModel::computingPass(const eval_context_type &eval_context, bool no_tmp_te
evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
- computePossiblySingularNormalization(contemporaneous_jacobian, cutoff == 0);
+ computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue);
diff --git a/StaticModel.hh b/StaticModel.hh
index 5b82ab54..1dbfb07d 100644
--- a/StaticModel.hh
+++ b/StaticModel.hh
@@ -74,7 +74,7 @@ private:
void computeTemporaryTermsOrdered();
//! Write derivative code of an equation w.r. to a variable
- void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const;
+ void compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_type &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_type &map_idx) const;
--
GitLab