Remove constraints [wip]

ilpy/decl.pxd

@@ -93,7 +93,8 @@ cdef extern from "impl/solvers/SolverBackend.h":
         ) except +
         void setObjective(Objective&)
         void setConstraints(Constraints&)
-        void addConstraint(Constraint&)
+        int addConstraint(Constraint&)
+        void removeConstraint(int)
         void setTimeout(double)
         void setOptimalityGap(double, bool)
         void setNumThreads(unsigned int)

ilpy/impl/solvers/GurobiBackend.cpp

@@ -157,7 +157,7 @@ GurobiBackend::setConstraints(const Constraints &constraints)
 	GRB_CHECK(GRBupdatemodel(_model));
 }
 
-void
+intptr_t
 GurobiBackend::addConstraint(const Constraint& constraint) {
 
 	// set the linear coefficients
@@ -210,8 +210,13 @@ GurobiBackend::addConstraint(const Constraint& constraint) {
 	delete[] qrows;
 	delete[] qcols;
 	delete[] qvals;
+
+	return 0;
 }
 
+void GurobiBackend::removeConstraint(intptr_t constraintId){
+
+}
 bool
 GurobiBackend::solve(Solution& x, std::string& msg) {
 

ilpy/impl/solvers/GurobiBackend.h

@@ -53,7 +53,9 @@ class GurobiBackend : public SolverBackend {
 
 	void setConstraints(const Constraints& constraints);
 
-	void addConstraint(const Constraint& constraint);
+	intptr_t addConstraint(const Constraint& constraint);
+
+	void removeConstraint(intptr_t constraintId);
 
 	void setTimeout(double timeout) { _timeout = timeout; }
 

ilpy/impl/solvers/ScipBackend.cpp

@@ -137,7 +137,7 @@ ScipBackend::setConstraints(const Constraints& constraints) {
 	}
 }
 
-void
+intptr_t
 ScipBackend::addConstraint(const Constraint& constraint) {
 
 	// create a list of variables and their coefficients
@@ -189,6 +189,18 @@ ScipBackend::addConstraint(const Constraint& constraint) {
 	// we do not release the constraint here
 	// so that we can remove the constraints later in freeConstraints()
 	// SCIP_CALL_ABORT(SCIPreleaseCons(_scip, &c));
+
+	// Cast the pointer to the new SCIP constraint to intptr_t and return it
+	return reinterpret_cast<intptr_t>(c);
+}
+
+void ScipBackend::removeConstraint(intptr_t constraintId){
+	// Cast the intptr_t back to a SCIP_CONS pointer
+	SCIP_CONS *cons = reinterpret_cast<SCIP_CONS*>(constraintId);
+	// Remove the constraint from the model
+	SCIP_CALL_ABORT(SCIPdelCons(_scip, cons));
+	// Release the constraint
+	SCIP_CALL_ABORT(SCIPreleaseCons(_scip, &cons));
 }
 
 void

ilpy/impl/solvers/ScipBackend.h

@@ -52,7 +52,9 @@ class ScipBackend : public SolverBackend {
 
 	void setConstraints(const Constraints& constraints);
 
-	void addConstraint(const Constraint& constraint);
+	intptr_t addConstraint(const Constraint& constraint);
+
+	void removeConstraint(intptr_t constraintId);
 
 	void setTimeout(double timeout);
 

ilpy/impl/solvers/SolverBackend.h

@@ -6,10 +6,10 @@
 #include "Solution.h"
 #include "VariableType.h"
 
-class SolverBackend {
+class SolverBackend
+{
 
 public:
-
 	virtual ~SolverBackend() {}
 
 	/**
@@ -20,74 +20,76 @@ class SolverBackend {
 	 *                     Binary).
 	 */
 	virtual void initialize(
-			unsigned int numVariables,
-			VariableType variableType) = 0;
+		unsigned int numVariables,
+		VariableType variableType) = 0;
 
 	/**
 	 * Initialise the linear solver for the given type of variables.
 	 *
 	 * @param numVariables
 	 *             The number of variables in the problem.
-	 * 
+	 *
 	 * @param defaultVariableType
-	 *             The default type of the variables (Continuous, Integer, 
+	 *             The default type of the variables (Continuous, Integer,
 	 *             Binary).
 	 *
 	 * @param specialVariableTypes
-	 *             A map of variable numbers to variable types to override the 
+	 *             A map of variable numbers to variable types to override the
 	 *             default.
 	 */
 	virtual void initialize(
-			unsigned int                                numVariables,
-			VariableType                                defaultVariableType,
-			const std::map<unsigned int, VariableType>& specialVariableTypes) = 0;
+		unsigned int numVariables,
+		VariableType defaultVariableType,
+		const std::map<unsigned int, VariableType> &specialVariableTypes) = 0;
 
 	/**
 	 * Set the objective.
 	 *
 	 * @param objective A linear objective.
 	 */
-	virtual void setObjective(const Objective& objective) = 0;
+	virtual void setObjective(const Objective &objective) = 0;
 
 	/**
 	 * Set the linear (in)equality constraints.
 	 *
 	 * @param constraints A set of linear constraints.
 	 */
-	virtual void setConstraints(const Constraints& constraints) = 0;
+	virtual void setConstraints(const Constraints &constraints) = 0;
 
 	/**
 	 * Add a single constraint.
 	 *
 	 * @param constraint A linear constraints.
 	 */
-	virtual void addConstraint(const Constraint& constraint) = 0;
+	virtual intptr_t addConstraint(const Constraint &constraint) = 0;
+
+	virtual void removeConstraint(intptr_t constraintId) = 0;
 
 	/**
 	 * Set a timeout in seconds for subsequent solve calls.
 	 */
 	virtual void setTimeout(double timeout) = 0;
 
 	/**
-	 * Set the solver's optimality gap. The solver will terminate with an 
-	 * "optimal" solution as soon as the gap between the upper and lower bound 
+	 * Set the solver's optimality gap. The solver will terminate with an
+	 * "optimal" solution as soon as the gap between the upper and lower bound
 	 * is less than the given value times the upper bound.
 	 *
 	 * @param gap
 	 *             The optimality gap.
 	 *
 	 * @param absolute
-	 *             If set to true, a solution is considered optimal if the gap 
-	 *             between the upper and lower bound is smaller than the given 
+	 *             If set to true, a solution is considered optimal if the gap
+	 *             between the upper and lower bound is smaller than the given
 	 *             value.
 	 */
-	virtual void setOptimalityGap(double gap, bool absolute=false) = 0;
+	virtual void setOptimalityGap(double gap, bool absolute = false) = 0;
 
 	/**
 	 * Set the number of threads the solver can use.
 	 *
 	 * @param numThreads
-	 *             The number of threads to use. Defaults to 0, which leaves the 
+	 *             The number of threads to use. Defaults to 0, which leaves the
 	 *             decision to the solver.
 	 */
 	virtual void setNumThreads(unsigned int numThreads) = 0;
@@ -98,7 +100,7 @@ class SolverBackend {
 	 * @param verbose
 	 *             If set to true, verbose logging is enabled.
 	 */
-        virtual void setVerbose(bool verbose) = 0;
+	virtual void setVerbose(bool verbose) = 0;
 
 	/**
 	 * Solve the problem.
@@ -108,8 +110,7 @@ class SolverBackend {
 	 * @param message A status message from the solver.
 	 * @return true, if the optimal value was found.
 	 */
-	virtual bool solve(Solution& solution, std::string& message) = 0;
+	virtual bool solve(Solution &solution, std::string &message) = 0;
 };
 
 #endif // INFERENCE_SOLVER_BACKEND_H__
-

ilpy/wrapper.pyi

@@ -117,6 +117,7 @@ class Solver:
     def set_objective(self, objective: Objective | Expression) -> None: ...
     def set_constraints(self, constraints: Constraints) -> None: ...
     def add_constraint(self, constraint: Constraint | Expression) -> None: ...
+    def remove_constraint(self, constraint: Constraint | Expression) -> None: ...
     def set_timeout(self, timeout: float) -> None: ...
     def set_optimality_gap(self, gap: float, absolute: bool = False) -> None: ...
     def set_num_threads(self, num_threads: int) -> None: ...

ilpy/wrapper.pyx

@@ -4,6 +4,7 @@ from typing import TYPE_CHECKING
 from libc.stdint cimport uint32_t
 from libcpp.memory cimport shared_ptr
 from libcpp.map cimport map as cppmap
+from libc.stdint cimport uintptr_t
 from libcpp.string cimport string
 from cython.operator cimport dereference as deref
 from . cimport decl
@@ -245,6 +246,7 @@ cdef class Solver:
 
     cdef shared_ptr[decl.SolverBackend] p
     cdef unsigned int num_variables
+    cdef dict _constraint_map
 
     def __cinit__(
             self,
@@ -259,6 +261,7 @@ cdef class Solver:
                 vtypes[k] = v
         self.p = factory.createSolverBackend(preference)
         self.num_variables = num_variables
+        self._constraint_map = {}
         deref(self.p).initialize(num_variables, default_variable_type, vtypes)
 
     def set_objective(self, objective: Objective | Expression):
@@ -278,7 +281,14 @@ cdef class Solver:
             const = constraint.as_constraint()
         else:
             const = constraint
-        deref(self.p).addConstraint(const.p[0])
+        backend_id = <uintptr_t>deref(self.p).addConstraint(const.p[0])
+        self._constraint_map[id(constraint)] = backend_id
+
+    def remove_constraint(self, constraint: Constraint | Expression):
+        if id(constraint) not in self._constraint_map:
+            raise ValueError("Constraint not found")
+        backend_id = <uintptr_t>self._constraint_map[id(constraint)]
+        deref(self.p).removeConstraint(backend_id)
 
     def set_timeout(self, timeout):
         deref(self.p).setTimeout(timeout)

tests/test_solvers.py

@@ -228,3 +228,29 @@ def test_solve_twice(preference: ilpy.Preference) -> None:
     solver.set_constraints(c1)
     solution = solver.solve()
     assert list(solution) != [7, 3]
+
+
+@pytest.mark.parametrize("preference", [ilpy.Scip])
+def test_remove_constraint(preference: ilpy.Preference) -> None:
+    solver = ilpy.Solver(2, ilpy.VariableType.Integer, preference=preference)
+    x = ilpy.Variable("x", index=0)
+    y = ilpy.Variable("y", index=1)
+    solver.set_objective((x + y).as_objective(ilpy.Minimize))
+    c1 = y == x + 1
+    c2 = y == 3 - x
+    c3 = y == 5 - x  # incompatible with c2
+    solver.add_constraint(c1)
+    solver.add_constraint(c2)
+    solution = solver.solve()
+    assert list(solution) == [1, 2]
+
+    # FIXME: if we add an unsatisfiable constraint and solve in the middle,
+    # it appears to leave the solver in an unrecoverable state.
+    # solver.add_constraint(c3)
+    # solution = solver.solve()
+    # assert list(solution) == [0, 0]  # optimal solution is now infeasible
+
+    solver.remove_constraint(c2)
+    solver.add_constraint(c3)
+    solution = solver.solve()
+    assert list(solution) == [2, 3]