From d9acf90f7bd6013e4b2644831235c3cd39ed3c0a Mon Sep 17 00:00:00 2001
From: Leonardo de Moura <leonardo@microsoft.com>
Date: Fri, 8 Jan 2016 15:59:45 -0800
Subject: [PATCH] refactor(util/lp): use Lean exception

---
 src/util/lp/lar_solver.h            |  11 +-
 src/util/lp/lp_core_solver_base.h   |   2 +-
 src/util/lp/lp_primal_core_solver.h |   7 +-
 src/util/lp/lp_primal_simplex.h     |   9 +-
 src/util/lp/lp_settings.h           |   8 +-
 src/util/lp/lp_solver.h             |  21 +-
 src/util/lp/numeric_pair.h          | 418 ++++++++++++++--------------
 src/util/lp/static_matrix.h         |  18 +-
 8 files changed, 235 insertions(+), 259 deletions(-)

diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h
index 25b43ca54..3f50a7996 100644
--- a/src/util/lp/lar_solver.h
+++ b/src/util/lp/lar_solver.h
@@ -261,7 +261,7 @@ class lar_solver {
         }
         break;
         default:
-            throw "unexpected";
+            lean_unreachable();
         }
     }
 
@@ -338,7 +338,7 @@ class lar_solver {
         case at_upper_bound: return upper_bound[j];
         case free_of_bounds: return zero_of_type<V>();
         default:
-            throw "unexpected type";
+            lean_unreachable();
         }
     }
 
@@ -421,10 +421,8 @@ public:
         case GE: return left_side_val >= constr.m_right_side;
         case GT: return left_side_val >= constr.m_right_side;
         case EQ: return left_side_val == constr.m_right_side;
-
         default:
-            throw "unexpected";
-            return false;
+            lean_unreachable();
         }
     }
 
@@ -833,8 +831,7 @@ public:
             return abs(left_side_val - norm_constr.m_right_side);
 
         default:
-            throw "unexpected";
-            return numeric_traits<mpq>::zero();
+            lean_unreachable();
         }
     }
 
diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h
index ed044931f..5090d6288 100644
--- a/src/util/lp/lp_core_solver_base.h
+++ b/src/util/lp/lp_core_solver_base.h
@@ -748,7 +748,7 @@ public:
         case upper_bound:
             return at_upper_bound;
         default:
-            throw "unexpected column type";
+            lean_unreachable();
         }
     }
 };
diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h
index 940a6758b..9a4c29ca5 100644
--- a/src/util/lp/lp_primal_core_solver.h
+++ b/src/util/lp/lp_primal_core_solver.h
@@ -1159,8 +1159,7 @@ public: // todo : move public lower ; it is for debug only
             slope_at_entering += delta;
             break;
         default:
-            throw 1;
-            lean_assert(false); // such a type does not exist
+            lean_unreachable();
         }
     }
 
@@ -1222,9 +1221,7 @@ public: // todo : move public lower ; it is for debug only
             return true;
             break;
         default:
-            lean_assert(false);
-            throw "cannot be here";
-            break;
+            lean_unreachable();
         }
     }
 
diff --git a/src/util/lp/lp_primal_simplex.h b/src/util/lp/lp_primal_simplex.h
index 1a5071cdb..fe16910c8 100644
--- a/src/util/lp/lp_primal_simplex.h
+++ b/src/util/lp/lp_primal_simplex.h
@@ -340,8 +340,7 @@ public:
         for (auto it : this->m_columns) {
             auto sol_it = solution.find(it.second->get_name());
             if (sol_it == solution.end()) {
-                std::cout << "cannot find column " << it.first << " in solution " << std::endl;
-                throw;
+                throw exception(sstream() << "cannot find column " << it.first << " in solutio");
             }
 
             if (!it.second->bounds_hold(sol_it->second)) {
@@ -356,8 +355,7 @@ public:
     T get_row_value(unsigned i, std::unordered_map<std::string, T> const & solution, bool print) {
         auto it = this->m_A_values.find(i);
         if (it == this->m_A_values.end()) {
-            std::cout << "cannot find row " << i << std::endl;
-            throw "get_row_value";
+            throw exception(sstream() << "cannot find row " << i);
         }
         T ret = numeric_traits<T>::zero();
         for (auto & pair : it->second) {
@@ -414,8 +412,7 @@ public:
             }
             return true;;
         }
-        std::cout << "throw in row_constraint_holds " << std::endl;
-        throw "wrong case";
+        lean_unreachable();
     }
 
     bool row_constraints_hold(std::unordered_map<std::string, T> const & solution) {
diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h
index 20cf3f940..3d6ca8b52 100644
--- a/src/util/lp/lp_settings.h
+++ b/src/util/lp/lp_settings.h
@@ -35,7 +35,7 @@ inline std::string column_type_to_string(column_type t) {
     case free_column:
         return std::string("free_column");
     default:
-        throw "unexpected";
+        lean_unreachable();
     }
 }
 
@@ -70,8 +70,7 @@ inline const char* lp_status_to_string(lp_status status) {
     case EMPTY: return "EMPTY";
     case UNSTABLE: return "UNSTABLE";
     default:
-        std::cout << "throwing in lp_status_to_string " << std::endl;
-        throw "lp_status_to_string";
+        lean_unreachable();
     }
 }
 
@@ -85,8 +84,7 @@ inline lp_status lp_status_from_string(std::string status) {
     if (status == "TIME_EXHAUSTED") return lp_status::TIME_EXHAUSTED;
     if (status == "ITERATIONS_EXHAUSTED") return lp_status::ITERATIONS_EXHAUSTED;
     if (status == "EMPTY") return lp_status::EMPTY;
-    std::cout << "unexpected status " << status << std::endl;
-    throw "wrong status in a test file";
+    lean_unreachable();
 }
 
 enum non_basic_column_value_position { at_low_bound, at_upper_bound, at_fixed, free_of_bounds };
diff --git a/src/util/lp/lp_solver.h b/src/util/lp/lp_solver.h
index 9bb4ac04f..7fe50160c 100644
--- a/src/util/lp/lp_solver.h
+++ b/src/util/lp/lp_solver.h
@@ -6,15 +6,16 @@
 */
 
 #pragma once
-
-#include "util/lp/lp_settings.h"
 #include <string>
 #include <unordered_map>
+#include <algorithm>
+#include <vector>
+#include "util/exception.h"
+#include "util/sstream.h"
+#include "util/lp/lp_settings.h"
 #include "util/lp/column_info.h"
 #include "util/lp/static_matrix.h"
-#include <algorithm>
 #include "util/lp/lp_core_solver_base.h"
-#include <vector>
 #include "util/lp/scaler.h"
 namespace lean {
 enum lp_relation  {
@@ -109,8 +110,7 @@ public:
     T get_column_value_by_name(std::string name) const {
         auto it = m_names_to_columns.find(name);
         if (it == m_names_to_columns.end()) {
-            std::cout << "throwing in get_column_value_by_name" << std::endl;
-            throw "get_column_value_by_name";
+            throw exception(sstream() << "get_column_value_by_name " << name);
         }
         return get_column_value(it -> second);
     }
@@ -435,11 +435,9 @@ protected:
         case lp_relation::Less_or_equal:
             return row_le_is_obsolete(row, row_index);
         }
-
-        throw "unexpected relation";
+        lean_unreachable();
     }
 
-
     void remove_fixed_or_zero_columns() {
         for (auto & i_row : m_A_values) {
             remove_fixed_or_zero_columns_from_row(i_row.first, i_row.second);
@@ -515,8 +513,7 @@ protected:
         for (auto & row : m_A_values) {
             for (auto & col : row.second) {
                 if (col.second == numeric_traits<T>::zero() || m_columns[col.first]->is_fixed()) {
-                    std::cout << "found fixed column " << std::endl;
-                    throw "map_external_columns_to_core_solver_columns";
+                    throw exception("found fixed column");
                 }
                 unsigned j = col.first;
                 auto j_place = m_external_columns_to_core_solver_columns.find(j);
@@ -615,7 +612,7 @@ protected:
 
         auto row = this->m_A_values.find(i);
         if (row == this->m_A_values.end()) {
-            throw "cannot find row ";
+            throw exception("cannot find row");
         }
         for (auto col : row->second) {
             ret += col.second * this->m_columns[col.first]->get_shift();
diff --git a/src/util/lp/numeric_pair.h b/src/util/lp/numeric_pair.h
index 9eb0ad5b4..2efd55116 100644
--- a/src/util/lp/numeric_pair.h
+++ b/src/util/lp/numeric_pair.h
@@ -4,7 +4,6 @@
 
   Author: Lev Nachmanson
 */
-
 #pragma once
 #include <string>
 #include <algorithm>
@@ -12,228 +11,227 @@
 #include "util/numerics/double.h"
 
 namespace lean {
-    template <typename T>
-    std::string T_to_string(const T & t); // forward definition
+template <typename T>
+std::string T_to_string(const T & t); // forward definition
 
-    template <typename X, typename Y>
-    struct convert_struct {
-        static X convert(const Y & y){ return X(y);}
-        static bool is_epsilon_small(const X & x,  const double & y) { return abs(get_double(x)) < y; }
-        static bool below_bound_numeric(const X &, const X &, const Y &) { throw "don't call";}
-        static bool above_bound_numeric(const X &, const X &, const Y &) { throw "don't call";}
-    };
+template <typename X, typename Y>
+struct convert_struct {
+    static X convert(const Y & y){ return X(y);}
+    static bool is_epsilon_small(const X & x,  const double & y) { return abs(get_double(x)) < y; }
+    static bool below_bound_numeric(const X &, const X &, const Y &) { lean_unreachable(); }
+    static bool above_bound_numeric(const X &, const X &, const Y &) { lean_unreachable(); }
+};
 
 
-    template <>
-    struct convert_struct<double, mpq> {
-        static double convert(const mpq & q) {return q.get_double();}
-    };
+template <>
+struct convert_struct<double, mpq> {
+    static double convert(const mpq & q) {return q.get_double();}
+};
 
-    template <typename T>
-    struct numeric_pair {
-        T x;
-        T y;
-        // empty constructor
-        numeric_pair() {}
-        // another constructor
+template <typename T>
+struct numeric_pair {
+    T x;
+    T y;
+    // empty constructor
+    numeric_pair() {}
+    // another constructor
 
-        numeric_pair(T xp, T yp) : x(xp), y(yp) {}
-        template <typename X, typename Y> numeric_pair(X xp, Y yp) : numeric_pair(convert_struct<T, X>::convert(xp), convert_struct<T, Y>::convert(yp)) {}
+    numeric_pair(T xp, T yp) : x(xp), y(yp) {}
+    template <typename X, typename Y> numeric_pair(X xp, Y yp) : numeric_pair(convert_struct<T, X>::convert(xp), convert_struct<T, Y>::convert(yp)) {}
 
-        bool operator<(const numeric_pair& a) const {
-            return x < a.x || (x == a.x && y < a.y);
-        }
-
-        bool operator>(const numeric_pair& a) const {
-            return x > a.x || (x == a.x && y > a.y);
-        }
-
-        bool operator==(const numeric_pair& a) const  {
-            return a.x == x &&  a.y == y;
-        }
-
-        bool operator!=(const numeric_pair& a) const  {
-            return !(*this == a);
-        }
-
-        bool operator<=(const numeric_pair& a) const  {
-            return *this < a || *this == a;
-        }
-
-        bool operator>=(const numeric_pair& a) const  {
-            return *this > a || a == *this;
-        }
-
-        numeric_pair operator*(const T & a) const {
-            return numeric_pair(a * x, a * y);
-        }
-
-        numeric_pair operator/(const T & a) const {
-            T a_as_T(a);
-            return numeric_pair(x / a_as_T, y / a_as_T);
-        }
-
-        numeric_pair operator/(const numeric_pair &) const {
-            throw "should not be called";
-        }
-
-
-        numeric_pair operator+(const numeric_pair & a) const  {
-            return numeric_pair(a.x + x, a.y + y);
-        }
-
-        numeric_pair operator*(const numeric_pair & /*a*/) const  {
-            throw "should not be called";
-        }
-
-
-        numeric_pair&  operator+=(const numeric_pair & a) {
-            x += a.x;
-            y += a.y;
-            return *this;
-        }
-
-        numeric_pair&  operator-=(const numeric_pair & a) {
-            x -= a.x;
-            y -= a.y;
-            return *this;
-        }
-
-        numeric_pair&  operator/=(const T & a) {
-            x /= a;
-            y /= a;
-            return *this;
-        }
-
-        numeric_pair&  operator*=(const T & a) {
-            x *= a;
-            y *= a;
-            return *this;
-        }
-
-        numeric_pair operator-(const numeric_pair & a) const {
-            return numeric_pair(x - a.x, y - a.y);
-        }
-
-        numeric_pair operator-() const {
-            return numeric_pair(-x, -y);
-        }
-
-        static bool precize() { return numeric_traits<T>::precize();}
-
-        std::string to_string() const { return std::string("(") + T_to_string(x) + ", "  + T_to_string(y) + ")"; }
-    };
-
-
-    template <typename T>
-    std::ostream& operator<<(std::ostream& os, numeric_pair<T> const & obj) {
-        os << obj.to_string();
-        return os;
+    bool operator<(const numeric_pair& a) const {
+        return x < a.x || (x == a.x && y < a.y);
     }
 
-    template <typename T, typename X>
-    numeric_pair<T> operator*(const X & a, const numeric_pair<T> & r) {
-        return numeric_pair<T>(a * r.x, a * r.y);
+    bool operator>(const numeric_pair& a) const {
+        return x > a.x || (x == a.x && y > a.y);
     }
 
-    template <typename T, typename X>
-    numeric_pair<T> operator*(const numeric_pair<T> & r, const X & a) {
-        return numeric_pair<T>(a * r.x, a * r.y);
+    bool operator==(const numeric_pair& a) const  {
+        return a.x == x &&  a.y == y;
+    }
+
+    bool operator!=(const numeric_pair& a) const  {
+        return !(*this == a);
+    }
+
+    bool operator<=(const numeric_pair& a) const  {
+        return *this < a || *this == a;
+    }
+
+    bool operator>=(const numeric_pair& a) const  {
+        return *this > a || a == *this;
+    }
+
+    numeric_pair operator*(const T & a) const {
+        return numeric_pair(a * x, a * y);
+    }
+
+    numeric_pair operator/(const T & a) const {
+        T a_as_T(a);
+        return numeric_pair(x / a_as_T, y / a_as_T);
+    }
+
+    numeric_pair operator/(const numeric_pair &) const {
+        lean_unreachable();
     }
 
 
+    numeric_pair operator+(const numeric_pair & a) const  {
+        return numeric_pair(a.x + x, a.y + y);
+    }
+
+    numeric_pair operator*(const numeric_pair & /*a*/) const  {
+        lean_unreachable();
+    }
+
+    numeric_pair&  operator+=(const numeric_pair & a) {
+        x += a.x;
+        y += a.y;
+        return *this;
+    }
+
+    numeric_pair&  operator-=(const numeric_pair & a) {
+        x -= a.x;
+        y -= a.y;
+        return *this;
+    }
+
+    numeric_pair&  operator/=(const T & a) {
+        x /= a;
+        y /= a;
+        return *this;
+    }
+
+    numeric_pair&  operator*=(const T & a) {
+        x *= a;
+        y *= a;
+        return *this;
+    }
+
+    numeric_pair operator-(const numeric_pair & a) const {
+        return numeric_pair(x - a.x, y - a.y);
+    }
+
+    numeric_pair operator-() const {
+        return numeric_pair(-x, -y);
+    }
+
+    static bool precize() { return numeric_traits<T>::precize();}
+
+    std::string to_string() const { return std::string("(") + T_to_string(x) + ", "  + T_to_string(y) + ")"; }
+};
 
 
-    // template <numeric_pair, typename T>  bool precise() { return numeric_traits<T>::precise();}
-    template <typename T> double get_double(const numeric_pair<T> & ) { throw "do not call"; }
-    template <typename T>
-    class numeric_traits<numeric_pair<T>> {
-    public:
-        static bool precise() { return numeric_traits<T>::precise();}
-        static numeric_pair<T> zero() { return numeric_pair<T>(numeric_traits<T>::zero(), numeric_traits<T>::zero()); }
-        static bool is_zero(const numeric_pair<T> & v) { return numeric_traits<T>::is_zero(v.x) && numeric_traits<T>::is_zero(v.y); }
-        static double get_double(const numeric_pair<T> & v){ return numeric_traits<T>::get_double(v.x); } // just return the double of the first coordinate
-        static double one() { throw "do not call"; }
-    };
-    template <>
-    struct convert_struct<double, numeric_pair<double>> {
-        static double convert(const numeric_pair<double> & q) {return q.x;}
-    };
-
-    template <typename X> bool is_epsilon_small(const X & v, const double& eps);   // forward definition { return convert_struct<X, double>::is_epsilon_small(v, eps);}
-
-    template <typename T>
-    struct convert_struct<numeric_pair<T>, double> {
-        static numeric_pair<T> convert(const double & q) {
-            return numeric_pair<T>(convert_struct<T, double>::convert(q), zero_of_type<T>());
-        }
-        static bool is_epsilon_small(const numeric_pair<T> & p, const double & eps) {
-            return convert_struct<T, double>::is_epsilon_small(p.x, eps) && convert_struct<T, double>::is_epsilon_small(p.y, eps);
-        }
-        static bool below_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
-            throw "do not call";
-        }
-        static bool above_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
-            throw "do not call";
-        }
-    };
-
-    template <>
-    struct convert_struct<numeric_pair<double>, double> {
-        static numeric_pair<double> convert(const double & q) {
-            return numeric_pair<double>(q, 0.0);
-        }
-        static bool is_epsilon_small(const numeric_pair<double> & p, const double & eps) {
-            return std::abs(p.x) < eps && std::abs(p.y) < eps;
-        }
-
-        static int compare_on_coord(const double & x, const double & bound, const double eps) {
-            lean_assert(eps > 0);
-            if (bound == 0) return (x < - eps)? -1: (x > eps? 1 : 0); // it is an important special case
-            double relative = (bound > 0)? - eps: eps;
-            return (x < bound * (1.0 + relative) - eps)? -1 : ((x > bound * (1.0 - relative) + eps)? 1 : 0);
-        }
-
-        static bool below_bound_numeric(const numeric_pair<double> & x, const numeric_pair<double> & bound, const double & eps) {
-            int r = compare_on_coord(x.x, bound.x, eps);
-            if (r == 1) return false;
-            if (r == -1) return true;
-            // the first coordinates are almost the same
-            lean_assert(r == 0);
-            return compare_on_coord(x.y, bound.y, eps) == -1;
-        }
-
-        static bool above_bound_numeric(const numeric_pair<double> & x, const numeric_pair<double> & bound, const double & eps) {
-            int r = compare_on_coord(x.x, bound.x, eps);
-            if (r == -1) return false;
-            if (r ==  1) return true;
-            // the first coordinates are almost the same
-            lean_assert(r == 0);
-            return compare_on_coord(x.y, bound.y, eps) == 1;
-        }
-    };
-
-    template <>
-    struct convert_struct<double, double> {
-        static bool is_epsilon_small(const double& x, const double & eps) {
-            return x < eps && x > -eps;
-        }
-        static double convert(const double & y){ return y;}
-        static bool below_bound_numeric(const double & x, const double & bound, const double & eps) {
-            lean_assert(eps > 0);
-            if (bound == 0) return x < - eps;
-            double relative = (bound > 0)? - eps: eps;
-            return x < bound * (1.0 + relative) - eps;
-        }
-        static bool above_bound_numeric(const double & x, const double & bound, const double & eps) {
-            lean_assert(eps > 0);
-            if (bound == 0) return x > eps;
-            double relative = (bound > 0)?  eps: - eps;
-            return x > bound * (1.0 + relative) + eps;
-        }
-    };
-
-    template <typename X> bool is_epsilon_small(const X & v, const double &eps) { return convert_struct<X, double>::is_epsilon_small(v, eps);}
-    template <typename X> bool below_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::below_bound_numeric(x, bound, eps);}
-    template <typename X> bool above_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::above_bound_numeric(x, bound, eps);}
+template <typename T>
+std::ostream& operator<<(std::ostream& os, numeric_pair<T> const & obj) {
+    os << obj.to_string();
+    return os;
+}
+
+template <typename T, typename X>
+numeric_pair<T> operator*(const X & a, const numeric_pair<T> & r) {
+    return numeric_pair<T>(a * r.x, a * r.y);
+}
+
+template <typename T, typename X>
+numeric_pair<T> operator*(const numeric_pair<T> & r, const X & a) {
+    return numeric_pair<T>(a * r.x, a * r.y);
+}
+
+
+
+
+// template <numeric_pair, typename T>  bool precise() { return numeric_traits<T>::precise();}
+template <typename T> double get_double(const numeric_pair<T> & ) { lean_unreachable(); }
+template <typename T>
+class numeric_traits<numeric_pair<T>> {
+  public:
+    static bool precise() { return numeric_traits<T>::precise();}
+    static numeric_pair<T> zero() { return numeric_pair<T>(numeric_traits<T>::zero(), numeric_traits<T>::zero()); }
+    static bool is_zero(const numeric_pair<T> & v) { return numeric_traits<T>::is_zero(v.x) && numeric_traits<T>::is_zero(v.y); }
+    static double get_double(const numeric_pair<T> & v){ return numeric_traits<T>::get_double(v.x); } // just return the double of the first coordinate
+    static double one() { lean_unreachable(); }
+};
+template <>
+struct convert_struct<double, numeric_pair<double>> {
+    static double convert(const numeric_pair<double> & q) {return q.x;}
+};
+
+template <typename X> bool is_epsilon_small(const X & v, const double& eps);   // forward definition { return convert_struct<X, double>::is_epsilon_small(v, eps);}
+
+template <typename T>
+struct convert_struct<numeric_pair<T>, double> {
+    static numeric_pair<T> convert(const double & q) {
+        return numeric_pair<T>(convert_struct<T, double>::convert(q), zero_of_type<T>());
+    }
+    static bool is_epsilon_small(const numeric_pair<T> & p, const double & eps) {
+        return convert_struct<T, double>::is_epsilon_small(p.x, eps) && convert_struct<T, double>::is_epsilon_small(p.y, eps);
+    }
+    static bool below_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
+        lean_unreachable();
+    }
+    static bool above_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
+        lean_unreachable();
+    }
+};
+
+template <>
+struct convert_struct<numeric_pair<double>, double> {
+    static numeric_pair<double> convert(const double & q) {
+        return numeric_pair<double>(q, 0.0);
+    }
+    static bool is_epsilon_small(const numeric_pair<double> & p, const double & eps) {
+        return std::abs(p.x) < eps && std::abs(p.y) < eps;
+    }
+
+    static int compare_on_coord(const double & x, const double & bound, const double eps) {
+        lean_assert(eps > 0);
+        if (bound == 0) return (x < - eps)? -1: (x > eps? 1 : 0); // it is an important special case
+        double relative = (bound > 0)? - eps: eps;
+        return (x < bound * (1.0 + relative) - eps)? -1 : ((x > bound * (1.0 - relative) + eps)? 1 : 0);
+    }
+
+    static bool below_bound_numeric(const numeric_pair<double> & x, const numeric_pair<double> & bound, const double & eps) {
+        int r = compare_on_coord(x.x, bound.x, eps);
+        if (r == 1) return false;
+        if (r == -1) return true;
+        // the first coordinates are almost the same
+        lean_assert(r == 0);
+        return compare_on_coord(x.y, bound.y, eps) == -1;
+    }
+
+    static bool above_bound_numeric(const numeric_pair<double> & x, const numeric_pair<double> & bound, const double & eps) {
+        int r = compare_on_coord(x.x, bound.x, eps);
+        if (r == -1) return false;
+        if (r ==  1) return true;
+        // the first coordinates are almost the same
+        lean_assert(r == 0);
+        return compare_on_coord(x.y, bound.y, eps) == 1;
+    }
+};
+
+template <>
+struct convert_struct<double, double> {
+    static bool is_epsilon_small(const double& x, const double & eps) {
+        return x < eps && x > -eps;
+    }
+    static double convert(const double & y){ return y;}
+    static bool below_bound_numeric(const double & x, const double & bound, const double & eps) {
+        lean_assert(eps > 0);
+        if (bound == 0) return x < - eps;
+        double relative = (bound > 0)? - eps: eps;
+        return x < bound * (1.0 + relative) - eps;
+    }
+    static bool above_bound_numeric(const double & x, const double & bound, const double & eps) {
+        lean_assert(eps > 0);
+        if (bound == 0) return x > eps;
+        double relative = (bound > 0)?  eps: - eps;
+        return x > bound * (1.0 + relative) + eps;
+    }
+};
+
+template <typename X> bool is_epsilon_small(const X & v, const double &eps) { return convert_struct<X, double>::is_epsilon_small(v, eps);}
+template <typename X> bool below_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::below_bound_numeric(x, bound, eps);}
+template <typename X> bool above_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::above_bound_numeric(x, bound, eps);}
 }
diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h
index f0fa7c95a..88e054e7e 100644
--- a/src/util/lp/static_matrix.h
+++ b/src/util/lp/static_matrix.h
@@ -7,6 +7,9 @@
 
 #pragma once
 #include <vector>
+#include <set>
+#include <unordered_map>
+#include <utility>
 #include "util/numerics/numeric_traits.h"
 #include "util/numerics/xnumeral.h"
 #include "util/numerics/mpq.h"
@@ -15,10 +18,7 @@
 #include "util/numerics/double.h"
 #include "util/numerics/float.h"
 #include "util/numerics/mpfp.h"
-#include <set>
 #include "util/lp/sparse_vector.h"
-#include <unordered_map>
-#include <utility>
 #include "util/lp/matrix_domain.h"
 #include "util/lp/indexed_vector.h"
 
@@ -48,9 +48,7 @@ private:
     T m_value;
 };
 
-
-
-    // each assignment for this matrix should be issued only once!!!
+// each assignment for this matrix should be issued only once!!!
 template <typename T, typename X>
 class static_matrix
 #ifdef LEAN_DEBUG
@@ -461,15 +459,9 @@ public:
     virtual void set_number_of_rows(unsigned /*m*/) { }
     virtual void set_number_of_columns(unsigned /*n*/) { }
 #endif
-    // const T & get_value_of_column_cell(column_cell<T> const & cc) const {
-    //     lean_assert(cc.m_i < m_rows.size());
-    //     lean_assert(cc.m_offset < m_rows[cc.m_i].size());
-    //     return m_rows[cc.m_i][cc.m_offset].get_val();
-    // }
 
     T get_max_val_in_row(unsigned i) const {
-        lean_assert(false);
-        throw 0;// not implemented
+        lean_unreachable();
     }
 
     T get_balance() const {