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