From 07269af6b2c42e336ba9573412bf5d7d3d643911 Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Tue, 12 Mar 2024 21:12:31 +0000
Subject: [PATCH 01/10] QUASS: added steepest edge pricing

---
 src/qpsolver/quass.cpp               |  3 +
 src/qpsolver/settings.hpp            |  2 +-
 src/qpsolver/steepestedgepricing.hpp | 99 ++++++++++++++++++++++++----
 3 files changed, 92 insertions(+), 12 deletions(-)

diff --git a/src/qpsolver/quass.cpp b/src/qpsolver/quass.cpp
index 91df50f965..e2621a3acb 100644
--- a/src/qpsolver/quass.cpp
+++ b/src/qpsolver/quass.cpp
@@ -171,6 +171,9 @@ static QpSolverStatus reduce(Runtime& rt, Basis& basis, const HighsInt newactive
 static std::unique_ptr<Pricing> getPricing(Runtime& runtime, Basis& basis,
                                     ReducedCosts& redcosts) {
   switch (runtime.settings.pricing) {
+    case PricingStrategy::SteepestEdge:
+      return std::unique_ptr<Pricing>(
+          new SteepestEdgePricing(runtime, basis, redcosts));
     case PricingStrategy::Devex:
       return std::unique_ptr<Pricing>(
           new DevexPricing(runtime, basis, redcosts));
diff --git a/src/qpsolver/settings.hpp b/src/qpsolver/settings.hpp
index 9108badf24..b31707f99e 100644
--- a/src/qpsolver/settings.hpp
+++ b/src/qpsolver/settings.hpp
@@ -6,7 +6,7 @@
 
 enum class RatiotestStrategy { TwoPass, Textbook };
 
-enum class PricingStrategy { DantzigWolfe, Devex };
+enum class PricingStrategy { SteepestEdge, DantzigWolfe, Devex };
 
 enum class OutputLevel { LIGHT, MEDIUM, HEAVY };
 
diff --git a/src/qpsolver/steepestedgepricing.hpp b/src/qpsolver/steepestedgepricing.hpp
index 4c30c692db..84a3064263 100644
--- a/src/qpsolver/steepestedgepricing.hpp
+++ b/src/qpsolver/steepestedgepricing.hpp
@@ -56,31 +56,108 @@ class SteepestEdgePricing : public Pricing {
       : runtime(rt),
         basis(bas),
         redcosts(rc),
-        weights(std::vector<double>(rt.instance.num_var, 1.0)) {};
+        weights(std::vector<double>(rt.instance.num_var, 1.0)) {
+          compute_exact_weights();
+        };
 
   HighsInt price(const Vector& x, const Vector& gradient) {
     HighsInt minidx = chooseconstrainttodrop(redcosts.getReducedCosts());
     return minidx;
   }
 
+
+  double compute_exact_weight(HighsInt i) {
+    Vector y_i = basis.btran(Vector::unit(runtime.instance.num_var, i));
+    return y_i.dot(y_i);
+  }
+
+  void compute_exact_weights() {
+    for (int i=0; i<runtime.instance.num_var; i++) {
+      weights[i] = compute_exact_weight(i);
+    }
+  }
+
+  bool check_weight(HighsInt i) {
+    double weight_is = weights[i];
+    double weight_comp = compute_exact_weight(i);
+    if(fabs(weight_comp - weight_is) > 10E-3) {
+      //printf("weights[%d] = %lf, should be %lf\n", i, weight_is, weight_comp);
+      return false;
+    }
+    return true;
+  }
+
+  bool check_weights() {
+
+    std::vector<int> correct_weights;
+    std::vector<int> incorrect_weights;
+    bool ret = true;
+    for (int i=0; i<runtime.instance.num_var; i++) {
+      bool correct = check_weight(i);
+        if (correct) {
+          correct_weights.push_back(i);
+        } else {
+          incorrect_weights.push_back(i);
+        }
+    }
+
+    printf("correct weights: ");
+    for (int& i : correct_weights) {
+      printf("%d ", i);
+    }
+    printf("\n");
+
+    printf("incorrect weights: ");
+    for (int& i : incorrect_weights) {
+      printf("%d ", i);
+    }
+    printf("\n");
+
+    return incorrect_weights.size() == 0;
+    }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     HighsInt rowindex_p = basis.getindexinfactor()[p];
+    //printf("Update weights, p = %d, rowindex = %d, q = %d\n", p, rowindex_p, q);
+
+    //if (!check_weights()) {
+    //  printf("weight check failed\n");
+    //  exit(1);
+    //}
+
+    Vector delta = basis.ftran(aq);
+
+    double old_weight_p_updated = weights[rowindex_p];
+    // exact weight coming in needs to come in before update.
+    double old_weight_p_computed = ep.dot(ep);
+
+    //if (fabs(old_weight_p_computed - old_weight_p_updated) >= 10E-3) {
+    //  printf("old weight[p] discrepancy: updated = %lf, computed=%lf\n", old_weight_p_updated, old_weight_p_computed);
+    //}
 
-    Vector v = basis.btran(aq);
+    double weight_p = old_weight_p_computed;
 
-    double weight_p = weights[rowindex_p];
+    
+    double t_p = aq.value[rowindex_p];
     for (HighsInt i = 0; i < runtime.instance.num_var; i++) {
-      if (i == rowindex_p) {
-        weights[i] = weight_p / (aq.value[rowindex_p] * aq.value[rowindex_p]);
-      } else {
-        weights[i] = weights[i] -
-                     2 * (aq.value[i] / aq.value[rowindex_p]) * (v.value[i]) +
-                     (aq.value[i] * aq.value[i]) /
-                         (aq.value[rowindex_p] * aq.value[rowindex_p]) *
-                         weight_p;
+      if (i != rowindex_p) {
+        double t_i = aq.value[i];
+        weights[i] = weights[i]
+                    - 2 * (t_i / t_p) * delta.value[i]
+                    + ((t_i * t_i) / (t_p * t_p)) * weight_p;
+        //printf("weights[%d] = %lf\n", i, weights[i]);
       }
     }
+    //Vector new_ep = basis.btran(Vector::unit(runtime.instance.num_var, rowindex_p)); 
+    //double computed_weight = new_ep.dot(new_ep);
+    double new_weight_p_updated = weight_p / (t_p * t_p);
+    
+    //if (fabs(updated_weight - computed_weight) > 10E-5) {
+    //  printf("updated weight %lf vs computed weight %lf. aq[p] = %lf\n", updated_weight, computed_weight, t_p);
+    //  printf("old weight = %lf, aq[p] = %lf, ^2 = %lf, new weight = %lf\n", weight_p, t_p, t_p*t_p, updated_weight);
+    //}
+    weights[rowindex_p] = new_weight_p_updated;
   }
 };
 

From 5b71f2e2857d676d605d7c96fe7f17cb86b6b0e3 Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Tue, 12 Mar 2024 21:51:22 +0000
Subject: [PATCH 02/10] QUASS: use simplex_primal_edge_weight_strategy in QP
 solver

---
 src/lp_data/Highs.cpp | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
index aab663f179..0f402e61bb 100644
--- a/src/lp_data/Highs.cpp
+++ b/src/lp_data/Highs.cpp
@@ -3406,6 +3406,20 @@ HighsStatus Highs::callSolveQp() {
   settings.iterationlimit = options_.simplex_iteration_limit;
   settings.lambda_zero_threshold = options_.dual_feasibility_tolerance;
 
+  switch(options_.simplex_primal_edge_weight_strategy) {
+    case 0:
+      settings.pricing = PricingStrategy::DantzigWolfe;
+      break;
+    case 1:
+      settings.pricing = PricingStrategy::Devex;
+      break;
+    case 2:
+      settings.pricing = PricingStrategy::SteepestEdge;
+      break;
+    default:
+      settings.pricing = PricingStrategy::Devex;
+  }
+
   // print header for QP solver output
   highsLogUser(options_.log_options, HighsLogType::kInfo,
                "Iteration, Runtime, ObjVal, NullspaceDim\n");

From ea6fb397c95e85b97dec3d53fcaf1db6a8924083 Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Tue, 12 Mar 2024 23:06:45 +0000
Subject: [PATCH 03/10] QUASS: reinvert in basis now triggers reinversion of
 all dependent structures (Cholesky factorization of reduced Hessian, reduced
 gradient, gradient, reduced costs, steepest edge weights)

---
 src/qpsolver/basis.cpp               |  3 +-
 src/qpsolver/basis.hpp               |  8 ++++-
 src/qpsolver/dantzigpricing.hpp      |  4 +++
 src/qpsolver/devexharrispricing.hpp  |  4 +++
 src/qpsolver/devexpricing.hpp        |  4 +++
 src/qpsolver/factor.hpp              | 53 ++++++++++++++--------------
 src/qpsolver/gradient.hpp            |  8 ++---
 src/qpsolver/pricing.hpp             |  1 +
 src/qpsolver/quass.cpp               | 20 +++++++++++
 src/qpsolver/reducedcosts.hpp        | 10 +++---
 src/qpsolver/reducedgradient.hpp     | 12 +++----
 src/qpsolver/settings.hpp            |  2 +-
 src/qpsolver/steepestedgepricing.hpp |  4 +++
 13 files changed, 89 insertions(+), 44 deletions(-)

diff --git a/src/qpsolver/basis.cpp b/src/qpsolver/basis.cpp
index 987a5f1caa..24541c2405 100644
--- a/src/qpsolver/basis.cpp
+++ b/src/qpsolver/basis.cpp
@@ -88,6 +88,7 @@ void Basis::rebuild() {
        i < activeconstraintidx.size() + nonactiveconstraintsidx.size(); i++) {
     constraintindexinbasisfactor[baseindex[i]] = i;
   }
+  reinversion_hint = false;
 }
 
 void Basis::report() {
@@ -166,7 +167,7 @@ void Basis::updatebasis(const Settings& settings, HighsInt newactivecon, HighsIn
 
   updatessinceinvert++;
   if (updatessinceinvert >= settings.reinvertfrequency || hint != 99999) {
-    rebuild();
+    reinversion_hint = true;
   }
   // since basis changed, buffered values are no longer valid
   buffered_p = -1;
diff --git a/src/qpsolver/basis.hpp b/src/qpsolver/basis.hpp
index aa40ed3dd1..96411bd0dd 100644
--- a/src/qpsolver/basis.hpp
+++ b/src/qpsolver/basis.hpp
@@ -73,7 +73,6 @@ class Basis {
   std::vector<HighsInt> constraintindexinbasisfactor;
 
   void build();
-  void rebuild();
 
   // buffer to avoid recreating vectors
   Vector buffer_column_aq;
@@ -85,10 +84,17 @@ class Basis {
   HVector row_ep;
   HVector col_aq;
 
+  bool reinversion_hint = false;
  public:
+
+
   Basis(Runtime& rt, std::vector<HighsInt> active,
         std::vector<BasisStatus> atlower, std::vector<HighsInt> inactive);
 
+  bool getreinversionhint() { return reinversion_hint; }
+
+  void rebuild();
+
   HighsInt getnupdatessinceinvert() { return updatessinceinvert; }
 
   HighsInt getnumactive() const { return activeconstraintidx.size(); };
diff --git a/src/qpsolver/dantzigpricing.hpp b/src/qpsolver/dantzigpricing.hpp
index e1f601a387..f043437cea 100644
--- a/src/qpsolver/dantzigpricing.hpp
+++ b/src/qpsolver/dantzigpricing.hpp
@@ -59,6 +59,10 @@ class DantzigPricing : public Pricing {
     return minidx;
   }
 
+  void recompute() {
+    // do nothing
+  }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     // does nothing
diff --git a/src/qpsolver/devexharrispricing.hpp b/src/qpsolver/devexharrispricing.hpp
index 3ee854106b..3a9b64126b 100644
--- a/src/qpsolver/devexharrispricing.hpp
+++ b/src/qpsolver/devexharrispricing.hpp
@@ -63,6 +63,10 @@ class DevexHarrisPricing : public Pricing {
     return minidx;
   }
 
+  void recompute() {
+    // do nothing
+  }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     HighsInt rowindex_p = basis.getindexinfactor()[p];
diff --git a/src/qpsolver/devexpricing.hpp b/src/qpsolver/devexpricing.hpp
index 744ca3f726..d3ba6a03dd 100644
--- a/src/qpsolver/devexpricing.hpp
+++ b/src/qpsolver/devexpricing.hpp
@@ -73,6 +73,10 @@ class DevexPricing : public Pricing {
     return minidx;
   }
 
+  void recompute() {
+    // do nothing
+  }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     HighsInt rowindex_p = basis.getindexinfactor()[p];
diff --git a/src/qpsolver/factor.hpp b/src/qpsolver/factor.hpp
index 6892aeeeb0..40ef603c46 100644
--- a/src/qpsolver/factor.hpp
+++ b/src/qpsolver/factor.hpp
@@ -26,6 +26,33 @@ class CholeskyFactor {
   bool has_negative_eigenvalue = false;
   std::vector<double> a;
 
+  void resize(HighsInt new_k_max) {
+    std::vector<double> L_old = L;
+    L.clear();
+    L.resize((new_k_max) * (new_k_max));
+    const HighsInt l_size = L.size();
+    // Driven by #958, changes made in following lines to avoid array
+    // bound error when new_k_max < current_k_max
+    HighsInt min_k_max = min(new_k_max, current_k_max);
+    for (HighsInt i = 0; i < min_k_max; i++) {
+      for (HighsInt j = 0; j < min_k_max; j++) {
+	assert(i * (new_k_max) + j < l_size);
+        L[i * (new_k_max) + j] = L_old[i * current_k_max + j];
+      }
+    }
+    current_k_max = new_k_max;
+  }
+
+ public:
+  CholeskyFactor(Runtime& rt, Basis& bas) : runtime(rt), basis(bas) {
+    uptodate = false;
+    current_k_max =
+        max(min((HighsInt)ceil(rt.instance.num_var / 16.0), (HighsInt)1000),
+            basis.getnuminactive());
+    L.resize(current_k_max * current_k_max);
+  }
+
+
   void recompute() {
     std::vector<std::vector<double>> orig;
     HighsInt dim_ns = basis.getinactive().size();
@@ -70,32 +97,6 @@ class CholeskyFactor {
     uptodate = true;
   }
 
-  void resize(HighsInt new_k_max) {
-    std::vector<double> L_old = L;
-    L.clear();
-    L.resize((new_k_max) * (new_k_max));
-    const HighsInt l_size = L.size();
-    // Driven by #958, changes made in following lines to avoid array
-    // bound error when new_k_max < current_k_max
-    HighsInt min_k_max = min(new_k_max, current_k_max);
-    for (HighsInt i = 0; i < min_k_max; i++) {
-      for (HighsInt j = 0; j < min_k_max; j++) {
-	assert(i * (new_k_max) + j < l_size);
-        L[i * (new_k_max) + j] = L_old[i * current_k_max + j];
-      }
-    }
-    current_k_max = new_k_max;
-  }
-
- public:
-  CholeskyFactor(Runtime& rt, Basis& bas) : runtime(rt), basis(bas) {
-    uptodate = false;
-    current_k_max =
-        max(min((HighsInt)ceil(rt.instance.num_var / 16.0), (HighsInt)1000),
-            basis.getnuminactive());
-    L.resize(current_k_max * current_k_max);
-  }
-
   QpSolverStatus expand(const Vector& yp, Vector& gyp, Vector& l, Vector& m) {
     if (!uptodate) {
       return QpSolverStatus::OK;
diff --git a/src/qpsolver/gradient.hpp b/src/qpsolver/gradient.hpp
index fea6a3616e..839c93ced3 100644
--- a/src/qpsolver/gradient.hpp
+++ b/src/qpsolver/gradient.hpp
@@ -11,6 +11,10 @@ class Gradient {
   bool uptodate;
   HighsInt numupdates = 0;
 
+ public:
+  Gradient(Runtime& rt)
+      : runtime(rt), gradient(Vector(rt.instance.num_var)), uptodate(false) {}
+
   void recompute() {
     runtime.instance.Q.vec_mat(runtime.primal, gradient);
     gradient += runtime.instance.c;
@@ -18,10 +22,6 @@ class Gradient {
     numupdates = 0;
   }
 
- public:
-  Gradient(Runtime& rt)
-      : runtime(rt), gradient(Vector(rt.instance.num_var)), uptodate(false) {}
-
   Vector& getGradient() {
     if (!uptodate ||
         numupdates >= runtime.settings.gradientrecomputefrequency) {
diff --git a/src/qpsolver/pricing.hpp b/src/qpsolver/pricing.hpp
index 293bc505d0..85fdae3bb4 100644
--- a/src/qpsolver/pricing.hpp
+++ b/src/qpsolver/pricing.hpp
@@ -8,6 +8,7 @@ class Pricing {
   virtual HighsInt price(const Vector& x, const Vector& gradient) = 0;
   virtual void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                               HighsInt q) = 0;
+  virtual void recompute() = 0;
   virtual ~Pricing() {}
 };
 
diff --git a/src/qpsolver/quass.cpp b/src/qpsolver/quass.cpp
index e2621a3acb..44006e1323 100644
--- a/src/qpsolver/quass.cpp
+++ b/src/qpsolver/quass.cpp
@@ -266,6 +266,20 @@ static double compute_dual_violation(Instance& instance, Vector& primal, Vector&
 }
 #endif
 
+bool check_reinvert_due(Basis& basis) {
+  // reinvert can be triggered by basis
+  return basis.getreinversionhint();
+}
+
+void reinvert(Basis& basis, CholeskyFactor& factor, Gradient& grad, ReducedCosts& rc, ReducedGradient& rg, std::unique_ptr<Pricing>& pricing) {
+  basis.rebuild();
+  factor.recompute();
+  grad.recompute();
+  rc.recompute();
+  rg.recompute();
+  pricing->recompute();
+}
+
 void Quass::solve(const Vector& x0, const Vector& ra, Basis& b0, HighsTimer& timer) {
 
   //feenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT & ~FE_UNDERFLOW);
@@ -332,6 +346,12 @@ void Quass::solve(const Vector& x0, const Vector& ra, Basis& b0, HighsTimer& tim
     }
     runtime.statistics.num_iterations++;
 
+    // REINVERSION
+    if (check_reinvert_due(basis)) {
+      printf("REINVERT\n");
+      reinvert(basis, factor, gradient, redcosts, redgrad, pricing);
+    }
+
     QpSolverStatus status;
 
     bool zero_curvature_direction = false;
diff --git a/src/qpsolver/reducedcosts.hpp b/src/qpsolver/reducedcosts.hpp
index 36a165a4d5..f3286a2c6c 100644
--- a/src/qpsolver/reducedcosts.hpp
+++ b/src/qpsolver/reducedcosts.hpp
@@ -14,11 +14,6 @@ class ReducedCosts {
   Vector reducedcosts;
   bool uptodate;
 
-  void recompute() {
-    basis.ftran(gradient.getGradient(), reducedcosts);
-    uptodate = true;
-  }
-
  public:
   ReducedCosts(Runtime& rt, Basis& bas, Gradient& grad)
       : basis(bas),
@@ -26,6 +21,11 @@ class ReducedCosts {
         reducedcosts(Vector(rt.instance.num_var)),
         uptodate(false) {}
 
+  void recompute() {
+    basis.ftran(gradient.getGradient(), reducedcosts);
+    uptodate = true;
+  }
+
   Vector& getReducedCosts() {
     if (!uptodate) {
       recompute();
diff --git a/src/qpsolver/reducedgradient.hpp b/src/qpsolver/reducedgradient.hpp
index 0ddd11f0de..b16c9bfb4a 100644
--- a/src/qpsolver/reducedgradient.hpp
+++ b/src/qpsolver/reducedgradient.hpp
@@ -11,12 +11,6 @@ class ReducedGradient {
   Gradient& gradient;
   Basis& basis;
 
-  void recompute() {
-    rg.dim = basis.getinactive().size();
-    basis.Ztprod(gradient.getGradient(), rg);
-    uptodate = true;
-  }
-
  public:
   ReducedGradient(Runtime& rt, Basis& bas, Gradient& grad)
       : rg(rt.instance.num_var), gradient(grad), basis(bas) {}
@@ -28,6 +22,12 @@ class ReducedGradient {
     return rg;
   }
 
+  void recompute() {
+    rg.dim = basis.getinactive().size();
+    basis.Ztprod(gradient.getGradient(), rg);
+    uptodate = true;
+  }
+
   void reduce(const Vector& buffer_d, const HighsInt maxabsd) {
     if (!uptodate) {
       return;
diff --git a/src/qpsolver/settings.hpp b/src/qpsolver/settings.hpp
index b31707f99e..e577b22963 100644
--- a/src/qpsolver/settings.hpp
+++ b/src/qpsolver/settings.hpp
@@ -36,7 +36,7 @@ struct Settings {
   HighsInt reportingfequency = 1;
   Eventhandler<Statistics&> endofiterationevent;
 
-  HighsInt reinvertfrequency = 100;
+  HighsInt reinvertfrequency = 1000;
   HighsInt gradientrecomputefrequency = 100;
   HighsInt reducedgradientrecomputefrequency =
       std::numeric_limits<HighsInt>::infinity();
diff --git a/src/qpsolver/steepestedgepricing.hpp b/src/qpsolver/steepestedgepricing.hpp
index 84a3064263..6a4fbcd86b 100644
--- a/src/qpsolver/steepestedgepricing.hpp
+++ b/src/qpsolver/steepestedgepricing.hpp
@@ -116,6 +116,10 @@ class SteepestEdgePricing : public Pricing {
     return incorrect_weights.size() == 0;
     }
 
+  void recompute() {
+    compute_exact_weights();
+  }
+
   void update_weights(const Vector& aq, const Vector& ep, HighsInt p,
                       HighsInt q) {
     HighsInt rowindex_p = basis.getindexinfactor()[p];

From d27d8fa8993c7ce6fed012a8606c0136b52dbf5a Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Wed, 13 Mar 2024 12:58:17 +0000
Subject: [PATCH 04/10] QUASS: fix rogue output, fix codestyle

---
 src/lp_data/Highs.cpp  | 2 +-
 src/qpsolver/quass.cpp | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
index 0f402e61bb..cdb3929a48 100644
--- a/src/lp_data/Highs.cpp
+++ b/src/lp_data/Highs.cpp
@@ -3406,7 +3406,7 @@ HighsStatus Highs::callSolveQp() {
   settings.iterationlimit = options_.simplex_iteration_limit;
   settings.lambda_zero_threshold = options_.dual_feasibility_tolerance;
 
-  switch(options_.simplex_primal_edge_weight_strategy) {
+  switch (options_.simplex_primal_edge_weight_strategy) {
     case 0:
       settings.pricing = PricingStrategy::DantzigWolfe;
       break;
diff --git a/src/qpsolver/quass.cpp b/src/qpsolver/quass.cpp
index 44006e1323..ce03649126 100644
--- a/src/qpsolver/quass.cpp
+++ b/src/qpsolver/quass.cpp
@@ -348,7 +348,6 @@ void Quass::solve(const Vector& x0, const Vector& ra, Basis& b0, HighsTimer& tim
 
     // REINVERSION
     if (check_reinvert_due(basis)) {
-      printf("REINVERT\n");
       reinvert(basis, factor, gradient, redcosts, redgrad, pricing);
     }
 

From 4e3c90b46a9913a05a57ef5be48857f0406d5224 Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Sun, 24 Mar 2024 15:22:49 +0000
Subject: [PATCH 05/10] QUASS fix bug where qp iteration count wasn't passed
 back to Highs

---
 src/qpsolver/a_asm.cpp   | 3 +--
 src/qpsolver/runtime.hpp | 7 ++++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/qpsolver/a_asm.cpp b/src/qpsolver/a_asm.cpp
index cd48273cde..6a289aad97 100644
--- a/src/qpsolver/a_asm.cpp
+++ b/src/qpsolver/a_asm.cpp
@@ -2,8 +2,7 @@
 #include "qpsolver/quass.hpp"
 
 QpAsmStatus solveqp_actual(Instance& instance, Settings& settings, QpHotstartInformation& startinfo, Statistics& stats, QpModelStatus& status, QpSolution& solution, HighsTimer& qp_timer) {
-  Runtime rt(instance);
-  rt.statistics = stats;
+  Runtime rt(instance, stats);
   rt.settings = settings;
   Quass quass(rt);
 
diff --git a/src/qpsolver/runtime.hpp b/src/qpsolver/runtime.hpp
index 3a9cd361e2..9c694ba950 100644
--- a/src/qpsolver/runtime.hpp
+++ b/src/qpsolver/runtime.hpp
@@ -13,7 +13,7 @@ struct Runtime {
   Instance scaled;
   Instance perturbed;
   Settings settings;
-  Statistics statistics;
+  Statistics& statistics;
 
   Vector primal;
   Vector rowactivity;
@@ -24,14 +24,15 @@ struct Runtime {
   std::vector<BasisStatus> status_var;
   std::vector<BasisStatus> status_con;
 
-  Runtime(Instance& inst)
+  Runtime(Instance& inst, Statistics& stats)
       : instance(inst),
         primal(Vector(instance.num_var)),
         rowactivity(Vector(instance.num_con)),
         dualvar(instance.num_var),
         dualcon(instance.num_con),
         status_var(instance.num_var),
-        status_con(instance.num_con) {}
+        status_con(instance.num_con),
+        statistics(stats) {}
 };
 
 #endif

From c09c0101219a60baaac68328eefeef19c57941ff Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Sun, 24 Mar 2024 15:33:24 +0000
Subject: [PATCH 06/10] QUASS don't recompute steepest edge weights. check for
 large nullspace regularly

---
 src/qpsolver/quass.cpp | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/src/qpsolver/quass.cpp b/src/qpsolver/quass.cpp
index ce03649126..cb0d94442d 100644
--- a/src/qpsolver/quass.cpp
+++ b/src/qpsolver/quass.cpp
@@ -277,7 +277,7 @@ void reinvert(Basis& basis, CholeskyFactor& factor, Gradient& grad, ReducedCosts
   grad.recompute();
   rc.recompute();
   rg.recompute();
-  pricing->recompute();
+  //pricing->recompute();
 }
 
 void Quass::solve(const Vector& x0, const Vector& ra, Basis& b0, HighsTimer& timer) {
@@ -316,11 +316,6 @@ void Quass::solve(const Vector& x0, const Vector& ra, Basis& b0, HighsTimer& tim
 
   runtime.relaxed_for_ratiotest = ratiotest_relax_instance(runtime);
 
-  if (basis.getnuminactive() > 4000) {
-    printf("nullspace too large %" HIGHSINT_FORMAT "\n", basis.getnuminactive());
-    runtime.status = QpModelStatus::LARGE_NULLSPACE;
-    return;
-  }
 
 
   bool atfsep = basis.getnumactive() == runtime.instance.num_var;
@@ -337,6 +332,12 @@ void Quass::solve(const Vector& x0, const Vector& ra, Basis& b0, HighsTimer& tim
       break;
     }
 
+    if (basis.getnuminactive() > 4000) {
+      printf("nullspace too large %" HIGHSINT_FORMAT "\n", basis.getnuminactive());
+      runtime.status = QpModelStatus::LARGE_NULLSPACE;
+    return;
+  }
+
     // LOGGING
     if (runtime.statistics.num_iterations %
             runtime.settings.reportingfequency ==

From 7a799cca6c8b907ec440dfcb0c1f855014d606ea Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Sun, 24 Mar 2024 15:34:08 +0000
Subject: [PATCH 07/10] QUASS problems without constraints (just bounds) no
 longer recompute reduced hessian factorization at every iteration

---
 src/qpsolver/factor.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qpsolver/factor.hpp b/src/qpsolver/factor.hpp
index 40ef603c46..3df7bff334 100644
--- a/src/qpsolver/factor.hpp
+++ b/src/qpsolver/factor.hpp
@@ -204,7 +204,7 @@ class CholeskyFactor {
   }
 
   void solve(Vector& rhs) {
-    if (!uptodate || (numberofreduces >= runtime.instance.num_con / 2 &&
+    if (!uptodate || (numberofreduces >= runtime.instance.num_var / 2 &&
                       !has_negative_eigenvalue)) {
       recompute();
     }

From 7f741d83dcbd2bce5ade8146fcd87af2813766ce Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Sun, 24 Mar 2024 15:37:57 +0000
Subject: [PATCH 08/10] QUASS special phase for instances with no constraints
 (just bounds)

---
 cmake/sources-python.cmake           |  1 +
 cmake/sources.cmake                  |  1 +
 src/qpsolver/a_quass.cpp             | 17 ++++--
 src/qpsolver/feasibility.hpp         |  4 ++
 src/qpsolver/feasibility_bounded.hpp | 90 ++++++++++++++++++++++++++++
 src/qpsolver/settings.hpp            |  2 +-
 6 files changed, 110 insertions(+), 5 deletions(-)
 create mode 100644 src/qpsolver/feasibility_bounded.hpp

diff --git a/cmake/sources-python.cmake b/cmake/sources-python.cmake
index 2002418225..3d6b3ea684 100644
--- a/cmake/sources-python.cmake
+++ b/cmake/sources-python.cmake
@@ -369,6 +369,7 @@ set(highs_headers_python
     src/qpsolver/devexpricing.hpp
     src/qpsolver/eventhandler.hpp
     src/qpsolver/factor.hpp
+    src/qpsolver/feasibility_bounded.hpp
     src/qpsolver/feasibility_highs.hpp
     src/qpsolver/feasibility_quass.hpp
     src/qpsolver/feasibility.hpp
diff --git a/cmake/sources.cmake b/cmake/sources.cmake
index cdf1096991..f8011b86f1 100644
--- a/cmake/sources.cmake
+++ b/cmake/sources.cmake
@@ -372,6 +372,7 @@ set(highs_headers
     qpsolver/devexpricing.hpp
     qpsolver/eventhandler.hpp
     qpsolver/factor.hpp
+    qpsolver/feasibility_bounded.hpp
     qpsolver/feasibility_highs.hpp
     qpsolver/feasibility_quass.hpp
     qpsolver/feasibility.hpp
diff --git a/src/qpsolver/a_quass.cpp b/src/qpsolver/a_quass.cpp
index f89f4506de..f4b091cd9a 100644
--- a/src/qpsolver/a_quass.cpp
+++ b/src/qpsolver/a_quass.cpp
@@ -2,6 +2,7 @@
 #include "qpsolver/a_asm.hpp"
 
 #include "qpsolver/feasibility_highs.hpp"
+#include "qpsolver/feasibility_bounded.hpp"
 
 
 
@@ -26,10 +27,18 @@ QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats, Q
 
   // compute initial feasible point
   QpHotstartInformation startinfo(instance.num_var, instance.num_con);
-  computestartingpoint_highs(instance, settings, stats, modelstatus, startinfo, qp_timer);
-  if (modelstatus == QpModelStatus::INFEASIBLE) {
-    return QpAsmStatus::OK;
-  }
+  if (instance.num_con == 0 && instance.num_var <= 15000) {
+    computestartingpoint_bounded(instance, settings, stats, modelstatus, startinfo, qp_timer);
+    if (modelstatus == QpModelStatus::OPTIMAL) {
+      solution.primal = startinfo.primal;
+      return QpAsmStatus::OK;
+    }
+  } else  {
+    computestartingpoint_highs(instance, settings, stats, modelstatus, startinfo, qp_timer);
+    if (modelstatus == QpModelStatus::INFEASIBLE) {
+      return QpAsmStatus::OK;
+    }
+  } 
 
   // solve
   QpAsmStatus status = solveqp_actual(instance, settings, startinfo, stats, modelstatus, solution, qp_timer);
diff --git a/src/qpsolver/feasibility.hpp b/src/qpsolver/feasibility.hpp
index 91265bb97b..e393b3a7ec 100644
--- a/src/qpsolver/feasibility.hpp
+++ b/src/qpsolver/feasibility.hpp
@@ -3,6 +3,7 @@
 
 #include "runtime.hpp"
 #include "crashsolution.hpp"
+#include "feasibility_bounded.hpp"
 #include "feasibility_highs.hpp"
 #include "feasibility_quass.hpp"
 
@@ -16,6 +17,9 @@ void computestartingpoint(Runtime& runtime, QpHotstartInformation& result) {
         case Phase1Strategy::QUASS:
             computestartingpoint_quass(runtime, result);
             break;
+        case Phase1Strategy::BOUNDED:
+            computestartingpoint_bounded(runtime.instance, runtime.settings, runtime.statistics, runtime.status, result, qp_timer);
+            break;
     }
 }
 
diff --git a/src/qpsolver/feasibility_bounded.hpp b/src/qpsolver/feasibility_bounded.hpp
new file mode 100644
index 0000000000..3442540e89
--- /dev/null
+++ b/src/qpsolver/feasibility_bounded.hpp
@@ -0,0 +1,90 @@
+#ifndef __SRC_LIB_FEASIBILITYBOUNDED_HPP__
+#define __SRC_LIB_FEASIBILITYBOUNDED_HPP__
+
+#include "Highs.h"
+#include "qpsolver/a_asm.hpp"
+#include "qpsolver/crashsolution.hpp"
+
+static void computestartingpoint_bounded(Instance& instance, Settings& settings, Statistics& stats, QpModelStatus& modelstatus, QpHotstartInformation& result, HighsTimer& timer) {
+  // compute initial feasible point for problems with bounds only (no general linear constraints)
+
+  // compute  Qx + c = 0 --> x = Q^-1c
+  std::vector<double> L;
+  L.resize(instance.num_var * instance.num_var);
+
+  // compute cholesky factorization of Q
+  for (size_t col = 0; col < instance.num_var; col++) {
+    for (size_t idx = instance.Q.mat.start[col]; idx < instance.Q.mat.start[col+1]; idx++) {
+      double sum = 0;
+      size_t row = instance.Q.mat.index[idx]; 
+      if (row == col) {
+        for (size_t k = 0; k < row; k++)
+          sum += L[k * instance.num_var + row] * L[k * instance.num_var + row];
+        L[row * instance.num_var + row] = sqrt(instance.Q.mat.value[idx] - sum);
+      } else {
+        for (size_t k = 0; k < row; k++)
+          sum += (L[k * instance.num_var + col] * L[k * instance.num_var + row]);
+        L[row * instance.num_var + col] =
+            (instance.Q.mat.value[idx] - sum) / L[row * instance.num_var + row];
+      }
+    }
+  }
+
+  // solve for c
+  Vector res = instance.c;
+  for (HighsInt r = 0; r <res.dim; r++) {
+    for (HighsInt j = 0; j < r; j++) {
+      res.value[r] -= res.value[j] * L[j * instance.num_var + r];
+    }
+    res.value[r] /= L[r * instance.num_var + r];
+  }
+
+  for (HighsInt i = res.dim - 1; i >= 0; i--) {
+    double sum = 0.0;
+    for (HighsInt j = res.dim - 1; j > i; j--) {
+      sum += res.value[j] * L[i * instance.num_var + j];
+    }
+    res.value[i] = (res.value[i] - sum) / L[i * instance.num_var + i];
+  }
+
+  // project solution to bounds and collect active bounds
+  Vector x0(instance.num_var);
+  Vector ra(instance.num_con);
+  std::vector<HighsInt> initialactive;
+  std::vector<HighsInt> initialinactive;
+  std::vector<BasisStatus> atlower;
+
+  for (int i=0; i<instance.num_var; i++) {
+    if (res.value[i] <= instance.var_lo[i]) {
+        res.value[i] = instance.var_lo[i];
+        initialactive.push_back(i + instance.num_con);
+        atlower.push_back(BasisStatus::ActiveAtLower);
+    } else if (res.value[i] >= instance.var_up[i]) {
+        res.value[i] = instance.var_up[i];
+        initialactive.push_back(i + instance.num_con);
+        atlower.push_back(BasisStatus::ActiveAtUpper);
+    } else {
+        initialinactive.push_back(i + instance.num_con);
+    }
+    if (fabs(res.value[i]) > 10E-5) {
+      x0.value[i] = res.value[i];
+      x0.index[x0.num_nz++] = i;
+    }
+  }
+
+  // if no bounds are active, solution lies in the interior -> optimal
+  if (initialactive.size() == 0) {
+    modelstatus = QpModelStatus::OPTIMAL;
+  }
+
+  assert((HighsInt)(initialactive.size() + initialinactive.size()) ==
+         instance.num_var);
+
+  result.status = atlower;
+  result.active = initialactive;
+  result.inactive = initialinactive;
+  result.primal = x0;
+  result.rowact = ra;
+}
+
+#endif
diff --git a/src/qpsolver/settings.hpp b/src/qpsolver/settings.hpp
index e577b22963..5aa8080867 100644
--- a/src/qpsolver/settings.hpp
+++ b/src/qpsolver/settings.hpp
@@ -10,7 +10,7 @@ enum class PricingStrategy { SteepestEdge, DantzigWolfe, Devex };
 
 enum class OutputLevel { LIGHT, MEDIUM, HEAVY };
 
-enum class Phase1Strategy { HIGHS, QUASS };
+enum class Phase1Strategy { HIGHS, QUASS, BOUNDED };
 
 struct Settings {
   RatiotestStrategy ratiotest = RatiotestStrategy::TwoPass;

From 0cde51a1a1daf9b81bb5707f89ca0177d4b99d57 Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Sun, 24 Mar 2024 16:58:02 +0000
Subject: [PATCH 09/10] QUASS Fixed bug in phase1 for unconstrained (but
 bounded) problems

---
 src/qpsolver/feasibility_bounded.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qpsolver/feasibility_bounded.hpp b/src/qpsolver/feasibility_bounded.hpp
index 3442540e89..d51970b6cc 100644
--- a/src/qpsolver/feasibility_bounded.hpp
+++ b/src/qpsolver/feasibility_bounded.hpp
@@ -31,7 +31,7 @@ static void computestartingpoint_bounded(Instance& instance, Settings& settings,
   }
 
   // solve for c
-  Vector res = instance.c;
+  Vector res = -instance.c;
   for (HighsInt r = 0; r <res.dim; r++) {
     for (HighsInt j = 0; j < r; j++) {
       res.value[r] -= res.value[j] * L[j * instance.num_var + r];

From 69ad5d3fe2a42e668f5deb46c59d77dfd285560e Mon Sep 17 00:00:00 2001
From: feldmeier <michael.feldmeier.grad@gmail.com>
Date: Sun, 24 Mar 2024 18:15:35 +0000
Subject: [PATCH 10/10] QUASS detect unboundedness in special-case phase1

---
 src/qpsolver/a_quass.cpp             |  3 +++
 src/qpsolver/feasibility_bounded.hpp | 12 +++++++++++-
 2 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/src/qpsolver/a_quass.cpp b/src/qpsolver/a_quass.cpp
index f4b091cd9a..06eaa93b0d 100644
--- a/src/qpsolver/a_quass.cpp
+++ b/src/qpsolver/a_quass.cpp
@@ -33,6 +33,9 @@ QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats, Q
       solution.primal = startinfo.primal;
       return QpAsmStatus::OK;
     }
+    if (modelstatus == QpModelStatus::UNBOUNDED) {
+      return QpAsmStatus::OK;
+    }
   } else  {
     computestartingpoint_highs(instance, settings, stats, modelstatus, startinfo, qp_timer);
     if (modelstatus == QpModelStatus::INFEASIBLE) {
diff --git a/src/qpsolver/feasibility_bounded.hpp b/src/qpsolver/feasibility_bounded.hpp
index d51970b6cc..5aa32d203b 100644
--- a/src/qpsolver/feasibility_bounded.hpp
+++ b/src/qpsolver/feasibility_bounded.hpp
@@ -55,7 +55,17 @@ static void computestartingpoint_bounded(Instance& instance, Settings& settings,
   std::vector<BasisStatus> atlower;
 
   for (int i=0; i<instance.num_var; i++) {
-    if (res.value[i] <= instance.var_lo[i]) {
+    if (res.value[i] > 0.5/settings.hessianregularizationfactor 
+        && instance.var_up[i] == std::numeric_limits<double>::infinity()
+        && instance.c.value[i] < 0.0) {
+          modelstatus = QpModelStatus::UNBOUNDED;
+          return;
+    } else if (res.value[i] < 0.5/settings.hessianregularizationfactor 
+        && instance.var_lo[i] == std::numeric_limits<double>::infinity()
+        && instance.c.value[i] > 0.0) {
+          modelstatus = QpModelStatus::UNBOUNDED;
+          return;
+    } else if (res.value[i] <= instance.var_lo[i]) {
         res.value[i] = instance.var_lo[i];
         initialactive.push_back(i + instance.num_con);
         atlower.push_back(BasisStatus::ActiveAtLower);