From 86e0cb48451123a0be4ccdec2d8f50436c44e56d Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Wed, 12 Feb 2025 16:43:17 +0100
Subject: [PATCH 01/28] Skip redundant constraints in
 HPresolve::strengthenInequalities

---
 src/presolve/HPresolve.cpp | 41 ++++++++++++++++++++++----------------
 1 file changed, 24 insertions(+), 17 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index e74b0d4aaf..a1f018be78 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -5176,9 +5176,11 @@ HPresolve::Result HPresolve::strengthenInequalities(
     double scale;
 
     if (model->row_lower_[row] != -kHighsInf) {
+      // ax >= lb, treat as -ax <= -lb --> -ax + lb <= 0
       maxviolation = model->row_lower_[row];
       scale = -1.0;
     } else {
+      // ax <= ub --> ax - ub <= 0
       maxviolation = -model->row_upper_[row];
       scale = 1.0;
     }
@@ -5192,6 +5194,7 @@ HPresolve::Result HPresolve::strengthenInequalities(
     reducedcost.reserve(rowsize[row]);
     upper.reserve(rowsize[row]);
     indices.reserve(rowsize[row]);
+    positions.reserve(rowsize[row]);
     stack.reserve(rowsize[row]);
     stack.push_back(rowroot[row]);
 
@@ -5205,25 +5208,25 @@ HPresolve::Result HPresolve::strengthenInequalities(
       if (ARleft[pos] != -1) stack.push_back(ARleft[pos]);
 
       int8_t comp;
-      double weight;
-      double ub;
-      weight = Avalue[pos] * scale;
       HighsInt col = Acol[pos];
-      ub = model->col_upper_[col] - model->col_lower_[col];
+      double weight = Avalue[pos] * scale;
+      double ub = model->col_upper_[col] - model->col_lower_[col];
 
-      skiprow = ub == kHighsInf;
+      skiprow = (ub == kHighsInf) ||
+                (weight > 0 && model->col_upper_[col] == kHighsInf) ||
+                (weight < 0 && model->col_lower_[col] == -kHighsInf);
       if (skiprow) break;
 
+      // compute maximum violation
+      // scale =  1:  ax <=  ub --> violation =  ax - ub > 0
+      // scale = -1: -ax <= -lb --> violation = -ax + lb > 0
+      // this means that for scale = 1 we sum up an upper bound on constraint
+      // activity, and for scale = -1 we sum up a lower bound on constraint
+      // activity.
       if (weight > 0) {
-        skiprow = model->col_upper_[col] == kHighsInf;
-        if (skiprow) break;
-
         comp = 1;
         maxviolation += model->col_upper_[col] * weight;
       } else {
-        skiprow = model->col_lower_[col] == -kHighsInf;
-        if (skiprow) break;
-
         comp = -1;
         maxviolation += model->col_lower_[col] * weight;
         weight = -weight;
@@ -5253,6 +5256,12 @@ HPresolve::Result HPresolve::strengthenInequalities(
       continue;
     }
 
+    // maxviolation <= 0 implies that the constraint is redundant:
+    // scale =  1: upper bound on activity <= model->row_upper_[row]
+    // scale = -1: model->row_lower_[row]  <= lower bound on activity
+    // we skip redundant constraints.
+    if (maxviolation <= primal_feastol) continue;
+
     const double smallVal =
         std::max(100 * primal_feastol, primal_feastol * double(maxviolation));
     while (true) {
@@ -5288,16 +5297,14 @@ HPresolve::Result HPresolve::strengthenInequalities(
             return reducedcost[i1] < reducedcost[i2];
           });
 
-      HighsInt coverend = cover.size();
-
       double al = reducedcost[alpos];
-      coefs.resize(coverend);
+      coefs.resize(cover.size());
       double coverrhs =
           std::max(std::ceil(double(lambda / al - primal_feastol)), 1.0);
       HighsCDouble slackupper = -coverrhs;
 
       double step = kHighsInf;
-      for (HighsInt i = 0; i != coverend; ++i) {
+      for (size_t i = 0; i != cover.size(); ++i) {
         coefs[i] =
             std::ceil(std::min(reducedcost[cover[i]], double(lambda)) / al -
                       options->small_matrix_value);
@@ -5311,7 +5318,7 @@ HPresolve::Result HPresolve::strengthenInequalities(
       reducedcost.push_back(step);
       upper.push_back(double(slackupper));
 
-      for (HighsInt i = 0; i != coverend; ++i)
+      for (size_t i = 0; i != cover.size(); ++i)
         reducedcost[cover[i]] -= step * coefs[i];
 
       indices.erase(std::remove_if(indices.begin(), indices.end(),
@@ -5332,7 +5339,7 @@ HPresolve::Result HPresolve::strengthenInequalities(
                   indices.end());
     if (indices.empty()) continue;
 
-    if (scale == -1.0) {
+    if (scale < 0) {
       HighsCDouble lhs = model->row_lower_[row];
       for (HighsInt i : indices) {
         double coefdelta = double(reducedcost[i] - maxviolation);

From f5f697fb236e7bf1fda31514ef120696ea084ea4 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Wed, 12 Feb 2025 16:57:47 +0100
Subject: [PATCH 02/28] Add test

---
 check/TestMipSolver.cpp  |  11 +
 check/instances/2171.mps | 717 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 728 insertions(+)
 create mode 100644 check/instances/2171.mps

diff --git a/check/TestMipSolver.cpp b/check/TestMipSolver.cpp
index 080f3d124b..a53a4da2a1 100644
--- a/check/TestMipSolver.cpp
+++ b/check/TestMipSolver.cpp
@@ -807,3 +807,14 @@ TEST_CASE("issue-2122", "[highs_test_mip_solver]") {
   const double optimal_objective = -187612.944194;
   solve(highs, kHighsOnString, require_model_status, optimal_objective);
 }
+
+TEST_CASE("issue-2171", "[highs_test_mip_solver]") {
+  std::string filename = std::string(HIGHS_DIR) + "/check/instances/2171.mps";
+  Highs highs;
+  highs.setOptionValue("mip_rel_gap", 0);
+  highs.setOptionValue("mip_abs_gap", 0);
+  highs.readModel(filename);
+  const HighsModelStatus require_model_status = HighsModelStatus::kOptimal;
+  const double optimal_objective = -22375.7585461;
+  solve(highs, kHighsOnString, require_model_status, optimal_objective);
+}
diff --git a/check/instances/2171.mps b/check/instances/2171.mps
new file mode 100644
index 0000000000..2a3e3e96db
--- /dev/null
+++ b/check/instances/2171.mps
@@ -0,0 +1,717 @@
+*SENSE:Minimize
+NAME          seed
+ROWS
+ N  OBJ
+ G  _C1
+ L  _C2
+ G  _C3
+ L  _C4
+ G  _C5
+ G  _C6
+ L  _C7
+ G  _C8
+ G  _C9
+ L  _C10
+ G  _C11
+ G  _C12
+ L  _C13
+ G  _C14
+ G  _C15
+ G  _C16
+ L  _C17
+ G  _C18
+ L  _C19
+ G  _C20
+ G  _C21
+ G  _C22
+ L  _C23
+ L  _C24
+ G  _C25
+ G  _C26
+ G  _C27
+ L  _C28
+ L  _C29
+ L  _C30
+ L  _C31
+ L  _C32
+ L  _C33
+ G  _C34
+ L  _C35
+ L  _C36
+ L  _C37
+ G  _C38
+ L  _C39
+ G  _C40
+ G  _C41
+ G  _C42
+ L  _C43
+ L  _C44
+ G  _C45
+ L  _C46
+ L  _C47
+ L  _C48
+ G  _C49
+ L  _C50
+ G  _C51
+ L  _C52
+ L  _C53
+ G  _C54
+ L  _C55
+ L  _C56
+ L  _C57
+ L  _C58
+ G  _C59
+ L  _C60
+ L  _C61
+ G  _C62
+ G  _C63
+ L  _C64
+ G  _C65
+ G  _C66
+ L  _C67
+ G  _C68
+ L  _C69
+ G  _C70
+ L  _C71
+ L  _C72
+ G  _C73
+ G  _C74
+ L  _C75
+ L  _C76
+ G  _C77
+ L  _C78
+ L  _C79
+ L  _C80
+ L  _C81
+ L  _C82
+ G  _C83
+ G  _C84
+ L  _C85
+ L  _C86
+ L  _C87
+ L  _C88
+ L  _C89
+ L  _C90
+ G  _C91
+ G  _C92
+ G  _C93
+ L  _C94
+ G  _C95
+ L  _C96
+ L  _C97
+ L  _C98
+COLUMNS
+    MARK      'MARKER'                 'INTORG'
+    x0        _C1        1.340000000000e+01
+    x0        _C2       -2.760000000000e+00
+    x0        _C3        5.100000000000e+00
+    x0        _C4        1.152000000000e+01
+    x0        _C5        8.180000000000e+00
+    x0        _C6       -7.650000000000e+00
+    x0        _C7       -1.804000000000e+01
+    x0        _C8        1.074000000000e+01
+    x0        _C9       -1.290000000000e+00
+    x0        _C10       4.710000000000e+00
+    x0        _C11      -5.250000000000e+00
+    x0        _C12      -1.226000000000e+01
+    x0        _C13      -8.890000000000e+00
+    x0        _C14       1.504000000000e+01
+    x0        _C15       1.280000000000e+00
+    x0        _C16       1.159000000000e+01
+    x0        _C17       1.544000000000e+01
+    x0        _C18       1.772000000000e+01
+    x0        _C19       3.640000000000e+00
+    x0        _C20      -4.060000000000e+00
+    x0        _C21       4.460000000000e+00
+    x0        _C22       1.956000000000e+01
+    x0        _C23       1.569000000000e+01
+    x0        _C24       2.820000000000e+00
+    x0        _C25      -1.269000000000e+01
+    x0        _C26      -2.150000000000e+00
+    x0        _C27      -8.480000000000e+00
+    x0        _C28      -1.337000000000e+01
+    x0        _C29      -1.955000000000e+01
+    x0        _C30       1.436000000000e+01
+    x0        _C31       1.324000000000e+01
+    x0        _C32       1.875000000000e+01
+    x0        _C33      -5.800000000000e+00
+    x0        _C34      -1.870000000000e+01
+    x0        _C35       1.912000000000e+01
+    x0        _C36      -1.093000000000e+01
+    x0        _C37      -1.607000000000e+01
+    x0        _C38      -1.501000000000e+01
+    x0        _C39      -1.180000000000e+01
+    x0        _C40       7.850000000000e+00
+    x0        _C41      -2.990000000000e+00
+    x0        _C42      -2.560000000000e+00
+    x0        _C43      -3.160000000000e+00
+    x0        _C44      -1.572000000000e+01
+    x0        _C45      -1.987000000000e+01
+    x0        _C46      -2.710000000000e+00
+    x0        _C47      -1.862000000000e+01
+    x0        _C48       8.510000000000e+00
+    x0        _C49      -1.876000000000e+01
+    x0        _C50      -9.260000000000e+00
+    x0        _C51       6.340000000000e+00
+    x0        _C52      -8.050000000000e+00
+    x0        _C53      -3.200000000000e+00
+    x0        _C54      -1.610000000000e+01
+    x0        _C55       1.545000000000e+01
+    x0        _C56       1.140000000000e+00
+    x0        _C57       8.090000000000e+00
+    x0        _C58       1.640000000000e+00
+    x0        _C59      -4.870000000000e+00
+    x0        _C60       1.505000000000e+01
+    x0        _C61      -1.969000000000e+01
+    x0        _C62      -1.210000000000e+00
+    x0        _C63       4.000000000000e-02
+    x0        _C64      -3.370000000000e+00
+    x0        _C65      -1.870000000000e+01
+    x0        _C66       9.000000000000e-02
+    x0        _C67       4.380000000000e+00
+    x0        _C68       5.920000000000e+00
+    x0        _C69      -8.280000000000e+00
+    x0        _C70       1.511000000000e+01
+    x0        _C71       1.936000000000e+01
+    x0        _C72       8.850000000000e+00
+    x0        _C73      -9.110000000000e+00
+    x0        _C74      -1.670000000000e+01
+    x0        _C75       1.800000000000e+00
+    x0        _C76      -1.536000000000e+01
+    x0        _C77       1.426000000000e+01
+    x0        _C78      -1.222000000000e+01
+    x0        _C79      -5.660000000000e+00
+    x0        _C80       1.860000000000e+00
+    x0        _C81      -5.870000000000e+00
+    x0        _C82       1.893000000000e+01
+    x0        _C83       1.061000000000e+01
+    x0        _C84       1.028000000000e+01
+    x0        _C85       1.674000000000e+01
+    x0        _C86       1.297000000000e+01
+    x0        _C87      -1.448000000000e+01
+    x0        _C88      -2.460000000000e+00
+    x0        _C89       9.730000000000e+00
+    x0        _C90       2.240000000000e+00
+    x0        _C91      -1.416000000000e+01
+    x0        _C92       1.572000000000e+01
+    x0        _C93      -1.233000000000e+01
+    x0        _C94      -7.080000000000e+00
+    x0        _C95       7.940000000000e+00
+    x0        _C96      -9.240000000000e+00
+    x0        _C97       1.146000000000e+01
+    x0        _C98      -8.090000000000e+00
+    x0        OBJ        9.260000000000e+01
+    MARK      'MARKER'                 'INTEND'
+    MARK      'MARKER'                 'INTORG'
+    x1        _C1       -6.810000000000e+00
+    x1        _C2        6.840000000000e+00
+    x1        _C3        2.400000000000e+00
+    x1        _C4        1.833000000000e+01
+    x1        _C5       -5.200000000000e-01
+    x1        _C6        8.370000000000e+00
+    x1        _C7        1.711000000000e+01
+    x1        _C8       -1.300000000000e+01
+    x1        _C9       -1.417000000000e+01
+    x1        _C10      -5.700000000000e+00
+    x1        _C11      -1.775000000000e+01
+    x1        _C12      -5.930000000000e+00
+    x1        _C13      -1.298000000000e+01
+    x1        _C14      -3.140000000000e+00
+    x1        _C15      -1.700000000000e+00
+    x1        _C16      -1.325000000000e+01
+    x1        _C17       1.566000000000e+01
+    x1        _C18      -9.560000000000e+00
+    x1        _C19       1.712000000000e+01
+    x1        _C20      -6.100000000000e+00
+    x1        _C21      -8.400000000000e+00
+    x1        _C22      -1.632000000000e+01
+    x1        _C23       6.600000000000e+00
+    x1        _C24       1.699000000000e+01
+    x1        _C25      -1.220000000000e+01
+    x1        _C26      -7.300000000000e-01
+    x1        _C27      -1.624000000000e+01
+    x1        _C28       1.866000000000e+01
+    x1        _C29       1.603000000000e+01
+    x1        _C30       1.436000000000e+01
+    x1        _C31       1.488000000000e+01
+    x1        _C32       1.915000000000e+01
+    x1        _C33       1.328000000000e+01
+    x1        _C34      -1.010000000000e+01
+    x1        _C35       1.850000000000e+01
+    x1        _C36       4.220000000000e+00
+    x1        _C37       1.310000000000e+01
+    x1        _C38      -1.066000000000e+01
+    x1        _C39       1.475000000000e+01
+    x1        _C40       1.350000000000e+01
+    x1        _C41      -1.653000000000e+01
+    x1        _C42      -4.710000000000e+00
+    x1        _C43      -4.880000000000e+00
+    x1        _C44       1.442000000000e+01
+    x1        _C45      -5.590000000000e+00
+    x1        _C46       1.569000000000e+01
+    x1        _C47      -4.240000000000e+00
+    x1        _C48       9.880000000000e+00
+    x1        _C49      -7.340000000000e+00
+    x1        _C50       7.660000000000e+00
+    x1        _C51      -1.339000000000e+01
+    x1        _C52       9.800000000000e+00
+    x1        _C53      -1.900000000000e+00
+    x1        _C54       1.309000000000e+01
+    x1        _C55      -4.140000000000e+00
+    x1        _C56       7.800000000000e+00
+    x1        _C57       1.903000000000e+01
+    x1        _C58      -6.070000000000e+00
+    x1        _C59      -1.667000000000e+01
+    x1        _C60       1.449000000000e+01
+    x1        _C61       7.790000000000e+00
+    x1        _C62      -1.874000000000e+01
+    x1        _C63      -1.723000000000e+01
+    x1        _C64      -1.500000000000e+00
+    x1        _C65      -1.342000000000e+01
+    x1        _C66      -1.278000000000e+01
+    x1        _C67       4.790000000000e+00
+    x1        _C68      -5.100000000000e-01
+    x1        _C69       1.400000000000e+01
+    x1        _C70       1.222000000000e+01
+    x1        _C71       1.909000000000e+01
+    x1        _C72       1.730000000000e+01
+    x1        _C73      -6.140000000000e+00
+    x1        _C74      -6.570000000000e+00
+    x1        _C75       1.240000000000e+00
+    x1        _C76       1.273000000000e+01
+    x1        _C77      -1.753000000000e+01
+    x1        _C78       1.795000000000e+01
+    x1        _C79      -6.210000000000e+00
+    x1        _C80      -1.094000000000e+01
+    x1        _C81       1.400000000000e+01
+    x1        _C82      -5.790000000000e+00
+    x1        _C83      -1.975000000000e+01
+    x1        _C84      -2.800000000000e-01
+    x1        _C85       3.360000000000e+00
+    x1        _C86       4.710000000000e+00
+    x1        _C87      -3.000000000000e-01
+    x1        _C88      -6.790000000000e+00
+    x1        _C89      -5.440000000000e+00
+    x1        _C90      -7.650000000000e+00
+    x1        _C91      -2.580000000000e+00
+    x1        _C92      -4.600000000000e+00
+    x1        _C93       3.030000000000e+00
+    x1        _C94       1.149000000000e+01
+    x1        _C95      -1.040000000000e+01
+    x1        _C96       1.382000000000e+01
+    x1        _C97      -1.402000000000e+01
+    x1        _C98      -8.600000000000e+00
+    x1        OBJ        8.100000000000e-01
+    MARK      'MARKER'                 'INTEND'
+    x2        _C1       -4.140000000000e+00
+    x2        _C2       -1.465000000000e+01
+    x2        _C3        1.629000000000e+01
+    x2        _C4       -1.327000000000e+01
+    x2        _C5       -6.640000000000e+00
+    x2        _C6        1.569000000000e+01
+    x2        _C7        8.810000000000e+00
+    x2        _C8        1.003000000000e+01
+    x2        _C9        1.979000000000e+01
+    x2        _C10      -1.145000000000e+01
+    x2        _C11      -1.139000000000e+01
+    x2        _C12       6.140000000000e+00
+    x2        _C13      -1.324000000000e+01
+    x2        _C14       1.631000000000e+01
+    x2        _C15      -7.050000000000e+00
+    x2        _C16      -1.740000000000e+01
+    x2        _C17       1.532000000000e+01
+    x2        _C18       4.870000000000e+00
+    x2        _C19       1.695000000000e+01
+    x2        _C20       1.570000000000e+01
+    x2        _C21       8.200000000000e-01
+    x2        _C22       5.210000000000e+00
+    x2        _C23      -6.600000000000e+00
+    x2        _C24      -5.310000000000e+00
+    x2        _C25      -1.139000000000e+01
+    x2        _C26       1.026000000000e+01
+    x2        _C27      -1.507000000000e+01
+    x2        _C28       7.050000000000e+00
+    x2        _C29      -1.160000000000e+01
+    x2        _C30      -1.833000000000e+01
+    x2        _C31      -7.800000000000e+00
+    x2        _C32       1.505000000000e+01
+    x2        _C33       1.841000000000e+01
+    x2        _C34      -7.100000000000e+00
+    x2        _C35       9.190000000000e+00
+    x2        _C36      -1.640000000000e+01
+    x2        _C37      -6.230000000000e+00
+    x2        _C38      -2.170000000000e+00
+    x2        _C39       3.010000000000e+00
+    x2        _C40      -1.292000000000e+01
+    x2        _C41       8.450000000000e+00
+    x2        _C42       1.781000000000e+01
+    x2        _C43       1.038000000000e+01
+    x2        _C44       9.380000000000e+00
+    x2        _C45       1.762000000000e+01
+    x2        _C46      -1.538000000000e+01
+    x2        _C47       1.215000000000e+01
+    x2        _C48       1.950000000000e+00
+    x2        _C49      -1.264000000000e+01
+    x2        _C50      -1.285000000000e+01
+    x2        _C51      -2.720000000000e+00
+    x2        _C52      -1.383000000000e+01
+    x2        _C53      -1.822000000000e+01
+    x2        _C54       3.580000000000e+00
+    x2        _C55       1.912000000000e+01
+    x2        _C56      -1.760000000000e+01
+    x2        _C57      -1.684000000000e+01
+    x2        _C58      -7.450000000000e+00
+    x2        _C59       1.955000000000e+01
+    x2        _C60       7.200000000000e-01
+    x2        _C61      -4.250000000000e+00
+    x2        _C62       1.070000000000e+01
+    x2        _C63       1.980000000000e+01
+    x2        _C64       9.040000000000e+00
+    x2        _C65       1.647000000000e+01
+    x2        _C66      -2.550000000000e+00
+    x2        _C67      -1.448000000000e+01
+    x2        _C68       1.058000000000e+01
+    x2        _C69      -1.275000000000e+01
+    x2        _C70       1.840000000000e+00
+    x2        _C71       3.540000000000e+00
+    x2        _C72      -3.590000000000e+00
+    x2        _C73      -3.560000000000e+00
+    x2        _C74      -1.826000000000e+01
+    x2        _C75      -1.708000000000e+01
+    x2        _C76       1.986000000000e+01
+    x2        _C77       1.352000000000e+01
+    x2        _C78       4.930000000000e+00
+    x2        _C79      -2.280000000000e+00
+    x2        _C80       2.800000000000e-01
+    x2        _C81      -1.860000000000e+00
+    x2        _C82       8.470000000000e+00
+    x2        _C83       9.890000000000e+00
+    x2        _C84       8.050000000000e+00
+    x2        _C85       6.000000000000e-02
+    x2        _C86       2.000000000000e-02
+    x2        _C87       8.620000000000e+00
+    x2        _C88      -8.130000000000e+00
+    x2        _C89      -3.920000000000e+00
+    x2        _C90       7.840000000000e+00
+    x2        _C91       5.010000000000e+00
+    x2        _C92       3.430000000000e+00
+    x2        _C93      -2.870000000000e+00
+    x2        _C94       7.550000000000e+00
+    x2        _C95       1.557000000000e+01
+    x2        _C96      -2.900000000000e-01
+    x2        _C97      -3.310000000000e+00
+    x2        _C98      -9.320000000000e+00
+    x2        OBJ        7.437000000000e+01
+    MARK      'MARKER'                 'INTORG'
+    x3        _C1        1.764000000000e+01
+    x3        _C2        2.710000000000e+00
+    x3        _C3       -6.100000000000e+00
+    x3        _C4       -1.784000000000e+01
+    x3        _C5        1.790000000000e+01
+    x3        _C6       -1.309000000000e+01
+    x3        _C7       -1.333000000000e+01
+    x3        _C8       -5.710000000000e+00
+    x3        _C9       -2.000000000000e-02
+    x3        _C10      -1.067000000000e+01
+    x3        _C11       2.000000000000e+00
+    x3        _C12       8.630000000000e+00
+    x3        _C13      -1.980000000000e+00
+    x3        _C14       8.490000000000e+00
+    x3        _C15       1.025000000000e+01
+    x3        _C16       8.190000000000e+00
+    x3        _C17      -1.681000000000e+01
+    x3        _C18       1.614000000000e+01
+    x3        _C19       4.100000000000e-01
+    x3        _C20       2.520000000000e+00
+    x3        _C21       1.930000000000e+01
+    x3        _C22      -2.400000000000e-01
+    x3        _C23       8.530000000000e+00
+    x3        _C24      -1.052000000000e+01
+    x3        _C25       1.771000000000e+01
+    x3        _C26       7.330000000000e+00
+    x3        _C27       8.200000000000e-01
+    x3        _C28      -1.120000000000e+00
+    x3        _C29      -1.303000000000e+01
+    x3        _C30       3.880000000000e+00
+    x3        _C31      -1.736000000000e+01
+    x3        _C32      -5.860000000000e+00
+    x3        _C33      -5.260000000000e+00
+    x3        _C34       1.108000000000e+01
+    x3        _C35       6.720000000000e+00
+    x3        _C36      -5.150000000000e+00
+    x3        _C37      -1.021000000000e+01
+    x3        _C38      -1.606000000000e+01
+    x3        _C39      -5.980000000000e+00
+    x3        _C40       1.902000000000e+01
+    x3        _C41       1.190000000000e+00
+    x3        _C42      -7.440000000000e+00
+    x3        _C43      -7.170000000000e+00
+    x3        _C44      -1.461000000000e+01
+    x3        _C45       1.540000000000e+00
+    x3        _C46      -6.140000000000e+00
+    x3        _C47      -1.483000000000e+01
+    x3        _C48       4.790000000000e+00
+    x3        _C49      -3.340000000000e+00
+    x3        _C50      -1.911000000000e+01
+    x3        _C51      -3.940000000000e+00
+    x3        _C52      -2.690000000000e+00
+    x3        _C53      -7.390000000000e+00
+    x3        _C54       1.307000000000e+01
+    x3        _C55      -1.950000000000e+01
+    x3        _C56       1.062000000000e+01
+    x3        _C57       1.740000000000e+01
+    x3        _C58      -1.852000000000e+01
+    x3        _C59       1.547000000000e+01
+    x3        _C60       1.176000000000e+01
+    x3        _C61       1.710000000000e+00
+    x3        _C62       8.030000000000e+00
+    x3        _C63      -6.480000000000e+00
+    x3        _C64      -8.030000000000e+00
+    x3        _C65      -8.130000000000e+00
+    x3        _C66       1.476000000000e+01
+    x3        _C67      -6.810000000000e+00
+    x3        _C68       1.010000000000e+01
+    x3        _C69      -1.694000000000e+01
+    x3        _C70       1.149000000000e+01
+    x3        _C71       8.920000000000e+00
+    x3        _C72      -1.481000000000e+01
+    x3        _C73      -8.150000000000e+00
+    x3        _C74       8.730000000000e+00
+    x3        _C75      -1.625000000000e+01
+    x3        _C76      -8.980000000000e+00
+    x3        _C77      -2.020000000000e+00
+    x3        _C78      -3.920000000000e+00
+    x3        _C79      -8.410000000000e+00
+    x3        _C80      -3.510000000000e+00
+    x3        _C81      -8.280000000000e+00
+    x3        _C82      -1.620000000000e+01
+    x3        _C83       1.239000000000e+01
+    x3        _C84       1.398000000000e+01
+    x3        _C85       8.900000000000e+00
+    x3        _C86       6.660000000000e+00
+    x3        _C87      -3.880000000000e+00
+    x3        _C88       5.920000000000e+00
+    x3        _C89      -8.440000000000e+00
+    x3        _C90      -5.860000000000e+00
+    x3        _C91       1.549000000000e+01
+    x3        _C92       1.636000000000e+01
+    x3        _C93       1.766000000000e+01
+    x3        _C94      -1.244000000000e+01
+    x3        _C95       1.255000000000e+01
+    x3        _C96      -1.951000000000e+01
+    x3        _C97      -9.720000000000e+00
+    x3        _C98      -1.427000000000e+01
+    x3        OBJ       -7.574000000000e+01
+    MARK      'MARKER'                 'INTEND'
+    MARK      'MARKER'                 'INTORG'
+    x4        _C1        3.620000000000e+00
+    x4        _C2       -6.880000000000e+00
+    x4        _C3       -1.026000000000e+01
+    x4        _C4        1.855000000000e+01
+    x4        _C5       -6.390000000000e+00
+    x4        _C6       -1.488000000000e+01
+    x4        _C7       -1.490000000000e+01
+    x4        _C8       -2.370000000000e+00
+    x4        _C9       -1.908000000000e+01
+    x4        _C10       1.402000000000e+01
+    x4        _C11      -1.273000000000e+01
+    x4        _C12       3.720000000000e+00
+    x4        _C13       1.806000000000e+01
+    x4        _C14      -7.200000000000e+00
+    x4        _C15      -1.660000000000e+01
+    x4        _C16       1.260000000000e+00
+    x4        _C17      -4.590000000000e+00
+    x4        _C18       1.133000000000e+01
+    x4        _C19       1.300000000000e+00
+    x4        _C20       1.142000000000e+01
+    x4        _C21      -1.149000000000e+01
+    x4        _C22      -1.652000000000e+01
+    x4        _C23       1.326000000000e+01
+    x4        _C24      -7.800000000000e-01
+    x4        _C25       5.920000000000e+00
+    x4        _C26      -1.194000000000e+01
+    x4        _C27       8.400000000000e+00
+    x4        _C28      -5.660000000000e+00
+    x4        _C29      -8.350000000000e+00
+    x4        _C30       5.230000000000e+00
+    x4        _C31      -5.830000000000e+00
+    x4        _C32      -1.913000000000e+01
+    x4        _C33       3.900000000000e+00
+    x4        _C34      -1.280000000000e+01
+    x4        _C35       1.460000000000e+00
+    x4        _C36       1.867000000000e+01
+    x4        _C37       1.474000000000e+01
+    x4        _C38      -1.753000000000e+01
+    x4        _C39       5.830000000000e+00
+    x4        _C40      -5.830000000000e+00
+    x4        _C41       7.790000000000e+00
+    x4        _C42      -1.888000000000e+01
+    x4        _C43       1.361000000000e+01
+    x4        _C44       1.787000000000e+01
+    x4        _C45       5.900000000000e-01
+    x4        _C46      -1.435000000000e+01
+    x4        _C47       1.749000000000e+01
+    x4        _C48       9.080000000000e+00
+    x4        _C49       7.770000000000e+00
+    x4        _C50       1.183000000000e+01
+    x4        _C51      -7.610000000000e+00
+    x4        _C52      -5.480000000000e+00
+    x4        _C53      -7.040000000000e+00
+    x4        _C54      -8.760000000000e+00
+    x4        _C55      -9.500000000000e-01
+    x4        _C56      -2.230000000000e+00
+    x4        _C57       9.990000000000e+00
+    x4        _C58       1.582000000000e+01
+    x4        _C59      -1.104000000000e+01
+    x4        _C60      -9.090000000000e+00
+    x4        _C61       1.813000000000e+01
+    x4        _C62       1.051000000000e+01
+    x4        _C63       8.920000000000e+00
+    x4        _C64       1.691000000000e+01
+    x4        _C65       9.580000000000e+00
+    x4        _C66       9.180000000000e+00
+    x4        _C67       1.583000000000e+01
+    x4        _C68       2.490000000000e+00
+    x4        _C69       1.256000000000e+01
+    x4        _C70      -1.404000000000e+01
+    x4        _C71       4.130000000000e+00
+    x4        _C72       8.550000000000e+00
+    x4        _C73      -1.204000000000e+01
+    x4        _C74      -8.820000000000e+00
+    x4        _C75       5.400000000000e-01
+    x4        _C76       1.760000000000e+01
+    x4        _C77      -1.271000000000e+01
+    x4        _C78      -7.680000000000e+00
+    x4        _C79       1.494000000000e+01
+    x4        _C80       6.050000000000e+00
+    x4        _C81       1.763000000000e+01
+    x4        _C82       1.535000000000e+01
+    x4        _C83      -1.625000000000e+01
+    x4        _C84       6.870000000000e+00
+    x4        _C85      -1.560000000000e+00
+    x4        _C86       3.150000000000e+00
+    x4        _C87       1.113000000000e+01
+    x4        _C88       1.362000000000e+01
+    x4        _C89      -1.259000000000e+01
+    x4        _C90       1.124000000000e+01
+    x4        _C91      -9.740000000000e+00
+    x4        _C92       9.950000000000e+00
+    x4        _C93      -1.756000000000e+01
+    x4        _C94       1.072000000000e+01
+    x4        _C95       7.510000000000e+00
+    x4        _C96       9.230000000000e+00
+    x4        _C97      -7.070000000000e+00
+    x4        _C98       2.750000000000e+00
+    x4        OBJ       -1.143000000000e+01
+    MARK      'MARKER'                 'INTEND'
+RHS
+    RHS       _C1        7.704000000000e+01
+    RHS       _C2       -1.601000000000e+01
+    RHS       _C3        7.370000000000e+01
+    RHS       _C4        1.586000000000e+01
+    RHS       _C5       -3.274000000000e+01
+    RHS       _C6       -4.474000000000e+01
+    RHS       _C7       -8.043000000000e+01
+    RHS       _C8       -8.811000000000e+01
+    RHS       _C9        6.287000000000e+01
+    RHS       _C10      -3.800000000000e-01
+    RHS       _C11       5.304000000000e+01
+    RHS       _C12       2.398000000000e+01
+    RHS       _C13      -6.376000000000e+01
+    RHS       _C14      -1.804000000000e+01
+    RHS       _C15       4.328000000000e+01
+    RHS       _C16      -5.559000000000e+01
+    RHS       _C17      -9.775000000000e+01
+    RHS       _C18       7.473000000000e+01
+    RHS       _C19       5.981000000000e+01
+    RHS       _C20       8.278000000000e+01
+    RHS       _C21      -1.600000000000e-01
+    RHS       _C22      -7.803000000000e+01
+    RHS       _C23      -2.654000000000e+01
+    RHS       _C24      -1.725000000000e+01
+    RHS       _C25       2.388000000000e+01
+    RHS       _C26       7.512000000000e+01
+    RHS       _C27       7.503000000000e+01
+    RHS       _C28      -5.365000000000e+01
+    RHS       _C29       7.422000000000e+01
+    RHS       _C30       8.743000000000e+01
+    RHS       _C31      -9.644000000000e+01
+    RHS       _C32      -2.519000000000e+01
+    RHS       _C33       9.879000000000e+01
+    RHS       _C34       6.983000000000e+01
+    RHS       _C35       6.159000000000e+01
+    RHS       _C36       2.694000000000e+01
+    RHS       _C37       5.155000000000e+01
+    RHS       _C38       9.375000000000e+01
+    RHS       _C39      -1.479000000000e+01
+    RHS       _C40       7.601000000000e+01
+    RHS       _C41       3.679000000000e+01
+    RHS       _C42       6.707000000000e+01
+    RHS       _C43       1.250000000000e+00
+    RHS       _C44       3.029000000000e+01
+    RHS       _C45       2.800000000000e+00
+    RHS       _C46      -2.180000000000e+01
+    RHS       _C47      -7.917000000000e+01
+    RHS       _C48       1.299000000000e+01
+    RHS       _C49      -8.319000000000e+01
+    RHS       _C50       6.803000000000e+01
+    RHS       _C51      -7.105000000000e+01
+    RHS       _C52       4.473000000000e+01
+    RHS       _C53       3.513000000000e+01
+    RHS       _C54      -1.420000000000e+00
+    RHS       _C55      -3.600000000000e+00
+    RHS       _C56      -2.532000000000e+01
+    RHS       _C57       9.410000000000e+00
+    RHS       _C58       2.140000000000e+00
+    RHS       _C59       7.600000000000e+01
+    RHS       _C60      -3.150000000000e+01
+    RHS       _C61      -9.743000000000e+01
+    RHS       _C62       9.303000000000e+01
+    RHS       _C63       6.213000000000e+01
+    RHS       _C64       4.019000000000e+01
+    RHS       _C65       6.691000000000e+01
+    RHS       _C66       1.925000000000e+01
+    RHS       _C67       8.789000000000e+01
+    RHS       _C68       4.610000000000e+00
+    RHS       _C69       1.242000000000e+01
+    RHS       _C70       8.063000000000e+01
+    RHS       _C71       9.538000000000e+01
+    RHS       _C72      -1.581000000000e+01
+    RHS       _C73      -5.479000000000e+01
+    RHS       _C74       2.576000000000e+01
+    RHS       _C75      -6.193000000000e+01
+    RHS       _C76       8.649000000000e+01
+    RHS       _C77       8.911000000000e+01
+    RHS       _C78      -8.219000000000e+01
+    RHS       _C79      -5.545000000000e+01
+    RHS       _C80       2.457000000000e+01
+    RHS       _C81       1.380000000000e+00
+    RHS       _C82       9.627000000000e+01
+    RHS       _C83      -2.008000000000e+01
+    RHS       _C84       1.059000000000e+01
+    RHS       _C85      -6.302000000000e+01
+    RHS       _C86      -8.181000000000e+01
+    RHS       _C87       6.376000000000e+01
+    RHS       _C88      -1.916000000000e+01
+    RHS       _C89       7.153000000000e+01
+    RHS       _C90      -2.240000000000e+00
+    RHS       _C91      -5.183000000000e+01
+    RHS       _C92       2.430000000000e+00
+    RHS       _C93       3.911000000000e+01
+    RHS       _C94       2.501000000000e+01
+    RHS       _C95      -7.167000000000e+01
+    RHS       _C96      -6.749000000000e+01
+    RHS       _C97      -5.058000000000e+01
+    RHS       _C98      -3.153000000000e+01
+BOUNDS
+ LO BND       x0        -2.000000000000e+02
+ UP BND       x0         2.000000000000e+02
+ LO BND       x1        -2.000000000000e+02
+ UP BND       x1         2.000000000000e+02
+ LO BND       x2        -2.000000000000e+02
+ UP BND       x2         2.000000000000e+02
+ LO BND       x3        -2.000000000000e+02
+ UP BND       x3         2.000000000000e+02
+ LO BND       x4        -2.000000000000e+02
+ UP BND       x4         2.000000000000e+02
+ENDATA

From d2e1edb509dd857931425eed944b932cb621940c Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Thu, 13 Feb 2025 08:57:03 +0100
Subject: [PATCH 03/28] Simplify a little more

---
 src/presolve/HPresolve.cpp | 46 +++++++++++++++++---------------------
 1 file changed, 21 insertions(+), 25 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index a1f018be78..5465384afd 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -5201,20 +5201,20 @@ HPresolve::Result HPresolve::strengthenInequalities(
     bool skiprow = false;
 
     while (!stack.empty()) {
+      // pop element from stack
       HighsInt pos = stack.back();
       stack.pop_back();
 
+      // add non-zeros to stack
       if (ARright[pos] != -1) stack.push_back(ARright[pos]);
       if (ARleft[pos] != -1) stack.push_back(ARleft[pos]);
 
-      int8_t comp;
+      // get column index
       HighsInt col = Acol[pos];
-      double weight = Avalue[pos] * scale;
-      double ub = model->col_upper_[col] - model->col_lower_[col];
 
-      skiprow = (ub == kHighsInf) ||
-                (weight > 0 && model->col_upper_[col] == kHighsInf) ||
-                (weight < 0 && model->col_lower_[col] == -kHighsInf);
+      // skip row if a column bound is not finite
+      skiprow = model->col_lower_[col] == -kHighsInf ||
+                model->col_upper_[col] == kHighsInf;
       if (skiprow) break;
 
       // compute maximum violation
@@ -5223,6 +5223,9 @@ HPresolve::Result HPresolve::strengthenInequalities(
       // this means that for scale = 1 we sum up an upper bound on constraint
       // activity, and for scale = -1 we sum up a lower bound on constraint
       // activity.
+      int8_t comp;
+      double weight = Avalue[pos] * scale;
+      double ub = model->col_upper_[col] - model->col_lower_[col];
       if (weight > 0) {
         comp = 1;
         maxviolation += model->col_upper_[col] * weight;
@@ -5339,26 +5342,10 @@ HPresolve::Result HPresolve::strengthenInequalities(
                   indices.end());
     if (indices.empty()) continue;
 
-    if (scale < 0) {
-      HighsCDouble lhs = model->row_lower_[row];
+    auto updateNonZeros = [&](HighsInt row, HighsCDouble& rhs,
+                              HighsInt direction) {
       for (HighsInt i : indices) {
-        double coefdelta = double(reducedcost[i] - maxviolation);
-        HighsInt pos = positions[i];
-
-        if (complementation[i] == -1) {
-          lhs -= coefdelta * model->col_lower_[Acol[pos]];
-          addToMatrix(row, Acol[pos], -coefdelta);
-        } else {
-          lhs += coefdelta * model->col_upper_[Acol[pos]];
-          addToMatrix(row, Acol[pos], coefdelta);
-        }
-      }
-
-      model->row_lower_[row] = double(lhs);
-    } else {
-      HighsCDouble rhs = model->row_upper_[row];
-      for (HighsInt i : indices) {
-        double coefdelta = double(reducedcost[i] - maxviolation);
+        double coefdelta = direction * double(reducedcost[i] - maxviolation);
         HighsInt pos = positions[i];
 
         if (complementation[i] == -1) {
@@ -5369,7 +5356,16 @@ HPresolve::Result HPresolve::strengthenInequalities(
           addToMatrix(row, Acol[pos], -coefdelta);
         }
       }
+    };
 
+    // update / add non-zeros
+    if (scale < 0) {
+      HighsCDouble lhs = model->row_lower_[row];
+      updateNonZeros(row, lhs, HighsInt{-1});
+      model->row_lower_[row] = double(lhs);
+    } else {
+      HighsCDouble rhs = model->row_upper_[row];
+      updateNonZeros(row, rhs, HighsInt{1});
       model->row_upper_[row] = double(rhs);
     }
 

From f3b7952f068dfb8af8523badcfac90df56cd3d27 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Thu, 13 Feb 2025 10:18:56 +0100
Subject: [PATCH 04/28] Use row presolve

---
 src/presolve/HPresolve.cpp | 11 +++++++----
 src/presolve/HPresolve.h   |  3 ++-
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index 5465384afd..d87ad1b03d 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -4203,7 +4203,8 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
 
       if (mipsolver != nullptr) {
         HighsInt num_strengthened = -1;
-        HPRESOLVE_CHECKED_CALL(strengthenInequalities(num_strengthened));
+        HPRESOLVE_CHECKED_CALL(
+            strengthenInequalities(postsolve_stack, num_strengthened));
         assert(num_strengthened >= 0);
         if (num_strengthened > 0)
           highsLogDev(options->log_options, HighsLogType::kInfo,
@@ -5143,7 +5144,7 @@ HPresolve::Result HPresolve::removeDoubletonEquations(
 }
 
 HPresolve::Result HPresolve::strengthenInequalities(
-    HighsInt& num_strengthened) {
+    HighsPostsolveStack& postsolve_stack, HighsInt& num_strengthened) {
   std::vector<int8_t> complementation;
   std::vector<double> reducedcost;
   std::vector<double> upper;
@@ -5262,8 +5263,10 @@ HPresolve::Result HPresolve::strengthenInequalities(
     // maxviolation <= 0 implies that the constraint is redundant:
     // scale =  1: upper bound on activity <= model->row_upper_[row]
     // scale = -1: model->row_lower_[row]  <= lower bound on activity
-    // we skip redundant constraints.
-    if (maxviolation <= primal_feastol) continue;
+    if (maxviolation <= primal_feastol) {
+      HPRESOLVE_CHECKED_CALL(rowPresolve(postsolve_stack, row));
+      continue;
+    }
 
     const double smallVal =
         std::max(100 * primal_feastol, primal_feastol * double(maxviolation));
diff --git a/src/presolve/HPresolve.h b/src/presolve/HPresolve.h
index 13bf2f726f..2240f79981 100644
--- a/src/presolve/HPresolve.h
+++ b/src/presolve/HPresolve.h
@@ -360,7 +360,8 @@ class HPresolve {
 
   Result removeDoubletonEquations(HighsPostsolveStack& postsolve_stack);
 
-  Result strengthenInequalities(HighsInt& num_strenghtened);
+  Result strengthenInequalities(HighsPostsolveStack& postsolve_stack,
+                                HighsInt& num_strenghtened);
 
   HighsInt detectImpliedIntegers();
 

From e2957b663f7dd9345f4d0a7b07d9b6c358a093c4 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Thu, 13 Feb 2025 11:50:13 +0100
Subject: [PATCH 05/28] Minor changes to HPresolve::transformColumn

---
 src/presolve/HPresolve.cpp | 25 ++++++++++++++++++-------
 1 file changed, 18 insertions(+), 7 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index e74b0d4aaf..99956d14b9 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -2298,11 +2298,21 @@ void HPresolve::reinsertEquation(HighsInt row) {
 
 void HPresolve::transformColumn(HighsPostsolveStack& postsolve_stack,
                                 HighsInt col, double scale, double constant) {
+  // replace column x by x = scale * x' + constant
   if (mipsolver != nullptr)
     mipsolver->mipdata_->implications.columnTransformed(col, scale, constant);
 
   postsolve_stack.linearTransform(col, scale, constant);
 
+  // new variable x' has the following bounds:
+  // scale > 0 --> (lb - constant) / scale <= x' <= (ub - constant) / scale
+  // scale < 0 --> (ub - constant) / scale <= x' <= (lb - constant) / scale
+  // each matrix coefficient a will be replaced by scale * a and, therefore, the
+  // contributions to the constraint activities will change as follows:
+  // a * lb --> (a * scale) * (lb - constant) / scale = a * (lb - constant)
+  // a * ub --> (a * scale) * (ub - constant) / scale = a * (ub - constant).
+  // therefore, for now the scaling can be neglected and the bounds on
+  // constraint activities can be updated using the constant term.
   double oldLower = model->col_lower_[col];
   double oldUpper = model->col_upper_[col];
   model->col_upper_[col] -= constant;
@@ -2337,13 +2347,6 @@ void HPresolve::transformColumn(HighsPostsolveStack& postsolve_stack,
   model->col_upper_[col] *= boundScale;
   implColLower[col] *= boundScale;
   implColUpper[col] *= boundScale;
-  if (model->integrality_[col] != HighsVarType::kContinuous) {
-    // we rely on the integrality status being already updated to the newly
-    // scaled column by the caller, if necessary
-    model->col_upper_[col] =
-        std::floor(model->col_upper_[col] + primal_feastol);
-    model->col_lower_[col] = std::ceil(model->col_lower_[col] - primal_feastol);
-  }
 
   if (scale < 0) {
     std::swap(model->col_lower_[col], model->col_upper_[col]);
@@ -2366,6 +2369,14 @@ void HPresolve::transformColumn(HighsPostsolveStack& postsolve_stack,
       model->row_upper_[row] -= rowConstant;
   }
 
+  // use utility functions for rounding scaled bounds of integer-constrained
+  // variables if possible. this should not be done before the preceding bound
+  // updates (scaling and swaps) and matrix updates. we rely on the integrality
+  // status being already updated to the newly scaled column by the caller, if
+  // necessary.
+  changeColLower(col, model->col_lower_[col]);
+  changeColUpper(col, model->col_upper_[col]);
+
   markChangedCol(col);
 }
 

From 4c783d98eb208e2ef79aa686d365aab7385e89b0 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Thu, 13 Feb 2025 14:25:26 +0000
Subject: [PATCH 06/28] Corrected integer values for PAMI and SIP in
 documentation

---
 docs/src/options/definitions.md | 2 +-
 src/lp_data/HighsOptions.h      | 2 +-
 src/simplex/HEkk.cpp            | 8 ++++----
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/docs/src/options/definitions.md b/docs/src/options/definitions.md
index 92ac80c889..5ff9199ac6 100644
--- a/docs/src/options/definitions.md
+++ b/docs/src/options/definitions.md
@@ -122,7 +122,7 @@
 - Default: 0
 
 ## [simplex\_strategy](@id option-simplex_strategy)
-- Strategy for simplex solver 0 => Choose; 1 => Dual (serial); 2 => Dual (PAMI); 3 => Dual (SIP); 4 => Primal
+- Strategy for simplex solver 0 => Choose; 1 => Dual (serial); 2 => Dual (SIP); 3 => Dual (PAMI); 4 => Primal
 - Type: integer
 - Range: {0, 4}
 - Default: 1
diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h
index 1a0b927f63..18c919869f 100644
--- a/src/lp_data/HighsOptions.h
+++ b/src/lp_data/HighsOptions.h
@@ -773,7 +773,7 @@ class HighsOptions : public HighsOptionsStruct {
     record_int = new OptionRecordInt(
         "simplex_strategy",
         "Strategy for simplex solver 0 => Choose; 1 => Dual (serial); 2 => "
-        "Dual (PAMI); 3 => Dual (SIP); 4 => Primal",
+        "Dual (SIP); 3 => Dual (PAMI); 4 => Primal",
         advanced, &simplex_strategy, kSimplexStrategyMin, kSimplexStrategyDual,
         kSimplexStrategyMax);
     records.push_back(record_int);
diff --git a/src/simplex/HEkk.cpp b/src/simplex/HEkk.cpp
index ff70340069..95121d9536 100644
--- a/src/simplex/HEkk.cpp
+++ b/src/simplex/HEkk.cpp
@@ -1105,13 +1105,13 @@ HighsStatus HEkk::solve(const bool force_phase2) {
     if (simplex_strategy == kSimplexStrategyDualTasks) {
       highsLogUser(options_->log_options, HighsLogType::kInfo,
                    "Using EKK parallel dual simplex solver - SIP with "
-                   "concurrency of %" HIGHSINT_FORMAT "\n",
-                   info_.num_concurrency);
+                   "concurrency of %d\n",
+                   int(info_.num_concurrency));
     } else if (simplex_strategy == kSimplexStrategyDualMulti) {
       highsLogUser(options_->log_options, HighsLogType::kInfo,
                    "Using EKK parallel dual simplex solver - PAMI with "
-                   "concurrency of %" HIGHSINT_FORMAT "\n",
-                   info_.num_concurrency);
+                   "concurrency of %d\n",
+                   int(info_.num_concurrency));
     } else {
       highsLogUser(options_->log_options, HighsLogType::kInfo,
                    "Using EKK dual simplex solver - serial\n");

From e501a3025ca8fba9c4de723ae5c53c1a5b35b410 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Fri, 14 Feb 2025 08:50:19 +0100
Subject: [PATCH 07/28] Fix some typos along the way

---
 src/interfaces/highs_c_api.h   | 2 +-
 src/ipm/ipx/lp_solver.h        | 2 +-
 src/lp_data/Highs.cpp          | 2 +-
 src/mip/HighsMipSolverData.cpp | 2 +-
 src/qpsolver/a_asm.cpp         | 2 +-
 5 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/interfaces/highs_c_api.h b/src/interfaces/highs_c_api.h
index 02fe1930b4..2aa17cb6cd 100644
--- a/src/interfaces/highs_c_api.h
+++ b/src/interfaces/highs_c_api.h
@@ -1832,7 +1832,7 @@ HighsInt Highs_getColsByMask(const void* highs, const HighsInt* mask,
  * @param from_row      The first row for which to query data for.
  * @param to_row        The last row (inclusive) for which to query data for.
  * @param num_row       An integer to be populated with the number of rows got
- *                      from the smodel.
+ *                      from the model.
  * @param lower         An array of size [to_row - from_row + 1] for the row
  *                      lower bounds.
  * @param upper         An array of size [to_row - from_row + 1] for the row
diff --git a/src/ipm/ipx/lp_solver.h b/src/ipm/ipx/lp_solver.h
index a7a2da4bd1..2fb8097776 100644
--- a/src/ipm/ipx/lp_solver.h
+++ b/src/ipm/ipx/lp_solver.h
@@ -95,7 +95,7 @@ class LpSolver {
     Int GetBasicSolution(double* x, double* slack, double* y, double* z,
                          Int* cbasis, Int* vbasis) const;
 
-    // Returns/sets all paramters. Without calling SetParameters(), the solver
+    // Returns/sets all parameters. Without calling SetParameters(), the solver
     // uses the default values of a Parameters object.
     Parameters GetParameters() const;
     void SetParameters(Parameters new_parameters);
diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
index 3da151d73a..488da77a5e 100644
--- a/src/lp_data/Highs.cpp
+++ b/src/lp_data/Highs.cpp
@@ -1264,7 +1264,7 @@ HighsStatus Highs::solve() {
   const bool unconstrained_lp = incumbent_lp.a_matrix_.numNz() == 0;
   assert(incumbent_lp.num_row_ || unconstrained_lp);
   // Even if options_.solver == kHighsChooseString in isolation will,
-  // untimately lead to a choice between simplex and IPM, if a basis
+  // ultimately lead to a choice between simplex and IPM, if a basis
   // is available, simplex should surely be chosen.
   const bool solver_will_use_basis = options_.solver == kSimplexString ||
                                      options_.solver == kHighsChooseString;
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 26a15f0a4c..8a4ce2a08b 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -1280,7 +1280,7 @@ void HighsMipSolverData::performRestart() {
   }
 
   // Transform the reference of the objective limit and lower/upper
-  // bounds to the original model, since offset will generally changte
+  // bounds to the original model, since offset will generally change
   // in presolve. Bound changes are transitory, so no real gap change,
   // and no update to P-D integral is necessary
   upper_limit += mipsolver.model_->offset_;
diff --git a/src/qpsolver/a_asm.cpp b/src/qpsolver/a_asm.cpp
index 4fb9706c4e..4dc238e0be 100644
--- a/src/qpsolver/a_asm.cpp
+++ b/src/qpsolver/a_asm.cpp
@@ -47,7 +47,7 @@ std::string qpModelStatusToString(const QpModelStatus qp_model_status) {
     case QpModelStatus::kNotset:
       return "Not set";
     case QpModelStatus::kUndetermined:
-      return "Undertermined";
+      return "Undetermined";
     case QpModelStatus::kOptimal:
       return "Optimal";
     case QpModelStatus::kUnbounded:

From a8edba912fcb7238a849450a20fe78cf24b0acc3 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Fri, 14 Feb 2025 15:13:01 +0100
Subject: [PATCH 08/28] Check if variable is integer

---
 src/presolve/HPresolve.cpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index c6e1aa8ede..933cdf2d20 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -2374,8 +2374,10 @@ void HPresolve::transformColumn(HighsPostsolveStack& postsolve_stack,
   // updates (scaling and swaps) and matrix updates. we rely on the integrality
   // status being already updated to the newly scaled column by the caller, if
   // necessary.
-  changeColLower(col, model->col_lower_[col]);
-  changeColUpper(col, model->col_upper_[col]);
+  if (model->integrality_[col] != HighsVarType::kContinuous) {
+    changeColLower(col, model->col_lower_[col]);
+    changeColUpper(col, model->col_upper_[col]);
+  }
 
   markChangedCol(col);
 }

From c90a2283ea357f83a1613e613bdc21a0f0a18248 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Mon, 17 Feb 2025 09:06:46 +0100
Subject: [PATCH 09/28] Fix comment

---
 src/presolve/HPresolve.cpp | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index 933cdf2d20..e7cde4fdda 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -2369,10 +2369,11 @@ void HPresolve::transformColumn(HighsPostsolveStack& postsolve_stack,
       model->row_upper_[row] -= rowConstant;
   }
 
-  // use utility functions for rounding scaled bounds of integer-constrained
-  // variables if possible. this should not be done before the preceding bound
-  // updates (scaling and swaps) and matrix updates. we rely on the integrality
-  // status being already updated to the newly scaled column by the caller, if
+  // finally, use utility methods for rounding scaled bounds of
+  // integer-constrained variables and updating bounds on constraint activities
+  // accordingly. this should not be done before the preceding bound updates
+  // (scaling and swaps) and matrix updates. we rely on the integrality status
+  // being already updated to the newly scaled column by the caller, if
   // necessary.
   if (model->integrality_[col] != HighsVarType::kContinuous) {
     changeColLower(col, model->col_lower_[col]);

From 3a782e7492ae57cabccd63ce926923e43a726043 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Mon, 17 Feb 2025 09:22:36 +0100
Subject: [PATCH 10/28] Skip part of HPresolve::transformColumn() when constant
 is zero

---
 src/presolve/HPresolve.cpp | 46 +++++++++++++++++++++-----------------
 1 file changed, 25 insertions(+), 21 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index e7cde4fdda..d68e2b5da4 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -2313,28 +2313,32 @@ void HPresolve::transformColumn(HighsPostsolveStack& postsolve_stack,
   // a * ub --> (a * scale) * (ub - constant) / scale = a * (ub - constant).
   // therefore, for now the scaling can be neglected and the bounds on
   // constraint activities can be updated using the constant term.
-  double oldLower = model->col_lower_[col];
-  double oldUpper = model->col_upper_[col];
-  model->col_upper_[col] -= constant;
-  model->col_lower_[col] -= constant;
-
-  for (const HighsSliceNonzero& nonzero : getColumnVector(col)) {
-    impliedRowBounds.updatedVarLower(nonzero.index(), col, nonzero.value(),
-                                     oldLower);
-    impliedRowBounds.updatedVarUpper(nonzero.index(), col, nonzero.value(),
-                                     oldUpper);
-  }
-
-  double oldImplLower = implColLower[col];
-  double oldImplUpper = implColUpper[col];
-  implColLower[col] -= constant;
-  implColUpper[col] -= constant;
+  if (constant != 0.0) {
+    double oldLower = model->col_lower_[col];
+    double oldUpper = model->col_upper_[col];
+    model->col_upper_[col] -= constant;
+    model->col_lower_[col] -= constant;
+
+    for (const HighsSliceNonzero& nonzero : getColumnVector(col)) {
+      impliedRowBounds.updatedVarLower(nonzero.index(), col, nonzero.value(),
+                                       oldLower);
+      impliedRowBounds.updatedVarUpper(nonzero.index(), col, nonzero.value(),
+                                       oldUpper);
+    }
 
-  for (const HighsSliceNonzero& nonzero : getColumnVector(col)) {
-    impliedRowBounds.updatedImplVarLower(nonzero.index(), col, nonzero.value(),
-                                         oldImplLower, colLowerSource[col]);
-    impliedRowBounds.updatedImplVarUpper(nonzero.index(), col, nonzero.value(),
-                                         oldImplUpper, colUpperSource[col]);
+    double oldImplLower = implColLower[col];
+    double oldImplUpper = implColUpper[col];
+    implColLower[col] -= constant;
+    implColUpper[col] -= constant;
+
+    for (const HighsSliceNonzero& nonzero : getColumnVector(col)) {
+      impliedRowBounds.updatedImplVarLower(nonzero.index(), col,
+                                           nonzero.value(), oldImplLower,
+                                           colLowerSource[col]);
+      impliedRowBounds.updatedImplVarUpper(nonzero.index(), col,
+                                           nonzero.value(), oldImplUpper,
+                                           colUpperSource[col]);
+    }
   }
 
   // now apply the scaling, which does not change the contributions to the

From 9975061e0e6c57ec279f1cb86f4ba4549715d39c Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Mon, 17 Feb 2025 11:03:19 +0100
Subject: [PATCH 11/28] Simplify a little more

---
 src/presolve/HPresolve.cpp | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index d68e2b5da4..99e9ddb764 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -5302,12 +5302,13 @@ HPresolve::Result HPresolve::strengthenInequalities(
       cover.clear();
       cover.reserve(indices.size());
 
-      for (HighsInt i = indices.size() - 1; i >= 0; --i) {
-        double delta = upper[indices[i]] * reducedcost[indices[i]];
+      for (size_t i = indices.size(); i > 0; --i) {
+        HighsInt index = indices[i - 1];
+        double delta = upper[index] * reducedcost[index];
 
-        if (upper[indices[i]] <= 1000.0 && reducedcost[indices[i]] > smallVal &&
+        if (upper[index] <= 1000.0 && reducedcost[index] > smallVal &&
             lambda - delta <= smallVal)
-          cover.push_back(indices[i]);
+          cover.push_back(index);
         else
           lambda -= delta;
       }
@@ -5366,15 +5367,16 @@ HPresolve::Result HPresolve::strengthenInequalities(
     auto updateNonZeros = [&](HighsInt row, HighsCDouble& rhs,
                               HighsInt direction) {
       for (HighsInt i : indices) {
+        assert(Arow[positions[i]] == row);
         double coefdelta = direction * double(reducedcost[i] - maxviolation);
-        HighsInt pos = positions[i];
+        HighsInt col = Acol[positions[i]];
 
         if (complementation[i] == -1) {
-          rhs += coefdelta * model->col_lower_[Acol[pos]];
-          addToMatrix(row, Acol[pos], coefdelta);
+          rhs += coefdelta * model->col_lower_[col];
+          addToMatrix(row, col, coefdelta);
         } else {
-          rhs -= coefdelta * model->col_upper_[Acol[pos]];
-          addToMatrix(row, Acol[pos], -coefdelta);
+          rhs -= coefdelta * model->col_upper_[col];
+          addToMatrix(row, col, -coefdelta);
         }
       }
     };

From 1293a1da25f94199d9fabb818ac002842d563275 Mon Sep 17 00:00:00 2001
From: fwesselm <fwesselm@mathworks.com>
Date: Tue, 18 Feb 2025 11:25:28 +0100
Subject: [PATCH 12/28] Avoid rounding twice (before and in
 HPresolve::changeColLower or HPresolve::changeColUpper)

---
 src/presolve/HPresolve.cpp | 18 ++++++------------
 1 file changed, 6 insertions(+), 12 deletions(-)

diff --git a/src/presolve/HPresolve.cpp b/src/presolve/HPresolve.cpp
index 99e9ddb764..85ce8d1275 100644
--- a/src/presolve/HPresolve.cpp
+++ b/src/presolve/HPresolve.cpp
@@ -404,13 +404,9 @@ bool HPresolve::convertImpliedInteger(HighsInt col, HighsInt row,
       ++rowsizeImplInt[nonzero.index()];
   }
 
-  // round bounds
-  double ceilLower = std::ceil(model->col_lower_[col] - primal_feastol);
-  double floorUpper = std::floor(model->col_upper_[col] + primal_feastol);
-
-  // use tighter bounds
-  if (ceilLower > model->col_lower_[col]) changeColLower(col, ceilLower);
-  if (floorUpper < model->col_upper_[col]) changeColUpper(col, floorUpper);
+  // round and update bounds
+  changeColLower(col, model->col_lower_[col]);
+  changeColUpper(col, model->col_upper_[col]);
   return true;
 }
 
@@ -4033,11 +4029,9 @@ HPresolve::Result HPresolve::initialRowAndColPresolve(
   for (HighsInt col = 0; col != model->num_col_; ++col) {
     if (colDeleted[col]) continue;
     if (model->integrality_[col] != HighsVarType::kContinuous) {
-      double ceilLower = std::ceil(model->col_lower_[col] - primal_feastol);
-      double floorUpper = std::floor(model->col_upper_[col] + primal_feastol);
-
-      if (ceilLower > model->col_lower_[col]) changeColLower(col, ceilLower);
-      if (floorUpper < model->col_upper_[col]) changeColUpper(col, floorUpper);
+      // round and update bounds
+      changeColLower(col, model->col_lower_[col]);
+      changeColUpper(col, model->col_upper_[col]);
     }
     HPRESOLVE_CHECKED_CALL(colPresolve(postsolve_stack, col));
     changedColFlag[col] = false;

From 54c74fa5b664caae77d1a6f6be8ecee88741ac98 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 10:25:59 +0000
Subject: [PATCH 13/28] Now to extend use of
 HighsMipSolverData::setCallbackDataOut

---
 check/TestCallbacks.cpp           | 24 +++++++++++++++++++++++-
 src/highs_bindings.cpp            |  2 ++
 src/highspy/highs.py              |  4 ++++
 src/lp_data/HConst.h              |  3 ++-
 src/lp_data/HighsCallback.cpp     |  3 ++-
 src/lp_data/HighsCallbackStruct.h |  1 +
 src/mip/HighsMipSolver.cpp        |  6 ++++++
 src/mip/HighsMipSolverData.cpp    | 22 ++++++++++++++++++++++
 src/mip/HighsMipSolverData.h      |  3 +++
 9 files changed, 65 insertions(+), 3 deletions(-)

diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp
index d97bc0d271..9509cd0741 100644
--- a/check/TestCallbacks.cpp
+++ b/check/TestCallbacks.cpp
@@ -6,7 +6,7 @@
 #include "catch.hpp"
 #include "lp_data/HighsCallback.h"
 
-const bool dev_run = false;
+const bool dev_run = true;
 
 const double egout_optimal_objective = 568.1007;
 const double egout_objective_target = 610;
@@ -179,6 +179,16 @@ HighsCallbackFunctionType userMipCutPoolCallback =
       }
     };
 
+HighsCallbackFunctionType userkMipUserSolution =
+    [](int callback_type, const std::string& message,
+       const HighsCallbackDataOut* data_out, HighsCallbackDataIn* data_in,
+       void* user_callback_data) {
+      if (dev_run) {
+        printf("userkMipUserSolution: mip_primal_bound = %g)\n",
+               data_out->mip_primal_bound);
+      }
+    };
+
 std::function<void(int, const std::string&, const HighsCallbackDataOut*,
                    HighsCallbackDataIn*, void*)>
     userDataCallback = [](int callback_type, const std::string& message,
@@ -380,3 +390,15 @@ TEST_CASE("highs-callback-mip-cut-pool", "[highs-callback]") {
   highs.startCallback(kCallbackMipGetCutPool);
   highs.run();
 }
+
+TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
+  std::string filename = std::string(HIGHS_DIR) + "/check/instances/flugpl.mps";
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("presolve", kHighsOffString);
+  highs.readModel(filename);
+  //  MipData user_callback_data;
+  highs.setCallback(userkMipUserSolution);  //, p_user_callback_data);
+  highs.startCallback(kCallbackMipUserSolution);
+  highs.run();
+}
diff --git a/src/highs_bindings.cpp b/src/highs_bindings.cpp
index a650071503..5b03d426c6 100644
--- a/src/highs_bindings.cpp
+++ b/src/highs_bindings.cpp
@@ -1579,6 +1579,8 @@ PYBIND11_MODULE(_core, m, py::mod_gil_not_used()) {
              HighsCallbackType::kCallbackMipGetCutPool)
       .value("kCallbackMipDefineLazyConstraints",
              HighsCallbackType::kCallbackMipDefineLazyConstraints)
+      .value("kCallbackMipUserSolution",
+             HighsCallbackType::kCallbackMipUserSolution)
       .value("kCallbackMax", HighsCallbackType::kCallbackMax)
       .value("kNumCallbackType", HighsCallbackType::kNumCallbackType)
       .export_values();
diff --git a/src/highspy/highs.py b/src/highspy/highs.py
index af42594be9..a434c2ed7f 100644
--- a/src/highspy/highs.py
+++ b/src/highspy/highs.py
@@ -1303,6 +1303,10 @@ def cbMipGetCutPool(self):
     def cbMipDefineLazyConstraints(self):
         return self.callbacks[int(cb.HighsCallbackType.kCallbackMipDefineLazyConstraints)]
 
+    @property
+    def cbkMipUserSolution(self):
+        return self.callbacks[int(cb.HighsCallbackType.kCallbackMipUserSolution)]
+
     # callback setters are required for +=/-= syntax
     # e.g., h.cbLogging += my_callback
     @cbLogging.setter
diff --git a/src/lp_data/HConst.h b/src/lp_data/HConst.h
index 178c84332a..aa7d0a745e 100644
--- a/src/lp_data/HConst.h
+++ b/src/lp_data/HConst.h
@@ -225,7 +225,8 @@ enum HighsCallbackType : int {
   kCallbackMipInterrupt,              // 6
   kCallbackMipGetCutPool,             // 7
   kCallbackMipDefineLazyConstraints,  // 8
-  kCallbackMax = kCallbackMipDefineLazyConstraints,
+  kCallbackMipUserSolution,           // 9
+  kCallbackMax = kCallbackMipUserSolution,
   kNumCallbackType
 };
 
diff --git a/src/lp_data/HighsCallback.cpp b/src/lp_data/HighsCallback.cpp
index 99e148aa51..75a491f882 100644
--- a/src/lp_data/HighsCallback.cpp
+++ b/src/lp_data/HighsCallback.cpp
@@ -65,7 +65,8 @@ bool HighsCallback::callbackAction(const int callback_type,
       callback_type == kCallbackMipSolution ||
       callback_type == kCallbackMipLogging ||
       callback_type == kCallbackMipGetCutPool ||
-      callback_type == kCallbackMipDefineLazyConstraints)
+      callback_type == kCallbackMipDefineLazyConstraints ||
+      callback_type == kCallbackMipUserSolution)
     assert(!action);
   return action;
 }
diff --git a/src/lp_data/HighsCallbackStruct.h b/src/lp_data/HighsCallbackStruct.h
index 342193f4a5..dc411f1daf 100644
--- a/src/lp_data/HighsCallbackStruct.h
+++ b/src/lp_data/HighsCallbackStruct.h
@@ -46,6 +46,7 @@ typedef struct {
 
 typedef struct {
   int user_interrupt;
+  double* user_solution;
 } HighsCallbackDataIn;
 
 // Additional callback handling
diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index e1aa17c63d..60a5d86a01 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -242,6 +242,12 @@ void HighsMipSolver::run() {
   double lowerBoundLastCheck = mipdata_->lower_bound;
   analysis_.mipTimerStart(kMipClockSearch);
   while (search.hasNode()) {
+
+    // Possibly look for primal solution from the user
+    if (!submip && callback_->user_callback &&
+	callback_->active[kCallbackMipUserSolution])
+      mipdata_->callbackUserSolution(solution_objective_);
+    
     analysis_.mipTimerStart(kMipClockPerformAging1);
     mipdata_->conflictPool.performAging();
     analysis_.mipTimerStop(kMipClockPerformAging1);
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 868737a629..309edc299c 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -2482,6 +2482,22 @@ void HighsMipSolverData::limitsToBounds(double& dual_bound,
 // incumbent value (mipsolver.solution_objective_) is not right for
 // callback_type = kCallbackMipSolution
 
+void HighsMipSolverData::setCallbackDataOut(const double mipsolver_objective_value) {
+  double dual_bound;
+  double primal_bound;
+  double mip_rel_gap;
+  limitsToBounds(dual_bound, primal_bound, mip_rel_gap);
+  mipsolver.callback_->data_out.running_time = mipsolver.timer_.read();
+  mipsolver.callback_->data_out.objective_function_value =
+    mipsolver_objective_value;
+  mipsolver.callback_->data_out.mip_node_count = mipsolver.mipdata_->num_nodes;
+  mipsolver.callback_->data_out.mip_total_lp_iterations =
+    mipsolver.mipdata_->total_lp_iterations;
+  mipsolver.callback_->data_out.mip_primal_bound = primal_bound;
+  mipsolver.callback_->data_out.mip_dual_bound = dual_bound;
+  mipsolver.callback_->data_out.mip_gap = mip_rel_gap;
+}
+
 bool HighsMipSolverData::interruptFromCallbackWithData(
     const int callback_type, const double mipsolver_objective_value,
     const std::string message) const {
@@ -2504,6 +2520,12 @@ bool HighsMipSolverData::interruptFromCallbackWithData(
   return mipsolver.callback_->callbackAction(callback_type, message);
 }
 
+void HighsMipSolverData::callbackUserSolution(const double mipsolver_objective_value) {
+  setCallbackDataOut(mipsolver_objective_value);
+  printf("HighsMipSolverData::callbackUserSolution() mipsolver_objective_value = %g\n",  mipsolver.callback_->data_out.objective_function_value);
+  
+}
+
 static double possInfRelDiff(const double v0, const double v1,
                              const double den) {
   double rel_diff;
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index 8678751e45..347ce471e5 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -286,9 +286,12 @@ struct HighsMipSolverData {
   bool checkLimits(int64_t nodeOffset = 0) const;
   void limitsToBounds(double& dual_bound, double& primal_bound,
                       double& mip_rel_gap) const;
+  void setCallbackDataOut(const double mipsolver_objective_value);
   bool interruptFromCallbackWithData(const int callback_type,
                                      const double mipsolver_objective_value,
                                      const std::string message = "") const;
+  void callbackUserSolution(const double mipsolver_objective_value);
+  
 };
 
 #endif

From 234b4cc9b58f280fa290109c2eefb87a59e465ea Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 11:10:15 +0000
Subject: [PATCH 14/28] User passing optimal solution via MIP user solultion
 callback

---
 check/TestCallbacks.cpp        | 31 ++++++++++++++++++++++++-------
 src/lp_data/HighsCallback.cpp  |  1 +
 src/mip/HighsMipSolverData.cpp | 26 +++++++++-----------------
 src/mip/HighsMipSolverData.h   |  2 +-
 4 files changed, 35 insertions(+), 25 deletions(-)

diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp
index 9509cd0741..75c5ece748 100644
--- a/check/TestCallbacks.cpp
+++ b/check/TestCallbacks.cpp
@@ -31,6 +31,11 @@ struct MipData {
   HighsVarType* integrality;
 };
 
+struct UserMipSolution {
+  double optimal_objective_value;
+  double* optimal_solution; 
+};
+
 // Callback that saves message for comparison
 HighsCallbackFunctionType myLogCallback =
     [](int callback_type, const std::string& message,
@@ -183,9 +188,15 @@ HighsCallbackFunctionType userkMipUserSolution =
     [](int callback_type, const std::string& message,
        const HighsCallbackDataOut* data_out, HighsCallbackDataIn* data_in,
        void* user_callback_data) {
+      UserMipSolution callback_data = *(static_cast<UserMipSolution*>(user_callback_data));
       if (dev_run) {
-        printf("userkMipUserSolution: mip_primal_bound = %g)\n",
-               data_out->mip_primal_bound);
+	if (data_out->mip_primal_bound > callback_data.optimal_objective_value) {
+	  // If current objective value is not optimal, pass the
+	  // optimal solution as a user solution
+	  printf("userkMipUserSolution: %g = mip_primal_bound > optimal_objective_value = %g\n",
+		 data_out->mip_primal_bound, callback_data.optimal_objective_value);
+	  data_in->user_solution = callback_data.optimal_solution;
+	}
       }
     };
 
@@ -385,8 +396,7 @@ TEST_CASE("highs-callback-mip-cut-pool", "[highs-callback]") {
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
   highs.readModel(filename);
-  //  MipData user_callback_data;
-  highs.setCallback(userMipCutPoolCallback);  //, p_user_callback_data);
+  highs.setCallback(userMipCutPoolCallback);
   highs.startCallback(kCallbackMipGetCutPool);
   highs.run();
 }
@@ -395,10 +405,17 @@ TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
   std::string filename = std::string(HIGHS_DIR) + "/check/instances/flugpl.mps";
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
-  highs.setOptionValue("presolve", kHighsOffString);
   highs.readModel(filename);
-  //  MipData user_callback_data;
-  highs.setCallback(userkMipUserSolution);  //, p_user_callback_data);
+  highs.run();
+  std::vector<double> optimal_solution = highs.getSolution().col_value;
+  UserMipSolution user_callback_data;
+  user_callback_data.optimal_objective_value = highs.getInfo().objective_function_value;
+  user_callback_data.optimal_solution = optimal_solution.data();
+  void* p_user_callback_data = (void*)(&user_callback_data);
+
+  highs.clearSolver();
+  highs.setOptionValue("presolve", kHighsOffString);
+  highs.setCallback(userkMipUserSolution, p_user_callback_data); 
   highs.startCallback(kCallbackMipUserSolution);
   highs.run();
 }
diff --git a/src/lp_data/HighsCallback.cpp b/src/lp_data/HighsCallback.cpp
index 75a491f882..8ab9b9e256 100644
--- a/src/lp_data/HighsCallback.cpp
+++ b/src/lp_data/HighsCallback.cpp
@@ -28,6 +28,7 @@ void HighsCallback::clearHighsCallbackDataOut() {
 
 void HighsCallback::clearHighsCallbackDataIn() {
   this->data_in.user_interrupt = false;
+  this->data_in.user_solution = nullptr;
 }
 
 void HighsCallback::clear() {
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 309edc299c..f3f821a11b 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -2482,7 +2482,7 @@ void HighsMipSolverData::limitsToBounds(double& dual_bound,
 // incumbent value (mipsolver.solution_objective_) is not right for
 // callback_type = kCallbackMipSolution
 
-void HighsMipSolverData::setCallbackDataOut(const double mipsolver_objective_value) {
+void HighsMipSolverData::setCallbackDataOut(const double mipsolver_objective_value) const {
   double dual_bound;
   double primal_bound;
   double mip_rel_gap;
@@ -2503,27 +2503,19 @@ bool HighsMipSolverData::interruptFromCallbackWithData(
     const std::string message) const {
   if (!mipsolver.callback_->callbackActive(callback_type)) return false;
   assert(!mipsolver.submip);
-
-  double dual_bound;
-  double primal_bound;
-  double mip_rel_gap;
-  limitsToBounds(dual_bound, primal_bound, mip_rel_gap);
-  mipsolver.callback_->data_out.running_time = mipsolver.timer_.read();
-  mipsolver.callback_->data_out.objective_function_value =
-      mipsolver_objective_value;
-  mipsolver.callback_->data_out.mip_node_count = mipsolver.mipdata_->num_nodes;
-  mipsolver.callback_->data_out.mip_total_lp_iterations =
-      mipsolver.mipdata_->total_lp_iterations;
-  mipsolver.callback_->data_out.mip_primal_bound = primal_bound;
-  mipsolver.callback_->data_out.mip_dual_bound = dual_bound;
-  mipsolver.callback_->data_out.mip_gap = mip_rel_gap;
+  setCallbackDataOut(mipsolver_objective_value);
   return mipsolver.callback_->callbackAction(callback_type, message);
 }
 
 void HighsMipSolverData::callbackUserSolution(const double mipsolver_objective_value) {
   setCallbackDataOut(mipsolver_objective_value);
-  printf("HighsMipSolverData::callbackUserSolution() mipsolver_objective_value = %g\n",  mipsolver.callback_->data_out.objective_function_value);
-  
+  mipsolver.callback_->clearHighsCallbackDataIn();
+  const bool interrupt = mipsolver.callback_->callbackAction(kCallbackMipUserSolution, "MIP User solution");
+  assert(!interrupt);
+  if (mipsolver.callback_->data_in.user_solution) {
+    printf("HighsMipSolverData::callbackUserSolution() User solution has first value %g\n",
+	   mipsolver.callback_->data_in.user_solution[0]);
+  }
 }
 
 static double possInfRelDiff(const double v0, const double v1,
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index 347ce471e5..f70ac11929 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -286,7 +286,7 @@ struct HighsMipSolverData {
   bool checkLimits(int64_t nodeOffset = 0) const;
   void limitsToBounds(double& dual_bound, double& primal_bound,
                       double& mip_rel_gap) const;
-  void setCallbackDataOut(const double mipsolver_objective_value);
+  void setCallbackDataOut(const double mipsolver_objective_value) const;
   bool interruptFromCallbackWithData(const int callback_type,
                                      const double mipsolver_objective_value,
                                      const std::string message = "") const;

From 336a183ded5ccb5cb6f30d91073b49da49813411 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 12:54:30 +0000
Subject: [PATCH 15/28] Now to refactor
 HighsMipSolverData::transformNewIntegerFeasibleSolution

---
 check/TestCallbacks.cpp        | 5 +++--
 src/mip/HighsMipSolverData.cpp | 7 +++++++
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp
index 75c5ece748..57dfcbbd7b 100644
--- a/check/TestCallbacks.cpp
+++ b/check/TestCallbacks.cpp
@@ -402,7 +402,8 @@ TEST_CASE("highs-callback-mip-cut-pool", "[highs-callback]") {
 }
 
 TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
-  std::string filename = std::string(HIGHS_DIR) + "/check/instances/flugpl.mps";
+  const std::string model = "bell5";
+  const std::string filename = std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
   highs.readModel(filename);
@@ -414,7 +415,7 @@ TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
   void* p_user_callback_data = (void*)(&user_callback_data);
 
   highs.clearSolver();
-  highs.setOptionValue("presolve", kHighsOffString);
+  //  highs.setOptionValue("presolve", kHighsOffString);
   highs.setCallback(userkMipUserSolution, p_user_callback_data); 
   highs.startCallback(kCallbackMipUserSolution);
   highs.run();
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index f3f821a11b..55cb321392 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -2515,6 +2515,13 @@ void HighsMipSolverData::callbackUserSolution(const double mipsolver_objective_v
   if (mipsolver.callback_->data_in.user_solution) {
     printf("HighsMipSolverData::callbackUserSolution() User solution has first value %g\n",
 	   mipsolver.callback_->data_in.user_solution[0]);
+    printf("HighsMipSolverData::callbackUserSolution() original model has %d columns; model has %d columns\n",
+	   int(mipsolver.orig_model_->num_col_),
+	   int(mipsolver.model_->num_col_)); 
+    std::vector<double> user_solution(mipsolver.orig_model_->num_col_);
+    for (HighsInt iCol = 0; iCol < mipsolver.orig_model_->num_col_; iCol++)
+      user_solution[iCol] = mipsolver.callback_->data_in.user_solution[iCol];
+    transformNewIntegerFeasibleSolution(user_solution, false);
   }
 }
 

From e8080ec8431d65537507e48dcc8ec692879a9053 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 14:05:04 +0000
Subject: [PATCH 16/28] Now to use HighsMipSolver::solutionFeasible

---
 src/mip/HighsMipSolver.cpp     | 67 ++++++++++++++++++++++++++++++++++
 src/mip/HighsMipSolver.h       |  8 ++++
 src/mip/HighsMipSolverData.cpp | 11 +++---
 src/mip/HighsMipSolverData.h   |  3 +-
 4 files changed, 81 insertions(+), 8 deletions(-)

diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index 60a5d86a01..72fde7b3a5 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -879,3 +879,70 @@ void HighsMipSolver::callbackGetCutPool() const {
                            &callback_->data_out, &callback_->data_in,
                            callback_->user_callback_data);
 }
+
+bool HighsMipSolver::solutionFeasible(const HighsLp* lp,
+				      const std::vector<double>& col_value,
+				      const std::vector<double>* pass_row_value,
+				      double& bound_violation,
+				      double& row_violation,
+				      double& integrality_violation,
+				      double& obj) {
+  bound_violation = 0;
+  row_violation = 0;
+  integrality_violation = 0;
+  const double mip_feasibility_tolerance = options_mip_->mip_feasibility_tolerance;
+
+  HighsCDouble quad_obj = lp->offset_;
+  assert((HighsInt)col_value.size() == lp->num_col_);
+  for (HighsInt i = 0; i != lp->num_col_; ++i) {
+    const double value = col_value[i];
+    quad_obj += lp->col_cost_[i] * value;
+
+    if (lp->integrality_[i] == HighsVarType::kInteger) {
+      integrality_violation =
+	std::max(fractionality(value), integrality_violation);
+    }
+
+    const double lower = lp->col_lower_[i];
+    const double upper = lp->col_upper_[i];
+    double primal_infeasibility;
+    if (value < lower - mip_feasibility_tolerance) {
+      primal_infeasibility = lower - value;
+    } else if (value > upper + mip_feasibility_tolerance) {
+      primal_infeasibility = value - upper;
+    } else
+      continue;
+    
+    bound_violation = std::max(bound_violation, primal_infeasibility);
+  }
+
+  
+  std::vector<double>row_value;
+  if (!pass_row_value)
+    calculateRowValuesQuad(*lp, col_value, row_value);
+  const double* row_value_p = pass_row_value ? (*pass_row_value).data() : row_value.data();
+  assert(row_value_p);
+  
+  for (HighsInt i = 0; i != lp->num_row_; ++i) {
+    const double value = row_value_p[i];
+    const double lower = lp->row_lower_[i];
+    const double upper = lp->row_upper_[i];
+    
+    double primal_infeasibility;
+    if (value < lower - mip_feasibility_tolerance) {
+      primal_infeasibility = lower - value;
+    } else if (value > upper + mip_feasibility_tolerance) {
+      primal_infeasibility = value - upper;
+    } else
+      continue;
+
+    row_violation = std::max(row_violation, primal_infeasibility);
+  }
+  obj = double(quad_obj);
+ 
+  const bool feasible =
+    bound_violation <= mip_feasibility_tolerance &&
+    integrality_violation <= mip_feasibility_tolerance &&
+    row_violation <= mip_feasibility_tolerance;
+  return feasible;
+}
diff --git a/src/mip/HighsMipSolver.h b/src/mip/HighsMipSolver.h
index 2517b09a85..45cecd2c91 100644
--- a/src/mip/HighsMipSolver.h
+++ b/src/mip/HighsMipSolver.h
@@ -103,6 +103,14 @@ class HighsMipSolver {
   presolve::HighsPostsolveStack getPostsolveStack() const;
 
   void callbackGetCutPool() const;
+  bool solutionFeasible(const HighsLp* lp,
+			const std::vector<double>& col_value,
+			const std::vector<double>* pass_row_value,
+			double& bound_violation,
+			double& row_violation,
+			double& integrality_violation,
+			double& obj);
+
 };
 
 #endif
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 55cb321392..e9efff29bd 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -69,12 +69,9 @@ std::string HighsMipSolverData::solutionSourceToString(
   } else if (solution_source == kSolutionSourceTrivialP) {
     if (code) return "p";
     return "Trivial point";
-    //  } else if (solution_source == kSolutionSourceOpt1) {
-    //    if (code) return "1";
-    //    return "1-opt";
-    //  } else if (solution_source == kSolutionSourceOpt2) {
-    //    if (code) return "2";
-    //    return "2-opt";
+  } else if (solution_source == kSolutionSourceUserSolution) {
+    if (code) return "X";
+    return "User solution";
   } else if (solution_source == kSolutionSourceCleanup) {
     if (code) return " ";
     return "";
@@ -970,6 +967,8 @@ void HighsMipSolverData::runSetup() {
                  "\n");
 }
 
+//double HighsMipSolverData::untransformSolution(}
+
 double HighsMipSolverData::transformNewIntegerFeasibleSolution(
     const std::vector<double>& sol,
     const bool possibly_store_as_new_incumbent) {
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index f70ac11929..30e21234c9 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -57,8 +57,7 @@ enum MipSolutionSource : int {
   kSolutionSourceTrivialL,
   kSolutionSourceTrivialU,
   kSolutionSourceTrivialP,
-  //  kSolutionSourceOpt1,
-  //  kSolutionSourceOpt2,
+  kSolutionSourceUserSolution,
   kSolutionSourceCleanup,
   kSolutionSourceCount
 };

From 19ead6b7597df61e0157e86d76733fa02f0b0f0e Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 15:22:53 +0000
Subject: [PATCH 17/28] HighsMipSolver constructor now uses solutionFeasible;
 now to use this in transformNewIntegerFeasibleSolution

---
 src/mip/HighsMipSolver.cpp | 87 +++++++++++---------------------------
 1 file changed, 25 insertions(+), 62 deletions(-)

diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index 72fde7b3a5..39095c4eaf 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -53,51 +53,8 @@ HighsMipSolver::HighsMipSolver(HighsCallback& callback,
                            integral, feasible);
     assert(valid);
 #endif
-    bound_violation_ = 0;
-    row_violation_ = 0;
-    integrality_violation_ = 0;
-
-    HighsCDouble obj = orig_model_->offset_;
-    assert((HighsInt)solution.col_value.size() == orig_model_->num_col_);
-    for (HighsInt i = 0; i != orig_model_->num_col_; ++i) {
-      const double value = solution.col_value[i];
-      obj += orig_model_->col_cost_[i] * value;
-
-      if (orig_model_->integrality_[i] == HighsVarType::kInteger) {
-        integrality_violation_ =
-            std::max(fractionality(value), integrality_violation_);
-      }
-
-      const double lower = orig_model_->col_lower_[i];
-      const double upper = orig_model_->col_upper_[i];
-      double primal_infeasibility;
-      if (value < lower - options_mip_->mip_feasibility_tolerance) {
-        primal_infeasibility = lower - value;
-      } else if (value > upper + options_mip_->mip_feasibility_tolerance) {
-        primal_infeasibility = value - upper;
-      } else
-        continue;
-
-      bound_violation_ = std::max(bound_violation_, primal_infeasibility);
-    }
-
-    for (HighsInt i = 0; i != orig_model_->num_row_; ++i) {
-      const double value = solution.row_value[i];
-      const double lower = orig_model_->row_lower_[i];
-      const double upper = orig_model_->row_upper_[i];
-
-      double primal_infeasibility;
-      if (value < lower - options_mip_->mip_feasibility_tolerance) {
-        primal_infeasibility = lower - value;
-      } else if (value > upper + options_mip_->mip_feasibility_tolerance) {
-        primal_infeasibility = value - upper;
-      } else
-        continue;
-
-      row_violation_ = std::max(row_violation_, primal_infeasibility);
-    }
-
-    solution_objective_ = double(obj);
+    // Initial solution can be infeasible, but need to set values for violation and objective
+    solutionFeasible(orig_model_, solution.col_value, &solution.row_value, bound_violation_, row_violation_, integrality_violation_, solution_objective_);
     solution_ = solution.col_value;
   }
 }
@@ -917,26 +874,32 @@ bool HighsMipSolver::solutionFeasible(const HighsLp* lp,
   }
 
   
-  std::vector<double>row_value;
-  if (!pass_row_value)
-    calculateRowValuesQuad(*lp, col_value, row_value);
-  const double* row_value_p = pass_row_value ? (*pass_row_value).data() : row_value.data();
-  assert(row_value_p);
+  // Check row feasibility if there are a positive number of rows.
+  //
+  // If there are no rows and pass_row_value is nullptr, then
+  // row_value_p is also nullptr since row_value is not resized
+  if (lp->num_row_ > 0) {
+    std::vector<double>row_value;
+    if (!pass_row_value)
+      calculateRowValuesQuad(*lp, col_value, row_value);
+    const double* row_value_p = pass_row_value ? (*pass_row_value).data() : row_value.data();
+    assert(row_value_p);
   
-  for (HighsInt i = 0; i != lp->num_row_; ++i) {
-    const double value = row_value_p[i];
-    const double lower = lp->row_lower_[i];
-    const double upper = lp->row_upper_[i];
+    for (HighsInt i = 0; i != lp->num_row_; ++i) {
+      const double value = row_value_p[i];
+      const double lower = lp->row_lower_[i];
+      const double upper = lp->row_upper_[i];
     
-    double primal_infeasibility;
-    if (value < lower - mip_feasibility_tolerance) {
-      primal_infeasibility = lower - value;
-    } else if (value > upper + mip_feasibility_tolerance) {
-      primal_infeasibility = value - upper;
-    } else
-      continue;
+      double primal_infeasibility;
+      if (value < lower - mip_feasibility_tolerance) {
+	primal_infeasibility = lower - value;
+      } else if (value > upper + mip_feasibility_tolerance) {
+	primal_infeasibility = value - upper;
+      } else
+	continue;
 
-    row_violation = std::max(row_violation, primal_infeasibility);
+      row_violation = std::max(row_violation, primal_infeasibility);
+    }
   }
   obj = double(quad_obj);
  

From 18c6ac5c837a0808318a817e45473af9fd9d532a Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 15:52:52 +0000
Subject: [PATCH 18/28] Added call to mipsolver.solutionFeasible in
 transformNewIntegerFeasibleSolution

---
 src/mip/HighsMipSolverData.cpp | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index e9efff29bd..ff27ef7827 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -991,7 +991,12 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
   double bound_violation_ = 0;
   double row_violation_ = 0;
   double integrality_violation_ = 0;
+  double temp_mipsolver_objective_value = 0;
 
+  bool temp_feasible =
+    mipsolver.solutionFeasible(mipsolver.orig_model_, solution.col_value, &solution.row_value,
+		     bound_violation_, row_violation_, integrality_violation_,
+		     temp_mipsolver_objective_value);
   // Compute to quad precision the objective function value of the MIP
   // being solved - including the offset, and independent of objective
   // sense

From 3d9d7b7c0d2c23426cb3951a82b0f6d78952a5c5 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 16:17:45 +0000
Subject: [PATCH 19/28]  HighsMipSolver::solutionFeasible now returning
 HighsCDouble obj, and being used in transformNewIntegerFeasibleSolution

---
 check/TestCallbacks.cpp        |  31 +++---
 src/mip/HighsMipSolver.cpp     |  71 +++++++------
 src/mip/HighsMipSolver.h       |  12 +--
 src/mip/HighsMipSolverData.cpp | 180 +++++----------------------------
 src/mip/HighsMipSolverData.h   |   1 -
 5 files changed, 89 insertions(+), 206 deletions(-)

diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp
index 57dfcbbd7b..a2fb5e60da 100644
--- a/check/TestCallbacks.cpp
+++ b/check/TestCallbacks.cpp
@@ -33,7 +33,7 @@ struct MipData {
 
 struct UserMipSolution {
   double optimal_objective_value;
-  double* optimal_solution; 
+  double* optimal_solution;
 };
 
 // Callback that saves message for comparison
@@ -188,15 +188,20 @@ HighsCallbackFunctionType userkMipUserSolution =
     [](int callback_type, const std::string& message,
        const HighsCallbackDataOut* data_out, HighsCallbackDataIn* data_in,
        void* user_callback_data) {
-      UserMipSolution callback_data = *(static_cast<UserMipSolution*>(user_callback_data));
+      UserMipSolution callback_data =
+          *(static_cast<UserMipSolution*>(user_callback_data));
       if (dev_run) {
-	if (data_out->mip_primal_bound > callback_data.optimal_objective_value) {
-	  // If current objective value is not optimal, pass the
-	  // optimal solution as a user solution
-	  printf("userkMipUserSolution: %g = mip_primal_bound > optimal_objective_value = %g\n",
-		 data_out->mip_primal_bound, callback_data.optimal_objective_value);
-	  data_in->user_solution = callback_data.optimal_solution;
-	}
+        if (data_out->mip_primal_bound >
+            callback_data.optimal_objective_value) {
+          // If current objective value is not optimal, pass the
+          // optimal solution as a user solution
+          printf(
+              "userkMipUserSolution: %g = mip_primal_bound > "
+              "optimal_objective_value = %g\n",
+              data_out->mip_primal_bound,
+              callback_data.optimal_objective_value);
+          data_in->user_solution = callback_data.optimal_solution;
+        }
       }
     };
 
@@ -403,20 +408,22 @@ TEST_CASE("highs-callback-mip-cut-pool", "[highs-callback]") {
 
 TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
   const std::string model = "bell5";
-  const std::string filename = std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  const std::string filename =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
   highs.readModel(filename);
   highs.run();
   std::vector<double> optimal_solution = highs.getSolution().col_value;
   UserMipSolution user_callback_data;
-  user_callback_data.optimal_objective_value = highs.getInfo().objective_function_value;
+  user_callback_data.optimal_objective_value =
+      highs.getInfo().objective_function_value;
   user_callback_data.optimal_solution = optimal_solution.data();
   void* p_user_callback_data = (void*)(&user_callback_data);
 
   highs.clearSolver();
   //  highs.setOptionValue("presolve", kHighsOffString);
-  highs.setCallback(userkMipUserSolution, p_user_callback_data); 
+  highs.setCallback(userkMipUserSolution, p_user_callback_data);
   highs.startCallback(kCallbackMipUserSolution);
   highs.run();
 }
diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index 39095c4eaf..1d0fa6af2c 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -53,8 +53,13 @@ HighsMipSolver::HighsMipSolver(HighsCallback& callback,
                            integral, feasible);
     assert(valid);
 #endif
-    // Initial solution can be infeasible, but need to set values for violation and objective
-    solutionFeasible(orig_model_, solution.col_value, &solution.row_value, bound_violation_, row_violation_, integrality_violation_, solution_objective_);
+    // Initial solution can be infeasible, but need to set values for violation
+    // and objective
+    HighsCDouble quad_solution_objective_;
+    solutionFeasible(orig_model_, solution.col_value, &solution.row_value,
+                     bound_violation_, row_violation_, integrality_violation_,
+                     quad_solution_objective_);
+    solution_objective_ = double(quad_solution_objective_);
     solution_ = solution.col_value;
   }
 }
@@ -199,12 +204,11 @@ void HighsMipSolver::run() {
   double lowerBoundLastCheck = mipdata_->lower_bound;
   analysis_.mipTimerStart(kMipClockSearch);
   while (search.hasNode()) {
-
     // Possibly look for primal solution from the user
     if (!submip && callback_->user_callback &&
-	callback_->active[kCallbackMipUserSolution])
+        callback_->active[kCallbackMipUserSolution])
       mipdata_->callbackUserSolution(solution_objective_);
-    
+
     analysis_.mipTimerStart(kMipClockPerformAging1);
     mipdata_->conflictPool.performAging();
     analysis_.mipTimerStop(kMipClockPerformAging1);
@@ -837,27 +841,26 @@ void HighsMipSolver::callbackGetCutPool() const {
                            callback_->user_callback_data);
 }
 
-bool HighsMipSolver::solutionFeasible(const HighsLp* lp,
-				      const std::vector<double>& col_value,
-				      const std::vector<double>* pass_row_value,
-				      double& bound_violation,
-				      double& row_violation,
-				      double& integrality_violation,
-				      double& obj) {
+bool HighsMipSolver::solutionFeasible(
+    const HighsLp* lp, const std::vector<double>& col_value,
+    const std::vector<double>* pass_row_value, double& bound_violation,
+    double& row_violation, double& integrality_violation, HighsCDouble& obj) {
   bound_violation = 0;
   row_violation = 0;
   integrality_violation = 0;
-  const double mip_feasibility_tolerance = options_mip_->mip_feasibility_tolerance;
+  const double mip_feasibility_tolerance =
+      options_mip_->mip_feasibility_tolerance;
+
+  obj = lp->offset_;
 
-  HighsCDouble quad_obj = lp->offset_;
-  assert((HighsInt)col_value.size() == lp->num_col_);
+  if (kAllowDeveloperAssert) assert(HighsInt(col_value.size()) == lp->num_col_);
   for (HighsInt i = 0; i != lp->num_col_; ++i) {
     const double value = col_value[i];
-    quad_obj += lp->col_cost_[i] * value;
+    obj += lp->col_cost_[i] * value;
 
     if (lp->integrality_[i] == HighsVarType::kInteger) {
       integrality_violation =
-	std::max(fractionality(value), integrality_violation);
+          std::max(fractionality(value), integrality_violation);
     }
 
     const double lower = lp->col_lower_[i];
@@ -869,43 +872,45 @@ bool HighsMipSolver::solutionFeasible(const HighsLp* lp,
       primal_infeasibility = value - upper;
     } else
       continue;
-    
+
     bound_violation = std::max(bound_violation, primal_infeasibility);
   }
 
-  
   // Check row feasibility if there are a positive number of rows.
   //
   // If there are no rows and pass_row_value is nullptr, then
   // row_value_p is also nullptr since row_value is not resized
   if (lp->num_row_ > 0) {
-    std::vector<double>row_value;
-    if (!pass_row_value)
+    std::vector<double> row_value;
+    if (pass_row_value) {
+      if (kAllowDeveloperAssert)
+        assert(HighsInt((*pass_row_value).size()) == lp->num_col_);
+    } else {
       calculateRowValuesQuad(*lp, col_value, row_value);
-    const double* row_value_p = pass_row_value ? (*pass_row_value).data() : row_value.data();
+    }
+    const double* row_value_p =
+        pass_row_value ? (*pass_row_value).data() : row_value.data();
     assert(row_value_p);
-  
+
     for (HighsInt i = 0; i != lp->num_row_; ++i) {
       const double value = row_value_p[i];
       const double lower = lp->row_lower_[i];
       const double upper = lp->row_upper_[i];
-    
+
       double primal_infeasibility;
       if (value < lower - mip_feasibility_tolerance) {
-	primal_infeasibility = lower - value;
+        primal_infeasibility = lower - value;
       } else if (value > upper + mip_feasibility_tolerance) {
-	primal_infeasibility = value - upper;
+        primal_infeasibility = value - upper;
       } else
-	continue;
+        continue;
 
       row_violation = std::max(row_violation, primal_infeasibility);
     }
   }
-  obj = double(quad_obj);
- 
-  const bool feasible =
-    bound_violation <= mip_feasibility_tolerance &&
-    integrality_violation <= mip_feasibility_tolerance &&
-    row_violation <= mip_feasibility_tolerance;
+
+  const bool feasible = bound_violation <= mip_feasibility_tolerance &&
+                        integrality_violation <= mip_feasibility_tolerance &&
+                        row_violation <= mip_feasibility_tolerance;
   return feasible;
 }
diff --git a/src/mip/HighsMipSolver.h b/src/mip/HighsMipSolver.h
index 45cecd2c91..e1e9b380ca 100644
--- a/src/mip/HighsMipSolver.h
+++ b/src/mip/HighsMipSolver.h
@@ -103,14 +103,10 @@ class HighsMipSolver {
   presolve::HighsPostsolveStack getPostsolveStack() const;
 
   void callbackGetCutPool() const;
-  bool solutionFeasible(const HighsLp* lp,
-			const std::vector<double>& col_value,
-			const std::vector<double>* pass_row_value,
-			double& bound_violation,
-			double& row_violation,
-			double& integrality_violation,
-			double& obj);
-
+  bool solutionFeasible(const HighsLp* lp, const std::vector<double>& col_value,
+                        const std::vector<double>* pass_row_value,
+                        double& bound_violation, double& row_violation,
+                        double& integrality_violation, HighsCDouble& obj);
 };
 
 #endif
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index ff27ef7827..18bf9c8b19 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -967,7 +967,7 @@ void HighsMipSolverData::runSetup() {
                  "\n");
 }
 
-//double HighsMipSolverData::untransformSolution(}
+// double HighsMipSolverData::untransformSolution(}
 
 double HighsMipSolverData::transformNewIntegerFeasibleSolution(
     const std::vector<double>& sol,
@@ -979,10 +979,8 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
   postSolveStack.undoPrimal(*mipsolver.options_mip_, solution);
   // Determine the row values, as they aren't computed in primal
   // postsolve
-  HighsInt first_check_row =
-      -1;  // mipsolver.mipdata_->presolve.debugGetCheckRow();
   HighsStatus return_status =
-      calculateRowValuesQuad(*mipsolver.orig_model_, solution, first_check_row);
+      calculateRowValuesQuad(*mipsolver.orig_model_, solution);
   if (kAllowDeveloperAssert) assert(return_status == HighsStatus::kOk);
   bool allow_try_again = true;
 try_again:
@@ -991,94 +989,12 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
   double bound_violation_ = 0;
   double row_violation_ = 0;
   double integrality_violation_ = 0;
-  double temp_mipsolver_objective_value = 0;
-
-  bool temp_feasible =
-    mipsolver.solutionFeasible(mipsolver.orig_model_, solution.col_value, &solution.row_value,
-		     bound_violation_, row_violation_, integrality_violation_,
-		     temp_mipsolver_objective_value);
-  // Compute to quad precision the objective function value of the MIP
-  // being solved - including the offset, and independent of objective
-  // sense
-  //
-  HighsCDouble mipsolver_quad_precision_objective_value =
-      mipsolver.orig_model_->offset_;
-  if (kAllowDeveloperAssert)
-    assert((HighsInt)solution.col_value.size() ==
-           mipsolver.orig_model_->num_col_);
-  HighsInt check_col = -1;
-  HighsInt check_int = -1;
-  HighsInt check_row = -1;
-  const bool debug_report = false;
-  for (HighsInt i = 0; i != mipsolver.orig_model_->num_col_; ++i) {
-    const double value = solution.col_value[i];
-    mipsolver_quad_precision_objective_value +=
-        mipsolver.orig_model_->col_cost_[i] * value;
-
-    if (mipsolver.orig_model_->integrality_[i] == HighsVarType::kInteger) {
-      double integrality_infeasibility = fractionality(value);
-      if (integrality_infeasibility >
-          mipsolver.options_mip_->mip_feasibility_tolerance) {
-        if (debug_report)
-          printf("Col %d[%s] value %g has integrality infeasibility %g\n",
-                 int(i), mipsolver.orig_model_->col_names_[i].c_str(), value,
-                 integrality_infeasibility);
-        check_int = i;
-      }
-      integrality_violation_ =
-          std::max(integrality_infeasibility, integrality_violation_);
-    }
-
-    const double lower = mipsolver.orig_model_->col_lower_[i];
-    const double upper = mipsolver.orig_model_->col_upper_[i];
-    double primal_infeasibility = 0;
-    if (value < lower - mipsolver.options_mip_->mip_feasibility_tolerance) {
-      primal_infeasibility = lower - value;
-    } else if (value >
-               upper + mipsolver.options_mip_->mip_feasibility_tolerance) {
-      primal_infeasibility = value - upper;
-    } else
-      continue;
-    if (primal_infeasibility >
-        mipsolver.options_mip_->primal_feasibility_tolerance) {
-      if (debug_report)
-        printf("Col %d[%s] [%g, %g, %g] has infeasibility %g\n", int(i),
-               mipsolver.orig_model_->col_names_[i].c_str(), lower, value,
-               upper, primal_infeasibility);
-      check_col = i;
-    }
-    bound_violation_ = std::max(bound_violation_, primal_infeasibility);
-  }
-
-  for (HighsInt i = 0; i != mipsolver.orig_model_->num_row_; ++i) {
-    const double value = solution.row_value[i];
-    const double lower = mipsolver.orig_model_->row_lower_[i];
-    const double upper = mipsolver.orig_model_->row_upper_[i];
-    double primal_infeasibility;
-    if (value < lower - mipsolver.options_mip_->mip_feasibility_tolerance) {
-      primal_infeasibility = lower - value;
-    } else if (value >
-               upper + mipsolver.options_mip_->mip_feasibility_tolerance) {
-      primal_infeasibility = value - upper;
-    } else
-      continue;
-    if (primal_infeasibility >
-        mipsolver.options_mip_->primal_feasibility_tolerance) {
-      if (debug_report)
-        printf("Row %d[%s] [%g, %g, %g] has infeasibility %g\n", int(i),
-               mipsolver.orig_model_->row_names_[i].c_str(), lower, value,
-               upper, primal_infeasibility);
-      check_row = i;
-    }
-    row_violation_ = std::max(row_violation_, primal_infeasibility);
-  }
-
-  bool feasible =
-      bound_violation_ <= mipsolver.options_mip_->mip_feasibility_tolerance &&
-      integrality_violation_ <=
-          mipsolver.options_mip_->mip_feasibility_tolerance &&
-      row_violation_ <= mipsolver.options_mip_->mip_feasibility_tolerance;
-
+  HighsCDouble mipsolver_quad_objective_value = 0;
+  bool feasible = mipsolver.solutionFeasible(
+      mipsolver.orig_model_, solution.col_value, &solution.row_value,
+      bound_violation_, row_violation_, integrality_violation_,
+      mipsolver_quad_objective_value);
+  double mipsolver_objective_value = double(mipsolver_quad_objective_value);
   if (!feasible && allow_try_again) {
     // printf(
     //     "trying to repair sol that is violated by %.12g bounds, %.12g "
@@ -1124,10 +1040,6 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
     }
   }
 
-  // Get a double precision version of the objective function value of
-  // the MIP being solved
-  const double mipsolver_objective_value =
-      double(mipsolver_quad_precision_objective_value);
   // Possible MIP solution callback
   if (!mipsolver.submip && feasible && mipsolver.callback_->user_callback &&
       mipsolver.callback_->active[kCallbackMipSolution]) {
@@ -1160,51 +1072,11 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
               mipsolver.options_mip_->mip_feasibility_tolerance &&
           mipsolver.row_violation_ <=
               mipsolver.options_mip_->mip_feasibility_tolerance;
-      //    check_col = 37;//mipsolver.mipdata_->presolve.debugGetCheckCol();
-      //    check_row = 37;//mipsolver.mipdata_->presolve.debugGetCheckRow();
-      std::string check_col_data = "";
-      if (check_col >= 0) {
-        check_col_data = " (col " + std::to_string(check_col);
-        if (mipsolver.orig_model_->col_names_.size())
-          check_col_data +=
-              "[" + mipsolver.orig_model_->col_names_[check_col] + "]";
-        check_col_data += ")";
-      }
-      std::string check_int_data = "";
-      if (check_int >= 0) {
-        check_int_data = " (col " + std::to_string(check_int);
-        if (mipsolver.orig_model_->col_names_.size())
-          check_int_data +=
-              "[" + mipsolver.orig_model_->col_names_[check_int] + "]";
-        check_int_data += ")";
-      }
-      std::string check_row_data = "";
-      if (check_row >= 0) {
-        check_row_data = " (row " + std::to_string(check_row);
-        if (mipsolver.orig_model_->row_names_.size())
-          check_row_data +=
-              "[" + mipsolver.orig_model_->row_names_[check_row] + "]";
-        check_row_data += ")";
-      }
       highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kWarning,
                    "Solution with objective %g has untransformed violations: "
-                   "bound = %.4g%s; integrality = %.4g%s; row = %.4g%s\n",
+                   "bound = %.4g; integrality = %.4g; row = %.4g\n",
                    mipsolver_objective_value, bound_violation_,
-                   check_col_data.c_str(), integrality_violation_,
-                   check_int_data.c_str(), row_violation_,
-                   check_row_data.c_str());
-
-      const bool debug_repeat = false;  // true;//
-      if (debug_repeat) {
-        HighsSolution check_solution;
-        check_solution.col_value = sol;
-        check_solution.value_valid = true;
-        postSolveStack.undoPrimal(*mipsolver.options_mip_, check_solution,
-                                  check_col);
-        fflush(stdout);
-        if (kAllowDeveloperAssert) assert(111 == 999);
-      }
-
+                   integrality_violation_, row_violation_);
       if (!currentFeasible) {
         // if the current incumbent is non existent or also not feasible we
         // still store the new one
@@ -1221,11 +1093,9 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
   }
   // return the objective value in the transformed space
   if (mipsolver.orig_model_->sense_ == ObjSense::kMaximize)
-    return -double(mipsolver_quad_precision_objective_value +
-                   mipsolver.model_->offset_);
+    return -double(mipsolver_quad_objective_value + mipsolver.model_->offset_);
 
-  return double(mipsolver_quad_precision_objective_value -
-                mipsolver.model_->offset_);
+  return double(mipsolver_quad_objective_value - mipsolver.model_->offset_);
 }
 
 double HighsMipSolverData::percentageInactiveIntegers() const {
@@ -2486,17 +2356,18 @@ void HighsMipSolverData::limitsToBounds(double& dual_bound,
 // incumbent value (mipsolver.solution_objective_) is not right for
 // callback_type = kCallbackMipSolution
 
-void HighsMipSolverData::setCallbackDataOut(const double mipsolver_objective_value) const {
+void HighsMipSolverData::setCallbackDataOut(
+    const double mipsolver_objective_value) const {
   double dual_bound;
   double primal_bound;
   double mip_rel_gap;
   limitsToBounds(dual_bound, primal_bound, mip_rel_gap);
   mipsolver.callback_->data_out.running_time = mipsolver.timer_.read();
   mipsolver.callback_->data_out.objective_function_value =
-    mipsolver_objective_value;
+      mipsolver_objective_value;
   mipsolver.callback_->data_out.mip_node_count = mipsolver.mipdata_->num_nodes;
   mipsolver.callback_->data_out.mip_total_lp_iterations =
-    mipsolver.mipdata_->total_lp_iterations;
+      mipsolver.mipdata_->total_lp_iterations;
   mipsolver.callback_->data_out.mip_primal_bound = primal_bound;
   mipsolver.callback_->data_out.mip_dual_bound = dual_bound;
   mipsolver.callback_->data_out.mip_gap = mip_rel_gap;
@@ -2511,17 +2382,22 @@ bool HighsMipSolverData::interruptFromCallbackWithData(
   return mipsolver.callback_->callbackAction(callback_type, message);
 }
 
-void HighsMipSolverData::callbackUserSolution(const double mipsolver_objective_value) {
+void HighsMipSolverData::callbackUserSolution(
+    const double mipsolver_objective_value) {
   setCallbackDataOut(mipsolver_objective_value);
   mipsolver.callback_->clearHighsCallbackDataIn();
-  const bool interrupt = mipsolver.callback_->callbackAction(kCallbackMipUserSolution, "MIP User solution");
+  const bool interrupt = mipsolver.callback_->callbackAction(
+      kCallbackMipUserSolution, "MIP User solution");
   assert(!interrupt);
   if (mipsolver.callback_->data_in.user_solution) {
-    printf("HighsMipSolverData::callbackUserSolution() User solution has first value %g\n",
-	   mipsolver.callback_->data_in.user_solution[0]);
-    printf("HighsMipSolverData::callbackUserSolution() original model has %d columns; model has %d columns\n",
-	   int(mipsolver.orig_model_->num_col_),
-	   int(mipsolver.model_->num_col_)); 
+    printf(
+        "HighsMipSolverData::callbackUserSolution() User solution has first "
+        "value %g\n",
+        mipsolver.callback_->data_in.user_solution[0]);
+    printf(
+        "HighsMipSolverData::callbackUserSolution() original model has %d "
+        "columns; model has %d columns\n",
+        int(mipsolver.orig_model_->num_col_), int(mipsolver.model_->num_col_));
     std::vector<double> user_solution(mipsolver.orig_model_->num_col_);
     for (HighsInt iCol = 0; iCol < mipsolver.orig_model_->num_col_; iCol++)
       user_solution[iCol] = mipsolver.callback_->data_in.user_solution[iCol];
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index 30e21234c9..c9846b876b 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -290,7 +290,6 @@ struct HighsMipSolverData {
                                      const double mipsolver_objective_value,
                                      const std::string message = "") const;
   void callbackUserSolution(const double mipsolver_objective_value);
-  
 };
 
 #endif

From 9a5244cd6d2668113da2cad048b330eb51bea54c Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 16:41:06 +0000
Subject: [PATCH 20/28] Now using solutionFeasible in
 HighsMipSolverData::callbackUserSolution

---
 src/mip/HighsMipSolverData.cpp | 24 ++++++++++++++----------
 src/mip/HighsMipSolverData.h   |  1 +
 2 files changed, 15 insertions(+), 10 deletions(-)

diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 18bf9c8b19..43444d2d63 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -967,7 +967,10 @@ void HighsMipSolverData::runSetup() {
                  "\n");
 }
 
-// double HighsMipSolverData::untransformSolution(}
+void HighsMipSolverData::presolveSolution(const std::vector<double>& sol, std::vector<double>& presolved_sol) {
+
+
+}
 
 double HighsMipSolverData::transformNewIntegerFeasibleSolution(
     const std::vector<double>& sol,
@@ -2390,18 +2393,19 @@ void HighsMipSolverData::callbackUserSolution(
       kCallbackMipUserSolution, "MIP User solution");
   assert(!interrupt);
   if (mipsolver.callback_->data_in.user_solution) {
-    printf(
-        "HighsMipSolverData::callbackUserSolution() User solution has first "
-        "value %g\n",
-        mipsolver.callback_->data_in.user_solution[0]);
-    printf(
-        "HighsMipSolverData::callbackUserSolution() original model has %d "
-        "columns; model has %d columns\n",
-        int(mipsolver.orig_model_->num_col_), int(mipsolver.model_->num_col_));
     std::vector<double> user_solution(mipsolver.orig_model_->num_col_);
     for (HighsInt iCol = 0; iCol < mipsolver.orig_model_->num_col_; iCol++)
       user_solution[iCol] = mipsolver.callback_->data_in.user_solution[iCol];
-    transformNewIntegerFeasibleSolution(user_solution, false);
+    double bound_violation_ = 0;
+    double row_violation_ = 0;
+    double integrality_violation_ = 0;
+    HighsCDouble mipsolver_quad_objective_value = 0;
+    const bool feasible = mipsolver.solutionFeasible(
+         mipsolver.orig_model_, user_solution, nullptr,
+         bound_violation_, row_violation_, integrality_violation_,
+         mipsolver_quad_objective_value);
+    printf("HighsMipSolverData::callbackUserSolution() User solution is %s\n",
+	   feasible ? "Feasible" : "Not feasible");
   }
 }
 
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index c9846b876b..a6f6114ce1 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -251,6 +251,7 @@ struct HighsMipSolverData {
   void setupDomainPropagation();
   void saveReportMipSolution(const double new_upper_limit = -kHighsInf);
   void runSetup();
+  void presolveSolution(const std::vector<double>& sol, std::vector<double>& presolved_sol);
   double transformNewIntegerFeasibleSolution(
       const std::vector<double>& sol,
       const bool possibly_store_as_new_incumbent = true);

From d637c5f23b7833a09b7381e324a225cbdd91a0ba Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Wed, 19 Feb 2025 16:43:04 +0000
Subject: [PATCH 21/28] Unfortunately we need to presolve the solution!
 Formatted

---
 src/mip/HighsMipSolverData.cpp | 13 +++++--------
 src/mip/HighsMipSolverData.h   |  3 ++-
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 43444d2d63..802e13a67a 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -967,10 +967,8 @@ void HighsMipSolverData::runSetup() {
                  "\n");
 }
 
-void HighsMipSolverData::presolveSolution(const std::vector<double>& sol, std::vector<double>& presolved_sol) {
-
-
-}
+void HighsMipSolverData::presolveSolution(const std::vector<double>& sol,
+                                          std::vector<double>& presolved_sol) {}
 
 double HighsMipSolverData::transformNewIntegerFeasibleSolution(
     const std::vector<double>& sol,
@@ -2401,11 +2399,10 @@ void HighsMipSolverData::callbackUserSolution(
     double integrality_violation_ = 0;
     HighsCDouble mipsolver_quad_objective_value = 0;
     const bool feasible = mipsolver.solutionFeasible(
-         mipsolver.orig_model_, user_solution, nullptr,
-         bound_violation_, row_violation_, integrality_violation_,
-         mipsolver_quad_objective_value);
+        mipsolver.orig_model_, user_solution, nullptr, bound_violation_,
+        row_violation_, integrality_violation_, mipsolver_quad_objective_value);
     printf("HighsMipSolverData::callbackUserSolution() User solution is %s\n",
-	   feasible ? "Feasible" : "Not feasible");
+           feasible ? "Feasible" : "Not feasible");
   }
 }
 
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index a6f6114ce1..cc87ba8150 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -251,7 +251,8 @@ struct HighsMipSolverData {
   void setupDomainPropagation();
   void saveReportMipSolution(const double new_upper_limit = -kHighsInf);
   void runSetup();
-  void presolveSolution(const std::vector<double>& sol, std::vector<double>& presolved_sol);
+  void presolveSolution(const std::vector<double>& sol,
+                        std::vector<double>& presolved_sol);
   double transformNewIntegerFeasibleSolution(
       const std::vector<double>& sol,
       const bool possibly_store_as_new_incumbent = true);

From 089579c0506e1a704a45ca954b7c6d910fc550f2 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Thu, 20 Feb 2025 11:43:36 +0000
Subject: [PATCH 22/28] Now to add random number generator for user solution
 callback

---
 check/TestCallbacks.cpp        | 45 ++++++++++++++---------
 src/mip/HighsMipSolverData.cpp | 66 +++++++++++++++++++++++++---------
 src/mip/HighsMipSolverData.h   |  5 ++-
 3 files changed, 81 insertions(+), 35 deletions(-)

diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp
index a2fb5e60da..3746665099 100644
--- a/check/TestCallbacks.cpp
+++ b/check/TestCallbacks.cpp
@@ -188,20 +188,20 @@ HighsCallbackFunctionType userkMipUserSolution =
     [](int callback_type, const std::string& message,
        const HighsCallbackDataOut* data_out, HighsCallbackDataIn* data_in,
        void* user_callback_data) {
+      
       UserMipSolution callback_data =
           *(static_cast<UserMipSolution*>(user_callback_data));
-      if (dev_run) {
-        if (data_out->mip_primal_bound >
-            callback_data.optimal_objective_value) {
-          // If current objective value is not optimal, pass the
-          // optimal solution as a user solution
-          printf(
-              "userkMipUserSolution: %g = mip_primal_bound > "
-              "optimal_objective_value = %g\n",
-              data_out->mip_primal_bound,
-              callback_data.optimal_objective_value);
-          data_in->user_solution = callback_data.optimal_solution;
-        }
+      if (data_out->mip_primal_bound >
+	  callback_data.optimal_objective_value) {
+	
+	// If current objective value is not optimal, pass the
+	// optimal solution as a user solution
+	if (dev_run) 
+	  printf("userkMipUserSolution: %g = mip_primal_bound > "
+		 "optimal_objective_value = %g\n",
+		 data_out->mip_primal_bound,
+		 callback_data.optimal_objective_value);
+	data_in->user_solution = callback_data.optimal_solution;
       }
     };
 
@@ -407,23 +407,36 @@ TEST_CASE("highs-callback-mip-cut-pool", "[highs-callback]") {
 }
 
 TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
-  const std::string model = "bell5";
+  //  const std::string model = "flugpl";
+  //  const std::string model = "lseu";
+    const std::string model = "egout";
+  //  const std::string model = "gt2";
+  //  const std::string model = "rgn";
+  //  const std::string model = "bell5";
+  //  const std::string model = "sp150x300d";
+  //  const std::string model = "p0548";
+  // const std::string model = "dcmulti";
   const std::string filename =
       std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("mip_rel_gap", 0);
   highs.readModel(filename);
   highs.run();
   std::vector<double> optimal_solution = highs.getSolution().col_value;
+  double objective_function_value0 = highs.getInfo().objective_function_value;
+  highs.clearSolver();
+  
   UserMipSolution user_callback_data;
-  user_callback_data.optimal_objective_value =
-      highs.getInfo().objective_function_value;
+  user_callback_data.optimal_objective_value = objective_function_value0;
   user_callback_data.optimal_solution = optimal_solution.data();
   void* p_user_callback_data = (void*)(&user_callback_data);
 
-  highs.clearSolver();
   //  highs.setOptionValue("presolve", kHighsOffString);
   highs.setCallback(userkMipUserSolution, p_user_callback_data);
   highs.startCallback(kCallbackMipUserSolution);
   highs.run();
+  double objective_function_value1 = highs.getInfo().objective_function_value;
+  double objective_diff = std::fabs(objective_function_value1 - objective_function_value0)/std::max(1.0, std::fabs(objective_function_value0));
+  REQUIRE(objective_diff < 1e-12);
 }
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 802e13a67a..ff63d6b39a 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -967,9 +967,6 @@ void HighsMipSolverData::runSetup() {
                  "\n");
 }
 
-void HighsMipSolverData::presolveSolution(const std::vector<double>& sol,
-                                          std::vector<double>& presolved_sol) {}
-
 double HighsMipSolverData::transformNewIntegerFeasibleSolution(
     const std::vector<double>& sol,
     const bool possibly_store_as_new_incumbent) {
@@ -1276,12 +1273,14 @@ const std::vector<double>& HighsMipSolverData::getSolution() const {
 
 bool HighsMipSolverData::addIncumbent(const std::vector<double>& sol,
                                       double solobj, const int solution_source,
-                                      const bool print_display_line) {
+                                      const bool print_display_line,
+				      const bool is_user_solution) {
   const bool execute_mip_solution_callback =
-      !mipsolver.submip &&
-      (mipsolver.callback_->user_callback
-           ? mipsolver.callback_->active[kCallbackMipSolution]
-           : false);
+    !is_user_solution &&
+    !mipsolver.submip &&
+    (mipsolver.callback_->user_callback
+     ? mipsolver.callback_->active[kCallbackMipSolution]
+     : false);
   // Determine whether the potential new incumbent should be
   // transformed
   //
@@ -1297,10 +1296,7 @@ bool HighsMipSolverData::addIncumbent(const std::vector<double>& sol,
                                : 0;
 
   if (possibly_store_as_new_incumbent) {
-    // #1463 use pre-computed transformed_solobj
     solobj = transformed_solobj;
-    //    solobj = transformNewIntegerFeasibleSolution(sol);
-
     if (solobj >= upper_bound) return false;
 
     double prev_upper_bound = upper_bound;
@@ -1315,7 +1311,7 @@ bool HighsMipSolverData::addIncumbent(const std::vector<double>& sol,
     incumbent = sol;
     double new_upper_limit = computeNewUpperLimit(solobj, 0.0, 0.0);
 
-    if (!mipsolver.submip) saveReportMipSolution(new_upper_limit);
+    if (!is_user_solution && !mipsolver.submip) saveReportMipSolution(new_upper_limit);
     if (new_upper_limit < upper_limit) {
       ++numImprovingSols;
       upper_limit = new_upper_limit;
@@ -1770,6 +1766,11 @@ void HighsMipSolverData::evaluateRootNode() {
 
   printDisplayLine();
 
+  // Possibly look for primal solution from the user
+  if (!mipsolver.submip && mipsolver.callback_->user_callback &&
+      mipsolver.callback_->active[kCallbackMipUserSolution])
+    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+
   if (firstrootbasis.valid)
     lp.getLpSolver().setBasis(firstrootbasis,
                               "HighsMipSolverData::evaluateRootNode");
@@ -2002,6 +2003,12 @@ void HighsMipSolverData::evaluateRootNode() {
     rootlpsolobj = lp.getObjective();
     lp.setIterationLimit(std::max(10000, int(10 * avgrootlpiters)));
     if (ncuts == 0) break;
+
+    // Possibly look for primal solution from the user
+    if (!mipsolver.submip && mipsolver.callback_->user_callback &&
+	mipsolver.callback_->active[kCallbackMipUserSolution])
+      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+
   }
   mipsolver.analysis_.mipTimerStop(kMipClockSeparation);
   if (mipsolver.analysis_.analyse_mip_time) {
@@ -2053,6 +2060,11 @@ void HighsMipSolverData::evaluateRootNode() {
   }
 
   printDisplayLine();
+  // Possibly look for primal solution from the user
+  if (!mipsolver.submip && mipsolver.callback_->user_callback &&
+      mipsolver.callback_->active[kCallbackMipUserSolution])
+    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+
   // Possible cut extraction callback
   if (!mipsolver.submip && mipsolver.callback_->user_callback &&
       mipsolver.callback_->callbackActive(kCallbackMipGetCutPool))
@@ -2100,6 +2112,11 @@ void HighsMipSolverData::evaluateRootNode() {
       ++nseparounds;
 
       printDisplayLine();
+      // Possibly look for primal solution from the user
+      if (!mipsolver.submip && mipsolver.callback_->user_callback &&
+	  mipsolver.callback_->active[kCallbackMipUserSolution])
+	mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+
     }
 
     if (upper_limit != kHighsInf || mipsolver.submip) break;
@@ -2132,6 +2149,11 @@ void HighsMipSolverData::evaluateRootNode() {
 
     ++nseparounds;
     printDisplayLine();
+    // Possibly look for primal solution from the user
+    if (!mipsolver.submip && mipsolver.callback_->user_callback &&
+	mipsolver.callback_->active[kCallbackMipUserSolution])
+      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+
   }
 
   removeFixedIndices();
@@ -2397,12 +2419,24 @@ void HighsMipSolverData::callbackUserSolution(
     double bound_violation_ = 0;
     double row_violation_ = 0;
     double integrality_violation_ = 0;
-    HighsCDouble mipsolver_quad_objective_value = 0;
+    HighsCDouble user_solution_quad_objective_value = 0;
     const bool feasible = mipsolver.solutionFeasible(
         mipsolver.orig_model_, user_solution, nullptr, bound_violation_,
-        row_violation_, integrality_violation_, mipsolver_quad_objective_value);
-    printf("HighsMipSolverData::callbackUserSolution() User solution is %s\n",
-           feasible ? "Feasible" : "Not feasible");
+        row_violation_, integrality_violation_, user_solution_quad_objective_value);
+    double user_solution_objective_value = double(user_solution_quad_objective_value);
+    if (!feasible) {
+      highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kWarning,
+                   "User-supplied solution has with objective %g has violations: "
+                   "bound = %.4g; integrality = %.4g; row = %.4g\n",
+                   user_solution_objective_value, bound_violation_,
+                   integrality_violation_, row_violation_);
+      return;
+    }
+    std::vector<double> reduced_user_solution;
+    reduced_user_solution = postSolveStack.getReducedPrimalSolution(user_solution);
+    const bool print_display_line = true;
+    const bool is_user_solution = true;
+    addIncumbent(reduced_user_solution, user_solution_objective_value, kSolutionSourceUserSolution, print_display_line, is_user_solution);
   }
 }
 
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index cc87ba8150..8f77337994 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -251,8 +251,6 @@ struct HighsMipSolverData {
   void setupDomainPropagation();
   void saveReportMipSolution(const double new_upper_limit = -kHighsInf);
   void runSetup();
-  void presolveSolution(const std::vector<double>& sol,
-                        std::vector<double>& presolved_sol);
   double transformNewIntegerFeasibleSolution(
       const std::vector<double>& sol,
       const bool possibly_store_as_new_incumbent = true);
@@ -267,7 +265,8 @@ struct HighsMipSolverData {
   void evaluateRootNode();
   bool addIncumbent(const std::vector<double>& sol, double solobj,
                     const int solution_source,
-                    const bool print_display_line = true);
+                    const bool print_display_line = true,
+		    const bool is_user_solution = false);
 
   const std::vector<double>& getSolution() const;
 

From f0dd2579be6e4917049396bb9c93d9dee8c4cb10 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Thu, 20 Feb 2025 14:23:43 +0000
Subject: [PATCH 23/28] All 8 unit tests pass for call
 kUserMipSolutionCallbackOriginAfterSetup

---
 check/TestCallbacks.cpp           | 51 ++++++++++++++++++-------------
 src/lp_data/HighsCallback.cpp     |  9 ++++++
 src/lp_data/HighsCallback.h       | 10 ++++++
 src/lp_data/HighsCallbackStruct.h |  1 +
 src/mip/HighsMipSolver.cpp        |  9 +++++-
 src/mip/HighsMipSolverData.cpp    | 26 ++++++++++------
 src/mip/HighsMipSolverData.h      |  3 +-
 7 files changed, 76 insertions(+), 33 deletions(-)

diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp
index 3746665099..a906e6460f 100644
--- a/check/TestCallbacks.cpp
+++ b/check/TestCallbacks.cpp
@@ -34,6 +34,7 @@ struct MipData {
 struct UserMipSolution {
   double optimal_objective_value;
   double* optimal_solution;
+  HighsInt require_user_solution_callback_origin;
 };
 
 // Callback that saves message for comparison
@@ -407,36 +408,42 @@ TEST_CASE("highs-callback-mip-cut-pool", "[highs-callback]") {
 }
 
 TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
+  const std::vector<std::string> model = {"flugpl", "lseu", "egout", "gt2", "rgn", "bell5", "sp150x300d", "p0548", "dcmulti"};
   //  const std::string model = "flugpl";
   //  const std::string model = "lseu";
-    const std::string model = "egout";
+  //  const std::string model = "egout";
   //  const std::string model = "gt2";
-  //  const std::string model = "rgn";
+  // const std::string model = "rgn";
   //  const std::string model = "bell5";
   //  const std::string model = "sp150x300d";
   //  const std::string model = "p0548";
-  // const std::string model = "dcmulti";
-  const std::string filename =
-      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  //   const std::string model = "dcmulti";
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
   highs.setOptionValue("mip_rel_gap", 0);
-  highs.readModel(filename);
-  highs.run();
-  std::vector<double> optimal_solution = highs.getSolution().col_value;
-  double objective_function_value0 = highs.getInfo().objective_function_value;
-  highs.clearSolver();
+  HighsInt from_model = 1;
+  HighsInt to_model = HighsInt(model.size());
+  for (HighsInt iModel = from_model; iModel < to_model; iModel++) {
+    const std::string filename =
+      std::string(HIGHS_DIR) + "/check/instances/" + model[iModel] + ".mps";
+    highs.readModel(filename);
+    highs.run();
+    std::vector<double> optimal_solution = highs.getSolution().col_value;
+    double objective_function_value0 = highs.getInfo().objective_function_value;
+    highs.clearSolver();
   
-  UserMipSolution user_callback_data;
-  user_callback_data.optimal_objective_value = objective_function_value0;
-  user_callback_data.optimal_solution = optimal_solution.data();
-  void* p_user_callback_data = (void*)(&user_callback_data);
-
-  //  highs.setOptionValue("presolve", kHighsOffString);
-  highs.setCallback(userkMipUserSolution, p_user_callback_data);
-  highs.startCallback(kCallbackMipUserSolution);
-  highs.run();
-  double objective_function_value1 = highs.getInfo().objective_function_value;
-  double objective_diff = std::fabs(objective_function_value1 - objective_function_value0)/std::max(1.0, std::fabs(objective_function_value0));
-  REQUIRE(objective_diff < 1e-12);
+    UserMipSolution user_callback_data;
+    user_callback_data.optimal_objective_value = objective_function_value0;
+    user_callback_data.optimal_solution = optimal_solution.data();
+    void* p_user_callback_data = (void*)(&user_callback_data);
+
+    //  highs.setOptionValue("presolve", kHighsOffString);
+    highs.setCallback(userkMipUserSolution, p_user_callback_data);
+    highs.startCallback(kCallbackMipUserSolution);
+    highs.run();
+    highs.stopCallback(kCallbackMipUserSolution);
+    double objective_function_value1 = highs.getInfo().objective_function_value;
+    double objective_diff = std::fabs(objective_function_value1 - objective_function_value0)/std::max(1.0, std::fabs(objective_function_value0));
+    REQUIRE(objective_diff < 1e-12);
+  }
 }
diff --git a/src/lp_data/HighsCallback.cpp b/src/lp_data/HighsCallback.cpp
index 8ab9b9e256..556f6b38b7 100644
--- a/src/lp_data/HighsCallback.cpp
+++ b/src/lp_data/HighsCallback.cpp
@@ -24,6 +24,15 @@ void HighsCallback::clearHighsCallbackDataOut() {
   this->data_out.mip_dual_bound = -kHighsInf;
   this->data_out.mip_gap = -1;
   this->data_out.mip_solution = nullptr;
+  this->data_out.cutpool_num_col = 0;
+  this->data_out.cutpool_num_cut = 0;
+  this->data_out.cutpool_num_nz = 0;
+  this->data_out.cutpool_start = nullptr;
+  this->data_out.cutpool_index = nullptr;
+  this->data_out.cutpool_value = nullptr;
+  this->data_out.cutpool_lower = nullptr;
+  this->data_out.cutpool_upper = nullptr;
+  this->data_out.user_solution_callback_origin = 0;
 }
 
 void HighsCallback::clearHighsCallbackDataIn() {
diff --git a/src/lp_data/HighsCallback.h b/src/lp_data/HighsCallback.h
index 73d53f3926..ae13c07811 100644
--- a/src/lp_data/HighsCallback.h
+++ b/src/lp_data/HighsCallback.h
@@ -16,6 +16,16 @@
 #include "lp_data/HStruct.h"
 #include "lp_data/HighsCallbackStruct.h"
 
+enum userMipSolutionCallbackOrigin {
+  kUserMipSolutionCallbackOriginAfterSetup = 0,
+  kUserMipSolutionCallbackOriginBeforeDive,
+  kUserMipSolutionCallbackOriginEvaluateRootNode0,
+  kUserMipSolutionCallbackOriginEvaluateRootNode1,
+  kUserMipSolutionCallbackOriginEvaluateRootNode2,
+  kUserMipSolutionCallbackOriginEvaluateRootNode3,
+  kUserMipSolutionCallbackOriginEvaluateRootNode4
+};
+
 using HighsCallbackFunctionType =
     std::function<void(int, const std::string&, const HighsCallbackDataOut*,
                        HighsCallbackDataIn*, void*)>;
diff --git a/src/lp_data/HighsCallbackStruct.h b/src/lp_data/HighsCallbackStruct.h
index dc411f1daf..d7b488598c 100644
--- a/src/lp_data/HighsCallbackStruct.h
+++ b/src/lp_data/HighsCallbackStruct.h
@@ -42,6 +42,7 @@ typedef struct {
   double* cutpool_value;
   double* cutpool_lower;
   double* cutpool_upper;
+  HighsInt user_solution_callback_origin;
 } HighsCallbackDataOut;
 
 typedef struct {
diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index 1d0fa6af2c..d1e4aa0cd8 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -131,6 +131,12 @@ void HighsMipSolver::run() {
       cleanupSolve();
       return;
     }
+    // Possibly look for primal solution from the user
+    if (!submip && callback_->user_callback &&
+	callback_->active[kCallbackMipUserSolution])
+      mipdata_->callbackUserSolution(solution_objective_,
+				     kUserMipSolutionCallbackOriginAfterSetup);
+
     // Apply the trivial heuristics
     analysis_.mipTimerStart(kMipClockTrivialHeuristics);
     HighsModelStatus model_status = mipdata_->trivialHeuristics();
@@ -207,7 +213,8 @@ void HighsMipSolver::run() {
     // Possibly look for primal solution from the user
     if (!submip && callback_->user_callback &&
         callback_->active[kCallbackMipUserSolution])
-      mipdata_->callbackUserSolution(solution_objective_);
+      mipdata_->callbackUserSolution(solution_objective_,
+				     kUserMipSolutionCallbackOriginBeforeDive);
 
     analysis_.mipTimerStart(kMipClockPerformAging1);
     mipdata_->conflictPool.performAging();
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index ff63d6b39a..10c272a4ea 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -1511,6 +1511,7 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) {
     ub = mipsolver.options_mip_->objective_bound;
 
   auto print_lp_iters = convertToPrintString(total_lp_iterations);
+  HighsInt dynamic_constraints_in_lp = lp.numRows() > 0 ? lp.numRows() - lp.getNumModelRows() : 0;
   if (upper_bound != kHighsInf) {
     std::array<char, 22> gap_string = {};
     if (gap >= 9999.)
@@ -1536,7 +1537,7 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) {
         solutionSourceToString(solution_source).c_str(), print_nodes.data(),
         queue_nodes.data(), print_leaves.data(), explored, lb_string.data(),
         ub_string.data(), gap_string.data(), cutpool.getNumCuts(),
-        lp.numRows() - lp.getNumModelRows(), conflictPool.getNumConflicts(),
+        dynamic_constraints_in_lp, conflictPool.getNumConflicts(),
         print_lp_iters.data(), time_string.c_str());
   } else {
     std::array<char, 22> ub_string;
@@ -1557,7 +1558,7 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) {
         solutionSourceToString(solution_source).c_str(), print_nodes.data(),
         queue_nodes.data(), print_leaves.data(), explored, lb_string.data(),
         ub_string.data(), gap, cutpool.getNumCuts(),
-        lp.numRows() - lp.getNumModelRows(), conflictPool.getNumConflicts(),
+        dynamic_constraints_in_lp, conflictPool.getNumConflicts(),
         print_lp_iters.data(), time_string.c_str());
   }
   // Check that limitsToBounds yields the same values for the
@@ -1769,7 +1770,8 @@ void HighsMipSolverData::evaluateRootNode() {
   // Possibly look for primal solution from the user
   if (!mipsolver.submip && mipsolver.callback_->user_callback &&
       mipsolver.callback_->active[kCallbackMipUserSolution])
-    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
+					     kUserMipSolutionCallbackOriginEvaluateRootNode0);
 
   if (firstrootbasis.valid)
     lp.getLpSolver().setBasis(firstrootbasis,
@@ -2007,7 +2009,8 @@ void HighsMipSolverData::evaluateRootNode() {
     // Possibly look for primal solution from the user
     if (!mipsolver.submip && mipsolver.callback_->user_callback &&
 	mipsolver.callback_->active[kCallbackMipUserSolution])
-      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
+					       kUserMipSolutionCallbackOriginEvaluateRootNode1);
 
   }
   mipsolver.analysis_.mipTimerStop(kMipClockSeparation);
@@ -2063,7 +2066,8 @@ void HighsMipSolverData::evaluateRootNode() {
   // Possibly look for primal solution from the user
   if (!mipsolver.submip && mipsolver.callback_->user_callback &&
       mipsolver.callback_->active[kCallbackMipUserSolution])
-    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
+					     kUserMipSolutionCallbackOriginEvaluateRootNode2);
 
   // Possible cut extraction callback
   if (!mipsolver.submip && mipsolver.callback_->user_callback &&
@@ -2115,7 +2119,8 @@ void HighsMipSolverData::evaluateRootNode() {
       // Possibly look for primal solution from the user
       if (!mipsolver.submip && mipsolver.callback_->user_callback &&
 	  mipsolver.callback_->active[kCallbackMipUserSolution])
-	mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+	mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
+						 kUserMipSolutionCallbackOriginEvaluateRootNode3);
 
     }
 
@@ -2152,7 +2157,8 @@ void HighsMipSolverData::evaluateRootNode() {
     // Possibly look for primal solution from the user
     if (!mipsolver.submip && mipsolver.callback_->user_callback &&
 	mipsolver.callback_->active[kCallbackMipUserSolution])
-      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_);
+      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
+					       kUserMipSolutionCallbackOriginEvaluateRootNode4);
 
   }
 
@@ -2405,9 +2411,11 @@ bool HighsMipSolverData::interruptFromCallbackWithData(
   return mipsolver.callback_->callbackAction(callback_type, message);
 }
 
-void HighsMipSolverData::callbackUserSolution(
-    const double mipsolver_objective_value) {
+void HighsMipSolverData::callbackUserSolution(const double mipsolver_objective_value,
+					      const HighsInt user_solution_callback_origin) {
   setCallbackDataOut(mipsolver_objective_value);
+  mipsolver.callback_->data_out.user_solution_callback_origin = user_solution_callback_origin;
+
   mipsolver.callback_->clearHighsCallbackDataIn();
   const bool interrupt = mipsolver.callback_->callbackAction(
       kCallbackMipUserSolution, "MIP User solution");
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index 8f77337994..36a69bd5e4 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -290,7 +290,8 @@ struct HighsMipSolverData {
   bool interruptFromCallbackWithData(const int callback_type,
                                      const double mipsolver_objective_value,
                                      const std::string message = "") const;
-  void callbackUserSolution(const double mipsolver_objective_value);
+  void callbackUserSolution(const double mipsolver_objective_value,
+			    const HighsInt user_solution_callback_origin);
 };
 
 #endif

From 9ac2d41d1269745c4936af1b24c891ee46ce17b5 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Thu, 20 Feb 2025 15:35:54 +0000
Subject: [PATCH 24/28] Added user MIP solution callback, and tested

---
 check/TestCallbacks.cpp        | 59 +++++++++++----------
 check/TestEkk.cpp              |  1 +
 check/TestMipSolver.cpp        |  1 +
 check/TestSpecialLps.cpp       |  2 +
 src/mip/HighsMipSolver.cpp     |  6 +--
 src/mip/HighsMipSolverData.cpp | 94 +++++++++++++++++++---------------
 src/mip/HighsMipSolverData.h   |  4 +-
 7 files changed, 94 insertions(+), 73 deletions(-)

diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp
index a906e6460f..6a23fd51c1 100644
--- a/check/TestCallbacks.cpp
+++ b/check/TestCallbacks.cpp
@@ -1,3 +1,4 @@
+// #include <algorithm>
 #include <cstdio>
 #include <cstring>
 
@@ -6,7 +7,7 @@
 #include "catch.hpp"
 #include "lp_data/HighsCallback.h"
 
-const bool dev_run = true;
+const bool dev_run = false;  // true;//
 
 const double egout_optimal_objective = 568.1007;
 const double egout_objective_target = 610;
@@ -189,20 +190,23 @@ HighsCallbackFunctionType userkMipUserSolution =
     [](int callback_type, const std::string& message,
        const HighsCallbackDataOut* data_out, HighsCallbackDataIn* data_in,
        void* user_callback_data) {
-      
       UserMipSolution callback_data =
           *(static_cast<UserMipSolution*>(user_callback_data));
-      if (data_out->mip_primal_bound >
-	  callback_data.optimal_objective_value) {
-	
-	// If current objective value is not optimal, pass the
-	// optimal solution as a user solution
-	if (dev_run) 
-	  printf("userkMipUserSolution: %g = mip_primal_bound > "
-		 "optimal_objective_value = %g\n",
-		 data_out->mip_primal_bound,
-		 callback_data.optimal_objective_value);
-	data_in->user_solution = callback_data.optimal_solution;
+      if (data_out->user_solution_callback_origin ==
+          callback_data.require_user_solution_callback_origin) {
+        if (data_out->mip_primal_bound >
+            callback_data.optimal_objective_value) {
+          // If current objective value is not optimal, pass the
+          // optimal solution as a user solution
+          if (dev_run)
+            printf(
+                "userkMipUserSolution: origin = %d; %g = mip_primal_bound > "
+                "optimal_objective_value = %g\n",
+                int(data_out->user_solution_callback_origin),
+                data_out->mip_primal_bound,
+                callback_data.optimal_objective_value);
+          data_in->user_solution = callback_data.optimal_solution;
+        }
       }
     };
 
@@ -408,33 +412,32 @@ TEST_CASE("highs-callback-mip-cut-pool", "[highs-callback]") {
 }
 
 TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
-  const std::vector<std::string> model = {"flugpl", "lseu", "egout", "gt2", "rgn", "bell5", "sp150x300d", "p0548", "dcmulti"};
-  //  const std::string model = "flugpl";
-  //  const std::string model = "lseu";
-  //  const std::string model = "egout";
-  //  const std::string model = "gt2";
-  // const std::string model = "rgn";
-  //  const std::string model = "bell5";
-  //  const std::string model = "sp150x300d";
-  //  const std::string model = "p0548";
-  //   const std::string model = "dcmulti";
+  //  const std::vector<std::string> model = {"rgn", "flugpl", "gt2", "egout",
+  //  "bell5", "lseu", "sp150x300d"};//, "p0548", "dcmulti"}; const
+  //  std::vector<HighsInt> require_origin = {0, 1, 2, 3, 4, 5, 6};
+  const std::vector<std::string> model = {"p0548", "flugpl", "gt2", "egout",
+                                          "sp150x300d"};
+  const std::vector<HighsInt> require_origin = {0, 1, 2, 3, 4};  //, 4, 5, 6};
+  assert(model.size() == require_origin.size());
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
   highs.setOptionValue("mip_rel_gap", 0);
-  HighsInt from_model = 1;
+  HighsInt from_model = 0;
   HighsInt to_model = HighsInt(model.size());
   for (HighsInt iModel = from_model; iModel < to_model; iModel++) {
     const std::string filename =
-      std::string(HIGHS_DIR) + "/check/instances/" + model[iModel] + ".mps";
+        std::string(HIGHS_DIR) + "/check/instances/" + model[iModel] + ".mps";
     highs.readModel(filename);
     highs.run();
     std::vector<double> optimal_solution = highs.getSolution().col_value;
     double objective_function_value0 = highs.getInfo().objective_function_value;
     highs.clearSolver();
-  
+
     UserMipSolution user_callback_data;
     user_callback_data.optimal_objective_value = objective_function_value0;
     user_callback_data.optimal_solution = optimal_solution.data();
+    user_callback_data.require_user_solution_callback_origin =
+        require_origin[iModel];
     void* p_user_callback_data = (void*)(&user_callback_data);
 
     //  highs.setOptionValue("presolve", kHighsOffString);
@@ -443,7 +446,9 @@ TEST_CASE("highs-callback-mip-user-solution", "[highs-callback]") {
     highs.run();
     highs.stopCallback(kCallbackMipUserSolution);
     double objective_function_value1 = highs.getInfo().objective_function_value;
-    double objective_diff = std::fabs(objective_function_value1 - objective_function_value0)/std::max(1.0, std::fabs(objective_function_value0));
+    double objective_diff =
+        std::fabs(objective_function_value1 - objective_function_value0) /
+        std::max(1.0, std::fabs(objective_function_value0));
     REQUIRE(objective_diff < 1e-12);
   }
 }
diff --git a/check/TestEkk.cpp b/check/TestEkk.cpp
index 02a15eb47f..d5566b02ed 100644
--- a/check/TestEkk.cpp
+++ b/check/TestEkk.cpp
@@ -30,6 +30,7 @@ void ekk_solve(Highs& highs, std::string presolve,
   }
 
   REQUIRE(highs.resetOptions() == HighsStatus::kOk);
+  highs.setOptionValue("output_flag", dev_run);
 }
 
 void ekk_distillation(Highs& highs) {
diff --git a/check/TestMipSolver.cpp b/check/TestMipSolver.cpp
index 080f3d124b..6f46ab5cd1 100644
--- a/check/TestMipSolver.cpp
+++ b/check/TestMipSolver.cpp
@@ -800,6 +800,7 @@ void rowlessMIP(Highs& highs) {
 TEST_CASE("issue-2122", "[highs_test_mip_solver]") {
   std::string filename = std::string(HIGHS_DIR) + "/check/instances/2122.lp";
   Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
   highs.setOptionValue("mip_rel_gap", 0);
   highs.setOptionValue("mip_abs_gap", 0);
   highs.readModel(filename);
diff --git a/check/TestSpecialLps.cpp b/check/TestSpecialLps.cpp
index 66d99ed304..32ddd4fe3a 100644
--- a/check/TestSpecialLps.cpp
+++ b/check/TestSpecialLps.cpp
@@ -47,6 +47,7 @@ void solve(Highs& highs, std::string presolve, std::string solver,
     REQUIRE(iteration_count == require_iteration_count);
   }
   REQUIRE(highs.resetOptions() == HighsStatus::kOk);
+  highs.setOptionValue("output_flag", dev_run);
 }
 
 void distillation(Highs& highs) {
@@ -591,6 +592,7 @@ void singularStartingBasis(Highs& highs) {
                                     optimal_objective, dev_run));
 
   REQUIRE(highs.resetOptions() == HighsStatus::kOk);
+  highs.setOptionValue("output_flag", dev_run);
 
   special_lps.reportSolution(highs, dev_run);
 }
diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index d1e4aa0cd8..cb88c8e5d5 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -133,9 +133,9 @@ void HighsMipSolver::run() {
     }
     // Possibly look for primal solution from the user
     if (!submip && callback_->user_callback &&
-	callback_->active[kCallbackMipUserSolution])
+        callback_->active[kCallbackMipUserSolution])
       mipdata_->callbackUserSolution(solution_objective_,
-				     kUserMipSolutionCallbackOriginAfterSetup);
+                                     kUserMipSolutionCallbackOriginAfterSetup);
 
     // Apply the trivial heuristics
     analysis_.mipTimerStart(kMipClockTrivialHeuristics);
@@ -214,7 +214,7 @@ void HighsMipSolver::run() {
     if (!submip && callback_->user_callback &&
         callback_->active[kCallbackMipUserSolution])
       mipdata_->callbackUserSolution(solution_objective_,
-				     kUserMipSolutionCallbackOriginBeforeDive);
+                                     kUserMipSolutionCallbackOriginBeforeDive);
 
     analysis_.mipTimerStart(kMipClockPerformAging1);
     mipdata_->conflictPool.performAging();
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 10c272a4ea..79b81262fd 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -1274,13 +1274,12 @@ const std::vector<double>& HighsMipSolverData::getSolution() const {
 bool HighsMipSolverData::addIncumbent(const std::vector<double>& sol,
                                       double solobj, const int solution_source,
                                       const bool print_display_line,
-				      const bool is_user_solution) {
+                                      const bool is_user_solution) {
   const bool execute_mip_solution_callback =
-    !is_user_solution &&
-    !mipsolver.submip &&
-    (mipsolver.callback_->user_callback
-     ? mipsolver.callback_->active[kCallbackMipSolution]
-     : false);
+      !is_user_solution && !mipsolver.submip &&
+      (mipsolver.callback_->user_callback
+           ? mipsolver.callback_->active[kCallbackMipSolution]
+           : false);
   // Determine whether the potential new incumbent should be
   // transformed
   //
@@ -1311,7 +1310,8 @@ bool HighsMipSolverData::addIncumbent(const std::vector<double>& sol,
     incumbent = sol;
     double new_upper_limit = computeNewUpperLimit(solobj, 0.0, 0.0);
 
-    if (!is_user_solution && !mipsolver.submip) saveReportMipSolution(new_upper_limit);
+    if (!is_user_solution && !mipsolver.submip)
+      saveReportMipSolution(new_upper_limit);
     if (new_upper_limit < upper_limit) {
       ++numImprovingSols;
       upper_limit = new_upper_limit;
@@ -1511,7 +1511,8 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) {
     ub = mipsolver.options_mip_->objective_bound;
 
   auto print_lp_iters = convertToPrintString(total_lp_iterations);
-  HighsInt dynamic_constraints_in_lp = lp.numRows() > 0 ? lp.numRows() - lp.getNumModelRows() : 0;
+  HighsInt dynamic_constraints_in_lp =
+      lp.numRows() > 0 ? lp.numRows() - lp.getNumModelRows() : 0;
   if (upper_bound != kHighsInf) {
     std::array<char, 22> gap_string = {};
     if (gap >= 9999.)
@@ -1557,9 +1558,9 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) {
         // clang-format on
         solutionSourceToString(solution_source).c_str(), print_nodes.data(),
         queue_nodes.data(), print_leaves.data(), explored, lb_string.data(),
-        ub_string.data(), gap, cutpool.getNumCuts(),
-        dynamic_constraints_in_lp, conflictPool.getNumConflicts(),
-        print_lp_iters.data(), time_string.c_str());
+        ub_string.data(), gap, cutpool.getNumCuts(), dynamic_constraints_in_lp,
+        conflictPool.getNumConflicts(), print_lp_iters.data(),
+        time_string.c_str());
   }
   // Check that limitsToBounds yields the same values for the
   // dual_bound, primal_bound (modulo optimization sense) and
@@ -1770,8 +1771,9 @@ void HighsMipSolverData::evaluateRootNode() {
   // Possibly look for primal solution from the user
   if (!mipsolver.submip && mipsolver.callback_->user_callback &&
       mipsolver.callback_->active[kCallbackMipUserSolution])
-    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
-					     kUserMipSolutionCallbackOriginEvaluateRootNode0);
+    mipsolver.mipdata_->callbackUserSolution(
+        mipsolver.solution_objective_,
+        kUserMipSolutionCallbackOriginEvaluateRootNode0);
 
   if (firstrootbasis.valid)
     lp.getLpSolver().setBasis(firstrootbasis,
@@ -2008,10 +2010,10 @@ void HighsMipSolverData::evaluateRootNode() {
 
     // Possibly look for primal solution from the user
     if (!mipsolver.submip && mipsolver.callback_->user_callback &&
-	mipsolver.callback_->active[kCallbackMipUserSolution])
-      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
-					       kUserMipSolutionCallbackOriginEvaluateRootNode1);
-
+        mipsolver.callback_->active[kCallbackMipUserSolution])
+      mipsolver.mipdata_->callbackUserSolution(
+          mipsolver.solution_objective_,
+          kUserMipSolutionCallbackOriginEvaluateRootNode1);
   }
   mipsolver.analysis_.mipTimerStop(kMipClockSeparation);
   if (mipsolver.analysis_.analyse_mip_time) {
@@ -2066,8 +2068,9 @@ void HighsMipSolverData::evaluateRootNode() {
   // Possibly look for primal solution from the user
   if (!mipsolver.submip && mipsolver.callback_->user_callback &&
       mipsolver.callback_->active[kCallbackMipUserSolution])
-    mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
-					     kUserMipSolutionCallbackOriginEvaluateRootNode2);
+    mipsolver.mipdata_->callbackUserSolution(
+        mipsolver.solution_objective_,
+        kUserMipSolutionCallbackOriginEvaluateRootNode2);
 
   // Possible cut extraction callback
   if (!mipsolver.submip && mipsolver.callback_->user_callback &&
@@ -2118,10 +2121,10 @@ void HighsMipSolverData::evaluateRootNode() {
       printDisplayLine();
       // Possibly look for primal solution from the user
       if (!mipsolver.submip && mipsolver.callback_->user_callback &&
-	  mipsolver.callback_->active[kCallbackMipUserSolution])
-	mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
-						 kUserMipSolutionCallbackOriginEvaluateRootNode3);
-
+          mipsolver.callback_->active[kCallbackMipUserSolution])
+        mipsolver.mipdata_->callbackUserSolution(
+            mipsolver.solution_objective_,
+            kUserMipSolutionCallbackOriginEvaluateRootNode3);
     }
 
     if (upper_limit != kHighsInf || mipsolver.submip) break;
@@ -2154,14 +2157,15 @@ void HighsMipSolverData::evaluateRootNode() {
 
     ++nseparounds;
     printDisplayLine();
-    // Possibly look for primal solution from the user
-    if (!mipsolver.submip && mipsolver.callback_->user_callback &&
-	mipsolver.callback_->active[kCallbackMipUserSolution])
-      mipsolver.mipdata_->callbackUserSolution(mipsolver.solution_objective_,
-					       kUserMipSolutionCallbackOriginEvaluateRootNode4);
-
   }
 
+  // Possibly look for primal solution from the user
+  if (!mipsolver.submip && mipsolver.callback_->user_callback &&
+      mipsolver.callback_->active[kCallbackMipUserSolution])
+    mipsolver.mipdata_->callbackUserSolution(
+        mipsolver.solution_objective_,
+        kUserMipSolutionCallbackOriginEvaluateRootNode4);
+
   removeFixedIndices();
   if (lp.getLpSolver().getBasis().valid) lp.removeObsoleteRows();
   rootlpsolobj = lp.getObjective();
@@ -2411,10 +2415,12 @@ bool HighsMipSolverData::interruptFromCallbackWithData(
   return mipsolver.callback_->callbackAction(callback_type, message);
 }
 
-void HighsMipSolverData::callbackUserSolution(const double mipsolver_objective_value,
-					      const HighsInt user_solution_callback_origin) {
+void HighsMipSolverData::callbackUserSolution(
+    const double mipsolver_objective_value,
+    const HighsInt user_solution_callback_origin) {
   setCallbackDataOut(mipsolver_objective_value);
-  mipsolver.callback_->data_out.user_solution_callback_origin = user_solution_callback_origin;
+  mipsolver.callback_->data_out.user_solution_callback_origin =
+      user_solution_callback_origin;
 
   mipsolver.callback_->clearHighsCallbackDataIn();
   const bool interrupt = mipsolver.callback_->callbackAction(
@@ -2430,21 +2436,27 @@ void HighsMipSolverData::callbackUserSolution(const double mipsolver_objective_v
     HighsCDouble user_solution_quad_objective_value = 0;
     const bool feasible = mipsolver.solutionFeasible(
         mipsolver.orig_model_, user_solution, nullptr, bound_violation_,
-        row_violation_, integrality_violation_, user_solution_quad_objective_value);
-    double user_solution_objective_value = double(user_solution_quad_objective_value);
+        row_violation_, integrality_violation_,
+        user_solution_quad_objective_value);
+    double user_solution_objective_value =
+        double(user_solution_quad_objective_value);
     if (!feasible) {
-      highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kWarning,
-                   "User-supplied solution has with objective %g has violations: "
-                   "bound = %.4g; integrality = %.4g; row = %.4g\n",
-                   user_solution_objective_value, bound_violation_,
-                   integrality_violation_, row_violation_);
+      highsLogUser(
+          mipsolver.options_mip_->log_options, HighsLogType::kWarning,
+          "User-supplied solution has with objective %g has violations: "
+          "bound = %.4g; integrality = %.4g; row = %.4g\n",
+          user_solution_objective_value, bound_violation_,
+          integrality_violation_, row_violation_);
       return;
     }
     std::vector<double> reduced_user_solution;
-    reduced_user_solution = postSolveStack.getReducedPrimalSolution(user_solution);
+    reduced_user_solution =
+        postSolveStack.getReducedPrimalSolution(user_solution);
     const bool print_display_line = true;
     const bool is_user_solution = true;
-    addIncumbent(reduced_user_solution, user_solution_objective_value, kSolutionSourceUserSolution, print_display_line, is_user_solution);
+    addIncumbent(reduced_user_solution, user_solution_objective_value,
+                 kSolutionSourceUserSolution, print_display_line,
+                 is_user_solution);
   }
 }
 
diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h
index 36a69bd5e4..b4cd86d4bf 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -266,7 +266,7 @@ struct HighsMipSolverData {
   bool addIncumbent(const std::vector<double>& sol, double solobj,
                     const int solution_source,
                     const bool print_display_line = true,
-		    const bool is_user_solution = false);
+                    const bool is_user_solution = false);
 
   const std::vector<double>& getSolution() const;
 
@@ -291,7 +291,7 @@ struct HighsMipSolverData {
                                      const double mipsolver_objective_value,
                                      const std::string message = "") const;
   void callbackUserSolution(const double mipsolver_objective_value,
-			    const HighsInt user_solution_callback_origin);
+                            const HighsInt user_solution_callback_origin);
 };
 
 #endif

From c6f82516b1434ecaf7837faec44a6e578699e2d1 Mon Sep 17 00:00:00 2001
From: jajhall <jajhall@ed.ac.uk>
Date: Thu, 20 Feb 2025 15:47:41 +0000
Subject: [PATCH 25/28] Updated FEATURES.md

---
 FEATURES.md | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/FEATURES.md b/FEATURES.md
index f6ff74c98c..b3d184e3e8 100644
--- a/FEATURES.md
+++ b/FEATURES.md
@@ -1,5 +1,7 @@
 ## Build changes
 
+Added code coverage report
+
 ## Code changes
 
 Any LP offset is communicated to the IPM solver, and used in logging and primal/dual objective calculations. 
@@ -18,4 +20,14 @@ Primal and dual residual tolerances - applied following IPM or PDLP solution - n
 
 Highs::getCols (Highs::getRows) now runs in linear time if the internal constraint matrix is stored column-wise (row-wise). Added ensureColwise/Rowwise to the Highs class, the C API and highspy so that users can set the internal constraint matrix storage orientation
 
-When columns and rows are deleted from the incumbent LP after a basic solution has been found, HiGHS no longer invalidates the basis. Now it maintains the basic and nonbasic status of the remaining variables and constraints. When the model is re-solved, this information is used to construct a starting basis.
\ No newline at end of file
+When columns and rows are deleted from the incumbent LP after a basic solution has been found, HiGHS no longer invalidates the basis. Now it maintains the basic and nonbasic status of the remaining variables and constraints. When the model is re-solved, this information is used to construct a starting basis.
+
+Fixed bugs in presolve
+
+When running from the command lin, changes to default option values are reported
+
+Added callback to allow users to supply integer feasible solutions to the MIP solver during execution
+
+
+
+

From 05db5c231704d3b4f5ea840f5fcaff7b8a9807b8 Mon Sep 17 00:00:00 2001
From: JAJHall <jajhall@ed.ac.uk>
Date: Sun, 23 Feb 2025 13:39:39 +0000
Subject: [PATCH 26/28] Corrected option message for log_dev_level

---
 src/lp_data/HighsOptions.h | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h
index 18c919869f..f19d4acd8c 100644
--- a/src/lp_data/HighsOptions.h
+++ b/src/lp_data/HighsOptions.h
@@ -1157,11 +1157,12 @@ class HighsOptions : public HighsOptionsStruct {
     // Advanced options
     advanced = true;
 
-    record_int = new OptionRecordInt(
-        "log_dev_level",
-        "Output development messages: 0 => none; 1 => info; 2 => verbose",
-        advanced, &log_dev_level, kHighsLogDevLevelMin, kHighsLogDevLevelNone,
-        kHighsLogDevLevelMax);
+    record_int =
+        new OptionRecordInt("log_dev_level",
+                            "Output development messages: 0 => none; 1 => "
+                            "info; 2 => detailed; 3 => verbose",
+                            advanced, &log_dev_level, kHighsLogDevLevelMin,
+                            kHighsLogDevLevelNone, kHighsLogDevLevelMax);
     records.push_back(record_int);
 
     record_bool = new OptionRecordBool("log_githash", "Log the githash",

From 7a77f519ad28814757a7022f188d0a043226c41e Mon Sep 17 00:00:00 2001
From: JAJHall <jajhall@ed.ac.uk>
Date: Tue, 25 Feb 2025 09:42:07 +0000
Subject: [PATCH 27/28] Added clocks from parallel-MIP branch

---
 src/mip/HighsMipAnalysis.cpp |  1 +
 src/mip/HighsMipSolver.cpp   | 32 ++++++++++++++++++-----
 src/mip/MipTimer.h           | 49 ++++++++++++++++++++++++++++++------
 3 files changed, 69 insertions(+), 13 deletions(-)

diff --git a/src/mip/HighsMipAnalysis.cpp b/src/mip/HighsMipAnalysis.cpp
index 7a6ea2d486..e1b9f4ad2d 100644
--- a/src/mip/HighsMipAnalysis.cpp
+++ b/src/mip/HighsMipAnalysis.cpp
@@ -147,6 +147,7 @@ void HighsMipAnalysis::reportMipTimer() {
   mip_timer.reportMipPresolveClock(mip_clocks);
   mip_timer.reportMipSearchClock(mip_clocks);
   mip_timer.reportMipDiveClock(mip_clocks);
+  mip_timer.reportMipNodeSearchClock(mip_clocks);
   mip_timer.reportMipPrimalHeuristicsClock(mip_clocks);
   mip_timer.reportMipEvaluateRootNodeClock(mip_clocks);
   mip_timer.reportMipSeparationClock(mip_clocks);
diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index cb88c8e5d5..b72c4f7ec6 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -140,6 +140,7 @@ void HighsMipSolver::run() {
     // Apply the trivial heuristics
     analysis_.mipTimerStart(kMipClockTrivialHeuristics);
     HighsModelStatus model_status = mipdata_->trivialHeuristics();
+    analysis_.mipTimerStop(kMipClockTrivialHeuristics);
     if (modelstatus_ == HighsModelStatus::kNotset &&
         model_status == HighsModelStatus::kInfeasible) {
       // trivialHeuristics can spot trivial infeasibility, so act on it
@@ -237,10 +238,10 @@ void HighsMipSolver::run() {
     while (true) {
       // Possibly apply primal heuristics
       if (considerHeuristics && mipdata_->moreHeuristicsAllowed()) {
-        analysis_.mipTimerStart(kMipClockEvaluateNode);
+        analysis_.mipTimerStart(kMipClockEvaluateNode0);
         const HighsSearch::NodeResult evaluate_node_result =
             search.evaluateNode();
-        analysis_.mipTimerStop(kMipClockEvaluateNode);
+        analysis_.mipTimerStop(kMipClockEvaluateNode0);
 
         if (evaluate_node_result == HighsSearch::NodeResult::kSubOptimal) break;
 
@@ -321,9 +322,9 @@ void HighsMipSolver::run() {
     }  // while (true)
     analysis_.mipTimerStop(kMipClockDive);
 
-    analysis_.mipTimerStart(kMipClockOpenNodesToQueue);
+    analysis_.mipTimerStart(kMipClockOpenNodesToQueue0);
     search.openNodesToQueue(mipdata_->nodequeue);
-    analysis_.mipTimerStop(kMipClockOpenNodesToQueue);
+    analysis_.mipTimerStop(kMipClockOpenNodesToQueue0);
 
     search.flushStatistics();
 
@@ -525,13 +526,23 @@ void HighsMipSolver::run() {
       // we evaluate the node directly here instead of performing a dive
       // because we first want to check if the node is not fathomed due to
       // new global information before we perform separation rounds for the node
-      if (search.evaluateNode() == HighsSearch::NodeResult::kSubOptimal)
-        search.currentNodeToQueue(mipdata_->nodequeue);
+      analysis_.mipTimerStart(kMipClockEvaluateNode1);
+      const HighsSearch::NodeResult evaluate_node_result =
+	search.evaluateNode();
+      analysis_.mipTimerStop(kMipClockEvaluateNode1);
+      if (evaluate_node_result == HighsSearch::NodeResult::kSubOptimal) {
+	analysis_.mipTimerStart(kMipClockCurrentNodeToQueue);
+	search.currentNodeToQueue(mipdata_->nodequeue);
+	analysis_.mipTimerStop(kMipClockCurrentNodeToQueue);
+      }
 
       // if the node was pruned we remove it from the search and install the
       // next node from the queue
+        analysis_.mipTimerStart(kMipClockNodePrunedLoop);
       if (search.currentNodePruned()) {
+	analysis_.mipTimerStart(kMipClockSearchBacktrack);
         search.backtrack();
+	analysis_.mipTimerStop(kMipClockSearchBacktrack);
         ++mipdata_->num_leaves;
         ++mipdata_->num_nodes;
         search.flushStatistics();
@@ -553,6 +564,7 @@ void HighsMipSolver::run() {
             mipdata_->updatePrimalDualIntegral(
                 prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
                 mipdata_->upper_bound);
+	  analysis_.mipTimerStop(kMipClockNodePrunedLoop);
           break;
         }
 
@@ -588,19 +600,27 @@ void HighsMipSolver::run() {
           mipdata_->domain.clearChangedCols();
           mipdata_->removeFixedIndices();
         }
+	analysis_.mipTimerStop(kMipClockStoreBasis);
 
+	analysis_.mipTimerStop(kMipClockNodePrunedLoop);
         continue;
       }
+      analysis_.mipTimerStop(kMipClockNodePrunedLoop);
 
       // the node is still not fathomed, so perform separation
+      analysis_.mipTimerStart(kMipClockNodeSearchSeparation);
       sepa.separate(search.getLocalDomain());
+      analysis_.mipTimerStop(kMipClockNodeSearchSeparation);
 
       if (mipdata_->domain.infeasible()) {
         search.cutoffNode();
+	analysis_.mipTimerStart(kMipClockOpenNodesToQueue1);
         search.openNodesToQueue(mipdata_->nodequeue);
+	analysis_.mipTimerStop(kMipClockOpenNodesToQueue1);
         mipdata_->nodequeue.clear();
         mipdata_->pruned_treeweight = 1.0;
 
+	analysis_.mipTimerStart(kMipClockStoreBasis);
         double prev_lower_bound = mipdata_->lower_bound;
 
         mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound);
diff --git a/src/mip/MipTimer.h b/src/mip/MipTimer.h
index 6b04972612..e6e281d323 100644
--- a/src/mip/MipTimer.h
+++ b/src/mip/MipTimer.h
@@ -29,14 +29,14 @@ enum iClockMip {
   kMipClockProbingPresolve,
   kMipClockPerformAging1,
   kMipClockDive,
-  kMipClockOpenNodesToQueue,
+  kMipClockOpenNodesToQueue0,
   kMipClockDomainPropgate,
   kMipClockPruneInfeasibleNodes,
   kMipClockUpdateLocalDomain,
   kMipClockNodeSearch,
   //  kMipClock@,
   // Dive
-  kMipClockEvaluateNode,
+  kMipClockEvaluateNode0,
   kMipClockPrimalHeuristics,
   kMipClockTheDive,
   kMipClockBacktrackPlunge,
@@ -58,6 +58,16 @@ enum iClockMip {
   kMipClockRootSeparationRound,
   kMipClockSolveSubMipRootReducedCost,
 
+  // NodeSearch
+  kMipClockCurrentNodeToQueue,
+  kMipClockSearchBacktrack,
+  kMipClockNodePrunedLoop,
+  kMipClockOpenNodesToQueue1,
+  kMipClockEvaluateNode1,
+  kMipClockNodeSearchSeparation,
+  kMipClockStoreBasis,
+  //  kMipClock@,
+
   // Separation
   kMipClockSeparationRootSeparationRound,
   kMipClockSeparationFinishAnalyticCentreComputation,
@@ -108,8 +118,8 @@ class MipTimer {
     // Search - Should correspond to kMipClockSearch
     clock[kMipClockPerformAging1] = timer_pointer->clock_def("Perform aging 1");
     clock[kMipClockDive] = timer_pointer->clock_def("Dive");
-    clock[kMipClockOpenNodesToQueue] =
-        timer_pointer->clock_def("Open nodes to queue");
+    clock[kMipClockOpenNodesToQueue0] =
+        timer_pointer->clock_def("Open nodes to queue 0");
     clock[kMipClockDomainPropgate] =
         timer_pointer->clock_def("Domain propagate");
     clock[kMipClockPruneInfeasibleNodes] =
@@ -120,7 +130,7 @@ class MipTimer {
     //    clock[kMipClock@] = timer_pointer->clock_def("@");
 
     // Dive - Should correspond to kMipClockDive
-    clock[kMipClockEvaluateNode] = timer_pointer->clock_def("Evaluate node");
+    clock[kMipClockEvaluateNode0] = timer_pointer->clock_def("Evaluate node");
     clock[kMipClockPrimalHeuristics] =
         timer_pointer->clock_def("Primal heuristics");
     clock[kMipClockTheDive] = timer_pointer->clock_def("The dive");
@@ -167,6 +177,21 @@ class MipTimer {
     clock[kMipClockSeparationEvaluateRootLp] =
         timer_pointer->clock_def("Evaluate root LP - s.");
 
+    // Node search
+    clock[kMipClockCurrentNodeToQueue] =
+        timer_pointer->clock_def("Current node to queue");
+    clock[kMipClockSearchBacktrack] =
+        timer_pointer->clock_def("Search backtrack");
+    clock[kMipClockNodePrunedLoop] =
+        timer_pointer->clock_def("Pruned loop search");
+    clock[kMipClockOpenNodesToQueue1] =
+        timer_pointer->clock_def("Open nodes to queue 1");
+    clock[kMipClockEvaluateNode1] = timer_pointer->clock_def("Evaluate node 1");
+    clock[kMipClockNodeSearchSeparation] =
+        timer_pointer->clock_def("Node search separation");
+    clock[kMipClockStoreBasis] = timer_pointer->clock_def("Store basis");
+    //    clock[] = timer_pointer->clock_def("");
+
     // Evaluate LPs
     clock[kMipClockSimplexBasisSolveLp] =
         timer_pointer->clock_def("Solve LP - simplex basis");
@@ -278,7 +303,7 @@ class MipTimer {
   void reportMipSearchClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{
         kMipClockPerformAging1,        kMipClockDive,
-        kMipClockOpenNodesToQueue,     kMipClockDomainPropgate,
+        kMipClockOpenNodesToQueue0,    kMipClockDomainPropgate,
         kMipClockPruneInfeasibleNodes, kMipClockUpdateLocalDomain,
         kMipClockNodeSearch,
         //	kMipClock@
@@ -289,7 +314,7 @@ class MipTimer {
 
   void reportMipDiveClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{
-        kMipClockEvaluateNode, kMipClockPrimalHeuristics, kMipClockTheDive,
+        kMipClockEvaluateNode0, kMipClockPrimalHeuristics, kMipClockTheDive,
         kMipClockBacktrackPlunge, kMipClockPerformAging2};
     reportMipClockList("MipDive_", mip_clock_list, mip_timer_clock,
                        kMipClockDive, tolerance_percent_report);
@@ -329,6 +354,16 @@ class MipTimer {
                        kMipClockSeparation);  //, tolerance_percent_report);
   };
 
+  void reportMipNodeSearchClock(const HighsTimerClock& mip_timer_clock) {
+    const std::vector<HighsInt> mip_clock_list{
+        kMipClockCurrentNodeToQueue, kMipClockNodePrunedLoop,
+        //      kMipClockSearchBacktrack,
+        kMipClockOpenNodesToQueue1, kMipClockEvaluateNode1,
+        kMipClockNodeSearchSeparation, kMipClockStoreBasis};
+    reportMipClockList("MipNodeSearch", mip_clock_list, mip_timer_clock,
+                       kMipClockNodeSearch);  //, tolerance_percent_report);
+  };
+
   void csvMipClock(const std::string model_name,
                    const HighsTimerClock& mip_timer_clock, const bool header,
                    const bool end_line) {

From b450b05ea7941d6ec1fb402e307f06247bd7e8ed Mon Sep 17 00:00:00 2001
From: JAJHall <jajhall@ed.ac.uk>
Date: Tue, 25 Feb 2025 10:03:12 +0000
Subject: [PATCH 28/28] Added submip timer and reporting time if number of
 calls is positive when time tolerance is 0

---
 src/lp_data/HighsOptions.h        |  3 ++-
 src/mip/HighsMipAnalysis.cpp      |  1 +
 src/mip/HighsMipSolver.cpp        | 28 ++++++++++++++--------------
 src/mip/HighsPrimalHeuristics.cpp |  3 +++
 src/mip/MipTimer.h                | 16 ++++++++++++++--
 src/util/HighsTimer.h             |  6 +++++-
 6 files changed, 39 insertions(+), 18 deletions(-)

diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h
index f19d4acd8c..f24146bb48 100644
--- a/src/lp_data/HighsOptions.h
+++ b/src/lp_data/HighsOptions.h
@@ -766,7 +766,8 @@ class HighsOptions : public HighsOptionsStruct {
 
     record_int = new OptionRecordInt(
         "highs_analysis_level", "Analysis level in HiGHS", now_advanced,
-        &highs_analysis_level, kHighsAnalysisLevelMin, kHighsAnalysisLevelMin,
+        &highs_analysis_level, kHighsAnalysisLevelMin,
+        kHighsAnalysisLevelMipTime,  // kHighsAnalysisLevelMin,
         kHighsAnalysisLevelMax);
     records.push_back(record_int);
 
diff --git a/src/mip/HighsMipAnalysis.cpp b/src/mip/HighsMipAnalysis.cpp
index e1b9f4ad2d..87c0f57506 100644
--- a/src/mip/HighsMipAnalysis.cpp
+++ b/src/mip/HighsMipAnalysis.cpp
@@ -151,6 +151,7 @@ void HighsMipAnalysis::reportMipTimer() {
   mip_timer.reportMipPrimalHeuristicsClock(mip_clocks);
   mip_timer.reportMipEvaluateRootNodeClock(mip_clocks);
   mip_timer.reportMipSeparationClock(mip_clocks);
+  mip_timer.reportMipSubMipSolveClock(mip_clocks);
   mip_timer.csvMipClock(this->model_name, mip_clocks, true, false);
   reportMipSolveLpClock(true);
   mip_timer.csvMipClock(this->model_name, mip_clocks, false, false);
diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index b72c4f7ec6..1eba575330 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -148,7 +148,6 @@ void HighsMipSolver::run() {
       cleanupSolve();
       return;
     }
-    analysis_.mipTimerStop(kMipClockTrivialHeuristics);
     if (analysis_.analyse_mip_time && !submip)
       highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                    "MIP-Timing: %11.2g - starting evaluate root node\n",
@@ -528,21 +527,21 @@ void HighsMipSolver::run() {
       // new global information before we perform separation rounds for the node
       analysis_.mipTimerStart(kMipClockEvaluateNode1);
       const HighsSearch::NodeResult evaluate_node_result =
-	search.evaluateNode();
+          search.evaluateNode();
       analysis_.mipTimerStop(kMipClockEvaluateNode1);
       if (evaluate_node_result == HighsSearch::NodeResult::kSubOptimal) {
-	analysis_.mipTimerStart(kMipClockCurrentNodeToQueue);
-	search.currentNodeToQueue(mipdata_->nodequeue);
-	analysis_.mipTimerStop(kMipClockCurrentNodeToQueue);
+        analysis_.mipTimerStart(kMipClockCurrentNodeToQueue);
+        search.currentNodeToQueue(mipdata_->nodequeue);
+        analysis_.mipTimerStop(kMipClockCurrentNodeToQueue);
       }
 
       // if the node was pruned we remove it from the search and install the
       // next node from the queue
-        analysis_.mipTimerStart(kMipClockNodePrunedLoop);
+      analysis_.mipTimerStart(kMipClockNodePrunedLoop);
       if (search.currentNodePruned()) {
-	analysis_.mipTimerStart(kMipClockSearchBacktrack);
+        //	analysis_.mipTimerStart(kMipClockSearchBacktrack);
         search.backtrack();
-	analysis_.mipTimerStop(kMipClockSearchBacktrack);
+        //	analysis_.mipTimerStop(kMipClockSearchBacktrack);
         ++mipdata_->num_leaves;
         ++mipdata_->num_nodes;
         search.flushStatistics();
@@ -564,7 +563,7 @@ void HighsMipSolver::run() {
             mipdata_->updatePrimalDualIntegral(
                 prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
                 mipdata_->upper_bound);
-	  analysis_.mipTimerStop(kMipClockNodePrunedLoop);
+          analysis_.mipTimerStop(kMipClockNodePrunedLoop);
           break;
         }
 
@@ -573,6 +572,7 @@ void HighsMipSolver::run() {
           break;
         }
 
+        //	analysis_.mipTimerStart(kMipClockStoreBasis);
         double prev_lower_bound = mipdata_->lower_bound;
 
         mipdata_->lower_bound = std::min(
@@ -600,9 +600,9 @@ void HighsMipSolver::run() {
           mipdata_->domain.clearChangedCols();
           mipdata_->removeFixedIndices();
         }
-	analysis_.mipTimerStop(kMipClockStoreBasis);
+        //	analysis_.mipTimerStop(kMipClockStoreBasis);
 
-	analysis_.mipTimerStop(kMipClockNodePrunedLoop);
+        analysis_.mipTimerStop(kMipClockNodePrunedLoop);
         continue;
       }
       analysis_.mipTimerStop(kMipClockNodePrunedLoop);
@@ -614,13 +614,13 @@ void HighsMipSolver::run() {
 
       if (mipdata_->domain.infeasible()) {
         search.cutoffNode();
-	analysis_.mipTimerStart(kMipClockOpenNodesToQueue1);
+        analysis_.mipTimerStart(kMipClockOpenNodesToQueue1);
         search.openNodesToQueue(mipdata_->nodequeue);
-	analysis_.mipTimerStop(kMipClockOpenNodesToQueue1);
+        analysis_.mipTimerStop(kMipClockOpenNodesToQueue1);
         mipdata_->nodequeue.clear();
         mipdata_->pruned_treeweight = 1.0;
 
-	analysis_.mipTimerStart(kMipClockStoreBasis);
+        analysis_.mipTimerStart(kMipClockStoreBasis);
         double prev_lower_bound = mipdata_->lower_bound;
 
         mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound);
diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp
index fc36b59d62..43f4244a41 100644
--- a/src/mip/HighsPrimalHeuristics.cpp
+++ b/src/mip/HighsPrimalHeuristics.cpp
@@ -134,6 +134,8 @@ bool HighsPrimalHeuristics::solveSubMip(
   solution.value_valid = false;
   solution.dual_valid = false;
   // Create HighsMipSolver instance for sub-MIP
+  if (!mipsolver.submip)
+    mipsolver.analysis_.mipTimerStart(kMipClockSubMipSolve);
   HighsMipSolver submipsolver(*mipsolver.callback_, submipoptions, submip,
                               solution, true, mipsolver.submip_level + 1);
   submipsolver.rootbasis = &basis;
@@ -145,6 +147,7 @@ bool HighsPrimalHeuristics::solveSubMip(
   submipsolver.run();
   mipsolver.max_submip_level =
       std::max(submipsolver.max_submip_level + 1, mipsolver.max_submip_level);
+  if (!mipsolver.submip) mipsolver.analysis_.mipTimerStop(kMipClockSubMipSolve);
   if (submipsolver.mipdata_) {
     double numUnfixed = mipsolver.mipdata_->integral_cols.size() +
                         mipsolver.mipdata_->continuous_cols.size();
diff --git a/src/mip/MipTimer.h b/src/mip/MipTimer.h
index e6e281d323..218c0fbb4c 100644
--- a/src/mip/MipTimer.h
+++ b/src/mip/MipTimer.h
@@ -74,10 +74,14 @@ enum iClockMip {
   kMipClockSeparationCentralRounding,
   kMipClockSeparationEvaluateRootLp,
 
+  // LP solves
   kMipClockSimplexBasisSolveLp,
   kMipClockSimplexNoBasisSolveLp,
   kMipClockIpmSolveLp,
 
+  // Sub-MIP solves
+  kMipClockSubMipSolve,
+
   kMipClockSolveSubMipRENS,
   kMipClockSolveSubMipRINS,
 
@@ -192,6 +196,9 @@ class MipTimer {
     clock[kMipClockStoreBasis] = timer_pointer->clock_def("Store basis");
     //    clock[] = timer_pointer->clock_def("");
 
+    // Sub-MIP clock
+    clock[kMipClockSubMipSolve] = timer_pointer->clock_def("Sub-MIP solves");
+
     // Evaluate LPs
     clock[kMipClockSimplexBasisSolveLp] =
         timer_pointer->clock_def("Solve LP - simplex basis");
@@ -208,7 +215,6 @@ class MipTimer {
     clock[kMipClockProbingImplications] =
         timer_pointer->clock_def("Probing - implications");
     //    clock[] = timer_pointer->clock_def("");
-    //    clock[] = timer_pointer->clock_def("");
   }
 
   bool reportMipClockList(const char* grepStamp,
@@ -294,6 +300,12 @@ class MipTimer {
                        kMipClockTotal);  //, tolerance_percent_report);
   };
 
+  void reportMipSubMipSolveClock(const HighsTimerClock& mip_timer_clock) {
+    const std::vector<HighsInt> mip_clock_list{kMipClockSubMipSolve};
+    reportMipClockList("MipSlvLp", mip_clock_list, mip_timer_clock,
+                       kMipClockTotal);  //, tolerance_percent_report);
+  };
+
   void reportMipPresolveClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{kMipClockProbingPresolve};
     reportMipClockList("MipPrslv", mip_clock_list, mip_timer_clock,
@@ -359,7 +371,7 @@ class MipTimer {
         kMipClockCurrentNodeToQueue, kMipClockNodePrunedLoop,
         //      kMipClockSearchBacktrack,
         kMipClockOpenNodesToQueue1, kMipClockEvaluateNode1,
-        kMipClockNodeSearchSeparation, kMipClockStoreBasis};
+        kMipClockNodeSearchSeparation};  //, kMipClockStoreBasis};
     reportMipClockList("MipNodeSearch", mip_clock_list, mip_timer_clock,
                        kMipClockNodeSearch);  //, tolerance_percent_report);
   };
diff --git a/src/util/HighsTimer.h b/src/util/HighsTimer.h
index eb427524fd..bb54b33a0d 100644
--- a/src/util/HighsTimer.h
+++ b/src/util/HighsTimer.h
@@ -301,7 +301,11 @@ class HighsTimer {
       double time_per_call = 0;
       if (clock_num_call[iClock] > 0) {
         time_per_call = time / clock_num_call[iClock];
-        if (percent_sum_clock_times[i] >= tolerance_percent_report) {
+        const bool report_time =
+            tolerance_percent_report > 0
+                ? percent_sum_clock_times[i] >= tolerance_percent_report
+                : clock_num_call[iClock] > 0;
+        if (report_time) {
           printf("%s-time  %-32s: %11.4e (%5.1f%%", grep_stamp,
                  clock_names[iClock].c_str(), time, percent_run_highs);
           if (ideal_sum_time > 0) {