Scippy

SCIP

Solving Constraint Integer Programs

cons_abspower.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2017 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_abspower.c
17  * @brief Constraint handler for absolute power constraints \f$\textrm{lhs} \leq \textrm{sign}(x+a) |x+a|^n + c z \leq \textrm{rhs}\f$
18  * @author Stefan Vigerske
19  */
20 
21 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
22 
23 #include <assert.h>
24 #include <string.h>
25 #include <ctype.h>
26 
27 #include "scip/cons_abspower.h"
28 #include "scip/cons_nonlinear.h"
29 #include "scip/cons_indicator.h"
30 #include "scip/cons_quadratic.h"
31 #include "scip/cons_linear.h"
32 #include "scip/cons_varbound.h"
33 #include "scip/intervalarith.h"
34 #include "scip/heur_subnlp.h"
35 #include "scip/heur_trysol.h"
36 #include "scip/debug.h"
37 
38 /* constraint handler properties */
39 #define CONSHDLR_NAME "abspower"
40 #define CONSHDLR_DESC "constraint handler for absolute power constraints lhs <= sign(x+offset)abs(x+offset)^n + c*z <= rhs"
41 #define CONSHDLR_SEPAPRIORITY 0 /**< priority of the constraint handler for separation */
42 #define CONSHDLR_ENFOPRIORITY -30 /**< priority of the constraint handler for constraint enforcing */
43 #define CONSHDLR_CHECKPRIORITY -3500000 /**< priority of the constraint handler for checking feasibility */
44 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
45 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
46 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
47  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
48 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
49 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
50 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
51 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
52 
53 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_FAST | SCIP_PRESOLTIMING_MEDIUM
54 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_ALWAYS /**< when should the constraint handlers propagation routines be called? */
55 
56 #define QUADCONSUPGD_PRIORITY 50000 /**< priority of the constraint handler for upgrading of quadratic constraints */
57 #define NONLINCONSUPGD_PRIORITY 50000 /**< priority of the constraint handler for upgrading of nonlinear constraints and reformulating expression graph nodes */
58 
59 /*
60  * Local defines
61  */
62 
63 #define PROPVARTOL SCIPepsilon(scip) /**< tolerance to add to variable bounds in domain propagation */
64 #define PROPSIDETOL SCIPepsilon(scip) /**< tolerance to add to constraint sides in domain propagation */
65 #define INITLPMAXVARVAL 1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */
66 
67 /** power function type to be used by a constraint instead of the general pow */
68 #define DECL_MYPOW(x) SCIP_Real x (SCIP_Real base, SCIP_Real exponent)
69 
70 /** sign of a value (-1 or +1)
71  *
72  * 0.0 has sign +1
73  */
74 #define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
75 
76 
77 /*
78  * Data structures
79  */
80 
81 #define ROOTS_KNOWN 10 /**< up to which (integer) exponents precomputed roots have been stored */
82 
83 /** The positive root of the polynomial (n-1) y^n + n y^(n-1) - 1 is needed in separation.
84  * Here we store these roots for small integer values of n.
85  */
86 static
88  -1.0, /* no root for n=0 */
89  -1.0, /* no root for n=1 */
90  0.41421356237309504880, /* root for n=2 (-1+sqrt(2)) */
91  0.5, /* root for n=3 */
92  0.56042566045031785945, /* root for n=4 */
93  0.60582958618826802099, /* root for n=5 */
94  0.64146546982884663257, /* root for n=6 */
95  0.67033204760309682774, /* root for n=7 */
96  0.69428385661425826738, /* root for n=8 */
97  0.71453772716733489700, /* root for n=9 */
98  0.73192937842370733350 /* root for n=10 */
99 };
100 
101 /** constraint data for absolute power constraints */
102 struct SCIP_ConsData
103 {
104  SCIP_VAR* x; /**< variable x in sign(x+offset)|x+offset|^n term */
105  SCIP_VAR* z; /**< linear variable in constraint */
106  SCIP_Real exponent; /**< exponent n of |x+offset| */
107  SCIP_Real xoffset; /**< offset in x+offset */
108  SCIP_Real zcoef; /**< coefficient of linear variable z */
109  SCIP_Real lhs; /**< left hand side of constraint */
110  SCIP_Real rhs; /**< right hand side of constraint */
111 
112  SCIP_Real root; /**< root of polynomial */
113  DECL_MYPOW ((*power)); /**< function for computing power*/
114 
115  SCIP_Real lhsviol; /**< current (scaled) violation of left hand side */
116  SCIP_Real rhsviol; /**< current (scaled) violation of right hand side */
117 
118  int xeventfilterpos; /**< position of x var event in SCIP event filter */
119  int zeventfilterpos; /**< position of z var event in SCIP event filter */
120  unsigned int propvarbounds:1; /**< have variable bounds been propagated? */
121 
122  SCIP_NLROW* nlrow; /**< nonlinear row representation of constraint */
123 };
124 
125 /** constraint handler data */
126 struct SCIP_ConshdlrData
127 {
128  SCIP_Real mincutefficacysepa; /**< minimal efficacy of a cut in order to add it to relaxation during separation */
129  SCIP_Real mincutefficacyenfofac;/**< minimal target efficacy of a cut in order to add it to relaxation during enforcement as factor of feasibility tolerance (may be ignored) */
130  char scaling; /**< scaling method of constraints in feasibility check */
131  SCIP_Real cutmaxrange; /**< maximal coef range (maximal abs coef / minimal abs coef) of a cut in order to be added to LP */
132  SCIP_Bool projectrefpoint; /**< whether to project the reference point when linearizing a absolute power constraint in a convex region */
133  int preferzerobranch; /**< how much we prefer to branch on 0.0 first */
134  SCIP_Bool branchminconverror; /**< whether to compute branching point such that the convexification error is minimized after branching on 0.0 */
135  SCIP_Bool addvarboundcons; /**< should variable bound constraints be added? */
136  SCIP_Bool linfeasshift; /**< try linear feasibility shift heuristic in CONSCHECK */
137  SCIP_Bool dualpresolve; /**< should dual presolve be applied? */
138  SCIP_Bool sepainboundsonly; /**< should tangents only be generated in variable bounds during separation? */
139  SCIP_Real sepanlpmincont; /**< minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation */
140  SCIP_Bool enfocutsremovable; /**< are cuts added during enforcement removable from the LP in the same node? */
141 
142  SCIP_HEUR* subnlpheur; /**< a pointer to the subnlp heuristic */
143  SCIP_HEUR* trysolheur; /**< a pointer to the trysol heuristic */
144  SCIP_EVENTHDLR* eventhdlr; /**< our handler for bound change events on variable x */
145  SCIP_CONSHDLR* conshdlrindicator; /**< a pointer to the indicator constraint handler */
146  int newsoleventfilterpos;/**< filter position of new solution event handler, if catched */
147  SCIP_Bool comparedpairwise; /**< did we compare absolute power constraints pairwise in this run? */
148  SCIP_Bool sepanlp; /**< where linearization of the NLP relaxation solution added? */
149  SCIP_NODE* lastenfonode; /**< the node for which enforcement was called the last time (and some constraint was violated) */
150  int nenforounds; /**< counter on number of enforcement rounds for the current node */
151  unsigned int nsecantcuts; /**< number of secant cuts created so far */
152  unsigned int ncuts; /**< number of linearization cuts created so far */
153 };
154 
155 /*
156  * Propagation rules
157  */
158 
159 enum Proprule
160 {
161  PROPRULE_1, /**< left hand side and bounds on z -> lower bound on x */
162  PROPRULE_2, /**< left hand side and upper bound on x -> bound on z */
163  PROPRULE_3, /**< right hand side and bounds on z -> upper bound on x */
164  PROPRULE_4, /**< right hand side and lower bound on x -> bound on z */
165  PROPRULE_INVALID /**< propagation was applied without a specific propagation rule */
166 };
167 typedef enum Proprule PROPRULE;
169 /*
170  * Local methods
171  */
172 
173 /** power function for square, that should be faster than using pow(x, 2.0) */
174 static
176 {
177  assert(exponent == 2.0);
178  return base*base;
179 }
180 
181 /** process variable event */
182 static
183 SCIP_DECL_EVENTEXEC(processVarEvent)
184 {
185  SCIP_CONS* cons;
186 
187  assert(scip != NULL);
188  assert(event != NULL);
190 
191  cons = (SCIP_CONS*) eventdata;
192  assert(cons != NULL);
193 
194  assert(SCIPconsGetData(cons) != NULL);
195  assert(SCIPeventGetVar(event) == SCIPconsGetData(cons)->x || SCIPeventGetVar(event) == SCIPconsGetData(cons)->z);
196 
198 
199  return SCIP_OKAY;
200 } /*lint !e715*/
201 
202 /** catch variable bound tightening events */
203 static
205  SCIP* scip, /**< SCIP data structure */
206  SCIP_EVENTHDLR* eventhdlr, /**< event handler for variables */
207  SCIP_CONS* cons /**< constraint for which to catch bound change events */
208  )
209 {
210  SCIP_CONSDATA* consdata;
211  SCIP_EVENTTYPE eventtype;
212 
213  assert(scip != NULL);
214  assert(cons != NULL);
215  assert(eventhdlr != NULL);
216 
217  consdata = SCIPconsGetData(cons);
218  assert(consdata != NULL);
219 
220  /* if z is multiaggregated, then bound changes on x could not be propagated, so we do not need to catch them */
221  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR )
222  {
223  eventtype = SCIP_EVENTTYPE_DISABLED;
224  if( !SCIPisInfinity(scip, -consdata->lhs) )
225  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
226  if( !SCIPisInfinity(scip, consdata->rhs) )
227  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
228 
229  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->x, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, &consdata->xeventfilterpos) );
230 
231  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
232  }
233 
234  /* if x is multiaggregated, then bound changes on z could not be propagated, so we do not need to catch them */
235  if( SCIPvarGetStatus(consdata->x) != SCIP_VARSTATUS_MULTAGGR )
236  {
237  eventtype = SCIP_EVENTTYPE_DISABLED;
238  if( consdata->zcoef > 0.0 )
239  {
240  if( !SCIPisInfinity(scip, -consdata->lhs) )
241  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
242  if( !SCIPisInfinity(scip, consdata->rhs) )
243  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
244  }
245  else
246  {
247  if( !SCIPisInfinity(scip, -consdata->lhs) )
248  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
249  if( !SCIPisInfinity(scip, consdata->rhs) )
250  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
251  }
252 
253  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->z, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, &consdata->zeventfilterpos) );
254 
255  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
256  }
257 
258  return SCIP_OKAY;
259 }
260 
261 /** drop variable bound tightening events */
262 static
264  SCIP* scip, /**< SCIP data structure */
265  SCIP_EVENTHDLR* eventhdlr, /**< event handler for variables */
266  SCIP_CONS* cons /**< constraint for which to drop bound change events */
267  )
268 {
269  SCIP_CONSDATA* consdata;
270  SCIP_EVENTTYPE eventtype;
271 
272  assert(scip != NULL);
273  assert(cons != NULL);
274  assert(eventhdlr != NULL);
275 
276  consdata = SCIPconsGetData(cons);
277  assert(consdata != NULL);
278 
279  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR )
280  {
281  eventtype = SCIP_EVENTTYPE_DISABLED;
282  if( !SCIPisInfinity(scip, -consdata->lhs) )
283  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
284  if( !SCIPisInfinity(scip, consdata->rhs) )
285  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
286 
287  SCIP_CALL( SCIPdropVarEvent(scip, consdata->x, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, consdata->xeventfilterpos) );
288  consdata->xeventfilterpos = -1;
289  }
290 
291  if( SCIPvarGetStatus(consdata->x) != SCIP_VARSTATUS_MULTAGGR )
292  {
293  eventtype = SCIP_EVENTTYPE_DISABLED;
294  if( consdata->zcoef > 0.0 )
295  {
296  if( !SCIPisInfinity(scip, -consdata->lhs) )
297  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
298  if( !SCIPisInfinity(scip, consdata->rhs) )
299  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
300  }
301  else
302  {
303  if( !SCIPisInfinity(scip, -consdata->lhs) )
304  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
305  if( !SCIPisInfinity(scip, consdata->rhs) )
306  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
307  }
308 
309  SCIP_CALL( SCIPdropVarEvent(scip, consdata->z, eventtype, eventhdlr, (SCIP_EVENTDATA*)cons, consdata->zeventfilterpos) );
310  consdata->zeventfilterpos = -1;
311  }
312 
313  assert(consdata->xeventfilterpos == -1);
314  assert(consdata->zeventfilterpos == -1);
315 
316  return SCIP_OKAY;
317 }
318 
319 /** get key of hash element */
320 static
321 SCIP_DECL_HASHGETKEY(presolveFindDuplicatesGetKey)
322 {
323  return elem;
324 } /*lint !e715*/
325 
326 /** checks if two constraints have the same x variable, the same exponent, and either the same offset or the same linear variable and are both equality constraint */
327 static
328 SCIP_DECL_HASHKEYEQ(presolveFindDuplicatesKeyEQ)
329 {
330  SCIP_CONSDATA* consdata1;
331  SCIP_CONSDATA* consdata2;
332 
333  consdata1 = SCIPconsGetData((SCIP_CONS*)key1);
334  consdata2 = SCIPconsGetData((SCIP_CONS*)key2);
335  assert(consdata1 != NULL);
336  assert(consdata2 != NULL);
337 
338  if( consdata1->x != consdata2->x )
339  return FALSE;
340 
341  if( consdata1->exponent != consdata2->exponent ) /*lint !e777*/
342  return FALSE;
343 
344  if( consdata1->xoffset != consdata2->xoffset && consdata1->z != consdata2->z ) /*lint !e777*/
345  return FALSE;
346 
347  return TRUE;
348 } /*lint !e715*/
349 
350 /** get value of hash element when comparing on x */
351 static
352 SCIP_DECL_HASHKEYVAL(presolveFindDuplicatesKeyVal)
353 {
354  SCIP_CONSDATA* consdata;
355 
356  consdata = SCIPconsGetData((SCIP_CONS*)key);
357  assert(consdata != NULL);
358 
359  return SCIPhashTwo(SCIPvarGetIndex(consdata->x),
360  SCIPpositiveRealHashCode(consdata->exponent, 7));
361 } /*lint !e715*/
362 
363 /** checks if two constraints have the same z variable and the same exponent */
364 static
365 SCIP_DECL_HASHKEYEQ(presolveFindDuplicatesKeyEQ2)
366 {
367  SCIP_CONSDATA* consdata1;
368  SCIP_CONSDATA* consdata2;
369 
370  consdata1 = SCIPconsGetData((SCIP_CONS*)key1);
371  consdata2 = SCIPconsGetData((SCIP_CONS*)key2);
372  assert(consdata1 != NULL);
373  assert(consdata2 != NULL);
374 
375  if( consdata1->z != consdata2->z )
376  return FALSE;
377 
378  if( consdata1->exponent != consdata2->exponent ) /*lint !e777*/
379  return FALSE;
380 
381  return TRUE;
382 } /*lint !e715*/
383 
384 /** get value of hash element when comparing on z */
385 static
386 SCIP_DECL_HASHKEYVAL(presolveFindDuplicatesKeyVal2)
387 {
388  SCIP_CONSDATA* consdata;
389 
390  consdata = SCIPconsGetData((SCIP_CONS*)key);
391  assert(consdata != NULL);
392 
393  return SCIPhashTwo(SCIPvarGetIndex(consdata->z),
394  SCIPpositiveRealHashCode(consdata->exponent, 7));
395 } /*lint !e715*/
396 
397 /** upgrades a signpower constraint to a linear constraint if a second signpower constraint with same nonlinear term is available */
398 static
400  SCIP* scip, /**< SCIP data structure */
401  SCIP_CONS* cons1, /**< constraint to upgrade to a linear constraint */
402  SCIP_CONS* cons2, /**< constraint which defines a relation for x|x|^{n-1} */
403  SCIP_Bool* infeas, /**< buffer where to indicate if infeasibility has been detected */
404  int* nupgdconss, /**< buffer where to add number of upgraded conss */
405  int* ndelconss, /**< buffer where to add number of deleted conss */
406  int* naggrvars /**< buffer where to add number of aggregated variables */
407  )
408 {
409  SCIP_CONSDATA* consdata1;
410  SCIP_CONSDATA* consdata2;
411  SCIP_CONS* lincons;
412  SCIP_Real lhs;
413  SCIP_Real rhs;
414  SCIP_VAR* vars[2];
415  SCIP_Real coefs[2];
416 
417  assert(scip != NULL);
418  assert(cons1 != NULL);
419  assert(cons2 != NULL);
420  assert(infeas != NULL);
421  assert(nupgdconss != NULL);
422  assert(ndelconss != NULL);
423  assert(naggrvars != NULL);
424 
425  consdata1 = SCIPconsGetData(cons1);
426  consdata2 = SCIPconsGetData(cons2);
427  assert(consdata1 != NULL);
428  assert(consdata2 != NULL);
429 
430  assert(SCIPisEQ(scip, consdata2->lhs, consdata2->rhs));
431  assert(!SCIPisInfinity(scip, consdata2->lhs));
432  assert(consdata1->x == consdata2->x);
433  assert(consdata1->exponent == consdata2->exponent); /*lint !e777*/
434  assert(consdata1->xoffset == consdata2->xoffset); /*lint !e777*/
435 
436  lhs = consdata1->lhs;
437  if( !SCIPisInfinity(scip, -lhs) )
438  lhs -= consdata2->lhs;
439  rhs = consdata1->rhs;
440  if( !SCIPisInfinity(scip, rhs) )
441  rhs -= consdata2->lhs;
442 
443  vars[0] = consdata1->z;
444  vars[1] = consdata2->z;
445 
446  coefs[0] = consdata1->zcoef;
447  coefs[1] = -consdata2->zcoef;
448 
449  if( SCIPisEQ(scip, lhs, rhs) )
450  {
451  SCIP_Bool redundant;
452  SCIP_Bool aggregated;
453 
454  /* try aggregation */
455  SCIP_CALL( SCIPaggregateVars(scip, consdata1->z, consdata2->z, consdata1->zcoef, -consdata2->zcoef, rhs, infeas, &redundant, &aggregated) );
456 
457  /* if infeasibility has been detected, stop here */
458  if( *infeas )
459  return SCIP_OKAY;
460 
461  if( redundant )
462  {
463  /* if redundant is TRUE, then either the aggregation has been done, or it was redundant */
464  if( aggregated )
465  ++*naggrvars;
466 
467  ++*ndelconss;
468 
469  SCIP_CALL( SCIPdelCons(scip, cons1) );
470  return SCIP_OKAY;
471  }
472  }
473 
474  /* if aggregation did not succeed, then either because some variable is multi-aggregated or due to numerics or because lhs != rhs
475  * we then add a linear constraint instead
476  */
477  vars[0] = consdata1->z;
478  vars[1] = consdata2->z;
479  coefs[0] = consdata1->zcoef;
480  coefs[1] = -consdata2->zcoef;
481 
482  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons1), 2, vars, coefs, lhs, rhs,
486  SCIPconsIsStickingAtNode(cons1)) );
487  SCIP_CALL( SCIPaddCons(scip, lincons) );
488  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
489 
490  SCIP_CALL( SCIPdelCons(scip, cons1) );
491  ++*nupgdconss;
492 
493  return SCIP_OKAY;
494 }
495 
496 /** solves a system of two absolute power equations
497  * Given: (x+xoffset1)|x+xoffset1|^{exponent-1} + zcoef1 * z == rhs1
498  * and (x+xoffset2)|x+xoffset2|^{exponent-1} + zcoef2 * z == rhs2
499  * with xoffset1 != xoffset2 and zcoef1 * rhs2 == zcoef2 * rhs1 and exponent == 2,
500  * finds values for x and z that satisfy these equations, or reports infeasibility if no solution exists.
501  *
502  * Multiplying the second equation by -zcoef1/zcoef2 and adding it to the first one gives
503  * (x+xoffset1)|x+xoffset1| - zcoef1/zcoef2 (x+offset2)|x+offset2| == 0
504  *
505  * If zcoef1 == zcoef2, then there exists, due to monotonicity of x|x|, no x such that
506  * (x+xoffset1)|x+xoffset1| == (x+xoffset2)|x+xoffset2|.
507  *
508  * In general, for zcoef1 / zcoef2 > 0.0, we get
509  * x = (xoffset2 - xoffset1) / (sqrt(zcoef2 / zcoef1) - 1.0) - xoffset1,
510  * and for zcoef1 / zcoef2 < 0.0, we get
511  * x = (xoffset2 - xoffset1) / (-sqrt(-zcoef2 / zcoef1) - 1.0) - xoffset1.
512  *
513  * This then yields z = (rhs1 - (x+xoffset1)|x+xoffset1|) / zcoef1.
514  */
515 static
517  SCIP* scip, /**< SCIP data structure */
518  SCIP_Bool* infeas, /**< buffer to indicate if the system of equations has no solution */
519  SCIP_Real* xval, /**< buffer to store value of x in the solution, if any */
520  SCIP_Real* zval, /**< buffer to store value of z in the solution, if any */
521  SCIP_Real exponent, /**< exponent in absolute power equations */
522  SCIP_Real xoffset1, /**< offset for x in first absolute power equation */
523  SCIP_Real zcoef1, /**< coefficient of z in first absolute power equation */
524  SCIP_Real rhs1, /**< right-hand-side in first absolute power equation */
525  SCIP_Real xoffset2, /**< offset for x in second absolute power equation */
526  SCIP_Real zcoef2, /**< coefficient of z in second absolute power equation */
527  SCIP_Real rhs2 /**< right-hand-side in second absolute power equation */
528  )
529 {
530  assert(scip != NULL);
531  assert(infeas != NULL);
532  assert(xval != NULL);
533  assert(zval != NULL);
534  assert(exponent == 2.0);
535  assert(!SCIPisEQ(scip, xoffset1, xoffset2));
536  assert(SCIPisEQ(scip, zcoef1 * rhs2, zcoef2 * rhs1));
537  assert(zcoef1 != 0.0);
538  assert(zcoef2 != 0.0);
539 
540  if( xoffset2 < xoffset1 )
541  {
542  presolveFindDuplicatesSolveEquations(scip, infeas, xval, zval, exponent, xoffset2, zcoef2, rhs2, xoffset1, zcoef1, rhs1);
543  return;
544  }
545 
546  if( SCIPisEQ(scip, zcoef1, zcoef2) )
547  {
548  *infeas = TRUE;
549  return;
550  }
551 
552  *infeas = FALSE;
553 
554  if( SCIPisEQ(scip, zcoef1, -zcoef2) )
555  {
556  *xval = - (xoffset1 + xoffset2) / 2.0;
557  }
558  else if( zcoef2 * zcoef1 > 0.0 )
559  {
560  *xval = (xoffset2 - xoffset1) / (sqrt(zcoef2 / zcoef1) - 1.0) - xoffset1;
561  }
562  else
563  {
564  assert(zcoef2 * zcoef1 < 0.0);
565  *xval = (xoffset2 - xoffset1) / (-sqrt(-zcoef2 / zcoef1) - 1.0) - xoffset1;
566  }
567 
568  *zval = rhs1 - (*xval + xoffset1) * REALABS(*xval + xoffset1);
569  *zval /= zcoef1;
570 
571  assert(SCIPisFeasEQ(scip, (*xval + xoffset1) * REALABS(*xval + xoffset1) + zcoef1 * *zval, rhs1));
572  assert(SCIPisFeasEQ(scip, (*xval + xoffset2) * REALABS(*xval + xoffset2) + zcoef2 * *zval, rhs2));
573 }
574 
575 /** finds and removes duplicates in a set of absolute power constraints */
576 static
578  SCIP* scip, /**< SCIP data structure */
579  SCIP_CONSHDLR* conshdlr, /**< constraint handler for absolute power constraints */
580  SCIP_CONS** conss, /**< constraints */
581  int nconss, /**< number of constraints */
582  int* nupgdconss, /**< pointer where to add number of upgraded constraints */
583  int* ndelconss, /**< pointer where to add number of deleted constraints */
584  int* naddconss, /**< pointer where to add number of added constraints */
585  int* nfixedvars, /**< pointer where to add number of fixed variables */
586  int* naggrvars, /**< pointer where to add number of aggregated variables */
587  SCIP_Bool* success, /**< pointer to store whether a duplicate was found (and removed) */
588  SCIP_Bool* infeas /**< pointer to store whether infeasibility was detected */
589  )
590 {
591  SCIP_MULTIHASH* multihash;
592  SCIP_MULTIHASHLIST* multihashlist;
593  SCIP_CONSHDLRDATA* conshdlrdata;
594  int c;
595 
596  assert(scip != NULL);
597  assert(conshdlr != NULL);
598  assert(conss != NULL || nconss == 0);
599  assert(nupgdconss != NULL);
600  assert(ndelconss != NULL);
601  assert(naddconss != NULL);
602  assert(nfixedvars != NULL);
603  assert(naggrvars != NULL);
604  assert(success != NULL);
605  assert(infeas != NULL);
606 
607  *success = FALSE;
608  *infeas = FALSE;
609 
610  if( nconss <= 1 )
611  return SCIP_OKAY;
612 
613  conshdlrdata = SCIPconshdlrGetData(conshdlr);
614  assert(conshdlrdata != NULL);
615 
616  /* check all constraints in the given set for duplicates, dominance, or possible simplifications w.r.t. the x variable */
617 
619  presolveFindDuplicatesGetKey, presolveFindDuplicatesKeyEQ, presolveFindDuplicatesKeyVal, (void*)scip) );
620 
621  for( c = 0; c < nconss && !*infeas; ++c )
622  {
623  SCIP_CONS* cons0;
624  SCIP_CONS* cons1;
625 
626  cons0 = conss[c]; /*lint !e613*/
627 
628  assert(!SCIPconsIsModifiable(cons0)); /* absolute power constraints aren't modifiable */
629  assert(!SCIPconsIsLocal(cons0)); /* shouldn't have local constraints in presolve */
630  assert(SCIPconsIsActive(cons0)); /* shouldn't get inactive constraints here */
631 
632  multihashlist = NULL;
633 
634  do
635  {
636  SCIP_CONSDATA* consdata0;
637  SCIP_CONSDATA* consdata1;
638 
639  /* get constraint from current hash table with same x variable as cons0 and same exponent */
640  cons1 = (SCIP_CONS*)(SCIPmultihashRetrieveNext(multihash, &multihashlist, (void*)cons0));
641  if( cons1 == NULL )
642  {
643  /* processed all constraints like cons0 from hash table, so insert cons0 and go to conss[c+1] */
644  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
645  break;
646  }
647 
648  assert(cons0 != cons1);
649 
650  consdata0 = SCIPconsGetData(cons0);
651  consdata1 = SCIPconsGetData(cons1);
652  assert(consdata0 != NULL);
653  assert(consdata1 != NULL);
654 
655  SCIPdebugPrintCons(scip, cons0, NULL);
656  SCIPdebugPrintCons(scip, cons1, NULL);
657 
658  assert(consdata0->x == consdata1->x);
659  assert(consdata0->exponent == consdata1->exponent); /*lint !e777*/
660 
661  if( SCIPisEQ(scip, consdata0->xoffset, consdata1->xoffset) )
662  {
663  /* we have two constraints with the same (x+offset)|x+offset|^n term */
664 
665  /* if both constraints have the same functions; strengthen sides of cons1 and throw cons0 away */
666  if( consdata0->z == consdata1->z && SCIPisEQ(scip, consdata0->zcoef, consdata1->zcoef) )
667  {
668  /* check if side strenghtening would result in inconsistency */
669  if( SCIPisGT(scip, consdata0->lhs, consdata1->rhs) || SCIPisGT(scip, consdata1->lhs, consdata0->rhs) )
670  {
671  SCIPdebugMsg(scip, "<%s> and <%s> are contradictory; declare infeasibility\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
672  *infeas = TRUE;
673  break;
674  }
675 
676  SCIPdebugMsg(scip, "<%s> and <%s> are equivalent; dropping the first\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
677 
678  /* if a side of cons1 gets finite via merging with cons0, then this changes locks and events */
679  if( (SCIPisInfinity(scip, -consdata1->lhs) && !SCIPisInfinity(scip, -consdata0->lhs)) ||
680  ( SCIPisInfinity(scip, consdata1->rhs) && !SCIPisInfinity(scip, consdata0->rhs)) )
681  {
682  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons1) );
683  SCIP_CALL( SCIPunlockVarCons(scip, consdata1->x, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
684  if( consdata1->zcoef > 0.0 )
685  SCIP_CALL( SCIPunlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
686  else
687  SCIP_CALL( SCIPunlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, consdata1->rhs), !SCIPisInfinity(scip, -consdata1->lhs)) );
688 
689  consdata1->lhs = MAX(consdata0->lhs, consdata1->lhs);
690  consdata1->rhs = MIN(consdata0->rhs, consdata1->rhs);
691 
692  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons1) );
693  SCIP_CALL( SCIPlockVarCons(scip, consdata1->x, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
694  if( consdata1->zcoef > 0.0 )
695  SCIP_CALL( SCIPlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, -consdata1->lhs), !SCIPisInfinity(scip, consdata1->rhs)) );
696  else
697  SCIP_CALL( SCIPlockVarCons(scip, consdata1->z, cons1, !SCIPisInfinity(scip, consdata1->rhs), !SCIPisInfinity(scip, -consdata1->lhs)) );
698  }
699  else
700  {
701  consdata1->lhs = MAX(consdata0->lhs, consdata1->lhs);
702  consdata1->rhs = MIN(consdata0->rhs, consdata1->rhs);
703  }
704 
705  SCIP_CALL( SCIPdelCons(scip, cons0) );
706  ++*ndelconss;
707  *success = TRUE;
708 
709  break;
710  }
711 
712  /* if cons1 defines a linear expression for sign(x+offset)|x+offset|^n, use it to replace cons0 by a linear constraint */
713  if( SCIPisEQ(scip, consdata1->lhs, consdata1->rhs) )
714  {
715  SCIPdebugMsg(scip, "substitute <%s> in <%s> to make linear constraint\n", SCIPconsGetName(cons1), SCIPconsGetName(cons0));
716  SCIP_CALL( presolveFindDuplicatesUpgradeCons(scip, cons0, cons1, infeas, nupgdconss, ndelconss, naggrvars) );
717 
718  *success = TRUE;
719  break;
720  }
721 
722  /* if cons0 defines a linear expression for sign(x+offset)|x+offset|^n, use it to replace cons1 by a linear constraint */
723  if( SCIPisEQ(scip, consdata0->lhs, consdata0->rhs) )
724  {
725  SCIPdebugMsg(scip, "substitute <%s> in <%s> to make linear constraint\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
726  SCIP_CALL( presolveFindDuplicatesUpgradeCons(scip, cons1, cons0, infeas, nupgdconss, ndelconss, naggrvars) );
727 
728  SCIP_CALL( SCIPmultihashRemove(multihash, cons1) );
729  *success = TRUE;
730 
731  if( *infeas )
732  break;
733  }
734  else
735  {
736  /* introduce a new equality constraint for sign(x+offset)|x+offset|^n and use it to replace cons0 and cons1 */
737  /* @todo maybe we could be more clever by looking which constraint sides are finite */
738  SCIP_VAR* auxvar;
739  SCIP_CONS* auxcons;
740  char name[SCIP_MAXSTRLEN];
741  SCIP_VAR* vars[2];
742  SCIP_Real coefs[2];
743 
744  SCIPdebugMsg(scip, "introduce new auxvar for signpower(%s+%g, %g) to make <%s> and <%s> linear constraint\n", SCIPvarGetName(consdata0->x), consdata0->exponent, consdata0->xoffset, SCIPconsGetName(cons0), SCIPconsGetName(cons1));
745 
746  /* create auxiliary variable to represent sign(x+offset)|x+offset|^n */
747  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "auxvar_abspower%s_%g_%g", SCIPvarGetName(consdata0->x), consdata0->exponent, consdata0->xoffset);
748  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
749  TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
750  SCIP_CALL( SCIPaddVar(scip, auxvar) );
751 
752  /* create auxiliary constraint auxvar = sign(x+offset)|x+offset|^n
753  * as we introduced a new variable, the constraint that "defines" the value for this variable need to be enforced, that is, is not redundent
754  */
755  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "auxcons_abspower%s_%g_%g", SCIPvarGetName(consdata0->x), consdata0->exponent, consdata0->xoffset);
756  SCIP_CALL( SCIPcreateConsAbspower(scip, &auxcons, name, consdata0->x, auxvar, consdata0->exponent, consdata0->xoffset, -1.0, 0.0, 0.0,
757  SCIPconsIsInitial(cons0) || SCIPconsIsInitial(cons1),
758  SCIPconsIsSeparated(cons0) || SCIPconsIsSeparated(cons1),
759  TRUE,
760  TRUE,
762  FALSE,
763  FALSE,
764  SCIPconsIsDynamic(cons0) || SCIPconsIsDynamic(cons1),
765  SCIPconsIsRemovable(cons0) || SCIPconsIsRemovable(cons1),
767  ) );
768  SCIP_CALL( SCIPaddCons(scip, auxcons) );
769  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
770  ++*naddconss;
771 
772 #ifdef SCIP_DEBUG_SOLUTION
773  if( SCIPdebugIsMainscip(scip) )
774  {
775  SCIP_Real xval;
776 
777  SCIP_CALL( SCIPdebugGetSolVal(scip, consdata0->x, &xval) );
778  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SIGN(xval + consdata0->xoffset) * pow(REALABS(xval + consdata0->xoffset), consdata0->exponent)) );
779  }
780 #endif
781 
782  /* create linear constraint equivalent for cons0 */
783  vars[0] = auxvar;
784  vars[1] = consdata0->z;
785  coefs[0] = 1.0;
786  coefs[1] = consdata0->zcoef;
787  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons0), 2, vars, coefs, consdata0->lhs, consdata0->rhs,
791  SCIPconsIsStickingAtNode(cons0)) );
792  SCIP_CALL( SCIPaddCons(scip, auxcons) );
793  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
794  ++*nupgdconss;
795 
796  /* create linear constraint equivalent for cons1 */
797  vars[1] = consdata1->z;
798  coefs[1] = consdata1->zcoef;
799  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons1), 2, vars, coefs, consdata1->lhs, consdata1->rhs,
803  SCIPconsIsStickingAtNode(cons1)) );
804  SCIP_CALL( SCIPaddCons(scip, auxcons) );
805  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
806  ++*nupgdconss;
807 
808  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
809 
810  SCIP_CALL( SCIPdelCons(scip, cons0) );
811  SCIP_CALL( SCIPdelCons(scip, cons1) );
812  SCIP_CALL( SCIPmultihashRemove(multihash, cons1) );
813  *success = TRUE;
814 
815  break;
816  }
817  }
818  else if( consdata0->z == consdata1->z &&
819  consdata0->exponent == 2.0 &&
820  !SCIPisZero(scip, consdata0->zcoef) &&
821  !SCIPisZero(scip, consdata1->zcoef) &&
822  SCIPisEQ(scip, consdata0->lhs, consdata0->rhs) &&
823  SCIPisEQ(scip, consdata1->lhs, consdata1->rhs) &&
824  SCIPisEQ(scip, consdata0->lhs * consdata1->zcoef, consdata1->lhs * consdata0->zcoef) )
825  {
826  /* If we have two equality constraints with the same variables and the same exponent and compatible constants,
827  * then this system of equations should have either no or a single solution.
828  * Thus, we can report cutoff or fix the variables to this solution, and forget about the constraints.
829  * @todo think about inequalities, differing exponents, and exponents != 2
830  */
831 
832  SCIP_Real xval;
833  SCIP_Real zval;
834 
835  assert(consdata0->x == consdata1->x);
836  assert(consdata0->exponent == consdata1->exponent); /*lint !e777*/
837  assert(!SCIPisEQ(scip, consdata0->xoffset, consdata1->xoffset));
838 
839  presolveFindDuplicatesSolveEquations(scip, infeas, &xval, &zval,
840  consdata0->exponent,
841  consdata0->xoffset, consdata0->zcoef, consdata0->lhs,
842  consdata1->xoffset, consdata1->zcoef, consdata1->lhs);
843 
844  if( *infeas )
845  {
846  SCIPdebugMsg(scip, "infeasibility detected while solving the equations, no solution exists\n");
847  SCIPdebugPrintCons(scip, cons0, NULL);
848  SCIPdebugPrintCons(scip, cons1, NULL);
849  break;
850  }
851 
852  SCIPdebugMsg(scip, "fixing variables <%s>[%g, %g] to %g and <%s>[%g, %g] to %g due to equations\n",
853  SCIPvarGetName(consdata0->x), SCIPvarGetLbLocal(consdata0->x), SCIPvarGetUbLocal(consdata0->x), xval,
854  SCIPvarGetName(consdata0->z), SCIPvarGetLbLocal(consdata0->z), SCIPvarGetUbLocal(consdata0->z), zval);
855  SCIPdebugPrintCons(scip, cons0, NULL);
856  SCIPdebugPrintCons(scip, cons1, NULL);
857 
859  {
860  SCIP_Bool fixed;
861 
862  SCIP_CALL( SCIPfixVar(scip, consdata0->x, xval, infeas, &fixed) );
863  ++*ndelconss;
864 
865  if( fixed )
866  ++*nfixedvars;
867 
868  if( *infeas )
869  {
870  SCIPdebugMsg(scip, "infeasibility detected after fixing <%s>\n", SCIPvarGetName(consdata0->x));
871  break;
872  }
873  }
874  else
875  {
876  SCIP_CONS* lincons;
877  SCIP_Real one;
878 
879  one = 1.0;
880  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons0), 1, &consdata0->x, &one, xval, xval,
882  SCIP_CALL( SCIPaddCons(scip, lincons) );
883  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
884  ++*nupgdconss;
885  }
886 
888  {
889  SCIP_Bool fixed;
890 
891  SCIP_CALL( SCIPfixVar(scip, consdata0->z, zval, infeas, &fixed) );
892  ++*ndelconss;
893 
894  if( fixed )
895  ++*nfixedvars;
896 
897  if( *infeas )
898  {
899  SCIPdebugMsg(scip, "infeasibility detected after fixing <%s>\n", SCIPvarGetName(consdata0->z));
900  break;
901  }
902  }
903  else
904  {
905  SCIP_CONS* lincons;
906  SCIP_Real one;
907 
908  one = 1.0;
909  SCIP_CALL( SCIPcreateConsLinear(scip, &lincons, SCIPconsGetName(cons1), 1, &consdata0->z, &one, zval, zval,
911  SCIP_CALL( SCIPaddCons(scip, lincons) );
912  SCIP_CALL( SCIPreleaseCons(scip, &lincons) );
913  ++*nupgdconss;
914  }
915 
916  SCIP_CALL( SCIPdelCons(scip, cons0) );
917  SCIP_CALL( SCIPdelCons(scip, cons1) );
918  SCIP_CALL( SCIPmultihashRemove(multihash, cons1) );
919  *success = TRUE;
920 
921  break;
922  }
923 
924  if( multihashlist == NULL )
925  {
926  /* processed all constraints like cons0 from hash table, but cons0 could not be removed, so insert cons0 into hashmap and go to conss[c+1] */
927  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
928  break;
929  }
930  }
931  while( TRUE ); /*lint !e506*/
932  }
933 
934  /* free hash table */
935  SCIPmultihashFree(&multihash);
936 
937  if( *infeas )
938  return SCIP_OKAY;
939 
940 
941  /* check all constraints in the given set for duplicates, dominance, or possible simplifications w.r.t. the z variable */
942 
944  presolveFindDuplicatesGetKey, presolveFindDuplicatesKeyEQ2, presolveFindDuplicatesKeyVal2, (void*) scip) );
945 
946  for( c = 0; c < nconss && !*infeas; ++c )
947  {
948  SCIP_CONS* cons0;
949  SCIP_CONS* cons1;
950  SCIP_CONSDATA* consdata0;
951 
952  cons0 = conss[c]; /*lint !e613*/
953 
954  assert(!SCIPconsIsModifiable(cons0)); /* absolute power constraints aren't modifiable */
955  assert(!SCIPconsIsLocal(cons0)); /* shouldn't have local constraints in presolve */
956 
957  /* do not consider constraints that we have deleted in the above loop */
958  if( SCIPconsIsDeleted(cons0) )
959  continue;
960  assert(SCIPconsIsActive(cons0)); /* shouldn't get inactive constraints here */
961 
962  consdata0 = SCIPconsGetData(cons0);
963  assert(consdata0 != NULL);
964 
965  /* consider only equality constraints so far
966  * @todo do also something with inequalities
967  */
968  if( !SCIPisEQ(scip, consdata0->lhs, consdata0->rhs) )
969  continue;
970 
971  multihashlist = NULL;
972 
973  do
974  {
975  SCIP_CONSDATA* consdata1;
976 
977  /* get constraint from current hash table with same z variable as cons0 and same exponent */
978  cons1 = (SCIP_CONS*)(SCIPmultihashRetrieveNext(multihash, &multihashlist, (void*)cons0));
979  if( cons1 == NULL )
980  {
981  /* processed all constraints like cons0 from hash table, so insert cons0 and go to conss[c+1] */
982  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
983  break;
984  }
985 
986  assert(cons0 != cons1);
987  assert(!SCIPconsIsDeleted(cons1));
988 
989  consdata1 = SCIPconsGetData(cons1);
990  assert(consdata1 != NULL);
991 
992  SCIPdebugPrintCons(scip, cons0, NULL);
993  SCIPdebugPrintCons(scip, cons1, NULL);
994 
995  assert(consdata0->z == consdata1->z);
996  assert(consdata0->exponent == consdata1->exponent); /*lint !e777*/
997  assert(SCIPisEQ(scip, consdata1->lhs, consdata1->rhs));
998  assert(!SCIPisZero(scip, consdata1->zcoef));
999 
1000  if( SCIPisEQ(scip, consdata0->lhs*consdata1->zcoef, consdata1->lhs*consdata0->zcoef) )
1001  {
1002  /* have two absolute power equations with same z and compatible constants
1003  * we can then reduce this to one absolute power and one linear equation
1004  * -> x0 + xoffset0 = signpower(zcoef0/zcoef1, 1/exponent) (x1 + xoffset1)
1005  * -> keep cons1
1006  * the latter can be realized as an aggregation (if x0 and x1 are not multiaggregated) or linear constraint
1007  */
1008  SCIP_Bool redundant;
1009  SCIP_Bool aggregated;
1010  SCIP_Real coef;
1011  SCIP_Real rhs;
1012 
1013  SCIPdebugMsg(scip, "<%s> and <%s> can be reformulated to one abspower and one aggregation\n", SCIPconsGetName(cons0), SCIPconsGetName(cons1));
1014  SCIPdebugPrintCons(scip, cons0, NULL);
1015  SCIPdebugPrintCons(scip, cons1, NULL);
1016 
1017  if( consdata0->exponent == 2.0 )
1018  coef = SIGN(consdata0->zcoef / consdata1->zcoef) * sqrt(REALABS(consdata0->zcoef / consdata1->zcoef));
1019  else
1020  coef = SIGN(consdata0->zcoef / consdata1->zcoef) * pow(REALABS(consdata0->zcoef / consdata1->zcoef), 1.0/consdata0->exponent);
1021  rhs = coef * consdata1->xoffset - consdata0->xoffset;
1022 
1023  /* try aggregation */
1024  SCIP_CALL( SCIPaggregateVars(scip, consdata0->x, consdata1->x, 1.0, -coef, rhs, infeas, &redundant, &aggregated) );
1025  if( *infeas )
1026  {
1027  /* if infeasibility has been detected, stop here */
1028  break;
1029  }
1030  else if( redundant )
1031  {
1032  /* if redundant is TRUE, then either the aggregation has been done, or it was redundant */
1033  if( aggregated )
1034  ++*naggrvars;
1035 
1036  ++*ndelconss;
1037  }
1038  else
1039  {
1040  /* if aggregation did not succeed, then either because some variable is multi-aggregated or due to numerics
1041  * we then add a linear constraint instead
1042  */
1043  SCIP_CONS* auxcons;
1044  SCIP_VAR* vars[2];
1045  SCIP_Real coefs[2];
1046 
1047  vars[0] = consdata0->x;
1048  vars[1] = consdata1->x;
1049  coefs[0] = 1.0;
1050  coefs[1] = -coef;
1051 
1052  /* create linear constraint equivalent for cons0 */
1053  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons0), 2, vars, coefs, rhs, rhs,
1057  SCIPconsIsStickingAtNode(cons0)) );
1058  SCIP_CALL( SCIPaddCons(scip, auxcons) );
1059  SCIPdebugPrintCons(scip, auxcons, NULL);
1060  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
1061 
1062  ++*nupgdconss;
1063  }
1064  SCIP_CALL( SCIPdelCons(scip, cons0) );
1065 
1066  *success = TRUE;
1067  break;
1068  }
1069 
1070  if( multihashlist == NULL )
1071  {
1072  /* processed all constraints like cons0 from hash table, but cons0 could not be removed, so insert cons0 into hashmap and go to conss[c+1] */
1073  SCIP_CALL( SCIPmultihashInsert(multihash, (void*) cons0) );
1074  break;
1075  }
1076  }
1077  while( TRUE ); /*lint !e506*/
1078  }
1079 
1080  /* free hash table */
1081  SCIPmultihashFree(&multihash);
1082 
1083  return SCIP_OKAY;
1084 }
1085 
1086 /** fix variables not appearing in any other constraint
1087  *
1088  * @todo generalize to inequalities
1089  * @todo generalize to support discrete variables
1090  * @todo generalize to arbitrary exponents also if z is in objective
1091  */
1092 static
1094  SCIP* scip, /**< SCIP data structure */
1095  SCIP_CONS* cons, /**< constraint */
1096  SCIP_Bool* cutoff, /**< buffer to indicate whether a cutoff was detected */
1097  int* ndelconss, /**< buffer to increase with the number of deleted constraint */
1098  int* nfixedvars /**< buffer to increase with the number of fixed variables */
1099  )
1100 {
1101  SCIP_CONSDATA* consdata;
1102  SCIP_Bool lhsexists;
1103  SCIP_Bool rhsexists;
1104 
1105  assert(scip != NULL);
1106  assert(cons != NULL);
1107  assert(cutoff != NULL);
1108  assert(nfixedvars != NULL);
1109  assert(ndelconss != NULL);
1110 
1111  /* only process checked constraints (for which the locks are increased);
1112  * otherwise we would have to check for variables with nlocks == 0, and these are already processed by the
1113  * dualfix presolver
1114  */
1115  if( !SCIPconsIsChecked(cons) )
1116  return SCIP_OKAY;
1117 
1118  consdata = SCIPconsGetData(cons);
1119  assert(consdata != NULL);
1120 
1121  /* skip dual presolve if multiaggregated variables are present for now (bounds are not updated, difficult to fix) */
1122  if( SCIPvarGetStatus(consdata->x) == SCIP_VARSTATUS_MULTAGGR )
1123  return SCIP_OKAY;
1124  if( SCIPvarGetStatus(consdata->z) == SCIP_VARSTATUS_MULTAGGR )
1125  return SCIP_OKAY;
1126 
1127  /* skip dual presolve if discrete variables are present for now (more difficult to compute fixing value) */
1128  if( SCIPvarGetType(consdata->x) <= SCIP_VARTYPE_INTEGER )
1129  return SCIP_OKAY;
1130  if( SCIPvarGetType(consdata->z) <= SCIP_VARTYPE_INTEGER )
1131  return SCIP_OKAY;
1132 
1133  /* we assume that domain propagation has been run and fixed variables were removed if possible */
1134  assert(!SCIPconsIsMarkedPropagate(cons));
1135  assert(consdata->zcoef != 0.0);
1136 
1137  lhsexists = !SCIPisInfinity(scip, -consdata->lhs);
1138  rhsexists = !SCIPisInfinity(scip, consdata->rhs);
1139 
1140  if( SCIPvarGetNLocksDown(consdata->x) == (lhsexists ? 1 : 0) &&
1141  SCIPvarGetNLocksUp(consdata->x) == (rhsexists ? 1 : 0) &&
1142  (consdata->zcoef > 0.0 ? SCIPvarGetNLocksDown(consdata->z) : SCIPvarGetNLocksUp(consdata->z)) == (lhsexists ? 1 : 0) &&
1143  (consdata->zcoef > 0.0 ? SCIPvarGetNLocksUp(consdata->z) : SCIPvarGetNLocksDown(consdata->z)) == (rhsexists ? 1 : 0) )
1144  {
1145  /* x and z are only locked by cons, so we can fix them to an optimal solution of
1146  * min xobj * x + zobj * z
1147  * s.t. lhs <= sign(x+offset)*abs(x+offset)^exponent + zcoef * z <= rhs
1148  * xlb <= x <= xub
1149  * zlb <= z <= zub
1150  */
1151  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) )
1152  {
1153  /* much simpler case where we can substitute z:
1154  * min xobj * x + zobj/zcoef * (rhs - sign(x+offset)*abs(x+offset)^exponent)
1155  * s.t. xlb <= x <= xub
1156  *
1157  * Since domain propagation had been applied, we can assume that for any valid value for x,
1158  * also the corresponding z value is valid.
1159  */
1160  SCIP_Real xfix;
1161  SCIP_Real xlb;
1162  SCIP_Real xub;
1163  SCIP_Real zfix;
1164  SCIP_Bool fixed;
1165 
1166  xlb = SCIPvarGetLbGlobal(consdata->x);
1167  xub = SCIPvarGetUbGlobal(consdata->x);
1168 
1169  if( SCIPisZero(scip, SCIPvarGetObj(consdata->z)) )
1170  {
1171  /* even simpler case where objective is linear in x */
1172  if( SCIPisZero(scip, SCIPvarGetObj(consdata->x)) )
1173  {
1174  /* simplest case where objective is zero:
1175  * if zero is within bounds, fix to zero, otherwise
1176  * fix x to middle of bounds for numerical stability. */
1177  if(SCIPisLT(scip, xlb, 0.0) && SCIPisGT(scip, xub, 0.0))
1178  xfix = 0.0;
1179  else
1180  xfix = 0.5 * (xlb + xub);
1181  }
1182  else
1183  {
1184  /* fix x to best bound */
1185  xfix = SCIPvarGetBestBoundGlobal(consdata->x);
1186  }
1187  }
1188  else if( consdata->exponent == 2.0 )
1189  {
1190  /* consider cases x <= -offset and x >= -offset separately */
1191  SCIP_Real a;
1192  SCIP_Real b;
1193  SCIP_Real c;
1194  SCIP_Real cand;
1195  SCIP_Real xfixobjval;
1196 
1197  xfix = SCIP_INVALID;
1198  xfixobjval = SCIP_INVALID;
1199 
1200  if( SCIPisLT(scip, xlb, -consdata->xoffset) )
1201  {
1202  /* For x <= -offset, the objective is equivalent to
1203  * zobj/zcoef * x^2 + (xobj + 2 offset zobj/zcoef) * x + offset^2 * zobj/zcoef + other constant
1204  * <-> a * x^2 + b * x + c
1205  *
1206  * critical values for x are xlb, MIN(xub,-offset), and -b/(2*a)
1207  */
1208  a = SCIPvarGetObj(consdata->z) / consdata->zcoef;
1209  b = SCIPvarGetObj(consdata->x) + 2 * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1210  c = consdata->xoffset * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1211 
1212  if( a < 0.0 && SCIPisInfinity(scip, -xlb) )
1213  {
1214  /* if a < 0.0, then a*x^2 is unbounded for x -> -infinity, thus fix x to -infinity */
1215  xfix = -SCIPinfinity(scip);
1216  xfixobjval = -SCIPinfinity(scip);
1217  }
1218  else
1219  {
1220  /* initialize with value for x=xlb */
1221  xfix = xlb;
1222  xfixobjval = a * xlb * xlb + b * xlb + c;
1223 
1224  /* compare with value for x=MIN(-offset,xub) */
1225  cand = MIN(-consdata->xoffset, xub);
1226  if( xfixobjval > a * cand * cand + b * cand + c )
1227  {
1228  xfix = cand;
1229  xfixobjval = a * cand * cand + b * cand + c;
1230  }
1231 
1232  /* compare with value for x=-b/(2*a), if within bounds */
1233  cand = -b/(2.0*a);
1234  if( cand > xlb && cand < -consdata->xoffset && cand < xub && xfixobjval > -b*b/(4.0*a) + c )
1235  {
1236  xfix = cand;
1237  xfixobjval = -b*b/(4.0*a) + c;
1238  }
1239  }
1240  }
1241 
1242  if( SCIPisGT(scip, xub, -consdata->xoffset) )
1243  {
1244  /* For x >= -offset, the objective is equivalent to
1245  * -zobj/zcoef * x^2 + (xobj - 2 offset zobj/zcoef) * x - offset^2 * zobj/zcoef + constants
1246  * <-> a * x^2 + b * x + c
1247  *
1248  * critical values for x are xub, MAX(xlb,-offset), and -b/(2*a)
1249  */
1250  a = -SCIPvarGetObj(consdata->z) / consdata->zcoef;
1251  b = SCIPvarGetObj(consdata->x) - 2 * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1252  c = -consdata->xoffset * consdata->xoffset * SCIPvarGetObj(consdata->z) / consdata->zcoef;
1253 
1254  if( a < 0.0 && SCIPisInfinity(scip, xub) )
1255  {
1256  /* if a < 0.0, then a*x^2 is unbounded for x -> infinity, thus fix x to infinity */
1257  xfix = SCIPinfinity(scip);
1258  /* not needed: xfixobjval = SCIPinfinity(scip); */
1259  }
1260  else
1261  {
1262  if( xfix == SCIP_INVALID ) /*lint !e777*/
1263  {
1264  /* initialize with value for x=xub */
1265  xfix = xub;
1266  xfixobjval = a * xub * xub + b * xub + c;
1267  }
1268  else
1269  {
1270  /* compare with value for x=xub */
1271  cand = xub;
1272  if( xfixobjval > a * cand * cand + b * cand + c )
1273  {
1274  xfix = cand;
1275  xfixobjval = a * cand * cand + b * cand + c;
1276  }
1277  }
1278 
1279  /* compare with value for x=MAX(xlb,-offset) */
1280  cand = MAX(xlb, -consdata->xoffset);
1281  if( xfixobjval > a * cand * cand + b * cand + c )
1282  {
1283  xfix = cand;
1284  xfixobjval = a * cand * cand + b * cand + c;
1285  }
1286 
1287  /* compare with value for x=-b/(2*a), if within bounds */
1288  cand = -b/(2.0*a);
1289  if( cand > xlb && cand > -consdata->xoffset && cand < xub && xfixobjval > -b*b/(4.0*a) + c )
1290  {
1291  xfix = cand;
1292  /* not needed: xfixobjval = -b*b/(4.0*a) + c; */
1293  }
1294  }
1295  }
1296  assert(xfix != SCIP_INVALID); /*lint !e777*/
1297  assert(SCIPisInfinity(scip, -xlb) || SCIPisLE(scip, xlb, xfix));
1298  assert(SCIPisInfinity(scip, xub) || SCIPisGE(scip, xub, xfix));
1299  }
1300  else
1301  {
1302  /* skip dual presolve for exponents != 2 and z in objective for now */
1303  return SCIP_OKAY;
1304  }
1305 
1306  /* compute fixing value for z */
1307  if( SCIPisInfinity(scip, xfix) )
1308  {
1309  if( consdata->zcoef > 0.0 )
1310  {
1311  assert(SCIPisInfinity(scip, -SCIPvarGetLbGlobal(consdata->z)));
1312  zfix = -SCIPinfinity(scip);
1313  }
1314  else
1315  {
1316  assert(SCIPisInfinity(scip, SCIPvarGetUbGlobal(consdata->z)));
1317  zfix = SCIPinfinity(scip);
1318  }
1319  }
1320  else if( SCIPisInfinity(scip, -xfix) )
1321  {
1322  if( consdata->zcoef > 0.0 )
1323  {
1324  assert(SCIPisInfinity(scip, SCIPvarGetUbGlobal(consdata->z)));
1325  zfix = SCIPinfinity(scip);
1326  }
1327  else
1328  {
1329  assert(SCIPisInfinity(scip, -SCIPvarGetLbGlobal(consdata->z)));
1330  zfix = -SCIPinfinity(scip);
1331  }
1332  }
1333  else
1334  {
1335  SCIP_Real zlb;
1336  SCIP_Real zub;
1337 
1338  zlb = SCIPvarGetLbGlobal(consdata->z);
1339  zub = SCIPvarGetUbGlobal(consdata->z);
1340  zfix = consdata->rhs - SIGN(xfix + consdata->xoffset) * consdata->power(ABS(xfix + consdata->xoffset), consdata->exponent);
1341  zfix /= consdata->zcoef;
1342 
1343  /* project zfix into box, it should be at least very close */
1344  assert(SCIPisFeasLE(scip, zlb, zfix));
1345  assert(SCIPisFeasGE(scip, zub, zfix));
1346  zfix = MAX(zlb, MIN(zub, zfix));
1347  }
1348 
1349  /* fix variables according to x=xfix */
1350  SCIPdebugMsg(scip, "dual presolve fixes x=<%s>[%g,%g] to %g and z=<%s>[%g,%g] to %g in cons <%s>\n",
1351  SCIPvarGetName(consdata->x), xlb, xub, xfix,
1352  SCIPvarGetName(consdata->z), SCIPvarGetLbGlobal(consdata->z), SCIPvarGetUbGlobal(consdata->z), zfix,
1353  SCIPconsGetName(cons));
1354  SCIPdebugPrintCons(scip, cons, NULL);
1355 
1356  /* fix x */
1357  SCIP_CALL( SCIPfixVar(scip, consdata->x, xfix, cutoff, &fixed) );
1358  if( *cutoff )
1359  return SCIP_OKAY;
1360  if( fixed )
1361  ++*nfixedvars;
1362 
1363  /* fix z */
1364  SCIP_CALL( SCIPfixVar(scip, consdata->z, zfix, cutoff, &fixed) );
1365  if( *cutoff )
1366  return SCIP_OKAY;
1367  if( fixed )
1368  ++*nfixedvars;
1369 
1370  /* delete constraint */
1371  SCIP_CALL( SCIPdelCons(scip, cons) );
1372  ++*ndelconss;
1373  }
1374  }
1375 
1376  return SCIP_OKAY;
1377 }
1378 
1379 /** given a variable and an interval, tightens the local bounds of this variable to the given interval */
1380 static
1382  SCIP* scip, /**< SCIP data structure */
1383  SCIP_VAR* var, /**< variable which bounds to tighten */
1384  SCIP_INTERVAL bounds, /**< new bounds */
1385  SCIP_Bool force, /**< force tightening even if below bound strengthening tolerance */
1386  SCIP_CONS* cons, /**< constraint that is propagated */
1387  SCIP_RESULT* result, /**< pointer to store the result of the propagation call */
1388  int* nchgbds, /**< buffer where to add the number of changed bounds */
1389  int* nfixedvars, /**< buffer where to add the number of fixed variables, can be equal to nchgbds */
1390  int* naddconss /**< buffer where to add the number of added constraints, can be NULL if force is FALSE */
1391  )
1392 {
1393  SCIP_Bool infeas;
1394  SCIP_Bool tightened;
1395 
1396  assert(scip != NULL);
1397  assert(var != NULL);
1398  assert(cons != NULL);
1399  assert(result != NULL);
1400  assert(nchgbds != NULL);
1401  assert(nfixedvars != NULL);
1402 
1403  *result = SCIP_DIDNOTFIND;
1404 
1405  if( SCIPisInfinity(scip, SCIPintervalGetInf(bounds)) || SCIPisInfinity(scip, -SCIPintervalGetSup(bounds)) )
1406  {
1407  /* domain outside [-infty, +infty] -> declare as infeasible */
1408  *result = SCIP_CUTOFF;
1409  return SCIP_OKAY;
1410  }
1411 
1412  /* if variable is not multiaggregated (or aggregated to a multiaggregated), then try SCIPfixVar or SCIPtightenVarLb/Ub
1413  * otherwise, if bound tightening is forced, add a linear constraint
1414  * otherwise, forget about the bound tightening
1415  */
1417  {
1418  /* check if variable can be fixed */
1419  if( SCIPisEQ(scip, bounds.inf, bounds.sup) )
1420  {
1421  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) )
1422  {
1423  /* if variable not fixed yet, then do so now */
1424  SCIP_Real fixval;
1425 
1426  if( bounds.inf != bounds.sup ) /*lint !e777*/
1427  fixval = (bounds.inf + bounds.sup) / 2.0;
1428  else
1429  fixval = bounds.inf;
1430  SCIP_CALL( SCIPfixVar(scip, var, fixval, &infeas, &tightened) );
1431 
1432  if( infeas )
1433  {
1434  SCIPdebugMsg(scip, "found <%s> infeasible due to fixing variable <%s>\n", SCIPconsGetName(cons), SCIPvarGetName(var));
1435  *result = SCIP_CUTOFF;
1436  return SCIP_OKAY;
1437  }
1438  if( tightened )
1439  {
1440  SCIPdebugMsg(scip, "fixed variable <%s> in constraint <%s> to %g\n", SCIPvarGetName(var), SCIPconsGetName(cons), SCIPvarGetLbLocal(var));
1441  ++*nfixedvars;
1442  *result = SCIP_REDUCEDDOM;
1443  }
1444  }
1445  else
1446  {
1447  /* only check if new fixing value is consistent with variable bounds, otherwise cutoff */
1448  if( SCIPisLT(scip, bounds.sup, SCIPvarGetUbLocal(var)) || SCIPisGT(scip, bounds.inf, SCIPvarGetLbLocal(var)) )
1449  {
1450  SCIPdebugMsg(scip, "found <%s> infeasible due to fixing fixed variable <%s>[%.20g,%.20g] to [%.20g,%.20g]\n",
1451  SCIPconsGetName(cons), SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), bounds.inf, bounds.sup);
1452  *result = SCIP_CUTOFF;
1453  return SCIP_OKAY;
1454  }
1455  }
1456 
1457  return SCIP_OKAY;
1458  }
1459 
1460  /* check if lower bound can be tightened */
1461  if( SCIPintervalGetInf(bounds) > SCIPvarGetLbLocal(var) )
1462  {
1463  assert(!SCIPisInfinity(scip, -SCIPintervalGetInf(bounds)));
1464  SCIP_CALL( SCIPtightenVarLb(scip, var, SCIPintervalGetInf(bounds), force, &infeas, &tightened) );
1465  if( infeas )
1466  {
1467  SCIPdebugMsg(scip, "found %s infeasible due to domain propagation for variable %s in constraint %s\n", SCIPconsGetName(cons), SCIPvarGetName(var), SCIPconsGetName(cons));
1468  *result = SCIP_CUTOFF;
1469  return SCIP_OKAY;
1470  }
1471  if( tightened )
1472  {
1473  SCIPdebugMsg(scip, "tightened lower bound of variable %s in constraint %s to %g\n", SCIPvarGetName(var), SCIPconsGetName(cons), SCIPvarGetLbLocal(var));
1474  ++*nchgbds;
1475  *result = SCIP_REDUCEDDOM;
1476  }
1477  }
1478 
1479  /* check if upper bound can be tightened */
1480  if( SCIPintervalGetSup(bounds) < SCIPvarGetUbLocal(var) )
1481  {
1482  assert(!SCIPisInfinity(scip, SCIPintervalGetSup(bounds)));
1483  SCIP_CALL( SCIPtightenVarUb(scip, var, SCIPintervalGetSup(bounds), force, &infeas, &tightened) );
1484  if( infeas )
1485  {
1486  SCIPdebugMsg(scip, "found %s infeasible due to domain propagation for linear variable %s in constraint %s\n", SCIPconsGetName(cons), SCIPvarGetName(var), SCIPconsGetName(cons));
1487  *result = SCIP_CUTOFF;
1488  return SCIP_OKAY;
1489  }
1490  if( tightened )
1491  {
1492  SCIPdebugMsg(scip, "tightened upper bound of variable %s in constraint %s to %g\n", SCIPvarGetName(var), SCIPconsGetName(cons), SCIPvarGetUbLocal(var));
1493  ++*nchgbds;
1494  *result = SCIP_REDUCEDDOM;
1495  }
1496  }
1497  }
1498  else if( force && (SCIPisLT(scip, SCIPvarGetLbLocal(var), bounds.inf) || SCIPisGT(scip, SCIPvarGetUbLocal(var), bounds.sup)) )
1499  {
1500  /* add a linear constraint bounds.inf <= x <= bounds.sup */
1501  SCIP_CONS* auxcons;
1502  SCIP_Bool local;
1503  SCIP_Real one;
1504 
1505  assert(naddconss != NULL);
1506 
1507  /* we add constraint as local constraint if we are during probing or if we are during solve and not at the root node */
1508  local = SCIPinProbing(scip) || (SCIPgetStage(scip) == SCIP_STAGE_SOLVING && (SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0));
1509 
1510  one = 1.0;
1511  SCIP_CALL( SCIPcreateConsLinear(scip, &auxcons, SCIPconsGetName(cons), 1, &var, &one, bounds.inf, bounds.sup,
1513  SCIPconsIsChecked(cons), SCIPconsIsPropagated(cons), local,
1514  FALSE, FALSE, TRUE, FALSE) );
1515 
1516  if( local )
1517  {
1518  SCIP_CALL( SCIPaddConsLocal(scip, auxcons, NULL) );
1519  }
1520  else
1521  {
1522  SCIP_CALL( SCIPaddCons(scip, auxcons) );
1523  }
1524  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
1525 
1526  ++*naddconss;
1527  *result = SCIP_CONSADDED;
1528  }
1529 
1530  return SCIP_OKAY;
1531 }
1532 
1533 /** computes bounds on z in a absolute power constraints for given bounds on x */
1534 static
1535 void computeBoundsZ(
1536  SCIP* scip, /**< SCIP data structure */
1537  SCIP_CONS* cons, /**< constraint */
1538  SCIP_INTERVAL xbnds, /**< bounds on x that are to be propagated */
1539  SCIP_INTERVAL* zbnds /**< buffer to store corresponding bounds on z */
1540  )
1541 {
1542  SCIP_CONSDATA* consdata;
1543  SCIP_Real bnd;
1544  SCIP_Real x;
1545 
1546  assert(scip != NULL);
1547  assert(cons != NULL);
1548  assert(zbnds != NULL);
1549  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), xbnds));
1550 
1551  consdata = SCIPconsGetData(cons);
1552  assert(consdata != NULL);
1553 
1554  SCIPintervalSetEntire(SCIPinfinity(scip), zbnds);
1555 
1556  /* apply zcoef*z <= rhs - signedpow(xbnds.inf + offset, n) */
1557  if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, -xbnds.inf) )
1558  {
1559  x = xbnds.inf - PROPVARTOL + consdata->xoffset;
1560  bnd = consdata->rhs + PROPSIDETOL - SIGN(x) * consdata->power(REALABS(x), consdata->exponent);
1561 
1562  if( consdata->zcoef > 0.0 )
1563  zbnds->sup = bnd / consdata->zcoef;
1564  else
1565  zbnds->inf = bnd / consdata->zcoef;
1566  }
1567 
1568  /* apply zcoef*z >= lhs - signedpow(xbnds.sup + offset, n) */
1569  if( !SCIPisInfinity(scip, -consdata->lhs) && !SCIPisInfinity(scip, xbnds.sup) )
1570  {
1571  x = xbnds.sup + PROPVARTOL + consdata->xoffset;
1572  bnd = consdata->lhs - PROPSIDETOL - SIGN(x) * consdata->power(REALABS(x), consdata->exponent);
1573 
1574  if( consdata->zcoef > 0.0 )
1575  zbnds->inf = bnd / consdata->zcoef;
1576  else
1577  zbnds->sup = bnd / consdata->zcoef;
1578  }
1579 
1580  SCIPdebugMsg(scip, "given x = [%.20g, %.20g], computed z = [%.20g, %.20g] via", xbnds.inf, xbnds.sup, zbnds->inf, zbnds->sup);
1581  SCIPdebugPrintCons(scip, cons, NULL);
1582 
1583  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), *zbnds));
1584 }
1585 
1586 /** computes bounds on x in a absolute power constraints for given bounds on z */
1587 static
1588 void computeBoundsX(
1589  SCIP* scip, /**< SCIP data structure */
1590  SCIP_CONS* cons, /**< constraint */
1591  SCIP_INTERVAL zbnds, /**< bounds on x that are to be propagated */
1592  SCIP_INTERVAL* xbnds /**< buffer to store corresponding bounds on z */
1593  )
1594 {
1595  SCIP_CONSDATA* consdata;
1596  SCIP_Real bnd;
1597  SCIP_Real z;
1598 
1599  assert(scip != NULL);
1600  assert(cons != NULL);
1601  assert(xbnds != NULL);
1602  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), zbnds));
1603 
1604  consdata = SCIPconsGetData(cons);
1605  assert(consdata != NULL);
1606 
1607  SCIPintervalSetEntire(SCIPinfinity(scip), xbnds);
1608 
1609  /* apply signedpow(x+offset, n) <= rhs - (zcoef * zbnds).inf */
1610  z = (consdata->zcoef > 0.0 ? zbnds.inf : zbnds.sup);
1611  if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, REALABS(z)) )
1612  {
1613  bnd = consdata->rhs + PROPSIDETOL - consdata->zcoef * z + REALABS(consdata->zcoef) * PROPVARTOL;
1614  if( consdata->exponent == 2.0 )
1615  bnd = SIGN(bnd) * sqrt(REALABS(bnd));
1616  else
1617  bnd = SIGN(bnd) * pow(REALABS(bnd), 1.0/consdata->exponent);
1618  xbnds->sup = bnd - consdata->xoffset;
1619  }
1620 
1621  /* apply signedpow(x+offset, n) >= lhs - (zcoef * zbnds).sup */
1622  z = (consdata->zcoef > 0.0 ? zbnds.sup : zbnds.inf);
1623  if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, REALABS(z)) )
1624  {
1625  bnd = consdata->lhs - PROPSIDETOL - consdata->zcoef * z - REALABS(consdata->zcoef) * PROPVARTOL;
1626  if( consdata->exponent == 2.0 )
1627  bnd = SIGN(bnd) * sqrt(REALABS(bnd));
1628  else
1629  bnd = SIGN(bnd) * pow(REALABS(bnd), 1.0/consdata->exponent);
1630  xbnds->inf = bnd - consdata->xoffset;
1631  }
1632 
1633  SCIPdebugMsg(scip, "given z = [%.20g, %.20g], computed x = [%.20g, %.20g] via", zbnds.inf, zbnds.sup, xbnds->inf, xbnds->sup);
1634  SCIPdebugPrintCons(scip, cons, NULL);
1635 
1636  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), *xbnds));
1637 }
1638 
1639 /** checks if x or z is fixed and replaces them or deletes constraint */
1640 static
1642  SCIP* scip, /**< SCIP data structure */
1643  SCIP_CONSHDLR* conshdlr, /**< constraint handler for absolute power constraints */
1644  SCIP_CONS* cons, /**< constraint */
1645  int* ndelconss, /**< counter for number of deleted constraints */
1646  int* nupgdconss, /**< counter for number of upgraded constraints */
1647  int* nchgbds, /**< counter for number of variable bound changes */
1648  int* nfixedvars, /**< counter for number of variable fixations */
1649  SCIP_RESULT* result /**< to store result if we detect infeasibility or remove constraint */
1650  )
1651 {
1652  SCIP_CONSHDLRDATA* conshdlrdata;
1653  SCIP_CONSDATA* consdata;
1654  SCIP_Real scalar;
1655  SCIP_Real constant;
1656  SCIP_Real factor;
1657  SCIP_VAR* var;
1658 
1659  assert(scip != NULL);
1660  assert(cons != NULL);
1661  assert(ndelconss != NULL);
1662  assert(nupgdconss != NULL);
1663  assert(nchgbds != NULL);
1664  assert(nfixedvars != NULL);
1665 
1666  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1667  assert(conshdlrdata != NULL);
1668 
1669  consdata = SCIPconsGetData(cons);
1670  assert(consdata != NULL);
1671 
1672  *result = SCIP_DIDNOTFIND;
1673 
1674  if( !SCIPvarIsActive(consdata->x) && SCIPvarGetStatus(consdata->x) != SCIP_VARSTATUS_MULTAGGR )
1675  {
1676  /* replace x variable */
1677 
1678  /* get relation x = scalar * var + constant */
1679  var = consdata->x;
1680  scalar = 1.0;
1681  constant = 0.0;
1682  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &scalar, &constant) );
1683 
1684  if( scalar == 0.0 )
1685  {
1686  SCIP_INTERVAL xbnds;
1687  SCIP_INTERVAL zbnds;
1688  int naddconss;
1689 
1690  naddconss = 0;
1691 
1692  /* x has been fixed to constant */
1693  assert(SCIPisFeasEQ(scip, SCIPvarGetLbGlobal(consdata->x), constant));
1694  assert(SCIPisFeasEQ(scip, SCIPvarGetUbGlobal(consdata->x), constant));
1695 
1696  /* compute corresponding bounds on z */
1697  SCIPintervalSet(&xbnds, constant);
1698  computeBoundsZ(scip, cons, xbnds, &zbnds);
1699 
1700  SCIPdebugMsg(scip, "in cons <%s>: x = <%s> fixed to %g -> tighten <%s> to [%g, %g]\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->x), constant, SCIPvarGetName(consdata->z), zbnds.inf, zbnds.sup);
1701 
1702  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) )
1703  {
1704  /* if sides are equal, then we should either fix z, or declare infeasibility */
1705  if( SCIPisFeasLT(scip, SCIPvarGetUbGlobal(consdata->z), zbnds.inf) || SCIPisFeasGT(scip, SCIPvarGetLbGlobal(consdata->z), zbnds.sup) )
1706  {
1707  SCIPdebugMsg(scip, "bounds inconsistent -> cutoff\n");
1708  *result = SCIP_CUTOFF;
1709  return SCIP_OKAY;
1710  }
1711  else
1712  {
1713  /* compute fixing value for z as value corresponding to fixing of x, projected onto bounds of z */
1714  SCIP_Real zfix;
1715 
1716  zfix = consdata->rhs - SIGN(constant + consdata->xoffset) * consdata->power(REALABS(constant + consdata->xoffset), consdata->exponent);
1717  zfix /= consdata->zcoef;
1718  assert(SCIPisLE(scip, zbnds.inf, zfix));
1719  assert(SCIPisGE(scip, zbnds.sup, zfix));
1720  zfix = MIN(SCIPvarGetUbGlobal(consdata->z), MAX(SCIPvarGetLbGlobal(consdata->z), zfix)); /*lint !e666*/
1721 
1722  zbnds.inf = zfix;
1723  zbnds.sup = zfix;
1724  SCIP_CALL( tightenBounds(scip, consdata->z, zbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1725  }
1726  }
1727  else
1728  {
1729  /* tighten bounds on z accordingly */
1730  SCIP_CALL( tightenBounds(scip, consdata->z, zbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1731  }
1732 
1733  /* delete constraint */
1734  SCIP_CALL( SCIPdelCons(scip, cons) );
1735 
1736  /* if tightenBounds added a constraint (because z was multiaggregated), then count this as constraint upgrade, otherwise as constraint deletion */
1737  if( naddconss > 0 )
1738  ++*nupgdconss;
1739  else
1740  ++*ndelconss;
1741 
1742  return SCIP_OKAY;
1743  }
1744 
1745  SCIPdebugMsg(scip, "in cons <%s>: x = <%s> replaced by %g*<%s> + %g\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->x), scalar, SCIPvarGetName(var), constant);
1746 
1747  /* constraint will be divided by scalar*pow(|scalar|,exponent-1), if scalar is not 1.0 */
1748  if( scalar == 1.0 )
1749  factor = 1.0;
1750  else if( scalar > 0.0 )
1751  factor = consdata->power( scalar, consdata->exponent);
1752  else
1753  factor = -consdata->power(-scalar, consdata->exponent);
1754 
1755  /* aggregate only if this would not lead to a vanishing or infinite coefficient for z */
1756  if( !SCIPisZero(scip, consdata->zcoef / factor) && !SCIPisInfinity(scip, REALABS(consdata->zcoef / factor)) )
1757  {
1758  /* we drop here the events for both variables, because if x is replaced by a multiaggregated variable here, then we do not need to catch bound tightenings on z anymore */
1759  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1760  SCIP_CALL( SCIPunlockVarCons(scip, consdata->x, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1761 
1762  consdata->x = var;
1763  if( SCIPvarIsActive(consdata->x) )
1764  {
1765  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, consdata->x) );
1766  }
1767 
1768  /* add constant to offset */
1769  consdata->xoffset += constant;
1770 
1771  /* divide constraint by factor */
1772  if( scalar == 1.0 ) ;
1773  else if( scalar > 0.0 )
1774  {
1775  if( !SCIPisInfinity(scip, -consdata->lhs) )
1776  consdata->lhs /= factor;
1777  if( !SCIPisInfinity(scip, consdata->rhs) )
1778  consdata->rhs /= factor;
1779  consdata->zcoef /= factor;
1780  consdata->xoffset /= scalar;
1781  }
1782  else
1783  {
1784  SCIP_Real oldlhs;
1785 
1786  assert(scalar < 0.0);
1787  assert(factor < 0.0);
1788 
1789  oldlhs = consdata->lhs;
1790 
1791  if( !SCIPisInfinity(scip, consdata->rhs) )
1792  consdata->lhs = consdata->rhs / factor;
1793  else
1794  consdata->lhs = -SCIPinfinity(scip);
1795  if( !SCIPisInfinity(scip, -oldlhs) )
1796  consdata->rhs = oldlhs / factor;
1797  else
1798  consdata->rhs = SCIPinfinity(scip);
1799  consdata->zcoef /= factor;
1800  consdata->xoffset /= scalar;
1801  /* since we flip both constraint sides and the sign of zcoef, the events catched for z remain the same, so update necessary there */
1802  }
1803 
1804  SCIP_CALL( SCIPlockVarCons(scip, consdata->x, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1805  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1806 
1807  SCIPdebugPrintCons(scip, cons, NULL);
1808 
1809  /* rerun constraint comparison */
1810  conshdlrdata->comparedpairwise = FALSE;
1811  }
1812  else
1813  {
1814  SCIPwarningMessage(scip, "Skip resolving aggregation of variable <%s> in abspower constraint <%s> to avoid zcoef = %g\n",
1815  SCIPvarGetName(consdata->x), SCIPconsGetName(cons), consdata->zcoef / factor);
1816  }
1817  }
1818 
1819  if( !SCIPvarIsActive(consdata->z) && SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR )
1820  {
1821  /* replace z variable */
1822 
1823  /* get relation z = scalar * var + constant */
1824  var = consdata->z;
1825  scalar = 1.0;
1826  constant = 0.0;
1827  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &scalar, &constant) );
1828 
1829  if( scalar == 0.0 )
1830  {
1831  SCIP_INTERVAL xbnds;
1832  SCIP_INTERVAL zbnds;
1833  int naddconss;
1834 
1835  naddconss = 0;
1836 
1837  /* z has been fixed to constant */
1838  assert(SCIPisFeasEQ(scip, SCIPvarGetLbGlobal(consdata->z), constant));
1839  assert(SCIPisFeasEQ(scip, SCIPvarGetUbGlobal(consdata->z), constant));
1840 
1841  /* compute corresponding bounds on x */
1842  SCIPintervalSet(&zbnds, constant);
1843  computeBoundsX(scip, cons, zbnds, &xbnds);
1844 
1845  SCIPdebugMsg(scip, "in cons <%s>: z = <%s> fixed to %g -> tighten <%s> to [%g, %g]\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->z), constant, SCIPvarGetName(consdata->x), xbnds.inf, xbnds.sup);
1846 
1847  if( SCIPisEQ(scip, consdata->lhs, consdata->rhs) )
1848  {
1849  /* if sides are equal, then we should either fix x, or declare infeasibility */
1850  if( SCIPisFeasLT(scip, SCIPvarGetUbGlobal(consdata->x), xbnds.inf) || SCIPisFeasGT(scip, SCIPvarGetLbGlobal(consdata->x), xbnds.sup) )
1851  {
1852  SCIPdebugMsg(scip, "bounds inconsistent -> cutoff\n");
1853  *result = SCIP_CUTOFF;
1854  return SCIP_OKAY;
1855  }
1856  else
1857  {
1858  /* compute fixing value for x as value corresponding to fixing of z, projected onto bounds of x */
1859  SCIP_Real xfix;
1860 
1861  xfix = consdata->rhs - consdata->zcoef * constant;
1862  if( consdata->exponent == 2.0 )
1863  xfix = SIGN(xfix) * sqrt(REALABS(xfix)) - consdata->xoffset;
1864  else
1865  xfix = SIGN(xfix) * pow(REALABS(xfix), 1.0/consdata->exponent) - consdata->xoffset;
1866  assert(SCIPisLE(scip, xbnds.inf, xfix));
1867  assert(SCIPisGE(scip, xbnds.sup, xfix));
1868  xfix = MIN(SCIPvarGetUbGlobal(consdata->x), MAX(SCIPvarGetLbGlobal(consdata->x), xfix)); /*lint !e666*/
1869 
1870  xbnds.inf = xfix;
1871  xbnds.sup = xfix;
1872  SCIP_CALL( tightenBounds(scip, consdata->x, xbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1873  }
1874  }
1875  else
1876  {
1877  /* tighten bounds on x accordingly */
1878  SCIP_CALL( tightenBounds(scip, consdata->x, xbnds, TRUE, cons, result, nchgbds, nfixedvars, &naddconss) );
1879  }
1880 
1881  /* delete constraint */
1882  SCIP_CALL( SCIPdelCons(scip, cons) );
1883 
1884  /* if tightenBounds added a constraint (because x was multiaggregated), then count this as constraint upgrade, otherwise as constraint deletion */
1885  if( naddconss > 0 )
1886  ++*nupgdconss;
1887  else
1888  ++*ndelconss;
1889 
1890  return SCIP_OKAY;
1891  }
1892 
1893  SCIPdebugMsg(scip, "in cons <%s>: z = <%s> replaced by %g*<%s> + %g\n", SCIPconsGetName(cons), SCIPvarGetName(consdata->z), scalar, SCIPvarGetName(var), constant);
1894 
1895  /* we drop here the events for both variables, because if z is replaced by a multiaggregated variable here, then we do not need to catch bound tightenings on x anymore */
1896  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1897  if( consdata->zcoef > 0.0 )
1898  SCIP_CALL( SCIPunlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1899  else
1900  SCIP_CALL( SCIPunlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
1901 
1902  consdata->z = var;
1903  if( SCIPvarIsActive(consdata->z) )
1904  {
1905  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, consdata->z) );
1906  }
1907 
1908  /* substract constant from constraint sides */
1909  if( !SCIPisInfinity(scip, -consdata->lhs) )
1910  consdata->lhs -= consdata->zcoef * constant;
1911  if( !SCIPisInfinity(scip, consdata->rhs) )
1912  consdata->rhs -= consdata->zcoef * constant;
1913 
1914  /* multiply zcoef by scalar */
1915  consdata->zcoef *= scalar;
1916 
1917  if( consdata->zcoef > 0.0 )
1918  SCIP_CALL( SCIPlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
1919  else
1920  SCIP_CALL( SCIPlockVarCons(scip, consdata->z, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
1921  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, cons) );
1922 
1923  /* rerun constraint comparison */
1924  conshdlrdata->comparedpairwise = FALSE;
1925  }
1926 
1927  assert(SCIPvarIsActive(consdata->z) || SCIPvarGetStatus(consdata->z) == SCIP_VARSTATUS_MULTAGGR);
1928 
1929  return SCIP_OKAY;
1930 }
1931 
1932 /** gets maximal absolute value in gradient of quadratic function
1933  * thus, gives \f$max(n |x+offset|^{n-1}, |zcoef|)\f$.
1934  */
1935 static
1937  SCIP* scip, /**< SCIP data structure */
1938  SCIP_CONS* cons, /**< constraint */
1939  SCIP_SOL* sol /**< solution or NULL if LP solution should be used */
1940  )
1941 {
1942  SCIP_CONSDATA* consdata;
1943  SCIP_Real xval;
1944  SCIP_Real val;
1945 
1946  assert(scip != NULL);
1947  assert(cons != NULL);
1948 
1949  consdata = SCIPconsGetData(cons);
1950  assert(consdata != NULL);
1951 
1952  xval = SCIPgetSolVal(scip, sol, consdata->x);
1953  assert(!SCIPisInfinity(scip, REALABS(xval)));
1954 
1955  if( consdata->exponent == 2.0 )
1956  val = consdata->exponent * REALABS(xval + consdata->xoffset);
1957  else
1958  val = consdata->exponent * pow(REALABS(xval + consdata->xoffset), consdata->exponent - 1.0);
1959 
1960  return MAX(val, REALABS(consdata->zcoef)); /*lint !e666*/
1961 }
1962 
1963 /** computes violation of a constraint */
1964 static
1966  SCIP* scip, /**< SCIP data structure */
1967  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1968  SCIP_CONS* cons, /**< constraint */
1969  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
1970  SCIP_Real* viol, /**< pointer to store absolute (unscaled) constraint violation */
1971  SCIP_Bool* solviolbounds /**< buffer to store whether the solution violates bounds on x by more than feastol */
1972  )
1973 {
1974  SCIP_CONSHDLRDATA* conshdlrdata;
1975  SCIP_CONSDATA* consdata;
1976  SCIP_Real val;
1977  SCIP_Real xval;
1978  SCIP_Real zval;
1979 
1980  assert(scip != NULL);
1981  assert(conshdlr != NULL);
1982  assert(cons != NULL);
1983  assert(viol != NULL);
1984  assert(solviolbounds != NULL);
1985 
1986  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1987  assert(conshdlrdata != NULL);
1988 
1989  consdata = SCIPconsGetData(cons);
1990  assert(consdata != NULL);
1991 
1992  *solviolbounds = FALSE;
1993  xval = SCIPgetSolVal(scip, sol, consdata->x);
1994  zval = SCIPgetSolVal(scip, sol, consdata->z);
1995 
1996  if( SCIPisInfinity(scip, REALABS(xval)) )
1997  {
1998  consdata->lhsviol = (SCIPisInfinity(scip, -consdata->lhs) ? 0.0 : SCIPinfinity(scip));
1999  consdata->rhsviol = (SCIPisInfinity(scip, consdata->rhs) ? 0.0 : SCIPinfinity(scip));
2000 
2001  return SCIP_OKAY;
2002  }
2003 
2004  if( sol == NULL )
2005  {
2006  SCIP_Real lb;
2007  SCIP_Real ub;
2008 
2009  lb = SCIPvarGetLbLocal(consdata->x);
2010  ub = SCIPvarGetUbLocal(consdata->x);
2011 
2012  /* 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 be outside bounds */
2013  if( !SCIPisFeasGE(scip, xval, lb) || !SCIPisFeasLE(scip, xval, ub) )
2014  *solviolbounds = TRUE;
2015  else
2016  xval = MAX(lb, MIN(ub, xval)); /* project x onto local box if just slightly outside (or even not outside) */
2017  }
2018 
2019  xval += consdata->xoffset;
2020 
2021  val = SIGN(xval) * consdata->power(REALABS(xval), consdata->exponent);
2022  val += consdata->zcoef * zval;
2023 
2024  if( val < consdata->lhs && !SCIPisInfinity(scip, -consdata->lhs) )
2025  consdata->lhsviol = *viol = consdata->lhs - val;
2026  else
2027  consdata->lhsviol = 0.0;
2028 
2029  if( val > consdata->rhs && !SCIPisInfinity(scip, consdata->rhs) )
2030  consdata->rhsviol = *viol = val - consdata->rhs;
2031  else
2032  consdata->rhsviol = 0.0;
2033 
2034  switch( conshdlrdata->scaling )
2035  {
2036  case 'o' :
2037  /* no scaling */
2038  break;
2039 
2040  case 'g' :
2041  /* scale by sup-norm of gradient in current point */
2042  if( consdata->lhsviol > 0.0 || consdata->rhsviol > 0.0 )
2043  {
2044  SCIP_Real norm;
2045  norm = getGradientMaxElement(scip, cons, sol);
2046  if( norm > 1.0 )
2047  {
2048  consdata->lhsviol /= norm;
2049  consdata->rhsviol /= norm;
2050  }
2051  }
2052  break;
2053 
2054  case 's' :
2055  {
2056  SCIP_Real absval;
2057 
2058  /* scale by left/right hand side of constraint */
2059  if( consdata->lhsviol > 0.0 )
2060  {
2061  absval = REALABS(consdata->lhs);
2062  consdata->lhsviol /= MAX(1.0, absval);
2063  }
2064 
2065  if( consdata->rhsviol > 0.0 )
2066  {
2067  absval = REALABS(consdata->rhs);
2068  consdata->rhsviol /= MAX(1.0, absval);
2069  }
2070 
2071  break;
2072  }
2073 
2074  default :
2075  SCIPerrorMessage("Unknown scaling method '%c'.", conshdlrdata->scaling);
2076  SCIPABORT();
2077  return SCIP_INVALIDDATA; /*lint !e527*/
2078  }
2079 
2080  return SCIP_OKAY;
2081 }
2082 
2083 /** computes violation of a set of constraints */
2084 static
2086  SCIP* scip, /**< SCIP data structure */
2087  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2088  SCIP_CONS** conss, /**< constraints */
2089  int nconss, /**< number of constraints */
2090  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
2091  SCIP_Bool* solviolbounds, /**< buffer to store whether the solution violates bounds on x by more than feastol */
2092  SCIP_CONS** maxviolcon /**< buffer to store constraint with largest violation, or NULL if solution is feasible */
2093  )
2094 {
2095  SCIP_CONSDATA* consdata;
2096  SCIP_Real viol;
2097  SCIP_Real maxviol;
2098  SCIP_Bool solviolbounds1;
2099  int c;
2100 
2101  assert(scip != NULL);
2102  assert(conss != NULL || nconss == 0);
2103  assert(solviolbounds != NULL);
2104  assert(maxviolcon != NULL);
2105 
2106  *solviolbounds = FALSE;
2107  *maxviolcon = NULL;
2108 
2109  maxviol = 0.0;
2110 
2111  for( c = 0; c < nconss; ++c )
2112  {
2113  assert(conss != NULL);
2114  assert(conss[c] != NULL);
2115 
2116  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol, &viol, &solviolbounds1) );
2117  *solviolbounds |= solviolbounds1;
2118 
2119  consdata = SCIPconsGetData(conss[c]);
2120  assert(consdata != NULL);
2121 
2122  viol = MAX(consdata->lhsviol, consdata->rhsviol);
2123  if( viol > maxviol && SCIPisGT(scip, viol, SCIPfeastol(scip)) )
2124  {
2125  maxviol = viol;
2126  *maxviolcon = conss[c];
2127  }
2128  }
2129 
2130  return SCIP_OKAY;
2131 }
2132 
2133 /** proposes branching point for constraint */
2134 static
2136  SCIP* scip, /**< SCIP data structure */
2137  SCIP_CONS* cons, /**< constraint which variable to get branching point for */
2138  SCIP_SOL* sol, /**< solution to branch on (NULL for LP or pseudosol) */
2139  int preferzero, /**< how much we prefer branching on -xoffset (0, 1, or 2) if sign is not fixed */
2140  SCIP_Bool branchminconverror /**< whether to minimize convexification error if sign is fixed */
2141  )
2142 {
2143  SCIP_CONSDATA* consdata;
2144  SCIP_VAR* x;
2145  SCIP_Real xref;
2146  SCIP_Real zref;
2147  SCIP_Real xlb;
2148  SCIP_Real xub;
2149 
2150  assert(scip != NULL);
2151  assert(cons != NULL);
2152 
2153  consdata = SCIPconsGetData(cons);
2154  assert(consdata != NULL);
2155 
2156  x = consdata->x;
2157  xlb = SCIPvarGetLbLocal(x);
2158  xub = SCIPvarGetUbLocal(x);
2159 
2160  /* check if sign of x is not fixed yet */
2161  if( SCIPisLT(scip, xlb, -consdata->xoffset) && SCIPisGT(scip, xub, -consdata->xoffset) )
2162  {
2163  /* if preferzero is 0, just return SCIP_INVALID
2164  * if preferzero is 1, then propose -xoffset if branching on -xoffset would cut off solution in both child nodes, otherwise return SCIP_INVALID
2165  * if preferzero is >1, then always propose -xoffset
2166  */
2167  assert(preferzero >= 0);
2168 
2169  if( preferzero == 0 )
2170  return SCIP_INVALID;
2171 
2172  if( preferzero > 1 || SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) )
2173  return -consdata->xoffset;
2174 
2175  xlb += consdata->xoffset;
2176  xub += consdata->xoffset;
2177 
2178  xref = SCIPgetSolVal(scip, sol, x) + consdata->xoffset;
2179  zref = SCIPgetSolVal(scip, sol, consdata->z);
2180  if( SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
2181  {
2182  /* signpow(x,n,offset) + c*z <= 0 is violated
2183  * if we are close to or right of -offset, then branching on -offset gives a convex function on the right branch, this is good
2184  * otherwise if branching on -offset yields a violated secant cut in left branch, then current solution would be cutoff there, this is also still good
2185  */
2186  if( !SCIPisFeasNegative(scip, xref) || SCIPisFeasPositive(scip, -consdata->power(-xlb, consdata->exponent)*xref/xlb + consdata->zcoef * zref) )
2187  return -consdata->xoffset;
2188  return SCIP_INVALID;
2189  }
2190 
2191  assert(SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) );
2192  /* signpow(x,n) + c*z >= 0 is violated
2193  * if we are close to or left of zero, then branching on 0.0 gives a concave function on the left branch, this is good
2194  * otherwise if branching on 0.0 yields a violated secant cut in right branch, then current solution would be cutoff there, this is also still good
2195  */
2196  if( !SCIPisFeasPositive(scip, xref) || SCIPisFeasNegative(scip, -consdata->power(xub, consdata->exponent)*xref/xub + consdata->zcoef * zref) )
2197  return -consdata->xoffset;
2198  return SCIP_INVALID;
2199  }
2200 
2201  if( branchminconverror )
2202  {
2203  /* given x^n with xlb <= x <= xub, then the sum of the integrals between the function and its secant on the left and right branches are minimized
2204  * for branching on ( (ub^n - lb^n) / (n*(ub - lb)) ) ^ (1/(n-1))
2205  */
2206  if( SCIPisGE(scip, xlb, -consdata->xoffset) )
2207  {
2208  SCIP_Real ref;
2209  xlb = MAX(0.0, xlb + consdata->xoffset);
2210  xub = MAX(0.0, xub + consdata->xoffset);
2211 
2212  ref = (consdata->power(xub, consdata->exponent) - consdata->power(xlb, consdata->exponent)) / (consdata->exponent * (xub - xlb));
2213  ref = pow(ref, 1.0/(consdata->exponent-1.0));
2214  ref -= consdata->xoffset;
2215  assert(SCIPisGE(scip, ref, SCIPvarGetLbLocal(x)));
2216  assert(SCIPisLE(scip, ref, SCIPvarGetUbLocal(x)));
2217 
2218  return ref;
2219  }
2220  else
2221  {
2222  SCIP_Real ref;
2223 
2224  assert(SCIPisLE(scip, xub, -consdata->xoffset));
2225 
2226  xlb = MIN(0.0, xlb + consdata->xoffset);
2227  xub = MIN(0.0, xub + consdata->xoffset);
2228 
2229  ref = (consdata->power(-xlb, consdata->exponent) - consdata->power(-xub, consdata->exponent)) / (consdata->exponent * (-xlb + xub));
2230  ref = -pow(ref, 1.0/(consdata->exponent-1.0));
2231  ref -= consdata->xoffset;
2232  assert(SCIPisGE(scip, ref, SCIPvarGetLbLocal(x)));
2233  assert(SCIPisLE(scip, ref, SCIPvarGetUbLocal(x)));
2234 
2235  return ref;
2236  }
2237  }
2238 
2239  return SCIP_INVALID;
2240 }
2241 
2242 /** registers branching variable candidates
2243  * registers x for all violated absolute power constraints where x is not in convex region
2244  */
2245 static
2247  SCIP* scip, /**< SCIP data structure */
2248  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2249  SCIP_CONS** conss, /**< constraints to check */
2250  int nconss, /**< number of constraints to check */
2251  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
2252  int* nnotify /**< counter for number of notifications performed */
2253  )
2254 {
2255  SCIP_CONSHDLRDATA* conshdlrdata;
2256  SCIP_CONSDATA* consdata;
2257  SCIP_Bool onlynonfixedsign;
2258  int c;
2259 
2260  assert(scip != NULL);
2261  assert(conshdlr != NULL);
2262  assert(conss != NULL || nconss == 0);
2263 
2264  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2265  assert(conshdlrdata != NULL);
2266 
2267  *nnotify = 0;
2268 
2269  onlynonfixedsign = conshdlrdata->preferzerobranch == 3;
2270 
2271  do
2272  {
2273  for( c = 0; c < nconss; ++c )
2274  {
2275  assert(conss[c] != NULL); /*lint !e613*/
2276 
2277  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
2278  assert(consdata != NULL);
2279 
2280  SCIPdebugMsg(scip, "cons <%s> violation: %g %g\n", SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol); /*lint !e613*/
2281 
2282  /* skip variables which sign is already fixed, if we are only interested in variables with unfixed sign here */
2283  if( onlynonfixedsign &&
2284  ( !SCIPisLT(scip, SCIPvarGetLbLocal(consdata->x), -consdata->xoffset) ||
2285  !SCIPisGT(scip, SCIPvarGetUbLocal(consdata->x), consdata->xoffset)) )
2286  continue;
2287 
2288  /* if the value of x lies in a concave range (i.e., where a secant approximation is used), then register x as branching variable */
2289  if( (SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) && (SCIPisInfinity(scip, -SCIPvarGetLbLocal(consdata->x)) || SCIPgetSolVal(scip, sol, consdata->x) + consdata->xoffset <= -consdata->root * (SCIPvarGetLbLocal(consdata->x) + consdata->xoffset))) ||
2290  ( SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && (SCIPisInfinity(scip, SCIPvarGetUbLocal(consdata->x)) || SCIPgetSolVal(scip, sol, consdata->x) + consdata->xoffset >= -consdata->root * (SCIPvarGetUbLocal(consdata->x) + consdata->xoffset))) )
2291  {
2292  /* domain propagation should have removed constraints with fixed x, at least for violated constraints */
2293  assert(!SCIPisRelEQ(scip, SCIPvarGetLbLocal(consdata->x), SCIPvarGetUbLocal(consdata->x)));
2294 
2295  SCIPdebugMsg(scip, "register var <%s> in cons <%s> with violation %g %g\n", SCIPvarGetName(consdata->x), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol); /*lint !e613*/
2296  SCIP_CALL( SCIPaddExternBranchCand(scip, consdata->x, MAX(consdata->lhsviol, consdata->rhsviol), proposeBranchingPoint(scip, conss[c], sol, conshdlrdata->preferzerobranch, conshdlrdata->branchminconverror)) ); /*lint !e613*/
2297  ++*nnotify;
2298  }
2299  }
2300 
2301  if( onlynonfixedsign && *nnotify == 0 )
2302  {
2303  /* if we could not a variable in a violated constraint which sign is not already fixed, do another round where we consider all variables again */
2304  onlynonfixedsign = FALSE;
2305  continue;
2306  }
2307  break;
2308  }
2309  while( TRUE ); /*lint !e506 */
2310 
2311  return SCIP_OKAY; /*lint !e438*/
2312 }
2313 
2314 /** registers a variable from a violated constraint as branching candidate that has a large absolute value in the relaxation */
2315 static
2317  SCIP* scip, /**< SCIP data structure */
2318  SCIP_CONS** conss, /**< constraints */
2319  int nconss, /**< number of constraints */
2320  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
2321  SCIP_VAR** brvar /**< buffer to store branching variable */
2322  )
2323 {
2324  SCIP_CONSDATA* consdata;
2325  SCIP_Real val;
2326  SCIP_Real brvarval;
2327  int c;
2328 
2329  assert(scip != NULL);
2330  assert(conss != NULL || nconss == 0);
2331 
2332  *brvar = NULL;
2333  brvarval = -1.0;
2334 
2335  for( c = 0; c < nconss; ++c )
2336  {
2337  assert(conss != NULL);
2338  consdata = SCIPconsGetData(conss[c]);
2339  assert(consdata != NULL);
2340 
2341  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
2342  continue;
2343 
2344  val = SCIPgetSolVal(scip, sol, consdata->x) + consdata->xoffset;
2345  if( REALABS(val) > brvarval )
2346  {
2347  brvarval = ABS(val);
2348  *brvar = consdata->x;
2349  }
2350  }
2351 
2352  if( *brvar != NULL )
2353  {
2354  SCIP_CALL( SCIPaddExternBranchCand(scip, *brvar, brvarval, SCIP_INVALID) );
2355  }
2356 
2357  return SCIP_OKAY;
2358 }
2359 
2360 /** resolves a propagation on the given variable by supplying the variables needed for applying the corresponding
2361  * propagation rule (see propagateCons()):
2362  * see cons_varbound
2363  */
2364 static
2366  SCIP* scip, /**< SCIP data structure */
2367  SCIP_CONS* cons, /**< constraint that inferred the bound change */
2368  SCIP_VAR* infervar, /**< variable that was deduced */
2369  PROPRULE proprule, /**< propagation rule that deduced the bound change */
2370  SCIP_BOUNDTYPE boundtype, /**< the type of the changed bound (lower or upper bound) */
2371  SCIP_BDCHGIDX* bdchgidx /**< bound change index (time stamp of bound change), or NULL for current time */
2372  )
2373 {
2374  SCIP_CONSDATA* consdata;
2375 
2376  assert(scip != NULL);
2377  assert(cons != NULL);
2378  assert(infervar != NULL);
2379 
2380  consdata = SCIPconsGetData(cons);
2381  assert(consdata != NULL);
2382  assert(consdata->zcoef != 0.0);
2383 
2384  switch( proprule )
2385  {
2386  case PROPRULE_1:
2387  /* lhs <= sign(x+offset)|x+offset|^n + c*z: left hand side and bounds on z -> lower bound on x */
2388  assert(infervar == consdata->x);
2389  assert(boundtype == SCIP_BOUNDTYPE_LOWER);
2390  assert(!SCIPisInfinity(scip, -consdata->lhs));
2391  if( consdata->zcoef > 0.0 )
2392  {
2393  SCIP_CALL( SCIPaddConflictUb(scip, consdata->z, bdchgidx) );
2394  }
2395  else
2396  {
2397  SCIP_CALL( SCIPaddConflictLb(scip, consdata->z, bdchgidx) );
2398  }
2399  break;
2400 
2401  case PROPRULE_2:
2402  /* lhs <= sign(x+offset)|x+offset|^n + c*z: left hand side and upper bound on x -> bound on z */
2403  assert(infervar == consdata->z);
2404  assert(!SCIPisInfinity(scip, -consdata->lhs));
2405  SCIP_CALL( SCIPaddConflictUb(scip, consdata->x, bdchgidx) );
2406  break;
2407 
2408  case PROPRULE_3:
2409  /* sign(x+offset)|x+offset|^n + c*z <= rhs: right hand side and bounds on z -> upper bound on x */
2410  assert(infervar == consdata->x);
2411  assert(boundtype == SCIP_BOUNDTYPE_UPPER);
2412  assert(!SCIPisInfinity(scip, consdata->rhs));
2413  if( consdata->zcoef > 0.0 )
2414  {
2415  SCIP_CALL( SCIPaddConflictLb(scip, consdata->z, bdchgidx) );
2416  }
2417  else
2418  {
2419  SCIP_CALL( SCIPaddConflictUb(scip, consdata->z, bdchgidx) );
2420  }
2421  break;
2422 
2423  case PROPRULE_4:
2424  /* sign(x+offset)|x+offset|^n + c*z <= rhs: right hand side and lower bound on x -> bound on z */
2425  assert(infervar == consdata->z);
2426  assert(!SCIPisInfinity(scip, consdata->rhs));
2427  SCIP_CALL( SCIPaddConflictLb(scip, consdata->x, bdchgidx) );
2428  break;
2429 
2430  case PROPRULE_INVALID:
2431  default:
2432  SCIPerrorMessage("invalid inference information %d in absolute power constraint <%s>\n", proprule, SCIPconsGetName(cons));
2433  return SCIP_INVALIDDATA;
2434  }
2435 
2436  return SCIP_OKAY;
2437 }
2438 
2439 /** analyze infeasibility */
2440 static
2442  SCIP* scip, /**< SCIP data structure */
2443  SCIP_CONS* cons, /**< variable bound constraint */
2444  SCIP_VAR* infervar, /**< variable that was deduced */
2445  PROPRULE proprule, /**< propagation rule that deduced the bound change */
2446  SCIP_BOUNDTYPE boundtype /**< the type of the changed bound (lower or upper bound) */
2447  )
2448 {
2449  /* conflict analysis can only be applied in solving stage and if it turned on */
2451  return SCIP_OKAY;
2452 
2453  /* initialize conflict analysis, and add all variables of infeasible constraint to conflict candidate queue */
2455 
2456  /* add the bound which got violated */
2457  if( boundtype == SCIP_BOUNDTYPE_LOWER )
2458  {
2459  SCIP_CALL( SCIPaddConflictUb(scip, infervar, NULL) );
2460  }
2461  else
2462  {
2463  assert(boundtype == SCIP_BOUNDTYPE_UPPER);
2464  SCIP_CALL( SCIPaddConflictLb(scip, infervar, NULL) );
2465  }
2466 
2467  /* add the reason for the violated of the bound */
2468  SCIP_CALL( resolvePropagation(scip, cons, infervar, proprule, boundtype, NULL) );
2469 
2470  /* analyze the conflict */
2471  SCIP_CALL( SCIPanalyzeConflictCons(scip, cons, NULL) );
2472 
2473  return SCIP_OKAY;
2474 }
2475 
2476 /** propagation method for absolute power constraint
2477  * SCIPinferVarXbCons to allow for repropagation
2478  */
2479 static
2481  SCIP* scip, /**< SCIP data structure */
2482  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2483  SCIP_CONS* cons, /**< variable bound constraint */
2484  SCIP_Bool canaddcons, /**< are we allowed to add a linear constraint when enforcing bounds for a multiaggregated variable? */
2485  SCIP_Bool* cutoff, /**< pointer to store whether the node can be cut off */
2486  int* nchgbds, /**< pointer to count number of bound changes */
2487  int* naddconss /**< pointer to count number of added constraints */
2488  )
2489 {
2490  SCIP_CONSDATA* consdata;
2491  SCIP_Real xlb;
2492  SCIP_Real xub;
2493  SCIP_Real zlb;
2494  SCIP_Real zub;
2495  SCIP_Real newlb;
2496  SCIP_Real newub;
2497  SCIP_Bool tightened;
2498  SCIP_Bool tightenedround;
2499  SCIP_Real minact;
2500  SCIP_Real maxact;
2501 
2502  assert(conshdlr != NULL);
2503  assert(cutoff != NULL);
2504  assert(nchgbds != NULL);
2505  assert(naddconss != NULL);
2506 
2507  consdata = SCIPconsGetData(cons);
2508  assert(consdata != NULL);
2509 
2510  SCIPdebugMsg(scip, "propagating absolute power constraint <%s>\n", SCIPconsGetName(cons));
2511 
2512  *cutoff = FALSE;
2513 
2514  /* get current bounds of variables */
2515  xlb = SCIPvarGetLbLocal(consdata->x);
2516  xub = SCIPvarGetUbLocal(consdata->x);
2517  zlb = SCIPvarGetLbLocal(consdata->z);
2518  zub = SCIPvarGetUbLocal(consdata->z);
2519 
2520  /* if some bound is not tightened, tighten bounds of variables as long as possible */
2521  tightenedround = SCIPconsIsMarkedPropagate(cons);
2522  while( tightenedround )
2523  {
2524  tightenedround = FALSE;
2525 
2526  /* propagate left hand side inequality: lhs <= (x+offset)*|x+offset|^n + c*z */
2527  if( !SCIPisInfinity(scip, -consdata->lhs) )
2528  {
2529  assert(!*cutoff);
2530 
2531  /* propagate bounds on x (if not multiaggregated):
2532  * (1) left hand side and bounds on z -> lower bound on x
2533  */
2534  if( SCIPvarIsActive(SCIPvarGetProbvar(consdata->x)) && (!SCIPisFeasEQ(scip, zlb, zub) || !SCIPisInfinity(scip, REALABS(zlb))) )
2535  {
2536  /* if z is fixed, first compute new lower bound on x without tolerances
2537  * if that is feasible, project new lower bound onto current bounds
2538  * otherwise, recompute with tolerances and continue as usual
2539  * do this only if variable is not essentially fixed to value of infinity
2540  */
2541  if( SCIPisFeasEQ(scip, zlb, zub) && !SCIPisInfinity(scip, zub) )
2542  {
2543  assert(!SCIPisInfinity(scip, -zlb));
2544 
2545  newlb = consdata->lhs - consdata->zcoef * (consdata->zcoef > 0.0 ? zub : zlb);
2546 
2547  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2548  if( consdata->exponent == 2.0 )
2549  newlb = SIGN(newlb) * sqrt(ABS(newlb));
2550  else
2551  newlb = SIGN(newlb) * pow(ABS(newlb), 1.0/consdata->exponent);
2552  newlb -= consdata->xoffset;
2553 
2554  if( SCIPisFeasGT(scip, newlb, xub) )
2555  {
2556  /* if new lower bound for x would yield cutoff, recompute with tolerances */
2557  newlb = consdata->lhs - PROPSIDETOL - consdata->zcoef * (consdata->zcoef > 0.0 ? (zub + PROPVARTOL) : (zlb - PROPVARTOL));
2558 
2559  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2560  if( consdata->exponent == 2.0 )
2561  newlb = SIGN(newlb) * sqrt(ABS(newlb));
2562  else
2563  newlb = SIGN(newlb) * pow(ABS(newlb), 1.0/consdata->exponent);
2564  newlb -= consdata->xoffset;
2565  }
2566  else
2567  {
2568  /* project new lower bound onto current bounds */
2569  newlb = MIN(newlb, xub);
2570  }
2571  }
2572  else
2573  {
2574  if( consdata->zcoef > 0.0 )
2575  {
2576  if( !SCIPisInfinity(scip, zub) )
2577  newlb = consdata->lhs - PROPSIDETOL - consdata->zcoef * (zub + PROPVARTOL);
2578  else
2579  newlb = -SCIPinfinity(scip);
2580  }
2581  else
2582  {
2583  if( !SCIPisInfinity(scip, -zlb) )
2584  newlb = consdata->lhs - PROPSIDETOL - consdata->zcoef * (zlb - PROPVARTOL);
2585  else
2586  newlb = -SCIPinfinity(scip);
2587  }
2588 
2589  if( !SCIPisInfinity(scip, -newlb) )
2590  {
2591  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2592  if( consdata->exponent == 2.0 )
2593  newlb = SIGN(newlb) * sqrt(ABS(newlb));
2594  else
2595  newlb = SIGN(newlb) * pow(ABS(newlb), 1.0/consdata->exponent);
2596  newlb -= consdata->xoffset;
2597  }
2598  }
2599 
2600  if( SCIPisInfinity(scip, newlb) )
2601  {
2602  /* we cannot fix a variable to +infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2603  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->x), xlb, xub, newlb, xub);
2604 
2605  *cutoff = TRUE;
2606 
2607  /* analyze infeasibility */
2608  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_1, SCIP_BOUNDTYPE_LOWER) );
2609  break;
2610  }
2611 
2612  if( !SCIPisInfinity(scip, -newlb) )
2613  {
2614  if( SCIPisLbBetter(scip, newlb, xlb, xub) )
2615  {
2616  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2617  SCIPvarGetName(consdata->x), xlb, xub, newlb, xub);
2618  SCIP_CALL( SCIPinferVarLbCons(scip, consdata->x, newlb, cons, (int)PROPRULE_1, FALSE, cutoff, &tightened) );
2619 
2620  if( *cutoff )
2621  {
2622  assert(SCIPisInfinity(scip, newlb) || SCIPisGT(scip, newlb, SCIPvarGetUbLocal(consdata->x)));
2623 
2624  /* analyze infeasibility */
2625  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_1, SCIP_BOUNDTYPE_LOWER) );
2626  break;
2627  }
2628 
2629  if( tightened )
2630  {
2631  tightenedround = TRUE;
2632  (*nchgbds)++;
2633  }
2634  xlb = SCIPvarGetLbLocal(consdata->x);
2635  }
2636  }
2637  }
2638 
2639  assert(!*cutoff);
2640 
2641  /* propagate bounds on z:
2642  * (2) left hand side and upper bound on x -> bound on z
2643  */
2644  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR && !SCIPisInfinity(scip, xub) ) /* cannot change bounds of multaggr vars */
2645  {
2646  SCIP_Real newbd;
2647 
2648  /* if x is fixed, first compute new bound on z without tolerances
2649  * if that is feasible, project new bound onto current bounds
2650  * otherwise, recompute with tolerances and continue as usual
2651  */
2652  if( SCIPisFeasEQ(scip, xlb, xub) )
2653  {
2654  newbd = xub + consdata->xoffset;
2655  newbd = consdata->lhs - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2656  newbd /= consdata->zcoef;
2657 
2658  if( SCIPisInfinity(scip, newbd) )
2659  newbd = SCIPinfinity(scip);
2660  else if( SCIPisInfinity(scip, -newbd) )
2661  newbd = -SCIPinfinity(scip);
2662 
2663  if( (consdata->zcoef > 0.0 && SCIPisFeasGT(scip, newbd, zub)) || (consdata->zcoef < 0.0 && SCIPisFeasLT(scip, newbd, zlb)) )
2664  {
2665  /* if infeasible, recompute with tolerances */
2666  newbd = xub + PROPVARTOL + consdata->xoffset;
2667  newbd = consdata->lhs - PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2668  newbd /= consdata->zcoef;
2669  }
2670  else
2671  {
2672  /* project onto current bounds of z */
2673  newbd = MIN(zub, MAX(zlb, newbd) );
2674  }
2675  }
2676  else
2677  {
2678  newbd = xub + PROPVARTOL + consdata->xoffset;
2679  newbd = consdata->lhs - PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2680  newbd /= consdata->zcoef;
2681  }
2682 
2683  if( consdata->zcoef > 0.0 )
2684  {
2685  newlb = newbd;
2686 
2687  if( SCIPisInfinity(scip, newlb) )
2688  {
2689  /* we cannot fix a variable to +infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2690  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
2691 
2692  *cutoff = TRUE;
2693 
2694  /* analyze infeasibility */
2695  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_LOWER) );
2696  break;
2697  }
2698 
2699  if( SCIPisLbBetter(scip, newlb, zlb, zub) )
2700  {
2701  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2702  SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
2703  SCIP_CALL( SCIPinferVarLbCons(scip, consdata->z, newlb, cons, (int)PROPRULE_2, FALSE, cutoff, &tightened) );
2704 
2705  if( *cutoff )
2706  {
2707  assert(SCIPisInfinity(scip, newlb) || SCIPisGT(scip, newlb, SCIPvarGetUbLocal(consdata->z)));
2708 
2709  /* analyze infeasibility */
2710  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_LOWER) );
2711  break;
2712  }
2713 
2714  if( tightened )
2715  {
2716  tightenedround = TRUE;
2717  (*nchgbds)++;
2718  }
2719  zlb = SCIPvarGetLbLocal(consdata->z);
2720  }
2721  }
2722  else
2723  {
2724  newub = newbd;
2725 
2726  if( SCIPisInfinity(scip, -newub) )
2727  {
2728  /* we cannot fix a variable to -infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2729  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
2730 
2731  *cutoff = TRUE;
2732 
2733  /* analyze infeasibility */
2734  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_UPPER) );
2735  break;
2736  }
2737 
2738  if( SCIPisUbBetter(scip, newub, zlb, zub) )
2739  {
2740  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2741  SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
2742  SCIP_CALL( SCIPinferVarUbCons(scip, consdata->z, newub, cons, (int)PROPRULE_2, FALSE, cutoff, &tightened) );
2743 
2744  if( *cutoff )
2745  {
2746  assert(SCIPisInfinity(scip, -newub) || SCIPisLT(scip, newub, SCIPvarGetLbLocal(consdata->z)));
2747 
2748  /* analyze infeasibility */
2749  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_2, SCIP_BOUNDTYPE_UPPER) );
2750  break;
2751  }
2752 
2753  if( tightened )
2754  {
2755  tightenedround = TRUE;
2756  (*nchgbds)++;
2757  }
2758  zub = SCIPvarGetUbLocal(consdata->z);
2759  }
2760  }
2761  }
2762  }
2763 
2764  assert(!*cutoff);
2765 
2766  /* propagate right hand side inequality: sign(x+offset)|x+offset|^n + c*z <= rhs */
2767  if( !SCIPisInfinity(scip, consdata->rhs) )
2768  {
2769  /* propagate bounds on x:
2770  * (3) right hand side and bounds on z -> upper bound on x
2771  */
2772  if( SCIPvarIsActive(SCIPvarGetProbvar(consdata->x)) && (!SCIPisFeasEQ(scip, zlb, zub) || !SCIPisInfinity(scip, REALABS(zlb))) ) /* cannot change bounds of multaggr or fixed vars */
2773  {
2774  /* if z is fixed, first compute new upper bound on x without tolerances
2775  * if that is feasible, project new upper bound onto current bounds
2776  * otherwise, recompute with tolerances and continue as usual
2777  * do this only if variable is not essentially fixed to value of infinity
2778  */
2779  if( SCIPisFeasEQ(scip, zlb, zub) && !SCIPisInfinity(scip, zub) )
2780  {
2781  assert(!SCIPisInfinity(scip, -zlb));
2782 
2783  newub = consdata->rhs - consdata->zcoef * (consdata->zcoef > 0.0 ? zlb : zub);
2784 
2785  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2786  if( consdata->exponent == 2.0 )
2787  newub = SIGN(newub) * sqrt(ABS(newub));
2788  else
2789  newub = SIGN(newub) * pow(ABS(newub), 1.0/consdata->exponent);
2790  newub -= consdata->xoffset;
2791 
2792  if( SCIPisFeasLT(scip, newub, xlb) )
2793  {
2794  /* if new lower bound for x would yield cutoff, recompute with tolerances */
2795  newub = consdata->rhs + PROPSIDETOL - consdata->zcoef * (consdata->zcoef > 0.0 ? (zlb - PROPVARTOL) : (zub + PROPVARTOL));
2796 
2797  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2798  if( consdata->exponent == 2.0 )
2799  newub = SIGN(newub) * sqrt(ABS(newub));
2800  else
2801  newub = SIGN(newub) * pow(ABS(newub), 1.0/consdata->exponent);
2802  newub -= consdata->xoffset;
2803  }
2804  else
2805  {
2806  /* project new upper bound onto current bounds */
2807  newub = MAX(newub, xlb);
2808  }
2809  }
2810  else
2811  {
2812  if( consdata->zcoef > 0.0 )
2813  {
2814  if( !SCIPisInfinity(scip, -zlb) )
2815  newub = consdata->rhs + PROPSIDETOL - consdata->zcoef * (zlb - PROPVARTOL);
2816  else
2817  newub = SCIPinfinity(scip);
2818  }
2819  else
2820  {
2821  if( !SCIPisInfinity(scip, zub) )
2822  newub = consdata->rhs + PROPSIDETOL - consdata->zcoef * (zub + PROPVARTOL);
2823  else
2824  newub = SCIPinfinity(scip);
2825  }
2826  if( !SCIPisInfinity(scip, newub) )
2827  {
2828  /* invert sign(x+offset)|x+offset|^(n-1) = y -> x = sign(y)|y|^(1/n) - offset */
2829  if( consdata->exponent == 2.0 )
2830  newub = SIGN(newub) * sqrt(ABS(newub));
2831  else
2832  newub = SIGN(newub) * pow(ABS(newub), 1.0/consdata->exponent);
2833  newub -= consdata->xoffset;
2834  }
2835  }
2836 
2837  if( SCIPisInfinity(scip, -newub) )
2838  {
2839  /* we cannot fix a variable to -infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2840  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->x), xlb, xub, xlb, newub);
2841 
2842  *cutoff = TRUE;
2843 
2844  /* analyze infeasibility */
2845  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_3, SCIP_BOUNDTYPE_UPPER) );
2846  break;
2847  }
2848 
2849  if( !SCIPisInfinity(scip, newub) )
2850  {
2851  if( SCIPisUbBetter(scip, newub, xlb, xub) )
2852  {
2853  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2854  SCIPvarGetName(consdata->x), xlb, xub, xlb, newub);
2855  SCIP_CALL( SCIPinferVarUbCons(scip, consdata->x, newub, cons, (int)PROPRULE_3, FALSE, cutoff, &tightened) );
2856 
2857  if( *cutoff )
2858  {
2859  assert(SCIPisInfinity(scip, -newub) || SCIPisLT(scip, newub, SCIPvarGetLbLocal(consdata->x)));
2860 
2861  /* analyze infeasibility */
2862  SCIP_CALL( analyzeConflict(scip, cons, consdata->x, PROPRULE_3, SCIP_BOUNDTYPE_UPPER) );
2863  break;
2864  }
2865 
2866  if( tightened )
2867  {
2868  tightenedround = TRUE;
2869  (*nchgbds)++;
2870  }
2871  xub = SCIPvarGetUbLocal(consdata->x);
2872  }
2873  }
2874  }
2875 
2876  assert(!*cutoff);
2877 
2878  /* propagate bounds on z:
2879  * (4) right hand side and lower bound on x -> bound on z
2880  */
2881  if( SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR && !SCIPisInfinity(scip, -xlb) ) /* cannot change bounds of multaggr vars */
2882  {
2883  SCIP_Real newbd;
2884 
2885  /* if x is fixed, first compute new bound on z without tolerances
2886  * if that is feasible, project new bound onto current bounds
2887  * otherwise, recompute with tolerances and continue as usual
2888  */
2889  if( SCIPisFeasEQ(scip, xlb, xub) )
2890  {
2891  newbd = xlb + consdata->xoffset;
2892  newbd = consdata->rhs - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2893  newbd /= consdata->zcoef;
2894 
2895  if( SCIPisInfinity(scip, newbd) )
2896  newbd = SCIPinfinity(scip);
2897  else if( SCIPisInfinity(scip, -newbd) )
2898  newbd = -SCIPinfinity(scip);
2899 
2900  if( (consdata->zcoef > 0.0 && SCIPisFeasLT(scip, newbd, zlb)) || (consdata->zcoef < 0.0 && SCIPisFeasGT(scip, newbd, zub)) )
2901  {
2902  /* if infeasible, recompute with tolerances */
2903  newbd = xlb - PROPVARTOL + consdata->xoffset;
2904  newbd = consdata->rhs + PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2905  newbd /= consdata->zcoef;
2906  }
2907  else
2908  {
2909  /* project onto current bounds of z */
2910  newbd = MIN(zub, MAX(zlb, newbd) );
2911  }
2912  }
2913  else
2914  {
2915  newbd = xlb - PROPVARTOL + consdata->xoffset;
2916  newbd = consdata->rhs + PROPSIDETOL - SIGN(newbd) * consdata->power(REALABS(newbd), consdata->exponent);
2917  newbd /= consdata->zcoef;
2918  }
2919 
2920  if( consdata->zcoef > 0.0 )
2921  {
2922  newub = newbd;
2923 
2924  if( SCIPisInfinity(scip, -newub) )
2925  {
2926  /* we cannot fix a variable to -infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2927  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
2928 
2929  *cutoff = TRUE;
2930 
2931  /* analyze infeasibility */
2932  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_UPPER) );
2933  break;
2934  }
2935 
2936  if( SCIPisUbBetter(scip, newub, zlb, zub) )
2937  {
2938  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2939  SCIPvarGetName(consdata->z), zlb, zub, zlb, newub);
2940  SCIP_CALL( SCIPinferVarUbCons(scip, consdata->z, newub, cons, (int)PROPRULE_4, FALSE, cutoff, &tightened) );
2941 
2942  if( *cutoff )
2943  {
2944  assert(SCIPisInfinity(scip, -newub) || SCIPisLT(scip, newub, SCIPvarGetLbLocal(consdata->z)));
2945 
2946  /* analyze infeasibility */
2947  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_UPPER) );
2948  break;
2949  }
2950 
2951  if( tightened )
2952  {
2953  tightenedround = TRUE;
2954  (*nchgbds)++;
2955  }
2956  zub = SCIPvarGetUbLocal(consdata->z);
2957  }
2958  }
2959  else
2960  {
2961  newlb = newbd;
2962 
2963  if( SCIPisInfinity(scip, newlb) )
2964  {
2965  /* we cannot fix a variable to +infinity, so let's report cutoff (there is no solution within SCIPs limitations to infinity) */
2966  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g] -> cutoff\n", SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
2967 
2968  *cutoff = TRUE;
2969 
2970  /* analyze infeasibility */
2971  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_LOWER) );
2972  break;
2973  }
2974 
2975  if( SCIPisLbBetter(scip, newlb, zlb, zub) )
2976  {
2977  SCIPdebugMsg(scip, " -> tighten <%s>[%.15g,%.15g] -> [%.15g,%.15g]\n",
2978  SCIPvarGetName(consdata->z), zlb, zub, newlb, zub);
2979  SCIP_CALL( SCIPinferVarLbCons(scip, consdata->z, newlb, cons, (int)PROPRULE_4, FALSE, cutoff, &tightened) );
2980 
2981  if( *cutoff )
2982  {
2983  assert(SCIPisInfinity(scip, newlb) || SCIPisGT(scip, newlb, SCIPvarGetUbLocal(consdata->z)));
2984 
2985  /* analyze infeasibility */
2986  SCIP_CALL( analyzeConflict(scip, cons, consdata->z, PROPRULE_4, SCIP_BOUNDTYPE_LOWER) );
2987  break;
2988  }
2989 
2990  if( tightened )
2991  {
2992  tightenedround = TRUE;
2993  (*nchgbds)++;
2994  }
2995  zlb = SCIPvarGetLbLocal(consdata->z);
2996  }
2997  }
2998  }
2999  }
3000 
3001  assert(!*cutoff);
3002  }
3003 
3004  /* mark the constraint propagated */
3005  SCIP_CALL( SCIPunmarkConsPropagate(scip, cons) );
3006 
3007  if( *cutoff )
3008  return SCIP_OKAY;
3009 
3010  /* check for redundancy */
3011  if( !SCIPisInfinity(scip, -xlb) && !SCIPisInfinity(scip, consdata->zcoef > 0.0 ? -zlb : zub) )
3012  minact = SIGN(xlb + consdata->xoffset) * consdata->power(REALABS(xlb + consdata->xoffset), consdata->exponent) + consdata->zcoef * (consdata->zcoef > 0.0 ? zlb : zub);
3013  else
3014  minact = -SCIPinfinity(scip);
3015 
3016  if( !SCIPisInfinity(scip, xub) && !SCIPisInfinity(scip, consdata->zcoef > 0.0 ? zub : -zlb) )
3017  maxact = SIGN(xub + consdata->xoffset) * consdata->power(REALABS(xub + consdata->xoffset), consdata->exponent) + consdata->zcoef * (consdata->zcoef > 0.0 ? zub : zlb);
3018  else
3019  maxact = SCIPinfinity(scip);
3020 
3021  if( SCIPisFeasGE(scip, minact, consdata->lhs) && SCIPisFeasLE(scip, maxact, consdata->rhs) )
3022  {
3023  SCIPdebugMsg(scip, "absolute power constraint <%s> is redundant: <%s>[%.15g,%.15g], <%s>[%.15g,%.15g]\n",
3024  SCIPconsGetName(cons),
3025  SCIPvarGetName(consdata->x), SCIPvarGetLbLocal(consdata->x), SCIPvarGetUbLocal(consdata->x),
3026  SCIPvarGetName(consdata->z), SCIPvarGetLbLocal(consdata->z), SCIPvarGetUbLocal(consdata->z));
3027 
3028  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
3029 
3030  return SCIP_OKAY;
3031  }
3032 
3033  /* delete constraint if x has been fixed */
3034  if( SCIPisRelEQ(scip, xlb, xub) && (SCIPvarIsActive(consdata->z) || canaddcons) )
3035  {
3036  SCIP_RESULT tightenresult;
3037  SCIP_INTERVAL xbnds;
3038  SCIP_INTERVAL zbnds;
3039 
3040  SCIPdebugMsg(scip, "x-variable in constraint <%s> is fixed: x = <%s>[%.15g,%.15g], z = <%s>[%.15g,%.15g]\n",
3041  SCIPconsGetName(cons), SCIPvarGetName(consdata->x), xlb, xub, SCIPvarGetName(consdata->z), zlb, zub);
3042 
3043  SCIPintervalSetBounds(&xbnds, MIN(xlb, xub), MAX(xlb, xub));
3044  computeBoundsZ(scip, cons, xbnds, &zbnds);
3045 
3046  /* in difference to the loop above, here we enforce a possible bound tightening on z, and may add a linear cons if z is multiaggregated */
3047  SCIP_CALL( tightenBounds(scip, consdata->z, zbnds, TRUE, cons, &tightenresult, nchgbds, nchgbds, naddconss) );
3048  if( tightenresult == SCIP_CUTOFF )
3049  *cutoff = TRUE;
3050 
3051  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
3052 
3053  return SCIP_OKAY;
3054  }
3055 
3056  /* delete constraint if z has been fixed */
3057  if( SCIPisRelEQ(scip, zlb, zub) && (SCIPvarIsActive(consdata->x) || canaddcons) )
3058  {
3059  SCIP_RESULT tightenresult;
3060  SCIP_INTERVAL xbnds;
3061  SCIP_INTERVAL zbnds;
3062 
3063  SCIPdebugMsg(scip, "z-variable in constraint <%s> is fixed: x = <%s>[%.15g,%.15g], z = <%s>[%.15g,%.15g]\n",
3064  SCIPconsGetName(cons), SCIPvarGetName(consdata->x), xlb, xub, SCIPvarGetName(consdata->z), zlb, zub);
3065 
3066  SCIPintervalSetBounds(&zbnds, MIN(zlb, zub), MAX(zlb, zub));
3067  computeBoundsX(scip, cons, zbnds, &xbnds);
3068 
3069  /* in difference to the loop above, here we enforce a possible bound tightening on x, and may add a linear cons if x is multiaggregated */
3070  SCIP_CALL( tightenBounds(scip, consdata->x, xbnds, TRUE, cons, &tightenresult, nchgbds, nchgbds, naddconss) );
3071  if( tightenresult == SCIP_CUTOFF )
3072  *cutoff = TRUE;
3073 
3074  SCIP_CALL( SCIPdelConsLocal(scip, cons) );
3075 
3076  return SCIP_OKAY;
3077  }
3078 
3079  return SCIP_OKAY;
3080 }
3081 
3082 /** notifies SCIP about a variable bound lhs <= x + c*y <= rhs */
3083 static
3085  SCIP* scip, /**< SCIP data structure */
3086  SCIP_CONS* cons, /**< absolute power constraint this variable bound is derived form */
3087  SCIP_Bool addcons, /**< should the variable bound be added as constraint to SCIP? */
3088  SCIP_VAR* var, /**< variable x for which we want to add a variable bound */
3089  SCIP_VAR* vbdvar, /**< variable y which makes the bound a variable bound */
3090  SCIP_Real vbdcoef, /**< coefficient c of bounding variable vbdvar */
3091  SCIP_Real lhs, /**< left hand side of varbound constraint */
3092  SCIP_Real rhs, /**< right hand side of varbound constraint */
3093  SCIP_Bool* infeas, /**< pointer to store whether an infeasibility was detected */
3094  int* nbdchgs, /**< pointer where to add number of performed bound changes */
3095  int* naddconss /**< pointer where to add number of added constraints */
3096  )
3097 {
3098  int nbdchgs_local;
3099 
3100  assert(scip != NULL);
3101  assert(cons != NULL);
3102  assert(var != NULL);
3103  assert(vbdvar != NULL);
3104  assert(!SCIPisZero(scip, vbdcoef));
3105  assert(!SCIPisInfinity(scip, ABS(vbdcoef)));
3106  assert(infeas != NULL);
3107 
3108  *infeas = FALSE;
3109 
3110  /* make sure vbdvar is active, so we can search for it in SCIPvarGetVxbdVars() */
3111  if( !SCIPvarIsActive(vbdvar) )
3112  {
3113  SCIP_Real constant;
3114 
3115  constant = 0.0;
3116  SCIP_CALL( SCIPgetProbvarSum(scip, &vbdvar, &vbdcoef, &constant) );
3117  if( !SCIPvarIsActive(vbdvar) || (vbdcoef == 0.0) )
3118  return SCIP_OKAY;
3119 
3120  if( !SCIPisInfinity(scip, -lhs) )
3121  lhs -= constant;
3122  if( !SCIPisInfinity(scip, rhs) )
3123  rhs -= constant;
3124  }
3125 
3126  /* vbdvar should be a non-fixed binary variable */
3127  assert(SCIPvarIsIntegral(vbdvar));
3128  assert(SCIPisZero(scip, SCIPvarGetLbGlobal(vbdvar)));
3129  assert(SCIPisEQ(scip, SCIPvarGetUbGlobal(vbdvar), 1.0));
3130 
3131  SCIPdebugMsg(scip, "-> %g <= <%s> + %g*<%s> <= %g\n", lhs, SCIPvarGetName(var), vbdcoef, SCIPvarGetName(vbdvar), rhs);
3132 
3133  if( addcons && SCIPvarGetStatus(var) != SCIP_VARSTATUS_MULTAGGR )
3134  {
3135  SCIP_CONS* vbdcons;
3136  char name[SCIP_MAXSTRLEN];
3137 
3138  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_vbnd", SCIPconsGetName(cons));
3139 
3140  SCIP_CALL( SCIPcreateConsVarbound(scip, &vbdcons, name, var, vbdvar, vbdcoef, lhs, rhs,
3142  SCIP_CALL( SCIPaddCons(scip, vbdcons) );
3143  SCIP_CALL( SCIPreleaseCons(scip, &vbdcons) );
3144 
3145  ++*naddconss;
3146 
3147  return SCIP_OKAY;
3148  }
3149 
3150 
3151  if( !SCIPisInfinity(scip, -lhs) )
3152  {
3153  SCIP_CALL( SCIPaddVarVlb(scip, var, vbdvar, -vbdcoef, lhs, infeas, &nbdchgs_local) );
3154  if( *infeas )
3155  return SCIP_OKAY;
3156  *nbdchgs += nbdchgs_local;
3157  }
3158 
3159  if( !SCIPisInfinity(scip, rhs) )
3160  {
3161  SCIP_CALL( SCIPaddVarVub(scip, var, vbdvar, -vbdcoef, rhs, infeas, &nbdchgs_local) );
3162  if( *infeas )
3163  return SCIP_OKAY;
3164  *nbdchgs += nbdchgs_local;
3165  }
3166 
3167  return SCIP_OKAY;
3168 }
3169 
3170 /** propagates varbounds of variables
3171  * Let f(x) = sign(x+offset)|x+offset|^n, f^{-1}(y) = sign(y)|y|^(1/n) - offset.
3172  * Thus, constraint is lhs <= f(x) + c*z <= rhs.
3173  *
3174  * Given a variable bound constraint x <= a*y + b with y a binary variable, one obtains
3175  * y = 0 -> f(x) <= f(b) -> lhs <= f(b) + c*z
3176  * y = 1 -> f(x) <= f(a+b) -> lhs <= f(a+b) + c*z
3177  * => lhs <= f(b) + y * (f(a+b)-f(b)) + c*z
3178  *
3179  * Given a variable bound constraint x >= a*y + b with y a binary variable, one obtains analogously
3180  * f(b) + y * (f(a+b)-f(b)) + c*z <= rhs
3181  *
3182  * Given a variable bound constraint c*z <= a*y + b with y a binary variable, one obtains
3183  * y = 0 -> lhs <= f(x) + b -> x >= f^{-1}(lhs - b)
3184  * y = 1 -> lhs <= f(x) + a+b -> x >= f^{-1}(lhs - (a+b))
3185  * => x >= f^{-1}(lhs - b) + y * (f^{-1}(lhs - (a+b)) - f^{-1}(lhs - b))
3186  *
3187  * Given a variable bound constraint c*z >= a*y + b with y a binary variable, one obtains analogously
3188  * x <= f^{-1}(rhs - b) + y * (f^{-1}(rhs - (a+b)) - f^{-1}(rhs - b))
3189  */
3190 static
3192  SCIP* scip, /**< SCIP data structure */
3193  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3194  SCIP_CONS* cons, /**< absolute power constraint */
3195  SCIP_Bool* infeas, /**< pointer to store whether an infeasibility was detected */
3196  int* nbdchgs, /**< pointer where to add number of performed bound changes */
3197  int* naddconss /**< pointer where to add number of added constraints */
3198  )
3199 {
3200  SCIP_CONSHDLRDATA* conshdlrdata;
3201  SCIP_CONSDATA* consdata;
3202  SCIP_VAR* y;
3203  SCIP_Real a;
3204  SCIP_Real b;
3205  SCIP_Real fb;
3206  SCIP_Real fab;
3207  SCIP_Real vbcoef;
3208  SCIP_Real vbconst;
3209  int i;
3210 
3211  assert(scip != NULL);
3212  assert(conshdlr != NULL);
3213  assert(cons != NULL);
3214  assert(infeas != NULL);
3215  assert(nbdchgs != NULL);
3216  assert(naddconss != NULL);
3217 
3218  *infeas = FALSE;
3219 
3220  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3221  assert(conshdlrdata != NULL);
3222 
3223  consdata = SCIPconsGetData(cons);
3224  assert(consdata != NULL);
3225  assert(consdata->z != NULL);
3226 
3227  /* don't do anything if it looks like we have numerical troubles */
3228  if( SCIPisZero(scip, consdata->zcoef) )
3229  return SCIP_OKAY;
3230 
3231  if( !SCIPisInfinity(scip, -consdata->lhs) )
3232  {
3233  /* propagate varbounds x <= a*y+b onto z
3234  * lhs <= f(b) + y * (f(a+b)-f(b)) + c*z
3235  * -> c*z >= lhs-f(b) + y * (f(b)-f(a+b))
3236  */
3237  for( i = 0; i < SCIPvarGetNVubs(consdata->x); ++i )
3238  {
3239  y = SCIPvarGetVubVars(consdata->x)[i];
3240  a = SCIPvarGetVubCoefs(consdata->x)[i];
3241  b = SCIPvarGetVubConstants(consdata->x)[i];
3242 
3243  /* skip variable bound if y is not integer or its valid values are not {0,1}
3244  * @todo extend to arbitrary integer variables
3245  */
3246  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3247  continue;
3248 
3249  /* skip variable bound if coefficient is very small */
3250  if( SCIPisFeasZero(scip, consdata->power(a, consdata->exponent)) )
3251  continue;
3252 
3253  SCIPdebugMsg(scip, "propagate variable bound <%s> <= %g*<%s> + %g\n", SCIPvarGetName(consdata->x), a, SCIPvarGetName(y), b);
3254 
3255  fb = SIGN( b + consdata->xoffset) * consdata->power( b + consdata->xoffset, consdata->exponent); /* f( b) = sign( b) | b|^n */
3256  fab = SIGN(a+b + consdata->xoffset) * consdata->power(a+b + consdata->xoffset, consdata->exponent); /* f(a+b) = sign(a+b) |a+b|^n */
3257 
3258  vbcoef = (fb - fab) / consdata->zcoef;
3259  vbconst = (consdata->lhs - fb) / consdata->zcoef;
3260 
3261  if( consdata->zcoef > 0.0 )
3262  {
3263  /* add varbound z >= (lhs-f(b))/c + y * (f(b)-f(a+b))/c */
3264  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, vbconst, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3265  }
3266  else
3267  {
3268  /* add varbound z <= (lhs-f(b))/c + y * (f(b)-f(a+b))/c */
3269  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, -SCIPinfinity(scip), vbconst, infeas, nbdchgs, naddconss) );
3270  }
3271  if( *infeas )
3272  return SCIP_OKAY;
3273  }
3274  }
3275 
3276  /* propagate varbounds x >= a*y+b onto z
3277  * f(b) + y * (f(a+b)-f(b)) + c*z <= rhs
3278  * -> c*z <= rhs-f(b) + y * (f(b)-f(a+b))
3279  */
3280  if( !SCIPisInfinity(scip, consdata->rhs) )
3281  {
3282  for( i = 0; i < SCIPvarGetNVlbs(consdata->x); ++i )
3283  {
3284  y = SCIPvarGetVlbVars(consdata->x)[i];
3285  a = SCIPvarGetVlbCoefs(consdata->x)[i];
3286  b = SCIPvarGetVlbConstants(consdata->x)[i];
3287 
3288  /* skip variable bound if y is not integer or its valid values are not {0,1}
3289  * @todo extend to arbitrary integer variables
3290  */
3291  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3292  continue;
3293 
3294  /* skip variable bound if coefficient is very small */
3295  if( SCIPisFeasZero(scip, consdata->power(a, consdata->exponent)) )
3296  continue;
3297 
3298  SCIPdebugMsg(scip, "propagate variable bound <%s> >= %g*<%s> + %g\n", SCIPvarGetName(consdata->x), a, SCIPvarGetName(y), b);
3299 
3300  fb = SIGN( b + consdata->xoffset) * consdata->power( b + consdata->xoffset, consdata->exponent); /* f( b) = sign( b) | b|^n */
3301  fab = SIGN(a+b + consdata->xoffset) * consdata->power(a+b + consdata->xoffset, consdata->exponent); /* f(a+b) = sign(a+b) |a+b|^n */
3302 
3303  vbcoef = (fb - fab) / consdata->zcoef;
3304  vbconst = (consdata->rhs - fb) / consdata->zcoef;
3305 
3306  if( consdata->zcoef > 0.0 )
3307  {
3308  /* add varbound z <= (rhs-f(b))/c + y * (f(b)-f(a+b))/c */
3309  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, -SCIPinfinity(scip), vbconst, infeas, nbdchgs, naddconss) );
3310  }
3311  else
3312  {
3313  /* add varbound z >= (rhs-f(b))/c + y * (f(b)-f(a+b))/c */
3314  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->z, y, -vbcoef, vbconst, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3315  }
3316  if( *infeas )
3317  return SCIP_OKAY;
3318  }
3319  }
3320 
3321  /* propagate variable upper bounds on z onto x
3322  * c*z <= a*y+b -> x >= f^{-1}(lhs - b) + y * (f^{-1}(lhs - (a+b)) - f^{-1}(lhs - b))
3323  * c*z >= a*y+b -> x <= f^{-1}(rhs - b) + y * (f^{-1}(rhs - (a+b)) - f^{-1}(rhs - b))
3324  */
3325  if( (consdata->zcoef > 0.0 && !SCIPisInfinity(scip, -consdata->lhs)) ||
3326  ( consdata->zcoef < 0.0 && !SCIPisInfinity(scip, consdata->rhs)) )
3327  for( i = 0; i < SCIPvarGetNVubs(consdata->z); ++i )
3328  {
3329  y = SCIPvarGetVubVars(consdata->z)[i];
3330  a = SCIPvarGetVubCoefs(consdata->z)[i] * consdata->zcoef;
3331  b = SCIPvarGetVubConstants(consdata->z)[i] * consdata->zcoef;
3332 
3333  SCIPdebugMsg(scip, "propagate variable bound %g*<%s> %c= %g*<%s> + %g\n", consdata->zcoef, SCIPvarGetName(consdata->z), consdata->zcoef > 0 ? '<' : '>', a, SCIPvarGetName(y), b);
3334 
3335  /* skip variable bound if y is not integer or its valid values are not {0,1}
3336  * @todo extend to arbitrary integer variables
3337  */
3338  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3339  continue;
3340 
3341  if( consdata->zcoef > 0.0 )
3342  {
3343  fb = consdata->lhs - b;
3344  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3345  fab = consdata->lhs - (a+b);
3346  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3347  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, fb - consdata->xoffset, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3348  }
3349  else
3350  {
3351  fb = consdata->rhs - b;
3352  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3353  fab = consdata->rhs - (a+b);
3354  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3355  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, -SCIPinfinity(scip), fb - consdata->xoffset, infeas, nbdchgs, naddconss) );
3356  }
3357  if( *infeas )
3358  return SCIP_OKAY;
3359  }
3360 
3361  /* propagate variable lower bounds on z onto x
3362  * c*z <= a*y+b -> x >= f^{-1}(lhs - b) + y * (f^{-1}(lhs - (a+b)) - f^{-1}(lhs - b))
3363  * c*z >= a*y+b -> x <= f^{-1}(rhs - b) + y * (f^{-1}(rhs - (a+b)) - f^{-1}(rhs - b))
3364  */
3365  if( (consdata->zcoef < 0.0 && !SCIPisInfinity(scip, -consdata->lhs)) ||
3366  ( consdata->zcoef > 0.0 && !SCIPisInfinity(scip, consdata->rhs)) )
3367  for( i = 0; i < SCIPvarGetNVlbs(consdata->z); ++i )
3368  {
3369  y = SCIPvarGetVlbVars(consdata->z)[i];
3370  a = SCIPvarGetVlbCoefs(consdata->z)[i] * consdata->zcoef;
3371  b = SCIPvarGetVlbConstants(consdata->z)[i] * consdata->zcoef;
3372 
3373  SCIPdebugMsg(scip, "propagate variable bound %g*<%s> %c= %g*<%s> + %g\n", consdata->zcoef, SCIPvarGetName(consdata->z), consdata->zcoef > 0 ? '>' : '<', a, SCIPvarGetName(y), b);
3374 
3375  /* skip variable bound if y is not integer or its valid values are not {0,1}
3376  * @todo extend to arbitrary integer variables
3377  */
3378  if( !SCIPvarIsBinary(y) || SCIPvarGetLbGlobal(y) > 0.5 || SCIPvarGetUbGlobal(y) < 0.5 )
3379  continue;
3380 
3381  if( consdata->zcoef > 0.0 )
3382  {
3383  fb = consdata->rhs - b;
3384  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3385  fab = consdata->rhs - (a+b);
3386  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3387  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, -SCIPinfinity(scip), fb - consdata->xoffset, infeas, nbdchgs, naddconss) );
3388  }
3389  else
3390  {
3391  fb = consdata->lhs - b;
3392  fb = SIGN(fb) * pow(ABS(fb), 1.0/consdata->exponent);
3393  fab = consdata->lhs - (a+b);
3394  fab = SIGN(fab) * pow(ABS(fab), 1.0/consdata->exponent);
3395  SCIP_CALL( addVarbound(scip, cons, conshdlrdata->addvarboundcons, consdata->x, y, fb - fab, fb - consdata->xoffset, SCIPinfinity(scip), infeas, nbdchgs, naddconss) );
3396  }
3397  if( *infeas )
3398  return SCIP_OKAY;
3399  }
3400 
3401  return SCIP_OKAY;
3402 }
3403 
3404 /** computes linear underestimator for (x+offset)^n + c*z <= rhs by linearization in x
3405  *
3406  * the generated cut is xmul * n * (refpoint+offset)^(n-1) * x + c*z <= rhs + ((n-1)*refpoint-offset) * (refpoint+offset)^(n-1)
3407  */
3408 static
3410  SCIP* scip, /**< SCIP data structure */
3411  SCIP_ROW** row, /**< buffer to store row */
3412  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3413  SCIP_Real refpoint, /**< base point for linearization */
3414  SCIP_Real exponent, /**< exponent n in sign(x)abs(x)^n */
3415  SCIP_Real xoffset, /**< offset of x */
3416  SCIP_Real xmult, /**< multiplier for coefficient of x */
3417  SCIP_Real zcoef, /**< coefficient of z */
3418  SCIP_Real rhs, /**< right hand side */
3419  SCIP_VAR* x, /**< variable x */
3420  SCIP_VAR* z, /**< variable z */
3421  SCIP_Bool islocal /**< whether the cut is valid only locally */
3422  )
3423 {
3424  char rowname[SCIP_MAXSTRLEN];
3425  SCIP_CONSHDLRDATA* conshdlrdata;
3426  SCIP_Real tmp;
3427 
3428  assert(scip != NULL);
3429  assert(!SCIPisFeasNegative(scip, refpoint+xoffset));
3430  assert(!SCIPisInfinity(scip, refpoint));
3431 
3432  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3433  assert(conshdlrdata != NULL);
3434 
3435  if( refpoint < -xoffset )
3436  refpoint = -xoffset;
3437 
3438  tmp = exponent == 2.0 ? refpoint+xoffset : pow(refpoint+xoffset, exponent-1);
3439  if( SCIPisInfinity(scip, tmp) )
3440  {
3441  SCIPdebugMsg(scip, "skip linearization cut because (refpoint+offset)^(exponent-1) > infinity\n");
3442  *row = NULL;
3443  return SCIP_OKAY;
3444  }
3445 
3446  rhs += ((exponent-1)*refpoint-xoffset)*tmp; /* now rhs is the rhs of the cut */
3447  /* do not change the right hand side to a value > infinity (this would trigger an assertion in lp.c) */
3448  if( SCIPisInfinity(scip, rhs) )
3449  {
3450  SCIPdebugMsg(scip, "skip linearization cut because its rhs would be > infinity\n");
3451  *row = NULL;
3452  return SCIP_OKAY;
3453  }
3454 
3455  (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "signpowlinearizecut_%u", ++(conshdlrdata->ncuts));
3456 
3457  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, conshdlr, rowname, -SCIPinfinity(scip), rhs, islocal,
3458  FALSE /* modifiable */, TRUE /* removable */ ) );
3459 
3460  SCIP_CALL( SCIPaddVarToRow(scip, *row, x, xmult*exponent*tmp) );
3461  SCIP_CALL( SCIPaddVarToRow(scip, *row, z, zcoef) );
3462 
3463  return SCIP_OKAY;
3464 }
3465 
3466 /** computes linear underestimator for (x+xoffset)^n + c*z <= rhs by linearization in x
3467  *
3468  * the generated cut is xmul * n * (refpoint+offset)^(n-1) * x + c*z <= rhs + ((n-1)*refpoint-offset) * (refpoint+offset)^(n-1)
3469  * where refpoint is computed by projecting (xref, zref) onto the graph of (x+offset)^n w.r.t. euclidean norm
3470  *
3471  * Thus, the projection is computed by minimizing 1/2(x-xref)^2 + 1/2(((x+offset)^n-rhs)/(-c) - zref)^2.
3472  * I.e., we aim to find a root of
3473  * g(x) = x - xref + n/c (x+offset)^(n-1) (zref - rhs/c) + n/c^2 (x+offset)^(2n-1)
3474  * We do this numerically by executing up to five newton iterations. It is
3475  * g'(x) = 1 + n(n-1)/c (x+offset)^(n-2) (zref - rhs/c) + n(2n-1)/c^2 (x+offset)^(2n-2)
3476  */
3477 static
3479  SCIP* scip, /**< SCIP data structure */
3480  SCIP_ROW** row, /**< buffer to store row */
3481  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3482  SCIP_Real xref, /**< reference point for x */
3483  SCIP_Real zref, /**< reference point for z */
3484  SCIP_Real xmin, /**< minimal value x is allowed to take */
3485  SCIP_Real exponent, /**< exponent n in sign(x+offset)abs(x+offset)^n */
3486  SCIP_Real xoffset, /**< offset of x */
3487  SCIP_Real xmult, /**< multiplier for coefficient of x */
3488  SCIP_Real zcoef, /**< coefficient of z */
3489  SCIP_Real rhs, /**< right hand side */
3490  SCIP_VAR* x, /**< variable x */
3491  SCIP_VAR* z, /**< variable z */
3492  SCIP_Bool islocal /**< whether the cut is valid only locally */
3493  )
3494 {
3495  SCIP_Real tmp;
3496  SCIP_Real xproj;
3497  SCIP_Real gval;
3498  SCIP_Real gderiv;
3499  int iter;
3500 
3501  assert(scip != NULL);
3502  assert(!SCIPisFeasNegative(scip, xref+xoffset));
3503  assert(!SCIPisInfinity(scip, xref));
3504 
3505  if( xref < xmin )
3506  xref = xmin;
3507 
3508  xproj = xref;
3509  iter = 0;
3510  if( exponent == 2.0 )
3511  do
3512  {
3513  tmp = (xproj+xoffset) * (xproj+xoffset);
3514  gval = xproj - xref + 2*(xproj+xoffset) / zcoef * ((tmp-rhs)/zcoef + zref);
3515  if( !SCIPisFeasPositive(scip, ABS(gval)) )
3516  break;
3517 
3518  gderiv = 1 + 6 * tmp / (zcoef*zcoef) + 2 / zcoef * (zref - rhs/zcoef);
3519  xproj -= gval / gderiv;
3520 
3521  }
3522  while( ++iter <= 5 );
3523  else
3524  do
3525  {
3526  tmp = pow(xproj + xoffset, exponent-1);
3527  gval = xproj - xref + exponent / zcoef * (pow(xproj+xoffset, 2*exponent-1)/zcoef + tmp * (zref-rhs/zcoef));
3528  if( !SCIPisFeasPositive(scip, ABS(gval)) )
3529  break;
3530 
3531  gderiv = 1 + exponent / zcoef * ( (2*exponent-1)*tmp*tmp/zcoef + (exponent-1)*pow(xproj+xoffset, exponent-2) * (zref-rhs/zcoef) );
3532  xproj -= gval / gderiv;
3533 
3534  }
3535  while( ++iter <= 5 );
3536 
3537  if( xproj < xmin )
3538  xproj = xmin;
3539 
3540  SCIP_CALL( generateLinearizationCut(scip, row, conshdlr, xproj, exponent, xoffset, xmult, zcoef, rhs, x, z, islocal) );
3541 
3542  return SCIP_OKAY;
3543 }
3544 
3545 /** computes secant underestimator for sign(x+offset)abs(x+offset)^n + c*z <= rhs
3546  *
3547  * the generated cut is slope*xmult*x + c*z <= rhs + (-xlb-offset)^n + slope*xlb,
3548  * where slope = (sign(xub+offset)*abs(xub+offset)^n + (-xlb-offset)^n) / (xub - xlb).
3549  *
3550  * the cut is not generated if the given solution (or the LP solution) would not be cutoff
3551  */
3552 static
3554  SCIP* scip, /**< SCIP data structure */
3555  SCIP_ROW** row, /**< buffer to store row */
3556  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3557  SCIP_SOL* sol, /**< point we want to cut off, or NULL for LP solution */
3558  SCIP_Real xlb, /**< lower bound of x */
3559  SCIP_Real xub, /**< upper bound of x */
3560  SCIP_Real exponent, /**< exponent n in sign(x+offset)abs(x+offset)^n */
3561  SCIP_Real xoffset, /**< offset of x */
3562  DECL_MYPOW ((*mypow)), /**< function to use for computing power */
3563  SCIP_Real xmult, /**< multiplier for coefficient of x */
3564  SCIP_Real zcoef, /**< coefficient of z */
3565  SCIP_Real rhs, /**< right hand side */
3566  SCIP_VAR* x, /**< variable x */
3567  SCIP_VAR* z /**< variable z */
3568  )
3569 {
3570  char rowname[SCIP_MAXSTRLEN];
3571  SCIP_CONSHDLRDATA* conshdlrdata;
3572  SCIP_Real slope, tmp, val;
3573 
3574  assert(scip != NULL);
3575  assert(SCIPisLE(scip, xlb, xub));
3576  assert(!SCIPisPositive(scip, xlb+xoffset));
3577 
3578  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3579  assert(conshdlrdata != NULL);
3580 
3581  /* ignore constraints with fixed x (should be removed soon) */
3582  if( SCIPisRelEQ(scip, xlb, xub) )
3583  return SCIP_OKAY;
3584 
3585  if( xlb > -xoffset )
3586  xlb = -xoffset;
3587 
3588  tmp = mypow(-xlb-xoffset, exponent);
3589  slope = SIGN(xub+xoffset) * mypow(ABS(xub+xoffset), exponent) + tmp;
3590  slope /= xub - xlb;
3591 
3592  /* check if cut would violated solution, check that slope is not above value of infinity */
3593  val = -tmp + slope * (xmult * SCIPgetSolVal(scip, sol, x) - xlb) + zcoef * SCIPgetSolVal(scip, sol, z) - rhs;
3594  if( !SCIPisFeasPositive(scip, val) || SCIPisInfinity(scip, REALABS(slope)) )
3595  {
3596  *row = NULL;
3597  return SCIP_OKAY;
3598  }
3599 
3600  (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "signpowsecantcut_%u", ++(conshdlrdata->nsecantcuts));
3601 
3602  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, conshdlr, rowname, -SCIPinfinity(scip), SCIPinfinity(scip),
3603  SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0 /* local */, FALSE /* modifiable */, TRUE /* removable */ ) );
3604  SCIP_CALL( SCIPaddVarToRow(scip, *row, x, xmult*slope) );
3605  SCIP_CALL( SCIPaddVarToRow(scip, *row, z, zcoef) );
3606  SCIP_CALL( SCIPchgRowRhs(scip, *row, rhs + tmp + slope*xlb) );
3607 
3608  return SCIP_OKAY;
3609 }
3610 
3611 /** computes secant underestimator for sign(x+xoffset)abs(x+xoffset)^n + c*z <= rhs
3612  *
3613  * The generated cut is slope*xmult*x + c*z <= rhs + (-xlb-xoffset)^n + slope*xlb,
3614  * where slope = (sign(xub+xoffset)*abs(xub+xoffset)^n + (-xlb-xoffset)^n) / (xub - xlb).
3615  */
3616 static
3618  SCIP* scip, /**< SCIP data structure */
3619  SCIP_ROW** row, /**< buffer to store row */
3620  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3621  SCIP_Real xlb, /**< lower bound of x */
3622  SCIP_Real xub, /**< upper bound of x */
3623  SCIP_Real exponent, /**< exponent n in sign(x)abs(x)^n */
3624  SCIP_Real xoffset, /**< offset of x */
3625  DECL_MYPOW ((*mypow)), /**< function to use for computing power */
3626  SCIP_Real xmult, /**< multiplier for coefficient of x */
3627  SCIP_Real zcoef, /**< coefficient of z */
3628  SCIP_Real rhs, /**< right hand side */
3629  SCIP_VAR* x, /**< variable x */
3630  SCIP_VAR* z /**< variable z */
3631  )
3632 {
3633  SCIP_Real slope, tmp;
3634 
3635  assert(scip != NULL);
3636  assert(SCIPisLE(scip, xlb, xub));
3637  assert(!SCIPisPositive(scip, xlb + xoffset));
3638 
3639  /* ignore constraints with fixed x (should be removed soon) */
3640  if( SCIPisRelEQ(scip, xlb, xub) )
3641  return SCIP_OKAY;
3642 
3643  if( xlb > -xoffset )
3644  xlb = -xoffset;
3645 
3646  tmp = mypow(-xlb-xoffset, exponent);
3647  slope = SIGN(xub+xoffset) * mypow(ABS(xub+xoffset), exponent) + tmp;
3648  slope /= xub - xlb;
3649 
3650  if( SCIPisInfinity(scip, REALABS(slope)) )
3651  return SCIP_OKAY;
3652 
3653  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, conshdlr, "signpowcut", -SCIPinfinity(scip), SCIPinfinity(scip),
3654  SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0 /* local */, FALSE /* modifiable */, TRUE /* removable */ ) );
3655  SCIP_CALL( SCIPaddVarToRow(scip, *row, x, xmult*slope) );
3656  SCIP_CALL( SCIPaddVarToRow(scip, *row, z, zcoef) );
3657  SCIP_CALL( SCIPchgRowRhs(scip, *row, rhs + tmp + slope*xlb) );
3658 
3659  return SCIP_OKAY;
3660 }
3661 
3662 /** generates a cut
3663  * based on Liberti and Pantelides, Convex Envelopes of Monomials of Odd Degree, J. Global Optimization 25, 157-168, 2003, and previous publications
3664  */
3665 static
3667  SCIP* scip, /**< SCIP data structure */
3668  SCIP_CONS* cons, /**< constraint */
3669  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
3670  SCIP_ROW** row, /**< storage for cut */
3671  SCIP_Bool onlyinbounds /**< whether linearization is allowed only in variable bounds */
3672  )
3673 {
3674  SCIP_CONSHDLRDATA* conshdlrdata;
3675  SCIP_CONSDATA* consdata;
3676  SCIP_SIDETYPE violside;
3677  SCIP_Real c;
3678  SCIP_Real xlb;
3679  SCIP_Real xglb;
3680  SCIP_Real xub;
3681  SCIP_Real xval;
3682  SCIP_Real xoffset;
3683  SCIP_Real xmult;
3684  SCIP_Real zcoef;
3685  SCIP_Real rhs;
3686 
3687  assert(scip != NULL);
3688  assert(cons != NULL);
3689  assert(row != NULL);
3690 
3691  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
3692  assert(conshdlrdata != NULL);
3693 
3694  consdata = SCIPconsGetData(cons);
3695  assert(consdata != NULL);
3696 
3697  assert(SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) || SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)));
3698 
3699  violside = SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT;
3700  *row = NULL;
3701 
3702  SCIPdebugMsg(scip, "generate cut for constraint <%s> with violated side %d\n", SCIPconsGetName(cons), violside);
3703  SCIPdebugPrintCons(scip, cons, NULL);
3704  SCIPdebugMsg(scip, "xlb = %g xub = %g xval = %g\n", SCIPvarGetLbLocal(consdata->x), SCIPvarGetUbLocal(consdata->x), SCIPgetSolVal(scip, sol, consdata->x));
3705 
3706  if( violside == SCIP_SIDETYPE_RIGHT )
3707  {
3708  xglb = SCIPvarGetLbGlobal(consdata->x);
3709  xlb = SCIPvarGetLbLocal(consdata->x);
3710  xub = SCIPvarGetUbLocal(consdata->x);
3711  xval = SCIPgetSolVal(scip, sol, consdata->x);
3712  xoffset = consdata->xoffset;
3713  xmult = 1.0;
3714  zcoef = consdata->zcoef;
3715  rhs = consdata->rhs;
3716  }
3717  else
3718  {
3719  xglb = -SCIPvarGetUbGlobal(consdata->x);
3720  xlb = -SCIPvarGetUbLocal(consdata->x);
3721  xub = -SCIPvarGetLbLocal(consdata->x);
3722  xval = -SCIPgetSolVal(scip, sol, consdata->x);
3723  xoffset = -consdata->xoffset;
3724  xmult = -1.0;
3725  zcoef = -consdata->zcoef;
3726  rhs = -consdata->lhs;
3727  }
3728 
3729  if( SCIPisInfinity(scip, REALABS(xval)) )
3730  {
3731  SCIPdebugMsg(scip, "skip separation since x is at infinity\n");
3732  return SCIP_OKAY;
3733  }
3734 
3735  if( !SCIPisNegative(scip, xlb+xoffset) )
3736  {
3737  /* [xlb, xub] completely in positive orthant -> function is convex on whole domain */
3738  SCIP_Bool islocal;
3739 
3740  islocal = (!SCIPconsIsGlobal(cons) || SCIPisNegative(scip, xglb+xoffset)) && SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0;
3741  if( conshdlrdata->projectrefpoint && !onlyinbounds )
3742  {
3743  SCIP_CALL( generateLinearizationCutProject(scip, row, SCIPconsGetHdlr(cons), xval, SCIPgetSolVal(scip, sol, consdata->z), -xoffset, consdata->exponent,
3744  xoffset, xmult, zcoef, rhs, consdata->x, consdata->z, islocal) );
3745  }
3746  else if( !onlyinbounds )
3747  {
3748  SCIP_CALL( generateLinearizationCut(scip, row, SCIPconsGetHdlr(cons), xval, consdata->exponent, xoffset, xmult, zcoef, rhs,
3749  consdata->x, consdata->z, islocal) );
3750  }
3751  else
3752  {
3753  SCIP_CALL( generateLinearizationCut(scip, row, SCIPconsGetHdlr(cons), 2.0*xval > xlb + xub ? xub : xlb, consdata->exponent, xoffset, xmult, zcoef, rhs,
3754  consdata->x, consdata->z, islocal) );
3755  }
3756  }
3757  else if( !SCIPisPositive(scip, xub+xoffset) )
3758  {
3759  /* [xlb, xub] completely in negative orthant -> function is concave on whole domain */
3760  if( SCIPisInfinity(scip, -xlb) )
3761  return SCIP_OKAY;
3762  SCIP_CALL( generateSecantCut(scip, row, SCIPconsGetHdlr(cons), sol, xlb, xub, consdata->exponent, xoffset, consdata->power, xmult, zcoef, rhs, consdata->x, consdata->z) );
3763  }
3764  else if( (c = - consdata->root * (xlb+xoffset) - xoffset) > xub )
3765  {
3766  /* c is right of xub -> use secant */
3767  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) )
3768  return SCIP_OKAY;
3769  SCIP_CALL( generateSecantCut(scip, row, SCIPconsGetHdlr(cons), sol, xlb, xub, consdata->exponent, xoffset, consdata->power, xmult, zcoef, rhs, consdata->x, consdata->z) );
3770  }
3771  else if( xval >= c )
3772  {
3773  /* xval is right of c -> use linearization */
3774  if( conshdlrdata->projectrefpoint && !onlyinbounds )
3775  SCIP_CALL( generateLinearizationCutProject(scip, row, SCIPconsGetHdlr(cons), xval, SCIPgetSolVal(scip, sol, consdata->z), c, consdata->exponent,
3776  xoffset, xmult, zcoef, rhs, consdata->x, consdata->z, SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0) );
3777  else if( !onlyinbounds )
3778  SCIP_CALL( generateLinearizationCut(scip, row, SCIPconsGetHdlr(cons), xval, consdata->exponent, xoffset, xmult, zcoef, rhs,
3779  consdata->x, consdata->z, xval+xoffset < - consdata->root * (xglb+xoffset) && SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0) );
3780  else
3781  SCIP_CALL( generateLinearizationCut(scip, row, SCIPconsGetHdlr(cons), xub, consdata->exponent, xoffset, xmult, zcoef, rhs,
3782  consdata->x, consdata->z, xval+xoffset < - consdata->root * (xglb+xoffset) && SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) > 0) );
3783  }
3784  else
3785  {
3786  /* xval between xlb and c -> use secant */
3787  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, c) )
3788  return SCIP_OKAY;
3789  SCIP_CALL( generateSecantCut(scip, row, SCIPconsGetHdlr(cons), sol, xlb, c, consdata->exponent, xoffset, consdata->power, xmult, zcoef, rhs, consdata->x, consdata->z) );
3790  }
3791 
3792  /* check numerics */
3793  if( *row != NULL )
3794  {
3795  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) );
3796 
3797  /* check range of coefficients */
3798  SCIPdebugMsg(scip, " -> found cut rhs=%f, min=%f, max=%f range=%g\n",
3799  SCIProwGetRhs(*row),
3800  SCIPgetRowMinCoef(scip, *row), SCIPgetRowMaxCoef(scip, *row),
3801  SCIPgetRowMaxCoef(scip, *row)/SCIPgetRowMinCoef(scip, *row));
3802 
3803  if( SCIPisInfinity(scip, REALABS(SCIProwGetRhs(*row))) )
3804  {
3805  SCIPdebugMsg(scip, "skip cut for constraint <%s> because of very large right hand side: %g\n", SCIPconsGetName(cons), SCIProwGetRhs(*row));
3806  SCIP_CALL( SCIPreleaseRow(scip, row) );
3807  return SCIP_OKAY;
3808  }
3809 
3810  if( SCIPgetRowMaxCoef(scip, *row) / SCIPgetRowMinCoef(scip, *row) >= conshdlrdata->cutmaxrange )
3811  {
3812  SCIPdebugMsg(scip, "skip cut for constraint <%s> because of very large range: %g\n", SCIPconsGetName(cons), SCIPgetRowMaxCoef(scip, *row)/SCIPgetRowMinCoef(scip, *row));
3813  SCIP_CALL( SCIPreleaseRow(scip, row) );
3814  return SCIP_OKAY;
3815  }
3816  }
3817 
3818  return SCIP_OKAY;
3819 }
3820 
3821 /** tries to separate solution or LP solution by a linear cut
3822  * assumes that constraint violations have been computed
3823  */
3824 static
3826  SCIP* scip, /**< SCIP data structure */
3827  SCIP_CONSHDLR* conshdlr, /**< quadratic constraints handler */
3828  SCIP_CONS** conss, /**< constraints */
3829  int nconss, /**< number of constraints */
3830  int nusefulconss, /**< number of constraints that seem to be useful */
3831  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
3832  SCIP_Real minefficacy, /**< minimal efficacy of a cut if it should be added to the LP */
3833  SCIP_Bool inenforcement, /**< whether we are in constraint enforcement */
3834  SCIP_Bool onlyinbounds, /**< whether linearization is allowed only in variable bounds */
3835  SCIP_Bool* success, /**< result of separation: separated point (TRUE) or not (FALSE) */
3836  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
3837  SCIP_Real* bestefficacy /**< buffer to store best efficacy of a cut that was added to the LP, if found; or NULL if not of interest */
3838  )
3839 {
3840  SCIP_CONSHDLRDATA* conshdlrdata;
3841  SCIP_CONSDATA* consdata;
3842  SCIP_Real efficacy;
3843  SCIP_Real feasibility;
3844  SCIP_Real norm;
3845  SCIP_Bool convex;
3846  int c;
3847  SCIP_ROW* row;
3848 
3849  assert(scip != NULL);
3850  assert(conshdlr != NULL);
3851  assert(conss != NULL || nconss == 0);
3852  assert(success != NULL);
3853  assert(cutoff != NULL);
3854 
3855  *success = FALSE;
3856  *cutoff = FALSE;
3857 
3858  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3859  assert(conshdlrdata != NULL);
3860 
3861  if( bestefficacy != NULL )
3862  *bestefficacy = 0.0;
3863 
3864  for( c = 0; c < nconss && ! (*cutoff); ++c )
3865  {
3866  assert(conss[c] != NULL); /*lint !e613*/
3867 
3868  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
3869  assert(consdata != NULL);
3870 
3871  if( SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) || SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
3872  {
3873  /* try to generate a cut */
3874  SCIP_CALL( generateCut(scip, conss[c], sol, &row, onlyinbounds) ); /*lint !e613*/
3875  if( row == NULL ) /* failed to generate cut */
3876  continue;
3877 
3878  /* check if we separate in convex area */
3879  if( SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
3880  {
3881  convex = !SCIPisInfinity(scip, -SCIPvarGetLbLocal(consdata->x))
3882  && (!SCIPisNegative(scip, SCIPvarGetLbLocal(consdata->x)+consdata->xoffset)
3883  || SCIPgetSolVal(scip, sol, consdata->x)+consdata->xoffset >= -consdata->root*(SCIPvarGetLbLocal(consdata->x)+consdata->xoffset));
3884  }
3885  else
3886  {
3887  convex = !SCIPisInfinity(scip, SCIPvarGetUbLocal(consdata->x))
3888  && (!SCIPisPositive(scip, SCIPvarGetUbLocal(consdata->x)+consdata->xoffset)
3889  || SCIPgetSolVal(scip, sol, consdata->x)+consdata->xoffset <= -consdata->root*(SCIPvarGetUbLocal(consdata->x)+consdata->xoffset));
3890  }
3891 
3892  feasibility = SCIPgetRowSolFeasibility(scip, row, sol);
3893 
3894  switch( conshdlrdata->scaling )
3895  {
3896  case 'o' :
3897  efficacy = -feasibility;
3898  break;
3899 
3900  case 'g' :
3901  /* in difference to SCIPgetCutEfficacy, we scale by norm only if the norm is > 1.0 this avoid finding cuts
3902  * efficient which are only very slightly violated CPLEX does not seem to scale row coefficients up too
3903  * also we use infinity norm, since that seem to be the usual scaling strategy in LP solvers (equilibrium
3904  * scaling) */
3905  norm = SCIPgetRowMaxCoef(scip, row);
3906  efficacy = -feasibility / MAX(1.0, norm);
3907  break;
3908 
3909  case 's' :
3910  {
3911  SCIP_Real abslhs = REALABS(SCIProwGetLhs(row));
3912  SCIP_Real absrhs = REALABS(SCIProwGetRhs(row));
3913  SCIP_Real minval = MIN(abslhs, absrhs);
3914 
3915  efficacy = -feasibility / MAX(1.0, minval);
3916  break;
3917  }
3918 
3919  default:
3920  SCIPerrorMessage("Unknown scaling method '%c'.", conshdlrdata->scaling);
3921  SCIPABORT();
3922  return SCIP_INVALIDDATA; /*lint !e527*/
3923  }
3924 
3925  /* if cut is strong or it's weak but we are convex and desperate (speak, in enforcement), then add,
3926  * unless it corresponds to a bound change that is too weak (<eps) to be added
3927  */
3928  if( (efficacy > minefficacy || (inenforcement && convex && (SCIPgetRelaxFeastolFactor(scip) > 0.0 ? SCIPisPositive(scip, efficacy) : SCIPisFeasPositive(scip, efficacy)))) &&
3929  SCIPisCutApplicable(scip, row) )
3930  {
3931  SCIP_Bool infeasible;
3932 
3933  SCIP_CALL( SCIPaddCut(scip, sol, row, FALSE, &infeasible) );
3934  if ( infeasible )
3935  *cutoff = TRUE;
3936  else
3937  *success = TRUE;
3938  if( bestefficacy != NULL && efficacy > *bestefficacy )
3939  *bestefficacy = efficacy;
3940 
3941  /* notify indicator constraint handler about this cut */
3942  if( conshdlrdata->conshdlrindicator != NULL && !SCIProwIsLocal(row) )
3943  {
3944  SCIP_CALL( SCIPaddRowIndicator(scip, conshdlrdata->conshdlrindicator, row) );
3945  }
3946 
3947  /* mark row as not removable from LP for current node, if in enforcement */
3948  if( inenforcement && !conshdlrdata->enfocutsremovable )
3949  SCIPmarkRowNotRemovableLocal(scip, row);
3950  }
3951 
3952  SCIP_CALL( SCIPreleaseRow (scip, &row) );
3953  }
3954 
3955  /* enforce only useful constraints
3956  * others are only checked and enforced if we are still feasible or have not found a separating cut yet
3957  */
3958  if( c >= nusefulconss && *success )
3959  break;
3960  }
3961 
3962  return SCIP_OKAY;
3963 }
3964 
3965 /** adds linearizations cuts for convex constraints w.r.t. a given reference point to cutpool and sepastore
3966  * if separatedlpsol is not NULL, then a cut that separates the LP solution is added to the sepastore and is forced to enter the LP
3967  * if separatedlpsol is not NULL, but cut does not separate the LP solution, then it is added to the cutpool only
3968  * if separatedlpsol is NULL, then cut is added to cutpool only
3969  */
3970 static
3972  SCIP* scip, /**< SCIP data structure */
3973  SCIP_CONSHDLR* conshdlr, /**< quadratic constraints handler */
3974  SCIP_CONS** conss, /**< constraints */
3975  int nconss, /**< number of constraints */
3976  SCIP_SOL* ref, /**< reference point where to linearize, or NULL for LP solution */
3977  SCIP_Bool* separatedlpsol, /**< buffer to store whether a cut that separates the current LP solution was found and added to LP, or NULL if adding to cutpool only */
3978  SCIP_Real minefficacy /**< minimal efficacy of a cut when checking for separation of LP solution */
3979  )
3980 {
3981  SCIP_CONSDATA* consdata;
3982  SCIP_Bool addedtolp;
3983  SCIP_ROW* row;
3984  int c;
3985 
3986  assert(scip != NULL);
3987  assert(conshdlr != NULL);
3988  assert(conss != NULL || nconss == 0);
3989 
3990  if( separatedlpsol != NULL )
3991  *separatedlpsol = FALSE;
3992 
3993  for( c = 0; c < nconss; ++c )
3994  {
3995  assert(conss[c] != NULL); /*lint !e613*/
3996 
3997  if( SCIPconsIsLocal(conss[c]) ) /*lint !e613*/
3998  continue;
3999 
4000  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4001  assert(consdata != NULL);
4002 
4003  if( !SCIPisGT(scip, SCIPvarGetUbGlobal(consdata->x), -consdata->xoffset) && !SCIPisInfinity(scip, -consdata->lhs) )
4004  {
4005  /* constraint function is concave for x+offset <= 0.0, so can linearize w.r.t. lhs */
4006  consdata->lhsviol = 1.0;
4007  consdata->rhsviol = 0.0;
4008  SCIP_CALL( generateCut(scip, conss[c], ref, &row, FALSE) ); /*lint !e613*/
4009  }
4010  else if( !SCIPisLT(scip, SCIPvarGetLbGlobal(consdata->x), -consdata->xoffset) && !SCIPisInfinity(scip, -consdata->rhs) )
4011  {
4012  /* constraint function is convex for x+offset >= 0.0, so can linearize w.r.t. rhs */
4013  consdata->lhsviol = 0.0;
4014  consdata->rhsviol = 1.0;
4015  SCIP_CALL( generateCut(scip, conss[c], ref, &row, FALSE) ); /*lint !e613*/
4016  }
4017  else
4018  {
4019  /* sign not fixed or nothing to linearize */
4020  continue;
4021  }
4022 
4023  if( row == NULL )
4024  continue;
4025 
4026  addedtolp = FALSE;
4027 
4028  assert(!SCIProwIsLocal(row));
4029 
4030  /* if caller wants, then check if cut separates LP solution and add to sepastore if so */
4031  if( separatedlpsol != NULL )
4032  {
4033  SCIP_CONSHDLRDATA* conshdlrdata;
4034  SCIP_Real efficacy;
4035  SCIP_Real norm;
4036 
4037  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4038  assert(conshdlrdata != NULL);
4039 
4040  efficacy = -SCIPgetRowLPFeasibility(scip, row);
4041  switch( conshdlrdata->scaling )
4042  {
4043  case 'o' :
4044  break;
4045 
4046  case 'g' :
4047  /* in difference to SCIPgetCutEfficacy, we scale by norm only if the norm is > 1.0 this avoid finding cuts
4048  * efficient which are only very slightly violated CPLEX does not seem to scale row coefficients up too
4049  * also we use infinity norm, since that seem to be the usual scaling strategy in LP solvers (equilibrium
4050  * scaling) */
4051  norm = SCIPgetRowMaxCoef(scip, row);
4052  efficacy /= MAX(1.0, norm);
4053  break;
4054 
4055  case 's' :
4056  {
4057  SCIP_Real abslhs = REALABS(SCIProwGetLhs(row));
4058  SCIP_Real absrhs = REALABS(SCIProwGetRhs(row));
4059  SCIP_Real minval = MIN(abslhs, absrhs);
4060 
4061  efficacy /= MAX(1.0, minval);
4062  break;
4063  }
4064 
4065  default:
4066  SCIPerrorMessage("Unknown scaling method '%c'.", conshdlrdata->scaling);
4067  SCIPABORT();
4068  return SCIP_INVALIDDATA; /*lint !e527*/
4069  }
4070 
4071  if( efficacy >= minefficacy )
4072  {
4073  SCIP_Bool infeasible;
4074 
4075  *separatedlpsol = TRUE;
4076  addedtolp = TRUE;
4077  SCIP_CALL( SCIPaddCut(scip, NULL, row, TRUE, &infeasible) );
4078  assert( ! infeasible );
4079  }
4080  }
4081 
4082  if( !addedtolp )
4083  {
4084  SCIP_CALL( SCIPaddPoolCut(scip, row) );
4085  }
4086 
4087  SCIP_CALL( SCIPreleaseRow(scip, &row) );
4088  }
4089 
4090  return SCIP_OKAY;
4091 }
4092 
4093 /** processes the event that a new primal solution has been found */
4094 static
4095 SCIP_DECL_EVENTEXEC(processNewSolutionEvent)
4097  SCIP_CONSHDLRDATA* conshdlrdata;
4098  SCIP_CONSHDLR* conshdlr;
4099  SCIP_CONS** conss;
4100  int nconss;
4101  SCIP_SOL* sol;
4102 
4103  assert(scip != NULL);
4104  assert(event != NULL);
4105  assert(eventdata != NULL);
4106  assert(eventhdlr != NULL);
4107 
4108  assert((SCIPeventGetType(event) & SCIP_EVENTTYPE_SOLFOUND) != 0);
4109 
4110  conshdlr = (SCIP_CONSHDLR*)eventdata;
4111 
4112  nconss = SCIPconshdlrGetNConss(conshdlr);
4113 
4114  if( nconss == 0 )
4115  return SCIP_OKAY;
4116 
4117  sol = SCIPeventGetSol(event);
4118  assert(sol != NULL);
4119 
4120  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4121  assert(conshdlrdata != NULL);
4122 
4123  /* we are only interested in solution coming from some heuristic other than trysol, but not from the tree
4124  * the reason for ignoring trysol solutions is that they may come from an NLP solve in sepalp, where we already added linearizations,
4125  * or are from the tree, but postprocessed via proposeFeasibleSolution
4126  */
4127  if( SCIPsolGetHeur(sol) == NULL || SCIPsolGetHeur(sol) == conshdlrdata->trysolheur )
4128  return SCIP_OKAY;
4129 
4130  conss = SCIPconshdlrGetConss(conshdlr);
4131  assert(conss != NULL);
4132 
4133  SCIPdebugMsg(scip, "catched new sol event %"SCIP_EVENTTYPE_FORMAT" from heur <%s>; have %d conss\n", SCIPeventGetType(event), SCIPheurGetName(SCIPsolGetHeur(sol)), nconss);
4134 
4135  SCIP_CALL( addLinearizationCuts(scip, conshdlr, conss, nconss, sol, NULL, 0.0) );
4136 
4137  return SCIP_OKAY;
4138 }
4139 
4140 /** given a solution, try to make absolute power constraints feasible by shifting the linear variable z and pass this solution to the trysol heuristic */
4141 static
4143  SCIP* scip, /**< SCIP data structure */
4144  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4145  SCIP_CONS** conss, /**< constraints to process */
4146  int nconss, /**< number of constraints */
4147  SCIP_SOL* sol /**< solution to process */
4148  )
4149 {
4150  SCIP_CONSDATA* consdata;
4151  SCIP_SOL* newsol;
4152  SCIP_Real xtermval;
4153  SCIP_Real zval;
4154  SCIP_Real viol;
4155  SCIP_Bool solviolbounds;
4156  int c;
4157 
4158  assert(scip != NULL);
4159  assert(conshdlr != NULL);
4160  assert(conss != NULL || nconss == 0);
4161 
4162  /* don't propose new solutions if not in presolve or solving */
4164  return SCIP_OKAY;
4165 
4166  if( sol != NULL )
4167  {
4168  SCIP_CALL( SCIPcreateSolCopy(scip, &newsol, sol) );
4169  }
4170  else
4171  {
4172  SCIP_CALL( SCIPcreateLPSol(scip, &newsol, NULL) );
4173  }
4174  SCIP_CALL( SCIPunlinkSol(scip, newsol) );
4175 
4176  for( c = 0; c < nconss; ++c )
4177  {
4178  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
4179  assert(consdata != NULL);
4180  assert(consdata->z != NULL);
4181  assert(consdata->zcoef != 0.0);
4182 
4183  /* recompute violation w.r.t. current solution */
4184  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], newsol, &viol, &solviolbounds) ); /*lint !e613*/
4185  assert(!solviolbounds);
4186 
4187  /* do nothing if constraint is satisfied */
4188  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
4189  continue;
4190 
4191  /* if violation is at infinity, give up */
4192  if( SCIPisInfinity(scip, MAX(consdata->lhsviol, consdata->rhsviol)) )
4193  break;
4194 
4195  /* @todo could also adjust x while keeping z fixed */
4196 
4197  /* if variable is multiaggregated, then cannot set its solution value, so give up */
4198  if( SCIPvarGetStatus(consdata->z) == SCIP_VARSTATUS_MULTAGGR )
4199  break;
4200 
4201  /* compute value of x-term */
4202  xtermval = SCIPgetSolVal(scip, newsol, consdata->x);
4203  xtermval += consdata->xoffset;
4204  xtermval = SIGN(xtermval) * consdata->power(ABS(xtermval), consdata->exponent);
4205 
4206  /* if left hand side is violated, try to set z such that lhs is active */
4207  if( SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) )
4208  {
4209  assert(!SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip))); /* should only have one side violated (otherwise some variable is at infinity) */
4210 
4211  zval = (consdata->lhs - xtermval)/consdata->zcoef;
4212  /* bad luck: z would get value outside of its domain */
4213  if( SCIPisInfinity(scip, REALABS(zval)) || SCIPisFeasLT(scip, zval, SCIPvarGetLbGlobal(consdata->z)) || SCIPisFeasGT(scip, zval, SCIPvarGetUbGlobal(consdata->z)) )
4214  break;
4215  SCIP_CALL( SCIPsetSolVal(scip, newsol, consdata->z, zval) );
4216  }
4217 
4218  /* if right hand side is violated, try to set z such that rhs is active */
4219  if( SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
4220  {
4221  zval = (consdata->rhs - xtermval)/consdata->zcoef;
4222  /* bad luck: z would get value outside of its domain */
4223  if( SCIPisInfinity(scip, REALABS(zval)) || SCIPisFeasLT(scip, zval, SCIPvarGetLbGlobal(consdata->z)) || SCIPisFeasGT(scip, zval, SCIPvarGetUbGlobal(consdata->z)) )
4224  break;
4225  SCIP_CALL( SCIPsetSolVal(scip, newsol, consdata->z, zval) );
4226  }
4227  }
4228 
4229  /* if we have a solution that should satisfy all absolute power constraints and has a better objective than the current upper bound, then pass it to the trysol heuristic */
4230  if( c == nconss )
4231  {
4232  SCIP_CONSHDLRDATA* conshdlrdata;
4233 
4234  SCIPdebugMsg(scip, "pass solution with objective val %g to trysol heuristic\n", SCIPgetSolTransObj(scip, newsol));
4235 
4236  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4237  assert(conshdlrdata != NULL);
4238  assert(conshdlrdata->trysolheur != NULL);
4239 
4240  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, newsol) );
4241  }
4242 
4243  SCIP_CALL( SCIPfreeSol(scip, &newsol) );
4244 
4245  return SCIP_OKAY;
4246 }
4247 
4248 /** create a nonlinear row representation of the constraint and stores them in consdata */
4249 static
4251  SCIP* scip, /**< SCIP data structure */
4252  SCIP_CONS* cons /**< absolute power constraint */
4253  )
4254 {
4255  SCIP_CONSDATA* consdata;
4256  SCIP_EXPRTREE* exprtree;
4257  SCIP_QUADELEM quadelem;
4258  SCIP_VAR* linvars[2];
4259  SCIP_Real lincoefs[2];
4260  SCIP_VAR* quadvar;
4261  SCIP_Real constant;
4262  SCIP_Bool expisint;
4263  int sign;
4264  int nlinvars;
4265  int nquadvars;
4266  int nquadelems;
4267  int n;
4268 
4269  assert(scip != NULL);
4270  assert(cons != NULL);
4271 
4272  consdata = SCIPconsGetData(cons);
4273  assert(consdata != NULL);
4274 
4275  if( consdata->nlrow != NULL )
4276  {
4277  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
4278  }
4279 
4280  nlinvars = 0;
4281  nquadvars = 0;
4282  nquadelems = 0;
4283  exprtree = NULL;
4284  constant = 0.0;
4285 
4286  /* check if sign of x is fixed */
4287  if( !SCIPisNegative(scip, SCIPvarGetLbGlobal(consdata->x)+consdata->xoffset) )
4288  sign = 1;
4289  else if( !SCIPisPositive(scip, SCIPvarGetUbGlobal(consdata->x)+consdata->xoffset) )
4290  sign = -1;
4291  else
4292  sign = 0;
4293 
4294  /* check if exponent is integral */
4295  expisint = SCIPisIntegral(scip, consdata->exponent);
4296  n = (int)SCIPround(scip, consdata->exponent);
4297 
4298  /* create quadelem or expression tree for nonlinear part sign(x+offset)abs(x+offset)^n */
4299  if( sign != 0 || (expisint && (n % 2 == 1)) )
4300  {
4301  /* sign is fixes or exponent is odd integer */
4302  if( expisint && n == 2 )
4303  {
4304  /* sign of x is clear and exponent is 2.0 -> generate quadratic, linear, and constant term for +/- (x+offset)^n */
4305  assert(sign == -1 || sign == 1);
4306  nquadelems = 1;
4307  quadelem.idx1 = 0;
4308  quadelem.idx2 = 0;
4309  quadelem.coef = (SCIP_Real)sign;
4310  nquadvars = 1;
4311  quadvar = consdata->x;
4312 
4313  if( consdata->xoffset != 0.0 )
4314  {
4315  linvars[0] = consdata->x;
4316  lincoefs[0] = sign * 2.0 * consdata->xoffset;
4317  nlinvars = 1;
4318  constant = sign * consdata->xoffset * consdata->xoffset;
4319  }
4320  }
4321  else
4322  {
4323  /* exponent is odd or sign of x is clear, generate expression tree for +/- (+/-(x+offset))^exponent */
4324  SCIP_EXPR* expr;
4325  SCIP_EXPR* expr2;
4326 
4327  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* x */
4328  if( consdata->xoffset != 0.0 )
4329  {
4330  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->xoffset) );
4331  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, expr2) ); /* x + offset */
4332  }
4333  if( sign == -1 && !expisint )
4334  {
4335  /* if exponent is not integer and x is negative, then negate */
4336  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, -1.0) ); /* -1 */
4337  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MUL, expr, expr2) ); /* -(x+offset) */
4338  }
4339  /* use intpower for integer exponent and realpower for fractional exponent */
4340  if( expisint )
4341  {
4342  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_INTPOWER, expr, n) ); /* (x+offset)^n */
4343  }
4344  else
4345  {
4346  assert(sign == 1 || sign == -1);
4347  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_REALPOWER, expr, consdata->exponent) ); /* abs(x+offset)^exponent */
4348  }
4349  /* if exponent is even integer, then negate result; if it's an odd integer, then intpower already takes care of correct sign */
4350  if( sign == -1 && !(expisint && n % 2 == 1) )
4351  {
4352  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, -1.0) ); /* -1 */
4353  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MUL, expr, expr2) ); /* -abs(x+offset)^exponent */
4354  }
4355  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, expr, 1, 0, NULL) );
4356  }
4357  }
4358  else
4359  {
4360  /* exponent is not odd integer and sign of x is not fixed -> generate expression tree for signpower(x+offset, n) */
4361  SCIP_EXPR* expr;
4362  SCIP_EXPR* expr2;
4363 
4364  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* x */
4365  if( consdata->xoffset != 0.0 )
4366  {
4367  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_CONST, consdata->xoffset) );
4368  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, expr2) ); /* x + offset */
4369  }
4370  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_SIGNPOWER, expr, (SCIP_Real)consdata->exponent) ); /* signpower(x+offset, n) */
4371 
4372  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &exprtree, expr, 1, 0, NULL) );
4373  }
4374  assert(exprtree != NULL || nquadelems > 0);
4375 
4376  /* tell expression tree, which is its variable */
4377  if( exprtree != NULL )
4378  {
4379  SCIP_CALL( SCIPexprtreeSetVars(exprtree, 1, &consdata->x) );
4380  }
4381 
4382  assert(nlinvars < 2);
4383  linvars[nlinvars] = consdata->z;
4384  lincoefs[nlinvars] = consdata->zcoef;
4385  ++nlinvars;
4386 
4387  /* create nlrow */
4388  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), constant,
4389  nlinvars, linvars, lincoefs,
4390  nquadvars, &quadvar, nquadelems, &quadelem,
4391  exprtree, consdata->lhs, consdata->rhs,
4393  ) );
4394 
4395  if( exprtree != NULL )
4396  {
4397  SCIP_CALL( SCIPexprtreeFree(&exprtree) );
4398  }
4399 
4400  return SCIP_OKAY;
4401 }
4402 
4403 /** upgrades a quadratic constraint where the quadratic part is only a single square term and the quadratic variable sign is fixed to a signpower constraint */
4404 static
4405 SCIP_DECL_QUADCONSUPGD(quadconsUpgdAbspower)
4406 { /*lint --e{715}*/
4407  SCIP_QUADVARTERM quadvarterm;
4408  SCIP_VAR* x;
4409  SCIP_VAR* z;
4410  SCIP_Real xoffset;
4411  SCIP_Real zcoef;
4412  SCIP_Real signpowcoef;
4413  SCIP_Real lhs;
4414  SCIP_Real rhs;
4415 
4416  *nupgdconss = 0;
4417 
4418  /* need at least one linear variable */
4419  if( SCIPgetNLinearVarsQuadratic(scip, cons) == 0 )
4420  return SCIP_OKAY;
4421 
4422  /* consider only quadratic constraints with a single square term */
4423  if( SCIPgetNQuadVarTermsQuadratic(scip, cons) != 1 )
4424  return SCIP_OKAY;
4425  assert(SCIPgetNBilinTermsQuadratic(scip, cons) == 0);
4426 
4427  quadvarterm = SCIPgetQuadVarTermsQuadratic(scip, cons)[0];
4428  if( SCIPisZero(scip, quadvarterm.sqrcoef) )
4429  return SCIP_OKAY;
4430 
4431  x = quadvarterm.var;
4432  xoffset = quadvarterm.lincoef / (2.0 * quadvarterm.sqrcoef);
4433 
4434  /* check that x has fixed sign */
4435  if( SCIPisNegative(scip, SCIPvarGetLbGlobal(x) + xoffset) && SCIPisPositive(scip, SCIPvarGetUbGlobal(x) + xoffset) )
4436  return SCIP_OKAY;
4437 
4438  /* check whether upgdconss array has enough space to store 1 or 2 constraints */
4439  if( SCIPgetNLinearVarsQuadratic(scip, cons) > 1 )
4440  *nupgdconss = -2;
4441  else
4442  *nupgdconss = -1;
4443  if( -*nupgdconss > upgdconsssize )
4444  return SCIP_OKAY;
4445 
4446  *nupgdconss = 0;
4447 
4448  SCIPdebugMsg(scip, "upgrade quadratic constraint <%s> to absolute power, x = [%g,%g], offset = %g\n", SCIPconsGetName(cons), SCIPvarGetLbGlobal(x), SCIPvarGetUbGlobal(x), xoffset);
4449  SCIPdebugPrintCons(scip, cons, NULL);
4450 
4451  lhs = SCIPgetLhsQuadratic(scip, cons);
4452  rhs = SCIPgetRhsQuadratic(scip, cons);
4453 
4454  /* get z and its coefficient */
4455  if( SCIPgetNLinearVarsQuadratic(scip, cons) > 1 )
4456  {
4457  /* create auxiliary variable and constraint for linear part, since we can handle only at most one variable in cons_signpower */
4458  char name[SCIP_MAXSTRLEN];
4459  SCIP_VAR* auxvar;
4460 
4461  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_linpart", SCIPconsGetName(cons));
4462  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
4464  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4465 
4466  SCIP_CALL( SCIPcreateConsLinear(scip, &upgdconss[0], name, SCIPgetNLinearVarsQuadratic(scip, cons),
4468  SCIPisInfinity(scip, -lhs) ? -SCIPinfinity(scip) : 0.0,
4469  SCIPisInfinity(scip, rhs) ? SCIPinfinity(scip) : 0.0,
4473  SCIPconsIsStickingAtNode(cons)) );
4474  SCIP_CALL( SCIPaddCoefLinear(scip, upgdconss[*nupgdconss], auxvar, -1.0) );
4475 
4476  z = auxvar;
4477  zcoef = 1.0;
4478 
4479  ++*nupgdconss;
4480 
4481  /* compute and set value of auxvar in debug solution */
4482 #ifdef SCIP_DEBUG_SOLUTION
4483  if( SCIPdebugIsMainscip(scip) )
4484  {
4485  SCIP_Real debugval;
4486  SCIP_Real debugvarval;
4487  int i;
4488 
4489  debugval = 0.0;
4490  for( i = 0; i < SCIPgetNLinearVarsQuadratic(scip, cons); ++i )
4491  {
4492  SCIP_CALL( SCIPdebugGetSolVal(scip, SCIPgetLinearVarsQuadratic(scip, cons)[i], &debugvarval) );
4493  debugval += SCIPgetCoefsLinearVarsQuadratic(scip, cons)[i] * debugvarval;
4494  }
4495 
4496  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
4497  }
4498 #endif
4499 
4500  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4501  }
4502  else
4503  {
4504  assert(SCIPgetNLinearVarsQuadratic(scip, cons) == 1);
4505  z = SCIPgetLinearVarsQuadratic(scip, cons)[0];
4506  zcoef = SCIPgetCoefsLinearVarsQuadratic(scip, cons)[0];
4507  }
4508 
4509  /* we now have lhs <= sqrcoef * (x + offset)^2 - sqrcoef * offset^2 + zcoef * z <= rhs */
4510 
4511  /* move sqrcoef * offset^2 into lhs and rhs */
4512  if( !SCIPisInfinity(scip, -lhs) )
4513  lhs += quadvarterm.sqrcoef * xoffset * xoffset;
4514  if( !SCIPisInfinity(scip, rhs) )
4515  rhs += quadvarterm.sqrcoef * xoffset * xoffset;
4516 
4517  /* divide by sqrcoef if x+offset > 0 and by -sqrcoef if < 0 */
4518  signpowcoef = quadvarterm.sqrcoef;
4519  if( SCIPisNegative(scip, SCIPvarGetLbGlobal(x) + xoffset) )
4520  signpowcoef = -signpowcoef;
4521  if( signpowcoef > 0.0 )
4522  {
4523  if( !SCIPisInfinity(scip, -lhs) )
4524  lhs /= signpowcoef;
4525  if( !SCIPisInfinity(scip, rhs) )
4526  rhs /= signpowcoef;
4527  }
4528  else
4529  {
4530  SCIP_Real newrhs;
4531 
4532  if( !SCIPisInfinity(scip, -lhs) )
4533  newrhs = lhs / signpowcoef;
4534  else
4535  newrhs = SCIPinfinity(scip);
4536  if( !SCIPisInfinity(scip, rhs) )
4537  lhs = rhs / signpowcoef;
4538  else
4539  lhs = -SCIPinfinity(scip);
4540  rhs = newrhs;
4541  }
4542  zcoef /= signpowcoef;
4543 
4544  /* create the absolute power constraint */
4545  SCIP_CALL( SCIPcreateConsAbspower(scip, &upgdconss[*nupgdconss], SCIPconsGetName(cons), x, z, 2.0,
4546  xoffset, zcoef, lhs, rhs,
4550  SCIPconsIsStickingAtNode(cons)) );
4551  SCIPdebugPrintCons(scip, upgdconss[*nupgdconss], NULL);
4552  ++*nupgdconss;
4553 
4554  return SCIP_OKAY;
4555 }
4556 
4557 /** tries to upgrade a nonlinear constraint into a absolute power constraint */
4558 static
4559 SCIP_DECL_NONLINCONSUPGD(nonlinconsUpgdAbspower)
4561  SCIP_EXPRGRAPH* exprgraph;
4562  SCIP_EXPRGRAPHNODE* node;
4563  SCIP_EXPRGRAPHNODE* child;
4564  SCIP_Real exponent;
4565  SCIP_VAR* x;
4566  SCIP_VAR* z;
4567  SCIP_Real signpowcoef;
4568  SCIP_Real zcoef;
4569  SCIP_Real xoffset;
4570  SCIP_Real constant;
4571  SCIP_Real lhs;
4572  SCIP_Real rhs;
4573 
4574  assert(nupgdconss != NULL);
4575  assert(upgdconss != NULL);
4576 
4577  *nupgdconss = 0;
4578 
4579  /* absolute power needs at least one linear variable (constraint is trivial, otherwise) */
4580  if( SCIPgetNLinearVarsNonlinear(scip, cons) == 0 )
4581  return SCIP_OKAY;
4582 
4583  node = SCIPgetExprgraphNodeNonlinear(scip, cons);
4584 
4585  /* no interest in linear constraints */
4586  if( node == NULL )
4587  return SCIP_OKAY;
4588 
4589  /* need exactly one argument */
4590  if( SCIPexprgraphGetNodeNChildren(node) != 1 )
4591  return SCIP_OKAY;
4592 
4593  constant = 0.0;
4594  signpowcoef = 1.0; /* coefficient of sign(x)abs(x)^n term, to be reformulated away... */
4595 
4596  child = SCIPexprgraphGetNodeChildren(node)[0];
4597 
4598  /* check if node expression fits to absolute power constraint */
4599  switch( SCIPexprgraphGetNodeOperator(node) )
4600  {
4601  case SCIP_EXPR_REALPOWER:
4602  /* realpower with exponent > 1.0 can always be signpower, since it assumes that argument is >= 0.0 */
4603  exponent = SCIPexprgraphGetNodeRealPowerExponent(node);
4604  if( exponent <= 1.0 )
4605  return SCIP_OKAY;
4606 
4607  assert(SCIPexprgraphGetNodeBounds(child).inf >= 0.0);
4608  break;
4609 
4610  case SCIP_EXPR_INTPOWER:
4611  {
4612  /* check if exponent > 1.0 and either odd or even with child having fixed sign */
4613  SCIP_INTERVAL childbounds;
4614 
4616  if( exponent <= 1.0 )
4617  return SCIP_OKAY;
4618 
4619  childbounds = SCIPexprgraphGetNodeBounds(child);
4620  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4621  return SCIP_OKAY;
4622 
4623  /* use x^exponent = -sign(x) |x|^exponent if exponent is even and x always negative */
4624  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 )
4625  signpowcoef = -1.0;
4626 
4627  break;
4628  }
4629 
4630  case SCIP_EXPR_SQUARE:
4631  {
4632  /* check if child has fixed sign */
4633  SCIP_INTERVAL childbounds;
4634 
4635  childbounds = SCIPexprgraphGetNodeBounds(child);
4636  if( childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4637  return SCIP_OKAY;
4638 
4639  /* use x^2 = -sign(x) |x|^2 if x is always negative */
4640  if( childbounds.inf < 0.0 )
4641  signpowcoef = -1.0;
4642 
4643  exponent = 2.0;
4644  break;
4645  }
4646 
4647  case SCIP_EXPR_SIGNPOWER:
4648  /* check if exponent > 1.0 */
4649  exponent = SCIPexprgraphGetNodeSignPowerExponent(node);
4650  if( exponent <= 1.0 )
4651  return SCIP_OKAY;
4652  break;
4653 
4654  case SCIP_EXPR_POLYNOMIAL:
4655  {
4656  SCIP_EXPRDATA_MONOMIAL* monomial;
4657  SCIP_INTERVAL childbounds;
4658 
4659  /* check if only one univariate monomial with exponent > 1.0 */
4660 
4661  /* if sum of univariate monomials, then this should have been taken care of by exprgraphnodeReformSignpower */
4663  return SCIP_OKAY;
4664  assert(SCIPexprgraphGetNodePolynomialNMonomials(node) == 1); /* assume simplified, i.e., no constant polynomial */
4665 
4666  monomial = SCIPexprgraphGetNodePolynomialMonomials(node)[0];
4667  assert(SCIPexprGetMonomialNFactors(monomial) == 1); /* since we have only one children and assume simplified */
4668 
4669  exponent = SCIPexprGetMonomialExponents(monomial)[0];
4670  if( exponent <= 1.0 )
4671  return SCIP_OKAY;
4672 
4673  /* if exponent is even integer and child has mixed sign, then cannot do
4674  * if exponent is even integer and child is always negative, then can do via multiplication by -1.0 */
4675  childbounds = SCIPexprgraphGetNodeBounds(child);
4676  if( SCIPisIntegral(scip, exponent) && ((int)SCIPround(scip, exponent) % 2 == 0) && childbounds.inf < 0.0 )
4677  {
4678  if( childbounds.sup > 0.0 )
4679  return SCIP_OKAY;
4680  signpowcoef = -1.0;
4681  }
4682 
4683  constant = SCIPexprgraphGetNodePolynomialConstant(node);
4684  signpowcoef *= SCIPexprGetMonomialCoef(monomial);
4685 
4686  break;
4687  }
4688 
4689  default:
4690  return SCIP_OKAY;
4691  } /*lint !e788*/
4692  assert(SCIPexprgraphGetNodeNChildren(node) == 1);
4693 
4694  /* check magnitue of coefficient of z in signpower constraint */
4695  zcoef = 1.0;
4696  if( SCIPgetNLinearVarsNonlinear(scip, cons) == 1 )
4697  zcoef = SCIPgetLinearCoefsNonlinear(scip, cons)[0];
4698  zcoef /= signpowcoef;
4700  {
4701  zcoef /= pow(REALABS(SCIPexprgraphGetNodeLinearCoefs(child)[0]), exponent);
4702  }
4703  if( SCIPisZero(scip, zcoef) )
4704  {
4705  SCIPdebugMsg(scip, "skip upgrade to signpower since |zcoef| = %g would be zero\n", zcoef);
4706  return SCIP_OKAY;
4707  }
4708 
4709  /* count how many constraints we need to add (use negative numbers, for convenience):
4710  * one constraint for absolute power,
4711  * plus one if we need to replace the linear part by single variable,
4712  * plus one if we need to replace the argument of absolute power by a single variable
4713  */
4714  *nupgdconss = -1;
4715 
4717  {
4718  /* if node has known curvature and we would add auxiliary var for child, then don't upgrade
4719  * it's not really necessary, but may introduce more numerical troubles
4720  * @todo maybe still do if child is linear?
4721  */
4723  {
4724  *nupgdconss = 0;
4725  return SCIP_OKAY;
4726  }
4727 
4728  --*nupgdconss;
4729  }
4730 
4731  if( SCIPgetNLinearVarsNonlinear(scip, cons) > 1 )
4732  --*nupgdconss;
4733 
4734  /* request larger upgdconss array */
4735  if( upgdconsssize < -*nupgdconss )
4736  return SCIP_OKAY;
4737 
4738  SCIPdebugMsg(scip, "upgrading constraint <%s>\n", SCIPconsGetName(cons));
4739 
4740  /* start counting at zero again */
4741  *nupgdconss = 0;
4742 
4743  exprgraph = SCIPgetExprgraphNonlinear(scip, SCIPconsGetHdlr(cons));
4744 
4745  lhs = SCIPgetLhsNonlinear(scip, cons);
4746  rhs = SCIPgetRhsNonlinear(scip, cons);
4747 
4748  /* get x and it's offset */
4750  {
4751  x = (SCIP_VAR*)SCIPexprgraphGetNodeVar(exprgraph, child);
4752  xoffset = 0.0;
4753  }
4755  {
4756  SCIP_Real xcoef;
4757 
4759  xcoef = SCIPexprgraphGetNodeLinearCoefs(child)[0];
4760  assert(!SCIPisZero(scip, xcoef));
4761 
4762  signpowcoef *= (xcoef < 0.0 ? -1.0 : 1.0) * pow(REALABS(xcoef), exponent);
4763  xoffset = SCIPexprgraphGetNodeLinearConstant(child) / xcoef;
4764  }
4765  else
4766  {
4767  /* reformulate by adding auxiliary variable and constraint for child */
4768  char name[SCIP_MAXSTRLEN];
4769  SCIP_INTERVAL bounds;
4770  SCIP_VAR* auxvar;
4771  SCIP_Real minusone;
4772 
4773  bounds = SCIPexprgraphGetNodeBounds(child);
4774  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_powerarg", SCIPconsGetName(cons));
4775 
4776  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s for node %p(%d,%d)\n", name, (void*)child, SCIPexprgraphGetNodeDepth(child), SCIPexprgraphGetNodePosition(child));
4777 
4778  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, SCIPintervalGetInf(bounds), SCIPintervalGetSup(bounds), 0.0,
4780  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4781 
4782  /* create new constraint child == auxvar
4783  * since signpower is monotonic, we need only child <= auxvar or child >= auxvar, if not both sides are finite, and depending on signpowcoef
4784  * i.e., we need child - auxvar <= 0.0 if rhs is finite and signpowcoef > 0.0 or lhs is finite and signpowcoef < 0.0
4785  * and we need 0.0 <= child - auxvar if lhs is finite and signpowcoef > 0.0 or rhs is finite and signpowcoef < 0.0
4786  */
4787  minusone = -1.0;
4788  assert(upgdconsssize > *nupgdconss);
4789  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &upgdconss[*nupgdconss], name, 1, &auxvar, &minusone, child,
4790  ((signpowcoef > 0.0 && !SCIPisInfinity(scip, -lhs)) || (signpowcoef < 0.0 && !SCIPisInfinity(scip, rhs))) ? 0.0 : -SCIPinfinity(scip),
4791  ((signpowcoef > 0.0 && !SCIPisInfinity(scip, rhs)) || (signpowcoef < 0.0 && !SCIPisInfinity(scip, -lhs))) ? 0.0 : SCIPinfinity(scip),
4792  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
4793  ++*nupgdconss;
4794 
4795  /* use auxvar to setup absolute power constraint */
4796  x = auxvar;
4797  xoffset = 0.0;
4798 
4799  /* compute and set value of auxvar in debug solution, if debugging is enabled */
4800  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SCIPexprgraphGetNodeVal(child)) ); /*lint !e506 !e774*/
4801 
4802  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4803  }
4804 
4805  /* get z and its coefficient */
4806  if( SCIPgetNLinearVarsNonlinear(scip, cons) > 1 )
4807  {
4808  /* create auxiliary variable and constraint for linear part, since we can handle only at most one variable in cons_signpower */
4809  char name[SCIP_MAXSTRLEN];
4810  SCIP_VAR* auxvar;
4811 
4812  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_linpart", SCIPconsGetName(cons));
4813  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
4815  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4816 
4817  assert(upgdconsssize > *nupgdconss);
4818  SCIP_CALL( SCIPcreateConsLinear(scip, &upgdconss[*nupgdconss], name, SCIPgetNLinearVarsNonlinear(scip, cons),
4820  SCIPisInfinity(scip, -lhs) ? -SCIPinfinity(scip) : 0.0,
4821  SCIPisInfinity(scip, rhs) ? SCIPinfinity(scip) : 0.0,
4825  SCIPconsIsStickingAtNode(cons)) );
4826  SCIP_CALL( SCIPaddCoefLinear(scip, upgdconss[*nupgdconss], auxvar, -1.0) );
4827 
4828  z = auxvar;
4829  zcoef = 1.0;
4830 
4831  ++*nupgdconss;
4832 
4833  /* compute and set value of auxvar in debug solution */
4834 #ifdef SCIP_DEBUG_SOLUTION
4835  if( SCIPdebugIsMainscip(scip) )
4836  {
4837  SCIP_Real debugval;
4838  SCIP_Real debugvarval;
4839  int i;
4840 
4841  debugval = 0.0;
4842  for( i = 0; i < SCIPgetNLinearVarsNonlinear(scip, cons); ++i )
4843  {
4844  SCIP_CALL( SCIPdebugGetSolVal(scip, SCIPgetLinearVarsNonlinear(scip, cons)[i], &debugvarval) );
4845  debugval += SCIPgetLinearCoefsNonlinear(scip, cons)[i] * debugvarval;
4846  }
4847 
4848  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
4849  }
4850 #endif
4851 
4852  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4853  }
4854  else
4855  {
4856  assert(SCIPgetNLinearVarsNonlinear(scip, cons) == 1);
4857  z = SCIPgetLinearVarsNonlinear(scip, cons)[0];
4858  zcoef = SCIPgetLinearCoefsNonlinear(scip, cons)[0];
4859  }
4860 
4861  if( constant != 0.0 )
4862  {
4863  if( !SCIPisInfinity(scip, -lhs) )
4864  lhs -= constant;
4865  if( !SCIPisInfinity(scip, rhs) )
4866  rhs -= constant;
4867  }
4868 
4869  /* divide absolute power constraint by signpowcoef */
4870  if( signpowcoef != 1.0 )
4871  {
4872  zcoef /= signpowcoef;
4873  if( signpowcoef < 0.0 )
4874  {
4875  SCIP_Real newrhs;
4876 
4877  newrhs = SCIPisInfinity(scip, -lhs) ? SCIPinfinity(scip) : lhs/signpowcoef;
4878  lhs = SCIPisInfinity(scip, rhs) ? -SCIPinfinity(scip) : rhs/signpowcoef;
4879  rhs = newrhs;
4880  }
4881  else
4882  {
4883  if( !SCIPisInfinity(scip, -lhs) )
4884  lhs /= signpowcoef;
4885  if( !SCIPisInfinity(scip, rhs) )
4886  rhs /= signpowcoef;
4887  }
4888  }
4889 
4890  /* finally setup a absolute power constraint */
4891 
4892  assert(*nupgdconss < upgdconsssize);
4893  SCIP_CALL( SCIPcreateConsAbspower(scip, &upgdconss[*nupgdconss], SCIPconsGetName(cons),
4894  x, z, exponent, xoffset, zcoef, lhs, rhs,
4898  SCIPconsIsStickingAtNode(cons)) );
4899  ++*nupgdconss;
4900 
4901  return SCIP_OKAY;
4902 }
4903 
4904 /** tries to reformulate a expression graph node via introducing a absolute power constraint
4905  * if node fits to absolute power and has indefinte curvature and has no nonlinear parents and has siblings, then replace by auxvar and absolute power constraint
4906  * if it still has nonlinear parents, then we wait to see if reformulation code move node into auxiliary constraint,
4907  * so we do not add unnessary auxiliary variables for something like an x^2 in an exp(x^2)
4908  * if it has no siblings, then we let the upgrading for nonlinear constraints take care of it,
4909  * since it may be able to upgrade the constraint as a whole and can take the constraint sides into account too (may need only <=/>= auxcons)
4910  */
4911 static
4912 SCIP_DECL_EXPRGRAPHNODEREFORM(exprgraphnodeReformAbspower)
4914  SCIP_EXPRGRAPHNODE* child;
4915  char name[SCIP_MAXSTRLEN];
4916  SCIP_CONS* cons;
4917  SCIP_Real exponent;
4918  SCIP_VAR* x;
4919  SCIP_VAR* z;
4920  SCIP_Real signpowcoef;
4921  SCIP_Real xoffset;
4922  SCIP_Real constant;
4923 
4924  assert(scip != NULL);
4925  assert(exprgraph != NULL);
4926  assert(node != NULL);
4927  assert(naddcons != NULL);
4928  assert(reformnode != NULL);
4929 
4930  *reformnode = NULL;
4931 
4933  return SCIP_OKAY;
4934 
4935  constant = 0.0;
4936  signpowcoef = 1.0; /* coefficient of sign(x)abs(x)^n term, to be move in from of z... */
4937 
4938  /* check if node expression fits to absolute power constraint */
4939  switch( SCIPexprgraphGetNodeOperator(node) )
4940  {
4941  case SCIP_EXPR_REALPOWER:
4942  /* realpower with exponent > 1.0 can always be absolute power, since it assumes that argument is >= 0.0
4943  * @todo we should also ensure that argument is >= 0.0
4944  */
4945  exponent = SCIPexprgraphGetNodeRealPowerExponent(node);
4946  if( exponent <= 1.0 )
4947  return SCIP_OKAY;
4948 
4949  assert(SCIPexprgraphGetNodeBounds(SCIPexprgraphGetNodeChildren(node)[0]).inf >= 0.0);
4950  break;
4951 
4952  case SCIP_EXPR_INTPOWER:
4953  {
4954  /* check if exponent > 1.0 and either odd or even with child having fixed sign */
4955  SCIP_INTERVAL childbounds;
4956 
4958  if( exponent <= 1.0 )
4959  return SCIP_OKAY;
4960 
4962  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4963  return SCIP_OKAY;
4964 
4965  /* use x^exponent = -sign(x) |x|^exponent if exponent is even and x always negative */
4966  if( (int)exponent % 2 == 0 && childbounds.inf < 0.0 )
4967  signpowcoef = -1.0;
4968 
4969  break;
4970  }
4971 
4972  case SCIP_EXPR_SQUARE:
4973  {
4974  /* check if child has fixed sign */
4975  SCIP_INTERVAL childbounds;
4976 
4978  if( childbounds.inf < 0.0 && childbounds.sup > 0.0 )
4979  return SCIP_OKAY;
4980 
4981  /* use x^2 = -sign(x) |x|^2 if x is always negative */
4982  if( childbounds.inf < 0.0 )
4983  signpowcoef = -1.0;
4984 
4985  exponent = 2.0;
4986  break;
4987  }
4988 
4989  case SCIP_EXPR_SIGNPOWER:
4990  /* check if exponent > 1.0 */
4991  exponent = SCIPexprgraphGetNodeSignPowerExponent(node);
4992  if( exponent <= 1.0 )
4993  return SCIP_OKAY;
4994  break;
4995 
4996  case SCIP_EXPR_POLYNOMIAL:
4997  {
4998  SCIP_EXPRDATA_MONOMIAL* monomial;
4999  SCIP_INTERVAL childbounds;
5000 
5001  /* check if only one univariate monomial with exponent > 1.0 */
5002  if( SCIPexprgraphGetNodeNChildren(node) > 1 )
5003  return SCIP_OKAY;
5004  assert(SCIPexprgraphGetNodeNChildren(node) == 1);
5005 
5007  return SCIP_OKAY;
5008  assert(SCIPexprgraphGetNodePolynomialNMonomials(node) == 1); /* assume simplified, i.e., no constant polynomial */
5009 
5010  monomial = SCIPexprgraphGetNodePolynomialMonomials(node)[0];
5011  assert(SCIPexprGetMonomialNFactors(monomial) == 1); /* since we have only one children and assume simplified */
5012 
5013  exponent = SCIPexprGetMonomialExponents(monomial)[0];
5014  if( exponent <= 1.0 )
5015  return SCIP_OKAY;
5016 
5017  /* if exponent is even integer and child has mixed sign, then cannot do
5018  * if exponent is even integer and child is always negative, then can do via multiplication by -1.0 */
5020  if( SCIPisIntegral(scip, exponent) && ((int)SCIPround(scip, exponent) % 2 == 0) && childbounds.inf < 0.0 )
5021  {
5022  if( childbounds.sup > 0.0 )
5023  return SCIP_OKAY;
5024  signpowcoef = -1.0;
5025  }
5026 
5027  constant = SCIPexprgraphGetNodePolynomialConstant(node);
5028  signpowcoef *= SCIPexprGetMonomialCoef(monomial);
5029 
5030  break;
5031  }
5032 
5033  default:
5034  return SCIP_OKAY;
5035  } /*lint !e788*/
5036  assert(SCIPexprgraphGetNodeNChildren(node) == 1);
5037 
5039  return SCIP_OKAY;
5040  if( !SCIPexprgraphHasNodeSibling(node) )
5041  return SCIP_OKAY;
5042 
5043  SCIPdebugMsg(scip, "reformulate node %p via signpower\n", (void*)node);
5044 
5045  /* get x and its offset */
5046  child = SCIPexprgraphGetNodeChildren(node)[0];
5048  {
5049  x = (SCIP_VAR*)SCIPexprgraphGetNodeVar(exprgraph, child);
5050  xoffset = 0.0;
5051  }
5053  {
5054  SCIP_Real xcoef;
5055 
5057  xcoef = SCIPexprgraphGetNodeLinearCoefs(child)[0];
5058  assert(!SCIPisZero(scip, xcoef));
5059 
5060  signpowcoef *= (xcoef < 0.0 ? -1.0 : 1.0) * pow(REALABS(xcoef), exponent);
5061  xoffset = SCIPexprgraphGetNodeLinearConstant(child) / xcoef;
5062  }
5063  else
5064  {
5065  /* reformulate by adding auxiliary variable and constraint for child */
5066  SCIP_INTERVAL bounds;
5067  SCIP_VAR* auxvar;
5068  SCIP_Real minusone;
5069 
5070  bounds = SCIPexprgraphGetNodeBounds(child);
5071  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%dsp", *naddcons);
5072 
5073  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s for node %p(%d,%d)\n", name, (void*)child, SCIPexprgraphGetNodeDepth(child), SCIPexprgraphGetNodePosition(child));
5074 
5075  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, SCIPintervalGetInf(bounds), SCIPintervalGetSup(bounds), 0.0, SCIP_VARTYPE_CONTINUOUS,
5076  TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
5077  SCIP_CALL( SCIPaddVar(scip, auxvar) );
5078 
5079  /* create new constraint child == auxvar */
5080  minusone = -1.0;
5081  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &cons, name, 1, &auxvar, &minusone, child, 0.0, 0.0,
5082  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5083  SCIP_CALL( SCIPaddCons(scip, cons) );
5084  ++*naddcons;
5085 
5086  /* use auxvar to setup signpower constraint */
5087  x = auxvar;
5088  xoffset = 0.0;
5089 
5090  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SCIPexprgraphGetNodeVal(child)) ); /*lint !e506 !e774*/
5091 
5092  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
5093  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
5094  }
5095 
5096  /* create auxiliary variable z and add to expression graph */
5097  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%dsp", *naddcons);
5098  SCIP_CALL( SCIPcreateVar(scip, &z, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS,
5099  TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
5100  SCIP_CALL( SCIPaddVar(scip, z) );
5101  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&z, reformnode) );
5102 
5103  /* setup a absolute power constraint */
5104  if( REALABS(signpowcoef) * SCIPfeastol(scip) < 1.0 )
5105  {
5106  /* if signpowcoef is not huge (<10^6), then put it into absolute power constraint */
5107  SCIP_CALL( SCIPcreateConsAbspower(scip, &cons, name,
5108  x, z, exponent, xoffset, -1.0/signpowcoef, -constant/signpowcoef, -constant/signpowcoef,
5109  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5110  SCIP_CALL( SCIPaddCons(scip, cons) );
5111  SCIPdebugPrintCons(scip, cons, NULL);
5112  ++*naddcons;
5113 
5114  /* compute value of z and reformnode and set in debug solution and expression graph, resp. */
5115 #ifdef SCIP_DEBUG_SOLUTION
5116  if( SCIPdebugIsMainscip(scip) )
5117  {
5118  SCIP_Real xval;
5119  SCIP_Real zval;
5120 
5121  SCIP_CALL( SCIPdebugGetSolVal(scip, x, &xval) );
5122  zval = signpowcoef * SIGN(xval + xoffset) * pow(REALABS(xval + xoffset), exponent) + constant;
5123 
5124  SCIP_CALL( SCIPdebugAddSolVal(scip, z, zval) );
5125  SCIPexprgraphSetVarNodeValue(*reformnode, zval);
5126  }
5127 #endif
5128  }
5129  else
5130  {
5131  /* if signpowcoef is huge, then avoid very small coefficient of z
5132  * instead create additional node on top of current reformnode */
5133  SCIP_EXPRGRAPHNODE* linnode;
5134 
5135  SCIP_CALL( SCIPcreateConsAbspower(scip, &cons, name,
5136  x, z, exponent, xoffset, -1.0, 0.0, 0.0,
5137  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
5138  SCIP_CALL( SCIPaddCons(scip, cons) );
5139  SCIPdebugPrintCons(scip, cons, NULL);
5140  ++*naddcons;
5141 
5142  /* compute value of z and reformnode and set in debug solution and expression graph, resp. */
5143 #ifdef SCIP_DEBUG_SOLUTION
5144  if( SCIPdebugIsMainscip(scip) )
5145  {
5146  SCIP_Real xval;
5147  SCIP_Real zval;
5148 
5149  SCIP_CALL( SCIPdebugGetSolVal(scip, x, &xval) );
5150  zval = SIGN(xval + xoffset) * pow(REALABS(xval + xoffset), exponent);
5151 
5152  SCIP_CALL( SCIPdebugAddSolVal(scip, z, zval) );
5153  SCIPexprgraphSetVarNodeValue(*reformnode, zval);
5154  }
5155 #endif
5156 
5157  SCIP_CALL( SCIPexprgraphCreateNodeLinear(SCIPblkmem(scip), &linnode, 1, &signpowcoef, constant) );
5158  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, linnode, -1, 1, reformnode) );
5159 
5160  *reformnode = linnode;
5161  }
5162 
5163  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
5164  SCIP_CALL( SCIPreleaseVar(scip, &z) );
5165 
5166  return SCIP_OKAY;
5167 }
5168 
5169 /** helper function to enforce constraints */
5170 static
5172  SCIP* scip, /**< SCIP data structure */
5173  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
5174  SCIP_CONS** conss, /**< constraints to process */
5175  int nconss, /**< number of constraints */
5176  int nusefulconss, /**< number of useful (non-obsolete) constraints to process */
5177  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
5178  SCIP_Bool solinfeasible, /**< was the solution already declared infeasible by a constraint handler? */
5179  SCIP_RESULT* result /**< pointer to store the result of the enforcing call */
5180  )
5181 {
5182  SCIP_CONSHDLRDATA* conshdlrdata;
5183  SCIP_CONS* maxviolcons;
5184  SCIP_CONSDATA* consdata;
5185  SCIP_Bool success;
5186  SCIP_Bool cutoff;
5187  SCIP_Bool solviolbounds;
5188  SCIP_Real minefficacy;
5189  SCIP_Real sepaefficacy;
5190  SCIP_Real leastpossibleefficacy;
5191  SCIP_Real maxviol;
5192  int nnotify;
5193  int c;
5194 
5195  assert(scip != NULL);
5196  assert(conshdlr != NULL);
5197  assert(conss != NULL || nconss == 0);
5198  assert(result != NULL);
5199 
5200  conshdlrdata = SCIPconshdlrGetData(conshdlr);
5201  assert(conshdlrdata != NULL);
5202 
5203  SCIP_CALL( computeViolations(scip, conshdlr, conss, nconss, sol, &solviolbounds, &maxviolcons) );
5204 
5205  if( maxviolcons == NULL )
5206  {
5207  *result = SCIP_FEASIBLE;
5208  return SCIP_OKAY;
5209  }
5210 
5211  *result = SCIP_INFEASIBLE;
5212 
5213  if( solviolbounds )
5214  {
5215  /* if LP solution violates variable bounds, then this should be because a row was added that
5216  * introduced this variable newly to the LP, in which case it gets value 0.0; the row should
5217  * have been added to resolve an infeasibility, so solinfeasible should be TRUE
5218  * see also issue #627
5219  */
5220  assert(solinfeasible);
5221  return SCIP_OKAY;
5222  }
5223 
5224  /* if we are above the 100'th enforcement round for this node, something is strange
5225  * (maybe the LP does not think that the cuts we add are violated, or we do ECP on a high-dimensional convex function)
5226  * in this case, check if some limit is hit or SCIP should stop for some other reason and terminate enforcement by creating a dummy node
5227  * (in optimized more, returning SCIP_INFEASIBLE in *result would be sufficient, but in debug mode this would give an assert in scip.c)
5228  * the reason to wait for 100 rounds is to avoid calls to SCIPisStopped in normal runs, which may be expensive
5229  * we only increment nenforounds until 101 to avoid an overflow
5230  */
5231  if( conshdlrdata->lastenfonode == SCIPgetCurrentNode(scip) )
5232  {
5233  if( conshdlrdata->nenforounds > 100 )
5234  {
5235  if( SCIPisStopped(scip) )
5236  {
5237  SCIP_NODE* child;
5238 
5239  SCIP_CALL( SCIPcreateChild(scip, &child, 1.0, SCIPnodeGetEstimate(SCIPgetCurrentNode(scip))) );
5240  *result = SCIP_BRANCHED;
5241 
5242  return SCIP_OKAY;
5243  }
5244  }
5245  else
5246  ++conshdlrdata->nenforounds;
5247  }
5248  else
5249  {
5250  conshdlrdata->lastenfonode = SCIPgetCurrentNode(scip);
5251  conshdlrdata->nenforounds = 0;
5252  }
5253 
5254  /* run domain propagation for violated constraints */
5255  for( c = 0; c < nconss; ++c )
5256  {
5257  int nchgbds;
5258  int naddconss;
5259 
5260  assert(conss[c] != NULL); /*lint !e613*/
5261 
5262  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
5263  assert(consdata != NULL);
5264 
5265  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
5266  continue;
5267 
5268  nchgbds = 0;
5269  naddconss = 0;
5270  SCIP_CALL( propagateCons(scip, conshdlr, conss[c], TRUE, &cutoff, &nchgbds, &naddconss) ); /*lint !e613*/
5271  if( cutoff )
5272  {
5273  *result = SCIP_CUTOFF;
5274  return SCIP_OKAY;
5275  }
5276  if( nchgbds )
5277  *result = SCIP_REDUCEDDOM;
5278  if( naddconss )
5279  *result = SCIP_CONSADDED;
5280  }
5281  if( *result == SCIP_REDUCEDDOM || *result == SCIP_CONSADDED )
5282  return SCIP_OKAY;
5283 
5284  consdata = SCIPconsGetData(maxviolcons);
5285  assert(consdata != NULL);
5286  maxviol = consdata->lhsviol + consdata->rhsviol;
5287  assert(SCIPisGT(scip, maxviol, SCIPfeastol(scip)));
5288 
5289  /* we would like a cut that is efficient enough that it is not redundant in the LP (>feastol)
5290  * however, if the maximal violation is very small, also the best cut efficacy cannot be large
5291  * thus, in the latter case, we are also happy if the efficacy is at least, say, 75% of the maximal violation
5292  * but in any case we need an efficacy that is at least feastol
5293  */
5294  minefficacy = MIN(0.75*maxviol, conshdlrdata->mincutefficacyenfofac * SCIPfeastol(scip)); /*lint !e666*/
5295  minefficacy = MAX(minefficacy, SCIPfeastol(scip)); /*lint !e666*/
5296  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, sol, minefficacy, TRUE, FALSE, &success,
5297  &cutoff, &sepaefficacy) );
5298  if( cutoff )
5299  {
5300  SCIPdebugMsg(scip, "separation detected cutoff.\n");
5301  *result = SCIP_CUTOFF;
5302  return SCIP_OKAY;
5303  }
5304  if( success )
5305  {
5306  SCIPdebugMsg(scip, "separation succeeded (bestefficacy = %g, minefficacy = %g)\n", sepaefficacy, minefficacy);
5307  *result = SCIP_SEPARATED;
5308  return SCIP_OKAY;
5309  }
5310  SCIPdebugMsg(scip, "separation failed (bestefficacy = %g < %g = minefficacy ); max viol: %g\n", sepaefficacy, minefficacy,
5311  maxviol);
5312 
5313  /* we are not feasible, the whole node is not infeasible, and we cannot find a reasonable cut
5314  * -> collect variables for branching
5315  */
5316  SCIP_CALL( registerBranchingCandidates(scip, conshdlr, conss, nconss, sol, &nnotify) );
5317 
5318  /* if sepastore can decrease LP feasibility tolerance, we can add cuts with efficacy in [eps, feastol] */
5319  leastpossibleefficacy = SCIPgetRelaxFeastolFactor(scip) > 0.0 ? SCIPepsilon(scip) : SCIPfeastol(scip);
5320  if( nnotify == 0 && !solinfeasible && minefficacy > leastpossibleefficacy )
5321  {
5322  /* fallback 1: we also have no branching candidates, so try to find a weak cut */
5323  SCIP_CALL( separatePoint(scip, conshdlr, conss, nconss, nusefulconss, sol, leastpossibleefficacy, TRUE, FALSE,
5324  &success, &cutoff, &sepaefficacy) );
5325  if( cutoff )
5326  {
5327  SCIPdebugMsg(scip, "separation detected cutoff.\n");
5328  *result = SCIP_CUTOFF;
5329  return SCIP_OKAY;
5330  }
5331  if( success )
5332  {
5333  *result = SCIP_SEPARATED;
5334  return SCIP_OKAY;
5335  }
5336  }
5337 
5338  if( nnotify == 0 && !solinfeasible )
5339  {
5340  /* fallback 2: separation probably failed because of numerical difficulties with a convex constraint;
5341  if noone declared solution infeasible yet and we had not even found a weak cut, try to resolve by branching */
5342  SCIP_VAR* brvar = NULL;
5343  SCIP_CALL( registerLargeRelaxValueVariableForBranching(scip, conss, nconss, sol, &brvar) );
5344  if( brvar == NULL )
5345  {
5346  SCIPwarningMessage(scip, "Could not find any branching variable candidate. Cutting off node. Max viol = %g.\n",
5347  SCIPconsGetData(maxviolcons)->lhsviol+SCIPconsGetData(maxviolcons)->rhsviol);
5348  *result = SCIP_CUTOFF;
5349  return SCIP_OKAY;
5350  }
5351  else
5352  {
5353  SCIPdebugMsg(scip, "Could not find any usual branching variable candidate. Proposed variable %s with LP value %g for branching. Max. viol. cons. <%s>: %g+%g\n",
5354  SCIPvarGetName(brvar), SCIPgetSolVal(scip, sol, brvar), SCIPconsGetName(maxviolcons),
5355  SCIPconsGetData(maxviolcons)->lhsviol, SCIPconsGetData(maxviolcons)->rhsviol);
5356  nnotify = 1;
5357  }
5358  }
5359 
5360  assert(*result == SCIP_INFEASIBLE && (solinfeasible || nnotify > 0));
5361  return SCIP_OKAY;
5362 }
5363 
5364 /*
5365  * Callback methods of constraint handler
5366  */
5367 
5368 /** copy method for constraint handler plugins (called when SCIP copies plugins) */
5369 static
5370 SCIP_DECL_CONSHDLRCOPY(conshdlrCopyAbspower)
5371 { /*lint --e{715}*/
5372  assert(scip != NULL);
5373  assert(conshdlr != NULL);
5374  assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
5375 
5376  /* call inclusion method of constraint handler */
5378 
5379  *valid = TRUE;
5380 
5381  return SCIP_OKAY;
5382 }
5383 
5384 /** destructor of constraint handler to free constraint handler data (called when SCIP is exiting) */
5385 static
5386 SCIP_DECL_CONSFREE(consFreeAbspower)
5387 { /*lint --e{715}*/
5388  SCIP_CONSHDLRDATA* conshdlrdata;
5389 
5390  assert(scip != NULL);
5391  assert(conshdlr != NULL);
5392 
5393  conshdlrdata = SCIPconshdlrGetData(conshdlr);
5394  assert(conshdlrdata != NULL);
539