highspy: updated highs_linear_expression to be immutable for issue #1888

FEATURES.md

@@ -14,5 +14,8 @@ Added options `write_presolved_model_to_file` and `write_presolved_model_file` s
 
 Added `Highs::feasibilityRelaxation` to solve the problem of minimizing a (possibly weighted) sum of (allowable) infeasibilities in an LP/MIP.
 
+Added Python utility `examples/plot_highs_log.py` (due to @Thell) to visualise progress of the MIP solver.
+
+
 
 

examples/call_highs_from_csharp.cs

@@ -6,17 +6,22 @@
 
 class Program {
    static void Main(string[] args) {
-      double[] cc = {1, -2};
-      double[] cl = {0, 0};
-      double[] cu = {10, 10};
-      double[] rl = {0, 0};
-      double[] ru = {2, 1};
-      int[] astart = {0, 2};
-      int[] aindex = {0, 1, 0, 1};
-      double[] avalue = {1, 2, 1, 3};
+       // Illustrate the solution of a QP, after first solving just the LP
+       //
+       // minimize x_2 + (1/2)(2x_1^2 - 2x_1x_3 + 0.2x_2^2 + 2x_3^2)
+       //
+       // subject to x_1 + x_2 + x_3 >= 1; x>=0
+       double[] cc = {0, 1, 0};
+      double[] cl = {0, 0, 0};
+      double[] cu = {1.0e30, 1.0e30, 1.0e30};
+      double[] rl = {1};
+      double[] ru = {1.0e30};
+      int[] astart = {0, 3};
+      int[] aindex = {0, 1, 2};
+      double[] avalue = {1, 1, 1};
       HighsObjectiveSense sense = HighsObjectiveSense.kMinimize;
       double offset = 0;
-      HighsMatrixFormat a_format = HighsMatrixFormat.kColwise;
+      HighsMatrixFormat a_format = HighsMatrixFormat.kRowwise;
 
       HighsModel model = new HighsModel(cc, cl, cu, rl, ru, astart, aindex, avalue, null, offset, a_format, sense);
 
@@ -31,6 +36,9 @@ static void Main(string[] args) {
       Console.WriteLine("Status: " + status);
       Console.WriteLine("Modelstatus: " + modelStatus);
 
+      for (int i=0; i<sol.colvalue.Length; i++) {
+         Console.WriteLine("Activity for col " + i + " = " + sol.colvalue[i]);
+      }
       for (int i=0; i<sol.rowvalue.Length; i++) {
          Console.WriteLine("Activity for row " + i + " = " + sol.rowvalue[i]);
       }
@@ -43,5 +51,33 @@ static void Main(string[] args) {
       for (int i=0; i<sol.colvalue.Length; i++) {
          Console.WriteLine("x" + i + " = " + sol.colvalue[i] + " is " + bas.colbasisstatus[i]);
       }
+       // Add the Hessian
+      int dim = 2;
+      int[] qstart = {0, 2, 3};
+      int[] qindex = {0, 1, 1};
+      double[] qvalue = {2, -1, 2};
+      HessianFormat q_format = HessianFormat.kTriangular;
+      HighsHessian hessian = new HighsHessian(dim, qstart, qindex, qvalue, q_format);
+      status = solver.passHessian(hessian);
+      status = solver.run();
+      sol = solver.getSolution();
+      modelStatus = solver.GetModelStatus();
+
+      Console.WriteLine("Status: " + status);
+      Console.WriteLine("Modelstatus: " + modelStatus);
+
+      for (int i=0; i<sol.colvalue.Length; i++) {
+         Console.WriteLine("Activity for col " + i + " = " + sol.colvalue[i]);
+      }
+      for (int i=0; i<sol.rowvalue.Length; i++) {
+         Console.WriteLine("Activity for row " + i + " = " + sol.rowvalue[i]);
+      }
+      for (int i=0; i<sol.coldual.Length; i++) {
+         Console.WriteLine("Reduced cost x[" + i + "] = " + sol.coldual[i]);
+      }
+      for (int i=0; i<sol.rowdual.Length; i++) {
+         Console.WriteLine("Dual value for row " + i + " = " + sol.rowdual[i]);
+      }
+
    }
-}
+}

examples/nqueens.py

@@ -11,17 +11,17 @@
 h = highspy.Highs()
 h.silent()
 
-x = np.reshape(h.addBinaries(N*N), (N, N))
+x = h.addBinaries(N, N)
 
-h.addConstrs(sum(x[i,:]) == 1 for i in range(N))    # each row has exactly one queen
-h.addConstrs(sum(x[:,j]) == 1 for j in range(N))    # each col has exactly one queen
+h.addConstrs(h.qsum(x[i,:]) == 1 for i in range(N))    # each row has exactly one queen
+h.addConstrs(h.qsum(x[:,j]) == 1 for j in range(N))    # each col has exactly one queen
 
 y = np.fliplr(x)
-h.addConstrs(x.diagonal(k).sum() <= 1 for k in range(-N + 1, N))   # each diagonal has at most one queen
-h.addConstrs(y.diagonal(k).sum() <= 1 for k in range(-N + 1, N))   # each 'reverse' diagonal has at most one queen
+h.addConstrs(h.qsum(x.diagonal(k)) <= 1 for k in range(-N + 1, N))   # each diagonal has at most one queen
+h.addConstrs(h.qsum(y.diagonal(k)) <= 1 for k in range(-N + 1, N))   # each 'reverse' diagonal has at most one queen
 
 h.solve()
-sol = np.array(h.vals(x))
+sol = h.vals(x)
 
 print('Queens:')
 

examples/plot_highs_log.py

@@ -0,0 +1,134 @@
+import re
+import matplotlib.pyplot as plt
+import numpy as np
+import math as math
+
+
+def parse_highs_log(log_file_path):
+    last_full_entry = []
+    current_entry = []
+    found_solution = False
+
+    with open(log_file_path, "r") as f:
+        for line in f:
+            if "Running HiGHS" in line:
+                if found_solution:
+                    last_full_entry = current_entry
+                current_entry = [line]
+                found_solution = False
+            else:
+                current_entry.append(line)
+                if "Writing the solution to" in line:
+                    found_solution = True
+
+    if not last_full_entry:
+        last_full_entry = current_entry
+
+    if not last_full_entry:
+        return None, None, None, None, None, None
+
+    time_values, best_bound_values, best_sol_values, in_queue_values, expl_values, gap_values = (
+        [],
+        [],
+        [],
+        [],
+        [],
+        [],
+    )
+    for line in last_full_entry:
+        match = re.search(r"\dk?\s+\d+\.\ds$", line)
+
+        if not match:
+            continue
+
+        tokens = line.split()
+        if len(tokens) == 13:
+            tokens = tokens[1:]
+        assert len(tokens) == 12, f"{line}"
+
+        in_queue_values.append(float(tokens[1]))  # InQueue
+        expl_values.append(float(tokens[3].replace("%", "")))  # Expl.%
+        best_bound_values.append(float(tokens[4].replace("inf", "nan")))  # Best Bound
+        best_sol_values.append(float(tokens[5].replace("inf", "nan")))  # Best Sol
+        gap_values.append(
+            float(tokens[6].replace("%", "").replace("inf", "nan").replace("Large", "nan"))
+        )  # Gap%
+        time_values.append(float(tokens[11].replace("s", "")))  # Time
+
+    return time_values, best_bound_values, best_sol_values, in_queue_values, expl_values, gap_values
+
+
+def plot_highs_log(
+    time_values, best_bound_values, best_sol_values, in_queue_values, expl_values, gap_values
+):
+    fig, ax1 = plt.subplots(figsize=(10, 6))
+
+    best_bound_colour = "blue"
+    best_solution_colour = "green"
+    in_queue_colour = "red"
+    explored_colour = "purple"
+    gap_colour = "orange"
+
+    # Plot Objective Bounds
+    ax1.plot(time_values, best_bound_values, label="Best Bound", color=best_bound_colour)
+    ax1.plot(time_values, best_sol_values, label="Best Solution", color=best_solution_colour)
+    ax1.set_xlabel("Time (seconds)")
+    ax1.set_ylabel("Objective Bounds", color=best_bound_colour, labelpad=15)
+    ax1.tick_params(axis="y", labelcolor=best_bound_colour)
+
+    # Limit y-axis to the range between min and max of the non-NaN values
+    valid_gap_index = next(i for i, gap in enumerate(gap_values) if not np.isnan(gap))
+    min_y = min(best_bound_values[valid_gap_index], best_sol_values[valid_gap_index])
+    max_y = max(best_bound_values[valid_gap_index], best_sol_values[valid_gap_index])
+    padding = (max_y - min_y) * 0.1
+    ax1.set_ylim(min_y - padding, max_y + padding)
+
+    # Add second y-axis for InQueue values
+    ax2 = ax1.twinx()
+    ax2.plot(time_values, in_queue_values, label="InQueue", color=in_queue_colour)
+#    ax2.set_ylabel("InQueue", color=in_queue_colour, loc="top", labelpad=12)
+    ax2.yaxis.label.set_rotation(0)
+    ax2.tick_params(axis="y", labelcolor=in_queue_colour)
+
+    # Add third y-axis for Explored % values (scaled)
+    ax3 = ax1.twinx()
+    ax3.spines["right"].set_position(("outward", 50))
+    ax3.plot(time_values, expl_values, label="Explored %", color=explored_colour)
+#    ax3.set_ylabel("Expl.%", color=explored_colour, loc="top", labelpad=10)
+    ax3.yaxis.label.set_rotation(0)
+    ax3.tick_params(axis="y", labelcolor=explored_colour)
+
+    # Add fourth y-axis for Gap % values (scaled)
+    ax4 = ax1.twinx()
+    ax4.spines["right"].set_position(("outward", 90))
+    ax4.plot(time_values, gap_values, label="Gap %", color=gap_colour, linestyle="--", linewidth=0.5)
+#    ax4.set_ylabel("Gap.%", color=gap_colour, loc="top", labelpad=22)
+    ax4.yaxis.label.set_rotation(0)
+    ax4.tick_params(axis="y", labelcolor=gap_colour)
+
+    # Plot vertical hash lines where Best Solution changes
+    for i in range(1, len(best_sol_values)):
+        # Print if change detected and not NaN
+        if (best_sol_values[i] != best_sol_values[i - 1]) and not(math.isnan(best_sol_values[i])):
+            ax1.axvline(x=time_values[i], color="grey", linestyle="--", linewidth=0.5)
+
+    # Shift plot area left to make room on the right for the three y-axis labels.
+    fig.subplots_adjust(left=0.08, right=0.85)
+
+    # Set up legend
+    fig.legend(loc="lower center", ncols=5)
+
+    # Show plot
+    plt.title("HiGHS MIP Log Analysis", pad=20)
+    plt.show()
+
+
+#log_file_path = "/path/to/your/logfile.log"
+log_file_path = "HiGHS.log"
+time_values, best_bound_values, best_sol_values, in_queue_values, expl_values, gap_values = (
+    parse_highs_log(log_file_path)
+)
+
+plot_highs_log(
+    time_values, best_bound_values, best_sol_values, in_queue_values, expl_values, gap_values
+)