Scippy

SCIP

Solving Constraint Integer Programs

cons_quadratic.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2020 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_quadratic.c
17  * @ingroup DEFPLUGINS_CONS
18  * @brief constraint handler for quadratic constraints \f$\textrm{lhs} \leq \sum_{i,j=1}^n a_{i,j} x_i x_j + \sum_{i=1}^n b_i x_i \leq \textrm{rhs}\f$
19  * @author Stefan Vigerske
20  * @author Benjamin Mueller
21  * @author Felipe Serrano
22  *
23  * @todo SCIP might fix linear variables on +/- infinity; remove them in presolve and take care later
24  * @todo round constraint sides to integers if all coefficients and variables are (impl.) integer
25  * @todo constraints in one variable should be replaced by linear or bounddisjunction constraint
26  * @todo check if some quadratic terms appear in several constraints and try to simplify (e.g., nous1)
27  * @todo skip separation in enfolp if for current LP (check LP id) was already separated
28  * @todo watch unbounded variables to enable/disable propagation
29  * @todo sort order in bilinvar1/bilinvar2 such that the var which is involved in more terms is in bilinvar1, and use this info propagate and AddLinearReform
30  * @todo underestimate for multivariate concave quadratic terms as in cons_nonlinear
31  */
32 
33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34 
35 #define SCIP_PRIVATE_ROWPREP
36 
37 #include "blockmemshell/memory.h"
38 #include <ctype.h>
39 #include "nlpi/nlpi.h"
40 #include "nlpi/nlpi_ipopt.h"
41 #include "nlpi/pub_expr.h"
42 #include "nlpi/type_expr.h"
43 #include "scip/cons_and.h"
44 #include "scip/cons_bounddisjunction.h"
45 #include "scip/cons_linear.h"
46 #include "scip/cons_nonlinear.h"
47 #include "scip/cons_quadratic.h"
48 #include "scip/cons_varbound.h"
49 #include "scip/debug.h"
50 #include "scip/heur_subnlp.h"
51 #include "scip/heur_trysol.h"
52 #include "scip/intervalarith.h"
53 #include "scip/pub_cons.h"
54 #include "scip/pub_event.h"
55 #include "scip/pub_heur.h"
56 #include "scip/pub_lp.h"
57 #include "scip/pub_message.h"
58 #include "scip/pub_misc.h"
59 #include "scip/pub_misc_sort.h"
60 #include "scip/pub_nlp.h"
61 #include "scip/pub_sol.h"
62 #include "scip/pub_tree.h"
63 #include "scip/pub_var.h"
64 #include "scip/scip_branch.h"
65 #include "scip/scip_cons.h"
66 #include "scip/scip_copy.h"
67 #include "scip/scip_cut.h"
68 #include "scip/scip_event.h"
69 #include "scip/scip_general.h"
70 #include "scip/scip_heur.h"
71 #include "scip/scip_lp.h"
72 #include "scip/scip_mem.h"
73 #include "scip/scip_message.h"
74 #include "scip/scip_nlp.h"
75 #include "scip/scip_nonlinear.h"
76 #include "scip/scip_numerics.h"
77 #include "scip/scip_param.h"
78 #include "scip/scip_prob.h"
79 #include "scip/scip_probing.h"
80 #include "scip/scip_sepa.h"
81 #include "scip/scip_sol.h"
82 #include "scip/scip_solve.h"
83 #include "scip/scip_solvingstats.h"
84 #include "scip/scip_tree.h"
85 #include "scip/scip_var.h"
86 #include <string.h>
87 
88 /* Inform compiler that this code accesses the floating-point environment, so that
89  * certain optimizations should be omitted (http://www.cplusplus.com/reference/cfenv/FENV_ACCESS/).
90  * Not supported by Clang (gives warning) and GCC (silently), at the moment.
91  */
92 #if defined(__INTEL_COMPILER) || defined(_MSC_VER)
93 #pragma fenv_access (on)
94 #elif defined __GNUC__
95 #pragma STDC FENV_ACCESS ON
96 #endif
97 
98 /* constraint handler properties */
99 #define CONSHDLR_NAME "quadratic"
100 #define CONSHDLR_DESC "quadratic constraints of the form lhs <= b' x + x' A x <= rhs"
101 #define CONSHDLR_SEPAPRIORITY 10 /**< priority of the constraint handler for separation */
102 #define CONSHDLR_ENFOPRIORITY -50 /**< priority of the constraint handler for constraint enforcing */
103 #define CONSHDLR_CHECKPRIORITY -4000000 /**< priority of the constraint handler for checking feasibility */
104 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
105 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
106 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
107  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
108 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
109 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
110 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
111 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
112 
113 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP /**< propagation timing mask of the constraint handler */
114 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_ALWAYS /**< presolving timing of the constraint handler (fast, medium, or exhaustive) */
115 
116 #define MAXDNOM 10000LL /**< maximal denominator for simple rational fixed values */
117 #define NONLINCONSUPGD_PRIORITY 40000 /**< priority of upgrading nonlinear constraints */
118 #define INITLPMAXVARVAL 1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */
119 
120 /* Activating this define enables reformulation of bilinear terms x*y with implications from x to y into linear terms.
121  * However, implications are not enforced by SCIP. Thus, if, e.g., the used implication was derived from this constraint and we then reformulate the constraint,
122  * then the implication may not be enforced in a solution.
123  * This issue need to be fixed before this feature can be enabled.
124  */
125 /* #define CHECKIMPLINBILINEAR */
126 
127 /* enable new propagation for bivariate quadratic terms */
128 #define PROPBILINNEW
129 
130 /* epsilon for differentiating between a boundary and interior point */
131 #define INTERIOR_EPS 1e-1
132 
133 /* scaling factor for gauge function */
134 #define GAUGESCALE 0.99999
135 
136 #define ROWPREP_SCALEUP_VIOLNONZERO (10.0*SCIPepsilon(scip)) /**< minimal violation for considering up-scaling of rowprep (we want to avoid upscaling very small violations) */
137 #define ROWPREP_SCALEUP_MINVIOLFACTOR 2.0 /**< scale up will target a violation of ~MINVIOLFACTOR*minviol, where minviol is given by caller */
138 #define ROWPREP_SCALEUP_MAXMINCOEF (1.0 / SCIPfeastol(scip)) /**< scale up only if min. coef is below this number (before scaling) */
139 #define ROWPREP_SCALEUP_MAXMAXCOEF SCIPgetHugeValue(scip) /**< scale up only if max. coef will not exceed this number by scaling */
140 #define ROWPREP_SCALEUP_MAXSIDE SCIPgetHugeValue(scip) /**< scale up only if side will not exceed this number by scaling */
141 #define ROWPREP_SCALEDOWN_MINMAXCOEF (1.0 / SCIPfeastol(scip)) /**< scale down if max. coef is at least this number (before scaling) */
142 #define ROWPREP_SCALEDOWN_MINCOEF SCIPfeastol(scip) /**< scale down only if min. coef does not drop below this number by scaling */
143 
144 /*
145  * Data structures
146  */
147 
148 /** eventdata for variable bound change events in quadratic constraints */
149 struct SCIP_QuadVarEventData
150 {
151  SCIP_CONS* cons; /**< the constraint */
152  int varidx; /**< the index of the variable which bound change is caught, positive for linear variables, negative for quadratic variables */
153  int filterpos; /**< position of eventdata in SCIP's event filter */
154 };
155 
156 /** Data of a quadratic constraint. */
157 struct SCIP_ConsData
158 {
159  SCIP_Real lhs; /**< left hand side of constraint */
160  SCIP_Real rhs; /**< right hand side of constraint */
161 
162  int nlinvars; /**< number of linear variables */
163  int linvarssize; /**< length of linear variable arrays */
164  SCIP_VAR** linvars; /**< linear variables */
165  SCIP_Real* lincoefs; /**< coefficients of linear variables */
166  SCIP_QUADVAREVENTDATA** lineventdata; /**< eventdata for bound change of linear variable */
167 
168  int nquadvars; /**< number of variables in quadratic terms */
169  int quadvarssize; /**< length of quadratic variable terms arrays */
170  SCIP_QUADVARTERM* quadvarterms; /**< array with quadratic variable terms */
171 
172  int nbilinterms; /**< number of bilinear terms */
173  int bilintermssize; /**< length of bilinear term arrays */
174  SCIP_BILINTERM* bilinterms; /**< bilinear terms array */
175  int* bilintermsidx; /**< unique index of each bilinear term xy in the bilinestimators array of the constraint handler data */
176 
177  SCIP_NLROW* nlrow; /**< a nonlinear row representation of this constraint */
178 
179  unsigned int linvarssorted:1; /**< are the linear variables already sorted? */
180  unsigned int linvarsmerged:1; /**< are equal linear variables already merged? */
181  unsigned int quadvarssorted:1; /**< are the quadratic variables already sorted? */
182  unsigned int quadvarsmerged:1; /**< are equal quadratic variables already merged? */
183  unsigned int bilinsorted:1; /**< are the bilinear terms already sorted? */
184  unsigned int bilinmerged:1; /**< are equal bilinear terms (and bilinear terms with zero coefficient) already merged? */
185 
186  unsigned int isconvex:1; /**< is quadratic function is convex ? */
187  unsigned int isconcave:1; /**< is quadratic function is concave ? */
188  unsigned int iscurvchecked:1; /**< is quadratic function checked on convexity or concavity ? */
189  unsigned int isremovedfixings:1; /**< did we removed fixed/aggr/multiaggr variables ? */
190  unsigned int ispropagated:1; /**< was the constraint propagated with respect to the current bounds ? */
191  unsigned int ispresolved:1; /**< did we checked for possibilities of upgrading or implicit integer variables ? */
192  unsigned int initialmerge:1; /**< did we perform an initial merge and clean in presolving yet ? */
193 #ifdef CHECKIMPLINBILINEAR
194  unsigned int isimpladded:1; /**< has there been an implication added for a binary variable in a bilinear term? */
195 #endif
196  unsigned int isgaugeavailable:1; /**< is the gauge function computed? */
197  unsigned int isedavailable:1; /**< is the eigen decomposition of A available? */
198 
199  SCIP_Real minlinactivity; /**< sum of minimal activities of all linear terms with finite minimal activity */
200  SCIP_Real maxlinactivity; /**< sum of maximal activities of all linear terms with finite maximal activity */
201  int minlinactivityinf; /**< number of linear terms with infinite minimal activity */
202  int maxlinactivityinf; /**< number of linear terms with infinity maximal activity */
203  SCIP_INTERVAL quadactivitybounds; /**< bounds on the activity of the quadratic term, if up to date, otherwise empty interval */
204  SCIP_Real activity; /**< activity of quadratic function w.r.t. current solution */
205  SCIP_Real lhsviol; /**< violation of lower bound by current solution (used temporarily inside constraint handler) */
206  SCIP_Real rhsviol; /**< violation of lower bound by current solution (used temporarily inside constraint handler) */
207 
208  int linvar_maydecrease; /**< index of a variable in linvars that may be decreased without making any other constraint infeasible, or -1 if none */
209  int linvar_mayincrease; /**< index of a variable in linvars that may be increased without making any other constraint infeasible, or -1 if none */
210 
211  SCIP_VAR** sepaquadvars; /**< variables corresponding to quadvarterms to use in separation, only available in solving stage */
212  int* sepabilinvar2pos; /**< position of second variable in bilinear terms to use in separation, only available in solving stage */
213  SCIP_Real lincoefsmin; /**< minimal absolute value of coefficients in linear part, only available in solving stage */
214  SCIP_Real lincoefsmax; /**< maximal absolute value of coefficients in linear part, only available in solving stage */
215 
216  SCIP_Real* factorleft; /**< coefficients of left factor if constraint function is factorable */
217  SCIP_Real* factorright; /**< coefficients of right factor if constraint function is factorable */
218 
219  SCIP_Real* gaugecoefs; /**< coefficients of the gauge function */
220  SCIP_Real gaugeconst; /**< constant of the gauge function */
221  SCIP_Real* interiorpoint; /**< interior point of the region defined by the convex function */
222  SCIP_Real interiorpointval; /**< function value at interior point */
223 
224  SCIP_Real* eigenvalues; /**< eigenvalues of A */
225  SCIP_Real* eigenvectors; /**< orthonormal eigenvectors of A; if A = P D P^T, then eigenvectors is P^T */
226  SCIP_Real* bp; /**< stores b * P where b are the linear coefficients of the quadratic vars */
227  SCIP_Real maxnonconvexity; /**< nonconvexity measure: estimate on largest absolute value of nonconvex eigenvalues */
228 
229  SCIP_Bool isdisaggregated; /**< has the constraint already been disaggregated? if might happen that more disaggreation would be potentially
230  possible, but we reached the maximum number of sparsity components during presolveDisaggregate() */
231 };
232 
233 /** quadratic constraint update method */
234 struct SCIP_QuadConsUpgrade
235 {
236  SCIP_DECL_QUADCONSUPGD((*quadconsupgd)); /**< method to call for upgrading quadratic constraint */
237  int priority; /**< priority of upgrading method */
238  SCIP_Bool active; /**< is upgrading enabled */
239 };
240 typedef struct SCIP_QuadConsUpgrade SCIP_QUADCONSUPGRADE; /**< quadratic constraint update method */
241 
242 /** structure to store everything needed for using linear inequalities to improve upon the McCormick relaxation */
243 struct BilinearEstimator
244 {
245  SCIP_VAR* x; /**< first variable */
246  SCIP_VAR* y; /**< second variable */
247  SCIP_Real inequnderest[6]; /**< at most two inequalities that can be used to underestimate xy; stored as (xcoef,ycoef,constant) with xcoef x <= ycoef y + constant */
248  SCIP_Real ineqoverest[6]; /**< at most two inequalities that can be used to overestimate xy; stored as (xcoef,ycoef,constant) with xcoef x <= ycoef y + constant */
249  SCIP_Real maxnonconvexity; /**< estimate on largest absolute value of nonconvex eigenvalues of all quadratic constraint containing xy */
250  int ninequnderest; /**< total number of inequalities for underestimating xy */
251  int nineqoverest; /**< total number of inequalities for overestimating xy */
252  int nunderest; /**< number of constraints that require to underestimate xy */
253  int noverest; /**< number of constraints that require to overestimate xy */
254 
255  SCIP_Real lastimprfac; /**< last achieved improvement factor */
256 };
257 typedef struct BilinearEstimator BILINESTIMATOR;
258 
259 /** constraint handler data */
260 struct SCIP_ConshdlrData
261 {
262  int replacebinaryprodlength; /**< length of linear term which when multiplied with a binary variable is replaced by an auxiliary variable and an equivalent linear formulation */
263  int empathy4and; /**< how much empathy we have for using the AND constraint handler: 0 avoid always; 1 use sometimes; 2 use as often as possible */
264  SCIP_Bool binreforminitial; /**< whether to make constraints added due to replacing products with binary variables initial */
265  SCIP_Bool binreformbinaryonly;/**< whether to consider only binary variables when reformulating products with binary variables */
266  SCIP_Real binreformmaxcoef; /**< factor on 1/feastol to limit coefficients and coef range in linear constraints created by binary reformulation */
267  SCIP_Real cutmaxrange; /**< maximal range (maximal coef / minimal coef) of a cut in order to be added to LP */
268  SCIP_Bool linearizeheursol; /**< whether linearizations of convex quadratic constraints should be added to cutpool when some heuristics finds a new solution */
269  SCIP_Bool checkcurvature; /**< whether functions should be checked for convexity/concavity */
270  SCIP_Bool checkfactorable; /**< whether functions should be checked to be factorable */
271  char checkquadvarlocks; /**< whether quadratic variables contained in a single constraint should be forced to be at their lower or upper bounds ('d'isable, change 't'ype, add 'b'ound disjunction) */
272  SCIP_Bool linfeasshift; /**< whether to make solutions in check feasible if possible */
273  int maxdisaggrsize; /**< maximum number of components when disaggregating a quadratic constraint (<= 1: off) */
274  char disaggrmergemethod; /**< method on merging blocks in disaggregation */
275  int maxproprounds; /**< limit on number of propagation rounds for a single constraint within one round of SCIP propagation during solve */
276  int maxproproundspresolve; /**< limit on number of propagation rounds for a single constraint within one presolving round */
277  SCIP_Real sepanlpmincont; /**< minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation */
278  SCIP_Bool enfocutsremovable; /**< are cuts added during enforcement removable from the LP in the same node? */
279  SCIP_Bool gaugecuts; /**< should convex quadratics generated strong cuts via gauge function? */
280  SCIP_Bool projectedcuts; /**< should convex quadratics generated strong cuts via projections? */
281  char interiorcomputation;/**< how the interior point should be computed: 'a'ny point per constraint, 'm'ost interior per constraint */
282  char branchscoring; /**< method to use to compute score of branching candidates */
283  int enfolplimit; /**< maximum number of enforcement round before declaring the LP relaxation
284  * infeasible (-1: no limit); WARNING: if this parameter is not set to -1,
285  * SCIP might declare sub-optimal solutions optimal or feasible instances
286  * infeasible; thus, the result returned by SCIP might be incorrect!
287  */
288  SCIP_HEUR* subnlpheur; /**< a pointer to the subnlp heuristic, if available */
289  SCIP_HEUR* trysolheur; /**< a pointer to the trysol heuristic, if available */
290  SCIP_EVENTHDLR* eventhdlr; /**< our handler for variable bound change events */
291  int newsoleventfilterpos; /**< filter position of new solution event handler, if caught */
292  SCIP_Bool sepanlp; /**< where linearization of the NLP relaxation solution added? */
293  SCIP_NODE* lastenfonode; /**< the node for which enforcement was called the last time (and some constraint was violated) */
294  int nenforounds; /**< counter on number of enforcement rounds for the current node */
295  SCIP_QUADCONSUPGRADE** quadconsupgrades; /**< quadratic constraint upgrade methods for specializing quadratic constraints */
296  int quadconsupgradessize; /**< size of quadconsupgrade array */
297  int nquadconsupgrades; /**< number of quadratic constraint upgrade methods */
298 
299  BILINESTIMATOR* bilinestimators; /**< array containing all required information for using stronger estimators for each bilinear term in all quadratic constraints */
300  int nbilinterms; /**< number of bilinear terms in all quadratic constraints */
301 
302  SCIP_Bool usebilinineqbranch; /**< should linear inequalities be considered when computing the branching scores for bilinear terms? */
303  SCIP_Bool storedbilinearterms; /**< did we already try to store all bilinear terms? */
304 
305  SCIP_Real minscorebilinterms; /**< minimal required score in order to use linear inequalities for tighter bilinear relaxations */
306  SCIP_Real mincurvcollectbilinterms;/**< minimal curvature of constraints to be considered when returning bilinear terms to other plugins */
307  int bilinineqmaxseparounds; /**< maximum number of separation rounds to use linear inequalities for the bilinear term relaxation in a local node */
308 };
309 
310 
311 /*
312  * local methods for managing quadratic constraint update methods
313  */
314 
315 
316 /** checks whether a quadratic constraint upgrade method has already be registered */
317 static
318 SCIP_Bool conshdlrdataHasUpgrade(
319  SCIP* scip, /**< SCIP data structure */
320  SCIP_CONSHDLRDATA* conshdlrdata, /**< constraint handler data */
321  SCIP_DECL_QUADCONSUPGD((*quadconsupgd)), /**< method to call for upgrading quadratic constraint */
322  const char* conshdlrname /**< name of the constraint handler */
323  )
324 {
325  int i;
326 
327  assert(scip != NULL);
328  assert(conshdlrdata != NULL);
329  assert(quadconsupgd != NULL);
330  assert(conshdlrname != NULL);
331 
332  for( i = conshdlrdata->nquadconsupgrades - 1; i >= 0; --i )
333  {
334  if( conshdlrdata->quadconsupgrades[i]->quadconsupgd == quadconsupgd )
335  {
336  SCIPwarningMessage(scip, "Try to add already known upgrade message for constraint handler <%s>.\n", conshdlrname);
337  return TRUE;
338  }
339  }
340 
341  return FALSE;
342 }
343 
344 /*
345  * Local methods
346  */
347 
348 /** translate from one value of infinity to another
349  *
350  * if val is >= infty1, then give infty2, else give val
351  */
352 #define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))
353 
354 /** catches variable bound change events on a linear variable in a quadratic constraint */
355 static
356 SCIP_RETCODE catchLinearVarEvents(
357  SCIP* scip, /**< SCIP data structure */
358  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
359  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
360  int linvarpos /**< position of variable in linear variables array */
361  )
362 {
363  SCIP_CONSDATA* consdata;
364  SCIP_QUADVAREVENTDATA* eventdata;
365  SCIP_EVENTTYPE eventtype;
366 
367  assert(scip != NULL);
368  assert(eventhdlr != NULL);
369  assert(cons != NULL);
370 
371  consdata = SCIPconsGetData(cons);
372  assert(consdata != NULL);
373 
374  assert(linvarpos >= 0);
375  assert(linvarpos < consdata->nlinvars);
376  assert(consdata->lineventdata != NULL);
377 
378  SCIP_CALL( SCIPallocBlockMemory(scip, &eventdata) );
379 
380  eventdata->cons = cons;
381  eventdata->varidx = linvarpos;
382 
383  eventtype = SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
384  if( !SCIPisInfinity(scip, consdata->rhs) )
385  {
386  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest
387  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
388  if( consdata->lincoefs[linvarpos] > 0.0 )
389  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
390  else
391  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
392  }
393  if( !SCIPisInfinity(scip, -consdata->lhs) )
394  {
395  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest
396  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
397  if( consdata->lincoefs[linvarpos] > 0.0 )
398  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
399  else
400  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
401  }
402 
403  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->linvars[linvarpos], eventtype, eventhdlr, (SCIP_EVENTDATA*)eventdata, &eventdata->filterpos) );
404 
405  consdata->lineventdata[linvarpos] = eventdata;
406 
407  /* invalidate activity information
408  * NOTE: It could happen that a constraint gets temporary deactivated and some variable bounds change. In this case
409  * we do not recognize those bound changes with the variable events and thus we have to recompute the activities.
410  */
411  consdata->minlinactivity = SCIP_INVALID;
412  consdata->maxlinactivity = SCIP_INVALID;
413  consdata->minlinactivityinf = -1;
414  consdata->maxlinactivityinf = -1;
415 
416  return SCIP_OKAY;
417 }
418 
419 /** drops variable bound change events on a linear variable in a quadratic constraint */
420 static
421 SCIP_RETCODE dropLinearVarEvents(
422  SCIP* scip, /**< SCIP data structure */
423  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
424  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
425  int linvarpos /**< position of variable in linear variables array */
426  )
427 {
428  SCIP_CONSDATA* consdata;
429  SCIP_EVENTTYPE eventtype;
430 
431  assert(scip != NULL);
432  assert(eventhdlr != NULL);
433  assert(cons != NULL);
434 
435  consdata = SCIPconsGetData(cons);
436  assert(consdata != NULL);
437 
438  assert(linvarpos >= 0);
439  assert(linvarpos < consdata->nlinvars);
440  assert(consdata->lineventdata != NULL);
441  assert(consdata->lineventdata[linvarpos] != NULL);
442  assert(consdata->lineventdata[linvarpos]->cons == cons);
443  assert(consdata->lineventdata[linvarpos]->varidx == linvarpos);
444  assert(consdata->lineventdata[linvarpos]->filterpos >= 0);
445 
446  eventtype = SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
447  if( !SCIPisInfinity(scip, consdata->rhs) )
448  {
449  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest
450  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
451  if( consdata->lincoefs[linvarpos] > 0.0 )
452  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
453  else
454  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
455  }
456  if( !SCIPisInfinity(scip, -consdata->lhs) )
457  {
458  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest
459  * since we also want to keep activities in consdata up-to-date, we also need to know when the corresponding bound is relaxed */
460  if( consdata->lincoefs[linvarpos] > 0.0 )
461  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
462  else
463  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
464  }
465 
466  SCIP_CALL( SCIPdropVarEvent(scip, consdata->linvars[linvarpos], eventtype, eventhdlr, (SCIP_EVENTDATA*)consdata->lineventdata[linvarpos], consdata->lineventdata[linvarpos]->filterpos) );
467 
468  SCIPfreeBlockMemory(scip, &consdata->lineventdata[linvarpos]); /*lint !e866 */
469 
470  return SCIP_OKAY;
471 }
472 
473 /** catches variable bound change events on a quadratic variable in a quadratic constraint */
474 static
475 SCIP_RETCODE catchQuadVarEvents(
476  SCIP* scip, /**< SCIP data structure */
477  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
478  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
479  int quadvarpos /**< position of variable in quadratic variables array */
480  )
481 {
482  SCIP_CONSDATA* consdata;
483  SCIP_QUADVAREVENTDATA* eventdata;
484  SCIP_EVENTTYPE eventtype;
485 
486  assert(scip != NULL);
487  assert(eventhdlr != NULL);
488  assert(cons != NULL);
489 
490  consdata = SCIPconsGetData(cons);
491  assert(consdata != NULL);
492 
493  assert(quadvarpos >= 0);
494  assert(quadvarpos < consdata->nquadvars);
495  assert(consdata->quadvarterms[quadvarpos].eventdata == NULL);
496 
497  SCIP_CALL( SCIPallocBlockMemory(scip, &eventdata) );
498 
499  eventtype = SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
500 #ifdef CHECKIMPLINBILINEAR
501  eventtype |= SCIP_EVENTTYPE_IMPLADDED;
502 #endif
503  eventdata->cons = cons;
504  eventdata->varidx = -quadvarpos-1;
505  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->quadvarterms[quadvarpos].var, eventtype, eventhdlr, (SCIP_EVENTDATA*)eventdata, &eventdata->filterpos) );
506 
507  consdata->quadvarterms[quadvarpos].eventdata = eventdata;
508 
509  /* invalidate activity information
510  * NOTE: It could happen that a constraint gets temporary deactivated and some variable bounds change. In this case
511  * we do not recognize those bound changes with the variable events and thus we have to recompute the activities.
512  */
513  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
514 
515  return SCIP_OKAY;
516 }
517 
518 /** catches variable bound change events on a quadratic variable in a quadratic constraint */
519 static
520 SCIP_RETCODE dropQuadVarEvents(
521  SCIP* scip, /**< SCIP data structure */
522  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
523  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
524  int quadvarpos /**< position of variable in quadratic variables array */
525  )
526 {
527  SCIP_CONSDATA* consdata;
528  SCIP_EVENTTYPE eventtype;
529 
530  assert(scip != NULL);
531  assert(eventhdlr != NULL);
532  assert(cons != NULL);
533 
534  consdata = SCIPconsGetData(cons);
535  assert(consdata != NULL);
536 
537  assert(quadvarpos >= 0);
538  assert(quadvarpos < consdata->nquadvars);
539  assert(consdata->quadvarterms[quadvarpos].eventdata != NULL);
540  assert(consdata->quadvarterms[quadvarpos].eventdata->cons == cons);
541  assert(consdata->quadvarterms[quadvarpos].eventdata->varidx == -quadvarpos-1);
542  assert(consdata->quadvarterms[quadvarpos].eventdata->filterpos >= 0);
543 
544  eventtype = SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED | SCIP_EVENTTYPE_GBDCHANGED;
545 #ifdef CHECKIMPLINBILINEAR
546  eventtype |= SCIP_EVENTTYPE_IMPLADDED;
547 #endif
548 
549  SCIP_CALL( SCIPdropVarEvent(scip, consdata->quadvarterms[quadvarpos].var, eventtype, eventhdlr, (SCIP_EVENTDATA*)consdata->quadvarterms[quadvarpos].eventdata, consdata->quadvarterms[quadvarpos].eventdata->filterpos) );
550 
551  SCIPfreeBlockMemory(scip, &consdata->quadvarterms[quadvarpos].eventdata);
552 
553  return SCIP_OKAY;
554 }
555 
556 /** catch variable events */
557 static
558 SCIP_RETCODE catchVarEvents(
559  SCIP* scip, /**< SCIP data structure */
560  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
561  SCIP_CONS* cons /**< constraint for which to catch bound change events */
562  )
563 {
564  SCIP_CONSDATA* consdata;
565  SCIP_VAR* var;
566  int i;
567 
568  assert(scip != NULL);
569  assert(cons != NULL);
570  assert(eventhdlr != NULL);
571 
572  consdata = SCIPconsGetData(cons);
573  assert(consdata != NULL);
574  assert(consdata->lineventdata == NULL);
575 
576  /* we will update isremovedfixings, so reset it to TRUE first */
577  consdata->isremovedfixings = TRUE;
578 
579  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize) );
580  for( i = 0; i < consdata->nlinvars; ++i )
581  {
582  SCIP_CALL( catchLinearVarEvents(scip, eventhdlr, cons, i) );
583 
584  var = consdata->linvars[i];
585  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
586  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
587  }
588 
589  for( i = 0; i < consdata->nquadvars; ++i )
590  {
591  assert(consdata->quadvarterms[i].eventdata == NULL);
592 
593  SCIP_CALL( catchQuadVarEvents(scip, eventhdlr, cons, i) );
594 
595  var = consdata->quadvarterms[i].var;
596  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
597  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
598  }
599 
600  consdata->ispropagated = FALSE;
601 
602  return SCIP_OKAY;
603 }
604 
605 /** drop variable events */
606 static
607 SCIP_RETCODE dropVarEvents(
608  SCIP* scip, /**< SCIP data structure */
609  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
610  SCIP_CONS* cons /**< constraint for which to drop bound change events */
611  )
612 {
613  SCIP_CONSDATA* consdata;
614  int i;
615 
616  assert(scip != NULL);
617  assert(eventhdlr != NULL);
618  assert(cons != NULL);
619 
620  consdata = SCIPconsGetData(cons);
621  assert(consdata != NULL);
622 
623  if( consdata->lineventdata != NULL )
624  {
625  for( i = 0; i < consdata->nlinvars; ++i )
626  {
627  if( consdata->lineventdata[i] != NULL )
628  {
629  SCIP_CALL( dropLinearVarEvents(scip, eventhdlr, cons, i) );
630  }
631  }
632  SCIPfreeBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize);
633  }
634 
635  for( i = 0; i < consdata->nquadvars; ++i )
636  {
637  if( consdata->quadvarterms[i].eventdata != NULL )
638  {
639  SCIP_CALL( dropQuadVarEvents(scip, eventhdlr, cons, i) );
640  }
641  }
642 
643  return SCIP_OKAY;
644 }
645 
646 /** locks a linear variable in a constraint */
647 static
648 SCIP_RETCODE lockLinearVariable(
649  SCIP* scip, /**< SCIP data structure */
650  SCIP_CONS* cons, /**< constraint where to lock a variable */
651  SCIP_VAR* var, /**< variable to lock */
652  SCIP_Real coef /**< coefficient of variable in constraint */
653  )
654 {
655  SCIP_CONSDATA* consdata;
656 
657  assert(scip != NULL);
658  assert(cons != NULL);
659  assert(var != NULL);
660  assert(coef != 0.0);
661 
662  consdata = SCIPconsGetData(cons);
663  assert(consdata != NULL);
664 
665  if( coef > 0.0 )
666  {
667  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
668  }
669  else
670  {
671  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
672  }
673 
674  return SCIP_OKAY;
675 }
676 
677 /** unlocks a linear variable in a constraint */
678 static
679 SCIP_RETCODE unlockLinearVariable(
680  SCIP* scip, /**< SCIP data structure */
681  SCIP_CONS* cons, /**< constraint where to unlock a variable */
682  SCIP_VAR* var, /**< variable to unlock */
683  SCIP_Real coef /**< coefficient of variable in constraint */
684  )
685 {
686  SCIP_CONSDATA* consdata;
687 
688  assert(scip != NULL);
689  assert(cons != NULL);
690  assert(var != NULL);
691  assert(coef != 0.0);
692 
693  consdata = SCIPconsGetData(cons);
694  assert(consdata != NULL);
695 
696  if( coef > 0.0 )
697  {
698  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
699  }
700  else
701  {
702  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
703  }
704 
705  return SCIP_OKAY;
706 }
707 
708 /** locks a quadratic variable in a constraint */
709 static
710 SCIP_RETCODE lockQuadraticVariable(
711  SCIP* scip, /**< SCIP data structure */
712  SCIP_CONS* cons, /**< constraint where to lock a variable */
713  SCIP_VAR* var /**< variable to lock */
714  )
715 {
716  SCIP_CALL( SCIPlockVarCons(scip, var, cons, TRUE, TRUE) );
717 
718  return SCIP_OKAY;
719 }
720 
721 /** unlocks a quadratic variable in a constraint */
722 static
723 SCIP_RETCODE unlockQuadraticVariable(
724  SCIP* scip, /**< SCIP data structure */
725  SCIP_CONS* cons, /**< constraint where to unlock a variable */
726  SCIP_VAR* var /**< variable to unlock */
727  )
728 {
729  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, TRUE, TRUE) );
730 
731  return SCIP_OKAY;
732 }
733 
734 /** computes the minimal and maximal activity for the linear part in a constraint data
735  *
736  * Only sums up terms that contribute finite values.
737  * Gives the number of terms that contribute infinite values.
738  * Only computes those activities where the corresponding side of the constraint is finite.
739  */
740 static
741 void consdataUpdateLinearActivity(
742  SCIP* scip, /**< SCIP data structure */
743  SCIP_CONSDATA* consdata, /**< constraint data */
744  SCIP_Real intervalinfty /**< infinity value used in interval operations */
745  )
746 { /*lint --e{666}*/
747  SCIP_ROUNDMODE prevroundmode;
748  int i;
749  SCIP_Real bnd;
750 
751  assert(scip != NULL);
752  assert(consdata != NULL);
753 
754  /* if variable bounds are not strictly consistent, then the activity update methods may yield inconsistent activities
755  * in this case, we also recompute the activities
756  */
757  if( consdata->minlinactivity != SCIP_INVALID && consdata->maxlinactivity != SCIP_INVALID && /*lint !e777 */
758  (consdata->minlinactivityinf > 0 || consdata->maxlinactivityinf > 0 || consdata->minlinactivity <= consdata->maxlinactivity) )
759  {
760  /* activities should be up-to-date */
761  assert(consdata->minlinactivityinf >= 0);
762  assert(consdata->maxlinactivityinf >= 0);
763  return;
764  }
765 
766  consdata->minlinactivityinf = 0;
767  consdata->maxlinactivityinf = 0;
768 
769  /* if lhs is -infinite, then we do not compute a maximal activity, so we set it to infinity
770  * if rhs is infinite, then we do not compute a minimal activity, so we set it to -infinity
771  */
772  consdata->minlinactivity = SCIPisInfinity(scip, consdata->rhs) ? -intervalinfty : 0.0;
773  consdata->maxlinactivity = SCIPisInfinity(scip, -consdata->lhs) ? intervalinfty : 0.0;
774 
775  if( consdata->nlinvars == 0 )
776  return;
777 
778  /* if the activities computed here should be still up-to-date after bound changes,
779  * variable events need to be caught */
780  assert(consdata->lineventdata != NULL);
781 
782  prevroundmode = SCIPintervalGetRoundingMode();
783 
784  if( !SCIPisInfinity(scip, consdata->rhs) )
785  {
786  /* compute minimal activity only if there is a finite right hand side */
787  SCIPintervalSetRoundingModeDownwards();
788 
789  for( i = 0; i < consdata->nlinvars; ++i )
790  {
791  assert(consdata->lineventdata[i] != NULL);
792  if( consdata->lincoefs[i] >= 0.0 )
793  {
794  bnd = MIN(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
795  if( SCIPisInfinity(scip, -bnd) )
796  {
797  ++consdata->minlinactivityinf;
798  continue;
799  }
800  assert(!SCIPisInfinity(scip, bnd)); /* do not like variables that are fixed at +infinity */
801  }
802  else
803  {
804  bnd = MAX(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
805  if( SCIPisInfinity(scip, bnd) )
806  {
807  ++consdata->minlinactivityinf;
808  continue;
809  }
810  assert(!SCIPisInfinity(scip, -bnd)); /* do not like variables that are fixed at -infinity */
811  }
812  consdata->minlinactivity += consdata->lincoefs[i] * bnd;
813  }
814  }
815 
816  if( !SCIPisInfinity(scip, -consdata->lhs) )
817  {
818  /* compute maximal activity only if there is a finite left hand side */
819  SCIPintervalSetRoundingModeUpwards();
820 
821  for( i = 0; i < consdata->nlinvars; ++i )
822  {
823  assert(consdata->lineventdata[i] != NULL);
824  if( consdata->lincoefs[i] >= 0.0 )
825  {
826  bnd = MAX(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
827  if( SCIPisInfinity(scip, bnd) )
828  {
829  ++consdata->maxlinactivityinf;
830  continue;
831  }
832  assert(!SCIPisInfinity(scip, -bnd)); /* do not like variables that are fixed at -infinity */
833  }
834  else
835  {
836  bnd = MIN(SCIPvarGetLbLocal(consdata->linvars[i]), SCIPvarGetUbLocal(consdata->linvars[i]));
837  if( SCIPisInfinity(scip, -bnd) )
838  {
839  ++consdata->maxlinactivityinf;
840  continue;
841  }
842  assert(!SCIPisInfinity(scip, bnd)); /* do not like variables that are fixed at +infinity */
843  }
844  consdata->maxlinactivity += consdata->lincoefs[i] * bnd;
845  }
846  }
847 
848  SCIPintervalSetRoundingMode(prevroundmode);
849 
850  assert(consdata->minlinactivityinf > 0 || consdata->maxlinactivityinf > 0 || consdata->minlinactivity <= consdata->maxlinactivity);
851 }
852 
853 /** update the linear activities after a change in the lower bound of a variable */
854 static
855 void consdataUpdateLinearActivityLbChange(
856  SCIP* scip, /**< SCIP data structure */
857  SCIP_CONSDATA* consdata, /**< constraint data */
858  SCIP_Real coef, /**< coefficient of variable in constraint */
859  SCIP_Real oldbnd, /**< previous lower bound of variable */
860  SCIP_Real newbnd /**< new lower bound of variable */
861  )
862 {
863  SCIP_ROUNDMODE prevroundmode;
864 
865  assert(scip != NULL);
866  assert(consdata != NULL);
867  /* we can't deal with lower bounds at infinity */
868  assert(!SCIPisInfinity(scip, oldbnd));
869  assert(!SCIPisInfinity(scip, newbnd));
870 
871  /* @todo since we check the linear activity for consistency later anyway, we may skip changing the rounding mode here */
872 
873  /* assume lhs <= a*x + y <= rhs, then the following bound changes can be deduced:
874  * a > 0: y <= rhs - a*lb(x), y >= lhs - a*ub(x)
875  * a < 0: y <= rhs - a*ub(x), y >= lhs - a*lb(x)
876  */
877 
878  if( coef > 0.0 )
879  {
880  /* we should only be called if rhs is finite */
881  assert(!SCIPisInfinity(scip, consdata->rhs));
882 
883  /* we have no min activities computed so far, so cannot update */
884  if( consdata->minlinactivity == SCIP_INVALID ) /*lint !e777 */
885  return;
886 
887  prevroundmode = SCIPintervalGetRoundingMode();
888  SCIPintervalSetRoundingModeDownwards();
889 
890  /* update min activity */
891  if( SCIPisInfinity(scip, -oldbnd) )
892  {
893  --consdata->minlinactivityinf;
894  assert(consdata->minlinactivityinf >= 0);
895  }
896  else
897  {
898  SCIP_Real minuscoef;
899  minuscoef = -coef;
900  consdata->minlinactivity += minuscoef * oldbnd;
901  }
902 
903  if( SCIPisInfinity(scip, -newbnd) )
904  {
905  ++consdata->minlinactivityinf;
906  }
907  else
908  {
909  consdata->minlinactivity += coef * newbnd;
910  }
911 
912  SCIPintervalSetRoundingMode(prevroundmode);
913  }
914  else
915  {
916  /* we should only be called if lhs is finite */
917  assert(!SCIPisInfinity(scip, -consdata->lhs));
918 
919  /* we have no max activities computed so far, so cannot update */
920  if( consdata->maxlinactivity == SCIP_INVALID ) /*lint !e777 */
921  return;
922 
923  prevroundmode = SCIPintervalGetRoundingMode();
924  SCIPintervalSetRoundingModeUpwards();
925 
926  /* update max activity */
927  if( SCIPisInfinity(scip, -oldbnd) )
928  {
929  --consdata->maxlinactivityinf;
930  assert(consdata->maxlinactivityinf >= 0);
931  }
932  else
933  {
934  SCIP_Real minuscoef;
935  minuscoef = -coef;
936  consdata->maxlinactivity += minuscoef * oldbnd;
937  }
938 
939  if( SCIPisInfinity(scip, -newbnd) )
940  {
941  ++consdata->maxlinactivityinf;
942  }
943  else
944  {
945  consdata->maxlinactivity += coef * newbnd;
946  }
947 
948  SCIPintervalSetRoundingMode(prevroundmode);
949  }
950 }
951 
952 /** update the linear activities after a change in the upper bound of a variable */
953 static
954 void consdataUpdateLinearActivityUbChange(
955  SCIP* scip, /**< SCIP data structure */
956  SCIP_CONSDATA* consdata, /**< constraint data */
957  SCIP_Real coef, /**< coefficient of variable in constraint */
958  SCIP_Real oldbnd, /**< previous lower bound of variable */
959  SCIP_Real newbnd /**< new lower bound of variable */
960  )
961 {
962  SCIP_ROUNDMODE prevroundmode;
963 
964  assert(scip != NULL);
965  assert(consdata != NULL);
966  /* we can't deal with upper bounds at -infinity */
967  assert(!SCIPisInfinity(scip, -oldbnd));
968  assert(!SCIPisInfinity(scip, -newbnd));
969 
970  /* @todo since we check the linear activity for consistency later anyway, we may skip changing the rounding mode here */
971 
972  /* assume lhs <= a*x + y <= rhs, then the following bound changes can be deduced:
973  * a > 0: y <= rhs - a*lb(x), y >= lhs - a*ub(x)
974  * a < 0: y <= rhs - a*ub(x), y >= lhs - a*lb(x)
975  */
976 
977  if( coef > 0.0 )
978  {
979  /* we should only be called if lhs is finite */
980  assert(!SCIPisInfinity(scip, -consdata->lhs));
981 
982  /* we have no max activities computed so far, so cannot update */
983  if( consdata->maxlinactivity == SCIP_INVALID ) /*lint !e777 */
984  return;
985 
986  prevroundmode = SCIPintervalGetRoundingMode();
987  SCIPintervalSetRoundingModeUpwards();
988 
989  /* update max activity */
990  if( SCIPisInfinity(scip, oldbnd) )
991  {
992  --consdata->maxlinactivityinf;
993  assert(consdata->maxlinactivityinf >= 0);
994  }
995  else
996  {
997  SCIP_Real minuscoef;
998  minuscoef = -coef;
999  consdata->maxlinactivity += minuscoef * oldbnd;
1000  }
1001 
1002  if( SCIPisInfinity(scip, newbnd) )
1003  {
1004  ++consdata->maxlinactivityinf;
1005  }
1006  else
1007  {
1008  consdata->maxlinactivity += coef * newbnd;
1009  }
1010 
1011  SCIPintervalSetRoundingMode(prevroundmode);
1012  }
1013  else
1014  {
1015  /* we should only be called if rhs is finite */
1016  assert(!SCIPisInfinity(scip, consdata->rhs));
1017 
1018  /* we have no min activities computed so far, so cannot update */
1019  if( consdata->minlinactivity == SCIP_INVALID ) /*lint !e777 */
1020  return;
1021 
1022  prevroundmode = SCIPintervalGetRoundingMode();
1023  SCIPintervalSetRoundingModeDownwards();
1024 
1025  /* update min activity */
1026  if( SCIPisInfinity(scip, oldbnd) )
1027  {
1028  --consdata->minlinactivityinf;
1029  assert(consdata->minlinactivityinf >= 0);
1030  }
1031  else
1032  {
1033  SCIP_Real minuscoef;
1034  minuscoef = -coef;
1035  consdata->minlinactivity += minuscoef * oldbnd;
1036  }
1037 
1038  if( SCIPisInfinity(scip, newbnd) )
1039  {
1040  ++consdata->minlinactivityinf;
1041  }
1042  else
1043  {
1044  consdata->minlinactivity += coef * newbnd;
1045  }
1046 
1047  SCIPintervalSetRoundingMode(prevroundmode);
1048  }
1049 }
1050 
1051 /** returns whether a quadratic variable domain can be reduced to its lower or upper bound; this is the case if the
1052  * quadratic variable is in just one single quadratic constraint and (sqrcoef > 0 and LHS = -infinity), or
1053  * (sqrcoef < 0 and RHS = +infinity) hold
1054  */
1055 static
1056 SCIP_Bool hasQuadvarHpProperty(
1057  SCIP* scip, /**< SCIP data structure */
1058  SCIP_CONSDATA* consdata, /**< constraint data */
1059  int idx /**< index of quadratic variable */
1060  )
1061 {
1062  SCIP_VAR* var;
1063  SCIP_Real quadcoef;
1064  SCIP_Bool haslhs;
1065  SCIP_Bool hasrhs;
1066 
1067  assert(scip != NULL);
1068  assert(consdata != NULL);
1069  assert(idx >= 0 && idx < consdata->nquadvars);
1070 
1071  var = consdata->quadvarterms[idx].var;
1072  assert(var != NULL);
1073 
1074  quadcoef = consdata->quadvarterms[idx].sqrcoef;
1075  haslhs = !SCIPisInfinity(scip, -consdata->lhs);
1076  hasrhs = !SCIPisInfinity(scip, consdata->rhs);
1077 
1078  return SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == 1
1079  && SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == 1 && SCIPisZero(scip, SCIPvarGetObj(var))
1080  && SCIPvarGetType(var) != SCIP_VARTYPE_BINARY && ((quadcoef < 0.0 && !haslhs) || (quadcoef > 0.0 && !hasrhs));
1081 }
1082 
1083 /** processes variable fixing or bound change event */
1084 static
1085 SCIP_DECL_EVENTEXEC(processVarEvent)
1086 {
1087  SCIP_CONS* cons;
1088  SCIP_CONSDATA* consdata;
1089  SCIP_EVENTTYPE eventtype;
1090  int varidx;
1091 
1092  assert(scip != NULL);
1093  assert(event != NULL);
1094  assert(eventdata != NULL);
1095  assert(eventhdlr != NULL);
1096 
1097  cons = ((SCIP_QUADVAREVENTDATA*)eventdata)->cons;
1098  assert(cons != NULL);
1099  consdata = SCIPconsGetData(cons);
1100  assert(consdata != NULL);
1101 
1102  varidx = ((SCIP_QUADVAREVENTDATA*)eventdata)->varidx;
1103  assert(varidx < 0 || varidx < consdata->nlinvars);
1104  assert(varidx >= 0 || -varidx-1 < consdata->nquadvars);
1105 
1106  eventtype = SCIPeventGetType(event);
1107 
1108  /* process local bound changes */
1109  if( eventtype & SCIP_EVENTTYPE_BOUNDCHANGED )
1110  {
1111  if( varidx < 0 )
1112  {
1113  /* mark activity bounds for quad term as not up to date anymore */
1114  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
1115  }
1116  else
1117  {
1118  /* update activity bounds for linear terms */
1119  if( eventtype & SCIP_EVENTTYPE_LBCHANGED )
1120  consdataUpdateLinearActivityLbChange(scip, consdata, consdata->lincoefs[varidx], SCIPeventGetOldbound(event), SCIPeventGetNewbound(event));
1121  else
1122  consdataUpdateLinearActivityUbChange(scip, consdata, consdata->lincoefs[varidx], SCIPeventGetOldbound(event), SCIPeventGetNewbound(event));
1123  }
1124 
1125  if( eventtype & SCIP_EVENTTYPE_BOUNDTIGHTENED )
1126  {
1127  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
1128  consdata->ispropagated = FALSE;
1129  }
1130  }
1131 
1132  /* process global bound changes */
1133  if( eventtype & SCIP_EVENTTYPE_GBDCHANGED )
1134  {
1135  SCIP_VAR* var;
1136 
1137  var = varidx < 0 ? consdata->quadvarterms[-varidx-1].var : consdata->linvars[varidx];
1138  assert(var != NULL);
1139 
1140  if( varidx < 0 )
1141  {
1142  SCIP_QUADVARTERM* quadvarterm;
1143 
1144  quadvarterm = &consdata->quadvarterms[-varidx-1];
1145 
1146  /* if an integer variable x with a x^2 is tightened to [0,1], then we can replace the x^2 by x, which is done in mergeAndCleanQuadVarTerms()
1147  * we currently do this only if the binary variable does not show up in any bilinear terms
1148  */
1149  if( SCIPgetStage(scip) < SCIP_STAGE_SOLVING && SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER && SCIPvarIsBinary(var) &&
1150  quadvarterm->sqrcoef != 0.0 && quadvarterm->nadjbilin == 0 )
1151  {
1152  consdata->quadvarsmerged = FALSE;
1153  consdata->initialmerge = FALSE;
1154  }
1155  }
1156 
1157  if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
1158  consdata->isremovedfixings = FALSE;
1159  }
1160 
1161  /* process variable fixing event */
1162  if( eventtype & SCIP_EVENTTYPE_VARFIXED )
1163  {
1164  consdata->isremovedfixings = FALSE;
1165  }
1166 
1167 #ifdef CHECKIMPLINBILINEAR
1168  if( eventtype & SCIP_EVENTTYPE_IMPLADDED )
1169  {
1170  assert(varidx < 0); /* we catch impladded events only for quadratic variables */
1171  /* if variable is binary (quite likely if an implication has been added) and occurs in a bilinear term, then mark that we should check implications */
1172  if( SCIPvarIsBinary(SCIPeventGetVar(event)) && consdata->quadvarterms[-varidx-1].nadjbilin > 0 )
1173  consdata->isimpladded = TRUE;
1174  }
1175 #endif
1176 
1177  return SCIP_OKAY;
1178 }
1179 
1180 /** ensures, that linear vars and coefs arrays can store at least num entries */
1181 static
1182 SCIP_RETCODE consdataEnsureLinearVarsSize(
1183  SCIP* scip, /**< SCIP data structure */
1184  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1185  int num /**< minimum number of entries to store */
1186  )
1187 {
1188  assert(scip != NULL);
1189  assert(consdata != NULL);
1190  assert(consdata->nlinvars <= consdata->linvarssize);
1191 
1192  if( num > consdata->linvarssize )
1193  {
1194  int newsize;
1195 
1196  newsize = SCIPcalcMemGrowSize(scip, num);
1197  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->linvars, consdata->linvarssize, newsize) );
1198  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->lincoefs, consdata->linvarssize, newsize) );
1199  if( consdata->lineventdata != NULL )
1200  {
1201  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize, newsize) );
1202  }
1203  consdata->linvarssize = newsize;
1204  }
1205  assert(num <= consdata->linvarssize);
1206 
1207  return SCIP_OKAY;
1208 }
1209 
1210 /** ensures, that quadratic variable terms array can store at least num entries */
1211 static
1212 SCIP_RETCODE consdataEnsureQuadVarTermsSize(
1213  SCIP* scip, /**< SCIP data structure */
1214  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1215  int num /**< minimum number of entries to store */
1216  )
1217 {
1218  assert(scip != NULL);
1219  assert(consdata != NULL);
1220  assert(consdata->nquadvars <= consdata->quadvarssize);
1221 
1222  if( num > consdata->quadvarssize )
1223  {
1224  int newsize;
1225 
1226  newsize = SCIPcalcMemGrowSize(scip, num);
1227  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->quadvarterms, consdata->quadvarssize, newsize) );
1228  consdata->quadvarssize = newsize;
1229  }
1230  assert(num <= consdata->quadvarssize);
1231 
1232  return SCIP_OKAY;
1233 }
1234 
1235 /** ensures, that adjacency array can store at least num entries */
1236 static
1237 SCIP_RETCODE consdataEnsureAdjBilinSize(
1238  SCIP* scip, /**< SCIP data structure */
1239  SCIP_QUADVARTERM* quadvarterm, /**< quadratic variable term */
1240  int num /**< minimum number of entries to store */
1241  )
1242 {
1243  assert(scip != NULL);
1244  assert(quadvarterm != NULL);
1245  assert(quadvarterm->nadjbilin <= quadvarterm->adjbilinsize);
1246 
1247  if( num > quadvarterm->adjbilinsize )
1248  {
1249  int newsize;
1250 
1251  newsize = SCIPcalcMemGrowSize(scip, num);
1252  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &quadvarterm->adjbilin, quadvarterm->adjbilinsize, newsize) );
1253  quadvarterm->adjbilinsize = newsize;
1254  }
1255  assert(num <= quadvarterm->adjbilinsize);
1256 
1257  return SCIP_OKAY;
1258 }
1259 
1260 /** ensures, that bilinear term arrays can store at least num entries */
1261 static
1262 SCIP_RETCODE consdataEnsureBilinSize(
1263  SCIP* scip, /**< SCIP data structure */
1264  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1265  int num /**< minimum number of entries to store */
1266  )
1267 {
1268  assert(scip != NULL);
1269  assert(consdata != NULL);
1270  assert(consdata->nbilinterms <= consdata->bilintermssize);
1271 
1272  if( num > consdata->bilintermssize )
1273  {
1274  int newsize;
1275 
1276  newsize = SCIPcalcMemGrowSize(scip, num);
1277  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->bilinterms, consdata->bilintermssize, newsize) );
1278  consdata->bilintermssize = newsize;
1279  }
1280  assert(num <= consdata->bilintermssize);
1281 
1282  return SCIP_OKAY;
1283 }
1284 
1285 /** creates empty constraint data structure */
1286 static
1287 SCIP_RETCODE consdataCreateEmpty(
1288  SCIP* scip, /**< SCIP data structure */
1289  SCIP_CONSDATA** consdata /**< a buffer to store pointer to new constraint data */
1290  )
1291 {
1292  assert(scip != NULL);
1293  assert(consdata != NULL);
1294 
1295  SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );
1296  BMSclearMemory(*consdata);
1297 
1298  (*consdata)->lhs = -SCIPinfinity(scip);
1299  (*consdata)->rhs = SCIPinfinity(scip);
1300 
1301  (*consdata)->linvarssorted = TRUE;
1302  (*consdata)->linvarsmerged = TRUE;
1303  (*consdata)->quadvarssorted = TRUE;
1304  (*consdata)->quadvarsmerged = TRUE;
1305  (*consdata)->bilinsorted = TRUE;
1306  (*consdata)->bilinmerged = TRUE;
1307 
1308  (*consdata)->isremovedfixings = TRUE;
1309  (*consdata)->ispropagated = TRUE;
1310  (*consdata)->initialmerge = FALSE;
1311 
1312  (*consdata)->linvar_maydecrease = -1;
1313  (*consdata)->linvar_mayincrease = -1;
1314 
1315  (*consdata)->minlinactivity = SCIP_INVALID;
1316  (*consdata)->maxlinactivity = SCIP_INVALID;
1317  (*consdata)->minlinactivityinf = -1;
1318  (*consdata)->maxlinactivityinf = -1;
1319 
1320  (*consdata)->isgaugeavailable = FALSE;
1321  (*consdata)->isedavailable = FALSE;
1322 
1323  return SCIP_OKAY;
1324 }
1325 
1326 /** creates constraint data structure */
1327 static
1328 SCIP_RETCODE consdataCreate(
1329  SCIP* scip, /**< SCIP data structure */
1330  SCIP_CONSDATA** consdata, /**< a buffer to store pointer to new constraint data */
1331  SCIP_Real lhs, /**< left hand side of constraint */
1332  SCIP_Real rhs, /**< right hand side of constraint */
1333  int nlinvars, /**< number of linear variables */
1334  SCIP_VAR** linvars, /**< array of linear variables */
1335  SCIP_Real* lincoefs, /**< array of coefficients of linear variables */
1336  int nquadvars, /**< number of quadratic variables */
1337  SCIP_QUADVARTERM* quadvarterms, /**< array of quadratic variable terms */
1338  int nbilinterms, /**< number of bilinear terms */
1339  SCIP_BILINTERM* bilinterms, /**< array of bilinear terms */
1340  SCIP_Bool capturevars /**< whether we should capture variables */
1341  )
1342 {
1343  int i;
1344 
1345  assert(scip != NULL);
1346  assert(consdata != NULL);
1347 
1348  assert(nlinvars == 0 || linvars != NULL);
1349  assert(nlinvars == 0 || lincoefs != NULL);
1350  assert(nquadvars == 0 || quadvarterms != NULL);
1351  assert(nbilinterms == 0 || bilinterms != NULL);
1352 
1353  SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );
1354  BMSclearMemory(*consdata);
1355 
1356  (*consdata)->minlinactivity = SCIP_INVALID;
1357  (*consdata)->maxlinactivity = SCIP_INVALID;
1358  (*consdata)->minlinactivityinf = -1;
1359  (*consdata)->maxlinactivityinf = -1;
1360 
1361  (*consdata)->lhs = lhs;
1362  (*consdata)->rhs = rhs;
1363 
1364  if( nlinvars > 0 )
1365  {
1366  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->linvars, linvars, nlinvars) );
1367  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->lincoefs, lincoefs, nlinvars) );
1368  (*consdata)->nlinvars = nlinvars;
1369  (*consdata)->linvarssize = nlinvars;
1370 
1371  if( capturevars )
1372  for( i = 0; i < nlinvars; ++i )
1373  {
1374  SCIP_CALL( SCIPcaptureVar(scip, linvars[i]) );
1375  }
1376  }
1377  else
1378  {
1379  (*consdata)->linvarssorted = TRUE;
1380  (*consdata)->linvarsmerged = TRUE;
1381  (*consdata)->minlinactivity = 0.0;
1382  (*consdata)->maxlinactivity = 0.0;
1383  (*consdata)->minlinactivityinf = 0;
1384  (*consdata)->maxlinactivityinf = 0;
1385  }
1386 
1387  if( nquadvars > 0 )
1388  {
1389  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->quadvarterms, quadvarterms, nquadvars) );
1390 
1391  for( i = 0; i < nquadvars; ++i )
1392  {
1393  (*consdata)->quadvarterms[i].eventdata = NULL;
1394  if( quadvarterms[i].nadjbilin )
1395  {
1396  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->quadvarterms[i].adjbilin, quadvarterms[i].adjbilin, quadvarterms[i].nadjbilin) );
1397  (*consdata)->quadvarterms[i].adjbilinsize = quadvarterms[i].nadjbilin;
1398  }
1399  else
1400  {
1401  assert((*consdata)->quadvarterms[i].nadjbilin == 0);
1402  (*consdata)->quadvarterms[i].adjbilin = NULL;
1403  (*consdata)->quadvarterms[i].adjbilinsize = 0;
1404  }
1405  if( capturevars )
1406  {
1407  SCIP_CALL( SCIPcaptureVar(scip, quadvarterms[i].var) );
1408  }
1409  }
1410 
1411  (*consdata)->nquadvars = nquadvars;
1412  (*consdata)->quadvarssize = nquadvars;
1413  SCIPintervalSetEmpty(&(*consdata)->quadactivitybounds);
1414  }
1415  else
1416  {
1417  (*consdata)->quadvarssorted = TRUE;
1418  (*consdata)->quadvarsmerged = TRUE;
1419  SCIPintervalSet(&(*consdata)->quadactivitybounds, 0.0);
1420  }
1421 
1422  if( nbilinterms > 0 )
1423  {
1424  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->bilinterms, bilinterms, nbilinterms) );
1425  (*consdata)->nbilinterms = nbilinterms;
1426  (*consdata)->bilintermssize = nbilinterms;
1427  }
1428  else
1429  {
1430  (*consdata)->bilinsorted = TRUE;
1431  (*consdata)->bilinmerged = TRUE;
1432  }
1433 
1434  (*consdata)->linvar_maydecrease = -1;
1435  (*consdata)->linvar_mayincrease = -1;
1436 
1437  (*consdata)->activity = SCIP_INVALID;
1438  (*consdata)->lhsviol = SCIPisInfinity(scip, -lhs) ? 0.0 : SCIP_INVALID;
1439  (*consdata)->rhsviol = SCIPisInfinity(scip, rhs) ? 0.0 : SCIP_INVALID;
1440 
1441  (*consdata)->isgaugeavailable = FALSE;
1442 
1443  return SCIP_OKAY;
1444 }
1445 
1446 /** frees constraint data structure */
1447 static
1448 SCIP_RETCODE consdataFree(
1449  SCIP* scip, /**< SCIP data structure */
1450  SCIP_CONSDATA** consdata /**< pointer to constraint data to free */
1451  )
1452 {
1453  int i;
1454 
1455  assert(scip != NULL);
1456  assert(consdata != NULL);
1457  assert(*consdata != NULL);
1458 
1459  /* free sepa arrays, may exists if constraint is deleted in solving stage */
1460  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->sepaquadvars, (*consdata)->nquadvars);
1461  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->sepabilinvar2pos, (*consdata)->nbilinterms);
1462 
1463  /* release linear variables and free linear part */
1464  if( (*consdata)->linvarssize > 0 )
1465  {
1466  for( i = 0; i < (*consdata)->nlinvars; ++i )
1467  {
1468  assert((*consdata)->lineventdata == NULL || (*consdata)->lineventdata[i] == NULL); /* variable events should have been dropped earlier */
1469  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->linvars[i]) );
1470  }
1471  SCIPfreeBlockMemoryArray(scip, &(*consdata)->linvars, (*consdata)->linvarssize);
1472  SCIPfreeBlockMemoryArray(scip, &(*consdata)->lincoefs, (*consdata)->linvarssize);
1473  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->lineventdata, (*consdata)->linvarssize);
1474  }
1475  assert((*consdata)->linvars == NULL);
1476  assert((*consdata)->lincoefs == NULL);
1477  assert((*consdata)->lineventdata == NULL);
1478 
1479  /* release quadratic variables and free quadratic variable term part */
1480  for( i = 0; i < (*consdata)->nquadvars; ++i )
1481  {
1482  assert((*consdata)->quadvarterms[i].eventdata == NULL); /* variable events should have been dropped earlier */
1483  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->quadvarterms[i].adjbilin, (*consdata)->quadvarterms[i].adjbilinsize);
1484  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->quadvarterms[i].var) );
1485  }
1486  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->quadvarterms, (*consdata)->quadvarssize);
1487 
1488  /* free bilinear terms */
1489  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->bilinterms, (*consdata)->bilintermssize);
1490 
1491  /* free nonlinear row representation */
1492  if( (*consdata)->nlrow != NULL )
1493  {
1494  SCIP_CALL( SCIPreleaseNlRow(scip, &(*consdata)->nlrow) );
1495  }
1496 
1497  /* free interior point information, may exists if constraint is deleted in solving stage */
1498  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->interiorpoint, (*consdata)->nquadvars);
1499  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->gaugecoefs, (*consdata)->nquadvars);
1500 
1501  /* free eigen decomposition information */
1502  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->eigenvalues, (*consdata)->nquadvars);
1503  if( (*consdata)->eigenvectors != NULL ) /* explicit check on NULL to avoid runtime warning if nquadvars^2 > int_max */
1504  SCIPfreeBlockMemoryArray(scip, &(*consdata)->eigenvectors, (int)((*consdata)->nquadvars*(*consdata)->nquadvars));
1505  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->bp, (*consdata)->nquadvars);
1506 
1507  /* free unique indices of bilinear terms array */
1508  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->bilintermsidx, (*consdata)->nbilinterms);
1509 
1510  SCIPfreeBlockMemory(scip, consdata);
1511  *consdata = NULL;
1512 
1513  return SCIP_OKAY;
1514 }
1515 
1516 /** sorts linear part of constraint data */
1517 static
1518 void consdataSortLinearVars(
1519  SCIP_CONSDATA* consdata /**< quadratic constraint data */
1520  )
1521 {
1522  assert(consdata != NULL);
1523 
1524  if( consdata->linvarssorted )
1525  return;
1526 
1527  if( consdata->nlinvars <= 1 )
1528  {
1529  consdata->linvarssorted = TRUE;
1530  return;
1531  }
1532 
1533  if( consdata->lineventdata == NULL )
1534  {
1535  SCIPsortPtrReal((void**)consdata->linvars, consdata->lincoefs, SCIPvarComp, consdata->nlinvars);
1536  }
1537  else
1538  {
1539  int i;
1540 
1541  SCIPsortPtrPtrReal((void**)consdata->linvars, (void**)consdata->lineventdata, consdata->lincoefs, SCIPvarComp, consdata->nlinvars);
1542 
1543  /* update variable indices in event data */
1544  for( i = 0; i < consdata->nlinvars; ++i )
1545  if( consdata->lineventdata[i] != NULL )
1546  consdata->lineventdata[i]->varidx = i;
1547  }
1548 
1549  consdata->linvarssorted = TRUE;
1550 }
1551 
1552 #ifdef SCIP_DISABLED_CODE /* no-one needs this routine currently */
1553 /** returns the position of variable in the linear coefficients array of a constraint, or -1 if not found */
1554 static
1555 int consdataFindLinearVar(
1556  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1557  SCIP_VAR* var /**< variable to search for */
1558  )
1559 {
1560  int pos;
1561 
1562  assert(consdata != NULL);
1563  assert(var != NULL);
1564 
1565  if( consdata->nlinvars == 0 )
1566  return -1;
1567 
1568  consdataSortLinearVars(consdata);
1569 
1570  if( !SCIPsortedvecFindPtr((void**)consdata->linvars, SCIPvarComp, (void*)var, consdata->nlinvars, &pos) )
1571  pos = -1;
1572 
1573  return pos;
1574 }
1575 #endif
1576 
1577 /** index comparison method for quadratic variable terms: compares two indices of the quadratic variable set in the quadratic constraint */
1578 static
1579 SCIP_DECL_SORTINDCOMP(quadVarTermComp)
1580 { /*lint --e{715}*/
1581  SCIP_CONSDATA* consdata = (SCIP_CONSDATA*)dataptr;
1582 
1583  assert(consdata != NULL);
1584  assert(0 <= ind1 && ind1 < consdata->nquadvars);
1585  assert(0 <= ind2 && ind2 < consdata->nquadvars);
1586 
1587  return SCIPvarCompare(consdata->quadvarterms[ind1].var, consdata->quadvarterms[ind2].var);
1588 }
1589 
1590 /** sorting of quadratic variable terms */
1591 static
1592 SCIP_RETCODE consdataSortQuadVarTerms(
1593  SCIP* scip, /**< SCIP data structure */
1594  SCIP_CONSDATA* consdata /**< quadratic constraint data */
1595  )
1596 {
1597  int* perm;
1598  int i;
1599  int nexti;
1600  int v;
1601  SCIP_QUADVARTERM quadterm;
1602 
1603  assert(scip != NULL);
1604  assert(consdata != NULL);
1605 
1606  if( consdata->quadvarssorted )
1607  return SCIP_OKAY;
1608 
1609  if( consdata->nquadvars == 0 )
1610  {
1611  consdata->quadvarssorted = TRUE;
1612  return SCIP_OKAY;
1613  }
1614 
1615  /* get temporary memory to store the sorted permutation */
1616  SCIP_CALL( SCIPallocBufferArray(scip, &perm, consdata->nquadvars) );
1617 
1618  /* call bubble sort */
1619  SCIPsort(perm, quadVarTermComp, (void*)consdata, consdata->nquadvars);
1620 
1621  /* permute the quadratic variable terms according to the resulting permutation */
1622  for( v = 0; v < consdata->nquadvars; ++v )
1623  {
1624  if( perm[v] != v )
1625  {
1626  quadterm = consdata->quadvarterms[v];
1627 
1628  i = v;
1629  do
1630  {
1631  assert(0 <= perm[i] && perm[i] < consdata->nquadvars);
1632  assert(perm[i] != i);
1633  consdata->quadvarterms[i] = consdata->quadvarterms[perm[i]];
1634  if( consdata->quadvarterms[i].eventdata != NULL )
1635  {
1636  consdata->quadvarterms[i].eventdata->varidx = -i-1;
1637  }
1638  nexti = perm[i];
1639  perm[i] = i;
1640  i = nexti;
1641  }
1642  while( perm[i] != v );
1643  consdata->quadvarterms[i] = quadterm;
1644  if( consdata->quadvarterms[i].eventdata != NULL )
1645  {
1646  consdata->quadvarterms[i].eventdata->varidx = -i-1;
1647  }
1648  perm[i] = i;
1649  }
1650  }
1651  consdata->quadvarssorted = TRUE;
1652 
1653  /* free temporary memory */
1654  SCIPfreeBufferArray(scip, &perm);
1655 
1656  return SCIP_OKAY;
1657 }
1658 
1659 /** returns the position of variable in the quadratic variable terms array of a constraint, or -1 if not found */
1660 static
1661 SCIP_RETCODE consdataFindQuadVarTerm(
1662  SCIP* scip, /**< SCIP data structure */
1663  SCIP_CONSDATA* consdata, /**< quadratic constraint data */
1664  SCIP_VAR* var, /**< variable to search for */
1665  int* pos /**< buffer where to store position of var in quadvarterms array, or -1 if not found */
1666  )
1667 {
1668  int left;
1669  int right;
1670  int cmpres;
1671 
1672  assert(consdata != NULL);
1673  assert(var != NULL);
1674  assert(pos != NULL);
1675 
1676  if( consdata->nquadvars == 0 )
1677  {
1678  *pos = -1;
1679  return SCIP_OKAY;
1680  }
1681 
1682  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
1683 
1684  left = 0;
1685  right = consdata->nquadvars - 1;
1686  while( left <= right )
1687  {
1688  int middle;
1689 
1690  middle = (left+right)/2;
1691  assert(0 <= middle && middle < consdata->nquadvars);
1692 
1693  cmpres = SCIPvarCompare(var, consdata->quadvarterms[middle].var);
1694 
1695  if( cmpres < 0 )
1696  right = middle - 1;
1697  else if( cmpres > 0 )
1698  left = middle + 1;
1699  else
1700  {
1701  *pos = middle;
1702  return SCIP_OKAY;
1703  }
1704  }
1705  assert(left == right+1);
1706 
1707  *pos = -1;
1708 
1709  return SCIP_OKAY;
1710 }
1711 
1712 /** index comparison method for bilinear terms: compares two index pairs of the bilinear term set in the quadratic constraint */
1713 static
1714 SCIP_DECL_SORTINDCOMP(bilinTermComp)
1715 { /*lint --e{715}*/
1716  SCIP_CONSDATA* consdata = (SCIP_CONSDATA*)dataptr;
1717  int var1cmp;
1718 
1719  assert(consdata != NULL);
1720  assert(0 <= ind1 && ind1 < consdata->nbilinterms);
1721  assert(0 <= ind2 && ind2 < consdata->nbilinterms);
1722 
1723  var1cmp = SCIPvarCompare(consdata->bilinterms[ind1].var1, consdata->bilinterms[ind2].var1);
1724  if( var1cmp != 0 )
1725  return var1cmp;
1726 
1727  return SCIPvarCompare(consdata->bilinterms[ind1].var2, consdata->bilinterms[ind2].var2);
1728 }
1729 
1730 #ifndef NDEBUG
1731 /** checks if all bilinear terms are sorted correctly */
1732 static
1733 SCIP_Bool consdataCheckBilinTermsSort(
1734  SCIP_CONSDATA* consdata
1735  )
1736 {
1737  int i;
1738 
1739  assert(consdata != NULL);
1740 
1741  /* nothing to check if the bilinear terms have not been sorted yet */
1742  if( !consdata->bilinsorted )
1743  return TRUE;
1744 
1745  for( i = 0; i < consdata->nbilinterms - 1; ++i )
1746  {
1747  if( bilinTermComp(consdata, i, i+1) > 0 )
1748  return FALSE;
1749  }
1750  return TRUE;
1751 }
1752 #endif
1753 
1754 /** sorting of bilinear terms */
1755 static
1756 SCIP_RETCODE consdataSortBilinTerms(
1757  SCIP* scip, /**< SCIP data structure */
1758  SCIP_CONSDATA* consdata /**< quadratic constraint data */
1759  )
1760 {
1761  int* perm;
1762  int* invperm;
1763  int i;
1764  int nexti;
1765  int v;
1766  SCIP_BILINTERM bilinterm;
1767 
1768  assert(scip != NULL);
1769  assert(consdata != NULL);
1770 
1771  if( consdata->bilinsorted )
1772  return SCIP_OKAY;
1773 
1774  if( consdata->nbilinterms == 0 )
1775  {
1776  consdata->bilinsorted = TRUE;
1777  return SCIP_OKAY;
1778  }
1779 
1780  /* get temporary memory to store the sorted permutation and the inverse permutation */
1781  SCIP_CALL( SCIPallocBufferArray(scip, &perm, consdata->nbilinterms) );
1782  SCIP_CALL( SCIPallocBufferArray(scip, &invperm, consdata->nbilinterms) );
1783 
1784  /* call bubble sort */
1785  SCIPsort(perm, bilinTermComp, (void*)consdata, consdata->nbilinterms);
1786 
1787  /* compute inverted permutation */
1788  for( v = 0; v < consdata->nbilinterms; ++v )
1789  {
1790  assert(0 <= perm[v] && perm[v] < consdata->nbilinterms);
1791  invperm[perm[v]] = v;
1792  }
1793 
1794  /* permute the bilinear terms according to the resulting permutation */
1795  for( v = 0; v < consdata->nbilinterms; ++v )
1796  {
1797  if( perm[v] != v )
1798  {
1799  bilinterm = consdata->bilinterms[v];
1800 
1801  i = v;
1802  do
1803  {
1804  assert(0 <= perm[i] && perm[i] < consdata->nbilinterms);
1805  assert(perm[i] != i);
1806  consdata->bilinterms[i] = consdata->bilinterms[perm[i]];
1807  nexti = perm[i];
1808  perm[i] = i;
1809  i = nexti;
1810  }
1811  while( perm[i] != v );
1812  consdata->bilinterms[i] = bilinterm;
1813  perm[i] = i;
1814  }
1815  }
1816 
1817  /* update the adjacency information in the quadratic variable terms */
1818  for( v = 0; v < consdata->nquadvars; ++v )
1819  for( i = 0; i < consdata->quadvarterms[v].nadjbilin; ++i )
1820  consdata->quadvarterms[v].adjbilin[i] = invperm[consdata->quadvarterms[v].adjbilin[i]];
1821 
1822  consdata->bilinsorted = TRUE;
1823  assert(consdataCheckBilinTermsSort(consdata));
1824 
1825  /* free temporary memory */
1826  SCIPfreeBufferArray(scip, &invperm);
1827  SCIPfreeBufferArray(scip, &perm);
1828 
1829  return SCIP_OKAY;
1830 }
1831 
1832 /** moves a linear variable from one position to another */
1833 static
1834 void consdataMoveLinearVar(
1835  SCIP_CONSDATA* consdata, /**< constraint data */
1836  int oldpos, /**< position of variable that shall be moved */
1837  int newpos /**< new position of variable */
1838  )
1839 {
1840  assert(consdata != NULL);
1841  assert(oldpos >= 0);
1842  assert(oldpos < consdata->nlinvars);
1843  assert(newpos >= 0);
1844  assert(newpos < consdata->linvarssize);
1845 
1846  if( newpos == oldpos )
1847  return;
1848 
1849  consdata->linvars [newpos] = consdata->linvars [oldpos];
1850  consdata->lincoefs[newpos] = consdata->lincoefs[oldpos];
1851 
1852  if( consdata->lineventdata != NULL )
1853  {
1854  assert(newpos >= consdata->nlinvars || consdata->lineventdata[newpos] == NULL);
1855 
1856  consdata->lineventdata[newpos] = consdata->lineventdata[oldpos];
1857  consdata->lineventdata[newpos]->varidx = newpos;
1858 
1859  consdata->lineventdata[oldpos] = NULL;
1860  }
1861 
1862  consdata->linvarssorted = FALSE;
1863 }
1864 
1865 /** moves a quadratic variable from one position to another */
1866 static
1867 void consdataMoveQuadVarTerm(
1868  SCIP_CONSDATA* consdata, /**< constraint data */
1869  int oldpos, /**< position of variable that shall be moved */
1870  int newpos /**< new position of variable */
1871  )
1872 {
1873  assert(consdata != NULL);
1874  assert(oldpos >= 0);
1875  assert(oldpos < consdata->nquadvars);
1876  assert(newpos >= 0);
1877  assert(newpos < consdata->quadvarssize);
1878 
1879  if( newpos == oldpos )
1880  return;
1881 
1882  assert(newpos >= consdata->nquadvars || consdata->quadvarterms[newpos].eventdata == NULL);
1883 
1884  consdata->quadvarterms[newpos] = consdata->quadvarterms[oldpos];
1885 
1886  if( consdata->quadvarterms[newpos].eventdata != NULL )
1887  {
1888  consdata->quadvarterms[newpos].eventdata->varidx = -newpos-1;
1889  consdata->quadvarterms[oldpos].eventdata = NULL;
1890  }
1891 
1892  consdata->quadvarssorted = FALSE;
1893 }
1894 
1895 /** adds linear coefficient in quadratic constraint */
1896 static
1897 SCIP_RETCODE addLinearCoef(
1898  SCIP* scip, /**< SCIP data structure */
1899  SCIP_CONS* cons, /**< quadratic constraint */
1900  SCIP_VAR* var, /**< variable of constraint entry */
1901  SCIP_Real coef /**< coefficient of constraint entry */
1902  )
1903 {
1904  SCIP_CONSDATA* consdata;
1905  SCIP_Bool transformed;
1906 
1907  assert(scip != NULL);
1908  assert(cons != NULL);
1909  assert(var != NULL);
1910 
1911  /* ignore coefficient if it is nearly zero */
1912  if( SCIPisZero(scip, coef) )
1913  return SCIP_OKAY;
1914 
1915  consdata = SCIPconsGetData(cons);
1916  assert(consdata != NULL);
1917 
1918  /* are we in the transformed problem? */
1919  transformed = SCIPconsIsTransformed(cons);
1920 
1921  /* always use transformed variables in transformed constraints */
1922  if( transformed )
1923  {
1924  SCIP_CALL( SCIPgetTransformedVar(scip, var, &var) );
1925  }
1926  assert(var != NULL);
1927  assert(transformed == SCIPvarIsTransformed(var));
1928 
1929  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars+1) );
1930  consdata->linvars [consdata->nlinvars] = var;
1931  consdata->lincoefs[consdata->nlinvars] = coef;
1932 
1933  ++consdata->nlinvars;
1934 
1935  /* catch variable events */
1936  if( SCIPconsIsEnabled(cons) )
1937  {
1938  SCIP_CONSHDLR* conshdlr;
1939  SCIP_CONSHDLRDATA* conshdlrdata;
1940 
1941  /* get event handler */
1942  conshdlr = SCIPconsGetHdlr(cons);
1943  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1944  assert(conshdlrdata != NULL);
1945  assert(conshdlrdata->eventhdlr != NULL);
1946 
1947  assert(consdata->lineventdata != NULL);
1948  consdata->lineventdata[consdata->nlinvars-1] = NULL;
1949 
1950  /* catch bound change events of variable */
1951  SCIP_CALL( catchLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, consdata->nlinvars-1) );
1952  }
1953 
1954  /* invalidate activity information */
1955  consdata->activity = SCIP_INVALID;
1956  consdata->minlinactivity = SCIP_INVALID;
1957  consdata->maxlinactivity = SCIP_INVALID;
1958  consdata->minlinactivityinf = -1;
1959  consdata->maxlinactivityinf = -1;
1960 
1961  /* invalidate nonlinear row */
1962  if( consdata->nlrow != NULL )
1963  {
1964  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1965  }
1966 
1967  /* install rounding locks for new variable */
1968  SCIP_CALL( lockLinearVariable(scip, cons, var, coef) );
1969 
1970  /* capture new variable */
1971  SCIP_CALL( SCIPcaptureVar(scip, var) );
1972 
1973  consdata->ispropagated = FALSE;
1974  consdata->ispresolved = FALSE;
1975  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
1976  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
1977  if( consdata->nlinvars == 1 )
1978  consdata->linvarssorted = TRUE;
1979  else
1980  consdata->linvarssorted = consdata->linvarssorted && (SCIPvarCompare(consdata->linvars[consdata->nlinvars-2], consdata->linvars[consdata->nlinvars-1]) == -1);
1981  /* always set too FALSE since the new linear variable should be checked if already existing as quad var term */
1982  consdata->linvarsmerged = FALSE;
1983 
1984  return SCIP_OKAY;
1985 }
1986 
1987 /** deletes linear coefficient at given position from quadratic constraint data */
1988 static
1989 SCIP_RETCODE delLinearCoefPos(
1990  SCIP* scip, /**< SCIP data structure */
1991  SCIP_CONS* cons, /**< quadratic constraint */
1992  int pos /**< position of coefficient to delete */
1993  )
1994 {
1995  SCIP_CONSDATA* consdata;
1996  SCIP_VAR* var;
1997  SCIP_Real coef;
1998 
1999  assert(scip != NULL);
2000  assert(cons != NULL);
2001 
2002  consdata = SCIPconsGetData(cons);
2003  assert(consdata != NULL);
2004  assert(0 <= pos && pos < consdata->nlinvars);
2005 
2006  var = consdata->linvars[pos];
2007  coef = consdata->lincoefs[pos];
2008  assert(var != NULL);
2009 
2010  /* remove rounding locks for deleted variable */
2011  SCIP_CALL( unlockLinearVariable(scip, cons, var, coef) );
2012 
2013  /* if we catch variable events, drop the events on the variable */
2014  if( consdata->lineventdata != NULL )
2015  {
2016  SCIP_CONSHDLR* conshdlr;
2017  SCIP_CONSHDLRDATA* conshdlrdata;
2018 
2019  /* get event handler */
2020  conshdlr = SCIPconsGetHdlr(cons);
2021  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2022  assert(conshdlrdata != NULL);
2023  assert(conshdlrdata->eventhdlr != NULL);
2024 
2025  /* drop bound change events of variable */
2026  SCIP_CALL( dropLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2027  }
2028 
2029  /* release variable */
2030  SCIP_CALL( SCIPreleaseVar(scip, &consdata->linvars[pos]) );
2031 
2032  /* move the last variable to the free slot */
2033  consdataMoveLinearVar(consdata, consdata->nlinvars-1, pos);
2034 
2035  --consdata->nlinvars;
2036 
2037  /* invalidate activity */
2038  consdata->activity = SCIP_INVALID;
2039  consdata->minlinactivity = SCIP_INVALID;
2040  consdata->maxlinactivity = SCIP_INVALID;
2041  consdata->minlinactivityinf = -1;
2042  consdata->maxlinactivityinf = -1;
2043 
2044  /* invalidate nonlinear row */
2045  if( consdata->nlrow != NULL )
2046  {
2047  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2048  }
2049 
2050  consdata->ispropagated = FALSE;
2051  consdata->ispresolved = FALSE;
2052 
2053  return SCIP_OKAY;
2054 }
2055 
2056 /** changes linear coefficient value at given position of quadratic constraint */
2057 static
2058 SCIP_RETCODE chgLinearCoefPos(
2059  SCIP* scip, /**< SCIP data structure */
2060  SCIP_CONS* cons, /**< quadratic constraint */
2061  int pos, /**< position of linear coefficient to change */
2062  SCIP_Real newcoef /**< new value of linear coefficient */
2063  )
2064 {
2065  SCIP_CONSHDLR* conshdlr;
2066  SCIP_CONSHDLRDATA* conshdlrdata;
2067  SCIP_CONSDATA* consdata;
2068  SCIP_VAR* var;
2069  SCIP_Real coef;
2070 
2071  assert(scip != NULL);
2072  assert(cons != NULL);
2073  assert(!SCIPconsIsLockedType(cons, SCIP_LOCKTYPE_CONFLICT));
2074  assert(!SCIPisZero(scip, newcoef));
2075 
2076  conshdlrdata = NULL;
2077 
2078  consdata = SCIPconsGetData(cons);
2079  assert(consdata != NULL);
2080  assert(0 <= pos);
2081  assert(pos < consdata->nlinvars);
2082  assert(!SCIPisZero(scip, newcoef));
2083 
2084  var = consdata->linvars[pos];
2085  coef = consdata->lincoefs[pos];
2086  assert(var != NULL);
2087  assert(SCIPconsIsTransformed(cons) == SCIPvarIsTransformed(var));
2088 
2089  /* invalidate activity */
2090  consdata->activity = SCIP_INVALID;
2091  consdata->minlinactivity = SCIP_INVALID;
2092  consdata->maxlinactivity = SCIP_INVALID;
2093  consdata->minlinactivityinf = -1;
2094  consdata->maxlinactivityinf = -1;
2095 
2096  /* invalidate nonlinear row */
2097  if( consdata->nlrow != NULL )
2098  {
2099  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2100  }
2101 
2102  /* if necessary, remove the rounding locks and event catching of the variable */
2103  if( newcoef * coef < 0.0 )
2104  {
2105  if( SCIPconsIsLocked(cons) )
2106  {
2107  assert(SCIPconsIsTransformed(cons));
2108 
2109  /* remove rounding locks for variable with old coefficient */
2110  SCIP_CALL( unlockLinearVariable(scip, cons, var, coef) );
2111  }
2112 
2113  if( consdata->lineventdata[pos] != NULL )
2114  {
2115  /* get event handler */
2116  conshdlr = SCIPconsGetHdlr(cons);
2117  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2118  assert(conshdlrdata != NULL);
2119  assert(conshdlrdata->eventhdlr != NULL);
2120 
2121  /* drop bound change events of variable */
2122  SCIP_CALL( dropLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2123  }
2124  }
2125 
2126  /* change the coefficient */
2127  consdata->lincoefs[pos] = newcoef;
2128 
2129  /* if necessary, install the rounding locks and event catching of the variable again */
2130  if( newcoef * coef < 0.0 )
2131  {
2132  if( SCIPconsIsLocked(cons) )
2133  {
2134  /* install rounding locks for variable with new coefficient */
2135  SCIP_CALL( lockLinearVariable(scip, cons, var, newcoef) );
2136  }
2137 
2138  if( conshdlrdata != NULL )
2139  {
2140  assert(SCIPconsIsEnabled(cons));
2141 
2142  /* catch bound change events of variable */
2143  SCIP_CALL( catchLinearVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2144  }
2145  }
2146 
2147  consdata->ispropagated = FALSE;
2148  consdata->ispresolved = FALSE;
2149 
2150  return SCIP_OKAY;
2151 }
2152 
2153 /** adds quadratic variable term to quadratic constraint */
2154 static
2155 SCIP_RETCODE addQuadVarTerm(
2156  SCIP* scip, /**< SCIP data structure */
2157  SCIP_CONS* cons, /**< quadratic constraint */
2158  SCIP_VAR* var, /**< variable to add */
2159  SCIP_Real lincoef, /**< linear coefficient of variable */
2160  SCIP_Real sqrcoef /**< square coefficient of variable */
2161  )
2162 {
2163  SCIP_CONSDATA* consdata;
2164  SCIP_Bool transformed;
2165  SCIP_QUADVARTERM* quadvarterm;
2166 
2167  assert(scip != NULL);
2168  assert(cons != NULL);
2169  assert(var != NULL);
2170 
2171  consdata = SCIPconsGetData(cons);
2172  assert(consdata != NULL);
2173 
2174  /* are we in the transformed problem? */
2175  transformed = SCIPconsIsTransformed(cons);
2176 
2177  /* always use transformed variables in transformed constraints */
2178  if( transformed )
2179  {
2180  SCIP_CALL( SCIPgetTransformedVar(scip, var, &var) );
2181  }
2182  assert(var != NULL);
2183  assert(transformed == SCIPvarIsTransformed(var));
2184 
2185  SCIP_CALL( consdataEnsureQuadVarTermsSize(scip, consdata, consdata->nquadvars+1) );
2186 
2187  quadvarterm = &consdata->quadvarterms[consdata->nquadvars];
2188  quadvarterm->var = var;
2189  quadvarterm->lincoef = lincoef;
2190  quadvarterm->sqrcoef = sqrcoef;
2191  quadvarterm->adjbilinsize = 0;
2192  quadvarterm->nadjbilin = 0;
2193  quadvarterm->adjbilin = NULL;
2194  quadvarterm->eventdata = NULL;
2195 
2196  ++consdata->nquadvars;
2197 
2198  /* capture variable */
2199  SCIP_CALL( SCIPcaptureVar(scip, var) );
2200 
2201  /* catch variable events, if we do so */
2202  if( SCIPconsIsEnabled(cons) )
2203  {
2204  SCIP_CONSHDLR* conshdlr;
2205  SCIP_CONSHDLRDATA* conshdlrdata;
2206 
2207  /* get event handler */
2208  conshdlr = SCIPconsGetHdlr(cons);
2209  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2210  assert(conshdlrdata != NULL);
2211  assert(conshdlrdata->eventhdlr != NULL);
2212 
2213  /* catch bound change events of variable */
2214  SCIP_CALL( catchQuadVarEvents(scip, conshdlrdata->eventhdlr, cons, consdata->nquadvars-1) );
2215  }
2216 
2217  /* invalidate activity information */
2218  consdata->activity = SCIP_INVALID;
2219  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2220 
2221  /* invalidate nonlinear row */
2222  if( consdata->nlrow != NULL )
2223  {
2224  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2225  }
2226 
2227  /* install rounding locks for new variable */
2228  SCIP_CALL( lockQuadraticVariable(scip, cons, var) );
2229 
2230  consdata->ispropagated = FALSE;
2231  consdata->ispresolved = FALSE;
2232  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
2233  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
2234  if( consdata->nquadvars == 1 )
2235  consdata->quadvarssorted = TRUE;
2236  else
2237  consdata->quadvarssorted = consdata->quadvarssorted &&
2238  (SCIPvarCompare(consdata->quadvarterms[consdata->nquadvars-2].var, consdata->quadvarterms[consdata->nquadvars-1].var) == -1);
2239  /* also set to FALSE if nquadvars == 1, since the new variable should be checked for linearity and other stuff in mergeAndClean ... */
2240  consdata->quadvarsmerged = FALSE;
2241 
2242  consdata->iscurvchecked = FALSE;
2243 
2244  return SCIP_OKAY;
2245 }
2246 
2247 /** deletes quadratic variable term at given position from quadratic constraint data */
2248 static
2249 SCIP_RETCODE delQuadVarTermPos(
2250  SCIP* scip, /**< SCIP data structure */
2251  SCIP_CONS* cons, /**< quadratic constraint */
2252  int pos /**< position of term to delete */
2253  )
2254 {
2255  SCIP_CONSDATA* consdata;
2256  SCIP_VAR* var;
2257 
2258  assert(scip != NULL);
2259  assert(cons != NULL);
2260 
2261  consdata = SCIPconsGetData(cons);
2262  assert(consdata != NULL);
2263  assert(0 <= pos && pos < consdata->nquadvars);
2264 
2265  var = consdata->quadvarterms[pos].var;
2266  assert(var != NULL);
2267  assert(consdata->quadvarterms[pos].nadjbilin == 0);
2268 
2269  /* remove rounding locks for deleted variable */
2270  SCIP_CALL( unlockQuadraticVariable(scip, cons, var) );
2271 
2272  /* if we catch variable events, drop the events on the variable */
2273  if( consdata->quadvarterms[pos].eventdata != NULL )
2274  {
2275  SCIP_CONSHDLR* conshdlr;
2276  SCIP_CONSHDLRDATA* conshdlrdata;
2277 
2278  /* get event handler */
2279  conshdlr = SCIPconsGetHdlr(cons);
2280  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2281  assert(conshdlrdata != NULL);
2282  assert(conshdlrdata->eventhdlr != NULL);
2283 
2284  /* drop bound change events of variable */
2285  SCIP_CALL( dropQuadVarEvents(scip, conshdlrdata->eventhdlr, cons, pos) );
2286  }
2287 
2288  /* release variable */
2289  SCIP_CALL( SCIPreleaseVar(scip, &consdata->quadvarterms[pos].var) );
2290 
2291  /* free adjacency array */
2292  SCIPfreeBlockMemoryArrayNull(scip, &consdata->quadvarterms[pos].adjbilin, consdata->quadvarterms[pos].adjbilinsize);
2293 
2294  /* move the last variable term to the free slot */
2295  consdataMoveQuadVarTerm(consdata, consdata->nquadvars-1, pos);
2296 
2297  --consdata->nquadvars;
2298 
2299  /* invalidate activity */
2300  consdata->activity = SCIP_INVALID;
2301  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2302 
2303  /* invalidate nonlinear row */
2304  if( consdata->nlrow != NULL )
2305  {
2306  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2307  }
2308 
2309  consdata->ispropagated = FALSE;
2310  consdata->ispresolved = FALSE;
2311  consdata->iscurvchecked = FALSE;
2312 
2313  return SCIP_OKAY;
2314 }
2315 
2316 /** replace variable in quadratic variable term at given position of quadratic constraint data
2317  *
2318  * Allows to replace x by coef*y+offset, thereby maintaining linear and square coefficients and bilinear terms.
2319  */
2320 static
2321 SCIP_RETCODE replaceQuadVarTermPos(
2322  SCIP* scip, /**< SCIP data structure */
2323  SCIP_CONS* cons, /**< quadratic constraint */
2324  int pos, /**< position of term to replace */
2325  SCIP_VAR* var, /**< new variable */
2326  SCIP_Real coef, /**< linear coefficient of new variable */
2327  SCIP_Real offset /**< offset of new variable */
2328  )
2329 {
2330  SCIP_CONSDATA* consdata;
2331  SCIP_QUADVARTERM* quadvarterm;
2332  SCIP_EVENTHDLR* eventhdlr;
2333  SCIP_BILINTERM* bilinterm;
2334  SCIP_Real constant;
2335 
2336  int i;
2337  SCIP_VAR* var2;
2338 
2339  consdata = SCIPconsGetData(cons);
2340  assert(consdata != NULL);
2341  assert(pos >= 0);
2342  assert(pos < consdata->nquadvars);
2343 
2344  quadvarterm = &consdata->quadvarterms[pos];
2345 
2346  /* remove rounding locks for old variable */
2347  SCIP_CALL( unlockQuadraticVariable(scip, cons, quadvarterm->var) );
2348 
2349  /* if we catch variable events, drop the events on the old variable */
2350  if( quadvarterm->eventdata != NULL )
2351  {
2352  SCIP_CONSHDLR* conshdlr;
2353  SCIP_CONSHDLRDATA* conshdlrdata;
2354 
2355  /* get event handler */
2356  conshdlr = SCIPconsGetHdlr(cons);
2357  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2358  assert(conshdlrdata != NULL);
2359  assert(conshdlrdata->eventhdlr != NULL);
2360 
2361  eventhdlr = conshdlrdata->eventhdlr;
2362 
2363  /* drop bound change events of variable */
2364  SCIP_CALL( dropQuadVarEvents(scip, eventhdlr, cons, pos) );
2365  }
2366  else
2367  {
2368  eventhdlr = NULL;
2369  }
2370 
2371  /* compute constant and put into lhs/rhs */
2372  constant = quadvarterm->lincoef * offset + quadvarterm->sqrcoef * offset * offset;
2373  if( constant != 0.0 )
2374  {
2375  /* maintain constant part */
2376  if( !SCIPisInfinity(scip, -consdata->lhs) )
2377  consdata->lhs -= constant;
2378  if( !SCIPisInfinity(scip, consdata->rhs) )
2379  consdata->rhs -= constant;
2380  }
2381 
2382  /* update linear and square coefficient */
2383  quadvarterm->lincoef *= coef;
2384  quadvarterm->lincoef += 2.0 * quadvarterm->sqrcoef * coef * offset;
2385  quadvarterm->sqrcoef *= coef * coef;
2386 
2387  /* update bilinear terms */
2388  for( i = 0; i < quadvarterm->nadjbilin; ++i )
2389  {
2390  bilinterm = &consdata->bilinterms[quadvarterm->adjbilin[i]];
2391 
2392  if( bilinterm->var1 == quadvarterm->var )
2393  {
2394  bilinterm->var1 = var;
2395  var2 = bilinterm->var2;
2396  }
2397  else
2398  {
2399  assert(bilinterm->var2 == quadvarterm->var);
2400  bilinterm->var2 = var;
2401  var2 = bilinterm->var1;
2402  }
2403 
2404  if( var == var2 )
2405  {
2406  /* looks like we actually have a square term here */
2407  quadvarterm->lincoef += bilinterm->coef * offset;
2408  quadvarterm->sqrcoef += bilinterm->coef * coef;
2409  /* deleting bilinear terms is expensive, since it requires updating adjacency information
2410  * thus, for now we just set the coefficient to 0.0 and clear in later when the bilinear terms are merged */
2411  bilinterm->coef = 0.0;
2412  continue;
2413  }
2414 
2415  /* swap var1 and var2 if they are in wrong order */
2416  if( SCIPvarCompare(bilinterm->var1, bilinterm->var2) > 0 )
2417  {
2418  SCIP_VAR* tmp;
2419  tmp = bilinterm->var1;
2420  bilinterm->var1 = bilinterm->var2;
2421  bilinterm->var2 = tmp;
2422  }
2423  assert(SCIPvarCompare(bilinterm->var1, bilinterm->var2) == -1);
2424 
2425  if( offset != 0.0 )
2426  {
2427  /* need to find var2 and add offset*bilinterm->coef to linear coefficient */
2428  int var2pos;
2429 
2430  var2pos = 0;
2431  while( consdata->quadvarterms[var2pos].var != var2 )
2432  {
2433  ++var2pos;
2434  assert(var2pos < consdata->nquadvars);
2435  }
2436 
2437  consdata->quadvarterms[var2pos].lincoef += bilinterm->coef * offset;
2438  }
2439 
2440  bilinterm->coef *= coef;
2441  }
2442 
2443  /* release old variable */
2444  SCIP_CALL( SCIPreleaseVar(scip, &quadvarterm->var) );
2445 
2446  /* set new variable */
2447  quadvarterm->var = var;
2448 
2449  /* capture new variable */
2450  SCIP_CALL( SCIPcaptureVar(scip, quadvarterm->var) );
2451 
2452  /* catch variable events, if we do so */
2453  if( eventhdlr != NULL )
2454  {
2455  assert(SCIPconsIsEnabled(cons));
2456 
2457  /* catch bound change events of variable */
2458  SCIP_CALL( catchQuadVarEvents(scip, eventhdlr, cons, pos) );
2459  }
2460 
2461  /* invalidate activity information */
2462  consdata->activity = SCIP_INVALID;
2463  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2464 
2465  /* invalidate nonlinear row */
2466  if( consdata->nlrow != NULL )
2467  {
2468  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2469  }
2470 
2471  /* install rounding locks for new variable */
2472  SCIP_CALL( lockQuadraticVariable(scip, cons, var) );
2473 
2474  consdata->isremovedfixings = consdata->isremovedfixings && SCIPvarIsActive(var)
2475  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
2476  consdata->quadvarssorted = (consdata->nquadvars == 1);
2477  consdata->quadvarsmerged = FALSE;
2478  consdata->bilinsorted &= (quadvarterm->nadjbilin == 0); /*lint !e514*/
2479  consdata->bilinmerged &= (quadvarterm->nadjbilin == 0); /*lint !e514*/
2480 
2481  consdata->ispropagated = FALSE;
2482  consdata->ispresolved = FALSE;
2483  consdata->iscurvchecked = FALSE;
2484 
2485  return SCIP_OKAY;
2486 }
2487 
2488 /** adds a bilinear term to quadratic constraint */
2489 static
2490 SCIP_RETCODE addBilinearTerm(
2491  SCIP* scip, /**< SCIP data structure */
2492  SCIP_CONS* cons, /**< quadratic constraint */
2493  int var1pos, /**< position of first variable in quadratic variables array */
2494  int var2pos, /**< position of second variable in quadratic variables array */
2495  SCIP_Real coef /**< coefficient of bilinear term */
2496  )
2497 {
2498  SCIP_CONSDATA* consdata;
2499  SCIP_BILINTERM* bilinterm;
2500 
2501  assert(scip != NULL);
2502  assert(cons != NULL);
2503 
2504  if( var1pos == var2pos )
2505  {
2506  SCIPerrorMessage("tried to add bilinear term where both variables are the same\n");
2507  return SCIP_INVALIDDATA;
2508  }
2509 
2510  consdata = SCIPconsGetData(cons);
2511  assert(consdata != NULL);
2512 
2513  /* check if the bilinear terms are sorted (disabled for big constraints as becoming expensive) */
2514  assert(consdata->nbilinterms > 10 || consdataCheckBilinTermsSort(consdata));
2515 
2516  assert(var1pos >= 0);
2517  assert(var1pos < consdata->nquadvars);
2518  assert(var2pos >= 0);
2519  assert(var2pos < consdata->nquadvars);
2520 
2521  SCIP_CALL( consdataEnsureBilinSize(scip, consdata, consdata->nbilinterms + 1) );
2522 
2523  bilinterm = &consdata->bilinterms[consdata->nbilinterms];
2524  if( SCIPvarCompare(consdata->quadvarterms[var1pos].var, consdata->quadvarterms[var2pos].var) < 0 )
2525  {
2526  bilinterm->var1 = consdata->quadvarterms[var1pos].var;
2527  bilinterm->var2 = consdata->quadvarterms[var2pos].var;
2528  }
2529  else
2530  {
2531  bilinterm->var1 = consdata->quadvarterms[var2pos].var;
2532  bilinterm->var2 = consdata->quadvarterms[var1pos].var;
2533  }
2534  bilinterm->coef = coef;
2535 
2536  if( bilinterm->var1 == bilinterm->var2 )
2537  {
2538  SCIPerrorMessage("tried to add bilinear term where both variables are the same, but appear at different positions in quadvarterms array\n");
2539  return SCIP_INVALIDDATA;
2540  }
2541  assert(SCIPvarCompare(bilinterm->var1, bilinterm->var2) == -1);
2542 
2543  SCIP_CALL( consdataEnsureAdjBilinSize(scip, &consdata->quadvarterms[var1pos], consdata->quadvarterms[var1pos].nadjbilin + 1) );
2544  SCIP_CALL( consdataEnsureAdjBilinSize(scip, &consdata->quadvarterms[var2pos], consdata->quadvarterms[var2pos].nadjbilin + 1) );
2545 
2546  consdata->quadvarterms[var1pos].adjbilin[consdata->quadvarterms[var1pos].nadjbilin] = consdata->nbilinterms;
2547  consdata->quadvarterms[var2pos].adjbilin[consdata->quadvarterms[var2pos].nadjbilin] = consdata->nbilinterms;
2548  ++consdata->quadvarterms[var1pos].nadjbilin;
2549  ++consdata->quadvarterms[var2pos].nadjbilin;
2550 
2551  ++consdata->nbilinterms;
2552 
2553  /* invalidate activity information */
2554  consdata->activity = SCIP_INVALID;
2555  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2556 
2557  /* invalidate nonlinear row */
2558  if( consdata->nlrow != NULL )
2559  {
2560  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2561  }
2562 
2563  consdata->ispropagated = FALSE;
2564  consdata->ispresolved = FALSE;
2565  if( consdata->nbilinterms == 1 )
2566  {
2567  consdata->bilinsorted = TRUE;
2568 
2569  /* we have to take care of the bilinear term in mergeAndCleanBilinearTerms() if the coefficient is zero */
2570  consdata->bilinmerged = !SCIPisZero(scip, consdata->bilinterms[0].coef);
2571  }
2572  else
2573  {
2574  consdata->bilinsorted = consdata->bilinsorted
2575  && (bilinTermComp(consdata, consdata->nbilinterms-2, consdata->nbilinterms-1) <= 0);
2576  consdata->bilinmerged = FALSE;
2577  }
2578 
2579  consdata->iscurvchecked = FALSE;
2580 
2581  /* check if the bilinear terms are sorted (disabled as expensive if big constraint) */
2582  assert(consdata->nbilinterms > 10 || consdataCheckBilinTermsSort(consdata));
2583 
2584  return SCIP_OKAY;
2585 }
2586 
2587 /** removes a set of bilinear terms and updates adjacency information in quad var terms
2588  *
2589  * Note: this function sorts the given array termposs.
2590  */
2591 static
2592 SCIP_RETCODE removeBilinearTermsPos(
2593  SCIP* scip, /**< SCIP data structure */
2594  SCIP_CONS* cons, /**< quadratic constraint */
2595  int nterms, /**< number of terms to delete */
2596  int* termposs /**< indices of terms to delete */
2597  )
2598 {
2599  SCIP_CONSDATA* consdata;
2600  int* newpos;
2601  int i;
2602  int j;
2603  int offset;
2604 
2605  assert(scip != NULL);
2606  assert(cons != NULL);
2607  assert(nterms == 0 || termposs != NULL);
2608 
2609  if( nterms == 0 )
2610  return SCIP_OKAY;
2611 
2612  consdata = SCIPconsGetData(cons);
2613  assert(consdata != NULL);
2614 
2615  SCIPsortInt(termposs, nterms);
2616 
2617  SCIP_CALL( SCIPallocBufferArray(scip, &newpos, consdata->nbilinterms) );
2618 
2619  i = 0;
2620  offset = 0;
2621  for( j = 0; j < consdata->nbilinterms; ++j )
2622  {
2623  /* if j'th term is deleted, increase offset and continue */
2624  if( i < nterms && j == termposs[i] )
2625  {
2626  ++offset;
2627  ++i;
2628  newpos[j] = -1;
2629  continue;
2630  }
2631 
2632  /* otherwise, move it forward and remember new position */
2633  if( offset > 0 )
2634  consdata->bilinterms[j-offset] = consdata->bilinterms[j];
2635  newpos[j] = j - offset;
2636  }
2637  assert(offset == nterms);
2638 
2639  /* update adjacency and activity information in quad var terms */
2640  for( i = 0; i < consdata->nquadvars; ++i )
2641  {
2642  offset = 0;
2643  for( j = 0; j < consdata->quadvarterms[i].nadjbilin; ++j )
2644  {
2645  assert(consdata->quadvarterms[i].adjbilin[j] < consdata->nbilinterms);
2646  if( newpos[consdata->quadvarterms[i].adjbilin[j]] == -1 )
2647  {
2648  /* corresponding bilinear term was deleted, thus increase offset */
2649  ++offset;
2650  }
2651  else
2652  {
2653  /* update index of j'th bilinear term and store at position j-offset */
2654  consdata->quadvarterms[i].adjbilin[j-offset] = newpos[consdata->quadvarterms[i].adjbilin[j]];
2655  }
2656  }
2657  consdata->quadvarterms[i].nadjbilin -= offset;
2658  /* some bilinear term was removed, so invalidate activity bounds */
2659  }
2660 
2661  consdata->nbilinterms -= nterms;
2662 
2663  SCIPfreeBufferArray(scip, &newpos);
2664 
2665  /* some quad vars may be linear now */
2666  consdata->quadvarsmerged = FALSE;
2667 
2668  consdata->ispropagated = FALSE;
2669  consdata->ispresolved = FALSE;
2670  consdata->iscurvchecked = FALSE;
2671 
2672  /* invalidate activity */
2673  consdata->activity = SCIP_INVALID;
2674  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2675 
2676  /* invalidate nonlinear row */
2677  if( consdata->nlrow != NULL )
2678  {
2679  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2680  }
2681 
2682  return SCIP_OKAY;
2683 }
2684 
2685 /** changes side of constraint and allow to change between finite and infinite
2686  *
2687  * takes care of updating events and locks of linear variables
2688  */
2689 static
2690 SCIP_RETCODE chgSideQuadratic(
2691  SCIP* scip, /**< SCIP data structure */
2692  SCIP_CONS* cons, /**< constraint */
2693  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
2694  SCIP_SIDETYPE side, /**< which side to change */
2695  SCIP_Real sideval /**< new value for side */
2696  )
2697 {
2698  SCIP_CONSDATA* consdata;
2699  int i;
2700 
2701  assert(scip != NULL);
2702  assert(cons != NULL);
2703  assert(!SCIPisInfinity(scip, side == SCIP_SIDETYPE_LEFT ? sideval : -sideval));
2704 
2705  consdata = SCIPconsGetData(cons);
2706  assert(consdata != NULL);
2707 
2708  /* if remaining finite or remaining infinite, then can just update the value */
2709  if( side == SCIP_SIDETYPE_LEFT )
2710  {
2711  if( SCIPisInfinity(scip, -consdata->lhs) == SCIPisInfinity(scip, -sideval) )
2712  {
2713  consdata->lhs = sideval;
2714  return SCIP_OKAY;
2715  }
2716  }
2717  else
2718  {
2719  if( SCIPisInfinity(scip, consdata->rhs) == SCIPisInfinity(scip, sideval) )
2720  {
2721  consdata->rhs = sideval;
2722  return SCIP_OKAY;
2723  }
2724  }
2725 
2726  /* catched boundchange events and locks for linear variables depends on whether side is finite, so first drop all */
2727  for( i = 0; i < consdata->nlinvars; ++i )
2728  {
2729  if( consdata->lineventdata != NULL && consdata->lineventdata[i] != NULL )
2730  {
2731  assert(SCIPconsIsEnabled(cons));
2732 
2733  SCIP_CALL( dropLinearVarEvents(scip, eventhdlr, cons, i) );
2734  }
2735 
2736  if( SCIPconsIsLocked(cons) )
2737  {
2738  assert(SCIPconsIsTransformed(cons));
2739 
2740  /* remove rounding locks for variable with old side */
2741  SCIP_CALL( unlockLinearVariable(scip, cons, consdata->linvars[i], consdata->lincoefs[i]) );
2742  }
2743  }
2744 
2745  if( side == SCIP_SIDETYPE_LEFT )
2746  consdata->lhs = sideval;
2747  else
2748  consdata->rhs = sideval;
2749 
2750  /* catch boundchange events and locks on variables again */
2751  for( i = 0; i < consdata->nlinvars; ++i )
2752  {
2753  if( consdata->lineventdata != NULL )
2754  {
2755  SCIP_CALL( catchLinearVarEvents(scip, eventhdlr, cons, i) );
2756  }
2757 
2758  if( SCIPconsIsLocked(cons) )
2759  {
2760  /* add rounding locks for variable with new side */
2761  SCIP_CALL( lockLinearVariable(scip, cons, consdata->linvars[i], consdata->lincoefs[i]) );
2762  }
2763  }
2764 
2765  return SCIP_OKAY;
2766 }
2767 
2768 /** merges quad var terms that correspond to the same variable and does additional cleanup
2769  *
2770  * If a quadratic variable terms is actually linear, makes a linear term out of it
2771  * also replaces squares of binary variables by the binary variables, i.e., adds sqrcoef to lincoef.
2772  */
2773 static
2774 SCIP_RETCODE mergeAndCleanQuadVarTerms(
2775  SCIP* scip, /**< SCIP data structure */
2776  SCIP_CONS* cons /**< quadratic constraint */
2777  )
2778 {
2779  SCIP_QUADVARTERM* quadvarterm;
2780  SCIP_CONSDATA* consdata;
2781  int i;
2782  int j;
2783 
2784  assert(scip != NULL);
2785  assert(cons != NULL);
2786 
2787  consdata = SCIPconsGetData(cons);
2788 
2789  if( consdata->quadvarsmerged )
2790  return SCIP_OKAY;
2791 
2792  if( consdata->nquadvars == 0 )
2793  {
2794  consdata->quadvarsmerged = TRUE;
2795  return SCIP_OKAY;
2796  }
2797 
2798  i = 0;
2799  while( i < consdata->nquadvars )
2800  {
2801  /* make sure quad var terms are sorted (do this in every round, since we may move variables around) */
2802  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
2803 
2804  quadvarterm = &consdata->quadvarterms[i];
2805 
2806  for( j = i+1; j < consdata->nquadvars && consdata->quadvarterms[j].var == quadvarterm->var; ++j )
2807  {
2808  /* add quad var term j to current term i */
2809  quadvarterm->lincoef += consdata->quadvarterms[j].lincoef;
2810  quadvarterm->sqrcoef += consdata->quadvarterms[j].sqrcoef;
2811  if( consdata->quadvarterms[j].nadjbilin > 0 )
2812  {
2813  /* move adjacency information from j to i */
2814  SCIP_CALL( consdataEnsureAdjBilinSize(scip, quadvarterm, quadvarterm->nadjbilin + consdata->quadvarterms[j].nadjbilin) );
2815  BMScopyMemoryArray(&quadvarterm->adjbilin[quadvarterm->nadjbilin], consdata->quadvarterms[j].adjbilin, consdata->quadvarterms[j].nadjbilin); /*lint !e866*/
2816  quadvarterm->nadjbilin += consdata->quadvarterms[j].nadjbilin;
2817  consdata->quadvarterms[j].nadjbilin = 0;
2818  }
2819  consdata->quadvarterms[j].lincoef = 0.0;
2820  consdata->quadvarterms[j].sqrcoef = 0.0;
2821  /* mark that activity information in quadvarterm is not up to date anymore */
2822  }
2823 
2824  /* remove quad var terms i+1..j-1 backwards */
2825  for( j = j-1; j > i; --j )
2826  {
2827  SCIP_CALL( delQuadVarTermPos(scip, cons, j) );
2828  }
2829 
2830  /* for binary variables, x^2 = x
2831  * however, we may destroy convexity of a quadratic term that involves also bilinear terms
2832  * thus, we do this step only if the variable does not appear in any bilinear term */
2833  if( quadvarterm->sqrcoef != 0.0 && SCIPvarIsBinary(quadvarterm->var) && quadvarterm->nadjbilin == 0 )
2834  {
2835  SCIPdebugMsg(scip, "replace square of binary variable by itself: <%s>^2 --> <%s>\n", SCIPvarGetName(quadvarterm->var), SCIPvarGetName(quadvarterm->var));
2836  quadvarterm->lincoef += quadvarterm->sqrcoef;
2837  quadvarterm->sqrcoef = 0.0;
2838 
2839  /* invalidate nonlinear row */
2840  if( consdata->nlrow != NULL )
2841  {
2842  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
2843  }
2844  }
2845 
2846  /* if its 0.0 or linear, get rid of it */
2847  if( SCIPisZero(scip, quadvarterm->sqrcoef) && quadvarterm->nadjbilin == 0 )
2848  {
2849  if( !SCIPisZero(scip, quadvarterm->lincoef) )
2850  {
2851  /* seem to be a linear term now, thus add as linear term */
2852  SCIP_CALL( addLinearCoef(scip, cons, quadvarterm->var, quadvarterm->lincoef) );
2853  }
2854  /* remove term at pos i */
2855  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
2856  }
2857  else
2858  {
2859  ++i;
2860  }
2861  }
2862 
2863  consdata->quadvarsmerged = TRUE;
2864  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2865 
2866  return SCIP_OKAY;
2867 }
2868 
2869 /** merges entries with same linear variable into one entry and cleans up entries with coefficient 0.0 */
2870 static
2871 SCIP_RETCODE mergeAndCleanLinearVars(
2872  SCIP* scip, /**< SCIP data structure */
2873  SCIP_CONS* cons /**< quadratic constraint */
2874  )
2875 {
2876  SCIP_CONSDATA* consdata;
2877  SCIP_Real newcoef;
2878  int i;
2879  int j;
2880  int qvarpos;
2881 
2882  assert(scip != NULL);
2883  assert(cons != NULL);
2884 
2885  consdata = SCIPconsGetData(cons);
2886 
2887  if( consdata->linvarsmerged )
2888  return SCIP_OKAY;
2889 
2890  if( consdata->nlinvars == 0 )
2891  {
2892  consdata->linvarsmerged = TRUE;
2893  return SCIP_OKAY;
2894  }
2895 
2896  i = 0;
2897  while( i < consdata->nlinvars )
2898  {
2899  /* make sure linear variables are sorted (do this in every round, since we may move variables around) */
2900  consdataSortLinearVars(consdata);
2901 
2902  /* sum up coefficients that correspond to variable i */
2903  newcoef = consdata->lincoefs[i];
2904  for( j = i+1; j < consdata->nlinvars && consdata->linvars[i] == consdata->linvars[j]; ++j )
2905  newcoef += consdata->lincoefs[j];
2906  /* delete the additional variables in backward order */
2907  for( j = j-1; j > i; --j )
2908  {
2909  SCIP_CALL( delLinearCoefPos(scip, cons, j) );
2910  }
2911 
2912  /* check if there is already a quadratic variable term with this variable */
2913  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, consdata->linvars[i], &qvarpos) );
2914  if( qvarpos >= 0)
2915  {
2916  /* add newcoef to linear coefficient of quadratic variable and mark linear variable as to delete */
2917  assert(qvarpos < consdata->nquadvars);
2918  assert(consdata->quadvarterms[qvarpos].var == consdata->linvars[i]);
2919  consdata->quadvarterms[qvarpos].lincoef += newcoef;
2920  newcoef = 0.0;
2921  SCIPintervalSetEmpty(&consdata->quadactivitybounds);
2922  }
2923 
2924  /* delete also entry at position i, if it became zero (or was zero before) */
2925  if( SCIPisZero(scip, newcoef) )
2926  {
2927  SCIP_CALL( delLinearCoefPos(scip, cons, i) );
2928  }
2929  else
2930  {
2931  SCIP_CALL( chgLinearCoefPos(scip, cons, i, newcoef) );
2932  ++i;
2933  }
2934  }
2935 
2936  consdata->linvarsmerged = TRUE;
2937 
2938  return SCIP_OKAY;
2939 }
2940 
2941 /** merges bilinear terms with same variables into a single term, removes bilinear terms with coefficient 0.0 */
2942 static
2943 SCIP_RETCODE mergeAndCleanBilinearTerms(
2944  SCIP* scip, /**< SCIP data structure */
2945  SCIP_CONS* cons /**< quadratic constraint */
2946  )
2947 {
2948  SCIP_CONSDATA* consdata;
2949  SCIP_BILINTERM* bilinterm;
2950  int i;
2951  int j;
2952  int* todelete;
2953  int ntodelete;
2954 
2955  assert(scip != NULL);
2956  assert(cons != NULL);
2957 
2958  consdata = SCIPconsGetData(cons);
2959 
2960  /* check if the bilinear terms are sorted */
2961  assert(consdataCheckBilinTermsSort(consdata));
2962 
2963  if( consdata->bilinmerged )
2964  return SCIP_OKAY;
2965 
2966  if( consdata->nbilinterms == 0 )
2967  {
2968  consdata->bilinmerged = TRUE;
2969  return SCIP_OKAY;
2970  }
2971 
2972  /* alloc memory for array of terms that need to be deleted finally */
2973  ntodelete = 0;
2974  SCIP_CALL( SCIPallocBufferArray(scip, &todelete, consdata->nbilinterms) );
2975 
2976  /* make sure bilinear terms are sorted */
2977  SCIP_CALL( consdataSortBilinTerms(scip, consdata) );
2978 
2979  i = 0;
2980  while( i < consdata->nbilinterms )
2981  {
2982  bilinterm = &consdata->bilinterms[i];
2983 
2984  /* sum up coefficients that correspond to same variables as term i */
2985  for( j = i+1; j < consdata->nbilinterms && bilinterm->var1 == consdata->bilinterms[j].var1 && bilinterm->var2 == consdata->bilinterms[j].var2; ++j )
2986  {
2987  bilinterm->coef += consdata->bilinterms[j].coef;
2988  todelete[ntodelete++] = j;
2989  }
2990 
2991  /* delete also entry at position i, if it became zero (or was zero before) */
2992  if( SCIPisZero(scip, bilinterm->coef) )
2993  {
2994  todelete[ntodelete++] = i;
2995  }
2996 
2997  /* continue with term after the current series */
2998  i = j;
2999  }
3000 
3001  /* delete bilinear terms */
3002  SCIP_CALL( removeBilinearTermsPos(scip, cons, ntodelete, todelete) );
3003 
3004  SCIPfreeBufferArray(scip, &todelete);
3005 
3006  consdata->bilinmerged = TRUE;
3007 
3008  /* check if the bilinear terms are sorted */
3009  assert(consdataCheckBilinTermsSort(consdata));
3010 
3011  return SCIP_OKAY;
3012 }
3013 
3014 /** removes fixes (or aggregated) variables from a quadratic constraint */
3015 static
3016 SCIP_RETCODE removeFixedVariables(
3017  SCIP* scip, /**< SCIP data structure */
3018  SCIP_EVENTHDLR* eventhdlr, /**< event handler for variable bound changes */
3019  SCIP_CONS* cons /**< quadratic constraint */
3020  )
3021 {
3022  SCIP_CONSDATA* consdata;
3023  SCIP_BILINTERM* bilinterm;
3024  SCIP_Real bilincoef;
3025  SCIP_Real coef;
3026  SCIP_Real offset;
3027  SCIP_VAR* var;
3028  SCIP_VAR* var2;
3029  int var2pos;
3030  int i;
3031  int j;
3032  int k;
3033 
3034  SCIP_Bool have_change;
3035 
3036  assert(scip != NULL);
3037  assert(cons != NULL);
3038 
3039  consdata = SCIPconsGetData(cons);
3040 
3041  have_change = FALSE;
3042  i = 0;
3043  while( i < consdata->nlinvars )
3044  {
3045  var = consdata->linvars[i];
3046 
3047  if( SCIPvarIsActive(var) && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
3048  {
3049  ++i;
3050  continue;
3051  }
3052 
3053  have_change = TRUE;
3054 
3055  coef = consdata->lincoefs[i];
3056  offset = 0.0;
3057 
3058  if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
3059  {
3060  offset = coef * (SCIPvarGetLbGlobal(var) + SCIPvarGetUbGlobal(var)) / 2.0;
3061  coef = 0.0;
3062  }
3063  else
3064  {
3065  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &coef, &offset) );
3066  }
3067 
3068  SCIPdebugMsg(scip, " linear term %g*<%s> is replaced by %g * <%s> + %g\n", consdata->lincoefs[i], SCIPvarGetName(consdata->linvars[i]),
3069  coef, SCIPvarGetName(var), offset);
3070 
3071  /* delete previous variable (this will move another variable to position i) */
3072  SCIP_CALL( delLinearCoefPos(scip, cons, i) );
3073 
3074  /* put constant part into bounds */
3075  if( offset != 0.0 )
3076  {
3077  if( !SCIPisInfinity(scip, -consdata->lhs) )
3078  {
3079  SCIP_CALL( chgSideQuadratic(scip, cons, eventhdlr, SCIP_SIDETYPE_LEFT, consdata->lhs - offset) );
3080  }
3081  if( !SCIPisInfinity(scip, consdata->rhs) )
3082  {
3083  SCIP_CALL( chgSideQuadratic(scip, cons, eventhdlr, SCIP_SIDETYPE_RIGHT, consdata->rhs - offset) );
3084  }
3085  }
3086 
3087  /* nothing left to do if variable had been fixed */
3088  if( coef == 0.0 )
3089  continue;
3090 
3091  /* if GetProbvar gave a linear variable, just add it
3092  * if it's a multilinear variable, add it's disaggregated variables */
3093  if( SCIPvarIsActive(var) )
3094  {
3095  SCIP_CALL( addLinearCoef(scip, cons, var, coef) );
3096  }
3097  else
3098  {
3099  int naggrs;
3100  SCIP_VAR** aggrvars;
3101  SCIP_Real* aggrscalars;
3102  SCIP_Real aggrconstant;
3103 
3104  assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_MULTAGGR);
3105 
3106  naggrs = SCIPvarGetMultaggrNVars(var);
3107  aggrvars = SCIPvarGetMultaggrVars(var);
3108  aggrscalars = SCIPvarGetMultaggrScalars(var);
3109  aggrconstant = SCIPvarGetMultaggrConstant(var);
3110 
3111  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars + naggrs) );
3112 
3113  for( j = 0; j < naggrs; ++j )
3114  {
3115  SCIP_CALL( addLinearCoef(scip, cons, aggrvars[j], coef * aggrscalars[j]) );
3116  }
3117 
3118  if( aggrconstant != 0.0 )
3119  {
3120  if( !SCIPisInfinity(scip, -consdata->lhs) )
3121  consdata->lhs -= coef * aggrconstant;
3122  if( !SCIPisInfinity(scip, consdata->rhs) )
3123  consdata->rhs -= coef * aggrconstant;
3124  }
3125  }
3126  }
3127 
3128  i = 0;
3129  while( i < consdata->nquadvars )
3130  {
3131  var = consdata->quadvarterms[i].var;
3132 
3133  if( SCIPvarIsActive(var) && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
3134  {
3135  ++i;
3136  continue;
3137  }
3138 
3139  have_change = TRUE;
3140 
3141  coef = 1.0;
3142  offset = 0.0;
3143 
3144  if( !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)) )
3145  {
3146  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &coef, &offset) );
3147  }
3148  else
3149  {
3150  coef = 0.0;
3151  offset = (SCIPvarGetLbGlobal(var) + SCIPvarGetUbGlobal(var)) / 2.0;
3152  }
3153 
3154  SCIPdebugMsg(scip, " quadratic variable <%s> with status %d is replaced by %g * <%s> + %g\n", SCIPvarGetName(consdata->quadvarterms[i].var),
3155  SCIPvarGetStatus(consdata->quadvarterms[i].var), coef, SCIPvarGetName(var), offset);
3156 
3157  /* handle fixed variable */
3158  if( coef == 0.0 )
3159  {
3160  /* if not fixed to 0.0, add to linear coefs of vars in bilinear terms, and deal with linear and square term as constant */
3161  if( offset != 0.0 )
3162  {
3163  for( j = 0; j < consdata->quadvarterms[i].nadjbilin; ++j )
3164  {
3165  bilinterm = &consdata->bilinterms[consdata->quadvarterms[i].adjbilin[j]];
3166 
3167  var2 = bilinterm->var1 == consdata->quadvarterms[i].var ? bilinterm->var2 : bilinterm->var1;
3168  assert(var2 != consdata->quadvarterms[i].var);
3169 
3170  var2pos = 0;
3171  while( consdata->quadvarterms[var2pos].var != var2 )
3172  {
3173  ++var2pos;
3174  assert(var2pos < consdata->nquadvars);
3175  }
3176  consdata->quadvarterms[var2pos].lincoef += bilinterm->coef * offset;
3177  }
3178 
3179  offset = consdata->quadvarterms[i].lincoef * offset + consdata->quadvarterms[i].sqrcoef * offset * offset;
3180  if( !SCIPisInfinity(scip, -consdata->lhs) )
3181  consdata->lhs -= offset;
3182  if( !SCIPisInfinity(scip, consdata->rhs) )
3183  consdata->rhs -= offset;
3184  }
3185 
3186  /* remove bilinear terms */
3187  SCIP_CALL( removeBilinearTermsPos(scip, cons, consdata->quadvarterms[i].nadjbilin, consdata->quadvarterms[i].adjbilin) );
3188 
3189  /* delete quad. var term i */
3190  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
3191 
3192  continue;
3193  }
3194 
3195  assert(var != NULL);
3196 
3197  /* if GetProbvar gave an active variable, replace the quad var term so that it uses the new variable */
3198  if( SCIPvarIsActive(var) )
3199  {
3200  /* replace x by coef*y+offset */
3201  SCIP_CALL( replaceQuadVarTermPos(scip, cons, i, var, coef, offset) );
3202 
3203  continue;
3204  }
3205  else
3206  {
3207  /* if GetProbVar gave a multi-aggregated variable, add new quad var terms and new bilinear terms
3208  * x is replaced by coef * (sum_i a_ix_i + b) + offset
3209  * lcoef * x + scoef * x^2 + bcoef * x * y ->
3210  * (b*coef + offset) * (lcoef + (b*coef + offset) * scoef)
3211  * + sum_i a_i*coef * (lcoef + 2 (b*coef + offset) * scoef) x_i
3212  * + sum_i (a_i*coef)^2 * scoef * x_i^2
3213  * + 2 sum_{i,j, i<j} (a_i a_j coef^2 scoef) x_i x_j
3214  * + bcoef * (b*coef + offset + coef * sum_i a_ix_i) y
3215  */
3216  int naggrs;
3217  SCIP_VAR** aggrvars; /* x_i */
3218  SCIP_Real* aggrscalars; /* a_i */
3219  SCIP_Real aggrconstant; /* b */
3220  int nquadtermsold;
3221 
3222  SCIP_Real lcoef;
3223  SCIP_Real scoef;
3224 
3225  assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_MULTAGGR);
3226 
3227  naggrs = SCIPvarGetMultaggrNVars(var);
3228  aggrvars = SCIPvarGetMultaggrVars(var);
3229  aggrscalars = SCIPvarGetMultaggrScalars(var);
3230  aggrconstant = SCIPvarGetMultaggrConstant(var);
3231 
3232  lcoef = consdata->quadvarterms[i].lincoef;
3233  scoef = consdata->quadvarterms[i].sqrcoef;
3234 
3235  nquadtermsold = consdata->nquadvars;
3236 
3237  SCIP_CALL( consdataEnsureQuadVarTermsSize(scip, consdata, consdata->nquadvars + naggrs) );
3238 
3239  /* take care of constant part */
3240  if( aggrconstant != 0.0 || offset != 0.0 )
3241  {
3242  SCIP_Real constant;
3243  constant = (aggrconstant * coef + offset) * (lcoef + (aggrconstant * coef + offset) * scoef);
3244  if( !SCIPisInfinity(scip, -consdata->lhs) )
3245  consdata->lhs -= constant;
3246  if( !SCIPisInfinity(scip, consdata->rhs) )
3247  consdata->rhs -= constant;
3248  }
3249 
3250  /* add x_i's with linear and square coefficients */
3251  for( j = 0; j < naggrs; ++j )
3252  {
3253  SCIP_CALL( addQuadVarTerm(scip, cons, aggrvars[j],
3254  coef * aggrscalars[j] * (lcoef + 2.0 * scoef * (coef * aggrconstant + offset)),
3255  coef * coef * aggrscalars[j] * aggrscalars[j] * scoef) );
3256  }
3257 
3258  /* ensure space for bilinear terms */
3259  SCIP_CALL( consdataEnsureBilinSize(scip, consdata, consdata->nquadvars + (scoef != 0.0 ? (naggrs * (naggrs-1))/2 : 0) + consdata->quadvarterms[j].nadjbilin * naggrs) );
3260 
3261  /* add x_j*x_k's */
3262  if( scoef != 0.0 )
3263  {
3264  for( j = 0; j < naggrs; ++j )
3265  for( k = 0; k < j; ++k )
3266  {
3267  assert(aggrvars[j] != aggrvars[k]);
3268  SCIP_CALL( addBilinearTerm(scip, cons, nquadtermsold + j, nquadtermsold + k,
3269  2.0 * aggrscalars[j] * aggrscalars[k] * coef * coef * scoef) );
3270  }
3271  }
3272 
3273  /* add x_i*y's */
3274  for( k = 0; k < consdata->quadvarterms[i].nadjbilin; ++k )
3275  {
3276  bilinterm = &consdata->bilinterms[consdata->quadvarterms[i].adjbilin[k]];
3277  bilincoef = bilinterm->coef; /* copy coef, as bilinterm pointer may become invalid by realloc in addBilinearTerm() below */
3278  var2 = (bilinterm->var1 == consdata->quadvarterms[i].var) ? bilinterm->var2 : bilinterm->var1;
3279  assert(var2 != consdata->quadvarterms[i].var);
3280 
3281  /* this is not efficient, but we cannot sort the quadratic terms here, since we currently iterate over them */
3282  var2pos = 0;
3283  while( consdata->quadvarterms[var2pos].var != var2 )
3284  {
3285  ++var2pos;
3286  assert(var2pos < consdata->nquadvars);
3287  }
3288 
3289  for( j = 0; j < naggrs; ++j )
3290  {
3291  if( aggrvars[j] == var2 )
3292  { /* x_i == y, so we have a square term here */
3293  consdata->quadvarterms[var2pos].sqrcoef += bilincoef * coef * aggrscalars[j];
3294  }
3295  else
3296  { /* x_i != y, so we need to add a bilinear term here */
3297  SCIP_CALL( addBilinearTerm(scip, cons, nquadtermsold + j, var2pos, bilincoef * coef * aggrscalars[j]) );
3298  }
3299  }
3300 
3301  consdata->quadvarterms[var2pos].lincoef += bilincoef * (aggrconstant * coef + offset);
3302  }
3303 
3304  /* remove bilinear terms */
3305  SCIP_CALL( removeBilinearTermsPos(scip, cons, consdata->quadvarterms[i].nadjbilin, consdata->quadvarterms[i].adjbilin) );
3306 
3307  /* delete quad. var term i */
3308  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
3309  }
3310  }
3311 
3312  consdata->isremovedfixings = TRUE;
3313 
3314  SCIPdebugMsg(scip, "removed fixations from <%s>\n -> ", SCIPconsGetName(cons));
3315  SCIPdebugPrintCons(scip, cons, NULL);
3316 
3317 #ifndef NDEBUG
3318  for( i = 0; i < consdata->nlinvars; ++i )
3319  assert(SCIPvarIsActive(consdata->linvars[i]));
3320 
3321  for( i = 0; i < consdata->nquadvars; ++i )
3322  assert(SCIPvarIsActive(consdata->quadvarterms[i].var));
3323 #endif
3324 
3325  if( !have_change )
3326  return SCIP_OKAY;
3327 
3328  /* some quadratic variable may have been replaced by an already existing linear variable
3329  * in this case, we want the linear variable to be removed, which happens in mergeAndCleanLinearVars
3330  */
3331  consdata->linvarsmerged = FALSE;
3332 
3333  SCIP_CALL( mergeAndCleanBilinearTerms(scip, cons) );
3334  SCIP_CALL( mergeAndCleanQuadVarTerms(scip, cons) );
3335  SCIP_CALL( mergeAndCleanLinearVars(scip, cons) );
3336 
3337 #ifndef NDEBUG
3338  for( i = 0; i < consdata->nbilinterms; ++i )
3339  {
3340  assert(consdata->bilinterms[i].var1 != consdata->bilinterms[i].var2);
3341  assert(consdata->bilinterms[i].coef != 0.0);
3342  assert(SCIPvarCompare(consdata->bilinterms[i].var1, consdata->bilinterms[i].var2) < 0);
3343  }
3344 #endif
3345 
3346  return SCIP_OKAY;
3347 }
3348 
3349 /** create a nonlinear row representation of the constraint and stores them in consdata */
3350 static
3351 SCIP_RETCODE createNlRow(
3352  SCIP* scip, /**< SCIP data structure */
3353  SCIP_CONS* cons /**< quadratic constraint */
3354  )
3355 {
3356  SCIP_CONSDATA* consdata;
3357  int nquadvars; /* number of variables in quadratic terms */
3358  SCIP_VAR** quadvars; /* variables in quadratic terms */
3359  int nquadelems; /* number of quadratic elements (square and bilinear terms) */
3360  SCIP_QUADELEM* quadelems; /* quadratic elements (square and bilinear terms) */
3361  int nquadlinterms; /* number of linear terms using variables that are in quadratic terms */
3362  SCIP_VAR** quadlinvars; /* variables of linear terms using variables that are in quadratic terms */
3363  SCIP_Real* quadlincoefs; /* coefficients of linear terms using variables that are in quadratic terms */
3364  int i;
3365  int idx1;
3366  int idx2;
3367  int lincnt;
3368  int elcnt;
3369  SCIP_VAR* lastvar;
3370  int lastvaridx;
3371  SCIP_EXPRCURV curvature;
3372 
3373  assert(scip != NULL);
3374  assert(cons != NULL);
3375 
3376  consdata = SCIPconsGetData(cons);
3377  assert(consdata != NULL);
3378 
3379  if( consdata->nlrow != NULL )
3380  {
3381  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
3382  }
3383 
3384  nquadvars = consdata->nquadvars;
3385  nquadelems = consdata->nbilinterms;
3386  nquadlinterms = 0;
3387  for( i = 0; i < nquadvars; ++i )
3388  {
3389  if( consdata->quadvarterms[i].sqrcoef != 0.0 )
3390  ++nquadelems;
3391  if( !SCIPisZero(scip, consdata->quadvarterms[i].lincoef) )
3392  ++nquadlinterms;
3393  }
3394 
3395  SCIP_CALL( SCIPallocBufferArray(scip, &quadvars, nquadvars) );
3396  SCIP_CALL( SCIPallocBufferArray(scip, &quadelems, nquadelems) );
3397  SCIP_CALL( SCIPallocBufferArray(scip, &quadlinvars, nquadlinterms) );
3398  SCIP_CALL( SCIPallocBufferArray(scip, &quadlincoefs, nquadlinterms) );
3399 
3400  lincnt = 0;
3401  elcnt = 0;
3402  for( i = 0; i < nquadvars; ++i )
3403  {
3404  quadvars[i] = consdata->quadvarterms[i].var;
3405 
3406  if( consdata->quadvarterms[i].sqrcoef != 0.0 )
3407  {
3408  assert(elcnt < nquadelems);
3409  quadelems[elcnt].idx1 = i;
3410  quadelems[elcnt].idx2 = i;
3411  quadelems[elcnt].coef = consdata->quadvarterms[i].sqrcoef;
3412  ++elcnt;
3413  }
3414 
3415  if( !SCIPisZero(scip, consdata->quadvarterms[i].lincoef) )
3416  {
3417  assert(lincnt < nquadlinterms);
3418  quadlinvars [lincnt] = consdata->quadvarterms[i].var;
3419  quadlincoefs[lincnt] = consdata->quadvarterms[i].lincoef;
3420  ++lincnt;
3421  }
3422  }
3423  assert(lincnt == nquadlinterms);
3424 
3425  /* bilinear terms are sorted first by first variable, then by second variable
3426  * thus, it makes sense to remember the index of the previous first variable for the case a series of bilinear terms with the same first var appears */
3427  lastvar = NULL;
3428  lastvaridx = -1;
3429  for( i = 0; i < consdata->nbilinterms; ++i )
3430  {
3431  if( lastvar == consdata->bilinterms[i].var1 )
3432  {
3433  assert(lastvaridx >= 0);
3434  assert(consdata->quadvarterms[lastvaridx].var == consdata->bilinterms[i].var1);
3435  }
3436  else
3437  {
3438  lastvar = consdata->bilinterms[i].var1;
3439  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, lastvar, &lastvaridx) );
3440  }
3441  idx1 = lastvaridx;
3442 
3443  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, consdata->bilinterms[i].var2, &idx2) );
3444 
3445  assert(elcnt < nquadelems);
3446  quadelems[elcnt].idx1 = MIN(idx1, idx2);
3447  quadelems[elcnt].idx2 = MAX(idx1, idx2);
3448  quadelems[elcnt].coef = consdata->bilinterms[i].coef;
3449  ++elcnt;
3450  }
3451  assert(elcnt == nquadelems);
3452 
3453  /* set curvature for the nonlinear row */
3454  if( consdata->isconcave && consdata->isconvex )
3455  {
3456  assert(consdata->nbilinterms == 0 && consdata->nquadvars == 0);
3457  curvature = SCIP_EXPRCURV_LINEAR;
3458  }
3459  else if( consdata->isconcave )
3460  curvature = SCIP_EXPRCURV_CONCAVE;
3461  else if( consdata->isconvex )
3462  curvature = SCIP_EXPRCURV_CONVEX;
3463  else
3464  curvature = SCIP_EXPRCURV_UNKNOWN;
3465 
3466  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
3467  consdata->nlinvars, consdata->linvars, consdata->lincoefs,
3468  nquadvars, quadvars, nquadelems, quadelems,
3469  NULL, consdata->lhs, consdata->rhs,
3470  curvature) );
3471 
3472  SCIP_CALL( SCIPaddLinearCoefsToNlRow(scip, consdata->nlrow, nquadlinterms, quadlinvars, quadlincoefs) );
3473 
3474  SCIPfreeBufferArray(scip, &quadlincoefs);
3475  SCIPfreeBufferArray(scip, &quadlinvars);
3476  SCIPfreeBufferArray(scip, &quadelems);
3477  SCIPfreeBufferArray(scip, &quadvars);
3478 
3479  return SCIP_OKAY;
3480 }
3481 
3482 /** solve constraint as presolving */
3483 static
3484 SCIP_RETCODE presolveSolve(
3485  SCIP* scip, /**< SCIP data structure */
3486  SCIP_CONS* cons, /**< constraint */
3487  SCIP_RESULT* result, /**< to store result of solve: cutoff, success, or do-not-find */
3488  SCIP_Bool* redundant, /**< to store whether constraint is redundant now (should be deleted) */
3489  int* naggrvars /**< counter on number of variable aggregations */
3490  )
3491 {
3492  SCIP_CONSDATA* consdata;
3493 
3494  assert(scip != NULL);
3495  assert(cons != NULL);
3496  assert(result != NULL);
3497  assert(redundant != NULL);
3498 
3499  *result = SCIP_DIDNOTFIND;
3500  *redundant = FALSE;
3501 
3502  consdata = SCIPconsGetData(cons);
3503  assert(consdata != NULL);
3504 
3505  /* if constraint is an equality with two variables, at least one of them binary,
3506  * and linear after fixing the binary, then we can aggregate the variables */
3507  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) && consdata->nlinvars == 0 && consdata->nquadvars == 2 &&
3508  ((SCIPvarIsBinary(consdata->quadvarterms[0].var) && consdata->quadvarterms[1].sqrcoef == 0.0) ||
3509  (SCIPvarIsBinary(consdata->quadvarterms[1].var) && consdata->quadvarterms[0].sqrcoef == 0.0)) )
3510  {
3511  SCIP_Bool infeasible;
3512  SCIP_Bool aggregated;
3513  SCIP_Real a;
3514  SCIP_Real b;
3515  SCIP_Real c;
3516  SCIP_VAR* x;
3517  SCIP_VAR* y;
3518  int binvaridx;
3519 
3520  /* constraint is a*(x+x^2) + b*y + c*x*y = rhs, with x binary variable
3521  * x = 0 -> b*y == rhs
3522  * x = 1 -> (b+c)*y == rhs - a
3523  *
3524  * if b != 0 and b+c != 0, then y = (rhs-a)/(b+c) * x + rhs/b * (1-x) = ((rhs-a)/(b+c) - rhs/b) * x + rhs/b
3525  */
3526 
3527  binvaridx = (SCIPvarIsBinary(consdata->quadvarterms[0].var) && consdata->quadvarterms[1].sqrcoef == 0.0) ? 0 : 1;
3528 
3529  x = consdata->quadvarterms[binvaridx].var;
3530  a = consdata->quadvarterms[binvaridx].sqrcoef + consdata->quadvarterms[binvaridx].lincoef;
3531 
3532  y = consdata->quadvarterms[1-binvaridx].var;
3533  b = consdata->quadvarterms[1-binvaridx].lincoef;
3534 
3535  assert(consdata->nbilinterms <= 1); /* should actually be 1, since constraint is otherwise linear */
3536  c = (consdata->nbilinterms == 1) ? consdata->bilinterms[0].coef : 0.0;
3537 
3538  if( !SCIPisZero(scip, b) && !SCIPisZero(scip, b+c) )
3539  {
3540  SCIPdebugMsg(scip, "<%s> = 0 -> %g*<%s> = %g and <%s> = 1 -> %g*<%s> = %g\n", SCIPvarGetName(x), b, SCIPvarGetName(y), consdata->rhs,
3541  SCIPvarGetName(x), b+c, SCIPvarGetName(y), consdata->rhs - a);
3542  SCIPdebugMsg(scip, "=> attempt aggregation <%s> = %g*<%s> + %g\n", SCIPvarGetName(y), (consdata->rhs-a)/(b+c) - consdata->rhs/b,
3543  SCIPvarGetName(x), consdata->rhs/b);
3544 
3545  SCIP_CALL( SCIPaggregateVars(scip, x, y, (consdata->rhs-a)/(b+c) - consdata->rhs/b, -1.0, -consdata->rhs/b, &infeasible, redundant, &aggregated) );
3546  if( infeasible )
3547  *result = SCIP_CUTOFF;
3548  else if( *redundant || aggregated )
3549  {
3550  /* aggregated (or were already aggregated), so constraint is now redundant */
3551  *result = SCIP_SUCCESS;
3552  *redundant = TRUE;
3553 
3554  if( aggregated )
3555  ++*naggrvars;
3556  }
3557  }
3558 
3559  /* @todo if b is 0 or b+c is 0, or lhs != rhs, then could replace by varbound constraint */
3560  }
3561 
3562  return SCIP_OKAY;
3563 }
3564 
3565 
3566 /** reformulates products of binary variables as AND constraint
3567  *
3568  * For a product x*y, with x and y binary variables, the product is replaced by a new auxiliary variable z and the constraint z = {x and y} is added.
3569  */
3570 static
3571 SCIP_RETCODE presolveTryAddAND(
3572  SCIP* scip, /**< SCIP data structure */
3573  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3574  SCIP_CONS* cons, /**< constraint */
3575  int* naddconss /**< buffer where to add the number of AND constraints added */
3576  )
3577 {
3578  SCIP_CONSHDLRDATA* conshdlrdata;
3579  SCIP_CONSDATA* consdata;
3580  char name[SCIP_MAXSTRLEN];
3581  SCIP_VAR* vars[2];
3582  SCIP_VAR* auxvar;
3583  SCIP_CONS* andcons;
3584  int i;
3585  int ntodelete;
3586  int* todelete;
3587 
3588  assert(scip != NULL);
3589  assert(conshdlr != NULL);
3590  assert(cons != NULL);
3591  assert(naddconss != NULL);
3592 
3593  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3594  assert(conshdlrdata != NULL);
3595 
3596  /* if no binary variables, then we will find nothing to reformulate here
3597  * (note that this does not count in integer variables with {0,1} bounds...)
3598  */
3599  if( SCIPgetNBinVars(scip) == 0 )
3600  return SCIP_OKAY;
3601 
3602  /* if user does not like AND very much, then return */
3603  if( conshdlrdata->empathy4and < 2 )
3604  return SCIP_OKAY;
3605 
3606  consdata = SCIPconsGetData(cons);
3607  assert(consdata != NULL);
3608 
3609  if( consdata->nbilinterms == 0 )
3610  return SCIP_OKAY;
3611 
3612  /* get array to store indices of bilinear terms that shall be deleted */
3613  SCIP_CALL( SCIPallocBufferArray(scip, &todelete, consdata->nbilinterms) );
3614  ntodelete = 0;
3615 
3616  for( i = 0; i < consdata->nbilinterms; ++i )
3617  {
3618  vars[0] = consdata->bilinterms[i].var1;
3619  if( !SCIPvarIsBinary(vars[0]) )
3620  continue;
3621 
3622  vars[1] = consdata->bilinterms[i].var2;
3623  if( !SCIPvarIsBinary(vars[1]) )
3624  continue;
3625 
3626  /* create auxiliary variable */
3627  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_%s", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]), SCIPconsGetName(cons));
3628  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, 0.0, 1.0, 0.0, SCIP_VARTYPE_BINARY,
3629  SCIPvarIsInitial(vars[0]) || SCIPvarIsInitial(vars[1]), SCIPvarIsRemovable(vars[0]) && SCIPvarIsRemovable(vars[1]), NULL, NULL, NULL, NULL, NULL) );
3630  SCIP_CALL( SCIPaddVar(scip, auxvar) );
3631 #ifdef WITH_DEBUG_SOLUTION
3632  if( SCIPdebugIsMainscip(scip) )
3633  {
3634  SCIP_Real var0val;
3635  SCIP_Real var1val;
3636  SCIP_CALL( SCIPdebugGetSolVal(scip, vars[0], &var0val) );
3637  SCIP_CALL( SCIPdebugGetSolVal(scip, vars[1], &var1val) );
3638  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, var0val * var1val) );
3639  }
3640 #endif
3641 
3642  /* create AND-constraint auxvar = x and y, need to be enforced as not redundant */
3643  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%sAND%s", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]));
3644  SCIP_CALL( SCIPcreateConsAnd(scip, &andcons, name, auxvar, 2, vars,
3645  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
3646  SCIPconsIsSeparated(cons), TRUE, TRUE,
3647  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
3648  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
3649  SCIP_CALL( SCIPaddCons(scip, andcons) );
3650  SCIPdebugMsg(scip, "added AND constraint: ");
3651  SCIPdebugPrintCons(scip, andcons, NULL);
3652  SCIP_CALL( SCIPreleaseCons(scip, &andcons) );
3653  ++*naddconss;
3654 
3655  /* add bilincoef * auxvar to linear terms */
3656  SCIP_CALL( addLinearCoef(scip, cons, auxvar, consdata->bilinterms[i].coef) );
3657  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
3658 
3659  /* remember that we have to delete this bilinear term */
3660  assert(ntodelete < consdata->nbilinterms);
3661  todelete[ntodelete++] = i;
3662  }
3663 
3664  /* remove bilinear terms that have been replaced */
3665  SCIP_CALL( removeBilinearTermsPos(scip, cons, ntodelete, todelete) );
3666  SCIPfreeBufferArray(scip, &todelete);
3667 
3668  return SCIP_OKAY;
3669 }
3670 
3671 /** gets bounds of variable y if x takes a certain value; checks whether x = xval has implications on y */
3672 static
3673 SCIP_RETCODE getImpliedBounds(
3674  SCIP* scip, /**< SCIP data structure */
3675  SCIP_VAR* x, /**< variable which implications to check */
3676  SCIP_Bool xval, /**< value of x to check for (TRUE for 1, FALSE for 0) */
3677  SCIP_VAR* y, /**< variable to check if bounds can be reduced */
3678  SCIP_INTERVAL* resultant /**< buffer to store bounds on y */
3679  )
3680 {
3681  SCIP_VAR** implvars;
3682  SCIP_BOUNDTYPE* impltypes;
3683  SCIP_Real* implbounds;
3684  int nimpls;
3685  int pos;
3686 
3687  assert(scip != NULL);
3688  assert(x != NULL);
3689  assert(y != NULL);
3690  assert(resultant != NULL);
3691 
3692  SCIPintervalSetBounds(resultant, MIN(SCIPvarGetLbGlobal(y), SCIPvarGetUbGlobal(y)), MAX(SCIPvarGetLbGlobal(y), SCIPvarGetUbGlobal(y))); /*lint !e666 */
3693 
3694  if( !SCIPvarIsBinary(x) || !SCIPvarIsActive(x) )
3695  return SCIP_OKAY;
3696 
3697  /* check in cliques for binary to binary implications */
3698  if( SCIPvarIsBinary(y) )
3699  {
3700  resultant->inf = MAX(resultant->inf, MIN(resultant->sup, 0.0));
3701  resultant->sup = MIN(resultant->sup, MAX(resultant->inf, 1.0));
3702 
3703  if( SCIPhaveVarsCommonClique(scip, x, xval, y, TRUE, FALSE) )
3704  {
3705  resultant->sup = MIN(resultant->sup, MAX(resultant->inf, 0.0));
3706  }
3707  else if( SCIPhaveVarsCommonClique(scip, x, xval, y, FALSE, FALSE) )
3708  {
3709  resultant->inf = MAX(resultant->inf, MIN(resultant->sup, 1.0));
3710  }
3711 
3712  return SCIP_OKAY;
3713  }
3714 
3715  /* analyze implications for x = xval */
3716  nimpls = SCIPvarGetNImpls(x, xval);
3717  if( nimpls == 0 )
3718  return SCIP_OKAY;
3719 
3720  implvars = SCIPvarGetImplVars (x, xval);
3721  impltypes = SCIPvarGetImplTypes (x, xval);
3722  implbounds = SCIPvarGetImplBounds(x, xval);
3723 
3724  assert(implvars != NULL);
3725  assert(impltypes != NULL);
3726  assert(implbounds != NULL);
3727 
3728  /* find implications */
3729  if( !SCIPsortedvecFindPtr((void**)implvars, SCIPvarComp, (void*)y, nimpls, &pos) )
3730  return SCIP_OKAY;
3731 
3732  /* if there are several implications on y, go to the first one */
3733  while( pos > 0 && implvars[pos-1] == y )
3734  --pos;
3735 
3736  /* update implied lower and upper bounds on y
3737  * but make sure that resultant will not be empty, due to tolerances
3738  */
3739  while( pos < nimpls && implvars[pos] == y )
3740  {
3741  if( impltypes[pos] == SCIP_BOUNDTYPE_LOWER )
3742  resultant->inf = MAX(resultant->inf, MIN(resultant->sup, implbounds[pos]));
3743  else
3744  resultant->sup = MIN(resultant->sup, MAX(resultant->inf, implbounds[pos]));
3745  ++pos;
3746  }
3747 
3748  assert(resultant->sup >= resultant->inf);
3749 
3750  return SCIP_OKAY;
3751 }
3752 
3753 /** Reformulates products of binary times bounded continuous variables as system of linear inequalities (plus auxiliary variable).
3754  *
3755  * For a product x*y, with y a binary variable and x a continuous variable with finite bounds,
3756  * an auxiliary variable z and the inequalities \f$ x^L y \leq z \leq x^U y \f$ and \f$ x - (1-y) x^U \leq z \leq x - (1-y) x^L \f$ are added.
3757  *
3758  * If x is a linear term consisting of more than one variable, it is split up in groups of linear terms of length at most maxnrvar.
3759  * For each product of linear term of length at most maxnrvar with y, an auxiliary z and linear inequalities are added.
3760  *
3761  * If y is a binary variable, the AND constraint \f$ z = x \wedge y \f$ may be added instead of linear constraints.
3762  */
3763 static
3764 SCIP_RETCODE presolveTryAddLinearReform(
3765  SCIP* scip, /**< SCIP data structure */
3766  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3767  SCIP_CONS* cons, /**< constraint */
3768  int* naddconss /**< buffer where to add the number of auxiliary constraints added */
3769  )
3770 { /*lint --e{666} */
3771  SCIP_CONSHDLRDATA* conshdlrdata;
3772  SCIP_CONSDATA* consdata;
3773  SCIP_VAR** xvars;
3774  SCIP_Real* xcoef;
3775  SCIP_INTERVAL xbndszero;
3776  SCIP_INTERVAL xbndsone;
3777  SCIP_INTERVAL act0;
3778  SCIP_INTERVAL act1;
3779  int nxvars;
3780  SCIP_VAR* y;
3781  SCIP_VAR* bvar;
3782  char name[SCIP_MAXSTRLEN];
3783  int nbilinterms;
3784  SCIP_VAR* auxvar;
3785  SCIP_CONS* auxcons;
3786  int i;
3787  int j;
3788  int k;
3789  int bilinidx;
3790  SCIP_Real bilincoef;
3791  SCIP_Real mincoef;
3792  SCIP_Real maxcoef;
3793  int* todelete;
3794  int ntodelete;
3795  int maxnrvar;
3796  SCIP_Bool integral;
3797  SCIP_Longint gcd;
3798  SCIP_Bool auxvarinitial;
3799  SCIP_Bool auxvarremovable;
3800 
3801  assert(scip != NULL);
3802  assert(conshdlr != NULL);
3803  assert(cons != NULL);
3804  assert(naddconss != NULL);
3805 
3806  /* if no binary variables, then we will find nothing to reformulate here
3807  * (note that this does not count in integer variables with {0,1} bounds...)
3808  */
3809  if( SCIPgetNBinVars(scip) == 0 )
3810  return SCIP_OKAY;
3811 
3812  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3813  assert(conshdlrdata != NULL);
3814 
3815  maxnrvar = conshdlrdata->replacebinaryprodlength;
3816  if( maxnrvar == 0 )
3817  return SCIP_OKAY;
3818 
3819  consdata = SCIPconsGetData(cons);
3820  assert(consdata != NULL);
3821 
3822  xvars = NULL;
3823  xcoef = NULL;
3824  todelete = NULL;
3825  gcd = 0;
3826 
3827  for( i = 0; i < consdata->nquadvars; ++i )
3828  {
3829  y = consdata->quadvarterms[i].var;
3830  if( !SCIPvarIsBinary(y) )
3831  continue;
3832 
3833  nbilinterms = consdata->quadvarterms[i].nadjbilin;
3834  if( nbilinterms == 0 )
3835  continue;
3836 
3837  SCIP_CALL( SCIPreallocBufferArray(scip, &xvars, MIN(maxnrvar, nbilinterms)+2) ); /* add 2 for later use when creating linear constraints */
3838  SCIP_CALL( SCIPreallocBufferArray(scip, &xcoef, MIN(maxnrvar, nbilinterms)+2) );
3839 
3840  /* alloc array to store indices of bilinear terms that shall be deleted */
3841  SCIP_CALL( SCIPreallocBufferArray(scip, &todelete, nbilinterms) );
3842  ntodelete = 0;
3843 
3844  auxvarinitial = SCIPvarIsInitial(y);
3845  auxvarremovable = SCIPvarIsRemovable(y);
3846 
3847  /* setup a list of bounded variables x_i with coefficients a_i that are multiplied with binary y: y*(sum_i a_i*x_i)
3848  * and compute range of sum_i a_i*x_i for the cases y = 0 and y = 1
3849  * we may need several rounds if maxnrvar < nbilinterms
3850  */
3851  j = 0;
3852  do
3853  {
3854  nxvars = 0;
3855  SCIPintervalSet(&xbndszero, 0.0);
3856  SCIPintervalSet(&xbndsone, 0.0);
3857 
3858  mincoef = SCIPinfinity(scip);
3859  maxcoef = 0.0;
3860  integral = TRUE;
3861 
3862  /* collect at most maxnrvar variables for x term */
3863  for( ; j < nbilinterms && nxvars < maxnrvar; ++j )
3864  {
3865  bilinidx = consdata->quadvarterms[i].adjbilin[j];
3866  assert(bilinidx >= 0);
3867  assert(bilinidx < consdata->nbilinterms);
3868 
3869  bvar = consdata->bilinterms[bilinidx].var1;
3870  if( bvar == y )
3871  bvar = consdata->bilinterms[bilinidx].var2;
3872  assert(bvar != y);
3873 
3874  /* skip products with unbounded variables */
3875  if( SCIPisInfinity(scip, -SCIPvarGetLbGlobal(bvar)) || SCIPisInfinity(scip, SCIPvarGetUbGlobal(bvar)) )
3876  {
3877  SCIPdebugMsg(scip, "skip reform of <%s><%s> due to unbounded second variable [%g,%g]\n",
3878  SCIPvarGetName(y), SCIPvarGetName(bvar), SCIPvarGetLbGlobal(bvar), SCIPvarGetUbGlobal(bvar));
3879  continue;
3880  }
3881 
3882  /* skip products with non-binary variables if binreformbinaryonly is set */
3883  if( conshdlrdata->binreformbinaryonly && !SCIPvarIsBinary(bvar) )
3884  {
3885  SCIPdebugMsg(scip, "skip reform of <%s><%s> because second variable is not binary\n",
3886  SCIPvarGetName(y), SCIPvarGetName(bvar));
3887  continue;
3888  }
3889 
3890  bilincoef = consdata->bilinterms[bilinidx].coef;
3891  assert(bilincoef != 0.0);
3892 
3893  /* get activity of bilincoef * x if y = 0 */
3894  SCIP_CALL( getImpliedBounds(scip, y, FALSE, bvar, &act0) );
3895  SCIPintervalMulScalar(SCIPinfinity(scip), &act0, act0, bilincoef);
3896 
3897  /* get activity of bilincoef * x if y = 1 */
3898  SCIP_CALL( getImpliedBounds(scip, y, TRUE, bvar, &act1) );
3899  SCIPintervalMulScalar(SCIPinfinity(scip), &act1, act1, bilincoef);
3900 
3901  /* skip products that give rise to very large coefficients (big big-M's) */
3902  if( SCIPfeastol(scip) * REALABS(act0.inf) >= conshdlrdata->binreformmaxcoef || SCIPfeastol(scip) * REALABS(act0.sup) >= conshdlrdata->binreformmaxcoef )
3903  {
3904  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity [%g,%g] for <%s> = 0.0\n",
3905  bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar), SCIPintervalGetInf(act0), SCIPintervalGetSup(act0), SCIPvarGetName(y));
3906  continue;
3907  }
3908  if( SCIPfeastol(scip) * REALABS(act1.inf) >= conshdlrdata->binreformmaxcoef || SCIPfeastol(scip) * REALABS(act1.sup) >= conshdlrdata->binreformmaxcoef )
3909  {
3910  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity [%g,%g] for <%s> = 1.0\n",
3911  bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar), SCIPintervalGetInf(act1), SCIPintervalGetSup(act1), SCIPvarGetName(y));
3912  continue;
3913  }
3914  if( !SCIPisZero(scip, MIN(REALABS(act0.inf), REALABS(act0.sup))) &&
3915  SCIPfeastol(scip) * MAX(REALABS(act0.inf), REALABS(act0.sup)) / MIN(REALABS(act0.inf), REALABS(act0.sup)) >= conshdlrdata->binreformmaxcoef )
3916  {
3917  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity ratio %g for <%s> = 0.0\n", bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar),
3918  MAX(REALABS(act0.inf), REALABS(act0.sup)) / MIN(REALABS(act0.inf), REALABS(act0.sup)), SCIPvarGetName(y));
3919  continue;
3920  }
3921  if( !SCIPisZero(scip, MIN(REALABS(act1.inf), REALABS(act1.sup))) &&
3922  SCIPfeastol(scip) * MAX(REALABS(act1.inf), REALABS(act1.sup)) / MIN(REALABS(act1.inf), REALABS(act1.sup)) >= conshdlrdata->binreformmaxcoef )
3923  {
3924  SCIPdebugMsg(scip, "skip reform of %g<%s><%s> due to huge activity ratio %g for <%s> = 0.0\n", bilincoef, SCIPvarGetName(y), SCIPvarGetName(bvar),
3925  MAX(REALABS(act1.inf), REALABS(act1.sup)) / MIN(REALABS(act1.inf), REALABS(act1.sup)), SCIPvarGetName(y));
3926  continue;
3927  }
3928 
3929  /* add bvar to x term */
3930  xvars[nxvars] = bvar;
3931  xcoef[nxvars] = bilincoef;
3932  ++nxvars;
3933 
3934  /* update bounds on x term */
3935  SCIPintervalAdd(SCIPinfinity(scip), &xbndszero, xbndszero, act0);
3936  SCIPintervalAdd(SCIPinfinity(scip), &xbndsone, xbndsone, act1);
3937 
3938  if( REALABS(bilincoef) < mincoef )
3939  mincoef = ABS(bilincoef);
3940  if( REALABS(bilincoef) > maxcoef )
3941  maxcoef = ABS(bilincoef);
3942 
3943  /* update whether all coefficients will be integral and if so, compute their gcd */
3944  integral &= (SCIPvarGetType(bvar) < SCIP_VARTYPE_CONTINUOUS) && SCIPisIntegral(scip, bilincoef); /*lint !e514 */
3945  if( integral )
3946  {
3947  if( nxvars == 1 )
3948  gcd = (SCIP_Longint)SCIPround(scip, REALABS(bilincoef));
3949  else
3950  gcd = SCIPcalcGreComDiv(gcd, (SCIP_Longint)SCIPround(scip, REALABS(bilincoef)));
3951  }
3952 
3953  /* if bvar is initial, then also the auxiliary variable should be initial
3954  * if bvar is not removable, then also the auxiliary variable should not be removable
3955  */
3956  auxvarinitial |= SCIPvarIsInitial(bvar);
3957  auxvarremovable &= SCIPvarIsRemovable(bvar);
3958 
3959  /* remember that we have to remove this bilinear term later */
3960  assert(ntodelete < nbilinterms);
3961  todelete[ntodelete++] = bilinidx;
3962  }
3963 
3964  if( nxvars == 0 ) /* all (remaining) x_j seem to be unbounded */
3965  break;
3966 
3967  assert(!SCIPisInfinity(scip, -SCIPintervalGetInf(xbndszero)));
3968  assert(!SCIPisInfinity(scip, SCIPintervalGetSup(xbndszero)));
3969  assert(!SCIPisInfinity(scip, -SCIPintervalGetInf(xbndsone)));
3970  assert(!SCIPisInfinity(scip, SCIPintervalGetSup(xbndsone)));
3971 
3972 #ifdef SCIP_DEBUG
3973  if( SCIPintervalGetInf(xbndszero) != SCIPintervalGetInf(xbndsone) || /*lint !e777*/
3974  +SCIPintervalGetSup(xbndszero) != SCIPintervalGetSup(xbndsone) ) /*lint !e777*/
3975  {
3976  SCIPdebugMsg(scip, "got different bounds for y = 0: [%g, %g] and y = 1: [%g, %g]\n", xbndszero.inf, xbndszero.sup, xbndsone.inf, xbndsone.sup);
3977  }
3978 #endif
3979 
3980  if( nxvars == 1 && conshdlrdata->empathy4and >= 1 && SCIPvarIsBinary(xvars[0]) )
3981  {
3982  /* product of two binary variables, replace by auxvar and AND constraint */
3983  /* add auxiliary variable z */
3984  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
3985  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, 0.0, 1.0, 0.0, SCIP_VARTYPE_IMPLINT,
3986  auxvarinitial, auxvarremovable, NULL, NULL, NULL, NULL, NULL) );
3987  SCIP_CALL( SCIPaddVar(scip, auxvar) );
3988 
3989 #ifdef WITH_DEBUG_SOLUTION
3990  if( SCIPdebugIsMainscip(scip) )
3991  {
3992  SCIP_Real var0val;
3993  SCIP_Real var1val;
3994  SCIP_CALL( SCIPdebugGetSolVal(scip, xvars[0], &var0val) );
3995  SCIP_CALL( SCIPdebugGetSolVal(scip, y, &var1val) );
3996  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, var0val * var1val) );
3997  }
3998 #endif
3999 
4000  /* add constraint z = x and y; need to be enforced, as it is not redundant w.r.t. existing constraints */
4001  xvars[1] = y;
4002  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%sAND%s_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4003  SCIP_CALL( SCIPcreateConsAnd(scip, &auxcons, name, auxvar, 2, xvars,
4004  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
4005  SCIPconsIsSeparated(cons), TRUE, TRUE,
4006  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4007  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4008  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4009  SCIPdebugMsg(scip, "added AND constraint: ");
4010  SCIPdebugPrintCons(scip, auxcons, NULL);
4011  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4012  ++*naddconss;
4013 
4014  /* add linear term coef*auxvar */
4015  SCIP_CALL( addLinearCoef(scip, cons, auxvar, xcoef[0]) );
4016 
4017  /* forget about auxvar */
4018  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4019  }
4020  else
4021  {
4022  /* product of binary variable with more than one binary or with continuous variables or with binary and user
4023  * did not like AND -> replace by auxvar and linear constraints */
4024  SCIP_Real scale;
4025 
4026  /* scale auxiliary constraint by some nice value,
4027  * if all coefficients are integral, take a value that preserves integrality (-> gcd), so we can make the auxiliary variable impl. integer
4028  */
4029  if( integral )
4030  {
4031  scale = (SCIP_Real)gcd;
4032  assert(scale >= 1.0);
4033  }
4034  else if( nxvars == 1 )
4035  {
4036  /* scaling by the only coefficient gives auxiliary variable = x * y, which thus will be implicit integral provided y is not continuous */
4037  assert(mincoef == maxcoef); /*lint !e777 */
4038  scale = mincoef;
4039  integral = SCIPvarGetType(xvars[0]) < SCIP_VARTYPE_CONTINUOUS;
4040  }
4041  else
4042  {
4043  scale = 1.0;
4044  if( maxcoef < 0.5 )
4045  scale = maxcoef;
4046  if( mincoef > 2.0 )
4047  scale = mincoef;
4048  if( scale != 1.0 )
4049  scale = SCIPselectSimpleValue(scale / 2.0, 1.5 * scale, MAXDNOM);
4050  }
4051  assert(scale > 0.0);
4052  assert(!SCIPisInfinity(scip, scale));
4053 
4054  /* if x-term is always negative for y = 1, negate scale so we get a positive auxiliary variable; maybe this is better sometimes? */
4055  if( !SCIPisPositive(scip, SCIPintervalGetSup(xbndsone)) )
4056  scale = -scale;
4057 
4058  SCIPdebugMsg(scip, "binary reformulation using scale %g, nxvars = %d, integral = %u\n", scale, nxvars, integral);
4059  if( scale != 1.0 )
4060  {
4061  SCIPintervalDivScalar(SCIPinfinity(scip), &xbndszero, xbndszero, scale);
4062  SCIPintervalDivScalar(SCIPinfinity(scip), &xbndsone, xbndsone, scale);
4063  for( k = 0; k < nxvars; ++k )
4064  xcoef[k] /= scale;
4065  }
4066 
4067  /* add auxiliary variable z */
4068  if( nxvars == 1 )
4069  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4070  else
4071  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "prod%s_%s_more_%s", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4072  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, MIN(0., SCIPintervalGetInf(xbndsone)), MAX(0., SCIPintervalGetSup(xbndsone)),
4073  0.0, integral ? SCIP_VARTYPE_IMPLINT : SCIP_VARTYPE_CONTINUOUS,
4074  auxvarinitial, auxvarremovable, NULL, NULL, NULL, NULL, NULL) );
4075  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4076 
4077  /* compute value of auxvar in debug solution */
4078 #ifdef WITH_DEBUG_SOLUTION
4079  if( SCIPdebugIsMainscip(scip) )
4080  {
4081  SCIP_Real debugval;
4082  SCIP_Real varval;
4083 
4084  SCIP_CALL( SCIPdebugGetSolVal(scip, y, &varval) );
4085  if( SCIPisZero(scip, varval) )
4086  {
4087  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, 0.0) );
4088  }
4089  else
4090  {
4091  assert(SCIPisEQ(scip, varval, 1.0));
4092 
4093  debugval = 0.0;
4094  for( k = 0; k < nxvars; ++k )
4095  {
4096  SCIP_CALL( SCIPdebugGetSolVal(scip, xvars[k], &varval) );
4097  debugval += xcoef[k] * varval;
4098  }
4099  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
4100  }
4101  }
4102 #endif
4103 
4104  /* add auxiliary constraints
4105  * it seems to be advantageous to make the varbound constraints initial and the linear constraints not initial
4106  * maybe because it is more likely that a binary variable takes value 0 instead of 1, and thus the varbound constraints
4107  * are more often active, compared to the linear constraints added below
4108  * also, the varbound constraints are more sparse than the linear cons
4109  */
4110  if( SCIPisNegative(scip, SCIPintervalGetInf(xbndsone)) )
4111  {
4112  /* add 0 <= z - xbndsone.inf * y constraint (as varbound constraint), need to be enforced as not redundant */
4113  if( nxvars == 1 )
4114  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_1", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4115  else
4116  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_1", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4117  SCIP_CALL( SCIPcreateConsVarbound(scip, &auxcons, name, auxvar, y, -SCIPintervalGetInf(xbndsone), 0.0, SCIPinfinity(scip),
4118  SCIPconsIsInitial(cons) /*&& conshdlrdata->binreforminitial*/,
4119  SCIPconsIsSeparated(cons), TRUE, TRUE,
4120  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4121  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4122  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4123  SCIPdebugMsg(scip, "added varbound constraint: ");
4124  SCIPdebugPrintCons(scip, auxcons, NULL);
4125  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4126  ++*naddconss;
4127  }
4128  if( SCIPisPositive(scip, SCIPintervalGetSup(xbndsone)) )
4129  {
4130  /* add z - xbndsone.sup * y <= 0 constraint (as varbound constraint), need to be enforced as not redundant */
4131  if( nxvars == 1 )
4132  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_2", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4133  else
4134  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_2", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4135  SCIP_CALL( SCIPcreateConsVarbound(scip, &auxcons, name, auxvar, y, -SCIPintervalGetSup(xbndsone), -SCIPinfinity(scip), 0.0,
4136  SCIPconsIsInitial(cons) /*&& conshdlrdata->binreforminitial*/,
4137  SCIPconsIsSeparated(cons), TRUE, TRUE,
4138  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4139  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4140  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4141  SCIPdebugMsg(scip, "added varbound constraint: ");
4142  SCIPdebugPrintCons(scip, auxcons, NULL);
4143  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4144  ++*naddconss;
4145  }
4146 
4147  /* add xbndszero.inf <= sum_i a_i*x_i + xbndszero.inf * y - z constraint, need to be enforced as not redundant */
4148  xvars[nxvars] = y;
4149  xvars[nxvars+1] = auxvar;
4150  xcoef[nxvars] = SCIPintervalGetInf(xbndszero);
4151  xcoef[nxvars+1] = -1;
4152 
4153  if( nxvars == 1 )
4154  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_3", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4155  else
4156  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_3", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4157  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, name, nxvars+2, xvars, xcoef, SCIPintervalGetInf(xbndszero), SCIPinfinity(scip),
4158  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
4159  SCIPconsIsSeparated(cons), TRUE, TRUE,
4160  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4161  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4162  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4163  SCIPdebugMsg(scip, "added linear constraint: ");
4164  SCIPdebugPrintCons(scip, auxcons, NULL);
4165  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4166  ++*naddconss;
4167 
4168  /* add sum_i a_i*x_i + xbndszero.sup * y - z <= xbndszero.sup constraint, need to be enforced as not redundant */
4169  xcoef[nxvars] = SCIPintervalGetSup(xbndszero);
4170 
4171  if( nxvars == 1 )
4172  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s_%s_4", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4173  else
4174  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "linreform%s*%s*more_%s_4", SCIPvarGetName(y), SCIPvarGetName(xvars[0]), SCIPconsGetName(cons));
4175  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, name, nxvars+2, xvars, xcoef, -SCIPinfinity(scip), SCIPintervalGetSup(xbndszero),
4176  SCIPconsIsInitial(cons) && conshdlrdata->binreforminitial,
4177  SCIPconsIsSeparated(cons), TRUE, TRUE,
4178  SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4179  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons), SCIPconsIsStickingAtNode(cons)) );
4180  SCIP_CALL( SCIPaddCons(scip, auxcons) );
4181  SCIPdebugMsg(scip, "added linear constraint: ");
4182  SCIPdebugPrintCons(scip, auxcons, NULL);
4183  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
4184  ++*naddconss;
4185 
4186  /* add linear term scale*auxvar to this constraint */
4187  SCIP_CALL( addLinearCoef(scip, cons, auxvar, scale) );
4188 
4189  /* forget about auxvar */
4190  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4191  }
4192  }
4193  while( j < nbilinterms );
4194 
4195  /* remove bilinear terms that have been replaced */
4196  SCIP_CALL( removeBilinearTermsPos(scip, cons, ntodelete, todelete) );
4197  }
4198  SCIPdebugMsg(scip, "resulting quadratic constraint: ");
4199  SCIPdebugPrintCons(scip, cons, NULL);
4200 
4201  SCIPfreeBufferArrayNull(scip, &todelete);
4202  SCIPfreeBufferArrayNull(scip, &xcoef);
4203  SCIPfreeBufferArrayNull(scip, &xvars);
4204 
4205  return SCIP_OKAY;
4206 }
4207 
4208 /** tries to automatically convert a quadratic constraint (or a part of it) into a more specific and more specialized constraint */
4209 static
4210 SCIP_RETCODE presolveUpgrade(
4211  SCIP* scip, /**< SCIP data structure */
4212  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
4213  SCIP_CONS* cons, /**< source constraint to try to convert */
4214  SCIP_Bool* upgraded, /**< buffer to store whether constraint was upgraded */
4215  int* nupgdconss, /**< buffer to increase if constraint was upgraded */
4216  int* naddconss, /**< buffer to increase with number of additional constraints created during upgrade */
4217  SCIP_PRESOLTIMING presoltiming /**< current presolving timing */
4218  )
4219 {
4220  SCIP_CONSHDLRDATA* conshdlrdata;
4221  SCIP_CONSDATA* consdata;
4222  SCIP_VAR* var;
4223  SCIP_Real lincoef;
4224  SCIP_Real quadcoef;
4225  SCIP_Real lb;
4226  SCIP_Real ub;
4227  int nbinlin;
4228  int nbinquad;
4229  int nintlin;
4230  int nintquad;
4231  int nimpllin;
4232  int nimplquad;
4233  int ncontlin;
4234  int ncontquad;
4235  SCIP_Bool integral;
4236  int i;
4237  int j;
4238  SCIP_CONS** upgdconss;
4239  int upgdconsssize;
4240  int nupgdconss_;
4241 
4242  assert(scip != NULL);
4243  assert(conshdlr != NULL);
4244  assert(cons != NULL);
4245  assert(!SCIPconsIsModifiable(cons));
4246  assert(upgraded != NULL);
4247  assert(nupgdconss != NULL);
4248  assert(naddconss != NULL);
4249 
4250  *upgraded = FALSE;
4251 
4252  nupgdconss_ = 0;
4253 
4254  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4255  assert(conshdlrdata != NULL);
4256 
4257  /* if there are no upgrade methods, we can also stop */
4258  if( conshdlrdata->nquadconsupgrades == 0 )
4259  return SCIP_OKAY;
4260 
4261  upgdconsssize = 2;
4262  SCIP_CALL( SCIPallocBufferArray(scip, &upgdconss, upgdconsssize) );
4263 
4264  consdata = SCIPconsGetData(cons);
4265  assert(consdata != NULL);
4266 
4267  /* calculate some statistics on quadratic constraint */
4268  nbinlin = 0;
4269  nbinquad = 0;
4270  nintlin = 0;
4271  nintquad = 0;
4272  nimpllin = 0;
4273  nimplquad = 0;
4274  ncontlin = 0;
4275  ncontquad = 0;
4276  integral = TRUE;
4277  for( i = 0; i < consdata->nlinvars; ++i )
4278  {
4279  var = consdata->linvars[i];
4280  lincoef = consdata->lincoefs[i];
4281  lb = SCIPvarGetLbLocal(var);
4282  ub = SCIPvarGetUbLocal(var);
4283  assert(!SCIPisZero(scip, lincoef));
4284 
4285  switch( SCIPvarGetType(var) )
4286  {
4287  case SCIP_VARTYPE_BINARY:
4288  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4289  integral = integral && SCIPisIntegral(scip, lincoef);
4290  nbinlin++;
4291  break;
4292  case SCIP_VARTYPE_INTEGER:
4293  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4294  integral = integral && SCIPisIntegral(scip, lincoef);
4295  nintlin++;
4296  break;
4297  case SCIP_VARTYPE_IMPLINT:
4298  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4299  integral = integral && SCIPisIntegral(scip, lincoef);
4300  nimpllin++;
4301  break;
4302  case SCIP_VARTYPE_CONTINUOUS:
4303  integral = integral && SCIPisRelEQ(scip, lb, ub) && SCIPisIntegral(scip, lincoef * lb);
4304  ncontlin++;
4305  break;
4306  default:
4307  SCIPerrorMessage("unknown variable type\n");
4308  return SCIP_INVALIDDATA;
4309  }
4310  }
4311 
4312  for( i = 0; i < consdata->nquadvars; ++i )
4313  {
4314  var = consdata->quadvarterms[i].var;
4315  lincoef = consdata->quadvarterms[i].lincoef;
4316  quadcoef = consdata->quadvarterms[i].sqrcoef;
4317  lb = SCIPvarGetLbLocal(var);
4318  ub = SCIPvarGetUbLocal(var);
4319 
4320  switch( SCIPvarGetType(var) )
4321  {
4322  case SCIP_VARTYPE_BINARY:
4323  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4324  integral = integral && SCIPisIntegral(scip, lincoef) && SCIPisIntegral(scip, quadcoef);
4325  nbinquad++;
4326  break;
4327  case SCIP_VARTYPE_INTEGER:
4328  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4329  integral = integral && SCIPisIntegral(scip, lincoef) && SCIPisIntegral(scip, quadcoef);
4330  nintquad++;
4331  break;
4332  case SCIP_VARTYPE_IMPLINT:
4333  if( !SCIPisZero(scip, lb) || !SCIPisZero(scip, ub) )
4334  integral = integral && SCIPisIntegral(scip, lincoef) && SCIPisIntegral(scip, quadcoef);
4335  nimplquad++;
4336  break;
4337  case SCIP_VARTYPE_CONTINUOUS:
4338  integral = integral && SCIPisRelEQ(scip, lb, ub) && SCIPisIntegral(scip, lincoef * lb + quadcoef * lb * lb);
4339  ncontquad++;
4340  break;
4341  default:
4342  SCIPerrorMessage("unknown variable type\n");
4343  return SCIP_INVALIDDATA;
4344  }
4345  }
4346 
4347  if( integral )
4348  {
4349  for( i = 0; i < consdata->nbilinterms && integral; ++i )
4350  {
4351  if( SCIPvarGetType(consdata->bilinterms[i].var1) < SCIP_VARTYPE_CONTINUOUS && SCIPvarGetType(consdata->bilinterms[i].var2) < SCIP_VARTYPE_CONTINUOUS )
4352  integral = integral && SCIPisIntegral(scip, consdata->bilinterms[i].coef);
4353  else
4354  integral = FALSE;
4355  }
4356  }
4357 
4358  /* call the upgrading methods */
4359 
4360  SCIPdebugMsg(scip, "upgrading quadratic constraint <%s> (%d upgrade methods):\n",
4361  SCIPconsGetName(cons), conshdlrdata->nquadconsupgrades);
4362  SCIPdebugMsg(scip, " binlin=%d binquad=%d intlin=%d intquad=%d impllin=%d implquad=%d contlin=%d contquad=%d integral=%u\n",
4363  nbinlin, nbinquad, nintlin, nintquad, nimpllin, nimplquad, ncontlin, ncontquad, integral);
4364  SCIPdebugPrintCons(scip, cons, NULL);
4365 
4366  /* try all upgrading methods in priority order in case the upgrading step is enable */
4367  for( i = 0; i < conshdlrdata->nquadconsupgrades; ++i )
4368  {
4369  if( !conshdlrdata->quadconsupgrades[i]->active )
4370  continue;
4371 
4372  SCIP_CALL( conshdlrdata->quadconsupgrades[i]->quadconsupgd(scip, cons,
4373  nbinlin, nbinquad, nintlin, nintquad, nimpllin, nimplquad, ncontlin, ncontquad, integral,
4374  &nupgdconss_, upgdconss, upgdconsssize, presoltiming) );
4375 
4376  while( nupgdconss_ < 0 )
4377  {
4378  /* upgrade function requires more memory: resize upgdconss and call again */
4379  assert(-nupgdconss_ > upgdconsssize);
4380  upgdconsssize = -nupgdconss_;
4381  SCIP_CALL( SCIPreallocBufferArray(scip, &upgdconss, -nupgdconss_) );
4382 
4383  SCIP_CALL( conshdlrdata->quadconsupgrades[i]->quadconsupgd(scip, cons,
4384  nbinlin, nbinquad, nintlin, nintquad, nimpllin, nimplquad, ncontlin, ncontquad, integral,
4385  &nupgdconss_, upgdconss, upgdconsssize, presoltiming) );
4386 
4387  assert(nupgdconss_ != 0);
4388  }
4389 
4390  if( nupgdconss_ > 0 )
4391  {
4392  /* got upgrade */
4393  SCIPdebugPrintCons(scip, cons, NULL);
4394  SCIPdebugMsg(scip, " -> upgraded to %d constraints:\n", nupgdconss_);
4395 
4396  /* add the upgraded constraints to the problem and forget them */
4397  for( j = 0; j < nupgdconss_; ++j )
4398  {
4399  SCIPdebugMsgPrint(scip, "\t");
4400  SCIPdebugPrintCons(scip, upgdconss[j], NULL);
4401 
4402  SCIP_CALL( SCIPaddCons(scip, upgdconss[j]) ); /*lint !e613*/
4403  SCIP_CALL( SCIPreleaseCons(scip, &upgdconss[j]) ); /*lint !e613*/
4404  }
4405 
4406  /* count the first upgrade constraint as constraint upgrade and the remaining ones as added constraints */
4407  *nupgdconss += 1;
4408  *naddconss += nupgdconss_ - 1;
4409  *upgraded = TRUE;
4410 
4411  /* delete upgraded constraint */
4412  SCIPdebugMsg(scip, "delete constraint <%s> after upgrade\n", SCIPconsGetName(cons));
4413  SCIP_CALL( SCIPdelCons(scip, cons) );
4414 
4415  break;
4416  }
4417  }
4418 
4419  SCIPfreeBufferArray(scip, &upgdconss);
4420 
4421  return SCIP_OKAY;
4422 }
4423 
4424 /** helper function for presolveDisaggregate */
4425 static
4426 SCIP_RETCODE presolveDisaggregateMarkComponent(
4427  SCIP* scip, /**< SCIP data structure */
4428  SCIP_CONSDATA* consdata, /**< constraint data */
4429  int quadvaridx, /**< index of quadratic variable to mark */
4430  SCIP_HASHMAP* var2component, /**< variables to components mapping */
4431  int componentnr, /**< the component number to mark to */
4432  int* componentsize /**< buffer to store size of component (incremented by 1) */
4433  )
4434 {
4435  SCIP_QUADVARTERM* quadvarterm;
4436  SCIP_VAR* othervar;
4437  int othervaridx;
4438  int i;
4439 
4440  assert(consdata != NULL);
4441  assert(quadvaridx >= 0);
4442  assert(quadvaridx < consdata->nquadvars);
4443  assert(var2component != NULL);
4444  assert(componentnr >= 0);
4445 
4446  quadvarterm = &consdata->quadvarterms[quadvaridx];
4447 
4448  if( SCIPhashmapExists(var2component, quadvarterm->var) )
4449  {
4450  /* if we saw the variable before, then it should have the same component number */
4451  assert(SCIPhashmapGetImageInt(var2component, quadvarterm->var) == componentnr);
4452  return SCIP_OKAY;
4453  }
4454 
4455  /* assign component number to variable */
4456  SCIP_CALL( SCIPhashmapInsertInt(var2component, quadvarterm->var, componentnr) );
4457  ++*componentsize;
4458 
4459  /* assign same component number to all variables this variable is multiplied with */
4460  for( i = 0; i < quadvarterm->nadjbilin; ++i )
4461  {
4462  othervar = consdata->bilinterms[quadvarterm->adjbilin[i]].var1 == quadvarterm->var ?
4463  consdata->bilinterms[quadvarterm->adjbilin[i]].var2 : consdata->bilinterms[quadvarterm->adjbilin[i]].var1;
4464  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, othervar, &othervaridx) );
4465  assert(othervaridx >= 0);
4466  SCIP_CALL( presolveDisaggregateMarkComponent(scip, consdata, othervaridx, var2component, componentnr, componentsize) );
4467  }
4468 
4469  return SCIP_OKAY;
4470 }
4471 
4472 /** merges components in variables connectivity graph */
4473 static
4474 SCIP_RETCODE presolveDisaggregateMergeComponents(
4475  SCIP* scip, /**< SCIP data structure */
4476  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
4477  SCIP_HASHMAP* var2component, /**< variables to component mapping */
4478  int nvars, /**< number of variables */
4479  int* ncomponents, /**< number of components */
4480  int* componentssize /**< size of components */
4481  )
4482 {
4483  SCIP_CONSHDLRDATA* conshdlrdata;
4484  SCIP_HASHMAPENTRY* entry;
4485  int maxncomponents;
4486  int* oldcompidx;
4487  int* newcompidx;
4488  int i;
4489  int oldcomponent;
4490  int newcomponent;
4491 
4492  assert(scip != NULL);
4493  assert(conshdlr != NULL);
4494  assert(var2component != NULL);
4495  assert(ncomponents != NULL);
4496  assert(componentssize != NULL);
4497 
4498  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4499  assert(conshdlrdata != NULL);
4500 
4501  maxncomponents = conshdlrdata->maxdisaggrsize;
4502  assert(maxncomponents > 0);
4503 
4504  /* if already not too many components, then nothing to do */
4505  if( *ncomponents <= maxncomponents )
4506  return SCIP_OKAY;
4507 
4508  /*
4509  printf("component sizes before:");
4510  for( i = 0; i < *ncomponents; ++i )
4511  printf(" %d", componentssize[i]);
4512  printf("\n");
4513  */
4514 
4515  SCIP_CALL( SCIPallocBufferArray(scip, &oldcompidx, *ncomponents) );
4516  SCIP_CALL( SCIPallocBufferArray(scip, &newcompidx, *ncomponents) );
4517 
4518  for( i = 0; i < *ncomponents; ++i )
4519  oldcompidx[i] = i;
4520 
4521  switch( conshdlrdata->disaggrmergemethod )
4522  {
4523  case 's' :
4524  /* sort components by size, increasing order */
4525  SCIPsortIntInt(componentssize, oldcompidx, *ncomponents);
4526  break;
4527  case 'b' :
4528  case 'm' :
4529  /* sort components by size, decreasing order */
4530  SCIPsortDownIntInt(componentssize, oldcompidx, *ncomponents);
4531  break;
4532  default :
4533  SCIPerrorMessage("invalid value for constraints/quadratic/disaggrmergemethod parameter");
4534  return SCIP_PARAMETERWRONGVAL;
4535  }
4536 
4537  SCIPdebugMsg(scip, "%-30s: % 4d components of size % 4d to % 4d, median: % 4d\n", SCIPgetProbName(scip), *ncomponents, componentssize[0], componentssize[*ncomponents-1], componentssize[*ncomponents/2]);
4538 
4539  if( conshdlrdata->disaggrmergemethod == 'm' )
4540  {
4541  SCIP_Real targetsize;
4542  int count = 0;
4543 
4544  /* a minimal component size we should reach to have all components roughly the same size */
4545  targetsize = nvars / maxncomponents; /*lint !e653*/
4546  for( i = 0; i < *ncomponents; ++i )
4547  {
4548  newcompidx[oldcompidx[i]] = i;
4549  count += componentssize[i];
4550 
4551  /* fill with small components until we reach targetsize
4552  * Since targetsize might be fractional, we also add another component if
4553  * the number of variables remaining (=nvars-count) is larger than
4554  * what we expect to put into the remaining components (=targetsize * (maxncomponents - i-1)).
4555  * Thus, from time to time, a component is made larger than the targetsize to avoid
4556  * having to add much into the last component.
4557  */
4558  while( i < *ncomponents-1 && (componentssize[i] + componentssize[*ncomponents-1] <= targetsize ||
4559  nvars - count > targetsize * (maxncomponents - i)) )
4560  {
4561  /* map last (=smallest) component to component i */
4562  newcompidx[oldcompidx[*ncomponents-1]] = i;
4563 
4564  /* increase size of component i accordingly */
4565  componentssize[i] += componentssize[*ncomponents-1];
4566  count += componentssize[*ncomponents-1];
4567 
4568  /* forget about last component */
4569  --*ncomponents;
4570  }
4571  }
4572  assert(count == nvars);
4573  }
4574  else
4575  {
4576  /* get inverse permutation */
4577  for( i = 0; i < *ncomponents; ++i )
4578  newcompidx[oldcompidx[i]] = i;
4579  }
4580 
4581  /* assign new component numbers to variables, cutting off at maxncomponents */
4582  for( i = 0; i < SCIPhashmapGetNEntries(var2component); ++i )
4583  {
4584  entry = SCIPhashmapGetEntry(var2component, i);
4585  if( entry == NULL )
4586  continue;
4587 
4588  oldcomponent = (int)(size_t)SCIPhashmapEntryGetImage(entry);
4589 
4590  newcomponent = newcompidx[oldcomponent];
4591  if( newcomponent >= maxncomponents )
4592  {
4593  newcomponent = maxncomponents-1;
4594  ++componentssize[maxncomponents-1];
4595  }
4596 
4597  SCIPhashmapEntrySetImage(entry, (void*)(size_t)newcomponent); /*lint !e571*/
4598  }
4599  if( *ncomponents > maxncomponents )
4600  *ncomponents = maxncomponents;
4601 
4602  /*
4603  printf("component sizes after :");
4604  for( i = 0; i < *ncomponents; ++i )
4605  printf(" %d", componentssize[i]);
4606  printf("\n");
4607  */
4608 
4609  SCIPfreeBufferArray(scip, &newcompidx);
4610  SCIPfreeBufferArray(scip, &oldcompidx);
4611 
4612  return SCIP_OKAY;
4613 }
4614 
4615 /** compute the next highest power of 2 for a 32-bit argument
4616  *
4617  * Source: https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
4618  *
4619  * @note Returns 0 for v=0.
4620  */
4621 static
4622 unsigned int nextPowerOf2(
4623  unsigned int v /**< input */
4624  )
4625 {
4626  v--;
4627  v |= v >> 1;
4628  v |= v >> 2;
4629  v |= v >> 4;
4630  v |= v >> 8;
4631  v |= v >> 16;
4632  v++;
4633 
4634  return v;
4635 }
4636 
4637 
4638 /** for quadratic constraints that consists of a sum of quadratic terms, disaggregates the sum into a set of constraints by introducing auxiliary variables
4639  *
4640  * Assume the quadratic constraint can be written in the form
4641  * lhs <= b'x + sum_{k=1..p} q_k(x_k) <= rhs
4642  * where x_k denotes a subset of the variables in x and these subsets are pairwise disjunct
4643  * and q_k(.) is a quadratic form.
4644  * p is selected as large as possible, but to be <= conshdlrdata->maxdisaggrsize.
4645  *
4646  * Without additional scaling, the constraint is disaggregated into
4647  * lhs <= b'x + sum_k c_k z_k <= rhs
4648  * c_k z_k ~ q_k(x)
4649  * where "~" is either "<=", "==", or ">=", depending on whether lhs or rhs are infinite.
4650  * Further, c_k is chosen to be the maximal absolute value of the coefficients of the quadratic terms in q_k(x).
4651  * This is done to ensure that z_k takes values with a similar magnitute as the variables in x_k (better for separation).
4652  *
4653  * However, a solution of this disaggregated system can violate the original constraint by (p+1)*epsilon
4654  * (assuming unscaled violations are used, which is the default).
4655  * Therefore, all constraints are scaled by p+1:
4656  * (p+1)*lhs <= (p+1)*b'x + (p+1) * sum_k c_k z_k <= (p+1) * rhs
4657  * (p+1)*c_k z_k ~ (p+1)*q_k(x)
4658  */
4659 static
4660 SCIP_RETCODE presolveDisaggregate(
4661  SCIP* scip, /**< SCIP data structure */
4662  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
4663  SCIP_CONS* cons, /**< source constraint to try to convert */
4664  int* naddconss /**< pointer to counter of added constraints */
4665  )
4666 {
4667  SCIP_CONSDATA* consdata;
4668  SCIP_HASHMAP* var2component;
4669  int* componentssize;
4670  int ncomponents;
4671  int i;
4672  int comp;
4673  SCIP_CONS** auxconss;
4674  SCIP_VAR** auxvars;
4675  SCIP_Real* auxcoefs;
4676 #ifdef WITH_DEBUG_SOLUTION
4677  SCIP_Real* auxsolvals; /* value of auxiliary variable in debug solution */
4678 #endif
4679  SCIP_Real scale;
4680  char name[SCIP_MAXSTRLEN];
4681 
4682  assert(scip != NULL);
4683  assert(conshdlr != NULL);
4684  assert(cons != NULL);
4685  assert(naddconss != NULL);
4686 
4687  consdata = SCIPconsGetData(cons);
4688  assert(consdata != NULL);
4689 
4690  /* skip if constraint has been already disaggregated */
4691  if( consdata->isdisaggregated )
4692  return SCIP_OKAY;
4693 
4694  consdata->isdisaggregated = TRUE;
4695 
4696  /* make sure there are no quadratic variables without coefficients */
4697  SCIP_CALL( mergeAndCleanBilinearTerms(scip, cons) );
4698  SCIP_CALL( mergeAndCleanQuadVarTerms(scip, cons) );
4699 
4700  if( consdata->nquadvars <= 1 )
4701  return SCIP_OKAY;
4702 
4703  /* sort quadratic variable terms here, so we can later search in it without reordering the array */
4704  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
4705 
4706  /* check how many quadratic terms with non-overlapping variables we have
4707  * in other words, the number of components in the sparsity graph of the quadratic term matrix
4708  */
4709  ncomponents = 0;
4710  SCIP_CALL( SCIPhashmapCreate(&var2component, SCIPblkmem(scip), consdata->nquadvars) );
4711  SCIP_CALL( SCIPallocBufferArray(scip, &componentssize, consdata->nquadvars) );
4712  for( i = 0; i < consdata->nquadvars; ++i )
4713  {
4714  /* if variable was marked already, skip it */
4715  if( SCIPhashmapExists(var2component, (void*)consdata->quadvarterms[i].var) )
4716  continue;
4717 
4718  /* start a new component with variable i */
4719  componentssize[ncomponents] = 0;
4720  SCIP_CALL( presolveDisaggregateMarkComponent(scip, consdata, i, var2component, ncomponents, componentssize + ncomponents) );
4721  ++ncomponents;
4722  }
4723 
4724  assert(ncomponents >= 1);
4725 
4726  /* if there is only one component, we cannot disaggregate
4727  * @todo we could still split the constraint into several while keeping the number of variables sharing several constraints as small as possible
4728  */
4729  if( ncomponents == 1 )
4730  {
4731  SCIPhashmapFree(&var2component);
4732  SCIPfreeBufferArray(scip, &componentssize);
4733  return SCIP_OKAY;
4734  }
4735 
4736  /* merge some components, if necessary */
4737  SCIP_CALL( presolveDisaggregateMergeComponents(scip, conshdlr, var2component, consdata->nquadvars, &ncomponents, componentssize) );
4738 
4739  SCIPfreeBufferArray(scip, &componentssize);
4740 
4741  /* scale all new constraints (ncomponents+1 many) by ncomponents+1 (or its next power of 2), so violations sum up to at most epsilon */
4742  scale = nextPowerOf2((unsigned int)ncomponents + 1);
4743 
4744  SCIP_CALL( SCIPallocBufferArray(scip, &auxconss, ncomponents) );
4745  SCIP_CALL( SCIPallocBufferArray(scip, &auxvars, ncomponents) );
4746  SCIP_CALL( SCIPallocBufferArray(scip, &auxcoefs, ncomponents) );
4747 #ifdef WITH_DEBUG_SOLUTION
4748  SCIP_CALL( SCIPallocClearBufferArray(scip, &auxsolvals, ncomponents) );
4749 #endif
4750 
4751  /* create auxiliary variables and empty constraints for each component */
4752  for( comp = 0; comp < ncomponents; ++comp )
4753  {
4754  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_comp%d", SCIPconsGetName(cons), comp);
4755 
4756  SCIP_CALL( SCIPcreateVar(scip, &auxvars[comp], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
4757  SCIP_VARTYPE_CONTINUOUS, SCIPconsIsInitial(cons), FALSE, NULL, NULL, NULL, NULL, NULL) );
4758 
4759  SCIP_CALL( SCIPcreateConsQuadratic2(scip, &auxconss[comp], name, 0, NULL, NULL, 0, NULL, 0, NULL,
4760  (SCIPisInfinity(scip, -consdata->lhs) ? -SCIPinfinity(scip) : 0.0),
4761  (SCIPisInfinity(scip, consdata->rhs) ? SCIPinfinity(scip) : 0.0),
4762  SCIPconsIsInitial(cons), SCIPconsIsSeparated(cons), TRUE,
4763  TRUE, SCIPconsIsPropagated(cons), SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons),
4764  SCIPconsIsDynamic(cons), SCIPconsIsRemovable(cons)) );
4765 
4766  auxcoefs[comp] = SCIPinfinity(scip);
4767  }
4768 
4769  /* add quadratic variables to each component constraint
4770  * delete adjacency information */
4771  for( i = 0; i < consdata->nquadvars; ++i )
4772  {
4773  assert(SCIPhashmapExists(var2component, consdata->quadvarterms[i].var));
4774 
4775  comp = SCIPhashmapGetImageInt(var2component, consdata->quadvarterms[i].var);
4776  assert(comp >= 0);
4777  assert(comp < ncomponents);
4778 
4779  /* add variable term to corresponding constraint */
4780  SCIP_CALL( SCIPaddQuadVarQuadratic(scip, auxconss[comp], consdata->quadvarterms[i].var, scale * consdata->quadvarterms[i].lincoef, scale * consdata->quadvarterms[i].sqrcoef) );
4781 
4782  /* reduce coefficient of aux variable */
4783  if( !SCIPisZero(scip, consdata->quadvarterms[i].lincoef) && ABS(consdata->quadvarterms[i].lincoef) < auxcoefs[comp] )
4784  auxcoefs[comp] = REALABS(consdata->quadvarterms[i].lincoef);
4785  if( !SCIPisZero(scip, consdata->quadvarterms[i].sqrcoef) && ABS(consdata->quadvarterms[i].sqrcoef) < auxcoefs[comp] )
4786  auxcoefs[comp] = REALABS(consdata->quadvarterms[i].sqrcoef);
4787 
4788  SCIPfreeBlockMemoryArray(scip, &consdata->quadvarterms[i].adjbilin, consdata->quadvarterms[i].adjbilinsize);
4789  consdata->quadvarterms[i].nadjbilin = 0;
4790  consdata->quadvarterms[i].adjbilinsize = 0;
4791 
4792 #ifdef WITH_DEBUG_SOLUTION
4793  if( SCIPdebugIsMainscip(scip) )
4794  {
4795  SCIP_Real debugvarval;
4796 
4797  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata->quadvarterms[i].var, &debugvarval) );
4798  auxsolvals[comp] += consdata->quadvarterms[i].lincoef * debugvarval + consdata->quadvarterms[i].sqrcoef * debugvarval * debugvarval;
4799  }
4800 #endif
4801  }
4802 
4803  /* add bilinear terms to each component constraint */
4804  for( i = 0; i < consdata->nbilinterms; ++i )
4805  {
4806  assert(SCIPhashmapExists(var2component, consdata->bilinterms[i].var1));
4807  assert(SCIPhashmapExists(var2component, consdata->bilinterms[i].var2));
4808 
4809  comp = SCIPhashmapGetImageInt(var2component, consdata->bilinterms[i].var1);
4810  assert(comp == SCIPhashmapGetImageInt(var2component, consdata->bilinterms[i].var2));
4811  assert(!SCIPisZero(scip, consdata->bilinterms[i].coef));
4812 
4813  SCIP_CALL( SCIPaddBilinTermQuadratic(scip, auxconss[comp],
4814  consdata->bilinterms[i].var1, consdata->bilinterms[i].var2, scale * consdata->bilinterms[i].coef) );
4815 
4816  if( ABS(consdata->bilinterms[i].coef) < auxcoefs[comp] )
4817  auxcoefs[comp] = ABS(consdata->bilinterms[i].coef);
4818 
4819 #ifdef WITH_DEBUG_SOLUTION
4820  if( SCIPdebugIsMainscip(scip) )
4821  {
4822  SCIP_Real debugvarval1;
4823  SCIP_Real debugvarval2;
4824 
4825  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata->bilinterms[i].var1, &debugvarval1) );
4826  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata->bilinterms[i].var2, &debugvarval2) );
4827  auxsolvals[comp] += consdata->bilinterms[i].coef * debugvarval1 * debugvarval2;
4828  }
4829 #endif
4830  }
4831 
4832  /* forget about bilinear terms in cons */
4833  SCIPfreeBlockMemoryArray(scip, &consdata->bilinterms, consdata->bilintermssize);
4834  consdata->nbilinterms = 0;
4835  consdata->bilintermssize = 0;
4836 
4837  /* remove quadratic variable terms from cons */
4838  for( i = consdata->nquadvars - 1; i >= 0; --i )
4839  {
4840  SCIP_CALL( delQuadVarTermPos(scip, cons, i) );
4841  }
4842  assert(consdata->nquadvars == 0);
4843 
4844  /* scale remaining linear variables and sides by scale */
4845  for( i = 0; i < consdata->nlinvars; ++i )
4846  {
4847  SCIP_CALL( chgLinearCoefPos(scip, cons, i, scale * consdata->lincoefs[i]) );
4848  }
4849  if( !SCIPisInfinity(scip, -consdata->lhs) )
4850  {
4851  consdata->lhs *= scale;
4852  assert(!SCIPisInfinity(scip, -consdata->lhs) );
4853  }
4854  if( !SCIPisInfinity(scip, consdata->rhs) )
4855  {
4856  consdata->rhs *= scale;
4857  assert(!SCIPisInfinity(scip, consdata->rhs) );
4858  }
4859 
4860  /* add auxiliary variables to auxiliary constraints
4861  * add aux vars and constraints to SCIP
4862  * add aux vars to this constraint
4863  * set value of aux vars in debug solution, if any
4864  */
4865  SCIPdebugMsg(scip, "add %d constraints for disaggregation of quadratic constraint <%s>\n", ncomponents, SCIPconsGetName(cons));
4866  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars + ncomponents) );
4867  for( comp = 0; comp < ncomponents; ++comp )
4868  {
4869  SCIP_CONSDATA* auxconsdata;
4870 
4871  SCIP_CALL( SCIPaddLinearVarQuadratic(scip, auxconss[comp], auxvars[comp], -scale * auxcoefs[comp]) );
4872 
4873  SCIP_CALL( SCIPaddVar(scip, auxvars[comp]) );
4874 
4875  SCIP_CALL( SCIPaddCons(scip, auxconss[comp]) );
4876  SCIPdebugPrintCons(scip, auxconss[comp], NULL);
4877 
4878  SCIP_CALL( addLinearCoef(scip, cons, auxvars[comp], scale * auxcoefs[comp]) );
4879 
4880  /* mark that the constraint should not further be disaggregated */
4881  auxconsdata = SCIPconsGetData(auxconss[comp]);
4882  assert(auxconsdata != NULL);
4883  auxconsdata->isdisaggregated = TRUE;
4884 
4885 #ifdef WITH_DEBUG_SOLUTION
4886  if( SCIPdebugIsMainscip(scip) )
4887  {
4888  /* auxvar should take value from auxsolvals in debug solution, but we also scaled auxvar by auxcoefs[comp] */
4889  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvars[comp], auxsolvals[comp] / auxcoefs[comp]) );
4890  }
4891 #endif
4892 
4893  SCIP_CALL( SCIPreleaseCons(scip, &auxconss[comp]) );
4894  SCIP_CALL( SCIPreleaseVar(scip, &auxvars[comp]) );
4895  }
4896  *naddconss += ncomponents;
4897 
4898  SCIPdebugPrintCons(scip, cons, NULL);
4899 
4900  SCIPfreeBufferArray(scip, &auxconss);
4901  SCIPfreeBufferArray(scip, &auxvars);
4902  SCIPfreeBufferArray(scip, &auxcoefs);
4903 #ifdef WITH_DEBUG_SOLUTION
4904  SCIPfreeBufferArray(scip, &auxsolvals);
4905 #endif
4906  SCIPhashmapFree(&var2component);
4907 
4908  return SCIP_OKAY;
4909 }
4910 
4911 #ifdef CHECKIMPLINBILINEAR
4912 /** checks if there are bilinear terms x*y with a binary variable x and an implication x = {0,1} -> y = 0
4913  *
4914  * In this case, the bilinear term can be removed (x=0 case) or replaced by y (x=1 case).
4915  */
4916 static
4917 SCIP_RETCODE presolveApplyImplications(
4918  SCIP* scip, /**< SCIP data structure */
4919  SCIP_CONS* cons, /**< quadratic constraint */
4920  int* nbilinremoved /**< buffer to store number of removed bilinear terms */
4921  )
4922 {
4923  SCIP_CONSDATA* consdata;
4924  SCIP_VAR* x;
4925  SCIP_VAR* y;
4926  SCIP_INTERVAL implbnds;
4927  int i;
4928  int j;
4929  int k;
4930 
4931  assert(scip != NULL);
4932  assert(cons != NULL);
4933  assert(nbilinremoved != NULL);
4934 
4935  *nbilinremoved = 0;
4936 
4937  consdata = SCIPconsGetData(cons);
4938  assert(consdata != NULL);
4939 
4940  SCIPdebugMsg(scip, "apply implications in <%s>\n", SCIPconsGetName(cons));
4941 
4942  /* sort quadvarterms in case we need to search */
4943  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
4944 
4945  for( i = 0; i < consdata->nquadvars; ++i )
4946  {
4947  x = consdata->quadvarterms[i].var;
4948  assert(x != NULL);
4949 
4950  if( consdata->quadvarterms[i].nadjbilin == 0 )
4951  continue;
4952 
4953  if( !SCIPvarIsBinary(x) )
4954  continue;
4955 
4956  if( !SCIPvarIsActive(x) )
4957  continue;
4958 
4959  if( SCIPvarGetNImpls(x, TRUE) == 0 && SCIPvarGetNImpls(x, FALSE) == 0 )
4960  continue;
4961 
4962  for( j = 0; j < consdata->quadvarterms[i].nadjbilin; ++j )
4963  {
4964  k = consdata->quadvarterms[i].adjbilin[j];
4965  assert(k >= 0);
4966  assert(k < consdata->nbilinterms);
4967 
4968  if( consdata->bilinterms[k].coef == 0.0 )
4969  continue;
4970 
4971  y = consdata->bilinterms[k].var1 == x ? consdata->bilinterms[k].var2 : consdata->bilinterms[k].var1;
4972  assert(x != y);
4973 
4974  SCIP_CALL( getImpliedBounds(scip, x, TRUE, y, &implbnds) );
4975  if( SCIPisZero(scip, implbnds.inf) && SCIPisZero(scip, implbnds.sup) )
4976  {
4977  /* if x = 1 implies y = 0, then we can remove the bilinear term x*y, since it is always 0
4978  * we only set the coefficient to 0.0 here and mark the bilinterms as not merged */
4979  SCIPdebugMsg(scip, "remove bilinear term %g<%s><%s> from <%s> due to implication\n", consdata->bilinterms[k].coef, SCIPvarGetName(x), SCIPvarGetName(y), SCIPconsGetName(cons));
4980  consdata->bilinterms[k].coef = 0.0;
4981  consdata->bilinmerged = FALSE;
4982  ++*nbilinremoved;
4983  continue;
4984  }
4985 
4986  SCIP_CALL( getImpliedBounds(scip, x, FALSE, y, &implbnds) );
4987  if( SCIPisZero(scip, implbnds.inf) && SCIPisZero(scip, implbnds.sup) )
4988  {
4989  /* if x = 0 implies y = 0, then we can replace the bilinear term x*y by y
4990  * we only move the coefficient to the linear coef of y here and mark the bilinterms as not merged */
4991  SCIPdebugMsg(scip, "replace bilinear term %g<%s><%s> by %g<%s> in <%s> due to implication\n", consdata->bilinterms[k].coef, SCIPvarGetName(x), SCIPvarGetName(y), consdata->bilinterms[k].coef, SCIPvarGetName(y), SCIPconsGetName(cons));
4992  assert(consdata->quadvarssorted);
4993  SCIP_CALL( SCIPaddQuadVarLinearCoefQuadratic(scip, cons, y, consdata->bilinterms[k].coef) );
4994  consdata->bilinterms[k].coef = 0.0;
4995  consdata->bilinmerged = FALSE;
4996  ++*nbilinremoved;
4997  }
4998  }
4999  }
5000 
5001  if( *nbilinremoved > 0 )
5002  {
5003  SCIP_CALL( mergeAndCleanBilinearTerms(scip, cons) );
5004 
5005  /* invalidate nonlinear row */
5006  if( consdata->nlrow != NULL )
5007  {
5008  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
5009  }
5010 
5011  consdata->ispropagated = FALSE;
5012  consdata->ispresolved = FALSE;
5013  consdata->iscurvchecked = FALSE;
5014  }
5015 
5016  consdata->isimpladded = FALSE;
5017 
5018  return SCIP_OKAY;
5019 }
5020 #endif
5021 
5022 /** checks a quadratic constraint for convexity and/or concavity without checking multivariate functions */
5023 static
5024 void checkCurvatureEasy(
5025  SCIP* scip, /**< SCIP data structure */
5026  SCIP_CONS* cons, /**< quadratic constraint */
5027  SCIP_HASHMAP* assumevarfixed, /**< variables to be assumed to be fixed, or NULL */
5028  SCIP_Bool* determined, /**< pointer to store whether the curvature could be determined */
5029  SCIP_Bool checkmultivariate, /**< whether curvature will be checked later on for multivariate functions */
5030  SCIP_Bool* isconvex, /**< buffer to store whether found convex */
5031  SCIP_Bool* isconcave, /**< buffer to store whether found concave */
5032  SCIP_Real* maxnonconvexity /**< buffer to store "maximal nonconvexity" */
5033  )
5034 {
5035  SCIP_CONSDATA* consdata;
5036  int nquadvars;
5037 
5038  assert(scip != NULL);
5039  assert(cons != NULL);
5040  assert(determined != NULL);
5041  assert(isconvex != NULL);
5042  assert(isconcave != NULL);
5043  assert(maxnonconvexity != NULL);
5044 
5045  consdata = SCIPconsGetData(cons);
5046  assert(consdata != NULL);
5047 
5048  nquadvars = consdata->nquadvars;
5049  *determined = TRUE;
5050 
5051  SCIPdebugMsg(scip, "Checking curvature of constraint <%s> without multivariate functions\n", SCIPconsGetName(cons));
5052 
5053  *maxnonconvexity = 0.0;
5054  if( nquadvars == 1 )
5055  {
5056  assert(consdata->nbilinterms == 0);
5057  if( assumevarfixed == NULL || !SCIPhashmapExists(assumevarfixed, (void*)consdata->quadvarterms[0].var) )
5058  {
5059  *isconvex = !SCIPisNegative(scip, consdata->quadvarterms[0].sqrcoef);
5060  *isconcave = !SCIPisPositive(scip, consdata->quadvarterms[0].sqrcoef);
5061 
5062  if( !SCIPisInfinity(scip, -consdata->lhs) && consdata->quadvarterms[0].sqrcoef > 0.0 )
5063  *maxnonconvexity = consdata->quadvarterms[0].sqrcoef;
5064  if( !SCIPisInfinity(scip, consdata->rhs) && consdata->quadvarterms[0].sqrcoef < 0.0 )
5065  *maxnonconvexity = -consdata->quadvarterms[0].sqrcoef;
5066  }
5067  else
5068  {
5069  /* only variable should be assumed to be fixed, so we are linear */
5070  *isconvex = TRUE;
5071  *isconcave = TRUE;
5072  }
5073  }
5074  else if( nquadvars == 0 )
5075  {
5076  *isconvex = TRUE;
5077  *isconcave = TRUE;
5078  }
5079  else if( consdata->nbilinterms == 0 )
5080  {
5081  int v;
5082 
5083  *isconvex = TRUE;
5084  *isconcave = TRUE;
5085 
5086  for( v = nquadvars - 1; v >= 0; --v )
5087  {
5088  /* skip assumed-to-be-fixed variables */
5089  if( assumevarfixed != NULL && SCIPhashmapExists(assumevarfixed, (void*)consdata->quadvarterms[v].var) )
5090  continue;
5091 
5092  *isconvex = *isconvex && !SCIPisNegative(scip, consdata->quadvarterms[v].sqrcoef);
5093  *isconcave = *isconcave && !SCIPisPositive(scip, consdata->quadvarterms[v].sqrcoef);
5094 
5095  if( !SCIPisInfinity(scip, -consdata->lhs) && consdata->quadvarterms[v].sqrcoef > consdata->maxnonconvexity )
5096  *maxnonconvexity = consdata->quadvarterms[0].sqrcoef;
5097  if( !SCIPisInfinity(scip, consdata->rhs) && -consdata->quadvarterms[v].sqrcoef > consdata->maxnonconvexity )
5098  *maxnonconvexity = -consdata->quadvarterms[0].sqrcoef;
5099  }
5100  }
5101  else if( !checkmultivariate )
5102  {
5103  *isconvex = FALSE;
5104  *isconcave = FALSE;
5105  *maxnonconvexity = SCIPinfinity(scip);
5106  }
5107  else
5108  {
5109  *isconvex = FALSE;
5110  *isconcave = FALSE;
5111  *determined = FALSE;
5112  }
5113 }
5114 
5115 /** checks a quadratic constraint for convexity and/or concavity while checking multivariate functions */
5116 static
5117 SCIP_RETCODE checkCurvatureExpensive(
5118  SCIP* scip, /**< SCIP data structure */
5119  SCIP_CONS* cons, /**< quadratic constraint */
5120  SCIP_HASHMAP* assumevarfixed, /**< variables to be assumed to be fixed, or NULL */
5121  SCIP_Bool* isconvex, /**< buffer to store whether found convex */
5122  SCIP_Bool* isconcave, /**< buffer to store whether found concave */
5123  SCIP_Real* maxnonconvexity /**< buffer to store "maximal nonconvexity" */
5124  )
5125 {
5126  SCIP_CONSDATA* consdata;
5127  double* matrix;
5128  SCIP_HASHMAP* var2index;
5129  int i;
5130  int n;
5131  int nn;
5132  int row;
5133  int col;
5134  double* alleigval;
5135 
5136  assert(scip != NULL);
5137  assert(cons != NULL);
5138  assert(isconvex != NULL);
5139  assert(isconcave != NULL);
5140  assert(maxnonconvexity != NULL);
5141 
5142  consdata = SCIPconsGetData(cons);
5143  assert(consdata != NULL);
5144 
5145  n = consdata->nquadvars;
5146 
5147  SCIPdebugMsg(scip, "Checking curvature of constraint <%s> with multivariate functions\n", SCIPconsGetName(cons));
5148 
5149  if( n == 2 && assumevarfixed == NULL )
5150  {
5151  SCIP_Real tracehalf;
5152  SCIP_Real discriminantroot;
5153 
5154  /* compute eigenvalues by hand */
5155  assert(consdata->nbilinterms == 1);
5156  *isconvex =
5157  consdata->quadvarterms[0].sqrcoef >= 0 &&
5158  consdata->quadvarterms[1].sqrcoef >= 0 &&
5159  4 * consdata->quadvarterms[0].sqrcoef * consdata->quadvarterms[1].sqrcoef >= consdata->bilinterms[0].coef * consdata->bilinterms[0].coef;
5160  *isconcave =
5161  consdata->quadvarterms[0].sqrcoef <= 0 &&
5162  consdata->quadvarterms[1].sqrcoef <= 0 &&
5163  4 * consdata->quadvarterms[0].sqrcoef * consdata->quadvarterms[1].sqrcoef >= consdata->bilinterms[0].coef * consdata->bilinterms[0].coef;
5164 
5165  /* store largest eigenvalue causing nonconvexity according to sides */
5166  tracehalf = (consdata->quadvarterms[0].sqrcoef + consdata->quadvarterms[1].sqrcoef) / 2.0;
5167  discriminantroot = consdata->quadvarterms[0].sqrcoef * consdata->quadvarterms[1].sqrcoef - SQR(consdata->bilinterms[0].coef / 2.0);
5168  discriminantroot = SQR(tracehalf) - discriminantroot;
5169  assert(!SCIPisNegative(scip, discriminantroot));
5170  discriminantroot = SQRT(MAX(0.0, discriminantroot));
5171 
5172  *maxnonconvexity = 0.0;
5173  if( !SCIPisInfinity(scip, -consdata->lhs) )
5174  *maxnonconvexity = MAX(*maxnonconvexity, tracehalf + discriminantroot);
5175  if( !SCIPisInfinity(scip, consdata->rhs) )
5176  *maxnonconvexity = MAX(*maxnonconvexity, discriminantroot - tracehalf);
5177 
5178  return SCIP_OKAY;
5179  }
5180 
5181  /* do not check curvature if n is too large */
5182  nn = n * n;
5183  if( nn < 0 || (unsigned) (int) nn > UINT_MAX / sizeof(SCIP_Real) )
5184  {
5185  SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "cons_quadratic - n is too large to check the curvature\n");
5186  *isconvex = FALSE;
5187  *isconcave = FALSE;
5188  *maxnonconvexity = SCIPinfinity(scip);
5189  return SCIP_OKAY;
5190  }
5191 
5192  /* lower triangular of quadratic term matrix */
5193  SCIP_CALL( SCIPallocBufferArray(scip, &matrix, nn) );
5194  BMSclearMemoryArray(matrix, nn);
5195 
5196  *isconvex = TRUE;
5197  *isconcave = TRUE;
5198  *maxnonconvexity = 0.0;
5199 
5200  SCIP_CALL( SCIPhashmapCreate(&var2index, SCIPblkmem(scip), n) );
5201  for( i = 0; i < n; ++i )
5202  {
5203  /* skip variables that we assume are fixed -> 0 row/col in coef matrix */
5204  if( assumevarfixed != NULL && SCIPhashmapExists(assumevarfixed, (void*)consdata->quadvarterms[i].var) )
5205  continue;
5206 
5207  if( consdata->quadvarterms[i].nadjbilin > 0 )
5208  {
5209  SCIP_CALL( SCIPhashmapInsertInt(var2index, consdata->quadvarterms[i].var, i) );
5210  matrix[i*n + i] = consdata->quadvarterms[i].sqrcoef;
5211  }
5212  else
5213  {
5214  /* if pure square term, then update maximal nonconvex eigenvalue, as it will not be considered in lapack call below */
5215  if( !SCIPisInfinity(scip, -consdata->lhs) && consdata->quadvarterms[i].sqrcoef > *maxnonconvexity )
5216  *maxnonconvexity = consdata->quadvarterms[i].sqrcoef;
5217  if( !SCIPisInfinity(scip, consdata->rhs) && -consdata->quadvarterms[i].sqrcoef > *maxnonconvexity )
5218  *maxnonconvexity = -consdata->quadvarterms[i].sqrcoef;
5219  }
5220  /* nonzero elements on diagonal tell a lot about convexity/concavity */
5221  if( SCIPisNegative(scip, consdata->quadvarterms[i].sqrcoef) )
5222  *isconvex = FALSE;
5223  if( SCIPisPositive(scip, consdata->quadvarterms[i].sqrcoef) )
5224  *isconcave = FALSE;
5225  }
5226 
5227  /* skip lapack call, if we know already that we are indefinite
5228  * NOTE: this will leave out updating consdata->maxnonconvexity, so that it only provides a lower bound in this case
5229  */
5230  if( !*isconvex && !*isconcave )
5231  {
5232  SCIPfreeBufferArray(scip, &matrix);
5233  SCIPhashmapFree(&var2index);
5234  /* make sure that maxnonconvexity is strictly different from zero if nonconvex
5235  * TODO one could think about doing some eigenvalue estimation here (Gershgorin)
5236  */
5237  *maxnonconvexity = MAX(1000.0, *maxnonconvexity);
5238  return SCIP_OKAY;
5239  }
5240 
5241  if( SCIPisIpoptAvailableIpopt() )
5242  {
5243  for( i = 0; i < consdata->nbilinterms; ++i )
5244  {
5245  /* skip variables that we assume are fixed -> 0 row/col in coef matrix */
5246  if( assumevarfixed != NULL && (SCIPhashmapExists(assumevarfixed, (void*)consdata->bilinterms[i].var1) || SCIPhashmapExists(assumevarfixed, (void*)consdata->bilinterms[i].var2)) )
5247  continue;
5248 
5249  assert(SCIPhashmapExists(var2index, consdata->bilinterms[i].var1));
5250  assert(SCIPhashmapExists(var2index, consdata->bilinterms[i].var2));
5251  row = SCIPhashmapGetImageInt(var2index, consdata->bilinterms[i].var1);
5252  col = SCIPhashmapGetImageInt(var2index, consdata->bilinterms[i].var2);
5253  if( row < col )
5254  matrix[row * n + col] = consdata->bilinterms[i].coef/2;
5255  else
5256  matrix[col * n + row] = consdata->bilinterms[i].coef/2;
5257  }
5258 
5259  SCIP_CALL( SCIPallocBufferArray(scip, &alleigval, n) );
5260  /* @todo Can we compute only min and max eigen value?
5261  * @todo Can we estimate the numerical error?
5262  * @todo Trying a cholesky factorization may be much faster.
5263  */
5264  if( LapackDsyev(FALSE, n, matrix, alleigval) != SCIP_OKAY )
5265  {
5266  SCIPwarningMessage(scip, "Failed to compute eigenvalues of quadratic coefficient matrix of constraint %s. Assuming matrix is indefinite.\n", SCIPconsGetName(cons));
5267  *isconvex = FALSE;
5268  *isconcave = FALSE;
5269  }
5270  else
5271  {
5272  /* deconvexification reformulates a stricly convex quadratic function in binaries such that it becomes not-strictly convex
5273  * by adding the -lambda*(x^2-x) terms for lambda the smallest eigenvalue of the matrix
5274  * the result is still a convex form "but less so" (ref. papers by Guignard et.al.), but with hopefully tighter value for the continuous relaxation
5275  */
5276 #ifdef DECONVEXIFY
5277  SCIP_Bool allbinary;
5278  printf("cons <%s>[%g,%g] spectrum = [%g,%g]\n", SCIPconsGetName(cons), consdata->lhs, consdata->rhs, alleigval[0], alleigval[n-1]);
5279 #endif
5280  *isconvex &= !SCIPisNegative(scip, alleigval[0]); /*lint !e514*/
5281  *isconcave &= !SCIPisPositive(scip, alleigval[n-1]); /*lint !e514*/
5282 #ifdef DECONVEXIFY
5283  for( i = 0; i < consdata->nquadvars; ++i )
5284  if( !SCIPvarIsBinary(consdata->quadvarterms[i].var) )
5285  break;
5286  allbinary = i == consdata->nquadvars;
5287 
5288  if( !SCIPisInfinity(scip, consdata->rhs) && alleigval[0] > 0.1 && allbinary )
5289  {
5290  printf("deconvexify cons <%s> by shifting hessian by %g\n", SCIPconsGetName(cons), alleigval[0]);
5291  for( i = 0; i < consdata->nquadvars; ++i )
5292  {
5293  consdata->quadvarterms[i].sqrcoef -= alleigval[0];
5294  consdata->quadvarterms[i].lincoef += alleigval[0];
5295  }
5296  }
5297 
5298  if( !SCIPisInfinity(scip, consdata->lhs) && alleigval[n-1] < -0.1 && allbinary )
5299  {
5300  printf("deconcavify cons <%s> by shifting hessian by %g\n", SCIPconsGetName(cons), alleigval[n-1]);
5301  for( i = 0; i < consdata->nquadvars; ++i )
5302  {
5303  consdata->quadvarterms[i].sqrcoef -= alleigval[n-1];
5304  consdata->quadvarterms[i].lincoef += alleigval[n-1];
5305  }
5306  }
5307 #endif
5308  }
5309 
5310  /* update largest eigenvalue causing nonconvexity according to sides */
5311  if( !SCIPisInfinity(scip, -consdata->lhs) )
5312  *maxnonconvexity = MAX(*maxnonconvexity, alleigval[n-1]);
5313  if( !SCIPisInfinity(scip, consdata->rhs) )
5314  *maxnonconvexity = MAX(*maxnonconvexity, -alleigval[0]);
5315 
5316  SCIPfreeBufferArray(scip, &alleigval);
5317  }
5318  else
5319  {
5320  *isconvex = FALSE;
5321  *isconcave = FALSE;
5322  *maxnonconvexity = SCIPinfinity(scip);
5323  }
5324 
5325  SCIPhashmapFree(&var2index);
5326  SCIPfreeBufferArray(scip, &matrix);
5327 
5328  return SCIP_OKAY;
5329 }
5330 
5331 
5332 /** checks a quadratic constraint for convexity and/or concavity */
5333 static
5334 SCIP_RETCODE checkCurvature(
5335  SCIP* scip, /**< SCIP data structure */
5336  SCIP_CONS* cons, /**< quadratic constraint */
5337  SCIP_Bool checkmultivariate /**< whether curvature should also be checked for multivariate functions */
5338  )
5339 {
5340  SCIP_Bool determined;
5341  SCIP_Bool isconvex;
5342  SCIP_Bool isconcave;
5343  SCIP_CONSDATA* consdata;
5344 
5345  consdata = SCIPconsGetData(cons);
5346  assert(consdata != NULL);
5347 
5348  if( consdata->iscurvchecked )
5349  return SCIP_OKAY;
5350 
5351  checkCurvatureEasy(scip, cons, NULL, &determined, checkmultivariate, &isconvex, &isconcave,
5352  &consdata->maxnonconvexity);
5353  if( !determined && checkmultivariate )
5354  {
5355  SCIP_CALL( checkCurvatureExpensive(scip, cons, NULL, &isconvex, &isconcave,
5356  &consdata->maxnonconvexity) );
5357  }
5358 
5359  consdata->isconvex = isconvex;
5360  consdata->isconcave = isconcave;
5361  consdata->iscurvchecked = TRUE;
5362 
5363  return SCIP_OKAY;
5364 }
5365 
5366 /** check whether indefinite constraint function is factorable and store corresponding coefficients */
5367 static
5368 SCIP_RETCODE checkFactorable(
5369  SCIP* scip, /**< SCIP data structure */
5370  SCIP_CONS* cons /**< constraint */
5371  )
5372 {
5373  SCIP_BILINTERM* bilinterm;
5374  SCIP_CONSDATA* consdata;
5375  SCIP_Real* a;
5376  SCIP_Real* eigvals;
5377  SCIP_Real sigma1;
5378  SCIP_Real sigma2;
5379  SCIP_Bool success;
5380  int n;
5381  int i;
5382  int idx1;
5383  int idx2;
5384  int posidx;
5385  int negidx;
5386 
5387  assert(scip != NULL);
5388  assert(cons != NULL);
5389 
5390  consdata = SCIPconsGetData(cons);
5391  assert(consdata != NULL);
5392  assert(consdata->factorleft == NULL);
5393  assert(consdata->factorright == NULL);
5394 
5395  /* we don't need this if there are no bilinear terms */
5396  if( consdata->nbilinterms == 0 )
5397  return SCIP_OKAY;
5398 
5399  /* write constraint as lhs <= linear + x'^T A x' <= rhs where x' = (x,1) and
5400  * A = ( Q b/2 )
5401  * ( b^T/2 0 )
5402  * compute an eigenvalue factorization of A and check if there are one positive and one negative eigenvalue
5403  * if so, then let sigma1^2 and -sigma2^2 be these eigenvalues and v1 and v2 be the first two rows of the inverse eigenvector matrix
5404  * thus, x'^T A x' = sigma1^2 (v1^T x')^2 - sigma2^2 (v2^T x')^2
5405  * = (sigma1 (v1^T x') - sigma2 (v2^T x')) * (sigma1 (v1^T x') + sigma2 (v2^T x'))
5406  * we then store sigma1 v1^T - sigma2 v2^T as left factor coef, and sigma1 v1^T + sigma2 v2^T as right factor coef
5407  */
5408 
5409  /* if we already know that there are only positive or only negative eigenvalues, then don't try */
5410  if( consdata->iscurvchecked && (consdata->isconvex || consdata->isconcave) )
5411  return SCIP_OKAY;
5412 
5413  n = consdata->nquadvars + 1;
5414 
5415  /* @todo handle case n=3 explicitly */
5416 
5417  /* skip too large matrices */
5418  if( n > 50 )
5419  return SCIP_OKAY;
5420 
5421  /* need routine to compute eigenvalues/eigenvectors */
5422  if( !SCIPisIpoptAvailableIpopt() )
5423  return SCIP_OKAY;
5424 
5425  SCIP_CALL( consdataSortQuadVarTerms(scip, consdata) );
5426 
5427  SCIP_CALL( SCIPallocBufferArray(scip, &a, n*n) );
5428  BMSclearMemoryArray(a, n*n);
5429 
5430  /* set lower triangular entries of A corresponding to bilinear terms */
5431  for( i = 0; i < consdata->nbilinterms; ++i )
5432  {
5433  bilinterm = &consdata->bilinterms[i];
5434 
5435  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var1, &idx1) );
5436  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var2, &idx2) );
5437  assert(idx1 >= 0);
5438  assert(idx2 >= 0);
5439  assert(idx1 != idx2);
5440 
5441  a[MIN(idx1,idx2) * n + MAX(idx1,idx2)] = bilinterm->coef / 2.0;
5442  }
5443 
5444  /* set lower triangular entries of A corresponding to square and linear terms */
5445  for( i = 0; i < consdata->nquadvars; ++i )
5446  {
5447  a[i*n + i] = consdata->quadvarterms[i].sqrcoef;
5448  a[i*n + n-1] = consdata->quadvarterms[i].lincoef / 2.0;
5449  }
5450 
5451  SCIP_CALL( SCIPallocBufferArray(scip, &eigvals, n) );
5452  if( LapackDsyev(TRUE, n, a, eigvals) != SCIP_OKAY )
5453  {
5454  SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors of augmented quadratic form matrix for constraint <%s>.\n", SCIPconsGetName(cons));
5455  goto CLEANUP;
5456  }
5457 
5458  /* check if there is exactly one positive and one negative eigenvalue */
5459  posidx = -1;
5460  negidx = -1;
5461  for( i = 0; i < n; ++i )
5462  {
5463  if( SCIPisPositive(scip, eigvals[i]) )
5464  {
5465  if( posidx == -1 )
5466  posidx = i;
5467  else
5468  break;
5469  }
5470  else if( SCIPisNegative(scip, eigvals[i]) )
5471  {
5472  if( negidx == -1 )
5473  negidx = i;
5474  else
5475  break;
5476  }
5477  }
5478  if( i < n || posidx == -1 || negidx == -1 )
5479  {
5480  SCIPdebugMsg(scip, "Augmented quadratic form of constraint <%s> is not factorable.\n", SCIPconsGetName(cons));
5481  goto CLEANUP;
5482  }
5483  assert(SCIPisPositive(scip, eigvals[posidx]));
5484  assert(SCIPisNegative(scip, eigvals[negidx]));
5485 
5486  /* compute factorleft and factorright */
5487  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->factorleft, consdata->nquadvars + 1) );
5488  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->factorright, consdata->nquadvars + 1) );
5489 
5490  /* eigenvectors are stored in a, inverse eigenvector matrix is transposed of a
5491  * it seems that v1 and v2 are at &a[posidx*n] and &a[negidx*n]
5492  */
5493  sigma1 = sqrt( eigvals[posidx]);
5494  sigma2 = sqrt(-eigvals[negidx]);
5495  for( i = 0; i < n; ++i )
5496  {
5497  consdata->factorleft[i] = sigma1 * a[posidx * n + i] - sigma2 * a[negidx * n + i];
5498  consdata->factorright[i] = sigma1 * a[posidx * n + i] + sigma2 * a[negidx * n + i];
5499  /* set almost-zero elements to zero */
5500  if( SCIPisZero(scip, consdata->factorleft[i]) )
5501  consdata->factorleft[i] = 0.0;
5502  if( SCIPisZero(scip, consdata->factorright[i]) )
5503  consdata->factorright[i] = 0.0;
5504  }
5505 
5506 #ifdef SCIP_DEBUG
5507  SCIPdebugMsg(scip, "constraint <%s> has factorable quadratic form: (%g", SCIPconsGetName(cons), consdata->factorleft[n-1]);
5508  for( i = 0; i < consdata->nquadvars; ++i )
5509  {
5510  if( consdata->factorleft[i] != 0.0 )
5511  SCIPdebugMsgPrint(scip, " %+g<%s>", consdata->factorleft[i], SCIPvarGetName(consdata->quadvarterms[i].var));
5512  }
5513  SCIPdebugMsgPrint(scip, ") * (%g", consdata->factorright[n-1]);
5514  for( i = 0; i < consdata->nquadvars; ++i )
5515  {
5516  if( consdata->factorright[i] != 0.0 )
5517  SCIPdebugMsgPrint(scip, " %+g<%s>", consdata->factorright[i], SCIPvarGetName(consdata->quadvarterms[i].var));
5518  }
5519  SCIPdebugMsgPrint(scip, ")\n");
5520 #endif
5521 
5522  /* check whether factorleft * factorright^T is matrix of augmented quadratic form
5523  * we check here only the nonzero entries from the quadratic form
5524  */
5525  success = TRUE;
5526 
5527  /* check bilinear terms */
5528  for( i = 0; i < consdata->nbilinterms; ++i )
5529  {
5530  bilinterm = &consdata->bilinterms[i];
5531 
5532  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var1, &idx1) );
5533  SCIP_CALL( consdataFindQuadVarTerm(scip, consdata, bilinterm->var2, &idx2) );
5534 
5535  if( !SCIPisRelEQ(scip, consdata->factorleft[idx1] * consdata->factorright[idx2] + consdata->factorleft[idx2] * consdata->factorright[idx1], bilinterm->coef) )
5536  {
5537  success = FALSE;
5538  break;
5539  }
5540  }
5541 
5542  /* set lower triangular entries of A corresponding to square and linear terms */
5543  for( i = 0; i < consdata->nquadvars; ++i )
5544  {
5545  if( !SCIPisRelEQ(scip, consdata->factorleft[i] * consdata->factorright[i], consdata->quadvarterms[i].sqrcoef) )
5546  {
5547  success = FALSE;
5548  break;
5549  }
5550 
5551  if( !SCIPisRelEQ(scip, consdata->factorleft[n-1] * consdata->factorright[i] + consdata->factorleft[i] * consdata->factorright[n-1], consdata->quadvarterms[i].lincoef) )
5552  {
5553  success = FALSE;
5554  break;
5555  }
5556  }
5557 
5558  if( !success )
5559  {
5560  SCIPdebugMsg(scip, "Factorization not accurate enough. Dropping it.\n");
5561  SCIPfreeBlockMemoryArray(scip, &consdata->factorleft, consdata->nquadvars + 1);
5562  SCIPfreeBlockMemoryArray(scip, &consdata->factorright, consdata->nquadvars + 1);
5563  }
5564 
5565  CLEANUP:
5566  SCIPfreeBufferArray(scip, &eigvals);
5567  SCIPfreeBufferArray(scip, &a);
5568 
5569  return SCIP_OKAY;
5570 }
5571 
5572 /** computes activity and violation of a constraint
5573  *
5574  * If solution violates bounds by more than feastol, the violation is still computed, but *solviolbounds is set to TRUE
5575  */
5576 static
5577 SCIP_RETCODE computeViolation(
5578  SCIP* scip, /**< SCIP data structure */
5579  SCIP_CONS* cons, /**< constraint */
5580  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
5581  SCIP_Bool* solviolbounds /**< buffer to store whether quadratic variables in solution are outside their bounds by more than feastol */
5582  )
5583 { /*lint --e{666}*/
5584  SCIP_CONSDATA* consdata;
5585  SCIP_Real varval;
5586  SCIP_Real varval2;
5587  SCIP_Real absviol;
5588  SCIP_Real relviol;
5589  SCIP_VAR* var;
5590  SCIP_VAR* var2;
5591  int i;
5592  int j;
5593 
5594  assert(scip != NULL);
5595  assert(cons != NULL);
5596  assert(solviolbounds != NULL);
5597 
5598  consdata = SCIPconsGetData(cons);
5599  assert(consdata != NULL);
5600 
5601  *solviolbounds = FALSE;
5602  consdata->activity = 0.0;
5603  consdata->lhsviol = 0.0;
5604  consdata->rhsviol = 0.0;
5605 
5606  for( i = 0; i < consdata->nlinvars; ++i )
5607  {
5608  SCIP_Real activity;
5609 
5610  var = consdata->linvars[i];
5611  varval = SCIPgetSolVal(scip, sol, var);
5612  activity = consdata->lincoefs[i] * varval;
5613 
5614  /* the contribution of a variable with |varval| = +inf is +inf when activity > 0.0, -inf when activity < 0.0, and
5615  * 0.0 otherwise
5616  */
5617  if( SCIPisInfinity(scip, REALABS(varval)) )
5618  {
5619  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
5620  {
5621  consdata->activity = SCIPinfinity(scip);
5622  consdata->rhsviol = SCIPinfinity(scip);
5623  return SCIP_OKAY;
5624  }
5625 
5626  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
5627  {
5628  consdata->activity = -SCIPinfinity(scip);
5629  consdata->lhsviol = SCIPinfinity(scip);
5630  return SCIP_OKAY;
5631  }
5632  }
5633 
5634  consdata->activity += activity;
5635  }
5636 
5637  for( j = 0; j < consdata->nquadvars; ++j )
5638  {
5639  SCIP_Real activity;
5640 
5641  var = consdata->quadvarterms[j].var;
5642  varval = SCIPgetSolVal(scip, sol, var);
5643  activity = (consdata->quadvarterms[j].lincoef + consdata->quadvarterms[j].sqrcoef * varval) * varval;
5644 
5645  /* the contribution of a variable with |varval| = +inf is +inf when activity > 0.0, -inf when activity < 0.0, and
5646  * 0.0 otherwise
5647  */
5648  if( SCIPisInfinity(scip, REALABS(varval)) )
5649  {
5650  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
5651  {
5652  consdata->activity = SCIPinfinity(scip);
5653  consdata->rhsviol = SCIPinfinity(scip);
5654  return SCIP_OKAY;
5655  }
5656 
5657  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
5658  {
5659  consdata->activity = -SCIPinfinity(scip);
5660  consdata->lhsviol = SCIPinfinity(scip);
5661  return SCIP_OKAY;
5662  }
5663  }
5664 
5665  /* project onto local box, in case the LP solution is slightly outside the bounds (which is not our job to enforce) */
5666  if( sol == NULL )
5667  {
5668  /* with non-initial columns, variables can shortly be a column variable before entering the LP and have value 0.0 in this case, which might violated the variable bounds */
5669  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) && !SCIPisFeasGE(scip, varval, SCIPvarGetLbLocal(var))) ||
5670  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) && !SCIPisFeasLE(scip, varval, SCIPvarGetUbLocal(var))) )
5671  *solviolbounds = TRUE;
5672  else
5673  {
5674  varval = MAX(SCIPvarGetLbLocal(var), MIN(SCIPvarGetUbLocal(var), varval));
5675  activity = (consdata->quadvarterms[j].lincoef + consdata->quadvarterms[j].sqrcoef * varval) * varval;
5676  }
5677  }
5678 
5679  consdata->activity += activity;
5680  }
5681 
5682  for( j = 0; j < consdata->nbilinterms; ++j )
5683  {
5684  SCIP_Real activity;
5685 
5686  var = consdata->bilinterms[j].var1;
5687  var2 = consdata->bilinterms[j].var2;
5688  varval = SCIPgetSolVal(scip, sol, var);
5689  varval2 = SCIPgetSolVal(scip, sol, var2);
5690 
5691  /* project onto local box, in case the LP solution is slightly outside the bounds (which is not our job to enforce) */
5692  if( sol == NULL )
5693  {
5694  /* with non-initial columns, variables can shortly be a column variable before entering the LP and have value 0.0 in this case, which might violated the variable bounds */
5695  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) && !SCIPisFeasGE(scip, varval, SCIPvarGetLbLocal(var))) ||
5696  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) && !SCIPisFeasLE(scip, varval, SCIPvarGetUbLocal(var))) )
5697  *solviolbounds = TRUE;
5698  else
5699  varval = MAX(SCIPvarGetLbLocal(var), MIN(SCIPvarGetUbLocal(var), varval));
5700 
5701  /* with non-initial columns, variables can shortly be a column variable before entering the LP and have value 0.0 in this case, which might violated the variable bounds */
5702  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var2)) && !SCIPisFeasGE(scip, varval2, SCIPvarGetLbLocal(var2))) ||
5703  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var2)) && !SCIPisFeasLE(scip, varval2, SCIPvarGetUbLocal(var2))) )
5704  *solviolbounds = TRUE;
5705  else
5706  varval2 = MAX(SCIPvarGetLbLocal(var2), MIN(SCIPvarGetUbLocal(var2), varval2));
5707  }
5708 
5709  activity = consdata->bilinterms[j].coef * varval * varval2;
5710 
5711  /* consider var*var2 as a new variable and handle it as it would appear linearly */
5712  if( SCIPisInfinity(scip, REALABS(varval*varval2)) )
5713  {
5714  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
5715  {
5716  consdata->activity = SCIPinfinity(scip);
5717  consdata->rhsviol = SCIPinfinity(scip);
5718  return SCIP_OKAY;
5719  }
5720 
5721  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
5722  {
5723  consdata->activity = -SCIPinfinity(scip);
5724  consdata->lhsviol = SCIPinfinity(scip);
5725  return SCIP_OKAY;
5726  }
5727  }
5728 
5729  consdata->activity += activity;
5730  }
5731 
5732  absviol = 0.0;
5733  relviol = 0.0;
5734  /* compute absolute violation left hand side */
5735  if( consdata->activity < consdata->lhs && !SCIPisInfinity(scip, -consdata->lhs) )
5736  {
5737  consdata->lhsviol = consdata->lhs - consdata->activity;
5738  absviol = consdata->lhsviol;
5739  relviol = SCIPrelDiff(consdata->lhs, consdata->activity);
5740  }
5741  else
5742  consdata->lhsviol = 0.0;
5743 
5744  /* compute absolute violation right hand side */
5745  if( consdata->activity > consdata->rhs && !SCIPisInfinity(scip, consdata->rhs) )
5746  {
5747  consdata->rhsviol = consdata->activity - consdata->rhs;
5748  absviol = consdata->rhsviol;
5749  relviol = SCIPrelDiff(consdata->activity, consdata->rhs);
5750  }
5751  else
5752  consdata->rhsviol = 0.0;
5753 
5754  /* update absolute and relative violation of the solution */
5755  if( sol != NULL )
5756  SCIPupdateSolConsViolation(scip, sol, absviol, relviol);
5757 
5758  return SCIP_OKAY;
5759 }
5760 
5761 /** computes violation of a set of constraints */
5762 static
5763 SCIP_RETCODE computeViolations(
5764  SCIP* scip, /**< SCIP data structure */
5765  SCIP_CONS** conss, /**< constraints */
5766  int nconss, /**< number of constraints */
5767  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
5768  SCIP_Bool* solviolbounds, /**< buffer to store whether quadratic variables in solution are outside their bounds by more than feastol in some constraint */
5769  SCIP_CONS** maxviolcon /**< buffer to store constraint with largest violation, or NULL if solution is feasible */
5770  )
5771 {
5772  SCIP_CONSDATA* consdata;
5773  SCIP_Real viol;
5774  SCIP_Real maxviol;
5775  SCIP_Bool solviolbounds1;
5776  int c;
5777 
5778  assert(scip != NULL);
5779  assert(conss != NULL || nconss == 0);
5780  assert(solviolbounds != NULL);
5781  assert(maxviolcon != NULL);
5782 
5783  *solviolbounds = FALSE;
5784  *maxviolcon = NULL;
5785 
5786  maxviol = 0.0;
5787 
5788  for( c = 0; c < nconss; ++c )
5789  {
5790  assert(conss != NULL);
5791  assert(conss[c] != NULL);
5792 
5793  SCIP_CALL( computeViolation(scip, conss[c], sol, &solviolbounds1) );
5794  *solviolbounds |= solviolbounds1;
5795 
5796  consdata = SCIPconsGetData(conss[c]);
5797  assert(consdata != NULL);
5798 
5799  viol = MAX(consdata->lhsviol, consdata->rhsviol);
5800  if( viol > maxviol && SCIPisGT(scip, viol, SCIPfeastol(scip)) )
5801  {
5802  maxviol = viol;
5803  *maxviolcon = conss[c];
5804  }
5805  }
5806 
5807  return SCIP_OKAY;
5808 }
5809 
5810 
5811 /** index comparison method for bilinear terms */
5812 static
5813 SCIP_DECL_SORTINDCOMP(bilinTermComp2)
5814 { /*lint --e{715}*/
5815  SCIP_BILINTERM* bilinterms = (SCIP_BILINTERM*)dataptr;
5816  int var1cmp;
5817 
5818  assert(bilinterms != NULL);
5819 
5820  var1cmp = SCIPvarCompare(bilinterms[ind1].var1, bilinterms[ind2].var1);
5821  if( var1cmp != 0 )
5822  return var1cmp;
5823 
5824  return SCIPvarCompare(bilinterms[ind1].var2, bilinterms[ind2].var2);
5825 }
5826 
5827 /** volume comparison method for bilinear terms; prioritizes bilinear products with a larger volume */
5828 static
5829 SCIP_DECL_SORTINDCOMP(bilinTermCompVolume)
5830 { /*lint --e{715}*/
5831  SCIP_BILINTERM* bilinterms = (SCIP_BILINTERM*)dataptr;
5832  SCIP_Real vol1;
5833  SCIP_Real vol2;
5834 
5835  assert(bilinterms != NULL);
5836 
5837  vol1 = (SCIPvarGetUbLocal(bilinterms[ind1].var1) - SCIPvarGetLbLocal(bilinterms[ind1].var1))
5838  * (