Merge pull request #916 from scipopt/plot-pd-evolution

Primal-dual evolution event handler recipe

Author: Opt-Mucca

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 ### Added
+- Added primal_dual_evolution recipe and a plot recipe
 ### Fixed
 ### Changed
 ### Removed

examples/finished/__init__.py

@@ -0,0 +1,54 @@
+"""
+This example show how to retrieve the primal and dual solutions during the optimization process
+and plot them as a function of time. The model is about gas transportation and can be found in
+PySCIPOpt/tests/helpers/utils.py
+
+It makes use of the attach_primal_dual_evolution_eventhdlr recipe.
+
+Requires matplotlib, and may require PyQt6 to show the plot.
+"""
+
+from pyscipopt import Model
+
+def plot_primal_dual_evolution(model: Model):
+    try:
+        from matplotlib import pyplot as plt
+    except ImportError:
+        raise ImportError("matplotlib is required to plot the solution. Try running `pip install matplotlib` in the command line.\
+                          You may also need to install PyQt6 to show the plot.")
+
+    assert model.data["primal_log"], "Could not find any feasible solutions"
+    time_primal, val_primal = map(list,zip(*model.data["primal_log"]))
+    time_dual, val_dual = map(list,zip(*model.data["dual_log"]))
+
+    
+    if time_primal[-1] < time_dual[-1]:
+        time_primal.append(time_dual[-1])
+        val_primal.append(val_primal[-1])
+
+    if time_primal[-1] > time_dual[-1]:
+        time_dual.append(time_primal[-1])
+        val_dual.append(val_dual[-1])
+        
+    plt.plot(time_primal, val_primal, label="Primal bound")
+    plt.plot(time_dual, val_dual, label="Dual bound")
+
+    plt.legend(loc="best")
+    plt.show()
+
+if __name__=="__main__":
+    from pyscipopt.recipes.primal_dual_evolution import attach_primal_dual_evolution_eventhdlr
+    import os
+    import sys
+
+    # just a way to import files from different folders, not important
+    sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '../../tests/helpers')))
+    
+    from utils import gastrans_model
+    
+    model = gastrans_model()
+    model.data = {}
+    attach_primal_dual_evolution_eventhdlr(model)
+
+    model.optimize()
+    plot_primal_dual_evolution(model)

examples/finished/plot_primal_dual_evolution.py

@@ -0,0 +1,46 @@
+from pyscipopt import Model, Eventhdlr, SCIP_EVENTTYPE, Eventhdlr
+
+def attach_primal_dual_evolution_eventhdlr(model: Model):
+    """
+    Attaches an event handler to a given SCIP model that collects primal and dual solutions,
+    along with the solving time when they were found.
+    The data is saved in model.data["primal_log"] and model.data["dual_log"]. They consist of
+    a list of tuples, each tuple containing the solving time and the corresponding solution.
+
+    A usage example can be found in examples/finished/plot_primal_dual_evolution.py. The
+    example takes the information provided by this recipe and uses it to plot the evolution
+    of the dual and primal bounds over time. 
+    """
+    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_log"].append([self.model.getSolvingTime(), self.model.getPrimalbound()])
+            
+            if not self.model.data["dual_log"]:
+                self.model.data["dual_log"].append([self.model.getSolvingTime(), self.model.getDualbound()])
+            
+            if self.model.getObjectiveSense() == "minimize":
+                if self.model.isGT(self.model.getDualbound(), self.model.data["dual_log"][-1][1]):
+                    self.model.data["dual_log"].append([self.model.getSolvingTime(), self.model.getDualbound()])
+            else:
+                if self.model.isLT(self.model.getDualbound(), self.model.data["dual_log"][-1][1]):
+                    self.model.data["dual_log"].append([self.model.getSolvingTime(), self.model.getDualbound()])
+                    
+
+    if not hasattr(model, "data") or model.data==None:
+        model.data = {}
+
+    model.data["primal_log"] = []
+    model.data["dual_log"] = []
+    hdlr = GapEventhdlr()
+    model.includeEventhdlr(hdlr, "gapEventHandler", "Event handler which collects primal and dual solution evolution")
+
+    return model

src/pyscipopt/recipes/primal_dual_evolution.py

@@ -1,14 +1,6 @@
 from pyscipopt import Model, quicksum, SCIP_PARAMSETTING, exp, log, sqrt, sin
 from typing import List
 
-from pyscipopt.scip import is_memory_freed
-
-
-def is_optimized_mode():
-    model = Model()
-    return is_memory_freed()
-
-
 def random_mip_1(disable_sepa=True, disable_huer=True, disable_presolve=True, node_lim=2000, small=False):
     model = Model()
 
@@ -77,19 +69,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], knapsack_size = 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/data/readStatistics.stats

@@ -4,9 +4,9 @@
 import pytest
 
 from pyscipopt import Model, Heur, SCIP_RESULT, SCIP_PARAMSETTING, SCIP_HEURTIMING, SCIP_LPSOLSTAT
-from pyscipopt.scip import is_memory_freed
+from test_memory import is_optimized_mode
 
-from helpers.utils import random_mip_1, is_optimized_mode
+from helpers.utils import random_mip_1
 
 class MyHeur(Heur):
 

tests/helpers/utils.py

@@ -1,6 +1,5 @@
 import pytest
 from pyscipopt.scip import Model, is_memory_freed, print_memory_in_use
-from helpers.utils import is_optimized_mode
 
 def test_not_freed():
     if is_optimized_mode():
@@ -16,4 +15,8 @@ def test_freed():
     assert is_memory_freed()
 
 def test_print_memory_in_use():
-    print_memory_in_use()
+    print_memory_in_use()
+
+def is_optimized_mode():
+    model = Model()
+    return is_memory_freed()

tests/test_heur.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_memory.py

@@ -0,0 +1,28 @@
+from pyscipopt.recipes.primal_dual_evolution import attach_primal_dual_evolution_eventhdlr
+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 = attach_primal_dual_evolution_eventhdlr(model)
+
+    assert "test" in model.data
+    assert "primal_log" in model.data
+
+    model.optimize()
+
+    for i in range(1, len(model.data["primal_log"])):
+        if model.getObjectiveSense() == "minimize":
+            assert model.data["primal_log"][i][1] <= model.data["primal_log"][i-1][1]
+        else:
+            assert model.data["primal_log"][i][1] >= model.data["primal_log"][i-1][1]
+    
+    for i in range(1, len(model.data["dual_log"])):
+        if model.getObjectiveSense() == "minimize":
+            assert model.data["dual_log"][i][1] >= model.data["dual_log"][i-1][1]
+        else:
+            assert model.data["dual_log"][i][1] <= model.data["dual_log"][i-1][1]