nlhdlr_quadratic.c
Go to the documentation of this file.
33 * - a `QUADEXPRTERM` stores an expression `expr` that is known to appear in a nonlinear, quadratic term, that is
34 * `expr^2` or `expr*other_expr`. It stores its `sqrcoef` (that can be 0), its linear coef and all the bilinear expression
38/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
50#define INTERLOG(x) if( SCIPgetSubscipDepth(scip) == 0 && SCIPgetVerbLevel(scip) >= SCIP_VERBLEVEL_NORMAL ) { x }
74#define TABLE_EARLIEST_STAGE_QUADRATIC SCIP_STAGE_TRANSFORMED /**< output of the statistics table is only printed from this stage onwards */
101 SCIP_Real minquadfiniteact; /**< minimum activity of quadratic part where only terms with finite min
103 SCIP_Real maxquadfiniteact; /**< maximum activity of quadratic part where only terms with finite max
107 SCIP_INTERVAL* quadactivities; /**< activity of each quadratic term as defined in nlhdlrIntevalQuadratic */
122 int ncutsadded; /**< total number of cuts that where generated by separateQuadratic and actually added */
123 SCIP_Longint lastnodenumber; /**< number of last node for which cuts were (allowed to be) generated */
127 SCIP_Bool useintersectioncuts; /**< whether to use intersection cuts for quadratic constraints or not */
130 SCIP_Bool useminrep; /**< whether the minimal representation of the S-free set should be used (instead of the gauge) */
131 SCIP_Bool useboundsasrays; /**< use bounds of variables in quadratic as rays for intersection cuts */
135 SCIP_Real mincutviolation; /**< minimal cut violation the generated cuts must fulfill to be added to the LP */
136 SCIP_Real minviolation; /**< minimal violation the constraint must fulfill such that a cut can be generated */
137 int atwhichnodes; /**< determines at which nodes cut is used (if it's -1, it's used only at the root node,
141 SCIP_Bool ignorebadrayrestriction; /**< should cut be generated even with bad numerics when restricting to ray? */
143 SCIP_Bool trackmore; /**< for monoidal strengthening, should we track more statistics (more expensive) */
146 int ncouldimprovedcoef; /**< number of times a coefficient could improve but didn't because of numerics */
147 int nbadrayrestriction; /**< number of times a cut was aborted because of numerics when restricting to ray */
148 int nbadnonbasic; /**< number of times a cut was aborted because the nonbasic row was not nonbasic enough */
157 SCIP_Real monoidalimprovementsum; /**< sum of average improvement of a cut when using monoidal strengthening */
160 SCIP_Real currentavemonoidalimprovement;/**< average improvement of current cut when using monoidal strengthening */
198 SCIPinfoMessage(scip, file, "Quadratic Nlhdlr : %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %20s %10s %10s %10s \n", "GenCuts", "AddCuts", "CouldImpr", "NLargeRE",
199 "AbrtBadRay", "AbrtPosPhi", "AbrtNonBas", "NStrength", "NMonoidal", "AveCutcoef", "AveMonoidalImprov", "AveDensity", "AveEfficacy", "AveBCutsFrac");
210 SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncutsgenerated > 0 ? nlhdlrdata->cutcoefsum / nlhdlrdata->ncutsgenerated : 0.0);
211 SCIPinfoMessage(scip, file, " %20g", (nlhdlrdata->nmonoidal > 0 && nlhdlrdata->trackmore) ? nlhdlrdata->monoidalimprovementsum / nlhdlrdata->nmonoidal : -1.0);
212 SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncutsgenerated > 0 ? nlhdlrdata->densitysum / nlhdlrdata->ncutsgenerated : 0.0);
213 SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncutsgenerated > 0 ? nlhdlrdata->efficacysum / nlhdlrdata->ncutsgenerated : 0.0);
214 SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncalls > 0 ? nlhdlrdata->nboundcuts / nlhdlrdata->ncalls : 0.0);
238 SCIPinfoMessage(scip, NULL, "adding col %s to cut. %g <= col <= %g\n", SCIPvarGetName(SCIPcolGetVar(col)),
240 SCIPinfoMessage(scip, NULL, "col is active at %s. Value %.15f\n", SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ? "lower bound" :
282 SCIPinfoMessage(scip, NULL, "adding slack var row_%d to cut. %g <= row <= %g\n", SCIProwGetLPPos(row), SCIProwGetLhs(row), SCIProwGetRhs(row));
283 SCIPinfoMessage(scip, NULL, "row is active at %s = %.15f Activity %.15f\n", SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? "lhs" :
284 "rhs" , SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? SCIProwGetLhs(row) : SCIProwGetRhs(row),
313 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(rowcols[i]), -rowcoefs[i] * cutcoef) );
366 SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
423 * However, we want to store the the tableau row by columns. Thus, we need to know which of the basic vars `col` is.
463 SCIP_CALL( SCIPgetLPBInvARow(scip, basicvarpos2tableaurow[lppos], binvrow, binvarow, NULL, NULL) );
482/** stores the rows of the tableau corresponding to the basic variables in the quadratic expression
484 * Also return a map storing to which var the entry of a ray corresponds, i.e., if the tableau is
491 * The map maps k to the position of basicvar_k in the variables of the constraint assuming the variables are sorted as
529 SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
541 SCIP_CALL( storeDenseTableauRow(scip, col, basicvarpos2tableaurow, nrayentries, raylength, binvrow, binvarow,
554 SCIP_CALL( storeDenseTableauRow(scip, col, basicvarpos2tableaurow, nrayentries, raylength, binvrow, binvarow,
567 SCIP_CALL( storeDenseTableauRow(scip, col, basicvarpos2tableaurow, nrayentries, raylength, binvrow, binvarow,
607 SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->raysbegin, SCIPgetNLPCols(scip) + SCIPgetNLPRows(scip) + 1) );
608 SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->lpposray, SCIPgetNLPCols(scip) + SCIPgetNLPRows(scip)) );
686/** constructs map between the lppos of a variables and its position in the constraint assuming the constraint variables
706 SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
766 * TODO: in case of problems, an idea would be to scale the ray entries; compute the cut coef and scale it back down
768 * The problem is that if the cut coefficient is 1/t where lpsol + t*ray intersects the S-free set.
813 SCIPinfoMessage(scip, NULL, "entries of ray %d are between [%d, %d):\n", rays->nrays, rays->raysbegin[rays->nrays], *nnonz);
841 * Also, we store the rays as if every nonbasic variable was at lower (so that all rays moves to infinity)
854 * In constrast, the nonbasic part of the ray has a 1.0 for nonbasic at lower and a -1.0 for nonbasic at upper, i.e.
895 /* construct dense tableau and map between ray entries and position of corresponding var in quad cons */
933 SCIP_CALL( insertRayEntries(scip, rays, densetableaucols, rayentry2conspos, raylength, nnonbasic,
937 SCIPinfoMessage(scip, NULL, "looked at ray of var %s with basestat %d, it has %d nonzeros\n-----------------\n",
943 SCIPdebugMsg(scip, "nonzero ray associated with variable <%s> has base status zero -> abort storing rays\n",
969 /* set factor to store basic entries of ray as = [-BinvL, BinvU]; basic status of rows are flipped! See lpi.h! */
978 SCIP_CALL( insertRayEntries(scip, rays, densetableaucols, rayentry2conspos, raylength, nnonbasic, -1, factor,
984 SCIPinfoMessage(scip, NULL, "looked at ray of row %d, it has %d nonzeros\n-----------------\n", i, nnonz - oldnnonz);
1033 * To do this transformation and later to compute the actual cut we need to compute and store some quantities.
1035 * - \f$I_0\f$, \f$I_+\f$, and \f$I_-\f$ be the index set of zero, positive, and negative eigenvalues, respectively
1048 * @note if the constraint is q(z) ≤ rhs, then the constant when writing the constraint as quad ≤ 0 is c - rhs.
1049 * @note if the quadratic constraint we are separating is q(z) ≥ lhs, then we multiply by -1.
1086 SCIPexprGetQuadraticData(qexpr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, &eigenvalues,
1126 vdotzlp += SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr)) * eigenvectors[offset + j];
1156 *wzlp += (sidefactor * lincoefs[i]) * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(linexprs[i]));
1191 /* rays are sorted; the first entries correspond to the quad vars, so we stop after first nonquad var entry */
1230 SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, NULL, &lincoefs, &nquadexprs, NULL, NULL, NULL);
1245 /* rays are sorted; last entries correspond to the lin vars, so we stop after first quad var entry */
1251 lincoefs[rayidx[i] - nquadexprs]) * raycoefs[i], lincoefs[rayidx[i] - nquadexprs] ,raycoefs[i]);
1259/** computes the dot product of v_i and the current ray as well as of v_i and the apex where v_i
1280 SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
1320 * The restriction of the function representing the maximal S-free set to zlp + t * ray has the form
1325 * In case 4, it computes the coefficients for both pieces, in addition to coefficients needed to evaluate the condition
1363 SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
1421 /* In theory, the function at 0 must be negative. Because of bad numerics this might not always hold, so we abort
1435 SCIPinfoMessage(scip, NULL, "Restriction yields case 2: a,b,c,d,e %g %g %g %g %g\n", coefs1234a[0], coefs1234a[1], coefs1234a[2],
1453 * The restriction of the function representing the maximal S-free set to zlp + t * ray has the form
1458 * In case 4, it computes the coefficients for both pieces, in addition to coefficients needed to evaluate the condition
1478 SCIP_Real* coefs4b, /**< buffer to store A, B, C, D, and E of case 4b (or NULL if not needed) */
1479 SCIP_Real* coefscondition, /**< buffer to store data to evaluate condition to decide case 4a or 4b */
1498 SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
1545 vdotray = computeEigenvecDotRay(&eigenvectors[i * nquadexprs], nquadexprs, raycoefs, rayidx, raynnonz);
1563 printf("Positive eigenvalue: computing D: v^T ray %g, v^T( zlp + b/theta ) %g and theta %g \n", vdotray, dot, (sidefactor * eigenvalues[i]));
1574 printf("Negative eigenvalue: computing A: v^T ray %g, and theta %g \n", vdotray, (sidefactor * eigenvalues[i]));
1593 /* In theory, the function at 0 must be negative. Because of bad numerics this might not always hold, so we abort
1659 * coefficients of condition: stores -numerator of x_{r+1}/ norm xhat, w(ray), and numerator of y_{s+1} at zlp
1671 SCIPinfoMessage(scip, NULL, "Restriction yields case 1,2 or 3: a,b,c,d,e %g %g %g %g %g\n", coefs1234a[0], coefs1234a[1], coefs1234a[2],
1676 SCIPinfoMessage(scip, NULL, "Restriction yields\n Case 4a: a,b,c,d,e %g %g %g %g %g\n", coefs1234a[0],
1678 SCIPinfoMessage(scip, NULL, " Case 4b: a,b,c,d,e %g %g %g %g %g\n", coefs4b[0], coefs4b[1], coefs4b[2],
1680 SCIPinfoMessage(scip, NULL, " Condition: xextra/e, wray, yextra %g %g %g g\n", coefscondition[0],
1694 /*assert(4.0 * (*a) * (*c) >= SQR( *b ) ); *//* the function is defined everywhere, so minimum of radicand must be nonnegative */
1758 * \f[ -\lambda_{r+1} \Vert \hat y(zlp + tsol\, ray)\Vert + \hat y_{s+1}(zlp + tsol\, ray) \leq 0\f]
1761 * \f[ -num(\hat x_{r+1}(zlp)) \sqrt{A t^2 + B t + C} / E + w(ray) \cdot t + num(\hat y_{s+1}(zlp)) \leq 0\f]
1772 return (coefscondition[0] * sqrt(coefs4a[0] * SQR( tsol ) + coefs4a[1] * tsol + coefs4a[2]) + coefscondition[1] *
1800 printf("%d: lb,ub %.10f, %.10f. curr = %g -> phi at curr %g -> phi at lb %g \n", i, lb, ub, curr, phival, evalPhiAtRay(lb, a, b, c, d, e));
1837 /* there is an intersection point if and only if sqrt(A) > D: here we are beliving in math, this might cause
1848 /* SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar finds all x such that a x^2 + b x >= c and x in bounds.
1849 * it is known that if tsol is the root we are looking for, then gamma(zlp + t * ray) <= 0 between 0 and tsol, thus
1852 SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a - d * d, b - 2.0 * d *
1856 sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : SCIPintervalGetInf(result);
1867 /* check that solution is acceptable, ideally it should be <= 0, however when it is positive, we trigger a binary
1868 * search to make it negative. This binary search might return a solution point that is not at accurately 0 as the
1869 * one obtained from the function above. Thus, it might fail to satisfy the condition of case 4b in some cases, e.g.,
1875 printf("interval solution returned %g -> phival = %g, believe it\n", sol, evalPhiAtRay(sol, a, b, c, d, e));
1882 /* perform a binary search to make it negative: this might correct a wrong infinity (e.g. crudeoil_lee1_05) */
1892/** The maximal S-free set is \f$\gamma(z) \leq 0\f$; we find the intersection point of the ray `ray` starting from zlp with the
1907 * It can be shown (given the special properties of \f$\gamma\f$) that the smallest positive root of each function of the form
1911 * & \sqrt{a t^2 + b t + c} - (d t + e)) (\sqrt{a t^2 + b t + c} + (d t + e)) = 0 \\ \Leftrightarrow
1917 * If there is no solution, then the second piece can't have a solution (first piece ≥ second piece for all t)
1971 /* not on 4a --> then the intersection point is whatever 4b says: as phi4a >= phi4b, the solution of phi4b should
1972 * always be larger (but shouldn't be equal at this point given that isCase4a failed, and the condition function
1973 * evaluates to 0 when phi4a == phi4b) than the solution of phi4a; However, because of numerics (or limits in the
1981 if( sol4b < sol1234a && evalPhiAtRay(1.1 * sol1234a, coefs4b[0], coefs4b[1], coefs4b[2], coefs4b[3], coefs4b[4]) <=
2002 /* check at phi at 0 is negative (note; this could be checked before restricting to the ray) also, if this
2030 INTERLOG(printf("Bad numerics 1 2 3 or 4a: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min); )
2064/** computes the coefficients a, b, c defining the quadratic function defining set S restricted to the line
2067 * The solution to the monoidal strengthening problem is then given by the smallest root of the function
2092 SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2117/** check if ray was in strip by checking if the point in the monoid corresponding to the cutcoef we just found
2141 SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2154 num += sidefactor * eigenvalues[i] * dot * (cutcoef * (vzlp[i] - vapex[i]) + vray[i] + vapex[i]);
2195 SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a, b, -c, bounds);
2196 sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : SCIPintervalGetInf(result);
2204 SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a, -b, -c, bounds);
2205 sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : -SCIPintervalGetInf(result);
2209 /* check if that solution is close enough or if we need to improve it more with binary search */
2278 SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2310 /* if denom = 0, the S-free set is just the strip, so we don't have an apex -> monoidal not possible */
2366 SCIPexprGetQuadraticData(nlhdlrexprdata->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
2386 /* check if ray is in strip. If not, monoidal is possible and cutcoef is the strengthened cut coef */
2400/** sparsify intersection cut by replacing non-basic variables with their bounds if their coefficient allows it */
2433 /* if the variable is at its lower or upper bound and the coefficient has the correct sign, we can
2451/** computes intersection cut cuts off sol (because solution sol violates the quadratic constraint cons)
2468 SCIP_Real* interpoints, /**< array to store intersection points for all rays or NULL if nothing
2497 /* if we use monoidal and we are in the right case for it, compute the apex of the S-free set */
2502 SCIPexprGetQuadraticData(nlhdlrexprdata->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
2537 SCIP_CALL( computeMonoidalStrengthCoef(scip, nlhdlrexprdata, rays->lpposray[i], &rays->rays[rays->raysbegin[i]],
2538 &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] - rays->raysbegin[i], vb, vzlp, wcoefs, kappa,
2552 rays->raysbegin[i], vb, vzlp, wcoefs, wzlp, kappa, coefs1234a, coefs4b, coefscondition, success) );
2564 interpoint = computeIntersectionPoint(scip, nlhdlrdata, iscase4, coefs1234a, coefs4b, coefscondition);
2574 /* if we are only computing the "normal" cut coefficient to track statistics (and we have been successful
2629 assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
2648 assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
2705 /* if the pointers look at different variables (or one already arrieved at the end), only one pointer can move
2724 /* if both pointers look at the same variable, just compute the difference and move both pointers */
2784/** checks if the ray alpha * ray_i + (1 - alpha) * ray_j is in the recession cone of the S-free set. To do so,
2818 newraynnonz = (rays->raysbegin[i + 1] - rays->raysbegin[i]) + (rays->raysbegin[j + 1] - rays->raysbegin[j]);
2823 combineRays(&rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] -
2829 SCIP_CALL( computeRestrictionToRay(scip, nlhdlrexprdata, sidefactor, iscase4, newraycoefs, newrayidx,
2835 /* check if restriction to "new" ray is numerically nasty. If so, treat the corresponding rho as if phi is
2840 interpoint = computeIntersectionPoint(scip, nlhdlrdata, iscase4, coefs1234a, coefs4b, coefscondition);
2870 SCIP_Real* interpoints, /**< array to store intersection points for all rays or NULL if nothing
2878 /* go through all rays not in the recession cone and compute the largest negative steplength possible. The
2897 if( raysAreDependent(scip, &rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]],
2899 &rays->raysidx[rays->raysbegin[idx]], rays->raysbegin[idx + 1] - rays->raysbegin[idx], &coef) )
2908 * Since we know that we can only use alpha < maxalpha, we don't need to do the whole binary search
2909 * for every ray i. We only need to search the intervall [0, maxalpha]. Thereby, we start by checking
2930 SCIP_CALL( rayInRecessionCone(scip, nlhdlrdata, nlhdlrexprdata, rays, idx, i, sidefactor, iscase4, vb,
2947 /* now we found the best convex combination which we use to derive the corresponding coef. If alpha = 0, we
2972/** computes intersection cut using negative edge extension to strengthen rays that do not intersect
3010 /* compute all intersection points and store them in interpoints; build not-stregthened intersection cut */
3011 SCIP_CALL( computeIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3021 /* go through all intersection points that are equal to infinity -> these correspond to the rays which are in the
3022 * recession cone of C, i.e. the rays for which we (possibly) can compute a negative steplength */
3037 SCIP_CALL( findRho(scip, nlhdlrdata, nlhdlrexprdata, rays, i, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3059 assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
3062 SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef :
3074 assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
3076 SCIP_CALL( addColToCut(scip, rowprep, sol, SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef :
3108 if( SCIPisInfinity(scip, SCIPvarGetLbLocal(var)) && SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) )
3129/** This function finds vertex (w.r.t. bounds of variables appearing in the quadratic) that is closest to the current
3158 SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
3191 rays->lpposray[i + nquadexprs] = SCIPcolGetLPPos(SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(linexprs[i])));
3207 SCIP_CALL( setVarToNearestBound(scip, sol, vertex, auxvar, &rays->rays[nquadexprs + nlinexprs], success) );
3302/** generates intersection cut that cuts off sol (which violates the quadratic constraint cons) */
3339 SCIPinfoMessage(scip, NULL, "Generating intersection cut for quadratic expr %p aka", (void*)expr);
3414 SCIP_CALL( intercutsComputeCommonQuantities(scip, nlhdlrexprdata, auxvar, sidefactor, soltoseparate,
3437 SCIP_CALL( computeStrengthenedIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs,
3445 SCIP_CALL( computeIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3464 * It is propagable, if a variable (aka child expr) appears at least twice, which is the case if at least two of the following hold:
3488 if( (lincoef != 0.0) + (sqrcoef != 0.0) + nadjbilin >= 2 ) /*lint !e514*/ /* actually MIN(2, nadjbilin), but we check >= 2 */
3497 * A term is propagable, if its variable (aka child expr) appears at least twice, which is the case if at least two of the following hold:
3515 return (lincoef != 0.0) + (sqrcoef != 0.0) + nadjbilin >= 2; /*lint !e514*/ /* actually MIN(2, nadjbilin), but we check >= 2 */
3518/** solves a quadratic equation \f$ a\, \text{expr}^2 + b\, \text{expr} \in \text{rhs} \f$ (with \f$b\f$ an interval)
3542 SCIPinfoMessage(scip, NULL, "Propagating <expr> by solving a <expr>^2 + b <expr> in rhs, where <expr> is: ");
3560 SCIPintervalSolveUnivariateQuadExpression(SCIP_INTERVAL_INFINITY, &newrange, a, b, rhs, exprbounds);
3571/** solves a linear equation \f$ b\, \text{expr} \in \text{rhs} \f$ (with \f$b\f$ a scalar) and reduces bounds on `expr` or deduces infeasibility if possible */
3590 SCIPinfoMessage(scip, NULL, "Propagating <expr> by solving %g <expr> in [%g, %g], where <expr> is: ", b, rhs.inf, rhs.sup);
3651 * if c = 0, then the function is monotone which means the maximum is also at one of the boundaries
3653 * if a < 0, then the function is concave. The function then has a maximum if and only if there is a point with derivative 0,
3654 * that is, iff -a/x^2 - c = 0 has a solution; i.e. if -a/c >= 0, i.e. (using a<0 and c != 0), c > 0.
3669 * the (restricted) maximum is at a boundary (we could even say at which boundary, but that doesn't save much)
3674 /* the maximum at sqrt(-a/c) is -2*sqrt(-a*c), so we compute an upper bound for that by computing a lower bound for 2*sqrt(-a*c) */
3680 /* if the interval containing sqrt(-a/c) is contained in dom, then we can return -negunresmax */
3684 /* now what is left is the case where we cannot say for sure whether sqrt(-a/c) is contained in dom or not
3692/** computes the range of rhs/x - coef * x for x in exprdom; this is used for the propagation of bilinear terms
3694 * If 0 is in the exprdom, we set range to \f$\mathbb{R}\f$ (even though this is not quite correct, it is correct for the
3696 * TODO: maybe check before calling it whether 0 is in the domain and then just avoid calling it
3698 * If rhs is [A,B] and x > 0, then we want the min of A/x - coef*x and max of B/x - coef*x for x in [exprdom].
3699 * If rhs is [A,B] and x < 0, then we want the min of B/x - coef*x and max of A/x - coef*x for x in [exprdom].
3700 * However, this is the same as min of -B/x + coef*x and max of -A/x + coef*x for x in -[exprdom].
3746 SCIP_Bool* infeasible, /**< buffer to store whether an exps' bounds were propagated to an empty interval */
3761 oldboundslin[i] = SCIPexprGetActivity(linexprs[i]); /* TODO use SCIPgetExprBoundsNonlinear(scip, linexprs[i]) ? */
3768 /* SCIP is more conservative with what constitutes a reduction than interval arithmetic so we follow SCIP */
3772 SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, linexprs[i], newboundslin[i], infeasible, nreductions) );
3797 SCIPexprGetQuadraticData((*nlhdlrexprdata)->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
3812 * We define a _propagable_ quadratic expression as a quadratic expression whose termwise propagation does not yield the
3813 * best propagation. In other words, is a quadratic expression that suffers from the dependency problem.
3815 * Specifically, a propagable quadratic expression is a sum expression such that there is at least one expr that appears
3816 * at least twice (because of simplification, this means it appears in a quadratic terms and somewhere else).
3817 * For example: \f$x^2 + y^2\f$ is not a propagable quadratic expression; \f$x^2 + x\f$ is a propagable quadratic expression;
3820 * Furthermore, we distinguish between propagable and non-propagable terms. A term is propagable if any of the expressions
3821 * involved in it appear somewhere else. For example, \f$xy + z^2 + z\f$ is a propagable quadratic, the term \f$xy\f$ is
3822 * non-propagable, and \f$z^2\f$ is propagable. For propagation, non-propagable terms are handled as if they were linear
3823 * terms, that is, we do not use the activity of \f$x\f$ and \f$y\f$ to compute the activity of \f$xy\f$ but rather we use directly
3824 * the activity of \f$xy\f$. Similarly, we do not backward propagate to \f$x\f$ and \f$y\f$ (the product expr handler will do this),
3825 * but we backward propagate to \f$x*y\f$. More technically, we register \f$xy\f$ for its activity usage, rather than\f$x\f$ and \f$y\f$.
3827 * For propagation, we store the quadratic in our data structure in the following way: We count how often a variable
3828 * appears. Then, a bilinear product expr_i * expr_j is stored as expr_i * expr_j if # expr_i appears > # expr_j
3829 * appears. When # expr_i appears = # expr_j appears, it then it will be stored as expr_i * expr_j if and only if
3830 * expr_i < expr_j, where '<' is the expression order (see \ref EXPR_ORDER "Ordering Rules" in \ref scip_expr.h).
3831 * Heuristically, this should be useful for propagation. The intuition is that by factoring out the variable that
3834 * Simple convex quadratics like \f$x^2 + y^2\f$ are ignored since the default nlhdlr will take care of them.
3837 * @note Common subexpressions are also assumed to have been identified, the hashing will fail otherwise!
3841 * - expr < expr * other_expr: u*v < w holds if and only if v < w (OR8), but here w = u < v, since expr comes before
3851 * It also implies that x^-2 < x^-1, but since, so far, we do not interpret x^-2 as (x^-1)^2, it is not a problem.
3907 SCIPdebugMsg(scip, "expr %p is not propagable and in presolving -> abort detect\n", (void*)expr);
3911 /* if we do not use intersection cuts and are not propagable, then we do not want to handle it at all;
3912 * if not propagable, then we need to check the curvature to decide if we want to generate intersection cuts
3943 SCIPexprGetQuadraticData(expr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, &nbilin, NULL, NULL);
3959 SCIPinfoMessage(scip, NULL, "quadterm %d propagable, using %p, unbounded=%d\n", i, (void*)argexpr, nbilin >
3966 * we should make use nlhdlrs in pow or product for this term, so we register usage of the square or product
3980 SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable square, using %p\n", i, (void*)sqrexpr);
3985 /* we have expr1 * other_expr or other_expr * expr1; know that expr1 is non propagable, but to decide if
3986 * we want the bounds of expr1 or of the product expr1 * other_expr (or other_expr * expr1), we have to
3988 * frequency), we can deduce that other_expr doesn't appear anywhere else (i.e. is non propagable) if the
3989 * product is of the form expr1 * other_expr; however, if we see other_expr * expr1 we need to find
4003 SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable product, using %p\n", i, (void*)prodexpr);
4009 /* check if other_expr is propagable in which case we need the bounds of expr1; otherwise we just need
4012 * TODO this should be done faster, maybe store pos1 in bilinexprterm or store quadexprterm's in bilinexprterm
4024 SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable alien product, using %p\n", i, (void*)argexpr);
4054 assert(SCIPgetStage(scip) >= SCIP_STAGE_INITSOLVE); /* separation should only be required in (init)solving stage */
4056 /* check if we can do something more: check curvature of quadratic function stored in nlexprdata
4057 * this is currently only used to decide whether we want to separate, so it can be skipped if in presolve
4060 SCIP_CALL( SCIPcomputeExprQuadraticCurvature(scip, expr, &nlexprdata->curvature, NULL, nlhdlrdata->useintersectioncuts) );
4066 if( nlhdlrdata->useintersectioncuts && eigenvalues != NULL && (*enforcing & SCIP_NLHDLR_METHOD_SEPABELOW) ==
4072 if( nlhdlrdata->useintersectioncuts && eigenvalues != NULL && (*enforcing & SCIP_NLHDLR_METHOD_SEPAABOVE) == FALSE &&
4093 SCIPexprGetQuadraticData(expr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
4106 SCIPdebugMsg(scip, "expr %p is quadratic and propagable -> propagate and separate\n", (void*)expr);
4118 SCIPdebugMsg(scip, "expr is %s in the original variables\n", nlexprdata->curvature == SCIP_EXPRCURV_CONCAVE ? "concave" : "convex");
4155 SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, &nbilinexprs, NULL, NULL);
4183 *auxvalue += coef * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr1)) * SCIPgetSolVal(scip, sol,
4232 /* right now can use interesction cuts only if a basic LP solution is at hand; TODO: in principle we can do something
4235 if( sol != NULL || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL || !SCIPisLPSolBasic(scip) )
4244 if( (nlhdlrdata->atwhichnodes == -1 && depth != 0) || (nlhdlrdata->atwhichnodes != -1 && depth % nlhdlrdata->atwhichnodes != 0) )
4250 /* do not add more than ncutslimitroot cuts in root node and ncutslimit cuts in the non-root nodes */
4257 /*else if( (depth > 0 && nlhdlrdata->ncutsadded - nlhdlrdata->lastncuts >= nlhdlrdata->ncutslimit) || (depth == 0 &&
4293 /* we can't build an intersection cut when the expr is the root of some constraint and also a subexpression of
4301 /* if we are the root of a constraint and we are feasible w.r.t our auxiliary variables, that is, auxvalue is
4304 if( cons == nlhdlrexprdata->cons && ((overestimate && (SCIPgetLhsNonlinear(cons)) - auxvalue < SCIPfeastol(scip)) ||
4333 SCIP_CALL( generateIntercut(scip, expr, nlhdlrdata, nlhdlrexprdata, cons, sol, rowprep, overestimate, &success) );
4350 SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, sol, nlhdlrdata->mincutviolation, &violation, &success) );
4354 /* if cut looks good (numerics ok and cutting off solution), then turn into row and add to sepastore */
4364 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%s_intersection_quadratic%p_lp%" SCIP_LONGINT_FORMAT,
4372 printf(" -> found maxquad-free cut <%s>: act=%f, lhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n\n",
4379 /*printf("SCIP DEPTH %d got a cut with violation %g, efficacy %g and r/e %g\n", SCIPgetSubscipDepth(scip),
4380 * violation, SCIPgetCutEfficacy(scip, NULL, row), SCIPgetRowMaxCoef(scip, row) / SCIPgetRowMinCoef(scip, row) /
4384 if( ! nlhdlrdata->ignorehighre || SCIPgetRowMaxCoef(scip, row) / SCIPgetRowMinCoef(scip, row) / SCIPgetCutEfficacy(scip, NULL, row) < 1e9 )
4403 nlhdlrdata->densitysum += (SCIP_Real) SCIProwprepGetNVars(rowprep) / (SCIP_Real) SCIPgetNVars(scip);
4426 * Interval arithmetic suffices when no variable appears twice, however this is seldom the case, so we try
4430 * where \f$q_l = a_l \text{expr}_l^2 + c_l \text{expr}_l + \sum_{i \in P_l} b_{il} \text{expr}_i \text{expr}_l\f$
4431 * 2. build interval quadratic functions, i.e., \f$a x^2 + b x\f$ where \f$b\f$ is an interval, i.e.,
4434 * \f$\min/\max a_l \text{expr}_l^2 + \text{expr}_l [\sum_{i \in P_l} b_{il} \text{expr}_i + c_l] : \text{expr}_l \in [\text{expr}_l]\f$
4437 * 1. The \f$l\f$-th quadratic expr (expressions that appear quadratically) is associated with \f$q_l\f$.
4438 * 2. `nlhdlrdata->quadactivities[l]` is the activity of \f$q_l\f$ as computed in the description above.
4439 * 3. The \f$q_l\f$ of a quadratic term might be empty, in which case `nlhdlrdata->quadactivities[l]` is [0,0].\n
4440 * For example, consider \f$x^2 + xy\f$. There are two quadratic expressions, \f$x\f$ and \f$y\f$.
4441 * The \f$q\f$ associated to \f$x\f$ is \f$x^2 + xy\f$, while the \f$q\f$ associated to \f$y\f$ is empty.
4445 * @note The order matters! If \f$\text{expr}_i\, \text{expr}_l\f$ is a term in the quadratic, then \f$i\f$ is *not* in \f$P_l\f$
4465 SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, NULL, NULL);
4488 SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &nlhdlrexprdata->linactivity, nlhdlrexprdata->linactivity, linterminterval);
4520 SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, &sqrexpr);
4524 /* term is not propagable, i.e., the exprs involved in term only appear once; thus use the activity of the
4525 * quadratic term directly and not the activity of the exprs involed in the term. See also documentation of
4549 SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>, where <expr> is: ", sqrcoef);
4564 /* the quadratic expression expr1 appears only as expr1 * expr2, so its 'q' is expr1 * expr2 */
4577 SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>, where <expr> is: ", prodcoef);
4583 /* the quadratic expression expr1 appears as expr2 * expr1, thus its 'q' is empty, see also the Notes
4596 SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, NULL);
4635 SCIPinfoMessage(scip, NULL, " [%g,%g]\n", SCIPexprGetActivity(expr2).inf, SCIPexprGetActivity(expr2).sup);
4639 /* TODO: under which assumptions do we know that we just need to compute min or max? its probably the locks that give some information here */
4653 SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>^2 + [%g,%g] <expr>, where <expr> is: ", sqrcoef, b.inf, b.sup);
4662 SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &nlhdlrexprdata->quadactivity, nlhdlrexprdata->quadactivity, nlhdlrexprdata->quadactivities[i]);
4693 SCIPdebugMsg(scip, "Activity of quadratic part is [%g, %g]\n", nlhdlrexprdata->quadactivity.inf, nlhdlrexprdata->quadactivity.sup);
4697 SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, nlhdlrexprdata->linactivity, nlhdlrexprdata->quadactivity);
4699 nlhdlrexprdata->activitiestag = SCIPgetCurBoundsTagNonlinear(SCIPfindConshdlr(scip, "nonlinear"));
4706 * @note the implemented technique is a proxy for solving the problem min/max{ x_i : quad expr in [quad expr] }
4724 SCIPdebugMsg(scip, "Reverse propagation of quadratic expr given bounds = [%g,%g]\n", bounds.inf, bounds.sup);
4744 * if the activity stored in expr is more recent than the partial activities stored in this nlhdlrexprdata,
4749 SCIP_CALL( nlhdlrIntevalQuadratic(scip, nlhdlr, expr, nlhdlrexprdata, &quadactivity, NULL, NULL) );
4752 SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, NULL, NULL);
4754 /* propagate linear part in rhs = expr's interval - quadratic activity; first, reconstruct the quadratic activity */
4756 nlhdlrexprdata->nneginfinityquadact > 0 ? -SCIP_INTERVAL_INFINITY : nlhdlrexprdata->minquadfiniteact,
4757 nlhdlrexprdata->nposinfinityquadact > 0 ? SCIP_INTERVAL_INFINITY : nlhdlrexprdata->maxquadfiniteact);
4761 SCIP_CALL( reversePropagateLinearExpr(scip, linexprs, nlinexprs, lincoefs, constant, rhs, infeasible, nreductions) );
4767 /* propagate quadratic part in expr's interval - linear activity, where linear activity was computed in INTEVAL.
4771 * - for each expression expr_i, write the quadratic expression as a_i expr^2_i + expr_i ( \sum_{j \in J_i} b_ij
4773 * - compute the interval b = [\sum_{j \in J_i} b_ij expr_j + c_i], where J_i are all the indices j such that the
4775 * - use some technique (like the one in nlhdlrIntevalQuadratic), to evaluate the activity of rest_i = [quadratic
4779 * However, this might be expensive, especially computing rest_i. Hence, we implement a simpler version.
4780 * - we use the same partition as in nlhdlrIntevalQuadratic for the bilinear terms. This way, b = [\sum_{j \in P_i}
4781 * b_ij expr_j + c_i], where P_i is the set of indices j such that expr_i * expr_j appears in that order
4782 * - we evaluate the activity of rest_i as sum_{k \neq i} [\min q_k, \max q_k] where q_k = a_k expr_k^2 + [\sum_{j
4783 * \in P_k} b_jk expr_j + c_k] expr_k. The intervals [\min q_k, \max q_k] were already computed in
4786 * A downside of the above is that we might not deduce any bounds for variables that appear less often. For example,
4787 * consider x^2 + x * y + x * z + y * z + z. This quadratic gets partitioned as (x^2 + x*y + x*z) + (z*y + z). The
4788 * first parenthesis is interpreted as a function of x, while the second one as a function of z.
4789 * To also get bounds on y, after reverse propagating x in x^2 + x*y + x*z \in rhs, we rewrite this as y + z \in rhs/x -
4792 * \sum_{j \in J_i} b_ij expr_j in ([expr activity] - quadratic expression in expr_k for k \neq i - c_i) / expr_i - a_i expr_i,
4793 * compute an interval for the right hand side (see computeRangeForBilinearProp) and use that to propagate the
4797 * The idea of that technique was to borrow a bilinear term expr_k expr_l when propagating expr_l and the quadratic
4799 * Since in P_l we only consider the indices of expressions that appear multiplying expr_l as _second_ factor, we
4801 * The problem is that the contribution of b_kl * expr_k * expr_l to rest_i is not just [b_kl * expr_k * expr_l], but
4802 * rather quadactivities[k] (= max/min of a_k expr_k^2 + expr_k * [c_k + sum_i \in P_k b_ki expr_i]).
4804 * But, if expr_k only appears as expr_k * expr_l, then quadactivities[k] = [b_kl * expr_k * expr_l]. So this
4826 SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, &sqrexpr);
4831 * if [q_i].sup = +infinity and there is = 1 contributing +infinity -> rest_i.sup = maxquadfiniteact
4834 * if [q_i].sup = finite and there is = 0 contributing +infinity -> rest_i.sup = maxquadfiniteact - [q_i].sup
4846 rest_i.sup = nlhdlrexprdata->maxquadfiniteact - SCIPintervalGetSup(nlhdlrexprdata->quadactivities[i]);
4864 rest_i.inf = nlhdlrexprdata->minquadfiniteact - SCIPintervalGetInf(nlhdlrexprdata->quadactivities[i]);
4874#ifdef SCIP_DISABLED_CODE /* I (SV) added the following in cons_quadratic to fix/workaround some bug. Maybe we'll need this here, too? */
4876 * what we tried to do here is to remove the contribution of the i'th bilinear term (=bilinterm) to [minquadactivity,maxquadactivity] from rhs
4877 * however, quadactivity is computed differently (as x*(a1*y1+...+an*yn)) than q_i (a*ak*yk) and since interval arithmetics do overestimation,
4878 * it can happen that q_i is actually slightly larger than quadactivity, which results in rest_i being (slightly) empty
4879 * a proper fix could be to compute the quadactivity also as x*a1*y1+...+x*an*yn if sqrcoef=0, but due to taking
4912 /* qexpr only appears in a term of the form qexpr * other_expr (or other_expr * qexpr); we only care about
4913 * getting bounds for the product, thus we will compute these bounds when qexpr appears as qexpr *
4914 * other_expr; note that if it appears as other_expr * qexpr, then when we process other_expr bounds for the
4916 * TODO: we can actually avoid computing rhs_i in the case that qexpr is not propagable and it appears as
4982 /* if 0 is not in [expr_i], then propagate bilincoefs^T bilinexpr in rhs_i/expr_i - a_i expr_i - c_i */
5004 /* TODO FIXME: we are overestimating the number of reductions: an expr might be tightened many times! */
5060 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY,
5111 "determines at which nodes cut is used (if it's -1, it's used only at the root node, if it's n >= 0, it's used at every multiple of n",
constraint handler for nonlinear constraints specified by algebraic expressions
power and signed power expression handlers
product expression handler
sum expression handler
variable expression handler