ModelTree::computeNormalization(): throw an exception in case of normalization failure

It would previously return a boolean. The exception is more convenient for
producing a different error message in the case of the specialized algorithm
introduced in the next commit.
master
Sébastien Villemot 2023-04-24 17:47:29 +02:00
parent e761da71bd
commit d246f9f99a
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
2 changed files with 41 additions and 36 deletions

View File

@ -234,8 +234,8 @@ ModelTree::operator=(const ModelTree &m)
return *this;
}
bool
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
void
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian)
{
const int n {static_cast<int>(equations.size())};
@ -261,19 +261,11 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
// Check if all variables are normalized
if (auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
it != mate_map.begin() + n)
{
if (verbose)
cerr << "Could not normalize the " << modelClassName() << ". Variable "
<< symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
<< " is not in the maximum cardinality matching." << endl;
return false;
}
throw ModelNormalizationFailed {symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin())) };
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
endo2eq.resize(equations.size());
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), [=](vertex_descriptor_t i) { return i-n; });
return true;
}
bool
@ -307,9 +299,8 @@ ModelTree::computeNonSingularNormalization(const eval_context_t &eval_context)
double current_cutoff = 0.99999999;
const double cutoff_lower_limit = 1e-19;
bool found_normalization = false;
int last_suppressed = 0;
while (!found_normalization && current_cutoff > cutoff_lower_limit)
while (current_cutoff > cutoff_lower_limit)
{
// Drop elements below cutoff from normalized contemporaneous Jacobian
jacob_map_t normalized_contemporaneous_jacobian_above_cutoff;
@ -321,35 +312,44 @@ ModelTree::computeNonSingularNormalization(const eval_context_t &eval_context)
suppressed++;
if (suppressed != last_suppressed)
found_normalization = computeNormalization(normalized_contemporaneous_jacobian_above_cutoff, false);
try
{
computeNormalization(normalized_contemporaneous_jacobian_above_cutoff);
return true;
}
catch (ModelNormalizationFailed &e)
{
}
last_suppressed = suppressed;
if (!found_normalization)
{
current_cutoff /= 2;
// In this last case try to normalize with the complete jacobian
if (current_cutoff <= cutoff_lower_limit)
found_normalization = computeNormalization(normalized_contemporaneous_jacobian, false);
}
current_cutoff /= 2;
}
if (!found_normalization)
// Try to normalize with the complete numerical jacobian
try
{
computeNormalization(normalized_contemporaneous_jacobian);
return true;
}
catch (ModelNormalizationFailed &e)
{
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. */
auto symbolic_jacobian = computeSymbolicJacobian(true);
found_normalization = computeNormalization(symbolic_jacobian, true);
}
#ifdef DEBUG
for (size_t i {0}; i < equations.size(); i++)
cout << "Variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, i))
<< " associated to equation " << endo2eq[i] << " (" << equation_tags.getTagValueByEqnAndKey(endo2eq[i], "name") << ")" << endl;
#endif
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. */
auto symbolic_jacobian { computeSymbolicJacobian(true) };
try
{
computeNormalization(symbolic_jacobian);
return true;
}
catch (ModelNormalizationFailed &e)
{
cerr << "Could not normalize the " << modelClassName() << ". Variable "
<< e.unmatched_endo << " is not in the maximum cardinality matching." << endl;
}
/* NB: If normalization failed, an explanatory message has been printed by the last call
to computeNormalization(), which has verbose=true */
return found_normalization;
return false;
}
ModelTree::jacob_map_t

View File

@ -436,12 +436,17 @@ private:
// Compute {var,eq}_idx_orig2block from {var,eq}_idx_block2orig
void updateReverseVariableEquationOrderings();
struct ModelNormalizationFailed
{
string unmatched_endo; // Name of endogenous not in maximum cardinality matching
};
//! Compute the matching between endogenous and variable using the jacobian contemporaneous_jacobian
/*!
\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_t &contemporaneous_jacobian, bool verbose);
void computeNormalization(const jacob_map_t &contemporaneous_jacobian);
//! Try to compute the matching between endogenous and variable using a decreasing cutoff
/*!