Primal-dual evolution event handler recipe

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 ### Added
+- Added primal_dual_evolution recipe and a plot recipe
 - Expanded Statistics class to more problems.
 - Created Statistics class
 - Added parser to read .stats file

src/pyscipopt/recipes/primal_dual_evolution.py

@@ -0,0 +1,50 @@
+from pyscipopt import Model, Eventhdlr, SCIP_EVENTTYPE, Eventhdlr
+
+def get_primal_dual_evolution(model: Model):
+
+    class gapEventhdlr(Eventhdlr):
+
+        def eventinit(self): # we want to collect best primal solutions and best dual solutions
+            self.model.catchEvent(SCIP_EVENTTYPE.BESTSOLFOUND, self)
+            self.model.catchEvent(SCIP_EVENTTYPE.LPSOLVED, self)
+            self.model.catchEvent(SCIP_EVENTTYPE.NODESOLVED, self)
+
+
+        def eventexec(self, event):
+            # if a new best primal solution was found, we save when it was found and also its objective
+            if event.getType() == SCIP_EVENTTYPE.BESTSOLFOUND:
+                self.model.data["primal_solutions"].append((self.model.getSolvingTime(), self.model.getPrimalbound()))
+
+            if not self.model.data["dual_solutions"]:
+                self.model.data["dual_solutions"].append((self.model.getSolvingTime(), self.model.getDualbound()))
+
+            if self.model.getObjectiveSense() == "minimize":
+                if self.model.isGT(self.model.getDualbound(), self.model.data["dual_solutions"][-1][1]):
+                    self.model.data["dual_solutions"].append((self.model.getSolvingTime(), self.model.getDualbound()))
+            else:
+                if self.model.isLT(self.model.getDualbound(), self.model.data["dual_solutions"][-1][1]):
+                    self.model.data["dual_solutions"].append((self.model.getSolvingTime(), self.model.getDualbound()))
+
+    if not hasattr(model, "data"):
+        model.data = {}
+
+    model.data["primal_solutions"] = []
+    model.data["dual_solutions"] = []
+    hdlr = gapEventhdlr()
+    model.includeEventhdlr(hdlr, "gapEventHandler", "Event handler which collects primal and dual solution evolution")
+
+    return model
+
+def plot_primal_dual_evolution(model: Model):
+    try:
+        from matplotlib import pyplot as plt
+    except ImportError:
+        raise("matplotlib is required to plot the solution. Try running `pip install matplotlib` in the command line.")
+
+    time_primal, val_primal = zip(*model.data["primal_solutions"])
+    plt.plot(time_primal, val_primal, label="Primal bound")
+    time_dual, val_dual = zip(*model.data["dual_solutions"])
+    plt.plot(time_dual, val_dual, label="Dual bound")
+
+    plt.legend(loc="best")
+    plt.show()

tests/helpers/utils.py

@@ -77,19 +77,15 @@ def random_nlp_1():
     return model
 
 
-def knapsack_model(weights=[4, 2, 6, 3, 7, 5], costs=[7, 2, 5, 4, 3, 4]):
+def knapsack_model(weights=[4, 2, 6, 3, 7, 5], costs=[7, 2, 5, 4, 3, 4], knapsackSize = 15):
     # create solver instance
     s = Model("Knapsack")
-    s.hideOutput()
 
     # setting the objective sense to maximise
     s.setMaximize()
 
     assert len(weights) == len(costs)
 
-    # knapsack size
-    knapsackSize = 15
-
     # adding the knapsack variables
     knapsackVars = []
     varNames = []

tests/test_model.py

@@ -476,4 +476,4 @@ def test_getObjVal():
     assert m.getVal(x) == 0
 
     assert m.getObjVal() == 16 
-    assert m.getVal(x) == 0
+    assert m.getVal(x) == 0

tests/test_recipe_primal_dual_evolution.py

@@ -0,0 +1,36 @@
+from pyscipopt.recipes.primal_dual_evolution import get_primal_dual_evolution
+from helpers.utils import bin_packing_model
+
+def test_primal_dual_evolution():
+    from random import randint
+
+    model = bin_packing_model(sizes=[randint(1,40) for _ in range(120)],  capacity=50)
+    model.setParam("limits/time",5)
+
+    model.data = {"test": True}
+    model = get_primal_dual_evolution(model)
+
+    assert "test" in model.data
+    assert "primal_solutions" in model.data
+
+    model.optimize()
+
+    # these are required because the event handler doesn't capture the final state
+    model.data["primal_solutions"].append((model.getSolvingTime(), model.getPrimalbound()))
+    model.data["dual_solutions"].append((model.getSolvingTime(), model.getDualbound()))
+
+    for i in range(1, len(model.data["primal_solutions"])):
+        if model.getObjectiveSense() == "minimize":
+            assert model.data["primal_solutions"][i][1] <= model.data["primal_solutions"][i-1][1]
+        else:
+            assert model.data["primal_solutions"][i][1] >= model.data["primal_solutions"][i-1][1]
+
+    for i in range(1, len(model.data["dual_solutions"])):
+        if model.getObjectiveSense() == "minimize":
+            assert model.data["dual_solutions"][i][1] >= model.data["dual_solutions"][i-1][1]
+        else:
+            assert model.data["dual_solutions"][i][1] <= model.data["dual_solutions"][i-1][1]
+
+    # how to get a simple plot of the data
+    #from pyscipopt.recipes.primal_dual_evolution import plot_primal_dual_evolution
+    #plot_primal_dual_evolution(model)