From d186985fe06071bc70afbbbb16dbd8262fff4b6b Mon Sep 17 00:00:00 2001 From: Ted Turocy Date: Thu, 2 Jan 2025 13:05:34 +0000 Subject: [PATCH] More QuickSolver improvements --- src/solvers/enumpoly/efgpoly.cc | 4 +- src/solvers/enumpoly/gpartltr.imp | 25 ++++----- src/solvers/enumpoly/nfgpoly.cc | 4 +- src/solvers/enumpoly/quiksolv.cc | 2 +- src/solvers/enumpoly/quiksolv.h | 85 ++++++++++++------------------- src/solvers/enumpoly/quiksolv.imp | 82 +++++++++-------------------- 6 files changed, 76 insertions(+), 126 deletions(-) diff --git a/src/solvers/enumpoly/efgpoly.cc b/src/solvers/enumpoly/efgpoly.cc index ed603fd6b..25a692191 100644 --- a/src/solvers/enumpoly/efgpoly.cc +++ b/src/solvers/enumpoly/efgpoly.cc @@ -167,9 +167,9 @@ std::list> SolveSupport(const BehaviorSupportProfil bottoms = 0; tops = 1; - QuikSolv solver(equations); + QuickSolver solver(equations); try { - solver.FindCertainNumberOfRoots({bottoms, tops}, std::numeric_limits::max(), p_stopAfter); + solver.FindRoots({bottoms, tops}, p_stopAfter); } catch (const SingularMatrixException &) { p_isSingular = true; diff --git a/src/solvers/enumpoly/gpartltr.imp b/src/solvers/enumpoly/gpartltr.imp index a41784bba..3853b5f26 100644 --- a/src/solvers/enumpoly/gpartltr.imp +++ b/src/solvers/enumpoly/gpartltr.imp @@ -28,8 +28,6 @@ using namespace Gambit; // class: TreeOfPartials //--------------------------------------------------------------- - - /// Recursive generation of all partial derivatives of the root polynomial template void TreeOfPartials::TreeOfPartialsRECURSIVE(gTree> &t, gTreeNode> *n) const @@ -60,9 +58,10 @@ T TreeOfPartials::ValueOfPartialOfRootPoly(const int coord, const Vector & } template -T TreeOfPartials::MaximalNonconstantContributionRECURSIVE( - const gTreeNode> *n, const Vector &p, - const Vector &halvesoflengths, Vector &wrtos) const +T TreeOfPartials::MaximalNonconstantContributionRECURSIVE(const gTreeNode> *n, + const Vector &p, + const Vector &halvesoflengths, + Vector &wrtos) const { T answer = static_cast(0); @@ -109,7 +108,9 @@ template bool TreeOfPartials::PolyHasNoRootsIn(const Rectangle & template bool TreeOfPartials::MultiaffinePolyHasNoRootsIn(const Rectangle &r) const { - T sign = (this->RootNode()->GetData().Evaluate(r.Center()) > static_cast(0)) ? static_cast(1) : static_cast(-1); + T sign = (this->RootNode()->GetData().Evaluate(r.Center()) > static_cast(0)) + ? static_cast(1) + : static_cast(-1); Array zeros(Dmnsn()); std::fill(zeros.begin(), zeros.end(), 0); @@ -167,15 +168,15 @@ template bool TreeOfPartials::PolyEverywhereNegativeIn(const Rectan auto center = r.Center(); T constant = this->RootNode()->GetData().Evaluate(center); if (constant >= static_cast(0)) { - return false; + return false; } auto HalvesOfSideLengths = r.SideLengths() / 2; - return (this->MaximalNonconstantContribution(center, HalvesOfSideLengths) + constant < static_cast(0)); + return (this->MaximalNonconstantContribution(center, HalvesOfSideLengths) + constant < + static_cast(0)); } template -Matrix ListOfPartialTrees::DerivativeMatrix(const Vector &p, - const int NoEquations) const +Matrix ListOfPartialTrees::DerivativeMatrix(const Vector &p, const int NoEquations) const { Matrix answer(NoEquations, Dmnsn()); for (int i = 1; i <= NoEquations; i++) { @@ -186,8 +187,8 @@ Matrix ListOfPartialTrees::DerivativeMatrix(const Vector &p, return answer; } -template SquareMatrix -ListOfPartialTrees::SquareDerivativeMatrix(const Vector &p) const +template +SquareMatrix ListOfPartialTrees::SquareDerivativeMatrix(const Vector &p) const { SquareMatrix answer(Dmnsn()); for (int i = 1; i <= Dmnsn(); i++) { diff --git a/src/solvers/enumpoly/nfgpoly.cc b/src/solvers/enumpoly/nfgpoly.cc index 58ae0f165..f360777f7 100644 --- a/src/solvers/enumpoly/nfgpoly.cc +++ b/src/solvers/enumpoly/nfgpoly.cc @@ -116,10 +116,10 @@ EnumPolyStrategySupportSolve(const StrategySupportProfile &support, bool &is_sin Vector bottoms(Space.Dmnsn()), tops(Space.Dmnsn()); bottoms = 0; tops = 1; - QuikSolv solver(equations); + QuickSolver solver(equations); is_singular = false; try { - solver.FindCertainNumberOfRoots({bottoms, tops}, std::numeric_limits::max(), p_stopAfter); + solver.FindRoots({bottoms, tops}, p_stopAfter); } catch (const SingularMatrixException &) { is_singular = true; diff --git a/src/solvers/enumpoly/quiksolv.cc b/src/solvers/enumpoly/quiksolv.cc index a2df5fbfb..b72936694 100644 --- a/src/solvers/enumpoly/quiksolv.cc +++ b/src/solvers/enumpoly/quiksolv.cc @@ -22,4 +22,4 @@ #include "quiksolv.imp" -template class QuikSolv; +template class QuickSolver; diff --git a/src/solvers/enumpoly/quiksolv.h b/src/solvers/enumpoly/quiksolv.h index 5f72507df..85d3fe0b1 100644 --- a/src/solvers/enumpoly/quiksolv.h +++ b/src/solvers/enumpoly/quiksolv.h @@ -32,49 +32,34 @@ using namespace Gambit; -/* - The (optimistically named) class described in this file is a method -of finding the roots of a system of polynomials and inequalities, with -equal numbers of equations and unknowns, that lie inside a given -rectangle. The general idea is to first ask whether the Taylor's -series information at the center of the rectangle precludes the -existence of roots, and if it does not, whether Newton's method leads -to a root, and if it does, whether the Taylor's series information at -the root precludes the existence of another root. If the roots in the -rectangle are not resolved by these queries, the rectangle is -subdivided into 2^d subrectangles, and the process is repeated on -each. This continues until it has been shown that all roots have been -found, or a predetermined search depth is reached. The bound on depth -is necessary because the procedure will not terminate if there are -singular roots. -*/ - -/* - The main constructor for this takes a gPolyList. The list must -be at least as long as the dimension Dmnsn() of the space of the -system. The first Dmnsn() polynomials are interpreted as equations, -while remaining polynomials are interpreted as inequalities in the -sense that the polynomial is required to be nonnegative. -*/ - -/* - * The original implementation of this used a custom floating-point - * class as its parameter T, which implemented fuzzy comparisons. - * This has since been rewritten such that it uses regular floating - * point with explicit tolerances; this did introduce some subtle - * bugs originally and it is possible some still remain. - */ - -template class QuikSolv { +/// @brief Find the roots of a system of polynomials and inequalities, with +/// equal numbers of equations and unknowns, that lie inside a given +/// rectangle. +/// +/// The general idea is to first ask whether the Taylor's +/// series information at the center of the rectangle precludes the +/// existence of roots, and if it does not, whether Newton's method leads +/// to a root, and if it does, whether the Taylor's series information at +/// the root precludes the existence of another root. If the roots in the +/// rectangle are not resolved by these queries, the rectangle is +/// subdivided into 2^d subrectangles, and the process is repeated on +/// each. This continues until it has been shown that all roots have been +/// found, or a predetermined search depth is reached. The bound on depth +/// is necessary because the procedure will not terminate if there are +/// singular roots. + +/// The list of polynomials must +/// be at least as long as the dimension Dmnsn() of the space of the +/// system. The first Dmnsn() polynomials are interpreted as equations, +/// while remaining polynomials are interpreted as inequalities in the +/// sense that the polynomial is required to be non-negative. +template class QuickSolver { private: - gPolyList m_system; - gPolyList m_normalizedSystem; + gPolyList m_system, m_normalizedSystem; int NoEquations; int NoInequalities; ListOfPartialTrees TreesOfPartials; - bool HasBeenSolved; List> Roots; - bool isMultiaffine; RectArray Equation_i_uses_var_j; // Supporting routines for the constructors @@ -100,39 +85,35 @@ template class QuikSolv { const SquareMatrix &) const; // Combine the last two steps into a single query - bool NewtonRootIsOnlyInRct(const Rectangle &, Vector &) const; - // Recursive parts of recursive methods - - void FindRootsRecursion(List> *, const Rectangle &, const int &, - Array &, int &iterations, int depth, const int &, int *) const; + void FindRoots(List> &, const Rectangle &, const int &, Array &, + int &iterations, int depth, const int &) const; public: - explicit QuikSolv(const gPolyList &p_system) + explicit QuickSolver(const gPolyList &p_system) : m_system(p_system), m_normalizedSystem(p_system.AmbientSpace(), p_system.NormalizedList()), NoEquations(std::min(m_system.Dmnsn(), m_system.Length())), NoInequalities(std::max(m_system.Length() - m_system.Dmnsn(), 0)), - TreesOfPartials(m_normalizedSystem), HasBeenSolved(false), - isMultiaffine(m_system.IsMultiaffine()), Equation_i_uses_var_j(Eq_i_Uses_j()) + TreesOfPartials(m_normalizedSystem), Equation_i_uses_var_j(Eq_i_Uses_j()) { } - QuikSolv(const QuikSolv &) = delete; - ~QuikSolv() = default; + QuickSolver(const QuickSolver &) = delete; + ~QuickSolver() = default; - QuikSolv &operator=(const QuikSolv &) = delete; + QuickSolver &operator=(const QuickSolver &) = delete; // Information int Dmnsn() const { return m_system.Dmnsn(); } const List> &RootList() const { return Roots; } - bool IsMultiaffine() const { return isMultiaffine; } + bool IsMultiaffine() const { return m_system.IsMultiaffine(); } // Refines the accuracy of roots obtained from other algorithms Vector NewtonPolishOnce(const Vector &) const; Vector SlowNewtonPolishOnce(const Vector &) const; - // The grand calculation - returns true if successful - bool FindCertainNumberOfRoots(const Rectangle &, const int &, const int &); + // Find up to `max_roots` roots inside rectangle `r` + bool FindRoots(const Rectangle &r, int max_roots); }; #endif // QUIKSOLV_H diff --git a/src/solvers/enumpoly/quiksolv.imp b/src/solvers/enumpoly/quiksolv.imp index 5534bf27e..4ebc1cf33 100644 --- a/src/solvers/enumpoly/quiksolv.imp +++ b/src/solvers/enumpoly/quiksolv.imp @@ -32,12 +32,11 @@ using namespace Gambit; - //------------------------------------------------------- // Supporting Calculations for the Constructors //------------------------------------------------------- -template RectArray QuikSolv::Eq_i_Uses_j() const +template RectArray QuickSolver::Eq_i_Uses_j() const { RectArray answer(m_system.Length(), Dmnsn()); for (int i = 1; i <= m_system.Length(); i++) { @@ -58,8 +57,7 @@ template RectArray QuikSolv::Eq_i_Uses_j() const //--------------------------- template -bool QuikSolv::SystemHasNoRootsIn(const Rectangle &r, - Array &precedence) const +bool QuickSolver::SystemHasNoRootsIn(const Rectangle &r, Array &precedence) const { for (int i = 1; i <= m_system.Length(); i++) { @@ -125,8 +123,7 @@ static bool fuzzy_equals(const Vector &x, const Vector &y) } template -bool QuikSolv::NewtonRootInRectangle(const Rectangle &r, - Vector &point) const +bool QuickSolver::NewtonRootInRectangle(const Rectangle &r, Vector &point) const { Vector zero(Dmnsn()); zero = 0; @@ -203,9 +200,8 @@ bool QuikSolv::NewtonRootInRectangle(const Rectangle &r, //------------------------------------ template -double QuikSolv::MaxDistanceFromPointToVertexAfterTransformation( - const Rectangle &r, const Vector &p, - const SquareMatrix &M) const +double QuickSolver::MaxDistanceFromPointToVertexAfterTransformation( + const Rectangle &r, const Vector &p, const SquareMatrix &M) const { // A very early implementation of this method used a type gDouble which // implemented fuzzy comparisons. Adding the epsilon parameter here is @@ -242,8 +238,8 @@ double QuikSolv::MaxDistanceFromPointToVertexAfterTransformation( } template -bool QuikSolv::HasNoOtherRootsIn(const Rectangle &r, const Vector &p, - const SquareMatrix &M) const +bool QuickSolver::HasNoOtherRootsIn(const Rectangle &r, const Vector &p, + const SquareMatrix &M) const { // assert (NoEquations == System.Dmnsn()); gPolyList system1 = m_normalizedSystem.TranslateOfSystem(p); @@ -256,11 +252,9 @@ bool QuikSolv::HasNoOtherRootsIn(const Rectangle &r, const Vector -bool QuikSolv::NewtonRootIsOnlyInRct(const Rectangle &r, - Vector &point) const +bool QuickSolver::NewtonRootIsOnlyInRct(const Rectangle &r, Vector &point) const { if (NewtonRootInRectangle(r, point)) { SquareMatrix Df = TreesOfPartials.SquareDerivativeMatrix(point); @@ -269,12 +263,8 @@ bool QuikSolv::NewtonRootIsOnlyInRct(const Rectangle &r, return false; } -//------------------------------------------- -// Improve Accuracy of Root -//------------------------------------------- - template -Vector QuikSolv::NewtonPolishOnce(const Vector &point) const +Vector QuickSolver::NewtonPolishOnce(const Vector &point) const { Vector oldevals = TreesOfPartials.ValuesOfRootPolys(point, NoEquations); Matrix Df = TreesOfPartials.DerivativeMatrix(point, NoEquations); @@ -284,7 +274,7 @@ Vector QuikSolv::NewtonPolishOnce(const Vector &point) const } template -Vector QuikSolv::SlowNewtonPolishOnce(const Vector &point) const +Vector QuickSolver::SlowNewtonPolishOnce(const Vector &point) const { Vector oldevals = TreesOfPartials.ValuesOfRootPolys(point, NoEquations); Matrix Df = TreesOfPartials.DerivativeMatrix(point, NoEquations); @@ -300,50 +290,30 @@ Vector QuikSolv::SlowNewtonPolishOnce(const Vector &point) co } } -//------------------------------------------- -// The Central Calculation -//------------------------------------------- - -template -bool QuikSolv::FindCertainNumberOfRoots(const Rectangle &r, const int &max_iterations, - const int &max_no_roots) +template bool QuickSolver::FindRoots(const Rectangle &r, const int max_roots) { - auto *rootlistptr = new List>(); + const int max_iterations = 100000; + Roots = List>(); if (NoEquations == 0) { Vector answer(0); - rootlistptr->push_back(answer); - Roots = *rootlistptr; - HasBeenSolved = true; + Roots.push_back(answer); return true; } Array precedence(m_system.Length()); // Orders search for nonvanishing poly - for (int i = 1; i <= m_system.Length(); i++) { - precedence[i] = i; - } + std::iota(precedence.begin(), precedence.end(), 1); int iterations = 0; - - int *no_found = new int(0); - FindRootsRecursion(rootlistptr, r, max_iterations, precedence, iterations, 1, max_no_roots, - no_found); - - if (iterations < max_iterations) { - Roots = *rootlistptr; - HasBeenSolved = true; - return true; - } - - return false; + FindRoots(Roots, r, max_iterations, precedence, iterations, 1, max_roots); + return iterations < max_iterations; } template -void QuikSolv::FindRootsRecursion(List> *rootlistptr, - const Rectangle &r, const int &max_iterations, - Array &precedence, int &iterations, int depth, - const int &max_no_roots, int *roots_found) const +void QuickSolver::FindRoots(List> &rootlist, const Rectangle &r, + const int &max_iterations, Array &precedence, int &iterations, + int depth, const int &max_no_roots) const { // // TLT: In some cases, this recursive process apparently goes into an @@ -371,26 +341,24 @@ void QuikSolv::FindRootsRecursion(List> *rootlistptr, } bool already_found = false; - for (size_t i = 1; i <= rootlistptr->size(); i++) { - if (fuzzy_equals(point, (*rootlistptr)[i])) { + for (size_t i = 1; i <= rootlist.size(); i++) { + if (fuzzy_equals(point, rootlist[i])) { already_found = true; } } if (!already_found) { - rootlistptr->push_back(point); - (*roots_found)++; + rootlist.push_back(point); } return; } for (const auto &cell : r.Orthants()) { - if (max_no_roots == 0 || *roots_found < max_no_roots) { + if (max_no_roots == 0 || rootlist.size() < max_no_roots) { if (iterations >= max_iterations || depth == MAX_DEPTH) { return; } iterations++; - FindRootsRecursion(rootlistptr, cell, max_iterations, precedence, iterations, depth + 1, - max_no_roots, roots_found); + FindRoots(rootlist, cell, max_iterations, precedence, iterations, depth + 1, max_no_roots); } } }