Merge branch 'latest' of https://github.com/ERGO-Code/HiGHS into pres…

…olveBugTransform

CONTRIBUTING.md

@@ -12,7 +12,7 @@ Contact [Julian](https://github.com/jajhall) (General issues and solvers), [Ivet
 
 ## Improve the documentation
 
-The top level [documentation](https://ergo-code.github.io/HiGHS/) is created using [Docsy](https://www.docsy.dev/), with the files held on the [HiGHS repository](https://github.com/ERGO-Code/HiGHS/tree/docsy). If your change is small (like fixing typos, or one or two sentence corrections), the easiest way to do this is to fork the branch and create a pull request. (See *Contribute code to HiGHS* below for more on this.) If your change is larger, or touches multiple files, please raise an issue describing what you want to do.
+The top level [documentation](https://ergo-code.github.io/HiGHS/) is created using [Docsy](https://www.docsy.dev/), with the files held on the [HiGHS repository](https://github.com/ERGO-Code/HiGHS/tree/master/docs). If your change is small (like fixing typos, or one or two sentence corrections), the easiest way to do this is to fork the branch and create a pull request. (See *Contribute code to HiGHS* below for more on this.) If your change is larger, or touches multiple files, please raise an issue describing what you want to do.
 
 ## Raise an issue
 

app/RunHighs.cpp

@@ -62,6 +62,24 @@ int main(int argc, char** argv) {
       return (int)read_solution_status;
     }
   }
+  if (options.write_presolved_model_to_file) {
+    // Run presolve and write the presolved model to a file
+    HighsStatus status = highs.presolve();
+    if (status == HighsStatus::kError) return int(status);
+    HighsPresolveStatus model_presolve_status = highs.getModelPresolveStatus();
+    const bool ok_to_write =
+        model_presolve_status == HighsPresolveStatus::kNotReduced ||
+        model_presolve_status == HighsPresolveStatus::kReduced ||
+        model_presolve_status == HighsPresolveStatus::kReducedToEmpty ||
+        model_presolve_status == HighsPresolveStatus::kTimeout;
+    if (!ok_to_write) {
+      highsLogUser(log_options, HighsLogType::kInfo,
+                   "No presolved model to write to file\n");
+      return int(status);
+    }
+    status = highs.writePresolvedModel(options.write_presolved_model_file);
+    return int(status);
+  }
   // Solve the model
   HighsStatus run_status = highs.run();
   if (run_status == HighsStatus::kError) return int(run_status);

docs/src/index.md

@@ -59,6 +59,19 @@ The C interface cannot make use of the C++ structures and enums, and its methods
 For LPs, HiGHS has implementations of both the revised simplex and interior
 point methods. MIPs are solved by branch-and-cut, and QPs by active set.
 
+For LP, the novel features of the dual simplex solver are described in
+_Parallelizing the dual revised simplex method_, Q. Huangfu and
+J. A. J. Hall, Mathematical Programming Computation, 10 (1), 119-142,
+2018 [DOI:
+10.1007/s12532-017-0130-5](https://link.springer.com/article/10.1007/s12532-017-0130-5). For
+the interior point solver, the reference is _Implementation of an
+interior point method with basis preconditioning_, Mathematical
+Programming Computation, 12, 603-635, 2020. [DOI:
+10.1007/s12532-020-00181-8](https://link.springer.com/article/10.1007/s12532-020-00181-8). There
+are no specific references to the techniques used in the MIP or QP
+solvers.
+
+
 ## Citing HiGHS
 
 If you use HiGHS in an academic context, please cite the following article:

docs/src/terminology.md

@@ -139,3 +139,4 @@ relative to the primal bound is a better measure. When the gap reaches
 zero then the MIP is solved to optimality. However, it is often
 preferable to stop the MIP solver when the relative gap is below a
 specified tolerance.
+

examples/chip.py

@@ -9,23 +9,12 @@
 h = highspy.Highs()
 h.silent()
 
-varNames = list()
-varNames.append('Tables')
-varNames.append('Sets of chairs')
+items = ['Tables', 'Sets of chairs']
+x = h.addVariables(items, obj = [10, 20], name = items)
 
-x1 = h.addVariable(obj = 10, name = varNames[0])
-x2 = h.addVariable(obj = 25, name = varNames[1])
-
-vars = list()
-vars.append(x1)
-vars.append(x2)
-
-constrNames = list()
-constrNames.append('Assembly')
-constrNames.append('Finishing')
-
-h.addConstr(x1 + 2*x2 <=  80, name = constrNames[0])
-h.addConstr(x1 + 4*x2 <= 120, name = constrNames[1])
+constrNames = ['Assembly', 'Finishing']
+cons = h.addConstrs(x['Tables'] + 2*x['Sets of chairs'] <=  80, 
+                    x['Tables'] + 4*x['Sets of chairs'] <= 120, name = constrNames)
 
 h.setMaximize()
 
@@ -36,16 +25,17 @@
 
 h.solve()
 
-for var in vars:
+
+for n, var in x.items():
     print('Make', h.variableValue(var), h.variableName(var), ': Reduced cost', h.variableDual(var))
-print('Make', h.variableValues(vars), 'of', h.variableNames(vars))
+
+print('Make', h.variableValues(x.values()), 'of', h.variableNames(x.values()))
 print('Make', h.allVariableValues(), 'of', h.allVariableNames())
 
-for name in constrNames:
-    print('Constraint', name, 'has value', h.constrValue(name), 'and dual', h.constrDual(name))
-
-print('Constraints have values', h.constrValues(constrNames), 'and duals', h.constrDuals(constrNames))
-
+for c in cons:
+    print('Constraint', c.name, 'has value', h.constrValue(c), 'and dual', h.constrDual(c))
+
+print('Constraints have values', h.constrValues(cons), 'and duals', h.constrDuals(cons))
 print('Constraints have values', h.allConstrValues(), 'and duals', h.allConstrDuals())
 
 print('Optimal objective value is', h.getObjectiveValue())

examples/distillation.py

@@ -15,25 +15,13 @@
 h = highspy.Highs()
 h.silent()
 
-variableNames = list()
-variableNames.append('TypeA')
-variableNames.append('TypeB')
+variableNames = ['TypeA', 'TypeB']
+x = h.addVariables(variableNames, obj = [8, 10], name = variableNames[0])
 
-useTypeA = h.addVariable(obj =  8, name = variableNames[0])
-useTypeB = h.addVariable(obj = 10, name = variableNames[1])
-
-vars = list()
-vars.append(useTypeA)
-vars.append(useTypeB)
-
-constrNames = list()
-constrNames.append('Product1')
-constrNames.append('Product2')
-constrNames.append('Product3')
-
-h.addConstr(2*useTypeA + 2*useTypeB >=  7, name = constrNames[0])
-h.addConstr(3*useTypeA + 4*useTypeB >= 12, name = constrNames[1])
-h.addConstr(2*useTypeA +   useTypeB >=  6, name = constrNames[2])
+constrNames = ['Product1', 'Product2', 'Product3']
+cons = h.addConstrs(2*x['TypeA'] + 2*x['TypeB'] >=  7,
+                    3*x['TypeA'] + 4*x['TypeB'] >= 12,
+                    2*x['TypeA'] +   x['TypeB'] >=  6, name = constrNames)
 
 h.setMinimize()
 
@@ -47,46 +35,46 @@
 
 h.solve()
 
-for var in vars:
+for name, var in x.items():
     print('Use {0:.1f} of {1:s}: reduced cost {2:.6f}'.format(h.variableValue(var), h.variableName(var), h.variableDual(var)))
-print('Use', h.variableValues(vars), 'of', h.variableNames(vars))
+
+print('Use', h.variableValues(x.values()), 'of', h.variableNames(x.values()))
 print('Use', h.allVariableValues(), 'of', h.allVariableNames())
 
-for name in constrNames:
-    print(f"Constraint {name} has value {h.constrValue(name):{width}.{precision}} and dual  {h.constrDual(name):{width}.{precision}}")
-
-print('Constraints have values', h.constrValues(constrNames), 'and duals', h.constrDuals(constrNames))
+for c in cons:
+    print(f"Constraint {c.name} has value {h.constrValue(c):{width}.{precision}} and dual  {h.constrDual(c):{width}.{precision}}")
 
+print('Constraints have values', h.constrValues(cons), 'and duals', h.constrDuals(cons))
 print('Constraints have values', h.allConstrValues(), 'and duals', h.allConstrDuals())
 
-for var in vars:
+for var in x.values():
     print(f"Use {h.variableValue(var):{width}.{precision}} of {h.variableName(var)}")
 print(f"Optimal objective value is {h.getObjectiveValue():{width}.{precision}}")
 
 print()
 print('Solve as MIP')
 
-for var in vars:
+for var in x.values():
     h.setInteger(var)
 
 h.solve()
 
-for var in vars:
+for var in x.values():
     print(f"Use {h.variableValue(var):{width}.{precision}} of {h.variableName(var)}")
 print(f"Optimal objective value is {h.getObjectiveValue():{width}.{precision}}")
 
 print()
 print('Solve as LP with Gomory cut')
 
 # Make the variables continuous
-for var in vars:
+for var in x.values():
     h.setContinuous(var)
 
 # Add Gomory cut
-h.addConstr(useTypeA + useTypeB >= 4, name = "Gomory")
+h.addConstr(x['TypeA'] + x['TypeB'] >= 4, name = "Gomory")
 
 h.solve()
 
-for var in vars:
+for var in x.values():
     print(f"Use {h.variableValue(var):{width}.{precision}} of {h.variableName(var)}")
 print(f"Optimal objective value is {h.getObjectiveValue():{width}.{precision}}")

examples/minimal.py

@@ -5,7 +5,7 @@
 x1 = h.addVariable(lb = -h.inf)
 x2 = h.addVariable(lb = -h.inf)
 
-h.addConstr(x2 - x1 >= 2)
-h.addConstr(x1 + x2 >= 0)
+h.addConstrs(x2 - x1 >= 2, 
+             x1 + x2 >= 0)
 
 h.minimize(x2)

examples/network_flow.py

@@ -0,0 +1,37 @@
+# Example of a shortest path network flow in a graph
+# Shows integration of highspy with networkx
+
+import highspy
+import networkx as nx
+
+orig, dest = ('A', 'D')
+
+# create directed graph with edge weights (distances)
+G = nx.DiGraph()
+G.add_weighted_edges_from([('A', 'B', 2.0), ('B', 'C', 3.0), ('A', 'C', 1.5), ('B', 'D', 2.5), ('C', 'D', 1.0)])
+
+h = highspy.Highs()
+h.silent()
+
+x = h.addBinaries(G.edges, obj=nx.get_edge_attributes(G, 'weight'))
+
+# add flow conservation constraints
+#                        {  1  if n = orig
+#   sum(out) - sum(in) = { -1  if n = dest
+#                        {  0     otherwise
+rhs  = lambda n: 1 if n == orig else -1 if n == dest else 0
+flow = lambda E: sum((x[e] for e in E))
+
+h.addConstrs(flow(G.out_edges(n)) - flow(G.in_edges(n)) == rhs(n) for n in G.nodes)
+h.minimize()
+
+# Print the solution
+print('Shortest path from', orig, 'to', dest, 'is: ', end = '')
+sol = h.vals(x)
+
+n = orig
+while n != dest:
+    print(n, end=' ')
+    n = next(e[1] for e in G.out_edges(n) if sol[e] > 0.5)
+
+print(dest)

examples/nqueens.py

@@ -0,0 +1,29 @@
+# This is an example of the N-Queens problem, which is a classic combinatorial problem.
+# The problem is to place N queens on an N x N chessboard so that no two queens attack each other.
+#
+# We show how to model the problem as a MIP and solve it using highspy.
+# Using numpy can simplify the construction of the constraints (i.e., diagonal).
+
+import highspy
+import numpy as np
+
+N = 8
+h = highspy.Highs()
+h.silent()
+
+x = np.reshape(h.addBinaries(N*N), (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
+
+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.solve()
+sol = np.array(h.vals(x))
+
+print('Queens:')
+
+for i in range(N):
+    print(''.join('Q' if sol[i, j] > 0.5 else '*' for j in range(N)))

src/Highs.h

@@ -1374,7 +1374,7 @@ class Highs {
   HighsStatus openWriteFile(const string filename, const string method_name,
                             FILE*& file, HighsFileType& file_type) const;
 
-  void reportModel();
+  void reportModel(const HighsModel& model);
   void newHighsBasis();
   void forceHighsSolutionBasisSize();
   //