Merge branch 'master' into feature-gil

Author: Joao-Dionisio

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ## Unreleased
 ### Added
+- Added stage checks to presolve, freereoptsolve, freetransform
+- 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

@@ -2101,6 +2101,15 @@ cdef class Model:
     def freeTransform(self):
         """Frees all solution process data including presolving and
         transformed problem, only original problem is kept."""
+        if self.getStage() not in [SCIP_STAGE_INIT,
+                                 SCIP_STAGE_PROBLEM,
+                                 SCIP_STAGE_TRANSFORMED,
+                                 SCIP_STAGE_PRESOLVING,
+                                 SCIP_STAGE_PRESOLVED,
+                                 SCIP_STAGE_SOLVING,
+                                 SCIP_STAGE_SOLVED]:
+            raise Warning("method cannot be called in stage %i." % self.getStage())
+            
         self._modelvars = {
             var: value
             for var, value in self._modelvars.items()
@@ -6227,6 +6236,11 @@ cdef class Model:
 
     def presolve(self):
         """Presolve the problem."""
+        if self.getStage() not in [SCIP_STAGE_PROBLEM, SCIP_STAGE_TRANSFORMED,\
+                                SCIP_STAGE_PRESOLVING, SCIP_STAGE_PRESOLVED, \
+                                SCIP_STAGE_SOLVED]:
+            raise Warning("method cannot be called in stage %i." % self.getStage())
+
         PY_SCIP_CALL(SCIPpresolve(self._scip))
         self._bestSol = Solution.create(self._scip, SCIPgetBestSol(self._scip))
 
@@ -9029,6 +9043,15 @@ cdef class Model:
 
     def freeReoptSolve(self):
         """Frees all solution process data and prepares for reoptimization."""
+
+        if self.getStage() not in [SCIP_STAGE_INIT,
+                                 SCIP_STAGE_PROBLEM,
+                                 SCIP_STAGE_TRANSFORMED,
+                                 SCIP_STAGE_PRESOLVING,
+                                 SCIP_STAGE_PRESOLVED,
+                                 SCIP_STAGE_SOLVING,
+                                 SCIP_STAGE_SOLVED]:
+            raise Warning("method cannot be called in stage %i." % self.getStage())
         PY_SCIP_CALL(SCIPfreeReoptSolve(self._scip))
 
     def chgReoptObjective(self, coeffs, sense = 'minimize'):