First draft before holiday

Author: Opt-Mucca

highs/mip/HighsPathSeparator.cpp

@@ -58,15 +58,23 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation,
   }
 
   std::vector<HighsInt> numContinuous(lp.num_row_);
-
+  // Score will be used for deciding order in which rows get aggregated
+  // TODO: The second entry is currently a hacky way to store the norm
+  std::vector<std::pair<double, double>> rowscore(lp.num_row_,
+                                                  std::make_pair(0.0, 0.0));
   size_t maxAggrRowSize = 0;
   for (HighsInt col : mip.mipdata_->continuous_cols) {
     if (transLp.boundDistance(col) == 0.0) continue;
 
     maxAggrRowSize += lp.a_matrix_.start_[col + 1] - lp.a_matrix_.start_[col];
     for (HighsInt i = lp.a_matrix_.start_[col];
-         i != lp.a_matrix_.start_[col + 1]; ++i)
+         i != lp.a_matrix_.start_[col + 1]; ++i) {
       ++numContinuous[lp.a_matrix_.index_[i]];
+      // Add the fractional score of the row
+      rowscore[i].first +=
+          lp.a_matrix_.value_[i] * transLp.getFracVbEstimate(col);
+      rowscore[i].second += lp.a_matrix_.value_[i];
+    }
   }
 
   std::vector<std::pair<HighsInt, double>> colSubstitutions(

highs/mip/HighsTransformedLp.h

@@ -16,6 +16,7 @@
 
 #include <vector>
 
+#include "HighsLpRelaxation.h"
 #include "lp_data/HConst.h"
 #include "mip/HighsImplications.h"
 #include "util/HighsCDouble.h"
@@ -52,6 +53,30 @@ class HighsTransformedLp {
 
   double boundDistance(HighsInt col) const { return boundDist[col]; }
 
+  // TODO: Should this simple be calculated when we initialise the transLp?
+  double getFracVbEstimate(HighsInt col) const {
+    HighsInt vbCol = -1;
+    double vbCoef = 0.0;
+    if (ubDist[col] <= lbDist[col] && bestVlb[col].first != -1) {
+      vbCol = bestVlb[col].first;
+      vbCoef = bestVlb[col].second.coef;
+    } else if (lbDist[col] <= ubDist[col] && bestVub[col].first != -1) {
+      vbCol = bestVub[col].first;
+      vbCoef = bestVub[col].second.coef;
+    } else if (bestVlb[col].first != -1) {
+      vbCol = bestVlb[col].first;
+      vbCoef = bestVlb[col].second.coef;
+    } else if (bestVub[col].first != -1) {
+      vbCol = bestVub[col].first;
+      vbCoef = bestVub[col].second.coef;
+    }
+    if (vbCol == -1) return 0;
+    // TODO: This differs from the SCIP definition? Their frac can be > 1?
+    double val = vbCoef * lprelaxation.solutionValue(col);
+    double frac = std::max(val - std::floor(val), 0.0);
+    return std::min(frac, 1 - frac);
+  }
+
   bool transform(std::vector<double>& vals, std::vector<double>& upper,
                  std::vector<double>& solval, std::vector<HighsInt>& inds,
                  double& rhs, bool& integralPositive, bool preferVbds = false);