getTermsQuadractic correctly returns linear terms

CHANGELOG.md

@@ -2,8 +2,10 @@
 
 ## Unreleased
 ### Added
+- Added getVarFarkasCoef() (untested)
 ### Fixed
 - all fundamental callbacks now raise an error if not implemented
+- getTermsQuadratic() now correctly returns all linear terms
 ### Changed
 ### Removed
 

src/pyscipopt/scip.pxd

@@ -1001,6 +1001,7 @@ cdef extern from "scip/scip.h":
     SCIP_Real SCIPgetDualboundRoot(SCIP* scip)
     SCIP_Real SCIPgetVarRedcost(SCIP* scip, SCIP_VAR* var)
     SCIP_RETCODE SCIPgetDualSolVal(SCIP* scip, SCIP_CONS* cons, SCIP_Real* dualsolval, SCIP_Bool* boundconstraint)
+    SCIP_Real SCIPgetVarFarkasCoef(SCIP* scip, SCIP_VAR* var)
 
     # Reader plugin
     SCIP_RETCODE SCIPincludeReader(SCIP* scip,

src/pyscipopt/scip.pxi

@@ -316,7 +316,7 @@
     if rc == SCIP_OKAY:
         pass
     elif rc == SCIP_ERROR:
         raise Exception('SCIP: unspecified error!')
     elif rc == SCIP_NOMEMORY:
         raise MemoryError('SCIP: insufficient memory error!')
     elif rc == SCIP_READERROR:
@@ -335,7 +335,7 @@
         raise Exception('SCIP: method cannot be called at this time'
                             + ' in solution process!')
     elif rc == SCIP_INVALIDDATA:
         raise Exception('SCIP: error in input data!')
     elif rc == SCIP_INVALIDRESULT:
         raise Exception('SCIP: method returned an invalid result code!')
     elif rc == SCIP_PLUGINNOTFOUND:
@@ -5730,6 +5730,13 @@
 
         PyCons = Constraint.create(scip_cons)
 
+        # Store the original polynomial expression on the constraint so that
+        # helpers such as getTermsQuadratic can reconstruct full linear terms
+        # even if SCIP's internal quadratic representation does not expose
+        # all linear coefficients explicitly.
+        if PyCons.data is None:
+            PyCons.data = quadcons.expr
+
         return PyCons
 
     def _createConsNonlinear(self, cons, **kwargs):
@@ -6064,6 +6071,11 @@
 
         PY_SCIP_CALL(SCIPaddCons(self._scip, scip_cons))
         pycons = Constraint.create(scip_cons)
+        # propagate any problem data (such as the original Expr for
+        # expression-based constraints) from the temporary constraint
+        # created in createConsFromExpr to the final Constraint object
+        # that is returned to the user
+        pycons.data = (<Constraint>pycons_initial).data
         PY_SCIP_CALL(SCIPreleaseCons(self._scip, &scip_cons))
 
         return pycons
@@ -8085,8 +8097,17 @@
         Returns
         -------
         bilinterms : list of tuple
+            Triples ``(var1, var2, coef)`` for terms of the form
+            ``coef * var1 * var2`` with ``var1 != var2``.
         quadterms : list of tuple
+            Triples ``(var, sqrcoef, lincoef)`` corresponding to diagonal
+            quadratic terms of the form ``sqrcoef * var**2`` and the linear
+            coefficient ``lincoef`` associated with the same variable when it
+            also appears linearly in the quadratic part.
         linterms : list of tuple
+            Pairs ``(var, coef)`` for all variables with a nonzero linear
+            coefficient in the constraint, including variables that also
+            appear in quadratic or bilinear terms.
 
         """
         cdef SCIP_EXPR* expr
@@ -8117,15 +8138,36 @@
         assert self.checkQuadraticNonlinear(cons), "constraint is not quadratic"
 
         expr = SCIPgetExprNonlinear(cons.scip_cons)
-        SCIPexprGetQuadraticData(expr, NULL, &nlinvars, &linexprs, &lincoefs, &nquadterms, &nbilinterms, NULL, NULL)
+        SCIPexprGetQuadraticData(expr, NULL, &nlinvars, &linexprs, &lincoefs,
+                                 &nquadterms, &nbilinterms, NULL, NULL)
 
         linterms   = []
         bilinterms = []
         quadterms  = []
 
-        for termidx in range(nlinvars):
-            var = Variable.create(SCIPgetVarExprVar(linexprs[termidx]))
-            linterms.append((var, lincoefs[termidx]))
+        # First try to recover all linear coefficients from the original
+        # polynomial expression, if it has been stored on the Constraint.
+        if isinstance(cons.data, Expr):
+            lindict = {}
+            for term, coef in cons.data.terms.items():
+                if coef == 0.0:
+                    continue
+                if len(term) == 1:
+                    var = term[0]
+                    key = var.ptr()
+                    if key in lindict:
+                        _, oldcoef = lindict[key]
+                        lindict[key] = (var, oldcoef + coef)
+                    else:
+                        lindict[key] = (var, coef)
+            for _, (var, coef) in lindict.items():
+                linterms.append((var, coef))
+        else:
+            # use only the purely linear part as exposed by SCIP's quadratic representation
+            for termidx in range(nlinvars):
+                scipvar1 = SCIPgetVarExprVar(linexprs[termidx])
+                var = Variable.create(scipvar1)
+                linterms.append((var, lincoefs[termidx]))
 
         for termidx in range(nbilinterms):
             SCIPexprGetQuadraticBilinTerm(expr, termidx, &bilinterm1, &bilinterm2, &bilincoef, NULL, NULL)
@@ -8137,13 +8179,13 @@
                 bilinterms.append((var1,var2,bilincoef))
             else:
                 quadterms.append((var1,bilincoef,0.0))
-
         for termidx in range(nquadterms):
             SCIPexprGetQuadraticQuadTerm(expr, termidx, NULL, &lincoef, &sqrcoef, NULL, NULL, &sqrexpr)
             if sqrexpr == NULL:
                 continue
-            var = Variable.create(SCIPgetVarExprVar(sqrexpr))
-            quadterms.append((var,sqrcoef,lincoef))
+            scipvar1 = SCIPgetVarExprVar(sqrexpr)
+            var = Variable.create(scipvar1)
+            quadterms.append((var, sqrcoef, lincoef))
 
         return (bilinterms, quadterms, linterms)
 
@@ -8488,6 +8530,31 @@
 
         return _dualsol
 
+    def getVarFarkasCoef(self, Variable var):
+        """
+        Returns the Farkas coefficient of the variable in the current node's LP relaxation; the current node has to have an infeasible LP.
+
+        Parameters
+        ----------
+        var : Variable
+            variable to get the farkas coefficient of
+
+        Returns
+        -------
+        float
+
+        """
+        assert SCIPgetStatus(self._scip) == SCIP_STATUS_INFEASIBLE, "Method can only be called with an infeasible model."
+
+        farkas_coef = None
+        try:
+            farkas_coef = SCIPgetVarFarkasCoef(self._scip, var.scip_var)
+            if self.getObjectiveSense() == "maximize":
+                farkas_coef = -farkas_coef
+        except Exception:
+            raise Warning("no farkas coefficient available for variable " + var.name)
+        return farkas_coef
+
     def optimize(self):
         """Optimize the problem."""
         PY_SCIP_CALL(SCIPsolve(self._scip))

tests/test_nonlinear.py

@@ -288,6 +288,52 @@ def test_quad_coeffs():
     assert linterms[0][0].name == z.name
     assert linterms[0][1] == 4
 
+
+def test_quad_coeffs_mixed_linear_and_quadratic():
+
+    scip = Model()
+
+    var1 = scip.addVar(name="var1", vtype='C', lb=None)
+    var2 = scip.addVar(name="var2", vtype='C')
+    var3 = scip.addVar(name="var3", vtype='B')
+    var4 = scip.addVar(name="var4", vtype='B')
+
+    cons = scip.addCons(
+        8 * var4
+        + 4 * var3
+        - 5 * var2
+        + 6 * var3 ** 2
+        - 3 * var1 ** 2
+        + 2 * var2 * var1
+        + 7 * var1 * var3
+        == -2
+    )
+
+    bilinterms, quadterms, linterms = scip.getTermsQuadratic(cons)
+
+    # Linear part: getTermsQuadratic must report all linear coefficients,
+    # including those of variables that also appear quadratically or
+    # bilinearly.
+    lin_only = {v.name: c for (v, c) in linterms}
+    assert lin_only["var4"] == 8
+    assert lin_only["var3"] == 4
+    assert lin_only["var2"] == -5
+    # var1 has no linear component and must not appear in linterms
+    assert "var1" not in lin_only or lin_only.get("var1", 0.0) == 0.0
+
+    # For completeness, checking if the coefficients from reconstructing the full linear
+    # coefficients from both linterms and quadterms match
+    full_lin = {}
+    for v, c in linterms:
+        full_lin[v.name] = full_lin.get(v.name, 0.0) + c
+    for v, _, lincoef in quadterms:
+        if lincoef != 0.0:
+            full_lin[v.name] = full_lin.get(v.name, 0.0) + lincoef
+
+    assert full_lin["var4"] == 8
+    assert full_lin["var3"] == 4
+    assert full_lin["var2"] == -5
+
 def test_addExprNonLinear():
     m = Model()
     x = m.addVar("x", lb=0, ub=1, obj=10)