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