Merge pull request #99 from chrhansk/hotfix-solver-errors

Handle miscellaneous edge cases in MA57

docs/conf.py

@@ -14,7 +14,7 @@
 project = "pygradflow"
 copyright = "2023, Christoph Hansknecht"
 author = "Christoph Hansknecht"
-release = "0.5.6"
+release = "0.5.7"
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

pygradflow/display.py

@@ -191,7 +191,11 @@ def iter_cols(problem):
             AttrColumn("Bound inf", 16, "{:16.8e}", IterateAttr("bound_violation"))
         )
 
-    cols.append(AttrColumn("Cons inf", 16, "{:16.8e}", IterateAttr("cons_violation")))
+    if problem.num_cons > 0:
+        cols.append(
+            AttrColumn("Cons inf", 16, "{:16.8e}", IterateAttr("cons_violation"))
+        )
+
     cols.append(AttrColumn("Dual inf", 16, "{:16.8e}", IterateAttr("stat_res")))
 
     return cols

pygradflow/linear_solver/ma57_solver.py

@@ -1,5 +1,6 @@
 from typing import Optional
 
+import numpy as np
 from pyomo.contrib.pynumero.linalg.base import LinearSolverStatus
 from pyomo.contrib.pynumero.linalg.ma57_interface import MA57
 
@@ -17,19 +18,52 @@ def __init__(self, mat, symmetric=False, report_rcond=False):
         # if report_rcond:
         #     self.solver.set_icntl(10, 1)
 
-        status = self.solver.do_symbolic_factorization(self.mat)
+        size = self.mat.shape[0]
 
-        if status.status != LinearSolverStatus.successful:
-            raise LinearSolverError("Failed to compute symbolic factorization")
+        if size == 0:
+            self.solver = None
+            return
 
-        status = self.solver.do_numeric_factorization(self.mat)
+        def symbolic():
+            return self.solver.do_symbolic_factorization(self.mat, raise_on_error=False)
 
-        if status.status != LinearSolverStatus.successful:
-            raise LinearSolverError("Failed to compute numeric factorization")
+        self._try_fact(symbolic)
+
+        def numeric():
+            return self.solver.do_numeric_factorization(self.mat, raise_on_error=False)
+
+        self._try_fact(numeric)
+
+    def _try_fact(self, func):
+        num_tries = 10
+        for _ in range(num_tries):
+            status = func()
+
+            if status.status == LinearSolverStatus.not_enough_memory:
+                self.solver.increase_memory_allocation()
+                continue
+
+            self._handle_fact_status(status)
+            break
+        else:
+            raise Exception("Failed to compute symbolic factorization")
+
+    def _handle_fact_status(self, status):
+        if status.status == LinearSolverStatus.successful:
+            return
+
+        if status.status == LinearSolverStatus.singular:
+            raise LinearSolverError("Matrix is singular")
+
+        raise Exception("Failed to compute factorization: %s", status)
 
     def solve(self, rhs, trans=False, initial_sol=None):
         from pyomo.contrib.pynumero.linalg.base import LinearSolverStatus
 
+        if self.solver is None:
+            assert rhs.size == 0
+            return np.zeros_like(rhs)
+
         assert not trans
 
         x, status = self.solver.do_back_solve(rhs)
@@ -40,6 +74,8 @@ def solve(self, rhs, trans=False, initial_sol=None):
         return x
 
     def num_neg_eigvals(self):
+        if self.solver is None:
+            return 0
         return self.solver.get_info(24)
 
     def rcond(self) -> Optional[float]:

pygradflow/runners/runner.py

@@ -246,30 +246,36 @@ def create_csv_row(self, args, instance, result):
             "size": instance.size,
         }
 
+        default_props = {
+            "iterations": 0,
+            "num_accepted_steps": 0,
+            "final_scaled_obj": 0.0,
+            "final_stat_res": 0.0,
+            "final_cons_violation": 0.0,
+            "dist_factor": 0.0,
+        }
+
         if result == "timeout":
             return {
                 **info,
                 "status": "timeout",
                 "total_time": args.time_limit,
-                "iterations": 0,
-                "num_accepted_steps": 0,
+                **default_props,
             }
 
         elif result == "error":
-            return {
-                **info,
-                "status": "error",
-                "total_time": 0.0,
-                "iterations": 0,
-                "num_accepted_steps": 0,
-            }
+            return {**info, "status": "error", "total_time": 0.0, **default_props}
         else:
             return {
                 **info,
                 "status": SolverStatus.short_name(result.status),
                 "total_time": result.total_time,
                 "iterations": result.iterations,
                 "num_accepted_steps": result.num_accepted_steps,
+                "final_scaled_obj": result.final_scaled_obj,
+                "final_stat_res": result.final_stat_res,
+                "final_cons_violation": result.final_cons_violation,
+                "dist_factor": result.dist_factor,
             }
 
     def solve(self, instances, args):
@@ -292,6 +298,10 @@ def solve(self, instances, args):
             "total_time",
             "iterations",
             "num_accepted_steps",
+            "final_scaled_obj",
+            "final_stat_res",
+            "final_cons_violation",
+            "dist_factor",
         ]
 
         with open(filename, "w") as output_file:

pygradflow/scale.py

@@ -5,6 +5,7 @@
 
 from pygradflow.params import Params, ScalingType
 from pygradflow.problem import Problem
+from pygradflow.util import sparse_zero
 
 
 def scale_symmetric(A):
@@ -244,12 +245,16 @@ def create_scaling(
         raise ValueError("Custom scaling requires explicit scaling")
 
     if scaling_type == ScalingType.Nominal:
-        cons_val = problem.cons(x0)
+        if problem.num_cons > 0:
+            cons_val = problem.cons(x0)
+        else:
+            cons_val = np.array([], dtype=x0.dtype)
         return Scaling.from_nominal_values(x0, cons_val)
 
-    cons_jac = None
     if problem.num_cons > 0:
         cons_jac = problem.cons_jac(x0)
+    else:
+        cons_jac = sparse_zero(shape=(0, problem.num_vars))
 
     if scaling_type == ScalingType.GradJac:
         obj_grad = problem.obj_grad(x0)

pygradflow/solver.py

@@ -379,5 +379,8 @@ def solve(
             num_accepted_steps=accepted_steps,
             total_time=total_time,
             dist_factor=dist_factor,
+            final_scaled_obj=iterate.obj,
+            final_stat_res=iterate.stat_res,
+            final_cons_violation=iterate.cons_violation,
             **result_props,
         )

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pygradflow"
-version = "0.5.6"
+version = "0.5.7"
 description = "PyGradFlow is a simple implementation of the sequential homotopy method to be used to solve general nonlinear programs."
 authors = ["Christoph Hansknecht <[email protected]>"]
 readme = "README.md"