More QuickSolver improvements

src/solvers/enumpoly/efgpoly.cc

@@ -167,9 +167,9 @@ std::list<MixedBehaviorProfile<double>> SolveSupport(const BehaviorSupportProfil
   bottoms = 0;
   tops = 1;
 
-  QuikSolv<double> solver(equations);
+  QuickSolver<double> solver(equations);
   try {
-    solver.FindCertainNumberOfRoots({bottoms, tops}, std::numeric_limits<int>::max(), p_stopAfter);
+    solver.FindRoots({bottoms, tops}, p_stopAfter);
   }
   catch (const SingularMatrixException &) {
     p_isSingular = true;

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 <class T>
 void TreeOfPartials<T>::TreeOfPartialsRECURSIVE(gTree<gPoly<T>> &t, gTreeNode<gPoly<T>> *n) const
@@ -60,9 +58,10 @@ T TreeOfPartials<T>::ValueOfPartialOfRootPoly(const int coord, const Vector<T> &
 }
 
 template <class T>
-T TreeOfPartials<T>::MaximalNonconstantContributionRECURSIVE(
-    const gTreeNode<gPoly<T>> *n, const Vector<T> &p,
-    const Vector<T> &halvesoflengths, Vector<int> &wrtos) const
+T TreeOfPartials<T>::MaximalNonconstantContributionRECURSIVE(const gTreeNode<gPoly<T>> *n,
+                                                             const Vector<T> &p,
+                                                             const Vector<T> &halvesoflengths,
+                                                             Vector<int> &wrtos) const
 {
   T answer = static_cast<T>(0);
 
@@ -109,7 +108,9 @@ template <class T> bool TreeOfPartials<T>::PolyHasNoRootsIn(const Rectangle<T> &
 
 template <class T> bool TreeOfPartials<T>::MultiaffinePolyHasNoRootsIn(const Rectangle<T> &r) const
 {
-  T sign = (this->RootNode()->GetData().Evaluate(r.Center()) > static_cast<T>(0)) ? static_cast<T>(1) : static_cast<T>(-1);
+  T sign = (this->RootNode()->GetData().Evaluate(r.Center()) > static_cast<T>(0))
+               ? static_cast<T>(1)
+               : static_cast<T>(-1);
 
   Array<int> zeros(Dmnsn());
   std::fill(zeros.begin(), zeros.end(), 0);
@@ -167,15 +168,15 @@ template <class T> bool TreeOfPartials<T>::PolyEverywhereNegativeIn(const Rectan
   auto center = r.Center();
   T constant = this->RootNode()->GetData().Evaluate(center);
   if (constant >= static_cast<T>(0)) {
-      return false;
+    return false;
   }
   auto HalvesOfSideLengths = r.SideLengths() / 2;
-  return (this->MaximalNonconstantContribution(center, HalvesOfSideLengths) + constant < static_cast<T>(0));
+  return (this->MaximalNonconstantContribution(center, HalvesOfSideLengths) + constant <
+          static_cast<T>(0));
 }
 
 template <class T>
-Matrix<T> ListOfPartialTrees<T>::DerivativeMatrix(const Vector<T> &p,
-                                                  const int NoEquations) const
+Matrix<T> ListOfPartialTrees<T>::DerivativeMatrix(const Vector<T> &p, const int NoEquations) const
 {
   Matrix<T> answer(NoEquations, Dmnsn());
   for (int i = 1; i <= NoEquations; i++) {
@@ -186,8 +187,8 @@ Matrix<T> ListOfPartialTrees<T>::DerivativeMatrix(const Vector<T> &p,
   return answer;
 }
 
-template <class T> SquareMatrix<T>
-ListOfPartialTrees<T>::SquareDerivativeMatrix(const Vector<T> &p) const
+template <class T>
+SquareMatrix<T> ListOfPartialTrees<T>::SquareDerivativeMatrix(const Vector<T> &p) const
 {
   SquareMatrix<T> answer(Dmnsn());
   for (int i = 1; i <= Dmnsn(); i++) {

src/solvers/enumpoly/nfgpoly.cc

@@ -116,10 +116,10 @@ EnumPolyStrategySupportSolve(const StrategySupportProfile &support, bool &is_sin
   Vector<double> bottoms(Space.Dmnsn()), tops(Space.Dmnsn());
   bottoms = 0;
   tops = 1;
-  QuikSolv<double> solver(equations);
+  QuickSolver<double> solver(equations);
   is_singular = false;
   try {
-    solver.FindCertainNumberOfRoots({bottoms, tops}, std::numeric_limits<int>::max(), p_stopAfter);
+    solver.FindRoots({bottoms, tops}, p_stopAfter);
   }
   catch (const SingularMatrixException &) {
     is_singular = true;

src/solvers/enumpoly/quiksolv.cc

@@ -22,4 +22,4 @@
 
 #include "quiksolv.imp"
 
-template class QuikSolv<double>;
+template class QuickSolver<double>;

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<T>.  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 T> 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 T> class QuickSolver {
 private:
-  gPolyList<T> m_system;
-  gPolyList<double> m_normalizedSystem;
+  gPolyList<T> m_system, m_normalizedSystem;
   int NoEquations;
   int NoInequalities;
   ListOfPartialTrees<double> TreesOfPartials;
-  bool HasBeenSolved;
   List<Vector<double>> Roots;
-  bool isMultiaffine;
   RectArray<bool> Equation_i_uses_var_j;
 
   // Supporting routines for the constructors
@@ -100,39 +85,35 @@ template <class T> class QuikSolv {
                          const SquareMatrix<double> &) const;
 
   // Combine the last two steps into a single query
-
   bool NewtonRootIsOnlyInRct(const Rectangle<double> &, Vector<double> &) const;
 
-  // Recursive parts of recursive methods
-
-  void FindRootsRecursion(List<Vector<double>> *, const Rectangle<double> &, const int &,
-                          Array<int> &, int &iterations, int depth, const int &, int *) const;
+  void FindRoots(List<Vector<double>> &, const Rectangle<double> &, const int &, Array<int> &,
+                 int &iterations, int depth, const int &) const;
 
 public:
-  explicit QuikSolv(const gPolyList<T> &p_system)
+  explicit QuickSolver(const gPolyList<T> &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<T> &) = delete;
-  ~QuikSolv() = default;
+  QuickSolver(const QuickSolver<T> &) = delete;
+  ~QuickSolver() = default;
 
-  QuikSolv<T> &operator=(const QuikSolv<T> &) = delete;
+  QuickSolver<T> &operator=(const QuickSolver<T> &) = delete;
 
   // Information
   int Dmnsn() const { return m_system.Dmnsn(); }
   const List<Vector<double>> &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<double> NewtonPolishOnce(const Vector<double> &) const;
   Vector<double> SlowNewtonPolishOnce(const Vector<double> &) const;
 
-  // The grand calculation - returns true if successful
-  bool FindCertainNumberOfRoots(const Rectangle<T> &, const int &, const int &);
+  // Find up to `max_roots` roots inside rectangle `r`
+  bool FindRoots(const Rectangle<T> &r, int max_roots);
 };
 
 #endif // QUIKSOLV_H