Scippy

SCIP

Solving Constraint Integer Programs

cons_bivariate.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-2018 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_bivariate.c
17  * @brief constraint handler for bivariate nonlinear constraints \f$\textrm{lhs} \leq f(x,y) + c z \leq \textrm{rhs}\f$
18  * @author Martin Ballerstein
19  * @author Dennis Michaels
20  * @author Stefan Vigerske
21  */
22 
23 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
24 
25 #include "blockmemshell/memory.h"
26 #include "nlpi/exprinterpret.h"
27 #include "nlpi/pub_expr.h"
29 #include "scip/cons_bivariate.h"
30 #include "scip/cons_nonlinear.h"
31 #include "scip/cons_quadratic.h"
32 #include "scip/debug.h"
33 #include "scip/heur_subnlp.h"
34 #include "scip/heur_trysol.h"
35 #include "scip/intervalarith.h"
36 #include "scip/pub_cons.h"
37 #include "scip/pub_event.h"
38 #include "scip/pub_heur.h"
39 #include "scip/pub_lp.h"
40 #include "scip/pub_message.h"
41 #include "scip/pub_misc.h"
42 #include "scip/pub_nlp.h"
43 #include "scip/pub_sol.h"
44 #include "scip/pub_tree.h"
45 #include "scip/pub_var.h"
46 #include "scip/scip_branch.h"
47 #include "scip/scip_cons.h"
48 #include "scip/scip_copy.h"
49 #include "scip/scip_cut.h"
50 #include "scip/scip_event.h"
51 #include "scip/scip_expr.h"
52 #include "scip/scip_general.h"
53 #include "scip/scip_heur.h"
54 #include "scip/scip_lp.h"
55 #include "scip/scip_mem.h"
56 #include "scip/scip_message.h"
57 #include "scip/scip_nlp.h"
58 #include "scip/scip_numerics.h"
59 #include "scip/scip_param.h"
60 #include "scip/scip_prob.h"
61 #include "scip/scip_probing.h"
62 #include "scip/scip_sepa.h"
63 #include "scip/scip_sol.h"
64 #include "scip/scip_solvingstats.h"
65 #include "scip/scip_tree.h"
66 #include "scip/scip_var.h"
67 
68 /* constraint handler properties */
69 #define CONSHDLR_NAME "bivariate"
70 #define CONSHDLR_DESC "constraint handler for constraints of the form lhs <= f(x,y) + c*z <= rhs where f(x,y) is a bivariate function"
71 #define CONSHDLR_SEPAPRIORITY 5 /**< priority of the constraint handler for separation */
72 #define CONSHDLR_ENFOPRIORITY -55 /**< priority of the constraint handler for constraint enforcing */
73 #define CONSHDLR_CHECKPRIORITY -3600000 /**< priority of the constraint handler for checking feasibility */
74 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
75 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
76 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
77  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
78 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
79 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
80 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
81 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
82 
83 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_FAST
84 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP
85 
86 #define INTERVALINFTY 1E+43 /**< value for infinity in interval operations */
87 #define NEWTONMAXITER 1000 /**< maximal number of iterations in newton method */
88 #define INITLPMAXVARVAL 1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */
89 
90 #define QUADCONSUPGD_PRIORITY 5000 /**< priority of the constraint handler for upgrading of quadratic constraints */
91 #define NONLINCONSUPGD_PRIORITY 10000 /**< priority of the constraint handler for upgrading of nonlinear constraints */
92 
93 /* activate the following define to get output on number of bivariate constraints for each convexity-type during INITSOL */
94 /* #define TYPESTATISTICS */
95 
96 /*
97  * Data structures
98  */
99 
100 /** data structure to cache data used for separation of convex-concave constraints */
101 struct SepaData_ConvexConcave
102 {
103  SCIP_Bool linearinx; /**< whether the function is linear in x */
104  SCIP_Bool lineariny; /**< whether the function is linear in y */
105  SCIP_EXPRTREE* f_yfixed; /**< expression tree for f(x,yfixed) */
106  SCIP_EXPRTREE* f_neg_swapped; /**< expression tree for -f(y,x) */
107  SCIP_EXPRTREE* f_neg_swapped_yfixed;/**< expression tree for -f(y,xfixed) */
108  SCIP_EXPRTREE* vred; /**< expression tree for vred to underestimate f(x,y) */
109  SCIP_EXPRTREE* vred_neg_swapped; /**< expression tree for vred to underestimate -f(y,x) */
110 };
111 /** data structure to cache data used for separation of convex-concave constraints */
112 typedef struct SepaData_ConvexConcave SEPADATA_CONVEXCONCAVE;
114 /** constraint data for bivariate constraints */
115 struct SCIP_ConsData
116 {
117  SCIP_EXPRTREE* f; /**< expression tree of bivariate function f(x,y) */
118  SCIP_BIVAR_CONVEXITY convextype; /**< kind of convexity of f(x,y) */
119  SCIP_VAR* z; /**< linear variable */
120  SCIP_Real zcoef; /**< coefficient of linear variable */
121  SCIP_Real lhs; /**< left hand side */
122  SCIP_Real rhs; /**< right hand side */
123 
124  SCIP_Real activity; /**< activity of bivariate function w.r.t. current solution */
125  SCIP_Real lhsviol; /**< violation of left hand side in current solution */
126  SCIP_Real rhsviol; /**< violation of left hand side in current solution */
127 
128  unsigned int mayincreasez:1; /**< whether z can be increased without harming other constraints */
129  unsigned int maydecreasez:1; /**< whether z can be decreased without harming other constraints */
130  int eventfilterpos; /**< position of z var events in SCIP event filter */
131 
132  SCIP_EXPRGRAPHNODE* exprgraphnode; /**< node in expression graph corresponding to bivariate function */
133 
134  SEPADATA_CONVEXCONCAVE sepaconvexconcave; /**< separation data for convex-concave constraints */
135 };
136 
137 /** constraint handler data */
138 struct SCIP_ConshdlrData
139 {
140  SCIP_EXPRINT* exprinterpreter; /**< expression interpreter (computer gradients and hessians) */
141 
142  SCIP_Real cutmaxrange; /**< maximal range (maximal coef / minimal coef) of a cut in order to be added to LP */
143  SCIP_Bool linfeasshift; /**< whether to make solutions in check feasible if possible */
144  int maxproprounds; /**< limit on number of propagation rounds for a single constraint within one round of SCIP propagation */
145  int ninitlprefpoints; /**< number of reference points in each direction where to compute linear support for envelope in LP initialization */
146  SCIP_Bool enfocutsremovable; /**< are cuts added during enforcement removable from the LP in the same node? */
147 
148  SCIP_EVENTHDLR* linvareventhdlr; /**< handler for linear variable bound change events */
149  SCIP_EVENTHDLR* nonlinvareventhdlr; /**< handler for nonlinear variable bound change events */
150  SCIP_HEUR* subnlpheur; /**< a pointer to the subNLP heuristic */
151  SCIP_HEUR* trysolheur; /**< a pointer to the TRYSOL heuristic, if available */
152  int newsoleventfilterpos;/**< filter position of new solution event handler, if catched */
153 
154  SCIP_EXPRGRAPH* exprgraph; /**< expression graph */
155  SCIP_Bool isremovedfixings; /**< whether variable fixations have been removed from the expression graph */
156  SCIP_Bool ispropagated; /**< whether the bounds on the variables in the expression graph have been propagated */
157  SCIP* scip; /**< SCIP data structure, needed in expression graph callbacks */
158 
159  SCIP_NODE* lastenfonode; /**< the node for which enforcement was called the last time (and some constraint was violated) */
160  int nenforounds; /**< counter on number of enforcement rounds for the current node */
161 };
162 
163 
164 /*
165  * Local methods
166  */
167 
168 /** translate from one value of infinity to another
169  *
170  * if val is >= infty1, then give infty2, else give val
171  */
172 #define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))
174 /** processes bound tightening event */
175 static
176 SCIP_DECL_EVENTEXEC(processLinearVarEvent)
177 {
178  SCIP_CONS* cons;
179 
180  assert(scip != NULL);
181  assert(event != NULL);
182  assert(eventdata != NULL);
183  assert(eventhdlr != NULL);
185 
186  cons = (SCIP_CONS*) eventdata;
187  assert(cons != NULL);
188 
190 
191  return SCIP_OKAY;
192 }
193 
194 /** catches variable bound change events on the linear variable in a bivariate constraint */
195 static
197  SCIP* scip, /**< SCIP data structure */
198  SCIP_CONS* cons /**< constraint for which to catch bound change events */
199  )
200 {
201  SCIP_CONSHDLRDATA* conshdlrdata;
202  SCIP_CONSDATA* consdata;
203  SCIP_EVENTTYPE eventtype;
204 
205  assert(scip != NULL);
206  assert(cons != NULL);
207  assert(SCIPconsIsEnabled(cons));
208  assert(SCIPconsIsTransformed(cons));
209 
210  assert(SCIPconsGetHdlr(cons) != NULL);
211  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
212  assert(conshdlrdata != NULL);
213  assert(conshdlrdata->linvareventhdlr != NULL);
214 
215  consdata = SCIPconsGetData(cons);
216  assert(consdata != NULL);
217 
218  if( consdata->z == NULL )
219  return SCIP_OKAY;
220  assert(consdata->eventfilterpos == -1);
221 
222  eventtype = SCIP_EVENTTYPE_DISABLED;
223  if( !SCIPisInfinity(scip, consdata->rhs) )
224  {
225  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest */
226  if( consdata->zcoef > 0.0 )
227  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
228  else
229  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
230  }
231  if( !SCIPisInfinity(scip, -consdata->lhs) )
232  {
233  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest */
234  if( consdata->zcoef > 0.0 )
235  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
236  else
237  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
238  }
239 
240  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->z, eventtype, conshdlrdata->linvareventhdlr, (SCIP_EVENTDATA*)cons, &consdata->eventfilterpos) );
241 
242  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
243 
244  return SCIP_OKAY;
245 }
246 
247 /** drops variable bound change events on the linear variable in a bivariate constraint */
248 static
250  SCIP* scip, /**< SCIP data structure */
251  SCIP_CONS* cons /**< constraint for which to catch bound change events */
252  )
253 {
254  SCIP_CONSHDLRDATA* conshdlrdata;
255  SCIP_CONSDATA* consdata;
256  SCIP_EVENTTYPE eventtype;
257 
258  assert(scip != NULL);
259  assert(cons != NULL);
260  assert(SCIPconsIsTransformed(cons));
261 
262  assert(SCIPconsGetHdlr(cons) != NULL);
263  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
264  assert(conshdlrdata != NULL);
265  assert(conshdlrdata->linvareventhdlr != NULL);
266 
267  consdata = SCIPconsGetData(cons);
268  assert(consdata != NULL);
269 
270  if( consdata->z == NULL )
271  return SCIP_OKAY;
272  assert(consdata->eventfilterpos >= 0);
273 
274  eventtype = SCIP_EVENTTYPE_DISABLED;
275  if( !SCIPisInfinity(scip, consdata->rhs) )
276  {
277  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest */
278  if( consdata->zcoef > 0.0 )
279  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
280  else
281  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
282  }
283  if( !SCIPisInfinity(scip, -consdata->lhs) )
284  {
285  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest */
286  if( consdata->zcoef > 0.0 )
287  eventtype |= SCIP_EVENTTYPE_UBTIGHTENED;
288  else
289  eventtype |= SCIP_EVENTTYPE_LBTIGHTENED;
290  }
291 
292  SCIP_CALL( SCIPdropVarEvent(scip, consdata->z, eventtype, conshdlrdata->linvareventhdlr, (SCIP_EVENTDATA*)cons, consdata->eventfilterpos) );
293  consdata->eventfilterpos = -1;
294 
295  return SCIP_OKAY;
296 }
297 
298 
299 /** processes bound change events for variables in expression graph */
300 static
301 SCIP_DECL_EVENTEXEC(processNonlinearVarEvent)
302 {
303  SCIP_CONSHDLRDATA* conshdlrdata;
304  SCIP_EVENTTYPE eventtype;
305 
306  assert(scip != NULL);
307  assert(event != NULL);
308  assert(eventdata != NULL);
309  assert(eventhdlr != NULL);
310 
311  conshdlrdata = (SCIP_CONSHDLRDATA*)SCIPeventhdlrGetData(eventhdlr);
312  assert(conshdlrdata != NULL);
313  assert(conshdlrdata->exprgraph != NULL);
314 
315  eventtype = SCIPeventGetType(event);
316  assert( eventtype & (SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED) );
317 
318  if( eventtype & SCIP_EVENTTYPE_BOUNDCHANGED )
319  {
320  SCIPdebugMsg(scip, "changed %s bound on expression graph variable <%s> from %g to %g\n",
321  (eventtype & SCIP_EVENTTYPE_LBCHANGED) ? "lower" : "upper",
323 
324  if( eventtype & SCIP_EVENTTYPE_BOUNDTIGHTENED )
325  conshdlrdata->ispropagated = FALSE;
326 
327  /* update variable bound in expression graph
328  * @todo should we add epsilon to variable range?
329  */
330  if( eventtype & SCIP_EVENTTYPE_LBCHANGED )
331  SCIPexprgraphSetVarNodeLb(conshdlrdata->exprgraph, (SCIP_EXPRGRAPHNODE*)eventdata,
332  -infty2infty(SCIPinfinity(scip), INTERVALINFTY, -SCIPeventGetNewbound(event))); /*lint !e666*/
333  else
334  SCIPexprgraphSetVarNodeUb(conshdlrdata->exprgraph, (SCIP_EXPRGRAPHNODE*)eventdata,
335  +infty2infty(SCIPinfinity(scip), INTERVALINFTY, SCIPeventGetNewbound(event))); /*lint !e666*/
336  }
337  else
338  {
339  assert(eventtype & SCIP_EVENTTYPE_VARFIXED);
340  conshdlrdata->isremovedfixings = FALSE;
341  }
342 
343  return SCIP_OKAY;
344 }
345 
346 /** callback method for variable addition in expression graph */
347 static
348 SCIP_DECL_EXPRGRAPHVARADDED( exprgraphVarAdded )
349 {
350  SCIP_CONSHDLRDATA* conshdlrdata;
351  SCIP_INTERVAL varbounds;
352  SCIP_VAR* var_;
353 
354  assert(exprgraph != NULL);
355  assert(var != NULL);
356  assert(varnode != NULL);
357 
358  var_ = (SCIP_VAR*)var;
359 
360  conshdlrdata = (SCIP_CONSHDLRDATA*)userdata;
361  assert(conshdlrdata != NULL);
362  assert(conshdlrdata->exprgraph == exprgraph);
363 
364  /* catch variable bound change events */
365  SCIP_CALL( SCIPcatchVarEvent(conshdlrdata->scip, (SCIP_VAR*)var, SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED, conshdlrdata->nonlinvareventhdlr, (SCIP_EVENTDATA*)varnode, NULL) );
366  SCIPdebugMessage("catch boundchange events on new expression graph variable <%s>\n", SCIPvarGetName(var_));
367 
368  /* set current bounds in expression graph */
369  SCIPintervalSetBounds(&varbounds,
370  -infty2infty(SCIPinfinity(conshdlrdata->scip), INTERVALINFTY, -MIN(SCIPvarGetLbLocal(var_), SCIPvarGetUbLocal(var_))), /*lint !e666*/
371  +infty2infty(SCIPinfinity(conshdlrdata->scip), INTERVALINFTY, MAX(SCIPvarGetLbLocal(var_), SCIPvarGetUbLocal(var_))) /*lint !e666*/
372  );
373  SCIPexprgraphSetVarNodeBounds(exprgraph, varnode, varbounds);
374 
375  SCIP_CALL( SCIPaddVarLocksType(conshdlrdata->scip, var_, SCIP_LOCKTYPE_MODEL, 1, 1) );
376  SCIPdebugMessage("increased up- and downlocks of variable <%s>\n", SCIPvarGetName(var_));
377 
378  conshdlrdata->isremovedfixings &= SCIPvarIsActive(var_);
379  conshdlrdata->ispropagated = FALSE;
380 
381  return SCIP_OKAY;
382 }
383 
384 /** callback method for variable removal in expression graph */
385 static
386 SCIP_DECL_EXPRGRAPHVARREMOVE( exprgraphVarRemove )
387 {
388  SCIP_CONSHDLRDATA* conshdlrdata;
389  SCIP_VAR* var_;
390 
391  assert(exprgraph != NULL);
392  assert(var != NULL);
393  assert(varnode != NULL);
394 
395  var_ = (SCIP_VAR*)var;
396 
397  conshdlrdata = (SCIP_CONSHDLRDATA*)userdata;
398  assert(conshdlrdata != NULL);
399  assert(conshdlrdata->exprgraph == exprgraph);
400 
401  SCIP_CALL( SCIPdropVarEvent(conshdlrdata->scip, var_, SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED, conshdlrdata->nonlinvareventhdlr, (SCIP_EVENTDATA*)varnode, -1) );
402  SCIPdebugMessage("drop boundchange events on expression graph variable <%s>\n", SCIPvarGetName(var_));
403 
404  SCIP_CALL( SCIPaddVarLocksType(conshdlrdata->scip, var_, SCIP_LOCKTYPE_MODEL, -1, -1) );
405  SCIPdebugMessage("decreased up- and downlocks of variable <%s>\n", SCIPvarGetName(var_));
406 
407  return SCIP_OKAY;
408 }
409 
410 /** locks linear variable in a constraint */
411 static
413  SCIP* scip, /**< SCIP data structure */
414  SCIP_CONS* cons, /**< constraint where to lock a variable */
415  SCIP_VAR* var, /**< variable to lock */
416  SCIP_Real coef /**< coefficient of variable in constraint */
417  )
418 {
419  SCIP_CONSDATA* consdata;
420 
421  assert(scip != NULL);
422  assert(cons != NULL);
423  assert(var != NULL);
424  assert(coef != 0.0);
425 
426  consdata = SCIPconsGetData(cons);
427  assert(consdata != NULL);
428 
429  if( coef > 0.0 )
430  {
431  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
432  }
433  else
434  {
435  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
436  }
437 
438  return SCIP_OKAY;
439 }
440 
441 /** unlocks linear variable in a constraint */
442 static
444  SCIP* scip, /**< SCIP data structure */
445  SCIP_CONS* cons, /**< constraint where to unlock a variable */
446  SCIP_VAR* var, /**< variable to unlock */
447  SCIP_Real coef /**< coefficient of variable in constraint */
448  )
449 {
450  SCIP_CONSDATA* consdata;
451 
452  assert(scip != NULL);
453  assert(cons != NULL);
454  assert(var != NULL);
455  assert(coef != 0.0);
456 
457  consdata = SCIPconsGetData(cons);
458  assert(consdata != NULL);
459 
460  if( coef > 0.0 )
461  {
462  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
463  }
464  else
465  {
466  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
467  }
468 
469  return SCIP_OKAY;
470 }
471 
472 /** resolves variable fixations and aggregations in a constraint */
473 static
475  SCIP* scip, /**< SCIP data structure */
476  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
477  SCIP_CONS* cons, /**< constraint where to remove fixed variables */
478  SCIP_Bool* ischanged, /**< buffer to store whether something was changed in the constraint */
479  SCIP_Bool* isupgraded /**< buffer to store whether the constraint has been upgraded (and deleted) */
480  )
481 {
482 #ifndef NDEBUG
483  SCIP_CONSHDLRDATA* conshdlrdata;
484 #endif
485  SCIP_CONSDATA* consdata;
486  SCIP_EXPR* substexpr[2];
487  SCIP_VAR* var;
488  SCIP_VAR* vars[2];
489  SCIP_Real coef;
490  SCIP_Real constant;
491  int i;
492 
493  assert(conshdlr != NULL);
494  assert(scip != NULL);
495  assert(cons != NULL);
496  assert(ischanged != NULL);
497  assert(isupgraded != NULL);
498 
499 #ifndef NDEBUG
500  conshdlrdata = SCIPconshdlrGetData(conshdlr);
501  assert(conshdlrdata != NULL);
502 #endif
503 
504  consdata = SCIPconsGetData(cons);
505  assert(consdata != NULL);
506  assert(consdata->f != NULL);
507 
508  *ischanged = FALSE;
509  *isupgraded = FALSE;
510 
511  if( consdata->z != NULL && !SCIPvarIsActive(consdata->z) && SCIPvarGetStatus(consdata->z) != SCIP_VARSTATUS_MULTAGGR )
512  {
513  /* replace z by active or multaggr. variable */
514 
515  /* drop events on z, unlock and release variable */
516  SCIP_CALL( dropLinearVarEvents(scip, cons) );
517  SCIP_CALL( unlockLinearVariable(scip, cons, consdata->z, consdata->zcoef) );
518 
519  /* replace by new variable, or NULL */
520  constant = 0.0;
521  SCIP_CALL( SCIPgetProbvarSum(scip, &consdata->z, &consdata->zcoef, &constant) );
522  if( consdata->zcoef == 0.0 )
523  consdata->z = NULL;
524  if( constant != 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
525  consdata->lhs -= constant;
526  if( constant != 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
527  consdata->rhs -= constant;
528 
529  if( consdata->z != NULL )
530  {
531  /* catch events on new z, lock and capture variable, mark as not to multaggr */
532  SCIP_CALL( catchLinearVarEvents(scip, cons) );
533  SCIP_CALL( lockLinearVariable(scip, cons, consdata->z, consdata->zcoef) );
534  if( SCIPvarIsActive(consdata->z) )
535  {
536  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, consdata->z) );
537  }
538  }
539 
540  *ischanged = TRUE;
541  }
542 
543  assert(SCIPexprtreeGetNVars(consdata->f) == 2);
544  vars[0] = SCIPexprtreeGetVars(consdata->f)[0];
545  vars[1] = SCIPexprtreeGetVars(consdata->f)[1];
546 
547  if( vars[0] == NULL || vars[1] == NULL )
548  return SCIP_INVALIDDATA;
549 
552  SCIPvarGetProbvar(vars[0]) == SCIPvarGetProbvar(vars[1]) )
553  {
554  /* if number of variable reduces, then upgrade to nonlinear constraint
555  * except if we are in the exit-presolving stage, where upgrading is not allowed
556  * in the latter case, we just do nothing, which may not be most efficient, but should still work
557  */
558  SCIP_EXPRTREE* tree;
559  SCIP_CONS* nlcons;
560 
562  return SCIP_OKAY;
563 
564  SCIP_CALL( SCIPexprtreeCopy(SCIPblkmem(scip), &tree, consdata->f) );
565 
566  for( i = 0; i < 2; ++i )
567  {
568  substexpr[i] = NULL;
569 
570  var = vars[i];
572  continue;
573 
574  coef = 1.0;
575  constant = 0.0;
576  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &coef, &constant) );
577 
578  if( coef == 0.0 )
579  {
580  /* replace var_i by constant in expression tree */
581  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &substexpr[i], SCIP_EXPR_CONST, constant) );
582  vars[i] = NULL;
583  }
584  else if( coef == 1.0 && constant == 0.0 )
585  {
586  /* do not need to change expression tree, just store new variable in tree */
587  substexpr[i] = NULL;
588  vars[i] = var;
589  }
590  else
591  {
592  /* replace var_i by coef * var_i + constant in expression tree */
593  SCIP_EXPR* child;
594 
595  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &child, SCIP_EXPR_VARIDX, i) );
596  SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &substexpr[i], 1, &child, &coef, constant) );
597  vars[i] = var;
598  }
599  }
600 
601  assert(substexpr[0] != NULL || substexpr[1] != NULL);
602 
603  SCIP_CALL( SCIPexprtreeSubstituteVars(tree, substexpr) );
604  if( substexpr[0] != NULL )
605  SCIPexprFreeDeep(SCIPblkmem(scip), &substexpr[0]);
606  if( substexpr[1] != NULL )
607  SCIPexprFreeDeep(SCIPblkmem(scip), &substexpr[1]);
608 
609  /* if variable 0 has been remove or is the same as variable 1, reindex 1 to 0 */
610  if( (vars[0] == NULL || vars[0] == vars[1]) && vars[1] != NULL )
611  {
612  int reindex[2];
613 
614  reindex[0] = 0;
615  reindex[1] = 0;
617  vars[0] = vars[1];
618  vars[1] = NULL;
619  }
620 
621  /* update variables array in tree */
622  assert(vars[1] == NULL || vars[0] != NULL);
623  SCIP_CALL( SCIPexprtreeSetVars(tree, vars[0] == NULL ? 0 : (vars[1] == NULL ? 1 : 2), vars) );
624 
625  SCIP_CALL( SCIPcreateConsNonlinear(scip, &nlcons, SCIPconsGetName(cons),
626  consdata->z != NULL ? 1 : 0, consdata->z != NULL ? &consdata->z : NULL, &consdata->zcoef,
627  1, &tree, NULL, consdata->lhs, consdata->rhs,
631  SCIPconsIsStickingAtNode(cons)) ); /*lint !e826*/
632  SCIP_CALL( SCIPaddCons(scip, nlcons) );
633  SCIPdebugMsg(scip, "upgraded to"); SCIPdebugPrintCons(scip, nlcons, NULL);
634  SCIP_CALL( SCIPreleaseCons(scip, &nlcons) );
635 
636  *isupgraded = TRUE;
637 
638  SCIP_CALL( SCIPexprtreeFree(&tree) );
639 
640  return SCIP_OKAY;
641  }
642 
643  for( i = 0; i < 2; ++i )
644  {
645  substexpr[i] = NULL;
646 
647  var = vars[i];
649  continue;
650 
651  coef = 1.0;
652  constant = 0.0;
653  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &coef, &constant) );
654  assert(coef != 0.0); /* fixed vars should have been handled above */
655 
656  if( coef == 1.0 && constant == 0.0 )
657  {
658  /* do not need to change expression tree, just store new variable in tree */
659  substexpr[i] = NULL;
660  vars[i] = var;
661  }
662  else
663  {
664  /* replace var_i by coef * var_i + constant in expression tree */
665  SCIP_EXPR* child;
666 
667  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &child, SCIP_EXPR_VARIDX, i) );
668  SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &substexpr[i], 1, &child, &coef, constant) );
669  vars[i] = var;
670  }
671 
672  /* update variables array in tree for next operation */
673  SCIP_CALL( SCIPexprtreeSetVars(consdata->f, 2, vars) );
674 
675  /* mark that variables in constraint should not be multiaggregated (bad for bound tightening and branching) */
676  if( SCIPvarIsActive(vars[0]) )
677  {
678  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, vars[0]) );
679  }
680  if( SCIPvarIsActive(vars[1]) )
681  {
682  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, vars[1]) );
683  }
684 
685  *ischanged = TRUE;
686  }
687 
688  /* update expression tree, if necessary */
689  if( substexpr[0] != NULL || substexpr[1] != NULL )
690  {
691  SCIP_CALL( SCIPexprtreeSubstituteVars(consdata->f, substexpr) );
692  if( substexpr[0] != NULL )
693  SCIPexprFreeDeep(SCIPblkmem(scip), &substexpr[0]);
694  if( substexpr[1] != NULL )
695  SCIPexprFreeDeep(SCIPblkmem(scip), &substexpr[1]);
696  }
697 
698  return SCIP_OKAY;
699 }
700 
701 /** removes fixed variables from expression graph */
702 static
704  SCIP* scip, /**< SCIP data structure */
705  SCIP_CONSHDLR* conshdlr /**< constraint handler */
706  )
707 {
708  SCIP_CONSHDLRDATA* conshdlrdata;
709  SCIP_VAR* var;
710  SCIP_VAR** vars;
711  SCIP_Real* coefs;
712  int nvars;
713  int varssize;
714  SCIP_Real constant;
715  int i;
716  int requsize;
717  SCIPdebug( int j );
718 
719  conshdlrdata = SCIPconshdlrGetData(conshdlr);
720  assert(conshdlrdata != NULL);
721  assert(conshdlrdata->exprgraph != NULL);
722 
723  if( conshdlrdata->isremovedfixings )
724  return SCIP_OKAY;
725 
726  varssize = 5;
727  SCIP_CALL( SCIPallocBufferArray(scip, &vars, varssize) );
728  SCIP_CALL( SCIPallocBufferArray(scip, &coefs, varssize) );
729 
730  i = 0;
731  while( i < SCIPexprgraphGetNVars(conshdlrdata->exprgraph) )
732  {
733  var = (SCIP_VAR*) SCIPexprgraphGetVars(conshdlrdata->exprgraph)[i];
734  if( SCIPvarIsActive(var) )
735  {
736  ++i;
737  continue;
738  }
739 
740  vars[0] = var;
741  coefs[0] = 1.0;
742  constant = 0.0;
743  nvars = 1;
744  SCIP_CALL( SCIPgetProbvarLinearSum(scip, vars, coefs, &nvars, varssize, &constant, &requsize, TRUE) );
745 
746  if( requsize > varssize )
747  {
748  SCIP_CALL( SCIPreallocBufferArray(scip, &vars, requsize) );
749  SCIP_CALL( SCIPreallocBufferArray(scip, &coefs, requsize) );
750  varssize = requsize;
751  SCIP_CALL( SCIPgetProbvarLinearSum(scip, vars, coefs, &nvars, varssize, &constant, &requsize, TRUE) );
752  assert(requsize <= varssize);
753  }
754 
755 #ifdef SCIP_DEBUG
756  SCIPdebugMsg(scip, "replace fixed variable <%s> by %g", SCIPvarGetName(var), constant);
757  for( j = 0; j < nvars; ++j )
758  {
759  SCIPdebugMsgPrint(scip, " %+g <%s>", coefs[j], SCIPvarGetName(vars[j]));
760  }
761  SCIPdebugMsgPrint(scip, "\n");
762 #endif
763 
764  SCIP_CALL( SCIPexprgraphReplaceVarByLinearSum(conshdlrdata->exprgraph, var, nvars, coefs, (void**)vars, constant) );
765 
766  i = 0;
767  }
768 
769  SCIPfreeBufferArray(scip, &vars);
770  SCIPfreeBufferArray(scip, &coefs);
771 
772  conshdlrdata->isremovedfixings = TRUE;
773 
774  return SCIP_OKAY;
775 }
776 
777 /** computes violation of a constraint */
778 static
780  SCIP* scip, /**< SCIP data structure */
781  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
782  SCIP_CONS* cons, /**< constraint */
783  SCIP_SOL* sol /**< solution or NULL if LP solution should be used */
784  )
785 { /*lint --e{666}*/
786  SCIP_CONSHDLRDATA* conshdlrdata;
787  SCIP_CONSDATA* consdata;
788  SCIP_Real xyvals[2];
789  SCIP_Real zval = 0.0;
790  SCIP_Real xlb;
791  SCIP_Real xub;
792  SCIP_Real ylb;
793  SCIP_Real yub;
794  SCIP_Real absviol;
795  SCIP_Real relviol;
796  SCIP_VAR* x;
797  SCIP_VAR* y;
798 
799  assert(scip != NULL);
800  assert(conshdlr != NULL);
801  assert(cons != NULL);
802 
803  conshdlrdata = SCIPconshdlrGetData(conshdlr);
804  assert(conshdlrdata != NULL);
805  assert(conshdlrdata->exprinterpreter != NULL);
806 
807  consdata = SCIPconsGetData(cons);
808  assert(consdata != NULL);
809 
810  if( SCIPexprtreeGetInterpreterData(consdata->f) == NULL )
811  {
812  SCIP_CALL( SCIPexprintCompile(conshdlrdata->exprinterpreter, consdata->f) );
813  }
814 
815  x = SCIPexprtreeGetVars(consdata->f)[0];
816  y = SCIPexprtreeGetVars(consdata->f)[1];
817 
818  xyvals[0] = SCIPgetSolVal(scip, sol, x);
819  xyvals[1] = SCIPgetSolVal(scip, sol, y);
820  if( consdata->z != NULL )
821  zval = SCIPgetSolVal(scip, sol, consdata->z);
822 
823  /* @todo proper handling of variables at infinity
824  * for now, just say infeasible and keep fingers crossed
825  */
826  if( SCIPisInfinity(scip, REALABS(xyvals[0])) )
827  {
828  consdata->lhsviol = consdata->rhsviol = SCIPinfinity(scip);
829  return SCIP_OKAY;
830  }
831 
832  if( SCIPisInfinity(scip, REALABS(xyvals[1])) )
833  {
834  consdata->lhsviol = consdata->rhsviol = SCIPinfinity(scip);
835  return SCIP_OKAY;
836  }
837 
838  /* project point onto box if from LP or very close to bounds to avoid eval error when function is not defined slightly outside bounds */
839  xlb = SCIPvarGetLbGlobal(x);
840  xub = SCIPvarGetUbGlobal(x);
841  ylb = SCIPvarGetLbGlobal(y);
842  yub = SCIPvarGetUbGlobal(y);
843  /* @todo handle case where variables are outside of bounds as in other constraint handlers, see also #627 */
844  if( sol == NULL )
845  {
846  assert(SCIPisFeasGE(scip, xyvals[0], xlb));
847  assert(SCIPisFeasLE(scip, xyvals[0], xub));
848  xyvals[0] = MAX(xlb, MIN(xub, xyvals[0]));
849 
850  assert(SCIPisFeasGE(scip, xyvals[1], ylb));
851  assert(SCIPisFeasLE(scip, xyvals[1], yub));
852  xyvals[1] = MAX(ylb, MIN(yub, xyvals[1]));
853 
854  if( consdata->z != NULL )
855  {
856  assert(SCIPisFeasGE(scip, zval, SCIPvarGetLbLocal(consdata->z)));
857  assert(SCIPisFeasLE(scip, zval, SCIPvarGetUbLocal(consdata->z)));
858  zval = MAX(SCIPvarGetLbLocal(consdata->z), MIN(SCIPvarGetUbLocal(consdata->z), zval));
859  }
860  }
861  else
862  {
863  if( SCIPisEQ(scip, xyvals[0], xlb) || SCIPisEQ(scip, xyvals[0], xub) )
864  xyvals[0] = MAX(xlb, MIN(xub, xyvals[0]));
865  if( SCIPisEQ(scip, xyvals[1], ylb) || SCIPisEQ(scip, xyvals[1], yub) )
866  xyvals[1] = MAX(ylb, MIN(yub, xyvals[1]));
867  }
868 
869  /* compute activity of constraint */
870  SCIP_CALL( SCIPexprintEval(conshdlrdata->exprinterpreter, consdata->f, xyvals, &consdata->activity) );
871 
872  /* point is outside the domain of f */
873  if( !SCIPisFinite(consdata->activity) )
874  {
875  consdata->lhsviol = consdata->rhsviol = SCIPinfinity(scip);
876  return SCIP_OKAY;
877  }
878 
879  if( consdata->z != NULL )
880  consdata->activity += consdata->zcoef * zval;
881 
882  /* compute violation of constraint sides */
883  absviol = 0.0;
884  relviol = 0.0;
885  if( consdata->activity < consdata->lhs && !SCIPisInfinity(scip, -consdata->lhs) )
886  {
887  consdata->lhsviol = consdata->lhs - consdata->activity;
888  absviol = consdata->lhsviol;
889  relviol = SCIPrelDiff(consdata->lhs, consdata->activity);
890  }
891  else
892  consdata->lhsviol = 0.0;
893 
894  if( consdata->activity > consdata->rhs && !SCIPisInfinity(scip, consdata->rhs) )
895  {
896  consdata->rhsviol = consdata->activity - consdata->rhs;
897  absviol = consdata->rhsviol;
898  relviol = SCIPrelDiff(consdata->activity, consdata->rhs);
899  }
900  else
901  consdata->rhsviol = 0.0;
902 
903  if( sol != NULL )
904  SCIPupdateSolConsViolation(scip, sol, absviol, relviol);
905 
906  return SCIP_OKAY;
907 }
908 
909 /** computes violation of a set of constraints */
910 static
912  SCIP* scip, /**< SCIP data structure */
913  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
914  SCIP_CONS** conss, /**< constraints */
915  int nconss, /**< number of constraints */
916  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
917  SCIP_CONS** maxviolcon /**< buffer to store constraint with largest violation, or NULL if solution is feasible */
918  )
919 {
920  SCIP_CONSDATA* consdata;
921  SCIP_Real viol;
922  SCIP_Real maxviol;
923  int c;
924 
925  assert(scip != NULL);
926  assert(conshdlr != NULL);
927  assert(conss != NULL || nconss == 0);
928  assert(maxviolcon != NULL);
929 
930  *maxviolcon = NULL;
931 
932  maxviol = 0.0;
933 
934  for( c = 0; c < nconss; ++c )
935  {
936  assert(conss != NULL);
937  assert(conss[c] != NULL);
938 
939  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol) );
940 
941  consdata = SCIPconsGetData(conss[c]);
942  assert(consdata != NULL);
943 
944  viol = MAX(consdata->lhsviol, consdata->rhsviol);
945  if( viol > maxviol && SCIPisGT(scip, viol, SCIPfeastol(scip)) )
946  {
947  maxviol = viol;
948  *maxviolcon = conss[c];
949  }
950  }
951 
952  return SCIP_OKAY;
953 }
954 
955 /** setup vred(s;x0,y0,ylb,yub) for a given f(x,y) for computing a convex-concave underestimator
956  * vred(s;x0,y0,ylb,yub) = (yub-y0)/(yub-ylb) f((yub-ylb)/(yub-y0)x0 - (y0-ylb)/(yub-y0)*s, ylb) + (y0-ylb)/(yub-ylb) f(s,yub)
957  */
958 static
960  SCIP* scip, /**< SCIP data structure */
961  SCIP_EXPRTREE** vred, /**< buffer where to store exprtree for vred */
962  SCIP_EXPRTREE* f /**< function f(x,y) for which vred should be setup */
963  )
964 {
965  SCIP_EXPR* subst[2];
966  SCIP_Real minusone;
967  SCIP_EXPR* e1;
968  SCIP_EXPR* e2;
969  SCIP_EXPR* e3;
970  SCIP_EXPR* e4;
971  SCIP_EXPR* e5;
972  SCIP_EXPR* e6;
973  SCIP_EXPR* arg1;
974  SCIP_EXPR* arg2;
975  SCIP_EXPR* vredexpr;
976 
977  assert(scip != NULL);
978  assert(vred != NULL);
979  assert(f != NULL);
980  assert(SCIPexprGetOperator(SCIPexprtreeGetRoot(f)) != SCIP_EXPR_VARIDX); /* substitute cannot substitute the root node, but f should not be a single variable anyway */
981 
982  /* setup vred(s;x0,y0,ylb,yub) for computing a convex-concave underestimator in the case where y is not at one of its bounds
983  * vred(s;x0,y0,ylb,yub) = (yub-y0)/(yub-ylb) f((yub-ylb)/(yub-y0)x0 - (y0-ylb)/(yub-y0)*s, ylb) + (y0-ylb)/(yub-ylb) f(s,yub)
984  */
985  /* create expression for x0(yub-ylb)/(yub-y0) */
986  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_PARAM, 2) ); /* ylb */
987  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_PARAM, 3) ); /* yub */
988  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e3, SCIP_EXPR_MINUS, e2, e1) ); /* yub-ylb */
989 
990  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_PARAM, 0) ); /* x0 */
991  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e3, SCIP_EXPR_MUL, e1, e3) ); /* x0(yub-ylb) */
992 
993  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_PARAM, 1) ); /* y0 */
994  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_PARAM, 3) ); /* yub */
995  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e4, SCIP_EXPR_MINUS, e2, e1) ); /* yub-y0 */
996 
997  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e5, SCIP_EXPR_DIV, e3, e4) ); /* x0(yub-ylb)/(yub-y0) */
998 
999  /* create expression for s(y0-ylb)/(yub-y0) */
1000  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_PARAM, 1) ); /* y0 */
1001  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_PARAM, 2) ); /* ylb */
1002  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e3, SCIP_EXPR_MINUS, e1, e2) ); /* y0-ylb */
1003 
1004  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_VARIDX, 0) ); /* s */
1005  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e3, SCIP_EXPR_MUL, e1, e3) ); /* s(y0-ylb) */
1006 
1007  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_PARAM, 1) ); /* y0 */
1008  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_PARAM, 3) ); /* yub */
1009  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e4, SCIP_EXPR_MINUS, e2, e1) ); /* yub-y0 */
1010 
1011  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e6, SCIP_EXPR_DIV, e3, e4) ); /* s(y0-ylb)/(yub-y0) */
1012 
1013  /* create expression for (yub-ylb)/(yub-y0)x0 - (y0-ylb)/(yub-y0)*s */
1014  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[0], SCIP_EXPR_MINUS, e5, e6) );
1015 
1016  /* create expression for ylb */
1017  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_PARAM, 2) );
1018 
1019  /* create expression for f((yub-ylb)/(yub-y0)x0 - (y0-ylb)/(yub-y0)*s, ylb) */
1021  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), arg1, subst) );
1022  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
1023  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
1024 
1025  /* create expression for f(s,yub) */
1027  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_PARAM, 3) );
1028  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), arg2, subst) );
1029  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
1030 
1031  /* create expression for (yub-y0)/(yub-ylb) */
1032  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_PARAM, 1) ); /* y0 */
1033  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_PARAM, 3) ); /* yub */
1034  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e3, SCIP_EXPR_MINUS, e2, e1) ); /* yub-y0 */
1035 
1036  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_PARAM, 2) ); /* ylb */
1037  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_PARAM, 3) ); /* yub */
1038  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e4, SCIP_EXPR_MINUS, e2, e1) ); /* yub-ylb */
1039 
1040  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e5, SCIP_EXPR_DIV, e3, e4) ); /* (yub-y0)/(yub-ylb) */
1041 
1042  /* create expression for 1 - (yub-y0)/(yub-ylb) */
1043  minusone = -1.0;
1044  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e1, e5) ); /* (yub-y0)/(yub-ylb) */
1045  SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &e6, 1, &e1, &minusone, 1.0) ); /* 1 - (yub-y0)/(yub-ylb) */
1046 
1047  /* create expression for vred = e5*arg1 + e6*arg2 */
1048  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_MUL, e5, arg1) );
1049  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_MUL, e6, arg2) );
1050  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &vredexpr, SCIP_EXPR_PLUS, e1, e2) );
1051 
1052  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), vred, vredexpr, 1, 4, NULL) );
1053 
1054  return SCIP_OKAY;
1055 }
1056 
1057 /** initializes separation data */
1058 static
1060  SCIP* scip, /**< SCIP data structure */
1061  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
1062  SCIP_CONS* cons /**< constraint */
1063  )
1064 {
1065  SCIP_CONSDATA* consdata;
1066 
1067  assert(scip != NULL);
1068  assert(exprinterpreter != NULL);
1069  assert(cons != NULL);
1070 
1071  consdata = SCIPconsGetData(cons);
1072  assert(consdata != NULL);
1073  assert(consdata->f != NULL);
1074 
1075  switch( consdata->convextype )
1076  {
1078  {
1079  SCIP_VAR** xy;
1080  SCIP_Real ref[2];
1081  SCIP_Bool sparsity[4];
1082 
1083  if( SCIPexprtreeGetInterpreterData(consdata->f) == NULL )
1084  {
1085  SCIP_CALL( SCIPexprintCompile(exprinterpreter, consdata->f) );
1086  }
1087 
1088  xy = SCIPexprtreeGetVars(consdata->f);
1089  assert(xy != NULL);
1090 
1091  /* check if the function is linear in x or y */
1092  ref[0] = MIN(MAX(SCIPvarGetLbLocal(xy[0]), 0.0), SCIPvarGetUbLocal(xy[0])); /*lint !e666*/
1093  ref[1] = MIN(MAX(SCIPvarGetLbLocal(xy[1]), 0.0), SCIPvarGetUbLocal(xy[1])); /*lint !e666*/
1094 
1095  SCIP_CALL( SCIPexprintHessianSparsityDense(exprinterpreter, consdata->f, ref, sparsity) );
1096 
1097  consdata->sepaconvexconcave.linearinx = !sparsity[0];
1098  consdata->sepaconvexconcave.lineariny = !sparsity[3];
1099 
1100  if( !consdata->sepaconvexconcave.linearinx && !SCIPisInfinity(scip, consdata->rhs) )
1101  {
1102  SCIP_EXPR* subst[2];
1103  SCIP_Real one;
1104 
1105  /* setup f(x,yfixed) for computing a convex-concave underestimator in the case where y is at one of its bounds */
1106  SCIP_CALL( SCIPexprtreeCopy(SCIPblkmem(scip), &consdata->sepaconvexconcave.f_yfixed, consdata->f) );
1107 
1108  /* x stays x, nothing to substitute
1109  * y is substituted by SCIP_EXPR_PARAM
1110  */
1111  subst[0] = NULL;
1112  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_PARAM, 0) );
1113 
1114  /* make y a parameter */
1115  SCIP_CALL( SCIPexprtreeSubstituteVars(consdata->sepaconvexconcave.f_yfixed, subst) );
1116 
1117  /* reset variables array to {x} and parameters array to {y} */
1118  one = 1.0;
1119  SCIP_CALL( SCIPexprtreeSetVars(consdata->sepaconvexconcave.f_yfixed, 1, &xy[0]) );
1120  SCIP_CALL( SCIPexprtreeSetParams(consdata->sepaconvexconcave.f_yfixed, 1, &one) );
1121 
1122  /* free subst[1] */
1123  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
1124 
1125  SCIP_CALL( SCIPexprintCompile(exprinterpreter, consdata->sepaconvexconcave.f_yfixed) );
1126 
1127  /* setup vred(s;x0,y0,ylb,yub) for computing a convex-concave underestimator in the case where y is not at one of its bounds
1128  * vred(s;x0,y0,ylb,yub) = (yub-y0)/(yub-ylb) f((yub-ylb)/(yub-y0)x0 - (y0-ylb)/(yub-y0)*s, ylb) + (y0-ylb)/(yub-ylb) f(s,yub)
1129  */
1130  SCIP_CALL( initSepaDataCreateVred(scip, &consdata->sepaconvexconcave.vred, consdata->f) );
1131  SCIP_CALL( SCIPexprintCompile(exprinterpreter, consdata->sepaconvexconcave.vred) );
1132  }
1133  else
1134  {
1135  consdata->sepaconvexconcave.f_yfixed = NULL;
1136  consdata->sepaconvexconcave.vred = NULL;
1137  }
1138 
1139  if( !consdata->sepaconvexconcave.lineariny && !SCIPisInfinity(scip, -consdata->lhs) )
1140  {
1141  /* if we have a left hand side and are not linear y in, then we may need to call
1142  * generateConvexConcaveUnderestimator for -f with swapped variables
1143  */
1144  SCIP_EXPR* minusf;
1145  SCIP_EXPR* fcopy;
1146  SCIP_VAR* vars[2];
1147  int reindex[2];
1148  SCIP_Real minusone;
1149  SCIP_Real one;
1150  SCIP_EXPR* subst[2];
1151 
1152  /* create expression for -f */
1153  minusone = -1.0;
1154  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &fcopy, SCIPexprtreeGetRoot(consdata->f)) );
1155  SCIP_CALL( SCIPexprCreateLinear(SCIPblkmem(scip), &minusf, 1, &fcopy, &minusone, 0.0) );
1156 
1157  /* reindex/swap variables */
1158  reindex[0] = 1;
1159  reindex[1] = 0;
1160  SCIPexprReindexVars(minusf, reindex);
1161 
1162  /* create expression tree for -f(y,x) */
1163  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &consdata->sepaconvexconcave.f_neg_swapped, minusf, 2, 0, NULL) );
1164 
1165  vars[0] = xy[1];
1166  vars[1] = xy[0];
1167  SCIP_CALL( SCIPexprtreeSetVars(consdata->sepaconvexconcave.f_neg_swapped, 2, vars) );
1168 
1169  SCIP_CALL( SCIPexprintCompile(exprinterpreter, consdata->sepaconvexconcave.f_neg_swapped) );
1170 
1171  /* setup -f(y, xfixed) for computing a convex-concave overestimator in the case where x is at on of it's bounds */
1172  SCIP_CALL( SCIPexprtreeCopy(SCIPblkmem(scip), &consdata->sepaconvexconcave.f_neg_swapped_yfixed, consdata->sepaconvexconcave.f_neg_swapped) );
1173 
1174  /* y stays y, nothing to substitute
1175  * x is substituted by SCIP_EXPR_PARAM
1176  */
1177  subst[0] = NULL;
1178  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_PARAM, 0) );
1179 
1180  /* make x a parameter */
1181  SCIP_CALL( SCIPexprtreeSubstituteVars(consdata->sepaconvexconcave.f_neg_swapped_yfixed, subst) );
1182 
1183  /* reset variables array to {y} and parameters array to {x} */
1184  one = 1.0;
1185  SCIP_CALL( SCIPexprtreeSetVars(consdata->sepaconvexconcave.f_neg_swapped_yfixed, 1, &xy[1]) );
1186  SCIP_CALL( SCIPexprtreeSetParams(consdata->sepaconvexconcave.f_neg_swapped_yfixed, 1, &one) );
1187 
1188  /* free subst[1] */
1189  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
1190 
1191  SCIP_CALL( SCIPexprintCompile(exprinterpreter, consdata->sepaconvexconcave.f_neg_swapped_yfixed) );
1192 
1193  /* setup vred(s;y0,x0,xlb,xub) for computing a convex-concave underestimator in the case where x is not at one of its bounds */
1194  SCIP_CALL( initSepaDataCreateVred(scip, &consdata->sepaconvexconcave.vred_neg_swapped, consdata->sepaconvexconcave.f_neg_swapped) );
1195  SCIP_CALL( SCIPexprintCompile(exprinterpreter, consdata->sepaconvexconcave.vred_neg_swapped) );
1196  }
1197  else
1198  {
1199  consdata->sepaconvexconcave.f_neg_swapped = NULL;
1200  consdata->sepaconvexconcave.f_neg_swapped_yfixed = NULL;
1201  consdata->sepaconvexconcave.vred_neg_swapped = NULL;
1202  }
1203 
1204  break;
1205  }
1206 
1207  default: ;
1208  } /*lint !e788*/
1209 
1210  return SCIP_OKAY;
1211 }
1212 
1213 /** frees separation data */
1214 static
1216  SCIP* scip, /**< SCIP data structure */
1217  SCIP_CONS* cons /**< constraint */
1218  )
1219 {
1220  SCIP_CONSDATA* consdata;
1221 
1222  assert(scip != NULL);
1223  assert(cons != NULL);
1224 
1225  consdata = SCIPconsGetData(cons);
1226  assert(consdata != NULL);
1227  assert(consdata->f != NULL);
1228 
1229  switch( consdata->convextype )
1230  {
1232  {
1233  if( consdata->sepaconvexconcave.f_yfixed != NULL )
1234  {
1235  SCIP_CALL( SCIPexprtreeFree(&consdata->sepaconvexconcave.f_yfixed) );
1236  }
1237  if( consdata->sepaconvexconcave.f_neg_swapped != NULL )
1238  {
1239  SCIP_CALL( SCIPexprtreeFree(&consdata->sepaconvexconcave.f_neg_swapped) );
1240  }
1241  if( consdata->sepaconvexconcave.f_neg_swapped_yfixed != NULL )
1242  {
1243  SCIP_CALL( SCIPexprtreeFree(&consdata->sepaconvexconcave.f_neg_swapped_yfixed) );
1244  }
1245  if( consdata->sepaconvexconcave.vred != NULL )
1246  {
1247  SCIP_CALL( SCIPexprtreeFree(&consdata->sepaconvexconcave.vred) );
1248  }
1249  if( consdata->sepaconvexconcave.vred_neg_swapped != NULL )
1250  {
1251  SCIP_CALL( SCIPexprtreeFree(&consdata->sepaconvexconcave.vred_neg_swapped) );
1252  }
1253  break;
1254  }
1255 
1256  default: ;
1257  } /*lint !e788*/
1258 
1259  return SCIP_OKAY;
1260 }
1261 
1262 /** perturbs a value w.r.t. bounds */
1263 static
1264 void perturb(
1265  SCIP_Real* val, /**< value to perturb on input; perturbed value on output */
1266  SCIP_Real lb, /**< lower bound */
1267  SCIP_Real ub, /**< upper bound */
1268  SCIP_Real amount /**< relative amount of perturbation */
1269  )
1270 {
1271  SCIP_Real range;
1272  SCIP_Real mid;
1273 
1274  assert(val != NULL);
1275 
1276  range = ub - lb;
1277  mid = 0.5 * (lb + ub);
1278 
1279  if( *val < mid )
1280  *val += MIN(1.0, amount * range);
1281  else
1282  *val -= MIN(1.0, amount * range);
1283 }
1284 
1285 /** solves an equation f'(s) = constant for a univariate convex or concave function f with respect to bounds on s
1286  * if there is no s between the bounds such that f'(s) = constant, then it returns the closest bound (and still claims success)
1287  */
1288 static
1290  SCIP* scip, /**< SCIP data structure */
1291  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
1292  SCIP_EXPRTREE* f, /**< expression tree for f(s) */
1293  SCIP_Real targetvalue, /**< target value for derivative */
1294  SCIP_Real lb, /**< lower bound on variable */
1295  SCIP_Real ub, /**< upper bound on variable */
1296  SCIP_Real* val, /**< buffer to store solution value */
1297  SCIP_Bool* success /**< buffer to indicate whether a solution has been found */
1298  )
1299 {
1300  SCIP_Real fval;
1301  SCIP_Real grad;
1302  SCIP_Real hess;
1303  SCIP_Real s;
1304  SCIP_Real nexts;
1305  SCIP_Real step;
1306  int iter;
1307 
1308  assert(scip != NULL);
1309  assert(exprinterpreter != NULL);
1310  assert(f != NULL);
1311  assert(SCIPexprtreeGetInterpreterData(f) != NULL);
1312  assert(SCIPexprtreeGetNVars(f) == 1);
1313  assert(val != NULL);
1314  assert(success != NULL);
1315 
1316  if( SCIPisEQ(scip, lb, ub) )
1317  {
1318  *val = lb;
1319  *success = TRUE;
1320  return SCIP_OKAY;
1321  }
1322 
1323  *success = FALSE;
1324 
1325  iter = 0;
1326 
1327  /* start at 0.0, projected onto interior of interval
1328  * we don't want to start at a bound, because we would not recognize if hessian is 0.0 then
1329  */
1330  s = MIN(MAX(0.0, lb), ub);
1331  perturb(&s, lb, ub, 0.1);
1332 
1333  while( ++iter < NEWTONMAXITER )
1334  {
1335  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, &s, TRUE, &fval, &grad) );
1336 
1337  /* SCIPdebugMsg(scip, "s = %.20g [%g,%g] f(s) = %g grad = %g\n", s, lb, ub, fval, grad); */
1338 
1339  if( !SCIPisFinite(grad) )
1340  {
1341  /* if f cannot be differentiated at s, perturb s to some other point close by
1342  * for that, we perturb by 0.1 * 2^{-iter}, if iter <= 65, otherwise by 1e-20
1343  * if that amount is too small to get a change in s, we increase by a factor of 2
1344  */
1345  SCIP_Real amount;
1346  SCIP_Real sold;
1347 
1348  sold = s;
1349  amount = iter <= 65 ? 0.1 / (1u<<iter) : 1e-20; /*lint !e790*/
1350  do
1351  {
1352  perturb(&s, lb, ub, amount);
1353  amount *= 2.0;
1354  } while( s == sold ); /*lint !e777*/
1355 
1356  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, &s, TRUE, &fval, &grad) );
1357 
1358  /* SCIPdebugMsg(scip, "s = %.20g [%g,%g] f(s) = %g grad = %g (perturbed by %g)\n", s, lb, ub, fval, grad, iter <= 65 ? 0.1 / (1<<iter) : 1e-20); */
1359 
1360  assert(SCIPisFinite(grad));
1361  }
1362 
1363  if( SCIPisRelEQ(scip, grad, targetvalue) )
1364  {
1365  /* if grad is targetvalue (w.r.t. epsilon), then we are done */
1366  *val = s;
1367  *success = TRUE;
1368  break;
1369  }
1370 
1371  /* coverity[callee_ptr_arith] */
1372  SCIP_CALL( SCIPexprintHessianDense(exprinterpreter, f, &s, FALSE, &fval, &hess) );
1373 
1374  /* SCIPdebugMsg(scip, "s = %.20g [%g,%g] f(s) = %g hess = %g\n", s, lb, ub, fval, hess); */
1375 
1376  if( !SCIPisFinite(hess) )
1377  {
1378  SCIP_Real smod;
1379  SCIP_Real smodval;
1380 
1381  /* if f cannot be two times differentiated at s, take the Hessian from another point close by */
1382  smod = s;
1383  perturb(&smod, lb, ub, 0.01);
1384  SCIP_CALL( SCIPexprintHessianDense(exprinterpreter, f, &smod, TRUE, &smodval, &hess) );
1385 
1386  assert(SCIPisFinite(hess));
1387  }
1388 
1389  /* next iterate would be s - (grad - targetvalue) / hess */
1390 
1391  if( SCIPisEQ(scip, s, lb) && (grad - targetvalue) * hess >= 0 )
1392  {
1393  /* if we are on the left boundary and would go left (or stay), then stop
1394  * (multiply instead of divide by hess for the case that hess is zero and since only the sign matters
1395  */
1396  *val = lb;
1397  *success = TRUE;
1398  break;
1399  }
1400 
1401  if( SCIPisEQ(scip, s, ub) && (grad - targetvalue) * hess <= 0 )
1402  {
1403  /* similar, if we are on the right boundary and would go right (or stay), then stop */
1404  *val = ub;
1405  *success = TRUE;
1406  break;
1407  }
1408 
1409  if( SCIPisZero(scip, hess) )
1410  {
1411  /* hmm, stationary point, don't know how to continue; thus, give up */
1412  break;
1413  }
1414 
1415  if( SCIPisZero(scip, (grad - targetvalue) / hess) && SCIPisFeasEQ(scip, grad, targetvalue) )
1416  {
1417  /* if grad is targetvalue (w.r.t. feastol) and step length would be almost 0, then we are also done */
1418  *val = s;
1419  *success = TRUE;
1420  break;
1421  }
1422 
1423  /* @todo we could also implement a damped Newton method if the step is too large */
1424  step = (grad - targetvalue) / hess;
1425  assert(step != 0.0);
1426 
1427  nexts = s - step;
1428  while( s == nexts ) /*lint !e777*/
1429  {
1430  /* if steplength is so tiny that there is no change in s, go by 1e-9 into given direction */
1431  step *= 2.0;
1432  nexts = s - step;
1433  }
1434  assert(nexts != s); /*lint !e777*/
1435  s = nexts;
1436 
1437  if( s < lb )
1438  s = lb;
1439  else if( s > ub )
1440  s = ub;
1441  }
1442 
1443  return SCIP_OKAY;
1444 }
1445 
1446 /** generates a cut for f(x,y) + c*z <= rhs with f(x,y) being convex or 1-convex with x or y fixed or convex-concave with y fixed
1447  * f(x0, y0) + <grad, (x,y)-(x0,y0)> + c*z <= rhs, where grad is gradient of f in (x0, y0)
1448  */
1449 static
1451  SCIP* scip, /**< SCIP data structure */
1452  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
1453  SCIP_CONS* cons, /**< constraint */
1454  SCIP_Real* x0y0, /**< value of x and y variables where to generate cut */
1455  SCIP_Bool newxy, /**< whether the last evaluation of f(x,y) with the expression interpreter was at (x0, y0) */
1456  SCIP_ROW** row /**< storage for cut */
1457  )
1458 {
1459  SCIP_VAR* x;
1460  SCIP_VAR* y;
1461  SCIP_CONSDATA* consdata;
1462  char rowname[SCIP_MAXSTRLEN];
1463  SCIP_Real fval;
1464  SCIP_Real fgrad[2];
1465  SCIP_Real rhs;
1466 
1467  assert(scip != NULL);
1468  assert(cons != NULL);
1469  assert(row != NULL);
1470 
1471  consdata = SCIPconsGetData(cons);
1472  assert(consdata != NULL);
1473  assert(!SCIPisInfinity(scip, consdata->rhs));
1474  assert(newxy || SCIPexprtreeGetInterpreterData(consdata->f) != NULL);
1475 
1476  /* compile expression if evaluated the first time; can only happen if newxy is FALSE */
1477  if( newxy && SCIPexprtreeGetInterpreterData(consdata->f) == NULL )
1478  {
1479  SCIP_CALL( SCIPexprintCompile(exprinterpreter, consdata->f) );
1480  }
1481 
1482  x = SCIPexprtreeGetVars(consdata->f)[0];
1483  y = SCIPexprtreeGetVars(consdata->f)[1];
1484 
1485  assert(consdata->convextype == SCIP_BIVAR_ALLCONVEX ||
1486  (consdata->convextype == SCIP_BIVAR_1CONVEX_INDEFINITE && (SCIPisEQ(scip, SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x)) || SCIPisEQ(scip, SCIPvarGetLbLocal(y), SCIPvarGetUbLocal(y)))) ||
1487  (consdata->convextype == SCIP_BIVAR_CONVEX_CONCAVE && SCIPisEQ(scip, SCIPvarGetLbLocal(y), SCIPvarGetUbLocal(y))) );
1488 
1489  /* compute f(x,y) and gradient of f in (x, y) */
1490  SCIP_CALL( SCIPexprintGrad(exprinterpreter, consdata->f, x0y0, newxy, &fval, fgrad) );
1491 
1492  if( !SCIPisFinite(fval) || !SCIPisFinite(fgrad[0]) || !SCIPisFinite(fgrad[1]) )
1493  {
1494  perturb(&x0y0[0], SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), 0.001);
1495  perturb(&x0y0[1], SCIPvarGetLbLocal(y), SCIPvarGetUbLocal(y), 0.001);
1496 
1497  SCIP_CALL( SCIPexprintGrad(exprinterpreter, consdata->f, x0y0, TRUE, &fval, fgrad) );
1498 
1499  if( !SCIPisFinite(fval) || !SCIPisFinite(fgrad[0]) || !SCIPisFinite(fgrad[1]) )
1500  {
1501  SCIPdebugMsg(scip, "could not evaluate f at given reference point and perturbed one");
1502  *row = NULL;
1503  return SCIP_OKAY;
1504  }
1505  }
1506 
1507  rhs = consdata->rhs - fval + fgrad[0] * x0y0[0] + fgrad[1] * x0y0[1];
1508 
1509  /* setup SCIP row */
1510  (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "%s_linearization_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
1511 
1512  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, SCIPconsGetHdlr(cons), rowname, -SCIPinfinity(scip), rhs, FALSE, FALSE /* modifiable */, TRUE /* removable */) );
1513 
1514  SCIP_CALL( SCIPaddVarsToRow(scip, *row, 2, SCIPexprtreeGetVars(consdata->f), fgrad) );
1515 
1516  if( consdata->z != NULL )
1517  SCIP_CALL( SCIPaddVarToRow(scip, *row, consdata->z, consdata->zcoef) );
1518 
1519  return SCIP_OKAY;
1520 }
1521 
1522 /** given a convex (concave, resp.) bivariate function, computes an over- (under-, resp.) estimating hyperplane
1523  * does not succeed if some variable is unbounded or both variables are fixed
1524  */
1525 static
1527  SCIP* scip, /**< SCIP data structure */
1528  SCIP_EXPRINT* exprinterpreter, /**< expression interpreter */
1529  SCIP_EXPRTREE* f, /**< bivariate function to compute under or overestimator for */
1530  SCIP_Bool doover, /**< whether to compute an overestimator (TRUE) or an underestimator (FALSE) */
1531  SCIP_Real* x0y0, /**< reference values for nonlinear variables */
1532  SCIP_Real* coefx, /**< coefficient of x in estimator */
1533  SCIP_Real* coefy, /**< coefficient of y in estimator */
1534  SCIP_Real* constant, /**< constant part of estimator */
1535  SCIP_Bool* success /**< pointer to indicate whether coefficients where successfully computed */
1536  )
1537 {
1538  SCIP_VAR* x;
1539  SCIP_VAR* y;
1540  SCIP_Real xlb;
1541  SCIP_Real xub;
1542  SCIP_Real ylb;
1543  SCIP_Real yub;
1544 
1545  SCIP_Real p1[2];
1546  SCIP_Real p2[2];
1547  SCIP_Real p3[2];
1548  SCIP_Real p4[2];
1549  SCIP_Real p1val;
1550  SCIP_Real p2val;
1551  SCIP_Real p3val;
1552  SCIP_Real p4val;
1553 
1554  SCIP_Real alpha;
1555  SCIP_Real beta;
1556  SCIP_Real gamma_;
1557  SCIP_Real delta;
1558 
1559  SCIP_Bool tryother;
1560 
1561  assert(scip != NULL);
1562  assert(exprinterpreter != NULL);
1563  assert(f != NULL);
1564  assert(x0y0 != NULL);
1565  assert(coefx != NULL);
1566  assert(coefy != NULL);
1567  assert(constant != NULL);
1568  assert(success != NULL);
1569 
1570  *success = FALSE;
1571 
1572  x = SCIPexprtreeGetVars(f)[0];
1573  y = SCIPexprtreeGetVars(f)[1];
1574 
1575  xlb = SCIPvarGetLbLocal(x);
1576  xub = SCIPvarGetUbLocal(x);
1577  ylb = SCIPvarGetLbLocal(y);
1578  yub = SCIPvarGetUbLocal(y);
1579 
1580  /* reference point should not be outside of bounds */
1581  assert(SCIPisLE(scip, xlb, x0y0[0]));
1582  assert(SCIPisGE(scip, xub, x0y0[0]));
1583  assert(SCIPisLE(scip, ylb, x0y0[1]));
1584  assert(SCIPisGE(scip, yub, x0y0[1]));
1585 
1586  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) || SCIPisInfinity(scip, -ylb) || SCIPisInfinity(scip, yub) )
1587  {
1588  SCIPdebugMsg(scip, "skip estimating hyperplane since <%s> or <%s> is unbounded\n", SCIPvarGetName(x), SCIPvarGetName(y));
1589  return SCIP_OKAY;
1590  }
1591 
1592  if( SCIPisEQ(scip, xlb, xub) && SCIPisEQ(scip, ylb, yub) )
1593  {
1594  SCIPdebugMsg(scip, "skip estimating hyperplane since both <%s> and <%s> are fixed\n", SCIPvarGetName(x), SCIPvarGetName(y));
1595  return SCIP_OKAY;
1596  }
1597 
1598  /* unten links */
1599  p1[0] = xlb;
1600  p1[1] = ylb;
1601 
1602  /* unten rechts */
1603  p2[0] = xub;
1604  p2[1] = ylb;
1605 
1606  /* oben rechts */
1607  p3[0] = xub;
1608  p3[1] = yub;
1609 
1610  /* oben links */
1611  p4[0] = xlb;
1612  p4[1] = yub;
1613 
1614  if( SCIPisEQ(scip, xlb, xub) )
1615  {
1616  /* secant between p1 and p4: p1val + [(p4val - p1val) / (yub - ylb)] * (y - ylb) */
1617  assert(!SCIPisEQ(scip, ylb, yub));
1618 
1619  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p1, &p1val) );
1620  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p4, &p4val) );
1621 
1622  if( !SCIPisFinite(p1val) || SCIPisInfinity(scip, REALABS(p1val)) || !SCIPisFinite(p4val) || SCIPisInfinity(scip, REALABS(p4val)) )
1623  {
1624  SCIPdebugMsg(scip, "skip hyperplane since function cannot be evaluated\n");
1625  return SCIP_OKAY;
1626  }
1627 
1628  *coefx = 0.0;
1629  *coefy = (p4val - p1val) / (yub - ylb);
1630  *constant = p1val - *coefy * ylb;
1631 
1632  *success = TRUE;
1633 
1634  return SCIP_OKAY;
1635  }
1636 
1637  if( SCIPisEQ(scip, ylb, yub) )
1638  {
1639  /* secant between p1 and p2: p1val + [(p2val - p1val) / (xub - xlb)] * (x - xlb) */
1640  assert(!SCIPisEQ(scip, xlb, xub));
1641 
1642  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p1, &p1val) );
1643  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p2, &p2val) );
1644 
1645  if( !SCIPisFinite(p1val) || SCIPisInfinity(scip, REALABS(p1val)) || !SCIPisFinite(p2val) || SCIPisInfinity(scip, REALABS(p2val)) )
1646  {
1647  SCIPdebugMsg(scip, "skip hyperplane since function cannot be evaluated\n");
1648  return SCIP_OKAY;
1649  }
1650 
1651  *coefx = (p2val - p1val) / (xub - xlb);
1652  *coefy = 0.0;
1653  *constant = p1val - *coefx * xlb;
1654 
1655  *success = TRUE;
1656 
1657  return SCIP_OKAY;
1658  }
1659 
1660  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p1, &p1val) );
1661  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p2, &p2val) );
1662  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p3, &p3val) );
1663  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, p4, &p4val) );
1664 
1665  /* if we want an underestimator, flip f(x,y), i.e., do as if we compute an overestimator for -f(x,y) */
1666  if( !doover )
1667  {
1668  p1val = -p1val;
1669  p2val = -p2val;
1670  p3val = -p3val;
1671  p4val = -p4val;
1672  }
1673 
1674  SCIPdebugMsg(scip, "p1 = (%g, %g), f(p1) = %g\n", p1[0], p1[1], p1val);
1675  SCIPdebugMsg(scip, "p2 = (%g, %g), f(p2) = %g\n", p2[0], p2[1], p2val);
1676  SCIPdebugMsg(scip, "p3 = (%g, %g), f(p3) = %g\n", p3[0], p3[1], p3val);
1677  SCIPdebugMsg(scip, "p4 = (%g, %g), f(p4) = %g\n", p4[0], p4[1], p4val);
1678 
1679  if( !SCIPisFinite(p1val) || SCIPisInfinity(scip, REALABS(p1val)) || !SCIPisFinite(p2val) || SCIPisInfinity(scip, REALABS(p2val)) ||
1680  ! SCIPisFinite(p3val) || SCIPisInfinity(scip, REALABS(p3val)) || !SCIPisFinite(p4val) || SCIPisInfinity(scip, REALABS(p4val)) )
1681  {
1682  SCIPdebugMsg(scip, "skip hyperplane since function cannot be evaluated\n");
1683  return SCIP_OKAY;
1684  }
1685 
1686  /* compute coefficients alpha, beta, gamma (>0), delta such that
1687  * alpha*x + beta*y + gamma*z = delta
1688  * is satisfied by at least three of the corner points (p1,f(p1)), ..., (p4,f(p4)) and
1689  * the fourth corner point lies below this hyperplane.
1690  * Since we assume that f is convex, we then know that all points (x,y,f(x,y)) are below this hyperplane, i.e.,
1691  * alpha*x + beta*y - delta <= -gamma * f(x,y),
1692  * or, equivalently,
1693  * -alpha/gamma*x - beta/gamma*y + delta/gamma >= f(x,y).
1694  */
1695 
1696  tryother = FALSE;
1697  if( x0y0[1] <= ylb + (yub - ylb)/(xub - xlb) * (x0y0[0] - xlb) )
1698  {
1699  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p2[0], p2[1], p2val, p3[0], p3[1], p3val, &alpha,
1700  &beta, &gamma_, &delta) );
1701 
1702  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p1[0] + beta * p1[1] + gamma_ * p1val, delta));
1703  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p2[0] + beta * p2[1] + gamma_ * p2val, delta));
1704  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p3[0] + beta * p3[1] + gamma_ * p3val, delta));
1705 
1706  /* if hyperplane through p1,p2,p3 does not overestimate f(p4), then it must be the other variant */
1707  if( SCIPisInfinity(scip, delta) || alpha * p4[0] + beta * p4[1] + gamma_ * p4val > delta )
1708  tryother = TRUE;
1709  }
1710  else
1711  {
1712  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p3[0], p3[1], p3val, p4[0], p4[1], p4val, &alpha,
1713  &beta, &gamma_, &delta) );
1714 
1715  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p1[0] + beta * p1[1] + gamma_ * p1val, delta));
1716  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p3[0] + beta * p3[1] + gamma_ * p3val, delta));
1717  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p4[0] + beta * p4[1] + gamma_ * p4val, delta));
1718 
1719  /* if hyperplane through p1,p3,p4 does not overestimate f(p2), then it must be the other variant */
1720  if( SCIPisInfinity(scip, delta) || alpha * p2[0] + beta * p2[1] + gamma_ * p2val > delta )
1721  tryother = TRUE;
1722  }
1723 
1724  if( tryother )
1725  {
1726  if( x0y0[1] <= yub + (ylb - yub)/(xub - xlb) * (x0y0[0] - xlb) )
1727  {
1728  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p2[0], p2[1], p2val, p4[0], p4[1], p4val,
1729  &alpha, &beta, &gamma_, &delta) );
1730 
1731  /* hyperplane should be above (p3,f(p3)) and other points should lie on hyperplane */
1732  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p1[0] + beta * p1[1] + gamma_ * p1val, delta));
1733  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p2[0] + beta * p2[1] + gamma_ * p2val, delta));
1734  assert(SCIPisInfinity(scip, delta) || SCIPisFeasLE(scip, alpha * p3[0] + beta * p3[1] + gamma_ * p3val, delta));
1735  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p4[0] + beta * p4[1] + gamma_ * p4val, delta));
1736  }
1737  else
1738  {
1739  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p2[0], p2[1], p2val, p3[0], p3[1], p3val, p4[0], p4[1], p4val,
1740  &alpha, &beta, &gamma_, &delta) );
1741 
1742  /* hyperplane should be above (p1,f(p1)) and other points should lie on hyperplane */
1743  assert(SCIPisInfinity(scip, delta) || SCIPisFeasLE(scip, alpha * p1[0] + beta * p1[1] + gamma_ * p1val, delta));
1744  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p2[0] + beta * p2[1] + gamma_ * p2val, delta));
1745  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p3[0] + beta * p3[1] + gamma_ * p3val, delta));
1746  assert(SCIPisInfinity(scip, delta) || SCIPisFeasEQ(scip, alpha * p4[0] + beta * p4[1] + gamma_ * p4val, delta));
1747  }
1748  }
1749 
1750  SCIPdebugMsg(scip, "alpha = %g, beta = %g, gamma = %g, delta = %g\n", alpha, beta, gamma_, delta);
1751 
1752  /* check if bad luck: should not happen if xlb != xub and ylb != yub and numerics are fine */
1753  if( SCIPisInfinity(scip, delta) || SCIPisZero(scip, gamma_) )
1754  return SCIP_OKAY;
1755  assert(!SCIPisNegative(scip, gamma_));
1756 
1757  /* flip hyperplane */
1758  if( !doover )
1759  gamma_ = -gamma_;
1760 
1761  *coefx = -alpha / gamma_;
1762  *coefy = -beta / gamma_;
1763  *constant = delta / gamma_;
1764 
1765  *success = TRUE;
1766 
1767  return SCIP_OKAY;
1768 }
1769 
1770 /** generates a cut for lhs <= f(x,y) + c*z with f(x,y) being convex */
1771 static
1773  SCIP* scip, /**< SCIP data structure */
1774  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
1775  SCIP_CONS* cons, /**< constraint */
1776  SCIP_Real* x0y0, /**< reference values for nonlinear variables */
1777  SCIP_ROW** row /**< storage for cut */
1778  )
1779 {
1780  SCIP_CONSDATA* consdata;
1781  SCIP_Real coefs[2];
1782  SCIP_Real constant = SCIP_INVALID;
1783  SCIP_Bool success;
1784 
1785  assert(scip != NULL);
1786  assert(cons != NULL);
1787  assert(row != NULL);
1788 
1789  *row = NULL;
1790 
1791  consdata = SCIPconsGetData(cons);
1792  assert(consdata != NULL);
1793 
1794  SCIP_CALL( generateEstimatingHyperplane(scip, exprinterpreter, consdata->f, TRUE, x0y0, &coefs[0], &coefs[1], &constant, &success) );
1795 
1796  if( success )
1797  {
1798  assert(!SCIPisInfinity(scip, -consdata->lhs));
1799  assert(SCIPisFinite(coefs[0]));
1800  assert(SCIPisFinite(coefs[1]));
1801  assert(SCIPisFinite(constant));
1802 
1803  SCIP_CALL( SCIPcreateRowCons(scip, row, SCIPconsGetHdlr(cons), "bivaroveresthyperplanecut", 0, NULL, NULL, consdata->lhs - constant, SCIPinfinity(scip), TRUE, FALSE, TRUE) );
1804 
1805  SCIP_CALL( SCIPaddVarsToRow(scip, *row, 2, SCIPexprtreeGetVars(consdata->f), coefs) );
1806  if( consdata->z != NULL )
1807  SCIP_CALL( SCIPaddVarToRow(scip, *row, consdata->z, consdata->zcoef) );
1808  }
1809  else
1810  {
1811  SCIPdebugMsg(scip, "failed to compute overestimator for all-convex constraint <%s>\n", SCIPconsGetName(cons));
1812  }
1813 
1814  return SCIP_OKAY;
1815 }
1816 
1817 /** generates a linear underestimator for f(x,y)
1818  * when the generators of the underestimating segment
1819  * are contained in y=ylb and y=yub.
1820  * Generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that
1821  * alpha * x + beta * y - delta <= gamma * f(x,y)
1822  */
1823 static
1825  SCIP* scip, /**< SCIP data structure */
1826  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
1827  SCIP_EXPRTREE* f, /**< function f(x,y) */
1828  SCIP_Real* xyref, /**< reference values for x and y */
1829  SCIP_Real cutcoeff[4], /**< cut coefficients alpha, beta, gamma, delta */
1830  SCIP_Real* convenvvalue, /**< function value of the convex envelope */
1831  SCIP_Bool* success /**< buffer to store whether coefficients were successfully computed */
1832  )
1833 {
1834  SCIP_VAR* x;
1835  SCIP_VAR* y;
1836  SCIP_Real xval;
1837  SCIP_Real xlb;
1838  SCIP_Real xub;
1839  SCIP_Real yval;
1840  SCIP_Real ylb;
1841  SCIP_Real yub;
1842 
1843  SCIP_Real t;
1844  SCIP_EXPR* vred;
1845  SCIP_EXPRTREE* vredtree;
1846  SCIP_EXPR* e1;
1847  SCIP_EXPR* e2;
1848  SCIP_EXPR* tmp;
1849  SCIP_EXPR* tmp2;
1850  SCIP_EXPR* subst[2];
1851 
1852  SCIP_Real sval;
1853  SCIP_Real slb;
1854  SCIP_Real sub;
1855  SCIP_Real rval;
1856 
1857  SCIP_Real frval;
1858  SCIP_Real fsval;
1859  SCIP_Real x0y0[2];
1860  SCIP_Real grad[2];
1861 
1862  assert(scip != NULL);
1863  assert(exprinterpreter != NULL);
1864  assert(f != NULL);
1865  assert(xyref != NULL);
1866  assert(success != NULL);
1867 
1868  x = SCIPexprtreeGetVars(f)[0];
1869  y = SCIPexprtreeGetVars(f)[1];
1870 
1871  xlb = SCIPvarGetLbLocal(x);
1872  xub = SCIPvarGetUbLocal(x);
1873 
1874  ylb = SCIPvarGetLbLocal(y);
1875  yub = SCIPvarGetUbLocal(y);
1876 
1877  xval = xyref[0];
1878  yval = xyref[1];
1879 
1880  *success = FALSE;
1881 
1882  /* check that variables are not unbounded or fixed and reference point is in interior */
1883  assert(!SCIPisInfinity(scip, -xlb));
1884  assert(!SCIPisInfinity(scip, xub));
1885  assert(!SCIPisInfinity(scip, -ylb));
1886  assert(!SCIPisInfinity(scip, yub));
1887  assert(!SCIPisEQ(scip,xlb,xub));
1888  assert(!SCIPisEQ(scip,ylb,yub));
1889  assert(!SCIPisEQ(scip,xlb,xval));
1890  assert(!SCIPisEQ(scip,xub,xval));
1891  assert(!SCIPisEQ(scip,ylb,yval));
1892  assert(!SCIPisEQ(scip,yub,yval));
1893 
1894  SCIPdebugMsg(scip, "f(%s, %s) = ", SCIPvarGetName(x), SCIPvarGetName(y));
1896  SCIPdebugMsgPrint(scip, "\n");
1897 
1898  t = (yub - yval) / (yub - ylb);
1899 
1900  /* construct v_red(s) := t f(1/t xval + (1-1/t) s, ylb) + (1-t) f(s, yub) */
1901 
1902  /* construct e1 := f(1/t xval + (1-1/t) s, ylb) */
1903  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e1, SCIPexprtreeGetRoot(f)) ); /* e1 = f(x,y) */
1904 
1905  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_VARIDX, 0) ); /* tmp = s */
1906  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp2, SCIP_EXPR_CONST, 1.0 - 1.0 / t) ); /* tmp2 = 1-1/t */
1907  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_MUL, tmp, tmp2) ); /* tmp = (1-1/t)*s */
1908  if( xval != 0.0 )
1909  {
1910  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp2, SCIP_EXPR_CONST, 1/t*xval) ); /* tmp2 = 1/t*xval */
1911  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_PLUS, tmp, tmp2) ); /* tmp = 1/t*xval + (1-1/t)*s */
1912  }
1913  subst[0] = tmp;
1914 
1915  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_CONST, ylb) ); /* tmp = ylb */
1916 
1917  assert(SCIPexprGetOperator(e1) != SCIP_EXPR_VARIDX); /* substitute cannot substitute the root node, but f should not be a single variable anyway */
1918  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e1, subst) ); /* e1 = f(1/t*xval + (1-1/t)*s, ylb) */
1919 
1920  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
1921  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
1922 
1923  /* construct e2 := f(s, yub) */
1924  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e2, SCIPexprtreeGetRoot(f)) ); /* e2 = f(x,y) */
1925 
1926  subst[0] = NULL;
1927 
1928  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_CONST, yub) );
1929 
1930  assert(SCIPexprGetOperator(e2) != SCIP_EXPR_VARIDX); /* substitute cannot substitute the root node, but f should not be a single variable anyway */
1931  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e2, subst) ); /* e2 = f(s,yub) */
1932 
1933  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
1934 
1935  /* construct vred := t * e1 + (1-t) * e2 */
1936  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, t) ); /* tmp = t */
1937  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e1, SCIP_EXPR_MUL, e1, tmp) ); /* e1 = t * f(1/t*xval+(1-1/t)*s,ylb) */
1938 
1939  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0 - t) ); /* tmp = 1 - t */
1940  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &e2, SCIP_EXPR_MUL, e2, tmp) ); /* e2 = (1-t) * f(s, yub) */
1941 
1942  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &vred, SCIP_EXPR_PLUS, e1, e2) );
1943  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &vredtree, vred, 1, 0, NULL) );
1944 
1945  SCIP_CALL( SCIPexprintCompile(exprinterpreter, vredtree) );
1946 
1947  /* compute bounds on s */
1948  slb = (yval - yub) / (ylb - yval) * (xval / t - xub);
1949  sub = (yval - yub) / (ylb - yval) * (xval / t - xlb);
1950  if( slb < xlb )
1951  slb = xlb;
1952  if( sub > xub )
1953  sub = xub;
1954 
1955  /* find s in [slb, sub] such that vred'(s) = 0 */
1956  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, vredtree, 0.0, slb, sub, &sval, success) );
1957 
1958  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
1959 
1960  if( *success == FALSE )
1961  {
1962  /* something went wrong when computing s */
1963  return SCIP_OKAY;
1964  }
1965 
1966  /* compute r from s */
1967  rval = 1.0 / t * xval + (1.0 - 1.0 / t) * sval;
1968  assert(SCIPisFeasGE(scip, rval, xlb));
1969  assert(SCIPisFeasLE(scip, rval, xub));
1970  rval = MAX(xlb, MIN(rval, xub));
1971 
1972  /* compute f(sval, yub) */
1973  x0y0[0] = sval;
1974  x0y0[1] = yub;
1975  SCIP_CALL( SCIPexprtreeEval(f, x0y0, &fsval) );
1976 
1977  /* compute f(rval, ylb) */
1978  x0y0[0] = rval;
1979  x0y0[1] = ylb;
1980  SCIP_CALL( SCIPexprtreeEval(f, x0y0, &frval) );
1981 
1982  if( !SCIPisEQ(scip, sval, xlb) && !SCIPisEQ(scip, sval, xub) )
1983  {
1984  x0y0[0] = sval;
1985  x0y0[1] = yub;
1986 
1987  /* compute f'(xbar, ybar) */
1988  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad) );
1989  }
1990  else if( !SCIPisEQ(scip, rval, xlb) && !SCIPisEQ(scip, rval, xub) )
1991  {
1992  x0y0[0] = rval;
1993  x0y0[1] = ylb;
1994 
1995  /* compute f'(xbar, ybar) */
1996  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad) );
1997  }
1998  else
1999  {
2000  /* rare case
2001  * both points (sval, yub) and (rval, ylb) should yield valid inequality
2002  * for now, just take the first one, if differentiable, otherwise second one */
2003  x0y0[0] = sval;
2004  x0y0[1] = yub;
2005 
2006  /* compute f'(xbar, ybar) */
2007  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad) );
2008 
2009  if( !SCIPisFinite(grad[0]) )
2010  {
2011  x0y0[0] = rval;
2012  x0y0[1] = ylb;
2013 
2014  /* compute f'(xbar, ybar) */
2015  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad) );
2016  }
2017  }
2018 
2019  /* compute vred(s) = t * f(rval, ylb) + (1-t) * f(s, yub) */
2020  /* SCIP_CALL( SCIPexprtreeEval(vredtree, &sval, &vredval) ); */
2021  *convenvvalue = t * frval + (1.0 - t) * fsval;
2022 
2023  SCIPdebugMsg(scip, "Parallel: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
2024  SCIPdebugMsg(scip, "Parallel: r=%g in [%g,%g], s=%g in [%g,%g], f(r,ylb)=%g, f(xlb,s)=%g\n",rval,xlb,xub,sval,ylb,yub,frval,fsval);
2025  SCIPdebugMsg(scip, "(r,ylb)=(%g,%g), (s,yub)=(%g,%g), vredval=%g\n",rval,ylb,sval,yub,*convenvvalue);
2026 
2027  if( !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
2028  {
2029  SCIPdebugMsg(scip, "f not differentiable in (x0,y0) w.r.t. x\n");
2030  return SCIP_OKAY;
2031  }
2032 
2033  /* compute cut coefficients */
2034  cutcoeff[0] = (yub - ylb) * grad[0];
2035  cutcoeff[1] = fsval - frval - (sval - rval) * grad[0];
2036  cutcoeff[2] = yub - ylb;
2037  cutcoeff[3] = cutcoeff[0] * xval + cutcoeff[1] * yval - cutcoeff[2] * *convenvvalue;
2038 
2039  SCIPdebugMsg(scip, "Parallel: cutcoeff[0]=%g, cutcoeff[1]=%g, cutcoeff[2]=1.0, cutcoeff[3]=%g\n",cutcoeff[0]/cutcoeff[2],cutcoeff[1]/cutcoeff[2],cutcoeff[3]/cutcoeff[2]);
2040 
2041  *success = TRUE;
2042 
2043  return SCIP_OKAY;
2044 }
2045 
2046 
2047 /** generates a linear underestimator for f(x,y)
2048  * with f(x,y) being convex in x and convex in y.
2049  * The segmenent connects orthogonal facets: Either (x=l_x,y=l_y)
2050  * or (x=u_x,y=u_y).
2051  * generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that
2052  * alpha * x + beta * y - delta <= gamma * f(x,y)
2053  */
2054 static
2056  SCIP* scip, /**< SCIP data structure */
2057  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
2058  SCIP_EXPRTREE* f, /**< function f(x,y) */
2059  SCIP_Real* xyref, /**< reference values for x and y */
2060  SCIP_Real cutcoeff[4], /**< cut coefficients alpha, beta, gamma, delta */
2061  SCIP_Real* convenvvalue, /**< function value of the convex envelope */
2062  SCIP_Bool* success /**< buffer to store whether coefficients were successfully computed */
2063  )
2064 {
2065  SCIP_VAR* x;
2066  SCIP_VAR* y;
2067  SCIP_Real xval;
2068  SCIP_Real xlb;
2069  SCIP_Real xub;
2070  SCIP_Real yval;
2071  SCIP_Real ylb;
2072  SCIP_Real yub;
2073 
2074  SCIP_Real x0y0[2];
2075 
2076  SCIP_EXPR* vred;
2077  SCIP_EXPRTREE* vredtree;
2078  SCIP_EXPR* e1;
2079  SCIP_EXPR* e2;
2080  SCIP_EXPR* tmp;
2081  SCIP_EXPR* expr;
2082  SCIP_EXPR* expr1;
2083  SCIP_EXPR* expr2;
2084  SCIP_EXPR* subst[2];
2085 
2086  SCIP_Real tval, tlb, tub;
2087  SCIP_Real sval;
2088  SCIP_Real rval;
2089 
2090  SCIP_Real frval,fsval;
2091  SCIP_Real grad_rval[2];
2092  SCIP_Real grad_sval[2];
2093 
2094  assert(scip != NULL);
2095  assert(exprinterpreter != NULL);
2096  assert(f != NULL);
2097  assert(convenvvalue != NULL);
2098  assert(success != NULL);
2099 
2100  x = SCIPexprtreeGetVars(f)[0];
2101  y = SCIPexprtreeGetVars(f)[1];
2102 
2103  xlb = SCIPvarGetLbLocal(x);
2104  xub = SCIPvarGetUbLocal(x);
2105 
2106  ylb = SCIPvarGetLbLocal(y);
2107  yub = SCIPvarGetUbLocal(y);
2108 
2109  xval = xyref[0];
2110  yval = xyref[1];
2111 
2112  /* check that variables are not unbounded or fixed and reference point is in interior */
2113  assert(!SCIPisInfinity(scip, -xlb));
2114  assert(!SCIPisInfinity(scip, xub));
2115  assert(!SCIPisInfinity(scip, -ylb));
2116  assert(!SCIPisInfinity(scip, yub));
2117  assert(!SCIPisEQ(scip,xlb,xub));
2118  assert(!SCIPisEQ(scip,ylb,yub));
2119  assert(!SCIPisEQ(scip,xlb,xval));
2120  assert(!SCIPisEQ(scip,xub,xval));
2121  assert(!SCIPisEQ(scip,ylb,yval));
2122  assert(!SCIPisEQ(scip,yub,yval));
2123 
2124  *success = FALSE;
2125 
2126  SCIPdebugMsg(scip, "f(%s, %s) = ", SCIPvarGetName(x), SCIPvarGetName(y));
2128  SCIPdebugMsgPrint(scip, "\n");
2129  SCIPdebugMsg(scip, "%s[%g,%g] = %g %s[%g,%g] = %g\n", SCIPvarGetName(x), xlb, xub, xval, SCIPvarGetName(y), ylb, yub, yval);
2130 
2131  /* check in which triangle the point (xval,yval) lies */
2132  if( yval <= (ylb-yub) / (xub-xlb) * (xval-xlb) + yub )
2133  {
2134  /* (xval,yval) lies in lower left triangle, i.e. region A_1 */
2135  /* construct v_red(t) := t f( xlb, (yval-(1-t)ylb)/t ) + (1-t)*f( (xval-xlb*t)/(1-t), ylb ) */
2136 
2137  /* construct e1 := f(xlb, ylb + (yval-ylb)/t) */
2138  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2139  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, yval-ylb) ); /* tmp = yval-ylb */
2140  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, tmp, expr) ); /* expr = (yval-ylb) / t */
2141  if( ylb != 0.0 )
2142  {
2143  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, ylb) ); /* tmp = ylb */
2144  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, tmp) ); /* expr = ylb + (yval-ylb) / t */
2145  }
2146  subst[1] = expr;
2147 
2148  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[0], SCIP_EXPR_CONST, xlb) ); /* subst[0] = xlb */
2149 
2150  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e1, SCIPexprtreeGetRoot(f)) ); /* e1 = f(x,y) */
2151  assert(SCIPexprGetOperator(e1) != SCIP_EXPR_VARIDX); /* expr substitute vars cannot substitute the root node, but f should not be a single variable anyway */
2152  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e1, subst) ); /* e1 = f(xlb, ylb + (yval-ylb)/t) */
2153 
2154  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2155  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2156 
2157  /* construct e2 := f((xval-xlb*t)/(1-t), ylb) */
2158  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_VARIDX, 0) ); /* expr1 = t */
2159  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2160  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MINUS, tmp, expr1) ); /* expr1 = 1-t */
2161 
2162  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_VARIDX, 0) ); /* expr2 = t */
2163  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xlb) ); /* tmp = xlb */
2164  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, expr2, tmp) ); /* expr2 = xlb * t */
2165  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xval) ); /* tmp = xval */
2166  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MINUS, tmp, expr2) ); /* expr2 = xval - xlb * t */
2167 
2168  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, expr2, expr1) ); /* expr = (xval-t*xlb)/(1-t) */
2169  subst[0] = expr;
2170 
2171  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_CONST, ylb) ); /* subst[0] = ylb */
2172 
2173  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e2, SCIPexprtreeGetRoot(f)) ); /* e2 = f(x,y) */
2174  assert(SCIPexprGetOperator(e2) != SCIP_EXPR_VARIDX); /* expr substitute vars cannot substitute the root node, but f should not be a single variable anyway */
2175  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e2, subst) ); /* e2 = f((xval-xlb*t)/(1-t), ylb) */
2176 
2177  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2178  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2179 
2180  /* construct vred := t * e1 + (1-t) * e2 */
2181  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2182  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MUL, expr, e1) ); /* expr1 = t * e1*/
2183 
2184  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2185  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2186  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MINUS, tmp, expr) ); /* expr = 1 - t */
2187  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, expr, e2) ); /* expr2 = (1-t) * e2 */
2188 
2189  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &vred, SCIP_EXPR_PLUS, expr1, expr2) );
2190  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &vredtree, vred, 1, 0, NULL) );
2191  SCIP_CALL( SCIPexprintCompile(exprinterpreter, vredtree) );
2192 
2193  /* compute bounds on t */
2194  tlb = (yval-ylb)/(yub-ylb);
2195  tub = (xub-xval)/(xub-xlb);
2196 
2197  /* find t in [lambalb, tub] such that vred'(t) = 0 */
2198  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, vredtree, 0.0, tlb, tub, &tval, success) );
2199 
2200  /* computing the cut coefficients */
2201  if( *success == FALSE )
2202  {
2203  /* something went wrong when computing s */
2204  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
2205  return SCIP_OKAY;
2206  }
2207 
2208  /* compute r and s from tval */
2209  rval = (yval-(1-tval)*ylb)/tval;
2210  rval = MAX(ylb, MIN(yub, rval));
2211  sval = (xval-xlb*tval)/(1-tval);
2212  sval = MAX(xlb, MIN(xub, sval));
2213 
2214  SCIPdebugMsg(scip, "LowerLeft: t[%g,%g] = %g -> r = %g, s = %g\n",tlb,tub,tval,rval,sval);
2215 
2216  /* compute vred(tval) */
2217  SCIP_CALL( SCIPexprtreeEval(vredtree, &tval, convenvvalue) );
2218 
2219  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
2220 
2221  /* compute f(s, ylb) and f'(s, ylb) */
2222  x0y0[0] = sval;
2223  x0y0[1] = ylb;
2224  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad_sval) );
2225 
2226  /* compute f(xlb, r) and f'(xlb,r) */
2227  x0y0[0] = xlb;
2228  x0y0[1] = rval;
2229  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad_rval) );
2230 
2231  /* generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that
2232  * alpha * x + beta * y - delta <= gamma * f(x,y)
2233  * cf. Section 2.5.2 Aux.prob. 2 case (ii)
2234  */
2235  if( !SCIPisEQ(scip, sval, xub) )
2236  {
2237  /* use the x-axis to determine the second direction */
2238  if( !SCIPisFinite(grad_sval[0]) || SCIPisInfinity(scip, REALABS(grad_sval[0])) )
2239  {
2240  *success = FALSE;
2241  return SCIP_OKAY;
2242  }
2243  cutcoeff[0] = (rval-ylb) * grad_sval[0];
2244  cutcoeff[1] = (sval-xlb) * grad_sval[0] + frval - fsval;
2245  cutcoeff[2] = rval-ylb;
2246  cutcoeff[3] = cutcoeff[0]*xlb+cutcoeff[1]*rval-cutcoeff[2]*frval;
2247  }
2248  else if( !SCIPisEQ(scip,rval,yub) )
2249  {
2250  /* use the y-axis to determine the second direction */
2251  if( !SCIPisFinite(grad_rval[1]) || SCIPisInfinity(scip, REALABS(grad_rval[1])) )
2252  {
2253  *success = FALSE;
2254  return SCIP_OKAY;
2255  }
2256  cutcoeff[0] = (rval-ylb)*grad_rval[1]+fsval-frval;
2257  cutcoeff[1] = (sval-xlb)*grad_rval[1];
2258  cutcoeff[2] = sval-xlb;
2259  cutcoeff[3] = cutcoeff[0]*xlb+cutcoeff[1]*rval-cutcoeff[2]*frval;
2260  }
2261  else
2262  {
2263  /* the point lies on the segment between (xlb,yub) and (xub,ylb) */
2264  if( !SCIPisFinite(grad_sval[0]) || !SCIPisFinite(grad_rval[0]) || SCIPisInfinity(scip, REALABS(MIN(grad_sval[0],grad_rval[0]))) )
2265  {
2266  /* FIXME maybe it is sufficient to have one of them finite, using that one for the MIN below? */
2267  *success = FALSE;
2268  return SCIP_OKAY;
2269  }
2270  cutcoeff[0] = (rval-ylb)* MIN(grad_sval[0],grad_rval[0]);
2271  cutcoeff[1] = (sval-xlb)* MIN(grad_sval[0],grad_rval[0])+frval-fsval;
2272  cutcoeff[2] = (rval-ylb);
2273  cutcoeff[3] = cutcoeff[0]*xlb+cutcoeff[1]*rval-cutcoeff[2]*frval;
2274  }
2275 
2276  SCIPdebugMsg(scip, "LowerLeft: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
2277  SCIPdebugMsg(scip, "LowerLeft: r=%g in [%g,%g], s=%g in [%g,%g], f(s,ylb)=%g, f(xlb,r)=%g\n",rval,xlb,xub,sval,ylb,yub,fsval,frval);
2278  SCIPdebugMsg(scip, "(s,ylb)=(%g,%g) (xlb,r)=(%g,%g) t=%g, vredval=%g\n",sval,ylb,xlb,rval,tval,*convenvvalue);
2279  SCIPdebugMsg(scip, "LowerLeft: cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=%g,cutcoeff[3]=%g\n",cutcoeff[0],cutcoeff[1],cutcoeff[2],cutcoeff[3]);
2280  }
2281  else
2282  {
2283  /* (xval,yval) lies in the upper right triangle, i.e region A_2 */
2284  /* construct v_red(t) := t f( xub, yub + (yval-yub)/t ) + (1-t)*f((xval-xub*t)/(1-t), yub) */
2285 
2286  /* construct e1 := f(xub, yub+(yval-yub)/t) */
2287  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t*/
2288  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, yval-yub) ); /* tmp = yval-yub*/
2289  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, tmp, expr) ); /* expr = (yval-yub) / t */
2290  if( yub != 0.0 )
2291  {
2292  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, yub) ); /* tmp = yub */
2293  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, tmp) ); /* expr = yub + (yval-yub)/t */
2294  }
2295  subst[1] = expr;
2296 
2297  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[0], SCIP_EXPR_CONST, xub) ); /* tmp = xub */
2298 
2299  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e1, SCIPexprtreeGetRoot(f)) ); /* e1 = f(x,y) */
2300  assert(SCIPexprGetOperator(e1) != SCIP_EXPR_VARIDX); /* cannot substitute root */
2301  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e1, subst) ); /* e1 = f(xub, yub + (yval-yub)/t) */
2302 
2303  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2304  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2305 
2306  /* construct e2 := f((xval-t*xub)/(1-t), yub) */
2307  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_VARIDX, 0) ); /* expr1 = t */
2308  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2309  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MINUS, tmp, expr1) ); /* expr1 = 1-t */
2310 
2311  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_VARIDX, 0) ); /* expr2 = t */
2312  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xub) ); /* tmp = xub */
2313  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, expr2, tmp) ); /* expr2 = xub * t */
2314  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xval) ); /* tmp = xval */
2315  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MINUS, tmp, expr2) ); /* expr2 = xval - xub * t */
2316 
2317  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, expr2, expr1) ); /* expr = (xval-t*xub)/(1-t) */
2318  subst[0] = expr;
2319 
2320  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_CONST, yub) ); /* tmp = yub */
2321 
2322  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e2, SCIPexprtreeGetRoot(f)) ); /* e2 = f(x,y) */
2323  assert(SCIPexprGetOperator(e2) != SCIP_EXPR_VARIDX); /* cannot substitute root */
2324  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e2, subst) ); /* e2 = f((xval-t*xub)/(1-t), yub) */
2325 
2326  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2327  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2328 
2329  /* construct vred := t * e1 + (1-t) * e2 */
2330  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2331  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MUL, e1, expr) ); /* expr1 = t * e1*/
2332 
2333  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2334  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2335  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MINUS, tmp, expr) ); /* expr = 1-t */
2336  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, e2, expr) ); /* expr2 = (1-t) * e2*/
2337 
2338  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &vred, SCIP_EXPR_PLUS, expr1, expr2) );
2339  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &vredtree, vred, 1, 0, NULL) );
2340  SCIP_CALL( SCIPexprintCompile(exprinterpreter, vredtree) );
2341 
2342  /* compute bounds on t */
2343  tlb = (yub-yval)/(yub-ylb);
2344  tub = (xval-xlb)/(xub-xlb);
2345 
2346  /* find t in [tlb, tub] such that vred'(t) = 0 */
2347  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, vredtree, 0.0, tlb, tub, &tval, success) );
2348 
2349  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
2350 
2351  if( *success == FALSE )
2352  {
2353  /* something went wrong when computing s */
2354  return SCIP_OKAY;
2355  }
2356 
2357  /* computing the cut coefficients */
2358 
2359  /* compute r and s from tval */
2360  rval = (yval-(1.0-tval)*yub)/tval;
2361  assert(SCIPisFeasGE(scip, rval, ylb));
2362  assert(SCIPisFeasLE(scip, rval, yub));
2363  rval = MAX(ylb, MIN(yub, rval));
2364 
2365  sval = (xval-xub*tval)/(1.0-tval);
2366  assert(SCIPisFeasGE(scip, sval, xlb));
2367  assert(SCIPisFeasLE(scip, sval, xub));
2368  sval = MAX(xlb, MIN(xub, sval));
2369 
2370  /* compute f(xub,r) and f'(xub,r) */
2371  x0y0[0] = xub;
2372  x0y0[1] = rval;
2373  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad_rval) );
2374 
2375  /* compute f(s,yub) and f'(s,yub) */
2376  x0y0[0] = sval;
2377  x0y0[1] = yub;
2378  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad_sval) );
2379 
2380  /* compute vred(tval) */
2381  *convenvvalue = tval * frval + (1.0-tval) * fsval;
2382 
2383  /* generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that
2384  * alpha * x + beta * y - delta <= gamma * f(x,y) */
2385 
2386  if( !SCIPisEQ(scip, sval, xlb) )
2387  {
2388  /* use the x-axis to determine the second direction */
2389  if( !SCIPisFinite(grad_sval[0]) || SCIPisInfinity(scip, REALABS(grad_sval[0])) )
2390  {
2391  *success = FALSE;
2392  return SCIP_OKAY;
2393  }
2394 
2395  cutcoeff[0] = (yub-rval)*grad_sval[0];
2396  cutcoeff[1] = (xub-sval)*grad_sval[0]+fsval-frval;
2397  cutcoeff[2] = yub-rval;
2398  cutcoeff[3] = cutcoeff[0]*sval+cutcoeff[1]*yub-cutcoeff[2]*fsval;
2399  }
2400  else if( !SCIPisEQ(scip,rval,ylb) )
2401  {
2402  /* use the y-axis to determine the second direction */
2403  if( !SCIPisFinite(grad_rval[1]) || SCIPisInfinity(scip, REALABS(grad_rval[1])) )
2404  {
2405  *success = FALSE;
2406  return SCIP_OKAY;
2407  }
2408  cutcoeff[0] = (yub-rval)*grad_rval[1]+frval-fsval;
2409  cutcoeff[1] = (xub-sval)*grad_rval[1];
2410  cutcoeff[2] = xub-sval;
2411  cutcoeff[3] = cutcoeff[0]*sval+cutcoeff[1]*yub-cutcoeff[2]*fsval;
2412  }
2413  else
2414  {
2415  /* the point lies on the segment between (xlb,yub) and (xub,ylb)
2416  * due to numerics, we get into this case here instead in the LowerLeft
2417  */
2418  assert(SCIPisFeasLE(scip, yval, (ylb-yub) / (xub-xlb) * (xval-xlb) + yub));
2419  if( !SCIPisFinite(grad_sval[0]) || !SCIPisFinite(grad_rval[0]) || SCIPisInfinity(scip, REALABS(MIN(grad_sval[0],grad_rval[0]))) )
2420  {
2421  /* FIXME maybe it is sufficient to have one of them finite, using that one for the MIN below? */
2422  *success = FALSE;
2423  return SCIP_OKAY;
2424  }
2425 
2426  cutcoeff[0] = (yub-rval)*MIN(grad_sval[0],grad_rval[0]);
2427  cutcoeff[1] = (xub-sval)*MIN(grad_sval[0],grad_rval[0])+fsval-frval;
2428  cutcoeff[2] = xub-sval;
2429  cutcoeff[3] = cutcoeff[0]*sval+cutcoeff[1]*yub-cutcoeff[2]*fsval;
2430  }
2431 
2432  SCIPdebugMsg(scip, "UpperRight: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
2433  SCIPdebugMsg(scip, "UpperRight: r=%g in [%g,%g], s=%g in [%g,%g], f(r,yub)=%g, f(xub,s)=%g\n",rval,xlb,xub,sval,ylb,yub,frval,fsval);
2434  SCIPdebugMsg(scip, "(s,yub)=(%g,%g) (xub,r)=(%g,%g) t=%g, vredval=%g\n",sval,yub,xub,rval,tval,*convenvvalue);
2435  SCIPdebugMsg(scip, "UpperRight: cutcoeff[0]=%g, cutcoeff[1]=%g, cutcoeff[2]=%g, cutcoeff[3]=%g\n",cutcoeff[0],cutcoeff[1],cutcoeff[2],cutcoeff[3]);
2436  }
2437 
2438  return SCIP_OKAY;
2439 }
2440 
2441 /** generates a linear underestimator for f(x,y) with f(x,y) being convex in x and convex in y
2442  * generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that
2443  * alpha * x + beta * y - delta <= gamma * f(x,y)
2444  */
2445 static
2447  SCIP* scip, /**< SCIP data structure */
2448  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
2449  SCIP_EXPRTREE* f, /**< function f(x,y) */
2450  SCIP_Real* xyref, /**< reference values for x and y */
2451  SCIP_Real cutcoeff[4], /**< cut coefficients alpha, beta, gamma, delta */
2452  SCIP_Real* convenvvalue, /**< function value of the convex envelope */
2453  SCIP_Bool* success /**< buffer to store whether coefficients were successfully computed */
2454  )
2455 {
2456  SCIP_VAR* x;
2457  SCIP_VAR* y;
2458  SCIP_Real xval;
2459  SCIP_Real xlb;
2460  SCIP_Real xub;
2461  SCIP_Real yval;
2462  SCIP_Real ylb;
2463  SCIP_Real yub;
2464  SCIP_Real x0y0[2];
2465 
2466  SCIP_EXPR* vred;
2467  SCIP_EXPRTREE* vredtree;
2468  SCIP_EXPR* e1;
2469  SCIP_EXPR* e2;
2470  SCIP_EXPR* tmp;
2471  SCIP_EXPR* expr;
2472  SCIP_EXPR* expr1;
2473  SCIP_EXPR* expr2;
2474  SCIP_EXPR* subst[2];
2475 
2476  SCIP_Real tval;
2477  SCIP_Real tlb;
2478  SCIP_Real tub;
2479  SCIP_Real sval;
2480  SCIP_Real rval;
2481 
2482  SCIP_Real frval;
2483  SCIP_Real fsval;
2484  SCIP_Real grad_rval[2];
2485  SCIP_Real grad_sval[2];
2486 
2487  assert(scip != NULL);
2488  assert(exprinterpreter != NULL);
2489  assert(f != NULL);
2490  assert(convenvvalue != NULL);
2491  assert(success != NULL);
2492 
2493  x = SCIPexprtreeGetVars(f)[0];
2494  y = SCIPexprtreeGetVars(f)[1];
2495 
2496  xlb = SCIPvarGetLbLocal(x);
2497  xub = SCIPvarGetUbLocal(x);
2498 
2499  ylb = SCIPvarGetLbLocal(y);
2500  yub = SCIPvarGetUbLocal(y);
2501 
2502  xval = xyref[0];
2503  yval = xyref[1];
2504 
2505  /* check that variables are not unbounded or fixed and reference point is in interior */
2506  assert(!SCIPisInfinity(scip, -xlb));
2507  assert(!SCIPisInfinity(scip, xub));
2508  assert(!SCIPisInfinity(scip, -ylb));
2509  assert(!SCIPisInfinity(scip, yub));
2510  assert(!SCIPisEQ(scip,xlb,xub));
2511  assert(!SCIPisEQ(scip,ylb,yub));
2512  assert(!SCIPisEQ(scip,xlb,xval));
2513  assert(!SCIPisEQ(scip,xub,xval));
2514  assert(!SCIPisEQ(scip,ylb,yval));
2515  assert(!SCIPisEQ(scip,yub,yval));
2516 
2517  *success = FALSE;
2518 
2519  SCIPdebugMsg(scip, "f(%s, %s) = ", SCIPvarGetName(x), SCIPvarGetName(y));
2521  SCIPdebugMsgPrint(scip, "\n");
2522 
2523  /* check in which triangle the point (xval,yval) lies */
2524  if( yval <= (yub-ylb)/(xub-xlb)*(xval-xlb)+ylb )
2525  {
2526  /* lower right triangle, i.e. region A_2 */
2527  /* construct v_red(t) := t f( xub+(xval-xub)/t, ylb ) + (1-t)*f( xub, (yval-ylb*t)/(1-t)) */
2528 
2529  /* construct e1:= f(xub+(xval-xub)/t, ylb) */
2530  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2531  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xval-xub) ); /* tmp = xval-xub */
2532  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, tmp, expr) ); /* expr = (xval-xub)/t */
2533  if( xub != 0.0 )
2534  {
2535  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xub) ); /* tmp = xub */
2536  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, tmp) ); /* expr = xub + (xval-xub)/t */
2537  }
2538  subst[0] = expr;
2539 
2540  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_CONST, ylb) ); /* subst[1] = ylb */
2541 
2542  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e1, SCIPexprtreeGetRoot(f)) ); /* e1 = f(x,y) */
2543  assert(SCIPexprGetOperator(e1) != SCIP_EXPR_VARIDX);
2544  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e1, subst) ); /* e1 = f(xub + (xval-xub)/t, ylb) */
2545 
2546  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2547  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2548 
2549  /* construct e2 := f(xub, (yval-t*ylb)/(1-t)) */
2550  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_VARIDX, 0) ); /* expr1 = t */
2551  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2552  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MINUS, tmp, expr1) ); /* expr1 = 1-t */
2553 
2554  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_VARIDX, 0) ); /* expr2 = t */
2555  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, ylb) ); /* tmp = ylb */
2556  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, expr2, tmp) ); /* expr2 = ylb * t */
2557  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, yval) ); /* tmp = yval */
2558  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MINUS, tmp, expr2) ); /* expr2 = yval - ylb * t */
2559 
2560  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, expr2, expr1) ); /* expr = (yval-t*ylb)/(1-t) */
2561  subst[1] = expr;
2562 
2563  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[0], SCIP_EXPR_CONST, xub) ); /* subst[0] = xub */
2564 
2565  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e2, SCIPexprtreeGetRoot(f)) ); /* e2 = f(x,y) */
2566  assert(SCIPexprGetOperator(e2) != SCIP_EXPR_VARIDX);
2567  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e2, subst) ); /* e2 = f(xub, (yval-t*ylb)/(1-t)) */
2568 
2569  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2570  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2571 
2572  /* construct vred := t * e1 + (1-t) * e2 */
2573  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2574  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MUL, e1, expr) ); /* expr1 = t * e1*/
2575 
2576  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2577  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2578  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MINUS, tmp, expr) ); /* expr = 1-t */
2579  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, e2, expr) ); /* expr2 = (1-t) * e2*/
2580 
2581  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &vred, SCIP_EXPR_PLUS, expr1, expr2) );
2582  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &vredtree, vred, 1, 0, NULL) );
2583  SCIP_CALL( SCIPexprintCompile(exprinterpreter, vredtree) );
2584 
2585  /* compute bounds on t */
2586  tlb = (xub-xval)/(xub-xlb);
2587  tub = (yub-yval)/(yub-ylb);
2588 
2589  /* find t in [tlb, tub] such that vred'(t) = 0 */
2590  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, vredtree, 0.0, tlb, tub, &tval, success) );
2591 
2592  if( *success == FALSE )
2593  {
2594  /* something went wrong when computing t */
2595  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
2596  return SCIP_OKAY;
2597  }
2598 
2599  /* computing the cut coefficients */
2600 
2601  /* compute r and s from tval */
2602  rval = xub+(xval-xub)/tval;
2603  rval = MAX(xlb, MIN(xub, rval));
2604  sval = (yval-tval*ylb)/(1-tval);
2605  sval = MAX(ylb, MIN(yub, sval));
2606 
2607  /* compute vred(tval) */
2608  SCIP_CALL( SCIPexprtreeEval(vredtree, &tval, convenvvalue) );
2609 
2610  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
2611 
2612  /* compute f(r, ylb) and f'(r, ylb) */
2613  x0y0[0] = rval;
2614  x0y0[1] = ylb;
2615  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad_rval) );
2616 
2617  /* compute f(xub, s) and f'(xub,s) */
2618  x0y0[0] = xub;
2619  x0y0[1] = sval;
2620  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad_sval) );
2621 
2622  /* generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that
2623  * alpha * x + beta * y - delta <= gamma * f(x,y) */
2624  if( !(SCIPisEQ(scip,rval,xlb)) )
2625  {
2626  /* take the slope along the x-axis and the slope between the points */
2627  if( !SCIPisFinite(grad_rval[0]) || SCIPisInfinity(scip, REALABS(grad_rval[0])) )
2628  {
2629  *success = FALSE;
2630  return SCIP_OKAY;
2631  }
2632  cutcoeff[0] = (sval-ylb)*grad_rval[0];
2633  cutcoeff[1] = (rval-xub)*grad_rval[0]-frval+fsval;
2634  cutcoeff[2] = sval-ylb;
2635  cutcoeff[3] = cutcoeff[0]*xub+cutcoeff[1]*sval-cutcoeff[2]*fsval;
2636  }
2637  else if( !(SCIPisEQ(scip,sval,yub)) )
2638  {
2639  /* take the slope along the y-axis and the slope between the points */
2640  if( !SCIPisFinite(grad_sval[1]) || SCIPisInfinity(scip, REALABS(grad_sval[1])) )
2641  {
2642  *success = FALSE;
2643  return SCIP_OKAY;
2644  }
2645  cutcoeff[0] = (ylb-sval)*grad_sval[1]-frval+fsval;
2646  cutcoeff[1] = (xub-rval)*grad_sval[1];
2647  cutcoeff[2] = xub-rval;
2648  cutcoeff[3] = cutcoeff[0]*xub+cutcoeff[1]*sval-cutcoeff[2]*fsval;
2649  }
2650  else
2651  {
2652  /* the point lies on the segment between (xlb,yub) and (xub,ylb) */
2653  if( !SCIPisFinite(grad_sval[0]) || !SCIPisFinite(grad_rval[0]) || SCIPisInfinity(scip, REALABS(MIN(grad_sval[0],grad_rval[0]))) )
2654  {
2655  /* FIXME maybe it is sufficient to have one of them finite, using that one for the MIN below? */
2656  *success = FALSE;
2657  return SCIP_OKAY;
2658  }
2659  cutcoeff[0] = (sval-ylb)*MIN(grad_sval[0],grad_rval[0]);
2660  cutcoeff[1] = (rval-xub)*MIN(grad_sval[0],grad_rval[0])+fsval-frval;
2661  cutcoeff[2] = sval-ylb;
2662  cutcoeff[3] = cutcoeff[0]*xub+cutcoeff[1]*sval-cutcoeff[2]*fsval;
2663  }
2664 
2665  SCIPdebugMsg(scip, "LowerRight: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
2666  SCIPdebugMsg(scip, "LowerRight: t=%g in [%g,%g], r=%g in [%g,%g], s=%g in [%g,%g]\n",tval,tlb,tub,rval,xlb,xub,sval,ylb,yub);
2667  SCIPdebugMsg(scip, "LowerRight: (r,ylb)=(%g,%g) (xub,sval)=(%g,%g) vredval=%g\n",rval,ylb,xub,sval,*convenvvalue);
2668  SCIPdebugMsg(scip, "LowerRight: cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=1.0,cutcoeff[3]=%g\n",cutcoeff[0]/cutcoeff[2],cutcoeff[1]/cutcoeff[2],cutcoeff[3]/cutcoeff[2]);
2669  }
2670  else
2671  {
2672  /* (xval,yval) lie in the upper left triangle, i.e. region A_1 */
2673  /* construct v_red(t) := t f( xlb+(xval-xlb)/t, yub ) + (1-t)*f( xlb, (yval-yub*t)/(1-t) ) */
2674 
2675  /* construct e1:= f(xlb+(xval-xlb)/t, yub) */
2676  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2677  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xval-xlb) ); /* tmp = xval-xlb */
2678  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, tmp, expr) ); /* expr = (xval-xlb)/lambda */
2679  if( xlb != 0.0 )
2680  {
2681  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, xlb) ); /* tmp = xlb */
2682  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_PLUS, expr, tmp) ); /* expr = xlb + (xval-xlb)/t */
2683  }
2684  subst[0] = expr;
2685 
2686  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_CONST, yub) ); /* subst[1] = yub */
2687 
2688  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e1, SCIPexprtreeGetRoot(f)) ); /* e1 = f(x,y) */
2689  assert(SCIPexprGetOperator(e1) != SCIP_EXPR_VARIDX);
2690  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e1, subst) ); /* e1 = f(xlb + (xval-xlb)/t, yub) */
2691 
2692  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2693  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2694 
2695  /* construct e2 := f(xlb, (yval-t*yub)/(1-t) ) */
2696  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_VARIDX, 0) ); /* expr1 = t */
2697  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2698  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MINUS, tmp, expr1) ); /* expr1 = 1-t */
2699 
2700  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_VARIDX, 0) ); /* expr2 = t */
2701  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, yub) ); /* tmp = yub */
2702  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, expr2, tmp) ); /* expr2 = yub * t */
2703  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, yval) ); /* tmp = yval */
2704  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MINUS, tmp, expr2) ); /* expr2 = yval - yub * t */
2705 
2706  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_DIV, expr2, expr1) ); /* expr = (yval-t*yub)/(1-t) */
2707  subst[1] = expr;
2708 
2709  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[0], SCIP_EXPR_CONST, xlb) ); /* subst[0] = xlb */
2710 
2711  SCIP_CALL( SCIPexprCopyDeep(SCIPblkmem(scip), &e2, SCIPexprtreeGetRoot(f)) ); /* e2 = f(x,y) */
2712  SCIP_CALL( SCIPexprSubstituteVars(SCIPblkmem(scip), e2, subst) ); /* e2 = f( xlb , (yval-t*yub)/(1-t) ) */
2713 
2714  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
2715  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
2716 
2717  /* construct vred := t * e1 + (1-t) * e2 */
2718  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2719  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr1, SCIP_EXPR_MUL, e1, expr) ); /* expr1 = t * e1*/
2720 
2721  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_VARIDX, 0) ); /* expr = t */
2722  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &tmp, SCIP_EXPR_CONST, 1.0) ); /* tmp = 1.0 */
2723  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr, SCIP_EXPR_MINUS, tmp, expr) ); /* expr = 1-t */
2724  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &expr2, SCIP_EXPR_MUL, e2, expr) ); /* expr2 = (1-t) * e2*/
2725 
2726  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &vred, SCIP_EXPR_PLUS, expr1, expr2) );
2727  SCIP_CALL( SCIPexprtreeCreate(SCIPblkmem(scip), &vredtree, vred, 1, 0, NULL) );
2728  SCIP_CALL( SCIPexprintCompile(exprinterpreter, vredtree) );
2729 
2730  /* compute bounds on lambda */
2731  tlb = (xval-xlb)/(xub-xlb);
2732  tub = (yval-ylb)/(yub-ylb);
2733 
2734  /* find t in [tlb, tub] such that vred'(t) = 0 */
2735  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, vredtree, 0.0, tlb, tub, &tval, success) );
2736 
2737  if( *success == FALSE )
2738  {
2739  /* something went wrong when computing s */
2740  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
2741  return SCIP_OKAY;
2742  }
2743 
2744  /* computing the cut coefficients */
2745 
2746  /* compute r and s from tval */
2747  rval = xlb+(xval-xlb)/tval;
2748  rval = MAX(xlb, MIN(xub, rval));
2749  sval = (yval-tval*yub)/(1-tval);
2750  sval = MAX(ylb, MIN(yub, sval));
2751 
2752  /* compute vred(tval) */
2753  SCIP_CALL( SCIPexprtreeEval(vredtree, &tval, convenvvalue) );
2754 
2755  SCIP_CALL( SCIPexprtreeFree(&vredtree) );
2756 
2757  /* compute f(r, yub) and f'(r, yub) */
2758  x0y0[0] = rval;
2759  x0y0[1] = yub;
2760  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad_rval) );
2761 
2762  /* compute f(xlb, s) and f'(xlb, s) */
2763  x0y0[0] = xlb;
2764  x0y0[1] = sval;
2765  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad_sval) );
2766 
2767  /* generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that
2768  * alpha * x + beta * y - delta <= gamma * f(x,y) */
2769  if( !SCIPisEQ(scip,rval,xub) )
2770  {
2771  /* take the slope along the x-axis and the slope between the points */
2772  if( !SCIPisFinite(grad_rval[0]) || SCIPisInfinity(scip, REALABS(grad_rval[0])) )
2773  {
2774  *success = FALSE;
2775  return SCIP_OKAY;
2776  }
2777  cutcoeff[0] = (yub-sval)*grad_rval[0];
2778  cutcoeff[1] = (xlb-rval)*grad_rval[0]-fsval+frval;
2779  cutcoeff[2] = yub-sval;
2780  cutcoeff[3] = cutcoeff[0]*xlb+cutcoeff[1]*sval-cutcoeff[2]*fsval;
2781  }
2782  else if( !SCIPisEQ(scip,sval,ylb) )
2783  {
2784  /* take the slope along the y-axis and the slope between the points */
2785  if( !SCIPisFinite(grad_sval[1]) || SCIPisInfinity(scip, REALABS(grad_sval[1])) )
2786  {
2787  *success = FALSE;
2788  return SCIP_OKAY;
2789  }
2790  cutcoeff[0] = (sval-yub)*grad_sval[1]-fsval+frval;
2791  cutcoeff[1] = (rval-xlb)*grad_sval[1];
2792  cutcoeff[2] = rval-xlb;
2793  cutcoeff[3] = cutcoeff[0]*xlb+cutcoeff[1]*sval-cutcoeff[2]*fsval;
2794  }
2795  else
2796  {
2797  /* the point lies on the segment between (xlb,yub) and (xub,ylb) */
2798  if( !SCIPisFinite(grad_sval[0]) || !SCIPisFinite(grad_rval[0]) || SCIPisInfinity(scip, REALABS(MIN(grad_rval[0],grad_sval[0]))) )
2799  {
2800  /* FIXME maybe it is sufficient to have one of them finite, using that one for the MIN below? */
2801  *success = FALSE;
2802  return SCIP_OKAY;
2803  }
2804  cutcoeff[0] = (yub-sval)*MIN(grad_rval[0],grad_sval[0]);
2805  cutcoeff[1] = (xlb-rval)*MIN(grad_rval[0],grad_sval[0])-fsval+frval;
2806  cutcoeff[2] = yub-sval;
2807  cutcoeff[3] = cutcoeff[0]*xlb+cutcoeff[1]*sval-cutcoeff[2]*fsval;
2808  }
2809 
2810  SCIPdebugMsg(scip, "UpperLeft: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
2811  SCIPdebugMsg(scip, "UpperLeft: r=%g in [%g,%g], s=%g in [%g,%g], f(r,yub)=%g, f(xlb,s)=%g\n",rval,xlb,xub,sval,ylb,yub,frval,fsval);
2812  SCIPdebugMsg(scip, "t=%g in [%g,%g], (r,yub)=(%g,%g) (xlb,sval)=(%g,%g) vredval=%g\n",tval,tlb,tub,rval,yub,xlb,sval,*convenvvalue);
2813  SCIPdebugMsg(scip, "UpperLeft: cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=1.0,cutcoeff[3]=%g\n",cutcoeff[0]/cutcoeff[2],cutcoeff[1]/cutcoeff[2],cutcoeff[3]/cutcoeff[2]);
2814  }
2815 
2816  return SCIP_OKAY;
2817 }
2818 
2819 
2820 /** generates a linear underestimator for f(x,y) with f(x,y) being STRICTLY convex in x and concave in y
2821  * generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that alpha * x + beta * y - delta <= gamma * f(x,y)
2822  */
2823 static
2825  SCIP* scip, /**< SCIP data structure */
2826  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
2827  SCIP_EXPRTREE* f, /**< function f(x,y) */
2828  SCIP_EXPRTREE* f_yfixed, /**< function f(x;y) with x variable and y parameter */
2829  SCIP_EXPRTREE* vred, /**< function vred(s;x0,y0,ylb,yub) */
2830  SCIP_Real xyref[2], /**< reference values for (x,y) */
2831  SCIP_Real cutcoeff[4], /**< cut coefficients alpha, beta, gamma, delta */
2832  SCIP_Real* convenvvalue, /**< function value of the convex envelope */
2833  SCIP_Bool* success /**< buffer to store whether coefficients were successfully computed */
2834  )
2835 {
2836  SCIP_VAR* x;
2837  SCIP_VAR* y;
2838  SCIP_Real xval;
2839  SCIP_Real xlb;
2840  SCIP_Real xub;
2841  SCIP_Real yval;
2842  SCIP_Real ylb;
2843  SCIP_Real yub;
2844 
2845  assert(scip != NULL);
2846  assert(exprinterpreter != NULL);
2847  assert(f != NULL);
2848  assert(success != NULL);
2849  assert(xyref != NULL);
2850 
2851  x = SCIPexprtreeGetVars(f)[0];
2852  y = SCIPexprtreeGetVars(f)[1];
2853 
2854  xlb = SCIPvarGetLbLocal(x);
2855  xub = SCIPvarGetUbLocal(x);
2856 
2857  ylb = SCIPvarGetLbLocal(y);
2858  yub = SCIPvarGetUbLocal(y);
2859 
2860  xval = xyref[0];
2861  yval = xyref[1];
2862 
2863  /* reference point should not be outside of bounds */
2864  assert(SCIPisLE(scip, xlb, xval));
2865  assert(SCIPisGE(scip, xub, xval));
2866  assert(SCIPisLE(scip, ylb, yval));
2867  assert(SCIPisGE(scip, yub, yval));
2868 
2869  *success = FALSE;
2870 
2871  if( SCIPisInfinity(scip, -ylb) || SCIPisInfinity(scip, yub) )
2872  {
2873  SCIPdebugMsg(scip, "skip convex-concave underestimator, since y is unbounded\n");
2874  return SCIP_OKAY;
2875  }
2876 
2877  SCIPdebugMsg(scip, "f(%s, %s) = ", SCIPvarGetName(x), SCIPvarGetName(y));
2879  SCIPdebugMsgPrint(scip, "\n");
2880 
2881  if( SCIPisEQ(scip, xlb, xub) )
2882  {
2883  /* x is fixed, so function is now concave -> generate secant between (x, ylb) and (x, yub) */
2884  SCIP_Real xy[2];
2885  SCIP_Real f_ylb;
2886  SCIP_Real f_yub;
2887  SCIP_Real slope;
2888 
2889  if( SCIPisEQ(scip, ylb, yub) )
2890  {
2891  SCIPdebugMsg(scip, "skip convex-concave underestimator, since both x and y are fixed\n");
2892  return SCIP_OKAY;
2893  }
2894 
2895  /* get f(x, ylb) */
2896  xy[0] = xlb;
2897  xy[1] = ylb;
2898  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, xy, &f_ylb) );
2899 
2900  if( !SCIPisFinite(f_ylb) )
2901  {
2902  SCIPdebugMsg(scip, "cannot evaluate function at (xlb, ylb)\n");
2903  return SCIP_OKAY;
2904  }
2905 
2906  /* get f(x, yub) */
2907  xy[1] = yub;
2908  SCIP_CALL( SCIPexprintEval(exprinterpreter, f, xy, &f_yub) );
2909 
2910  if( !SCIPisFinite(f_yub) )
2911  {
2912  SCIPdebugMsg(scip, "cannot evaluate function at (xlb, yub)\n");
2913  return SCIP_OKAY;
2914  }
2915 
2916  slope = (f_yub - f_ylb) / (yub - ylb);
2917 
2918  /* secant is f(x,ylb) + slope * (y - ylb) <= f(x,y)*/
2919 
2920  cutcoeff[0] = 0.0; /* coefficient of x == 0 */
2921  cutcoeff[1] = slope; /* coefficient of y == slope */
2922  cutcoeff[2] = 1.0; /* coefficient of f(x,y) == 1.0 */
2923  cutcoeff[3] = -(f_ylb - slope * ylb); /* constant == -(f(x,ylb) - slope * ylb) */
2924  *convenvvalue = f_ylb+slope*(yval-ylb);
2925 
2926  *success = TRUE;
2927  return SCIP_OKAY;
2928  }
2929 
2930  if( SCIPisEQ(scip, ylb, yub) )
2931  {
2932  /* y is fixed, so function is now convex -> linearize in (xval, ylb) */
2933  SCIP_Real xy[2];
2934  SCIP_Real grad[2];
2935  SCIP_Real fval;
2936 
2937  xy[0] = xval;
2938  xy[1] = ylb;
2939  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xy, TRUE, &fval, grad) );
2940 
2941  if( !SCIPisFinite(fval) || !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
2942  {
2943  perturb(&xval, xlb, xub, 0.001);
2944  xy[0] = xval;
2945 
2946  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xy, TRUE, &fval, grad) );
2947 
2948  if( !SCIPisFinite(fval) || !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
2949  {
2950  SCIPdebugMsg(scip, "cannot evaluate function or derivative in (xval,ylb), also after perturbation\n");
2951  return SCIP_OKAY;
2952  }
2953  }
2954 
2955  /* linearization is f(xval,ylb) + df/dx(xval,ylb) * (x - xval) <= f(x,y) */
2956 
2957  cutcoeff[0] = grad[0]; /* coefficient of x == gradient in x */
2958  cutcoeff[1] = 0.0; /* coefficient of y == 0 */
2959  cutcoeff[2] = 1.0; /* coefficient of f(x,y) == 1.0 */
2960  cutcoeff[3] = -(fval - grad[0] * xval); /* constant == -(f(xval,ylb) - grad * xval) */
2961  *convenvvalue = fval;
2962 
2963  *success = TRUE;
2964  return SCIP_OKAY;
2965  }
2966 
2967  /* compute coefficients of a valid underestimating hyperplane */
2968 
2969  if( SCIPisFeasEQ(scip, xlb, xval) || SCIPisFeasEQ(scip, xub, xval) )
2970  {
2971  /* x is at it's lower or upper bound */
2972  SCIP_Real x0y0[2];
2973  SCIP_Real gradylb[2];
2974  SCIP_Real gradyub[2];
2975  SCIP_Real fvalylb;
2976  SCIP_Real fvalyub;
2977 
2978  xval = SCIPisFeasEQ(scip, xlb, xval) ? xlb : xub;
2979 
2980  /* compute f'(xval, ylb) and f'(xval, yub) */
2981  x0y0[0] = xval;
2982  x0y0[1] = ylb;
2983  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fvalylb, gradylb) );
2984 
2985  x0y0[1] = yub;
2986  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fvalyub, gradyub) );
2987 
2988  if( !SCIPisFinite(gradylb[0]) || !SCIPisFinite(gradyub[0]) || !SCIPisFinite(fvalylb) || !SCIPisFinite(fvalyub) ||
2989  SCIPisInfinity(scip, REALABS(gradylb[0])) || SCIPisInfinity(scip, REALABS(gradyub[0])) )
2990  {
2991  /* move xval inside domain and continue below, hope this will work better */
2992  perturb(&xval, xlb, xub, 0.001);
2993  }
2994  else
2995  {
2996  /* setup cut coefficients */
2997  if( xval == xlb ) /*lint !e777*/
2998  cutcoeff[0] = (yub - ylb) * MIN(gradylb[0], gradyub[0]);/* coefficient of x */
2999  else
3000  cutcoeff[0] = (yub - ylb) * MAX(gradylb[0], gradyub[0]);/* coefficient of x */
3001  cutcoeff[1] = fvalyub - fvalylb; /* coefficient of y */
3002  cutcoeff[2] = yub - ylb; /* coefficient of f(x,y) */
3003  cutcoeff[3] = cutcoeff[0] * xval + cutcoeff[1] * ylb - cutcoeff[2] * fvalylb; /* constant */
3004  *convenvvalue = fvalylb;
3005 
3006  SCIPdebugMsg(scip, "alpha: %g, beta: %g, gamma: 1.0, delta: %g\n",
3007  cutcoeff[0]/cutcoeff[2], cutcoeff[1]/cutcoeff[2], cutcoeff[3]/cutcoeff[2]);
3008 
3009  *success = TRUE;
3010  return SCIP_OKAY;
3011  }
3012  }
3013 
3014  if( SCIPisFeasEQ(scip, ylb, yval) )
3015  {
3016  /* y is at it's lower bound */
3017  SCIP_Real x0y0[2];
3018  SCIP_Real grad[2];
3019  SCIP_Real xtilde;
3020  SCIP_Real fval, ftilde;
3021 
3022  /* these two cases should have been handled above */
3023  assert(!SCIPisEQ(scip, xlb, xval));
3024  assert(!SCIPisEQ(scip, xub, xval));
3025 
3026  assert(f_yfixed != NULL);
3027 
3028  /* compute f(xval, ylb) and f'(xval, ylb) */
3029  x0y0[0] = xval;
3030  x0y0[1] = ylb;
3031  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fval, grad) );
3032 
3033  if( !SCIPisFinite(fval) || !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
3034  {
3035  /* move yval inside domain and continue below, hope this will work better */
3036  perturb(&yval, ylb, yub, 0.001);
3037  }
3038  else
3039  {
3040  /* setup f(x,yub) */
3041  SCIPexprtreeSetParamVal(f_yfixed, 0, yub);
3042  SCIP_CALL( SCIPexprintNewParametrization(exprinterpreter, f_yfixed) );
3043 
3044  SCIPdebugMsg(scip, "f(x,yub) = ");
3046  SCIPdebugMsgPrint(scip, "\n");
3047 
3048  assert(SCIPexprtreeGetNVars(f_yfixed) == 1);
3049 
3050  /* find xtilde in [xlb, xub] such that f'(xtilde,yub) = f'(xval,ylb) */
3051  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, f_yfixed, grad[0], xlb, xub, &xtilde, success) );
3052 
3053  if( !*success )
3054  {
3055  SCIP_Real fxlb;
3056  SCIP_Real fxub;
3057 
3058  /* if we could not find an xtilde such that f'(xtilde,yub) = f'(xval,ylb), then probably because f'(x,yub) is constant
3059  * in this case, choose xtilde from {xlb, xub} such that it maximizes f'(xtilde, yub) - grad[0]*xtilde
3060  */
3061  /* coverity[callee_ptr_arith] */
3062  SCIP_CALL( SCIPexprintEval(exprinterpreter, f_yfixed, &xlb, &fxlb) );
3063  /* coverity[callee_ptr_arith] */
3064  SCIP_CALL( SCIPexprintEval(exprinterpreter, f_yfixed, &xub, &fxub) );
3065 
3066  SCIPdebugMsg(scip, "couldn't solve deriv equ, compare f(%g,%g) - %g*%g = %g and f(%g,%g) - %g*%g = %g\n",
3067  xlb, ylb, grad[0], xlb, fxlb - grad[0] * xlb,
3068  xub, ylb, grad[0], xub, fxub - grad[0] * xub);
3069 
3070  if( SCIPisFinite(fxlb) && SCIPisFinite(fxub) )
3071  {
3072  if( fxlb - grad[0] * xlb > fxub - grad[0] * xub )
3073  xtilde = xlb;
3074  else
3075  xtilde = xub;
3076  *success = TRUE;
3077  }
3078  else
3079  {
3080  /* move yval inside domain and continue below, hope this will work better */
3081  perturb(&yval, ylb, yub, 0.001);
3082  }
3083  }
3084 
3085  if( *success )
3086  {
3087  /* compute f(xtilde, yub) */
3088  /* coverity[callee_ptr_arith] */
3089  SCIP_CALL( SCIPexprintEval(exprinterpreter, f_yfixed, &xtilde, &ftilde) );
3090 
3091  SCIPdebugMsg(scip, "xtilde = %g, f(%g,%g) = %g\n", xtilde, xtilde, yub, ftilde);
3092 
3093  /* setup cut coefficients */
3094  cutcoeff[0] = (yub - ylb) * grad[0]; /* coefficient of x */
3095  cutcoeff[1] = ftilde - fval - grad[0] * (xtilde - xval); /* coefficient of y */
3096  cutcoeff[2] = yub - ylb; /* coefficient of f(x,y) */
3097  cutcoeff[3] = cutcoeff[0] * xval + cutcoeff[1] * ylb - cutcoeff[2] * fval; /* constant */
3098  *convenvvalue = fval;
3099 
3100  SCIPdebugMsg(scip, "alpha: %g, beta: %g, gamma: %g, delta: %g\n", cutcoeff[0], cutcoeff[1], cutcoeff[2], cutcoeff[3]);
3101 
3102  return SCIP_OKAY;
3103  }
3104  }
3105  }
3106 
3107  if( SCIPisFeasEQ(scip, yval, yub) )
3108  {
3109  /* y is at it's upper bound */
3110  SCIP_Real x0y0[2];
3111  SCIP_Real grad[2];
3112  SCIP_Real fval;
3113  SCIP_Real xtilde;
3114  SCIP_Real ftilde;
3115 
3116  assert(f_yfixed != NULL);
3117 
3118  /* compute f(xval, yub) and f'(xval, yub) */
3119  x0y0[0] = xval;
3120  x0y0[1] = yub;
3121  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fval, grad) );
3122 
3123  if( !SCIPisFinite(fval) || !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
3124  {
3125  /* move yval inside domain and continue below, hope this will work better */
3126  perturb(&yval, ylb, yub, 0.001);
3127  }
3128  else
3129  {
3130  /* setup f(x,ylb) */
3131  SCIPexprtreeSetParamVal(f_yfixed, 0, ylb);
3132  SCIP_CALL( SCIPexprintNewParametrization(exprinterpreter, f_yfixed) );
3133 
3134  assert(SCIPexprtreeGetNVars(f_yfixed) == 1);
3135 
3136  /* find xtilde in [xlb, xub] such that f'(x,ylb) = f'(xval,yub) */
3137  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, f_yfixed, grad[0], xlb, xub, &xtilde, success) );
3138 
3139  if( !*success )
3140  {
3141  SCIP_Real fxlb;
3142  SCIP_Real fxub;
3143 
3144  /* if we could not find an xtilde such that f'(xtilde,ylb) = f'(xval,yub), then probably because f'(x,ylb) is constant
3145  * in this case, choose xtilde from {xlb, xub} such that it maximizes f'(xtilde, yub) - grad[0]*xtilde
3146  */
3147  /* coverity[callee_ptr_arith] */
3148  SCIP_CALL( SCIPexprintEval(exprinterpreter, f_yfixed, &xlb, &fxlb) );
3149  /* coverity[callee_ptr_arith] */
3150  SCIP_CALL( SCIPexprintEval(exprinterpreter, f_yfixed, &xub, &fxub) );
3151 
3152  SCIPdebugMsg(scip, "couldn't solve deriv equ, compare f(%g,%g) - %g*%g = %g and f(%g,%g) - %g*%g = %g\n",
3153  xlb, yub, grad[0], xlb, fxlb - grad[0] * xlb,
3154  xub, yub, grad[0], xub, fxub - grad[0] * xub);
3155 
3156  if( SCIPisFinite(fxlb) && SCIPisFinite(fxub) )
3157  {
3158  if( fxlb - grad[0] * xlb < fxub - grad[0] * xub )
3159  xtilde = xlb;
3160  else
3161  xtilde = xub;
3162  *success = TRUE;
3163  }
3164  else
3165  {
3166  /* move yval inside domain and continue below, hope this will work better */
3167  perturb(&yval, ylb, yub, 0.001);
3168  }
3169  }
3170 
3171  if( *success )
3172  {
3173  /* compute f(xtilde, yub) */
3174  /* coverity[callee_ptr_arith] */
3175  SCIP_CALL( SCIPexprintEval(exprinterpreter, f_yfixed, &xtilde, &ftilde) );
3176 
3177  SCIPdebugMsg(scip, "xtilde = %g, f(%g,%g) = %g\n", xtilde, xtilde, ylb, ftilde);
3178 
3179  /* set up cut coefficients */
3180  cutcoeff[0] = (yub - ylb) * grad[0];
3181  cutcoeff[1] = grad[0] * (xtilde - xval) - ftilde + fval;
3182  cutcoeff[2] = yub - ylb;
3183  cutcoeff[3] = cutcoeff[0] * xval + cutcoeff[1] * yub - cutcoeff[2] * fval;
3184  *convenvvalue = fval;
3185 
3186  SCIPdebugMsg(scip, "alpha: %g, beta: %g, gamma: %g, delta: %g\n", cutcoeff[0], cutcoeff[1], cutcoeff[2], cutcoeff[3]);
3187 
3188  return SCIP_OKAY;
3189  }
3190  }
3191  }
3192 
3193  {
3194  /* x and y are somewhere between the bounds,
3195  * -> envelope is generated from f(x,y) in y=ylb and in y=yub
3196  */
3197  SCIP_Real paramvals[4];
3198 #ifdef SCIP_DEBUG
3199  const char* paramnames[4] = {"x0", "y0", "ylb", "yub"};
3200 #endif
3201  SCIP_Real t;
3202  SCIP_Real slb;
3203  SCIP_Real sub;
3204  SCIP_Real sval;
3205  SCIP_Real rval;
3206  SCIP_Real fsval;
3207  SCIP_Real frval;
3208  SCIP_Real grad[2];
3209  SCIP_Real x0y0[2];
3210 
3211  assert(vred != NULL);
3212 
3213  /* check that variables are not unbounded or fixed and reference point is in interior
3214  * @todo it should also work if x is unbounded, or? */
3215  /* assert(!SCIPisInfinity(scip, -xlb));
3216  assert(!SCIPisInfinity(scip, xub)); */
3217  assert(!SCIPisInfinity(scip, -ylb));
3218  assert(!SCIPisInfinity(scip, yub));
3219 
3220  /* update parameter values in vred */
3221  paramvals[0] = xval;
3222  paramvals[1] = yval;
3223  paramvals[2] = ylb;
3224  paramvals[3] = yub;
3225  SCIP_CALL( SCIPexprtreeSetParams(vred, 4, paramvals) );
3226  SCIP_CALL( SCIPexprintNewParametrization(exprinterpreter, vred) );
3227 
3228  SCIPdebugMsg(scip, "vred(s;x0,y0,ylb,yub) = ");
3229  SCIPdebug( SCIPexprtreePrint(vred, SCIPgetMessagehdlr(scip), NULL, NULL, paramnames) );
3230  SCIPdebugMsgPrint(scip, "\n");
3231 
3232  /* compute bounds on s */
3233  t = (yub - yval) / (yub - ylb);
3234  if( !SCIPisInfinity(scip, xub) )
3235  slb = (yval - yub) / (ylb - yval) * (xval / t - xub);
3236  else
3237  slb = -SCIPinfinity(scip);
3238  if( !SCIPisInfinity(scip, xlb) )
3239  sub = (yval - yub) / (ylb - yval) * (xval / t - xlb);
3240  else
3241  sub = SCIPinfinity(scip);
3242  if( slb < xlb )
3243  slb = xlb;
3244  if( sub > xub )
3245  sub = xub;
3246 
3247  /* find s in [slb, sub] such that vred'(s) = 0 */
3248  SCIP_CALL( solveDerivativeEquation(scip, exprinterpreter, vred, 0.0, slb, sub, &sval, success) );
3249  assert(!*success || !SCIPisInfinity(scip, REALABS(sval)));
3250 
3251  if( *success )
3252  {
3253  /* compute r from s */
3254  rval = xval / t + (1.0 - 1.0 / t) * sval;
3255  assert(SCIPisFeasGE(scip, rval, xlb));
3256  assert(SCIPisFeasLE(scip, rval, xub));
3257  rval = MAX(xlb, MIN(rval, xub));
3258 
3259  /* compute f(sval, yub) */
3260  x0y0[0] = sval;
3261  x0y0[1] = yub;
3262  SCIP_CALL( SCIPexprtreeEval(f, x0y0, &fsval) );
3263 
3264  /* compute f(rval, ylb) */
3265  x0y0[0] = rval;
3266  x0y0[1] = ylb;
3267  SCIP_CALL( SCIPexprtreeEval(f, x0y0, &frval) );
3268 
3269  if( !SCIPisEQ(scip, sval, xlb) && !SCIPisEQ(scip, sval, xub) )
3270  {
3271  x0y0[0] = sval;
3272  x0y0[1] = yub;
3273 
3274  /* compute f'(xbar, ybar) */
3275  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad) );
3276  }
3277  else if( !SCIPisEQ(scip, rval, xlb) && !SCIPisEQ(scip, rval, xub) )
3278  {
3279  x0y0[0] = rval;
3280  x0y0[1] = ylb;
3281 
3282  /* compute f'(xbar, ybar) */
3283  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad) );
3284  }
3285  else
3286  {
3287  /* rare case
3288  * both points (sval, yub) and (rval, ylb) should yield valid inequality
3289  * for now, just take the first one, if differentiable, otherwise second one
3290  */
3291  x0y0[0] = sval;
3292  x0y0[1] = yub;
3293 
3294  /* compute f'(xbar, ybar) */
3295  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fsval, grad) );
3296 
3297  if( !SCIPisFinite(grad[0]) )
3298  {
3299  x0y0[0] = rval;
3300  x0y0[1] = ylb;
3301 
3302  /* compute new f'(xbar, ybar) */
3303  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &frval, grad) );
3304  }
3305  }
3306 
3307  /* compute vred(s) = t * f(rval, ylb) + (1-t) * f(sval, yub) */
3308  *convenvvalue = t * frval + (1.0 - t) * fsval;
3309 
3310  SCIPdebugMsg(scip, "Parallel: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
3311  SCIPdebugMsg(scip, "Parallel: r=%g s=%g in [%g,%g], y in [%g,%g], f(r,ylb)=%g, f(xlb,s)=%g\n",rval,sval,xlb,xub,ylb,yub,frval,fsval);
3312  SCIPdebugMsg(scip, "(r,ylb)=(%g,%g), (s,yub)=(%g,%g), vredval=%g\n",rval,ylb,sval,yub,*convenvvalue);
3313 
3314  if( !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
3315  {
3316  SCIPdebugMsg(scip, "f not differentiable at (x0,y0) w.r.t. x\n");
3317  *success = FALSE;
3318  return SCIP_OKAY;
3319  }
3320 
3321  /* compute cut coefficients */
3322  cutcoeff[0] = (yub - ylb) * grad[0];
3323  cutcoeff[1] = fsval - frval - (sval - rval) * grad[0];
3324  cutcoeff[2] = yub - ylb;
3325  cutcoeff[3] = cutcoeff[0] * xval + cutcoeff[1] * yval - cutcoeff[2] * *convenvvalue;
3326 
3327  SCIPdebugMsg(scip, "Parallel: cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=1.0,cutcoeff[3]=%g\n",cutcoeff[0]/cutcoeff[2],cutcoeff[1]/cutcoeff[2],cutcoeff[3]/cutcoeff[2]);
3328  }
3329  }
3330 
3331  return SCIP_OKAY;
3332 }
3333 
3334 
3335 /** generates a cut for one side of lhs <= f(x,y) + c*z <= rhs with f(x,y) being convex in x and concave in y */
3336 static
3338  SCIP* scip, /**< SCIP data structure */
3339  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
3340  SCIP_CONS* cons, /**< constraint */
3341  SCIP_Real xyref[2], /**< reference values for nonlinear variables */
3342  SCIP_SIDETYPE violside, /**< for which side of constraint to find a cut */
3343  SCIP_ROW** row /**< storage for cut */
3344  )
3345 {
3346  SCIP_CONSDATA* consdata;
3347  SCIP_Real cutcoeff[4];
3348  SCIP_Real dummy;
3349  SCIP_Bool success;
3350  SCIP_Real coefs[2];
3351  char cutname[SCIP_MAXSTRLEN];
3352 
3353  assert(scip != NULL);
3354  assert(SCIPgetStage(scip) == SCIP_STAGE_SOLVING);
3355  assert(cons != NULL);
3356  assert(row != NULL);
3357 
3358  consdata = SCIPconsGetData(cons);
3359  assert(consdata != NULL);
3360  assert(consdata->f != NULL);
3361  assert(consdata->convextype == SCIP_BIVAR_CONVEX_CONCAVE);
3362 
3363  *row = NULL;
3364 
3365  SCIPdebugMsg(scip, "generate %sestimator for convex-concave constraint <%s>\n",
3366  (violside == SCIP_SIDETYPE_LEFT ? "over" : "under"), SCIPconsGetName(cons));
3367  SCIPdebugPrintCons(scip, cons, NULL);
3368 
3369  if( violside == SCIP_SIDETYPE_LEFT )
3370  {
3371  /* need overestimator */
3372  assert(!SCIPisInfinity(scip, -consdata->lhs));
3373 
3374  if( consdata->sepaconvexconcave.lineariny )
3375  {
3376  /* f is strictly convex in x and linear in y -> overestimator is polyhedral */
3377  SCIP_Real constant;
3378 
3379  SCIP_CALL( generateEstimatingHyperplane(scip, exprinterpreter, consdata->f, TRUE, xyref, &coefs[0], &coefs[1], &constant, &success) );
3380 
3381  if( success )
3382  {
3383  assert(SCIPisFinite(coefs[0]));
3384  assert(SCIPisFinite(coefs[1]));
3385  assert(SCIPisFinite(constant));
3386 
3387  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_overesthyperplanecut_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
3388  SCIP_CALL( SCIPcreateRowCons(scip, row, SCIPconsGetHdlr(cons), cutname, 0, NULL, NULL, consdata->lhs - constant, SCIPinfinity(scip), TRUE, FALSE, TRUE) );
3389 
3390  SCIP_CALL( SCIPaddVarsToRow(scip, *row, 2, SCIPexprtreeGetVars(consdata->f), coefs) );
3391  if( consdata->z != NULL )
3392  {
3393  SCIP_CALL( SCIPaddVarToRow(scip, *row, consdata->z, consdata->zcoef) );
3394  }
3395  }
3396  }
3397  else
3398  {
3399  SCIP_Real xyref_[2];
3400 
3401  /* f is strictly concave in y -> can compute overestimator by applying generateConvexConcaveUnderstimator on -f(y,x) */
3402  assert(consdata->sepaconvexconcave.f_neg_swapped != NULL);
3403 
3404  xyref_[0] = xyref[1];
3405  xyref_[1] = xyref[0];
3406  SCIP_CALL( generateConvexConcaveUnderestimator(scip, exprinterpreter, consdata->sepaconvexconcave.f_neg_swapped, consdata->sepaconvexconcave.f_neg_swapped_yfixed, consdata->sepaconvexconcave.vred_neg_swapped, xyref_, cutcoeff, &dummy, &success) );
3407 
3408  if( success )
3409  {
3410  assert(SCIPisFinite(cutcoeff[0]));
3411  assert(SCIPisFinite(cutcoeff[1]));
3412  assert(SCIPisFinite(cutcoeff[2]));
3413  assert(SCIPisFinite(cutcoeff[3]));
3414  assert(SCIPisPositive(scip, cutcoeff[2])); /* assert gamma > 0 */
3415 
3416  /* construct row from cut coefficients (alpha, beta, gamma, delta)
3417  * coefficients are such that alpha * y + beta * x - gamma * (-f(x,y)) <= delta,
3418  * i.e., gamma * f(x,y) <= delta - alpha * y - beta * x
3419  * -> lhs <= f(x,y) + c*z <= delta/gamma - alpha/gamma * y - beta/gamma * x + c*z
3420  */
3421  coefs[0] = -cutcoeff[1] / cutcoeff[2];
3422  coefs[1] = -cutcoeff[0] / cutcoeff[2];
3423  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_convexconcaveoverest_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
3424  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, SCIPconsGetHdlr(cons), cutname, consdata->lhs - cutcoeff[3]/cutcoeff[2], SCIPinfinity(scip),
3425  TRUE, FALSE /* modifiable */, TRUE /* removable */) );
3426  SCIP_CALL( SCIPaddVarsToRow(scip, *row, 2, SCIPexprtreeGetVars(consdata->f), coefs) );
3427  if( consdata->z != NULL )
3428  {
3429  SCIP_CALL( SCIPaddVarToRow(scip, *row, consdata->z, consdata->zcoef) );
3430  }
3431  }
3432  }
3433  }
3434  else
3435  {
3436  /* need underestimator */
3437  assert(violside == SCIP_SIDETYPE_RIGHT);
3438  assert(!SCIPisInfinity(scip, consdata->rhs));
3439 
3440  if( consdata->sepaconvexconcave.linearinx )
3441  {
3442  /* f is linear in x and strictly concave in y -> underestimator is polyhedral */
3443  SCIP_Real constant;
3444 
3445  SCIP_CALL( generateEstimatingHyperplane(scip, exprinterpreter, consdata->f, FALSE, xyref, &coefs[0], &coefs[1], &constant, &success) );
3446 
3447  if( success )
3448  {
3449  assert(SCIPisFinite(coefs[0]));
3450  assert(SCIPisFinite(coefs[1]));
3451  assert(SCIPisFinite(constant));
3452 
3453  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_underesthyperplanecut_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
3454  SCIP_CALL( SCIPcreateRowCons(scip, row, SCIPconsGetHdlr(cons), cutname, 0, NULL, NULL, -SCIPinfinity(scip), consdata->rhs - constant, TRUE, FALSE, TRUE) );
3455 
3456  SCIP_CALL( SCIPaddVarsToRow(scip, *row, 2, SCIPexprtreeGetVars(consdata->f), coefs) );
3457  if( consdata->z != NULL )
3458  {
3459  SCIP_CALL( SCIPaddVarToRow(scip, *row, consdata->z, consdata->zcoef) );
3460  }
3461  }
3462  }
3463  else
3464  {
3465  /* f is strictly convex in x -> can compute underestimator by applying generateConvexConcaveUnderstimator */
3466  assert(!consdata->sepaconvexconcave.linearinx); /* generateConvexConcaveUnderestimator assumes that if f is strictly convex in x */
3467 
3468  SCIP_CALL( generateConvexConcaveUnderestimator(scip, exprinterpreter, consdata->f, consdata->sepaconvexconcave.f_yfixed, consdata->sepaconvexconcave.vred, xyref, cutcoeff, &dummy, &success) );
3469 
3470  if( success )
3471  {
3472  assert(SCIPisFinite(cutcoeff[0]));
3473  assert(SCIPisFinite(cutcoeff[1]));
3474  assert(SCIPisFinite(cutcoeff[2]));
3475  assert(SCIPisFinite(cutcoeff[3]));
3476  assert(SCIPisPositive(scip, cutcoeff[2])); /* assert gamma > 0 */
3477 
3478  /* construct row from cut coefficients (alpha, beta, gamma, delta)
3479  * coefficients are such that alpha * x + beta * y - gamma * f(x,y) <= delta,
3480  * i.e., alpha/gamma * x + beta/gamma * y - delta/gamma <= f(x,y)
3481  * -> alpha/gamma * x + beta/gamma * y - delta/gamma + c*z <= f(x,y) + c*z <= rhs
3482  */
3483 
3484  coefs[0] = cutcoeff[0] / cutcoeff[2];
3485  coefs[1] = cutcoeff[1] / cutcoeff[2];
3486  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_convexconcaveunderest_%d", SCIPconsGetName(cons), SCIPgetNLPs(scip));
3487  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, SCIPconsGetHdlr(cons), cutname, -SCIPinfinity(scip), consdata->rhs + cutcoeff[3]/cutcoeff[2],
3488  TRUE, FALSE /* modifiable */, TRUE /* removable */) );
3489  SCIP_CALL( SCIPaddVarsToRow(scip, *row, 2, SCIPexprtreeGetVars(consdata->f), coefs) );
3490  if( consdata->z != NULL )
3491  {
3492  SCIP_CALL( SCIPaddVarToRow(scip, *row, consdata->z, consdata->zcoef) );
3493  }
3494  }
3495  }
3496  }
3497 
3498  return SCIP_OKAY;
3499 }
3500 
3501 
3502 /** computes an underestimating hyperplane for functions that are convex in x and y if the point to cut off lies on the boundary */
3503 static
3505  SCIP* scip, /**< SCIP data structure */
3506  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
3507  SCIP_EXPRTREE* f, /**< function f(x,y) */
3508  SCIP_Real xval, /**< current x value */
3509  SCIP_Real yval, /**< current y value */
3510  SCIP_Real xlb, /**< lower bound x */
3511  SCIP_Real xub, /**< upper bound x */
3512  SCIP_Real ylb, /**< lower bound y */
3513  SCIP_Real yub, /**< upper bound y */
3514  int min_max, /**< min=-1 max=1 */
3515  SCIP_Real cutcoeff[4], /**< returns the lifting coefficient*/
3516  SCIP_Real* convenvvalue, /**< value of the convex envelope at (xval,yval) */
3517  SCIP_Bool* success /**< buffer to indicate whether lifting was successful */
3518  )
3519 {
3520  int idx; /* indicates which variable is at the boundary */
3521 
3522  SCIP_Real mu;
3523  SCIP_Real fval;
3524  SCIP_Real grad[2];
3525 
3526  SCIP_Real x0y0[2];
3527  SCIP_Real f_lb;
3528  SCIP_Real f_ub;
3529  SCIP_Real grad_lb[2];
3530  SCIP_Real grad_ub[2];
3531 
3532  assert(SCIPisEQ(scip,xlb,xub) || SCIPisEQ(scip,ylb,yub));
3533  assert(success != NULL);
3534 
3535  *success = FALSE;
3536  idx = SCIPisEQ(scip, xlb, xub) ? 0 : 1;
3537 
3538  /* determine mu
3539  * if f is bivariate quadratic then f_x(xlb,yval) is linear in yval
3540  * thus the minimum is attained at the lower or the upper bound
3541  */
3542  x0y0[0] = xlb;
3543  x0y0[1] = ylb;
3544  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &f_lb, grad_lb) );
3545  if( !SCIPisFinite(grad_lb[idx]) )
3546  return SCIP_OKAY;
3547 
3548  x0y0[0] = xub;
3549  x0y0[1] = yub;
3550  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &f_ub, grad_ub) );
3551  if( !SCIPisFinite(grad_ub[idx]) )
3552  return SCIP_OKAY;
3553 
3554  /* if min_max=-1 choose min( grad_lb[idx], grad_ub[idx] )
3555  * if min_max= 1 choose max( grad_lb[idx], grad_ub[idx] )
3556  */
3557  if( min_max * (grad_lb[idx] - grad_ub[idx]) >= 0 )
3558  mu = grad_lb[idx];
3559  else
3560  mu = grad_ub[idx];
3561 
3562  /* determine coefficients for the hyperplane */
3563  x0y0[0] = xval;
3564  x0y0[1] = yval;
3565  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, x0y0, TRUE, &fval, grad) );
3566 
3567  if( idx == 0 )
3568  {
3569  if( !SCIPisFinite(grad[1]) || SCIPisInfinity(scip, REALABS(grad[1])) )
3570  return SCIP_OKAY;
3571  cutcoeff[0] = mu;
3572  cutcoeff[1] = grad[1];
3573  }
3574  else
3575  {
3576  assert(idx == 1);
3577  if( !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
3578  return SCIP_OKAY;
3579  cutcoeff[0] = grad[0];
3580  cutcoeff[1] = mu;
3581  }
3582  cutcoeff[2] = 1;
3583  cutcoeff[3] = -(fval-cutcoeff[0]*xval-cutcoeff[1]*yval);
3584  *convenvvalue = fval;
3585  *success = TRUE;
3586 
3587  return SCIP_OKAY;
3588 }
3589 
3590 /** generate a linear underestimator for f(x,y) with f(x,y) being convex in x and convex in y and the point to cut off lies on the boundary
3591  * generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that alpha * x + beta * y - delta <= gamma * f(x,y)
3592  */
3593 static
3595  SCIP* scip, /**< SCIP data structure */
3596  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
3597  SCIP_EXPRTREE* f, /**< function f(x,y) */
3598  SCIP_Real xyref[2], /**< reference values for x and y */
3599  SCIP_Real cutcoeff[4], /**< cut coefficients alpha, beta, gamma, delta */
3600  SCIP_Real* convenvvalue, /**< function value of the convex envelope */
3601  SCIP_Bool* success /**< buffer to store whether coefficients were successfully computed */
3602  )
3603 {
3604  SCIP_VAR* x;
3605  SCIP_VAR* y;
3606  SCIP_Real xval;
3607  SCIP_Real xlb;
3608  SCIP_Real xub;
3609  SCIP_Real yval;
3610  SCIP_Real ylb;
3611  SCIP_Real yub;
3612 
3613  assert(scip != NULL);
3614  assert(exprinterpreter != NULL);
3615  assert(f != NULL);
3616  assert(convenvvalue != NULL);
3617  assert(success != NULL);
3618 
3619  x = SCIPexprtreeGetVars(f)[0];
3620  y = SCIPexprtreeGetVars(f)[1];
3621 
3622  xlb = SCIPvarGetLbLocal(x);
3623  xub = SCIPvarGetUbLocal(x);
3624 
3625  ylb = SCIPvarGetLbLocal(y);
3626  yub = SCIPvarGetUbLocal(y);
3627 
3628  *success = FALSE;
3629 
3630  SCIPdebugMsg(scip, "f(%s, %s) = ", SCIPvarGetName(x), SCIPvarGetName(y));
3632  SCIPdebugMsgPrint(scip, "\n");
3633 
3634  xval = xyref[0];
3635  yval = xyref[1];
3636 
3637  SCIPdebugMsg(scip, "xval=%g in [%g,%g], yval=%g in [%g,%g]\n",xval,xlb,xub,yval,ylb,yub);
3638 
3639  if( SCIPisEQ(scip, ylb, yub) )
3640  {
3641  /* y is fixed, so function is now convex -> linearize in (xval, ylb) */
3642  SCIP_Real xy[2];
3643  SCIP_Real grad[2];
3644  SCIP_Real fval;
3645 
3646  xy[0] = xval;
3647  xy[1] = ylb;
3648  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xy, TRUE, &fval, grad) );
3649  if( !SCIPisFinite(grad[0]) || SCIPisInfinity(scip, REALABS(grad[0])) )
3650  return SCIP_OKAY;
3651 
3652  /* linearization is f(xval,ylb) + df/dx(xval,ylb) * (x - xval) <= f(x,y) */
3653 
3654  cutcoeff[0] = grad[0]; /* coefficient of x == gradient in x */
3655  cutcoeff[1] = 0.0; /* coefficient of y == 0 */
3656  cutcoeff[2] = 1.0; /* coefficient of f(x,y) == 1.0 */
3657  cutcoeff[3] = -(fval - grad[0] * xval); /* constant == -(f(xval,ylb) - grad * xval) */
3658 
3659  *success = TRUE;
3660  return SCIP_OKAY;
3661  }
3662 
3663  if( SCIPisEQ(scip, xlb, xub) )
3664  {
3665  /* x is fixed, so function is now convex -> linearize in (xlb, yval) */
3666  SCIP_Real xy[2];
3667  SCIP_Real grad[2];
3668  SCIP_Real fval;
3669 
3670  xy[0] = xlb;
3671  xy[1] = yval;
3672  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xy, TRUE, &fval, grad) );
3673  if( !SCIPisFinite(grad[1]) || SCIPisInfinity(scip, REALABS(grad[1])) )
3674  return SCIP_OKAY;
3675 
3676  /* linearization is f(xlb,yval) + df/dy(xlb,yval) * (y - yval) <= f(x,y) */
3677 
3678  cutcoeff[0] = 0.0; /* coefficient of x == 0.0 */
3679  cutcoeff[1] = grad[1]; /* coefficient of y == gradient in y */
3680  cutcoeff[2] = 1.0; /* coefficient of f(x,y) == 1.0 */
3681  cutcoeff[3] = -(fval - grad[1] * yval); /* constant == -(f(xlb,yval) - grad * yval) */
3682 
3683  *success = TRUE;
3684  return SCIP_OKAY;
3685  }
3686 
3687  /* check if the points lie on a boundary */
3688  if( SCIPisFeasEQ(scip, xlb, xval) )
3689  {
3690  /* apply a lifting and exploit that the function is convex in x and y
3691  * Idea: f(xlb,y) + mu (x-xlb) <= f(x,y)
3692  * determine mu with mu <= min_{x,y} ( f(x,y)-f(xlb,y) ) / (x-xlb)
3693  * f is convex in x: mu<= min_{y} f_x(xlb,y)
3694  *
3695  * mu (x-lb) + f_y(xlb,yval) * y <= f(x,y)
3696  */
3697  xval = xlb;
3698 
3699  SCIP_CALL( lifting(scip,exprinterpreter,f,xval,yval,xlb,xlb,ylb,yub,-1,cutcoeff,convenvvalue,success) );
3700 
3701  if( !*success )
3702  return SCIP_OKAY;
3703 
3704  SCIPdebugMsg(scip, "Boundary x=lb: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
3705  SCIPdebugMsg(scip, "convenvvalue = %g\n",*convenvvalue);
3706  SCIPdebugMsg(scip, "cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=%g,cutcoeff[3]=%g\n",
3707  cutcoeff[0],cutcoeff[1],cutcoeff[2],cutcoeff[3]);
3708 
3709  return SCIP_OKAY;
3710  }
3711 
3712  if( SCIPisFeasEQ(scip, ylb, yval) )
3713  {
3714  yval = ylb;
3715 
3716  SCIP_CALL( lifting(scip,exprinterpreter,f,xval,yval,xlb,xub,ylb,ylb,-1,cutcoeff,convenvvalue,success) );
3717 
3718  if( !*success )
3719  return SCIP_OKAY;
3720 
3721  SCIPdebugMsg(scip, "Boundary y=lb: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
3722  SCIPdebugMsg(scip, "convenvvalue = %g\n",*convenvvalue);
3723  SCIPdebugMsg(scip, "cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=%g,cutcoeff[3]=%g\n",
3724  cutcoeff[0],cutcoeff[1],cutcoeff[2],cutcoeff[3]);
3725 
3726  return SCIP_OKAY;
3727  }
3728 
3729  if( SCIPisFeasEQ(scip, xub, xval) )
3730  {
3731  /* apply a lifting and exploit that the function is convex in x and y
3732  * Idea: f(xlb,y) + mu (xub-x) <= f(x,y)
3733  * determine mu with mu <= min_{x,y} ( f(x,y)-f(xub,y) ) / (xub-x)
3734  * f is convex in x: -1 * mu >= min_{y} f_x(xub,y)
3735  *
3736  * mu (xub-x) + f_y(xub,yval) * y <= f(x,y)
3737  * -mu*x -mu*xub + f_y(xub,yval) * y <= f(x,y)
3738  */
3739  xval = xub;
3740 
3741  SCIP_CALL( lifting(scip,exprinterpreter,f,xval,yval,xub,xub,ylb,yub,1,cutcoeff,convenvvalue,success) );
3742 
3743  if( !*success )
3744  return SCIP_OKAY;
3745 
3746  SCIPdebugMsg(scip, "Boundary x=ub: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
3747  SCIPdebugMsg(scip, "convenvvalue = %g\n",*convenvvalue);
3748  SCIPdebugMsg(scip, "cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=%g,cutcoeff[3]=%g\n",
3749  cutcoeff[0],cutcoeff[1],cutcoeff[2],cutcoeff[3]);
3750 
3751  return SCIP_OKAY;
3752  }
3753 
3754  if( SCIPisFeasEQ(scip, yub, yval) )
3755  {
3756  yval = yub;
3757 
3758  SCIP_CALL( lifting(scip,exprinterpreter,f,xval,yval,xlb,xub,yub,yub,1,cutcoeff,convenvvalue,success) );
3759 
3760  if( !*success )
3761  return SCIP_OKAY;
3762 
3763  SCIPdebugMsg(scip, "Boundary y=ub: Cut of (xval,yval)=(%g,%g)\n",xval,yval);
3764  SCIPdebugMsg(scip, "convenvvalue = %g\n",*convenvvalue);
3765  SCIPdebugMsg(scip, "cutcoeff[0]=%g, cutcoeff[1]=%g,cutcoeff[2]=%g,cutcoeff[3]=%g\n",
3766  cutcoeff[0],cutcoeff[1],cutcoeff[2],cutcoeff[3]);
3767 
3768  return SCIP_OKAY;
3769  }
3770 
3771  /* (xval,yval) lies in the interior */
3772  SCIPerrorMessage("Tries to compute underestimator for a point at the boundary. But point is not on the boundary!\n");
3773  return SCIP_ERROR;
3774 }
3775 
3776 /** generates a linear underestimator for f(x,y) with f(x,y) being convex in x and convex in y but indefinite
3777  * This is for the case where the cone of the concave directions is (R_+ x R_-) union (R_\- x R_+).
3778  * We consider two cases:
3779  * a) the underestimating segmenent connects parallel facets
3780  * b) the underestimating segmenent connects orthogonal facets where
3781  * x=l_x, y=l_y and x=u_x, y=u_y
3782  * We ensure that the parallel facets are the horizontal with y=l_y and y=u_y
3783  * We compute the objective value of the two problems.
3784  * The smaller objective value corresponds to the convex envelope.
3785  * The supporting hyperplane is then constructed at the this point.
3786  */
3787 static
3789  SCIP* scip, /**< SCIP data structure */
3790  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
3791  SCIP_EXPRTREE* f, /**< function f(x,y) */
3792  SCIP_Real xyref[2], /**< reference values for x and y */
3793  SCIP_Real cutcoeff[4], /**< cut coefficients alpha, beta, gamma, delta */
3794  SCIP_Real* convenvvalue, /**< function value of the convex envelope */
3795  SCIP_Bool* success /**< buffer to store whether coefficients were successfully computed */
3796  )
3797 {
3798  SCIP_VAR* x;
3799  SCIP_VAR* y;
3800  SCIP_Real xlb;
3801  SCIP_Real xub;
3802  SCIP_Real ylb;
3803  SCIP_Real yub;
3804  SCIP_Real xub_ylb[2];
3805  SCIP_Real xlb_yub[2];
3806  SCIP_Real grad_xub_ylb[2];
3807  SCIP_Real grad_xlb_yub[2];
3808  SCIP_Real fval_xub_ylb;
3809  SCIP_Real fval_xlb_yub;
3810 
3811  SCIP_Real all_cutcoeff[2][4];
3812  SCIP_Real all_convenvvalue[2];
3813  SCIP_Bool all_success[2];
3814 
3815  SCIP_Real lowest;
3816  int lowestidx;
3817  int i;
3818 
3819  SCIP_EXPRTREE* fswapped;
3820  SCIP_VAR* vars[2];
3821  SCIP_Bool swapped;
3822  SCIP_Real swap_buffer;
3823  SCIP_EXPR* subst[2];
3824 
3825  assert(scip != NULL);
3826  assert(exprinterpreter != NULL);
3827  assert(f != NULL);
3828  assert(convenvvalue != NULL);
3829  assert(success != NULL);
3830 
3831  x = SCIPexprtreeGetVars(f)[0];
3832  y = SCIPexprtreeGetVars(f)[1];
3833 
3834  xlb = SCIPvarGetLbLocal(x);
3835  xub = SCIPvarGetUbLocal(x);
3836 
3837  ylb = SCIPvarGetLbLocal(y);
3838  yub = SCIPvarGetUbLocal(y);
3839 
3840  *success = FALSE;
3841 
3842  xub_ylb[0] = xub;
3843  xub_ylb[1] = ylb;
3844  xlb_yub[0] = xlb;
3845  xlb_yub[1] = yub;
3846 
3847  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xub_ylb, TRUE, &fval_xub_ylb, grad_xub_ylb) );
3848  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xlb_yub, TRUE, &fval_xlb_yub, grad_xlb_yub) );
3849 
3850  if( !SCIPisFinite(fval_xub_ylb) || SCIPisInfinity(scip, REALABS(fval_xub_ylb)) || !SCIPisFinite(fval_xlb_yub) || SCIPisInfinity(scip, REALABS(fval_xlb_yub)) )
3851  {
3852  SCIPdebugMsg(scip, "skip 1-convex underestimator since function cannot be evaluated\n");
3853  return SCIP_OKAY;
3854  }
3855 
3856  if( !SCIPisFinite(grad_xub_ylb[0]) || !SCIPisFinite(grad_xlb_yub[1]) )
3857  {
3858  SCIPdebugMsg(scip, "skip 1-convex underestimator since function cannot be differentiated\n");
3859  return SCIP_OKAY;
3860  }
3861 
3862  SCIPdebugMsg(scip, "f(%s, %s) = ", SCIPvarGetName(x), SCIPvarGetName(y));
3864  SCIPdebugMsgPrint(scip, "\n");
3865 
3866  SCIPdebugMsg(scip, "xval=%g in [%g,%g], yval=%g in [%g,%g]\n", xyref[0], xlb, xub, xyref[1], ylb, yub);
3867 
3868  /* assure (xub-xlb)*f_x(xub,ylb) - (yub-ylb)*f_y(xlb,yub) >= f(xub,ylb) - f(xlb,yub) */
3869  /* f_y(xlb,yub)*(ylb-yub)* + f(xlb,yub) >= f_x(xub,ylb)*(xub-xlb) + f(xub,ylb) */
3870  if( fval_xub_ylb-fval_xlb_yub <= (xub-xlb)*grad_xub_ylb[0]-(yub-ylb)*grad_xlb_yub[1] )
3871  {
3872  swapped = 0;
3873  }
3874  else
3875  {
3876  /* swap the variables */
3877  swapped = 1;
3878 
3879  vars[0] = SCIPexprtreeGetVars(f)[1];
3880  vars[1] = SCIPexprtreeGetVars(f)[0];
3881 
3882  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[0], SCIP_EXPR_VARIDX, 1) );
3883  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_VARIDX, 0) );
3884 
3885  SCIP_CALL( SCIPexprtreeCopy(SCIPblkmem(scip), &fswapped, f) );
3886  SCIP_CALL( SCIPexprtreeSubstituteVars(fswapped, subst) );
3887  SCIP_CALL( SCIPexprtreeSetVars(fswapped, 2, vars) );
3888  SCIP_CALL( SCIPexprintCompile(exprinterpreter, fswapped) );
3889 
3890  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
3891  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
3892  }
3893 
3894  if( swapped == 0 )
3895  {
3896  /* assume (xval,yval) lie in A1 (lower left triangle) or A2 (upper right triangle) */
3897  SCIP_CALL( generateOrthogonal_lx_ly_Underestimator(scip, exprinterpreter, f, xyref, all_cutcoeff[0], &all_convenvvalue[0], &all_success[0]) );
3898  /* assume (xval,yval) lie in A3 */
3899  SCIP_CALL( generateUnderestimatorParallelYFacets(scip, exprinterpreter, f, xyref, all_cutcoeff[1], &all_convenvvalue[1], &all_success[1]) );
3900  }
3901  else
3902  {
3903  SCIP_Real xyref_[2];
3904 
3905  assert(swapped == 1);
3906 
3907  xyref_[0] = xyref[1];
3908  xyref_[1] = xyref[0];
3909 
3910  /* assume (xval,yval) lie in A1 (lower left triangle) or A2 (upper right triangle) */
3911  SCIP_CALL( generateOrthogonal_lx_ly_Underestimator(scip, exprinterpreter, fswapped, xyref_, all_cutcoeff[0], &all_convenvvalue[0], &all_success[0]) ); /*lint !e644*/
3912  /* assume (xval,yval) lie in A3 */
3913  SCIP_CALL( generateUnderestimatorParallelYFacets(scip, exprinterpreter, fswapped, xyref_, all_cutcoeff[1], &all_convenvvalue[1], &all_success[1]) );
3914 
3915  /* swap back */
3916  swap_buffer = all_cutcoeff[0][0];
3917  all_cutcoeff[0][0] = all_cutcoeff[0][1];
3918  all_cutcoeff[0][1] = swap_buffer;
3919 
3920  swap_buffer = all_cutcoeff[1][0];
3921  all_cutcoeff[1][0] = all_cutcoeff[1][1];
3922  all_cutcoeff[1][1] = swap_buffer;
3923 
3924  SCIP_CALL( SCIPexprtreeFree(&fswapped) );
3925  }
3926 
3927  /* Select the underestimator with the lowest convex envelope */
3928  SCIPdebugMsg(scip, "\n");
3929  SCIPdebugMsg(scip, "Triangulation: convenvvalue=%g\n", all_convenvvalue[0]);
3930  SCIPdebugMsg(scip, "Parallel Y: convenvvalue=%g\n", all_convenvvalue[1]);
3931 
3932  lowest = SCIPinfinity(scip);
3933  lowestidx = -1;
3934 
3935  if( all_success[0] && all_success[1] )
3936  {
3937  *success = TRUE;
3938  for( i = 0; i < 2; ++i )
3939  {
3940  assert(SCIPisFinite(all_cutcoeff[i][0]));
3941  assert(SCIPisFinite(all_cutcoeff[i][1]));
3942  assert(SCIPisFinite(all_cutcoeff[i][2]));
3943  assert(SCIPisFinite(all_cutcoeff[i][3]));
3944 
3945  if( all_convenvvalue[i] < lowest )
3946  {
3947  /* if all_convenvvalue[0] == all_convenvalue[1], take all_convenvvalue[0] */
3948  lowest = all_convenvvalue[i];
3949  lowestidx = i;
3950  }
3951  }
3952  assert(lowestidx >= 0);
3953 
3954  *convenvvalue = all_convenvvalue[lowestidx];
3955  cutcoeff[0] = all_cutcoeff[lowestidx][0];
3956  cutcoeff[1] = all_cutcoeff[lowestidx][1];
3957  cutcoeff[2] = all_cutcoeff[lowestidx][2];
3958  cutcoeff[3] = all_cutcoeff[lowestidx][3];
3959  assert(SCIPisPositive(scip, cutcoeff[2])); /* assert gamma > 0 */
3960  }
3961  else
3962  {
3963  *success = FALSE;
3964  }
3965 
3966  return SCIP_OKAY;
3967 }
3968 
3969 
3970 /** generates a linear underestimator for f(x,y) with f(x,y) being convex in x and convex in y but indefinite
3971  * This is for the case where the cone of the concave directions is (R_+ x R_+) union (R_- x R_-).
3972  * We consider two cases:
3973  * a) the underestimating segmenent connects parallel facets
3974  * b) the underestimating segmenent connects orthogonal facets where
3975  * x=l_x, y=u_y and x=u_x, y=l_y
3976  * We ensure that the parallel facets are the horizontal with y=l_y and y=u_y
3977  * We compute the objective value of the two problems.
3978  * The smaller objective value corresponds to the convex envelope.
3979  * The supporting hyperplane is then constructed at the this point.
3980  * Generates coefficients cutcoeff = (alpha, beta, gamma, delta), such that alpha * x + beta * y - delta <= gamma * f(x,y)
3981  */
3982 static
3984  SCIP* scip, /**< SCIP data structure */
3985  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
3986  SCIP_EXPRTREE* f, /**< function f(x,y) */
3987  SCIP_Real xyref[2], /**< reference values for x and y */
3988  SCIP_Real cutcoeff[4], /**< cut coefficients alpha, beta, gamma, delta */
3989  SCIP_Real* convenvvalue, /**< function value of the convex envelope */
3990  SCIP_Bool* success /**< buffer to store whether coefficients were successfully computed */
3991  )
3992 {
3993  SCIP_VAR* x;
3994  SCIP_VAR* y;
3995  SCIP_Real xlb;
3996  SCIP_Real xub;
3997  SCIP_Real ylb;
3998  SCIP_Real yub;
3999  SCIP_Real xlb_ylb[2];
4000  SCIP_Real xub_yub[2];
4001  SCIP_Real grad_xlb_ylb[2];
4002  SCIP_Real grad_xub_yub[2];
4003  SCIP_Real fval_xlb_ylb;
4004  SCIP_Real fval_xub_yub;
4005 
4006  SCIP_Real all_cutcoeff[2][4];
4007  SCIP_Real all_convenvvalue[2];
4008  SCIP_Bool all_success[2];
4009 
4010  SCIP_Real lowest;
4011  int lowestidx;
4012  int i;
4013 
4014  SCIP_EXPRTREE* fswapped;
4015  SCIP_VAR* vars[2];
4016  SCIP_Bool swapped;
4017  SCIP_Real swap_buffer;
4018  SCIP_EXPR* subst[2];
4019 
4020  assert(scip != NULL);
4021  assert(exprinterpreter != NULL);
4022  assert(f != NULL);
4023  assert(convenvvalue != NULL);
4024  assert(success != NULL);
4025 
4026  x = SCIPexprtreeGetVars(f)[0];
4027  y = SCIPexprtreeGetVars(f)[1];
4028 
4029  xlb = SCIPvarGetLbLocal(x);
4030  xub = SCIPvarGetUbLocal(x);
4031 
4032  ylb = SCIPvarGetLbLocal(y);
4033  yub = SCIPvarGetUbLocal(y);
4034 
4035  *success = FALSE;
4036 
4037  SCIPdebugMsg(scip, "f(%s, %s) = ", SCIPvarGetName(x), SCIPvarGetName(y));
4039  SCIPdebugMsgPrint(scip, "\n");
4040 
4041  xlb_ylb[0] = xlb;
4042  xlb_ylb[1] = ylb;
4043  xub_yub[0] = xub;
4044  xub_yub[1] = yub;
4045 
4046  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xlb_ylb, TRUE, &fval_xlb_ylb, grad_xlb_ylb) );
4047  SCIP_CALL( SCIPexprintGrad(exprinterpreter, f, xub_yub, TRUE, &fval_xub_yub, grad_xub_yub) );
4048 
4049  if( !SCIPisFinite(fval_xlb_ylb) || SCIPisInfinity(scip, REALABS(fval_xlb_ylb)) || !SCIPisFinite(fval_xub_yub) || SCIPisInfinity(scip, REALABS(fval_xub_yub)) )
4050  {
4051  SCIPdebugMsg(scip, "skip 1-convex underestimator since function cannot be evaluated\n");
4052  return SCIP_OKAY;
4053  }
4054 
4055  if( !SCIPisFinite(grad_xlb_ylb[1]) || !SCIPisFinite(grad_xub_yub[0]) )
4056  {
4057  SCIPdebugMsg(scip, "skip 1-convex underestimator since function cannot be differentiated\n");
4058  return SCIP_OKAY;
4059  }
4060 
4061  SCIPdebugMsg(scip, "xval=%g in [%g,%g], yval=%g in [%g,%g]\n",xyref[0],xlb,xub,xyref[1],ylb,yub);
4062 
4063  /* assure f_y(xlb,ylb)*(yub-ylb)* + f(xlb,ylb) >= f_x(xub,yub)*(xlb-xub) + f(xub,yub) */
4064  if( SCIPisGE( scip, fval_xlb_ylb+(yub-ylb)*grad_xlb_ylb[1], fval_xub_yub+(xlb-xub)*grad_xub_yub[0] ) )
4065  {
4066  swapped = 0;
4067  }
4068  else
4069  {
4070  /* swap the variables */
4071  swapped = 1;
4072 
4073  vars[0] = SCIPexprtreeGetVars(f)[1];
4074  vars[1] = SCIPexprtreeGetVars(f)[0];
4075 
4076  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[0], SCIP_EXPR_VARIDX, 1) );
4077  SCIP_CALL( SCIPexprCreate(SCIPblkmem(scip), &subst[1], SCIP_EXPR_VARIDX, 0) );
4078 
4079  SCIP_CALL( SCIPexprtreeCopy(SCIPblkmem(scip), &fswapped, f) );
4080  SCIP_CALL( SCIPexprtreeSubstituteVars(fswapped, subst) );
4081  SCIP_CALL( SCIPexprtreeSetVars(fswapped, 2, vars) );
4082  SCIP_CALL( SCIPexprintCompile(exprinterpreter, fswapped) );
4083 
4084  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[0]);
4085  SCIPexprFreeDeep(SCIPblkmem(scip), &subst[1]);
4086  }
4087 
4088  if( swapped == 0 )
4089  {
4090  /* assume (xval,yval) lie in A1 (lower left triangle) or A2 (upper right triangle) */
4091  SCIP_CALL( generateOrthogonal_lx_uy_Underestimator(scip, exprinterpreter, f, xyref, all_cutcoeff[0], &all_convenvvalue[0], &all_success[0]) );
4092  /* assume (xval,yval) lie in A3*/
4093  SCIP_CALL( generateUnderestimatorParallelYFacets(scip, exprinterpreter, f, xyref, all_cutcoeff[1], &all_convenvvalue[1], &all_success[1]) );
4094  }
4095  else
4096  {
4097  SCIP_Real xyref_[2];
4098 
4099  assert(swapped == 1);
4100 
4101  xyref_[0] = xyref[1];
4102  xyref_[1] = xyref[0];
4103  /* assume (xval,yval) lie in A1 (upper left triangle) or A2 (lower left triangle) */
4104  SCIP_CALL( generateOrthogonal_lx_uy_Underestimator(scip, exprinterpreter, fswapped, xyref_, all_cutcoeff[0], &all_convenvvalue[0], &all_success[0]) ); /*lint !e644*/
4105  /* assume (xval,yval) lie in A3 */
4106  SCIP_CALL( generateUnderestimatorParallelYFacets(scip, exprinterpreter, fswapped, xyref_, all_cutcoeff[1], &all_convenvvalue[1], &all_success[1]) );
4107 
4108  /* swap back */
4109  swap_buffer = all_cutcoeff[0][0];
4110  all_cutcoeff[0][0] = all_cutcoeff[0][1];
4111  all_cutcoeff[0][1] = swap_buffer;
4112 
4113  swap_buffer = all_cutcoeff[1][0];
4114  all_cutcoeff[1][0] = all_cutcoeff[1][1];
4115  all_cutcoeff[1][1] = swap_buffer;
4116 
4117  SCIP_CALL( SCIPexprtreeFree(&fswapped) );
4118  }
4119 
4120  /* select the underestimator with the lowest convex envelope */
4121  SCIPdebugMsg(scip, "\n");
4122  SCIPdebugMsg(scip, "Triangulation: convenvvalue=%g\n", all_convenvvalue[0]);
4123  SCIPdebugMsg(scip, "Parallel Y: convenvvalue=%g\n", all_convenvvalue[1]);
4124 
4125  lowest = SCIPinfinity(scip);
4126  lowestidx = -1;
4127 
4128  if( all_success[0] && all_success[1] )
4129  {
4130  *success = TRUE;
4131  for( i = 0; i < 2; ++i )
4132  {
4133  assert(SCIPisFinite(all_cutcoeff[i][0]));
4134  assert(SCIPisFinite(all_cutcoeff[i][1]));
4135  assert(SCIPisFinite(all_cutcoeff[i][2]));
4136  assert(SCIPisFinite(all_cutcoeff[i][3]));
4137 
4138  /* if all_convenvvalue[0]==all_convenvalue[1], take all_convenvvalue[0] */
4139  if( all_convenvvalue[i] < lowest )
4140  {
4141  lowest = all_convenvvalue[i];
4142  lowestidx = i;
4143  }
4144  }
4145  assert(lowestidx >= 0);
4146 
4147  *convenvvalue = all_convenvvalue[lowestidx];
4148  cutcoeff[0] = all_cutcoeff[lowestidx][0];
4149  cutcoeff[1] = all_cutcoeff[lowestidx][1];
4150  cutcoeff[2] = all_cutcoeff[lowestidx][2];
4151  cutcoeff[3] = all_cutcoeff[lowestidx][3];
4152  assert(SCIPisPositive(scip, cutcoeff[2])); /* assert gamma > 0 */
4153  }
4154  else
4155  {
4156  *success = FALSE;
4157  }
4158 
4159  return SCIP_OKAY;
4160 }
4161 
4162 
4163 /** generates a linear underestimator for f(x,y) with f(x,y) being convex in x and convex in y but indefinite
4164  * generate coefficients cutcoeff = (alpha, beta, gamma, delta), such that alpha * x + beta * y - delta <= gamma * f(x,y)
4165  * 1. If the point lies on the boundary we apply the lifting technique.
4166  * 2. If the point lies in the interior we check the pattern of
4167  * the concave directions and compute the corresponding underestimators.
4168  */
4169 static
4171  SCIP* scip, /**< SCIP data structure */
4172  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
4173  SCIP_CONS* cons, /**< constraint */
4174  SCIP_Real* xyref, /**< reference values for x and y */
4175  SCIP_ROW** row /**< storage for cut */
4176  )
4177 {
4178  SCIP_CONSDATA* consdata;
4179  SCIP_EXPRTREE* f;
4180  SCIP_Real cutcoeff[4];
4181  SCIP_Bool success;
4182  SCIP_Real rhs;
4183  SCIP_Real convenvvalue;
4184 
4185  SCIP_VAR* x;
4186  SCIP_VAR* y;
4187  SCIP_Real xlb;
4188  SCIP_Real xub;
4189  SCIP_Real ylb;
4190  SCIP_Real yub;
4191  SCIP_Real xy_mid[2];
4192  SCIP_Real fval_mid;
4193  SCIP_Real hess[4];
4194 
4195  assert(scip != NULL);
4196  assert(cons != NULL);
4197  assert(row != NULL);
4198 
4199  consdata = SCIPconsGetData(cons);
4200  assert(consdata != NULL);
4201 
4202  assert(consdata->convextype == SCIP_BIVAR_1CONVEX_INDEFINITE);
4203 
4204  assert(!SCIPisInfinity(scip, consdata->rhs));
4205 
4206  f = consdata->f;
4207 
4208  x = SCIPexprtreeGetVars(f)[0];
4209  y = SCIPexprtreeGetVars(f)[1];
4210 
4211  xlb = SCIPvarGetLbLocal(x);
4212  xub = SCIPvarGetUbLocal(x);
4213 
4214  ylb = SCIPvarGetLbLocal(y);
4215  yub = SCIPvarGetUbLocal(y);
4216 
4217  xy_mid[0] = 0.5 * (xlb+xub);
4218  xy_mid[1] = 0.5 * (ylb+yub);
4219 
4220  /* assert that the bounds are finite */
4221  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) || SCIPisInfinity(scip, -ylb) || SCIPisInfinity(scip, yub) )
4222  {
4223  SCIPdebugMsg(scip, "skip underestimate for 1-convex indefinite constraint <%s> since <%s> or <%s> is unbounded\n", SCIPconsGetName(cons), SCIPvarGetName(x), SCIPvarGetName(y));
4224  return SCIP_OKAY;
4225  }
4226 
4227  success = FALSE;
4228  cutcoeff[0] = SCIP_INVALID;
4229  cutcoeff[1] = SCIP_INVALID;
4230  cutcoeff[2] = SCIP_INVALID;
4231  cutcoeff[3] = SCIP_INVALID;
4232 
4233  /* (xval,yval) lie on a boundary */
4234  if( SCIPisFeasEQ(scip,xyref[0],xlb) || SCIPisFeasEQ(scip,xyref[0],xub) || SCIPisFeasEQ(scip,xyref[1],ylb) || SCIPisFeasEQ(scip,xyref[1],yub) )
4235  {
4236  SCIP_CALL( generate1ConvexIndefiniteUnderestimatorAtBoundary(scip, exprinterpreter, f, xyref, cutcoeff, &convenvvalue, &success) );
4237 
4238  if( !success )
4239  {
4240  /* maybe f is not differentiable on boundary, so move reference point into interior
4241  * we do this here w.r.t. both coordinates
4242  */
4243  perturb(&xyref[0], xlb, xub, 0.001);
4244  perturb(&xyref[1], ylb, yub, 0.001);
4245  }
4246  }
4247 
4248  if( !success )
4249  {
4250  /* xyref lies in the interior */
4251  /* check the pattern of the concave directions */
4252  SCIP_CALL( SCIPexprintHessianDense(exprinterpreter, f, xy_mid, TRUE, &fval_mid, hess) );
4253  assert(SCIPisFinite(hess[1]));
4254 
4255  if( hess[1] > 0.0 )
4256  {
4257  /* Pattern A: (R>=0 x R<=0) union (R<=0 x R>=0)*/
4258  SCIPdebugMsg(scip, "Pattern A\n");
4259  SCIP_CALL( generate1ConvexIndefiniteUnderestimatorInTheInteriorPatternA(scip, exprinterpreter, f, xyref, cutcoeff, &convenvvalue, &success) );
4260  }
4261  else
4262  {
4263  /* Pattern B: (R>=0 x R>=0) union (R<=0 x R <=0)*/
4264  SCIPdebugMsg(scip, "Pattern B\n");
4265  SCIP_CALL( generate1ConvexIndefiniteUnderestimatorInTheInteriorPatternB(scip, exprinterpreter, f, xyref, cutcoeff, &convenvvalue, &success) );
4266  }
4267  }
4268 
4269  if( !success )
4270  {
4271  /* bad luck */
4272  *row = NULL;
4273  return SCIP_OKAY;
4274  }
4275 
4276  /* construct row from cut coefficients (alpha, beta, gamma, delta)
4277  * coefficients are such that alpha * x + beta * y - gamma * f(x,y) <= delta,
4278  * i.e., alpha/gamma * x + beta/gamma * y - delta/gamma <= f(x,y)
4279  * -> alpha/gamma * x + beta/gamma * y - delta/gamma + c*z <= f(x,y) + c*z <= rhs
4280  */
4281 
4282  assert(cutcoeff[0] != SCIP_INVALID); /*lint !e777*/
4283  assert(cutcoeff[1] != SCIP_INVALID); /*lint !e777*/
4284  assert(cutcoeff[2] != SCIP_INVALID); /*lint !e777*/
4285  assert(cutcoeff[3] != SCIP_INVALID); /*lint !e777*/
4286  assert(SCIPisFinite(cutcoeff[0]));
4287  assert(SCIPisFinite(cutcoeff[1]));
4288  assert(SCIPisFinite(cutcoeff[2]));
4289  assert(SCIPisFinite(cutcoeff[3]));
4290  assert(SCIPisPositive(scip, cutcoeff[2])); /* assert gamma > 0 */
4291 
4292  if( SCIPisInfinity(scip, REALABS(cutcoeff[0]/cutcoeff[2])) ||
4293  SCIPisInfinity( scip, REALABS(cutcoeff[1]/cutcoeff[2])) ||
4294  SCIPisInfinity( scip, REALABS(cutcoeff[3]/cutcoeff[2])) )
4295  {
4296  *row = NULL;
4297  return SCIP_OKAY;
4298  }
4299 
4300  rhs = consdata->rhs + cutcoeff[3]/cutcoeff[2];
4301  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, SCIPconsGetHdlr(cons), "1ConvexUnderest", -SCIPinfinity(scip), rhs,
4302  TRUE, FALSE /* modifiable */, TRUE /* removable */) );
4303  SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPexprtreeGetVars(consdata->f)[0], cutcoeff[0] / cutcoeff[2]) );
4304  SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPexprtreeGetVars(consdata->f)[1], cutcoeff[1] / cutcoeff[2]) );
4305  if( consdata->z != NULL )
4306  {
4307  SCIP_CALL( SCIPaddVarToRow(scip, *row, consdata->z, consdata->zcoef) );
4308  }
4309 
4310  return SCIP_OKAY;
4311 }
4312 
4313 /** generates a cut */
4314 static
4316  SCIP* scip, /**< SCIP data structure */
4317  SCIP_EXPRINT* exprinterpreter, /**< expressions interpreter */
4318  SCIP_CONS* cons, /**< constraint */
4319  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
4320  SCIP_SIDETYPE violside, /**< for which side of constraint we want to generate a cut */
4321  SCIP_Real cutmaxrange, /**< bound on cut coef range */
4322  SCIP_ROW** row /**< storage for cut */
4323  )
4324 {
4325  SCIP_CONSDATA* consdata;
4326  SCIP_VAR* x;
4327  SCIP_VAR* y;
4328  SCIP_Real x0y0[2];
4329 
4330  assert(scip != NULL);
4331  assert(cons != NULL);
4332  assert(row != NULL);
4333 
4334  consdata = SCIPconsGetData(cons);
4335  assert(consdata != NULL);
4336 
4337  *row = NULL;
4338 
4339  x = SCIPexprtreeGetVars(consdata->f)[0];
4340  y = SCIPexprtreeGetVars(consdata->f)[1];
4341 
4342  x0y0[0] = SCIPgetSolVal(scip, sol, x);
4343  x0y0[1] = SCIPgetSolVal(scip, sol, y);
4344 
4345  assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(x), x0y0[0]));
4346  assert(SCIPisFeasGE(scip, SCIPvarGetUbLocal(x), x0y0[0]));
4347  assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(y), x0y0[1]));
4348  assert(SCIPisFeasGE(scip, SCIPvarGetUbLocal(y), x0y0[1]));
4349 
4350  /* project into box */
4351  x0y0[0] = MIN(MAX(SCIPvarGetLbLocal(x),x0y0[0]),SCIPvarGetUbLocal(x)); /*lint !e666*/
4352  x0y0[1] = MIN(MAX(SCIPvarGetLbLocal(y),x0y0[1]),SCIPvarGetUbLocal(y)); /*lint !e666*/
4353 
4354  SCIPdebugMsgPrint(scip, "\n");
4355  SCIPdebugMsg(scip, "generate cut for constraint <%s> with %s hand side violated by %g\n", SCIPconsGetName(cons), violside == SCIP_SIDETYPE_LEFT ? "left" : "right", violside == SCIP_SIDETYPE_LEFT ? consdata->lhsviol : consdata->rhsviol);
4356  SCIPdebugMsg(scip, "convextype = %d\n",consdata->convextype);
4357  SCIPdebugMsg(scip, "%s = %g with bounds [%g, %g], %s = %g with bounds [%g, %g]",
4360  if( consdata->z != NULL )
4361  SCIPdebugMsgPrint(scip, ", %s = %g with bounds [%g, %g]", SCIPvarGetName(consdata->z), SCIPgetSolVal(scip, sol, consdata->z), SCIPvarGetLbLocal(consdata->z), SCIPvarGetUbLocal(consdata->z));
4362  SCIPdebugMsgPrint(scip, "\n");
4363  SCIPdebugPrintCons(scip, cons, NULL);
4364  SCIPdebugMsgPrint(scip, "\n");
4365 
4366  switch( consdata->convextype )
4367  {
4368  case SCIP_BIVAR_ALLCONVEX:
4369  {
4370  if( violside == SCIP_SIDETYPE_RIGHT )
4371  {
4372  /* rhs is violated */
4373  SCIP_CALL( generateLinearizationCut(scip, exprinterpreter, cons, x0y0, FALSE, row) );
4374  }
4375  else
4376  {
4377  /* lhs is violated */
4378  SCIP_CALL( generateOverestimatingHyperplaneCut(scip, exprinterpreter, cons, x0y0, row) );
4379  }
4380 
4381  break;
4382  }
4383 
4385  {
4386  SCIP_CALL( generateConvexConcaveEstimator(scip, exprinterpreter, cons, x0y0, violside, row) );
4387  break;
4388  }
4389 
4391  {
4392  if( violside == SCIP_SIDETYPE_RIGHT )
4393  {
4394  /* rhs is violated */
4395  SCIP_CALL( generate1ConvexIndefiniteUnderestimator(scip, exprinterpreter, cons, x0y0, row) );
4396  }
4397  else
4398  {
4399  /* lhs is violated */
4400  SCIP_CALL( generateOverestimatingHyperplaneCut(scip, exprinterpreter, cons, x0y0, row) );
4401  }
4402  break;
4403  }
4404  default:
4405  {
4406  SCIPdebugMsg(scip, "cut generation for convexity type not implemented\n");
4407  }
4408  } /*lint !e788*/
4409 
4410  if( *row == NULL )
4411  return SCIP_OKAY;
4412 
4413  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) );
4414 
4415  /* check numerics */
4416  {
4417  SCIP_Real mincoef;
4418  SCIP_Real maxcoef;
4419 
4420  mincoef = SCIPgetRowMinCoef(scip, *row);
4421  maxcoef = SCIPgetRowMaxCoef(scip, *row);
4422 
4423  while( maxcoef / mincoef > cutmaxrange )
4424  {
4425  SCIP_VAR* var;
4426  SCIP_Real coef;
4427  SCIP_Real constant;
4428  int j;
4429 
4430  /* if range of coefficients is bad, find very small coefficients and make them zero */
4431  SCIPdebugMsg(scip, "cut coefficients for constraint <%s> have very large range: mincoef = %g maxcoef = %g\n", SCIPconsGetName(cons), mincoef, maxcoef);
4432 
4433  /* if minimal coefficient is given by z, then give up (probably the maximal coefficient is the problem) */
4434  if( mincoef == consdata->zcoef ) /*lint !e777*/
4435  {
4436  SCIPdebugMsg(scip, "could not eliminate small coefficient, since it comes from linear part\n");
4437  break;
4438  }
4439 
4440  constant = 0.0;
4441  for( j = 0; j < SCIProwGetNNonz(*row); ++j )
4442  {
4443  coef = SCIProwGetVals(*row)[j];
4444  if( !SCIPisEQ(scip, REALABS(coef), mincoef) )
4445  continue;
4446 
4447  var = SCIPcolGetVar(SCIProwGetCols(*row)[j]);
4448  assert(var != NULL);
4449 
4450  /* try to eliminate coefficient with minimal absolute value by weakening cut and try again */
4451  if( ((coef > 0.0 && violside == SCIP_SIDETYPE_RIGHT) || (coef < 0.0 && violside == SCIP_SIDETYPE_LEFT)) && !SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) )
4452  {
4453  SCIPdebugMsg(scip, "eliminate coefficient %g for <%s> = %g [%g, %g]\n", coef, SCIPvarGetName(var), SCIPgetSolVal(scip, sol, var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
4454 
4455  constant += coef * (SCIProwIsLocal(*row) ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var));
4456  SCIP_CALL( SCIPaddVarToRow(scip, *row, var, -coef) );
4457  continue;
4458  }
4459 
4460  if( ((coef < 0.0 && violside == SCIP_SIDETYPE_RIGHT) || (coef > 0.0 && violside == SCIP_SIDETYPE_LEFT)) && !SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) )
4461  {
4462  SCIPdebugMsg(scip, "eliminate coefficient %g for <%s> = %g [%g, %g]\n", coef, SCIPvarGetName(var), SCIPgetSolVal(scip, sol, var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
4463 
4464  constant += coef * (SCIProwIsLocal(*row) ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var));
4465  SCIP_CALL( SCIPaddVarToRow(scip, *row, var, -coef) );
4466  continue;
4467  }
4468 
4469  break;
4470  }
4471 
4472  if( j < SCIProwGetNNonz(*row) )
4473  {
4474  SCIPdebugMsg(scip, "could not eliminate small coefficient\n");
4475  SCIP_CALL( SCIPreleaseRow(scip, row) );
4476  break;
4477  }
4478 
4479  if( violside == SCIP_SIDETYPE_LEFT )
4480  {
4481  SCIP_CALL( SCIPchgRowLhs(scip, *row, SCIProwGetLhs(*row) - constant) );
4482  }
4483  else
4484  {
4485  SCIP_CALL( SCIPchgRowRhs(scip, *row, SCIProwGetRhs(*row) - constant) );
4486  }
4487 
4488  /* update min/max coefficient */
4489  mincoef = SCIPgetRowMinCoef(scip, *row);
4490  maxcoef = SCIPgetRowMaxCoef(scip, *row);
4491  };
4492 
4493  /* avoid numerically very bad cuts */
4494  if( maxcoef / mincoef > cutmaxrange )
4495  {
4496  SCIPdebugMsg(scip, "drop row for constraint <%s> because range of coefficients is too large: mincoef = %g, maxcoef = %g -> range = %g\n",
4497  SCIPconsGetName(cons), mincoef, maxcoef, maxcoef / mincoef);
4498  }
4499 
4500  if( *row != NULL &&
4501  ( (violside == SCIP_SIDETYPE_LEFT && SCIPisInfinity(scip, -SCIProwGetLhs(*row))) ||
4502  (violside == SCIP_SIDETYPE_RIGHT && SCIPisInfinity(scip, SCIProwGetRhs(*row)))) )
4503  {
4504  SCIPdebugMsg(scip, "drop row for constraint <%s> because of very large side: %g\n", SCIPconsGetName(cons), violside == SCIP_SIDETYPE_LEFT ? -SCIProwGetLhs(*row) : SCIProwGetRhs(*row));
4505  SCIP_CALL( SCIPreleaseRow(scip, row) );
4506  }
4507  }
4508 
4509  return SCIP_OKAY;
4510 }
4511 
4512 /** returns whether one side of a constraint function is convex w.r.t. local bounds
4513  * i.e., if side == RIGHT, then returns whether constraint function is convex w.r.t. local bounds
4514  * and if side == LEFT, then returns whether constraint function is concave w.r.t. local bounds
4515  */
4516 static
4518  SCIP* scip, /**< SCIP data structure */
4519  SCIP_CONS* cons, /**< constraint */
4520  SCIP_SIDETYPE side /**< constraint side to consider */
4521  )
4522 {
4523  SCIP_CONSDATA* consdata;
4524  SCIP_VAR** xy;
4525 
4526  consdata = SCIPconsGetData(cons);
4527  assert(consdata != NULL);
4528  assert(consdata->f != NULL);
4529 
4530  switch( consdata->convextype )
4531  {
4532  case SCIP_BIVAR_ALLCONVEX:
4533  /* always convex w.r.t. right hand side and concave w.r.t. left hand side */
4534  return side == SCIP_SIDETYPE_RIGHT;
4535 
4537  {
4538  /* always not convex w.r.t. left hand side */
4539  if( side == SCIP_SIDETYPE_LEFT )
4540  return FALSE;
4541 
4542  xy = SCIPexprtreeGetVars(consdata->f);
4543  assert(xy != NULL);
4544 
4545  /* convex w.r.t. right hand side if one of the variables is fixed */
4546  return SCIPisEQ(scip, SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0])) ||
4547  SCIPisEQ(scip, SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1]));
4548  }
4549 
4551  {
4552  xy = SCIPexprtreeGetVars(consdata->f);
4553  assert(xy != NULL);
4554 
4555  /* convex w.r.t. right hand side if y is fixed and
4556  * convex w.r.t. left hand side if x is fixed */
4557  return (side == SCIP_SIDETYPE_RIGHT && SCIPisEQ(scip, SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1]))) ||
4558  (side == SCIP_SIDETYPE_LEFT && SCIPisEQ(scip, SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0])));
4559  }
4560 
4561  default:
4562  return FALSE;
4563  } /*lint !e788*/
4564 }
4565 
4566 #ifdef SCIP_DEBUG
4567 static
4568 void printEstimator(
4569  SCIP* scip, /**< SCIP data structure */
4570  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
4571  SCIP_CONS* cons, /**< constraint */
4572  SCIP_SIDETYPE side, /**< violated side of constraint */
4573  SCIP_ROW* row /**< row */
4574  )
4575 {
4576  SCIP_CONSDATA* consdata;
4577  const char* varnames[2] = {"x", "y"};
4578  SCIP_VAR* x;
4579  SCIP_VAR* y;
4580  int i;
4581 
4582  assert(scip != NULL);
4583  assert(cons != NULL);
4584  assert(row != NULL);
4585 
4586  consdata = SCIPconsGetData(cons);
4587  assert(consdata != NULL);
4588  x = SCIPexprtreeGetVars(consdata->f)[0];
4589  y = SCIPexprtreeGetVars(consdata->f)[1];
4590 
4591  SCIPinfoMessage(scip, NULL, "splot [%g:%g] [%g:%g] ", SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), SCIPvarGetLbLocal(y), SCIPvarGetUbLocal(y));
4592  SCIPexprtreePrint(consdata->f, SCIPgetMessagehdlr(scip), NULL, varnames, NULL);
4593  SCIPinfoMessage(scip, NULL, "%+g", side == SCIP_SIDETYPE_LEFT ? consdata->lhs : consdata->rhs);
4594 
4595  SCIPinfoMessage(scip, NULL, ", %g", SCIPisInfinity(scip, SCIProwGetRhs(row)) ? -SCIProwGetLhs(row) : -SCIProwGetRhs(row));
4596  for( i = 0; i < SCIProwGetNNonz(row); ++i )
4597  {
4598  SCIP_VAR* var;
4599 
4600  var = SCIPcolGetVar(SCIProwGetCols(row)[i]);
4601  if( var != x && var != y )
4602  continue;
4603 
4604  SCIPinfoMessage(scip, NULL, "%+g * %s", SCIProwGetVals(row)[i], var == x ? "x" : "y");
4605  }
4606 
4607  SCIPinfoMessage(scip, NULL, ", \"< echo '%g %g %g'\" with circles", SCIPgetSolVal(scip, sol, x), SCIPgetSolVal(scip, sol, y), consdata->activity);
4608 
4609  SCIPinfoMessage(scip, NULL, "\n");
4610 }
4611 #endif
4612 
4613 /** tries to separate solution or LP solution by a linear cut
4614  *
4615  * assumes that constraint violations have been computed
4616  */
4617 static
4619  SCIP* scip, /**< SCIP data structure */
4620  SCIP_CONSHDLR* conshdlr, /**< quadratic constraints handler */
4621  SCIP_CONS** conss, /**< constraints */
4622  int nconss, /**< number of constraints */
4623  int nusefulconss, /**< number of constraints that seem to be useful */
4624  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
4625  SCIP_Real minefficacy, /**< minimal efficacy of a cut if it should be added to the LP */
4626  SCIP_Bool inenforcement, /**< whether we are in constraint enforcement */
4627  SCIP_RESULT* result, /**< result of separation */
4628  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 */
4629  )
4630 {
4631  SCIP_CONSHDLRDATA* conshdlrdata;
4632  SCIP_CONSDATA* consdata;
4633  SCIP_SIDETYPE violside;
4634  SCIP_Real feasibility;
4635  SCIP_Real efficacy;
4636  int c;
4637  SCIP_ROW* row;
4638 
4639  assert(scip != NULL);
4640  assert(conshdlr != NULL);
4641  assert(conss != NULL || nconss == 0);
4642  assert(nusefulconss <= nconss);
4643  assert(result != NULL);
4644 
4645  *result = SCIP_FEASIBLE;
4646 
4647  if( bestefficacy != NULL )
4648  *bestefficacy = 0.0;
4649 
4650  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4651  assert(conshdlrdata != NULL);
4652 
4653  for( c = 0; c < nconss; ++c )
4654  {
4655  assert(conss != NULL);
4656  consdata = SCIPconsGetData(conss[c]);
4657  assert(consdata != NULL);
4658 
4659  if( SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) || SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
4660  {
4661  /* we are not feasible anymore */
4662  if( *result == SCIP_FEASIBLE )
4663  *result = SCIP_DIDNOTFIND;
4664 
4665  violside = SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT;
4666 
4667  /* generate cut */
4668  SCIP_CALL( generateCut(scip, conshdlrdata->exprinterpreter, conss[c], sol, violside, conshdlrdata->cutmaxrange, &row) );
4669  if( row == NULL ) /* failed to generate cut */
4670  continue;
4671 
4672  if( sol == NULL )
4673  feasibility = SCIPgetRowLPFeasibility(scip, row);
4674  else
4675  feasibility = SCIPgetRowSolFeasibility(scip, row, sol);
4676  efficacy = -feasibility;
4677 
4678  SCIPdebug( printEstimator(scip, sol, conss[c], violside, row) );
4679 
4680  /* if cut is strong enough or it's weak but we separate on a convex function and accept weak cuts there, add cut to SCIP */
4681  if( (SCIPisGT(scip, efficacy, minefficacy) ||
4682  (inenforcement && SCIPisGT(scip, efficacy, SCIPfeastol(scip)) && isConvexLocal(scip, conss[c], violside))) &&
4683  SCIPisCutApplicable(scip, row) )
4684  {
4685  SCIP_Bool infeasible;
4686 
4687  /* cut cuts off solution sufficiently */
4688  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
4689  if( infeasible )
4690  {
4691  SCIPdebugMsg(scip, "cut for constraint <%s> is infeasible -> cutoff.\n", SCIPconsGetName(conss[c]));
4692  *result = SCIP_CUTOFF;
4693  }
4694  else
4695  {
4696  SCIPdebugMsg(scip, "added cut with efficacy %g for constraint <%s> violated by %g\n", efficacy, SCIPconsGetName(conss[c]), MAX(consdata->lhsviol, consdata->rhsviol));
4697  *result = SCIP_SEPARATED;
4698  }
4699  if( bestefficacy != NULL && efficacy > *bestefficacy )
4700  *bestefficacy = efficacy;
4701 
4702  /* mark row as not removable from LP for current node, if in enforcement */
4703  if( inenforcement && !conshdlrdata->enfocutsremovable )
4704  SCIPmarkRowNotRemovableLocal(scip, row);
4705  }
4706  else
4707  {
4708  SCIPdebugMsg(scip, "abandon cut since efficacy %g is too small or not applicable\n", efficacy);
4709  }
4710 
4711  SCIP_CALL( SCIPreleaseRow(scip, &row) );
4712  }
4713 
4714  if( *result == SCIP_CUTOFF )
4715  break;
4716 
4717  /* enforce only useful constraints
4718  * others are only checked and enforced if we are still feasible or have not found a separating cut yet
4719  */
4720  if( c >= nusefulconss && *result == SCIP_FEASIBLE )
4721  break;
4722  }
4723 
4724  return SCIP_OKAY;
4725 }
4726 
4727 /** processes the event that a new primal solution has been found adds linearizations of all-convex constraints to the cutpool */
4728 static
4729 SCIP_DECL_EVENTEXEC(processNewSolutionEvent)
4731  SCIP_CONSHDLR* conshdlr;
4732  SCIP_CONSHDLRDATA* conshdlrdata;
4733  SCIP_CONS** conss;
4734  int nconss;
4735  SCIP_CONSDATA* consdata;
4736  int c;
4737  SCIP_SOL* sol;
4738  SCIP_ROW* row;
4739  SCIP_Real x0y0[2];
4740 
4741  assert(scip != NULL);
4742  assert(event != NULL);
4743  assert(eventdata != NULL);
4744  assert(eventhdlr != NULL);
4745 
4746  assert((SCIPeventGetType(event) & SCIP_EVENTTYPE_SOLFOUND) != 0);
4747 
4748  conshdlr = (SCIP_CONSHDLR*)eventdata;
4749 
4750  nconss = SCIPconshdlrGetNConss(conshdlr);
4751 
4752  if( nconss == 0 )
4753  return SCIP_OKAY;
4754 
4755  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4756  assert(conshdlrdata != NULL);
4757 
4758  sol = SCIPeventGetSol(event);
4759  assert(sol != NULL);
4760 
4761  /* we are only interested in solution coming from some heuristic other than trysol, but not from the tree
4762  * the reason for ignoring trysol solutions is that they may come from an NLP solve in sepalp, where we already added linearizations,
4763  * or are from the tree, but postprocessed via proposeFeasibleSolution
4764  */
4765  if( SCIPsolGetHeur(sol) == NULL || SCIPsolGetHeur(sol) == conshdlrdata->trysolheur )
4766  return SCIP_OKAY;
4767 
4768  conss = SCIPconshdlrGetConss(conshdlr);
4769  assert(conss != NULL);
4770 
4771  SCIPdebugMsg(scip, "catched new sol event %" SCIP_EVENTTYPE_FORMAT " from heur <%s>; have %d conss\n", SCIPeventGetType(event), SCIPheurGetName(SCIPsolGetHeur(sol)), nconss);
4772 
4773  row = NULL;
4774 
4775  for( c = 0; c < nconss; ++c )
4776  {
4777  if( SCIPconsIsLocal(conss[c]) )
4778  continue;
4779 
4780  consdata = SCIPconsGetData(conss[c]);
4781  assert(consdata != NULL);
4782 
4783  if( consdata->convextype == SCIP_BIVAR_ALLCONVEX && !SCIPisInfinity(scip, consdata->rhs) )
4784  {
4785  SCIP_CALL( SCIPgetSolVals(scip, sol, 2, SCIPexprtreeGetVars(consdata->f), x0y0) );
4786  SCIP_CALL( generateLinearizationCut(scip, conshdlrdata->exprinterpreter, conss[c], x0y0, TRUE, &row) );
4787  }
4788  else
4789  continue;
4790 
4791  if( row == NULL )
4792  continue;
4793 
4794  assert(!SCIProwIsLocal(row));
4795 
4796  SCIP_CALL( SCIPaddPoolCut(scip, row) );
4797  SCIP_CALL( SCIPreleaseRow(scip, &row) );
4798  }
4799 
4800  return SCIP_OKAY;
4801 }
4802 
4803 /** registers unfixed variables in nonlinear terms of violated constraints as external branching candidates
4804  * We score the variables by their gap between the convex envelope and the bivariate function in the current (x,y).
4805  * This value is given by the constraint violation, since we assume that cuts have been generated which support
4806  * the convex envelope in the LP.
4807  */
4808 static
4810  SCIP* scip, /**< SCIP data structure */
4811  SCIP_CONS** conss, /**< constraints to check */
4812  int nconss, /**< number of constraints to check */
4813  int* nnotify /**< counter for number of notifications performed */
4814  )
4815 {
4816  SCIP_CONSDATA* consdata;
4817  SCIP_VAR** xy;
4818  int c;
4819 
4820  assert(scip != NULL);
4821  assert(conss != NULL || nconss == 0);
4822 
4823  *nnotify = 0;
4824 
4825  for( c = 0; c < nconss; ++c )
4826  {
4827  assert(conss != NULL);
4828  consdata = SCIPconsGetData(conss[c]);
4829  assert(consdata != NULL);
4830  SCIPdebugMsg(scip, "cons <%s> violation: %g %g\n", SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol);
4831 
4832  xy = SCIPexprtreeGetVars(consdata->f);
4833  assert(xy != NULL);
4834 
4835  /* @todo prefer binary before continuous, prefer unbounded before bounded */
4836 
4837  switch( consdata->convextype )
4838  {
4840  {
4841  /* need to branch on the variable in which function is concave (or linear) */
4842  if( !SCIPisFeasZero(scip, consdata->lhsviol) )
4843  {
4844  /* regarding left hand side, we are concave in x and convex in y, so branch on x, if not fixed */
4845  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0])) )
4846  {
4847  SCIPdebugMsg(scip, "register variable x = <%s>[%g,%g] in convex-concave <%s> with violation %g %g\n", SCIPvarGetName(xy[0]), SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0]), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol);
4848  SCIP_CALL( SCIPaddExternBranchCand(scip, xy[0], consdata->lhsviol, SCIP_INVALID) );
4849  ++*nnotify;
4850  }
4851  }
4852  if( !SCIPisFeasZero(scip, consdata->rhsviol) )
4853  {
4854  /* regarding right hand side, we are convex in x and concave in y, so branch on y, if not fixed */
4855  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1])) )
4856  {
4857  SCIPdebugMsg(scip, "register variable y = <%s>[%g,%g] in convex-concave <%s> with violation %g %g\n", SCIPvarGetName(xy[1]), SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1]), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol);
4858  SCIP_CALL( SCIPaddExternBranchCand(scip, xy[1], consdata->lhsviol, SCIP_INVALID) );
4859  ++*nnotify;
4860  }
4861  }
4862  break;
4863  }
4864 
4866  {
4867  if( !SCIPisFeasZero(scip, consdata->rhsviol) )
4868  if( SCIPisEQ(scip, SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0])) || SCIPisEQ(scip, SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1])) )
4869  break;
4870 
4871  /* register both variables, if not fixed */
4872  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0])) )
4873  {
4874  SCIPdebugMsg(scip, "register variable x = <%s>[%g,%g] in 1-convex <%s> with violation %g %g\n", SCIPvarGetName(xy[0]), SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0]), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol);
4875  SCIP_CALL( SCIPaddExternBranchCand(scip, xy[0], consdata->lhsviol, SCIP_INVALID) );
4876  ++*nnotify;
4877  }
4878 
4879  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1])) )
4880  {
4881  SCIPdebugMsg(scip, "register variable y = <%s>[%g,%g] in 1-convex <%s> with violation %g %g\n", SCIPvarGetName(xy[1]), SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1]), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol);
4882  SCIP_CALL( SCIPaddExternBranchCand(scip, xy[1], consdata->lhsviol, SCIP_INVALID) );
4883  ++*nnotify;
4884  }
4885 
4886  break;
4887  }
4888 
4889  case SCIP_BIVAR_ALLCONVEX:
4890  {
4891  if( SCIPisFeasZero(scip, consdata->lhsviol) )
4892  continue;
4893  } /*lint -fallthrough*/
4894 
4895  default:
4896  {
4897  /* register both variables, if not fixed */
4898  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0])) )
4899  {
4900  SCIPdebugMsg(scip, "register variable x = <%s>[%g,%g] in allconvex <%s> with violation %g %g\n", SCIPvarGetName(xy[0]), SCIPvarGetLbLocal(xy[0]), SCIPvarGetUbLocal(xy[0]), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol);
4901  SCIP_CALL( SCIPaddExternBranchCand(scip, xy[0], consdata->lhsviol, SCIP_INVALID) );
4902  ++*nnotify;
4903  }
4904 
4905  if( !SCIPisEQ(scip, SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1])) )
4906  {
4907  SCIPdebugMsg(scip, "register variable y = <%s>[%g,%g] in allconvex <%s> with violation %g %g\n", SCIPvarGetName(xy[1]), SCIPvarGetLbLocal(xy[1]), SCIPvarGetUbLocal(xy[1]), SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol);
4908  SCIP_CALL( SCIPaddExternBranchCand(scip, xy[1], consdata->lhsviol, SCIP_INVALID) );
4909  ++*nnotify;
4910  }
4911  }
4912  } /*lint !e788*/
4913  }
4914 
4915  return SCIP_OKAY;
4916 }
4917 
4918 /** registers a nonlinear variable from a violated constraint as branching candidate that has a large absolute value in the relaxation */
4919 static
4921  SCIP* scip, /**< SCIP data structure */
4922  SCIP_CONS** conss, /**< constraints */
4923  int nconss, /**< number of constraints */
4924  SCIP_SOL* sol, /**< solution to enforce (NULL for the LP solution) */
4925  SCIP_VAR** brvar /**< buffer to store branching variable */
4926  )
4927 {
4928  SCIP_CONSDATA* consdata;
4929  SCIP_VAR* var;
4930  SCIP_Real val;
4931  SCIP_Real brvarval;
4932  int i;
4933  int c;
4934 
4935  assert(scip != NULL);
4936  assert(conss != NULL || nconss == 0);
4937 
4938  *brvar = NULL;
4939  brvarval = -1.0;
4940 
4941  for( c = 0; c < nconss; ++c )
4942  {
4943  assert(conss != NULL);
4944  consdata = SCIPconsGetData(conss[c]);
4945  assert(consdata != NULL);
4946  assert(consdata->f != NULL);
4947 
4948  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
4949  continue;
4950 
4951  for( i = 0; i < 2; ++i )
4952  {
4953  var = SCIPexprtreeGetVars(consdata->f)[i];
4954  /* do not propose fixed variables */
4955  if( SCIPisEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) )
4956  continue;
4957  val = SCIPgetSolVal(scip, sol, var);
4958  if( REALABS(val) > brvarval )
4959  {
4960  brvarval = REALABS(val);
4961  *brvar = var;
4962  }
4963  }
4964  }
4965 
4966  if( *brvar != NULL )
4967  {
4968  SCIP_CALL( SCIPaddExternBranchCand(scip, *brvar, brvarval, SCIP_INVALID) );
4969  }
4970 
4971  return SCIP_OKAY;
4972 }
4973 
4974 /** enforces violated bivariate constraints where both nonlinear variables can be assumed to be fixed
4975  * apply a bound change to the remaining linear variable, or recognizing infeasibility
4976  */
4977 static
4979  SCIP* scip, /**< SCIP data structure */
4980  SCIP_CONS** conss, /**< constraints */
4981  int nconss, /**< number of constraints */
4982  SCIP_Bool* reduceddom, /**< whether a domain has been reduced */
4983  SCIP_Bool* infeasible /**< whether we detected infeasibility */
4984  )
4985 {
4986  SCIP_CONSDATA* consdata;
4987  SCIP_INTERVAL nonlinact;
4988  SCIP_Real lhs;
4989  SCIP_Real rhs;
4990  int c;
4991 
4992  assert(scip != NULL);
4993  assert(conss != NULL || nconss == 0);
4994  assert(reduceddom != NULL);
4995  assert(infeasible != NULL);
4996 
4997  *reduceddom = FALSE;
4998  *infeasible = FALSE;
4999 
5000  for( c = 0; c < nconss; ++c )
5001  {
5002  assert(conss != NULL);
5003  consdata = SCIPconsGetData(conss[c]);
5004  assert(consdata != NULL);
5005 
5006  if( !SCIPisGT(scip, consdata->lhsviol, SCIPfeastol(scip)) && !SCIPisGT(scip, consdata->rhsviol, SCIPfeastol(scip)) )
5007  continue;
5008 
5009  /* get activity for f(x,y) */
5010  SCIP_CALL( SCIPevalExprtreeLocalBounds(scip, consdata->f, SCIPinfinity(scip), &nonlinact) );
5011  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), nonlinact));
5012 
5013  /* if all variables are fixed (at least up to epsilson), then the activity of the nonlinear part should be bounded */
5014  assert(!SCIPisInfinity(scip, -SCIPintervalGetInf(nonlinact)));
5015  assert(!SCIPisInfinity(scip, SCIPintervalGetSup(nonlinact)));
5016 
5017  if( !SCIPisInfinity(scip, -consdata->lhs) )
5018  lhs = consdata->lhs - SCIPintervalGetSup(nonlinact);
5019  else
5020  lhs = -SCIPinfinity(scip);
5021 
5022  if( !SCIPisInfinity(scip, consdata->rhs) )
5023  rhs = consdata->rhs - SCIPintervalGetInf(nonlinact);
5024  else
5025  rhs = SCIPinfinity(scip);
5026 
5027  if( consdata->z != NULL )
5028  {
5029  SCIP_Bool tightened;
5030  SCIP_Real coef;
5031 
5032  coef = consdata->zcoef;
5033  assert(!SCIPisZero(scip, coef));
5034 
5035  SCIPdebugMsg(scip, "Linear constraint with one variable: %g <= %g <%s> <= %g\n", lhs, coef, SCIPvarGetName(consdata->z), rhs);
5036 
5037  /* possibly correct lhs/rhs */
5038  if( coef >= 0.0 )
5039  {
5040  if( !SCIPisInfinity(scip, -lhs) )
5041  lhs /= coef;
5042  if( !SCIPisInfinity(scip, rhs) )
5043  rhs /= coef;
5044  }
5045  else
5046  {
5047  SCIP_Real h;
5048  h = rhs;
5049  if( !SCIPisInfinity(scip, -lhs) )
5050  rhs = lhs/coef;
5051  else
5052  rhs = SCIPinfinity(scip);
5053 
5054  if( !SCIPisInfinity(scip, h) )
5055  lhs = h/coef;
5056  else
5057  lhs = -SCIPinfinity(scip);
5058  }
5059  SCIPdebugMsg(scip, "Linear constraint is a bound: %g <= <%s> <= %g\n", lhs, SCIPvarGetName(consdata->z), rhs);
5060 
5061  if( !SCIPisInfinity(scip, -lhs) )
5062  {
5063  SCIP_CALL( SCIPtightenVarLb(scip, consdata->z, lhs, TRUE, infeasible, &tightened) );
5064  if( *infeasible )
5065  {
5066  SCIPdebugMsg(scip, "Lower bound leads to infeasibility.\n");
5067  return SCIP_OKAY;
5068  }
5069  if( tightened )
5070  {
5071  SCIPdebugMsg(scip, "Lower bound changed.\n");
5072  *reduceddom = TRUE;
5073  return SCIP_OKAY;
5074  }
5075  }
5076 
5077  if( !SCIPisInfinity(scip, rhs) )
5078  {
5079  SCIP_CALL( SCIPtightenVarUb(scip, consdata->z, rhs, TRUE, infeasible, &tightened) );
5080  if( *infeasible )
5081  {
5082  SCIPdebugMsg(scip, "Upper bound leads to infeasibility.\n");
5083  return SCIP_OKAY;
5084  }
5085  if( tightened )
5086  {
5087  SCIPdebugMsg(scip, "Upper bound changed.\n");
5088  *reduceddom = TRUE;
5089  return SCIP_OKAY;
5090  }
5091  }
5092  }
5093  else
5094  {
5095  /* no variable, thus check feasibility of lhs <= 0.0 <= rhs */
5096  *infeasible = SCIPisFeasGT(scip, lhs, 0.0) || SCIPisFeasLT(scip, rhs, 0.0);
5097  }
5098  }
5099 
5100  return SCIP_OKAY;
5101 }
5102 
5103 /** tightens bounds on a variable to given interval */
5104 static
5106  SCIP* scip, /**< SCIP data structure */
5107  SCIP_VAR* var, /**< variable which bounds to tighten */
5108  SCIP_INTERVAL bounds, /**< new bounds */
5109  SCIP_CONS* cons, /**< constraint that is propagated */
5110  SCIP_RESULT* result, /**< pointer where to update the result of the propagation call */
5111  int* nchgbds /**< buffer where to add the the number of changed bounds */
5112  )
5113 {
5114  SCIP_Bool infeas;
5115  SCIP_Bool tightened;
5116  SCIP_Real bnd;
5117 
5118  assert(scip != NULL);
5119  assert(var != NULL);
5120  assert(result != NULL);
5121  assert(*result == SCIP_DIDNOTFIND || *result == SCIP_REDUCEDDOM);
5122  assert(nchgbds != NULL);
5123 
5124  if( SCIPintervalIsPositiveInfinity(SCIPinfinity(scip), bounds) ||
5126  SCIPintervalIsEmpty(SCIPinfinity(scip), bounds) )
5127  {
5128  /* domain outside [-infty, +infty] or empty -> declare node infeasible */
5129  SCIPdebugMsg(scip, "found <%s> infeasible due to domain propagation for variable <%s>\n", cons != NULL ? SCIPconsGetName(cons) : "???", SCIPvarGetName(var)); /*lint !e585*/
5130  *result = SCIP_CUTOFF;
5131  return SCIP_OKAY;
5132  }
5133 
5135  {
5136  bnd = SCIPadjustedVarLb(scip, var, SCIPintervalGetInf(bounds));
5137  SCIP_CALL( SCIPtightenVarLb(scip, var, bnd, FALSE, &infeas, &tightened) );
5138  if( infeas )
5139  {
5140  SCIPdebugMsg(scip, "found <%s> infeasible due to domain propagation for variable <%s>\n", cons != NULL ? SCIPconsGetName(cons) : "???", SCIPvarGetName(var)); /*lint !e585*/
5141  *result = SCIP_CUTOFF;
5142  return SCIP_OKAY;
5143  }
5144  if( tightened )
5145  {
5146  SCIPdebugMsg(scip, "tightened lower bound of variable <%s> in constraint <%s> to %g\n", SCIPvarGetName(var), cons != NULL ? SCIPconsGetName(cons) : "???", SCIPvarGetLbLocal(var)); /*lint !e585*/
5147  ++*nchgbds;
5148  *result = SCIP_REDUCEDDOM;
5149  }
5150  }
5151 
5153  {
5154  bnd = SCIPadjustedVarLb(scip, var, SCIPintervalGetSup(bounds));
5155  SCIP_CALL( SCIPtightenVarUb(scip, var, bnd, FALSE, &infeas, &tightened) );
5156  if( infeas )
5157  {
5158  SCIPdebugMsg(scip, "found <%s> infeasible due to domain propagation for variable <%s>\n", cons != NULL ? SCIPconsGetName(cons) : "???", SCIPvarGetName(var)); /*lint !e585*/
5159  *result = SCIP_CUTOFF;
5160  return SCIP_OKAY;
5161  }
5162  if( tightened )
5163  {
5164  SCIPdebugMsg(scip, "tightened upper bound of variable <%s> in constraint <%s> to %g\n", SCIPvarGetName(var), cons != NULL ? SCIPconsGetName(cons) : "???", SCIPvarGetUbLocal(var)); /*lint !e585*/
5165  ++*nchgbds;
5166  *result = SCIP_REDUCEDDOM;
5167  }
5168  }
5169 
5170  return SCIP_OKAY;
5171 }
5172 
5173 /** tightens bounds of z in a single bivariate constraint
5174  * checks for redundancy and infeasibility
5175  */
5176 static
5178  SCIP* scip, /**< SCIP data structure */
5179  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
5180  SCIP_CONS* cons, /**< constraint to process */
5181  SCIP_RESULT* result, /**< pointer to store the result of the propagation call */
5182  int* nchgbds, /**< buffer where to add the the number of changed bounds */
5183  SCIP_Bool* redundant /**< buffer where to store whether constraint has been found to be redundant */
5184  )
5185 {
5186  SCIP_CONSHDLRDATA* conshdlrdata;
5187  SCIP_CONSDATA* consdata;
5188  SCIP_INTERVAL consbounds; /* left and right side of constraint */
5189  SCIP_INTERVAL ftermactivity; /* activity of f(x,y) */
5190  SCIP_INTERVAL ztermactivity; /* activity of c*z */
5191  SCIP_INTERVAL consactivity; /* activity of f(x,y) + c*z */
5192  SCIP_INTERVAL tmp;
5193  SCIP_Bool cutoff;
5194 
5195  assert(scip != NULL);
5196  assert(cons != NULL);
5197  assert(result != NULL);
5198  assert(nchgbds != NULL);
5199 
5200  conshdlrdata = SCIPconshdlrGetData(conshdlr);
5201  assert(conshdlrdata != NULL);
5202  assert(conshdlrdata->exprgraph != NULL);
5203 
5204  consdata = SCIPconsGetData(cons);
5205  assert(consdata != NULL);
5206  assert(consdata->exprgraphnode != NULL);
5207 
5208  *result = SCIP_DIDNOTRUN;
5209  *redundant = FALSE;
5210 
5211  /* extend interval by epsilon to avoid cutoff in forward propagation if constraint is only almost feasible */
5212  SCIPintervalSetBounds(&consbounds,
5213  -infty2infty(SCIPinfinity(scip), INTERVALINFTY, -consdata->lhs+SCIPepsilon(scip)), /*lint !e666*/
5214  +infty2infty(SCIPinfinity(scip), INTERVALINFTY, consdata->rhs+SCIPepsilon(scip)) ); /*lint !e666*/
5215 
5216  /* get activity for f(x,y) */
5217  ftermactivity = SCIPexprgraphGetNodeBounds(consdata->exprgraphnode);
5218  assert(!SCIPintervalIsEmpty(SCIPinfinity(scip), ftermactivity) );
5219 
5220  /* get activity for c*z */
5221  if( consdata->z != NULL )
5222  {
5223  SCIPintervalSetBounds(&ztermactivity,
5224  -infty2infty(SCIPinfinity(scip), INTERVALINFTY, -MIN(SCIPvarGetLbLocal(consdata->z), SCIPvarGetUbLocal(consdata->z))), /*lint !e666*/
5225  +infty2infty(SCIPinfinity(scip), INTERVALINFTY, MAX(SCIPvarGetLbLocal(consdata->z), SCIPvarGetUbLocal(consdata->z)))); /*lint !e666*/
5226  SCIPintervalMulScalar(INTERVALINFTY, &ztermactivity, ztermactivity, consdata->zcoef);
5227  }
5228  else
5229  {
5230  SCIPintervalSet(&ztermactivity, 0.0);
5231  }
5232 
5233  /* get activity for f(x,y)+c*z */
5234  SCIPintervalAdd(INTERVALINFTY, &consactivity, ftermactivity, ztermactivity);
5235 
5236  /* check redundancy */
5237  if( SCIPintervalIsSubsetEQ(INTERVALINFTY, consactivity, consbounds) )
5238  {
5239  SCIPdebugMsg(scip, "found constraint <%s> to be redundant: sides: [%g, %g], activity: [%g, %g]\n",
5240  SCIPconsGetName(cons), consdata->lhs, consdata->rhs, SCIPintervalGetInf(consactivity), SCIPintervalGetSup(consactivity));
5241  *redundant = TRUE;
5242  return SCIP_OKAY;
5243  }
5244 
5245  /* check infeasibility */
5246  if( SCIPintervalAreDisjoint(consbounds, consactivity) )
5247  {
5248  SCIPdebugMsg(scip, "found constraint <%s> to be infeasible; sides: [%g, %g], activity: [%g, %g], infeas: %g\n",
5249  SCIPconsGetName(cons), consdata->lhs, consdata->rhs, SCIPintervalGetInf(consactivity), SCIPintervalGetSup(consactivity),
5250  MAX(consdata->lhs - SCIPintervalGetSup(consactivity), SCIPintervalGetInf(consactivity) - consdata->rhs)); /*lint !e666*/
5251  *result = SCIP_CUTOFF;
5252  return SCIP_OKAY;
5253  }
5254 
5255  /* try to tighten bounds on z */
5256  if( consdata->z != NULL )
5257  {
5258  *result = SCIP_DIDNOTFIND;
5259 
5260  /* compute ([lhs, rhs] - f([xlb,xub], [ylb,yub])) / zcoef */
5261  SCIPintervalSub(INTERVALINFTY, &tmp, consbounds, ftermactivity);
5262  SCIPintervalDivScalar(INTERVALINFTY, &tmp, tmp, consdata->zcoef);
5263 
5264  SCIP_CALL( propagateBoundsTightenVar(scip, consdata->z, tmp, cons, result, nchgbds) );
5265 
5266  if( *result == SCIP_CUTOFF )
5267  return SCIP_OKAY;
5268 
5269  if( *result == SCIP_SUCCESS )
5270  {
5271  SCIPintervalSetBounds(&ztermactivity,
5272  -infty2infty(SCIPinfinity(scip), INTERVALINFTY, -MIN(SCIPvarGetLbLocal(consdata->z), SCIPvarGetUbLocal(consdata->z))), /*lint !e666*/
5273  +infty2infty(SCIPinfinity(scip), INTERVALINFTY, MAX(SCIPvarGetLbLocal(consdata->z), SCIPvarGetUbLocal(consdata->z)))); /*lint !e666*/
5274  SCIPintervalMulScalar(INTERVALINFTY, &ztermactivity, ztermactivity, consdata->zcoef);
5275  }
5276  }
5277