diff --git a/CHANGELOG.md b/CHANGELOG.md index 4280e772c..bcfe07c6b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/pyscipopt/scip.pxd b/src/pyscipopt/scip.pxd index b9bffc1d6..9f3110853 100644 --- a/src/pyscipopt/scip.pxd +++ b/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, diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi index 14741531d..d14371055 100644 --- a/src/pyscipopt/scip.pxi +++ b/src/pyscipopt/scip.pxi @@ -5730,6 +5730,13 @@ cdef class Model: 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 @@ cdef class Model: 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 = (pycons_initial).data PY_SCIP_CALL(SCIPreleaseCons(self._scip, &scip_cons)) return pycons @@ -8085,8 +8097,17 @@ cdef class Model: 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 @@ cdef class Model: 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 @@ cdef class Model: 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 @@ cdef class Model: 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)) diff --git a/tests/test_nonlinear.py b/tests/test_nonlinear.py index 383532f2e..2a5434bbd 100644 --- a/tests/test_nonlinear.py +++ b/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)