Scippy

SCIP

Solving Constraint Integer Programs

cons_nonlinear.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2019 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_nonlinear.c
17  * @brief constraint handler for nonlinear constraints \f$\textrm{lhs} \leq \sum_{i=1}^n a_ix_i + \sum_{j=1}^m c_jf_j(x) \leq \textrm{rhs}\f$
18  * @author Stefan Vigerske
19  * @author Ingmar Vierhaus (consparse)
20  */
21 
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23 
24 #include "blockmemshell/memory.h"
25 #include <ctype.h>
26 #include "lpi/lpi.h"
27 #include "lpi/type_lpi.h"
28 #include "nlpi/exprinterpret.h"
29 #include "nlpi/nlpi_ipopt.h"
30 #include "nlpi/pub_expr.h"
31 #include "nlpi/type_exprinterpret.h"
32 #include "nlpi/type_nlpi.h"
33 #include "scip/cons_linear.h"
34 #include "scip/cons_nonlinear.h"
35 #define SCIP_PRIVATE_ROWPREP
36 #include "scip/cons_quadratic.h" /* for SCIP_ROWPREP */
37 #include "scip/debug.h"
38 #include "scip/heur_subnlp.h"
39 #include "scip/heur_trysol.h"
40 #include "scip/intervalarith.h"
41 #include "scip/pub_cons.h"
42 #include "scip/pub_event.h"
43 #include "scip/pub_heur.h"
44 #include "scip/pub_lp.h"
45 #include "scip/pub_message.h"
46 #include "scip/pub_misc.h"
47 #include "scip/pub_misc_sort.h"
48 #include "scip/pub_nlp.h"
49 #include "scip/pub_sol.h"
50 #include "scip/pub_tree.h"
51 #include "scip/pub_var.h"
52 #include "scip/scip_branch.h"
53 #include "scip/scip_cons.h"
54 #include "scip/scip_copy.h"
55 #include "scip/scip_cut.h"
56 #include "scip/scip_event.h"
57 #include "scip/scip_expr.h"
58 #include "scip/scip_general.h"
59 #include "scip/scip_heur.h"
60 #include "scip/scip_lp.h"
61 #include "scip/scip_mem.h"
62 #include "scip/scip_message.h"
63 #include "scip/scip_nlp.h"
64 #include "scip/scip_numerics.h"
65 #include "scip/scip_param.h"
66 #include "scip/scip_prob.h"
67 #include "scip/scip_probing.h"
68 #include "scip/scip_sepa.h"
69 #include "scip/scip_sol.h"
70 #include "scip/scip_solvingstats.h"
71 #include "scip/scip_tree.h"
72 #include "scip/scip_var.h"
73 #include <string.h>
74 
75 
76 /* constraint handler properties */
77 #define CONSHDLR_NAME "nonlinear"
78 #define CONSHDLR_DESC "constraint handler for nonlinear constraints"
79 #define CONSHDLR_SEPAPRIORITY 10 /**< priority of the constraint handler for separation */
80 #define CONSHDLR_ENFOPRIORITY -60 /**< priority of the constraint handler for constraint enforcing */
81 #define CONSHDLR_CHECKPRIORITY -4000010 /**< priority of the constraint handler for checking feasibility */
82 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
83 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
84 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
85  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
86 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
87 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
88 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
89 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
90 
91 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP /**< propagation timing mask of the constraint handler */
92 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_ALWAYS /**< presolving timing of the constraint handler (fast, medium, or exhaustive) */
93 
94 #define INTERVALINFTY 1E+43 /**< value for infinity in interval operations */
95 #define BOUNDTIGHTENING_MINSTRENGTH 0.05/**< minimal required bound tightening strength in expression graph domain tightening for propagating bound change */
96 #define INITLPMAXVARVAL 1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */
97 
98 /*
99  * Data structures
100  */
101 
102 /** event data for linear variable bound change events */
103 struct LinVarEventData
104 {
105  SCIP_CONSHDLRDATA* conshdlrdata; /**< the constraint handler data */
106  SCIP_CONS* cons; /**< the constraint */
107  int varidx; /**< the index of the linear variable which bound change is catched */
108  int filterpos; /**< position of eventdata in SCIP's event filter */
109 };
110 typedef struct LinVarEventData LINVAREVENTDATA;
111 
112 /** constraint data for nonlinear constraints */
113 struct SCIP_ConsData
114 {
115  SCIP_Real lhs; /**< left hand side of constraint */
116  SCIP_Real rhs; /**< right hand side of constraint */
117 
118  int nlinvars; /**< number of linear variables */
119  int linvarssize; /**< length of linear variable arrays */
120  SCIP_VAR** linvars; /**< linear variables */
121  SCIP_Real* lincoefs; /**< coefficients of linear variables */
122  LINVAREVENTDATA** lineventdata; /**< eventdata for bound change of linear variable */
123 
124  int nexprtrees; /**< number of expression trees */
125  SCIP_Real* nonlincoefs; /**< coefficients of expression trees */
126  SCIP_EXPRTREE** exprtrees; /**< nonlinear part of constraint */
127  SCIP_EXPRCURV* curvatures; /**< curvature of each expression tree (taking nonlincoefs into account) */
128  SCIP_EXPRGRAPHNODE* exprgraphnode; /**< node in expression graph corresponding to expression tree of this constraint */
129  SCIP_EXPRCURV curvature; /**< curvature of complete nonlinear part, if checked */
130 
131  SCIP_NLROW* nlrow; /**< a nonlinear row representation of this constraint */
132 
133  unsigned int linvarssorted:1; /**< are the linear variables already sorted? */
134  unsigned int linvarsmerged:1; /**< are equal linear variables already merged? */
135 
136  unsigned int iscurvchecked:1; /**< is nonlinear function checked on convexity or concavity ? */
137  unsigned int isremovedfixingslin:1; /**< did we removed fixed/aggr/multiaggr variables in linear part? */
138  unsigned int ispresolved:1; /**< did we checked for possibilities of upgrading or implicit integer variables? */
139  unsigned int forcebackprop:1; /**< should we force to run the backward propagation on our subgraph in the exprgraph? */
140 
141  SCIP_Real minlinactivity; /**< sum of minimal activities of all linear terms with finite minimal activity */
142  SCIP_Real maxlinactivity; /**< sum of maximal activities of all linear terms with finite maximal activity */
143  int minlinactivityinf; /**< number of linear terms with infinite minimal activity */
144  int maxlinactivityinf; /**< number of linear terms with infinity maximal activity */
145  SCIP_Real activity; /**< activity of constraint function w.r.t. current solution */
146  SCIP_Real lhsviol; /**< violation of lower bound by current solution (used temporarily inside constraint handler) */
147  SCIP_Real rhsviol; /**< violation of lower bound by current solution (used temporarily inside constraint handler) */
148 
149  int linvar_maydecrease; /**< index of a variable in linvars that may be decreased without making any other constraint infeasible, or -1 if none */
150  int linvar_mayincrease; /**< index of a variable in linvars that may be increased without making any other constraint infeasible, or -1 if none */
151 
152  SCIP_Real lincoefsmin; /**< maximal absolute value of coefficients in linear part, only available in solving stage */
153  SCIP_Real lincoefsmax; /**< minimal absolute value of coefficients in linear part, only available in solving stage */
154  unsigned int ncuts; /**< number of cuts created for this constraint so far */
155 };
156 
157 /** nonlinear constraint update method */
158 struct SCIP_NlConsUpgrade
159 {
160  SCIP_DECL_NONLINCONSUPGD((*nlconsupgd)); /**< method to call for upgrading nonlinear constraint */
161  SCIP_DECL_EXPRGRAPHNODEREFORM((*nodereform));/**< method to call for reformulating an expression graph node */
162  int priority; /**< priority of upgrading method */
163  SCIP_Bool active; /**< is upgrading enabled */
164 };
165 typedef struct SCIP_NlConsUpgrade SCIP_NLCONSUPGRADE;
166 
167 /** constraint handler data */
168 struct SCIP_ConshdlrData
169 {
170  SCIP_EXPRINT* exprinterpreter; /**< expression interpreter to compute gradients */
171 
172  SCIP_Real cutmaxrange; /**< maximal range (maximal coef / minimal coef) of a cut in order to be added to LP */
173  SCIP_Bool linfeasshift; /**< whether to make solutions in check feasible if possible */
174  SCIP_Bool checkconvexexpensive;/**< whether to apply expensive curvature checking methods */
175  SCIP_Bool assumeconvex; /**< whether functions in inequalities should be assumed to be convex */
176  int maxproprounds; /**< limit on number of propagation rounds for a single constraint within one round of SCIP propagation */
177  SCIP_Bool reformulate; /**< whether to reformulate expression graph */
178  int maxexpansionexponent;/**< maximal exponent where still expanding non-monomial polynomials in expression simplification */
179  SCIP_Real sepanlpmincont; /**< minimal required fraction of continuous variables in problem to use solution of NLP relaxation in root for separation */
180  SCIP_Bool enfocutsremovable; /**< are cuts added during enforcement removable from the LP in the same node? */
181 
182  SCIP_HEUR* subnlpheur; /**< a pointer to the subNLP heuristic, if available */
183  SCIP_HEUR* trysolheur; /**< a pointer to the TRYSOL heuristic, if available */
184  SCIP_EVENTHDLR* linvareventhdlr; /**< our handler for linear variable bound change events */
185  SCIP_EVENTHDLR* nonlinvareventhdlr; /**< our handler for nonlinear variable bound change events */
186  int newsoleventfilterpos;/**< filter position of new solution event handler, if catched */
187 
188  SCIP_NLCONSUPGRADE** nlconsupgrades; /**< nonlinear constraint upgrade methods for specializing nonlinear constraints */
189  int nlconsupgradessize; /**< size of nlconsupgrade array */
190  int nnlconsupgrades; /**< number of nonlinear constraint upgrade methods */
191 
192  SCIP_EXPRGRAPH* exprgraph; /**< expression graph */
193  SCIP* scip; /**< SCIP pointer for use in expression graph callbacks */
194  unsigned int isremovedfixings:1; /**< have fixed variables been removed in the expression graph? */
195  unsigned int ispropagated:1; /**< have current bounds of linear variables in constraints and variables in expression graph been propagated? */
196  unsigned int isreformulated:1; /**< has expression graph been reformulated? */
197  unsigned int sepanlp:1; /**< has a linearization in the NLP relaxation been added? */
198  int naddedreformconss; /**< number of constraints added via reformulation */
199  SCIP_NODE* lastenfonode; /**< the node for which enforcement was called the last time (and some constraint was violated) */
200  int nenforounds; /**< counter on number of enforcement rounds for the current node */
201 };
202 
203 /*
204  * Local methods
205  */
206 
207 /** translate from one value of infinity to another
208  *
209  * if val is >= infty1, then give infty2, else give val
210  */
211 #define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))
212 
213 /* catches variable bound change events on a linear variable in a nonlinear constraint */
214 static
215 SCIP_RETCODE catchLinearVarEvents(
216  SCIP* scip, /**< SCIP data structure */
217  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
218  int linvarpos /**< position of variable in linear variables array */
219  )
220 {
221  SCIP_CONSHDLRDATA* conshdlrdata;
222  SCIP_CONSDATA* consdata;
223  LINVAREVENTDATA* eventdata;
224  SCIP_EVENTTYPE eventtype;
225 
226  assert(scip != NULL);
227  assert(cons != NULL);
228  assert(SCIPconsIsEnabled(cons));
229  assert(SCIPconsIsTransformed(cons));
230 
231  assert(SCIPconsGetHdlr(cons) != NULL);
232  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
233  assert(conshdlrdata != NULL);
234  assert(conshdlrdata->linvareventhdlr != NULL);
235 
236  consdata = SCIPconsGetData(cons);
237  assert(consdata != NULL);
238 
239  assert(linvarpos >= 0);
240  assert(linvarpos < consdata->nlinvars);
241 
242  SCIP_CALL( SCIPallocBlockMemory(scip, &eventdata) );
243 
244  eventtype = SCIP_EVENTTYPE_VARFIXED;
245  if( !SCIPisInfinity(scip, consdata->rhs) )
246  {
247  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest */
248  if( consdata->lincoefs[linvarpos] > 0.0 )
249  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
250  else
251  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
252  }
253  if( !SCIPisInfinity(scip, -consdata->lhs) )
254  {
255  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest */
256  if( consdata->lincoefs[linvarpos] > 0.0 )
257  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
258  else
259  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
260  }
261 
262  eventdata->conshdlrdata = conshdlrdata;
263  eventdata->cons = cons;
264  eventdata->varidx = linvarpos;
265  SCIP_CALL( SCIPcatchVarEvent(scip, consdata->linvars[linvarpos], eventtype, conshdlrdata->linvareventhdlr, (SCIP_EVENTDATA*)eventdata, &eventdata->filterpos) );
266 
267  /* ensure lineventdata array is existing */
268  if( consdata->lineventdata == NULL )
269  {
270  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize) );
271  }
272 
273  consdata->lineventdata[linvarpos] = eventdata;
274 
275  /* since bound changes were not catched before, a possibly stored linear activity may have become outdated, so set to invalid */
276  consdata->minlinactivity = SCIP_INVALID;
277  consdata->maxlinactivity = SCIP_INVALID;
278 
279  /* mark constraint for propagation */
280  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
281 
282  return SCIP_OKAY;
283 }
284 
285 /* drops variable bound change events on a linear variable in a nonlinear constraint */
286 static
287 SCIP_RETCODE dropLinearVarEvents(
288  SCIP* scip, /**< SCIP data structure */
289  SCIP_CONS* cons, /**< constraint for which to catch bound change events */
290  int linvarpos /**< position of variable in linear variables array */
291  )
292 {
293  SCIP_CONSHDLRDATA* conshdlrdata;
294  SCIP_CONSDATA* consdata;
295  SCIP_EVENTTYPE eventtype;
296 
297  assert(scip != NULL);
298  assert(cons != NULL);
299  assert(SCIPconsIsTransformed(cons));
300 
301  assert(SCIPconsGetHdlr(cons) != NULL);
302  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
303  assert(conshdlrdata != NULL);
304  assert(conshdlrdata->linvareventhdlr != NULL);
305 
306  consdata = SCIPconsGetData(cons);
307  assert(consdata != NULL);
308 
309  assert(linvarpos >= 0);
310  assert(linvarpos < consdata->nlinvars);
311  assert(consdata->lineventdata != NULL);
312  assert(consdata->lineventdata[linvarpos] != NULL);
313  assert(consdata->lineventdata[linvarpos]->cons == cons);
314  assert(consdata->lineventdata[linvarpos]->varidx == linvarpos);
315  assert(consdata->lineventdata[linvarpos]->filterpos >= 0);
316 
317  eventtype = SCIP_EVENTTYPE_VARFIXED;
318  if( !SCIPisInfinity(scip, consdata->rhs) )
319  {
320  /* if right hand side is finite, then a tightening in the lower bound of coef*linvar is of interest */
321  if( consdata->lincoefs[linvarpos] > 0.0 )
322  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
323  else
324  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
325  }
326  if( !SCIPisInfinity(scip, -consdata->lhs) )
327  {
328  /* if left hand side is finite, then a tightening in the upper bound of coef*linvar is of interest */
329  if( consdata->lincoefs[linvarpos] > 0.0 )
330  eventtype |= SCIP_EVENTTYPE_UBCHANGED;
331  else
332  eventtype |= SCIP_EVENTTYPE_LBCHANGED;
333  }
334 
335  SCIP_CALL( SCIPdropVarEvent(scip, consdata->linvars[linvarpos], eventtype, conshdlrdata->linvareventhdlr, (SCIP_EVENTDATA*)consdata->lineventdata[linvarpos], consdata->lineventdata[linvarpos]->filterpos) );
336 
337  SCIPfreeBlockMemory(scip, &consdata->lineventdata[linvarpos]); /*lint !e866*/
338 
339  return SCIP_OKAY;
340 }
341 
342 /** locks a linear variable in a constraint */
343 static
344 SCIP_RETCODE lockLinearVariable(
345  SCIP* scip, /**< SCIP data structure */
346  SCIP_CONS* cons, /**< constraint where to lock a variable */
347  SCIP_VAR* var, /**< variable to lock */
348  SCIP_Real coef /**< coefficient of variable in constraint */
349  )
350 {
351  SCIP_CONSDATA* consdata;
352 
353  assert(scip != NULL);
354  assert(cons != NULL);
355  assert(var != NULL);
356  assert(coef != 0.0);
357 
358  consdata = SCIPconsGetData(cons);
359  assert(consdata != NULL);
360 
361  if( coef > 0.0 )
362  {
363  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
364  }
365  else
366  {
367  SCIP_CALL( SCIPlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
368  }
369 
370  return SCIP_OKAY;
371 }
372 
373 /** unlocks a linear variable in a constraint */
374 static
375 SCIP_RETCODE unlockLinearVariable(
376  SCIP* scip, /**< SCIP data structure */
377  SCIP_CONS* cons, /**< constraint where to unlock a variable */
378  SCIP_VAR* var, /**< variable to unlock */
379  SCIP_Real coef /**< coefficient of variable in constraint */
380  )
381 {
382  SCIP_CONSDATA* consdata;
383 
384  assert(scip != NULL);
385  assert(cons != NULL);
386  assert(!SCIPconsIsLockedType(cons, SCIP_LOCKTYPE_CONFLICT));
387  assert(var != NULL);
388  assert(coef != 0.0);
389 
390  consdata = SCIPconsGetData(cons);
391  assert(consdata != NULL);
392 
393  if( coef > 0.0 )
394  {
395  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, -consdata->lhs), !SCIPisInfinity(scip, consdata->rhs)) );
396  }
397  else
398  {
399  SCIP_CALL( SCIPunlockVarCons(scip, var, cons, !SCIPisInfinity(scip, consdata->rhs), !SCIPisInfinity(scip, -consdata->lhs)) );
400  }
401 
402  return SCIP_OKAY;
403 }
404 
405 /** computes the minimal and maximal activity for the linear part in a constraint data
406  * only sums up terms that contribute finite values
407  * gives the number of terms that contribute infinite values
408  * only computes those activities where the corresponding side of the constraint is finite
409  */
410 static
411 void consdataUpdateLinearActivity(
412  SCIP* scip, /**< SCIP data structure */
413  SCIP_CONSDATA* consdata /**< constraint data */
414  )
415 { /*lint --e{666}*/
416  SCIP_ROUNDMODE prevroundmode;
417  int i;
418  SCIP_Real bnd;
419 
420  assert(scip != NULL);
421  assert(consdata != NULL);
422 
423  if( consdata->minlinactivity != SCIP_INVALID && consdata->maxlinactivity != SCIP_INVALID ) /*lint !e777*/
424  {
425  /* activities should be uptodate */
426  assert(consdata->minlinactivityinf >= 0);
427  assert(consdata->maxlinactivityinf >= 0);
428  return;
429  }
430 
431  consdata->minlinactivityinf = 0;
432  consdata->maxlinactivityinf = 0;
433 
434  /* if lhs is -infinite, then we do not compute a maximal activity, so we set it to infinity
435  * if rhs is infinite, then we do not compute a minimal activity, so we set it to -infinity
436  */
437  consdata->minlinactivity = SCIPisInfinity(scip, consdata->rhs) ? -INTERVALINFTY : 0.0;
438  consdata->maxlinactivity = SCIPisInfinity(scip, -consdata->lhs) ? INTERVALINFTY : 0.0;
439 
440  if( consdata->nlinvars == 0 )
441  return;
442 
443  /* if the activities computed here should be still uptodate after boundchanges,
444  * variable events need to be catched */
445  assert(consdata->lineventdata != NULL);
446 
447  prevroundmode = SCIPintervalGetRoundingMode();
448 
449  if( !SCIPisInfinity(scip, consdata->rhs) )
450  {
451  /* compute minimal activity only if there is a finite right hand side */
452  SCIPintervalSetRoundingModeDownwards();
453 
454  for( i = 0; i < consdata->nlinvars; ++i )
455  {
456  assert(SCIPvarGetLbLocal(consdata->linvars[i]) <= SCIPvarGetUbLocal(consdata->linvars[i]));
457  assert(consdata->lineventdata[i] != NULL);
458  if( consdata->lincoefs[i] >= 0.0 )
459  {
460  bnd = SCIPvarGetLbLocal(consdata->linvars[i]);
461  if( SCIPisInfinity(scip, -bnd) )
462  {
463  ++consdata->minlinactivityinf;
464  continue;
465  }
466  assert(!SCIPisInfinity(scip, bnd)); /* do not like variables that are fixed at +infinity */
467  }
468  else
469  {
470  bnd = SCIPvarGetUbLocal(consdata->linvars[i]);
471  if( SCIPisInfinity(scip, bnd) )
472  {
473  ++consdata->minlinactivityinf;
474  continue;
475  }
476  assert(!SCIPisInfinity(scip, -bnd)); /* do not like variables that are fixed at -infinity */
477  }
478  consdata->minlinactivity += consdata->lincoefs[i] * bnd;
479  }
480  }
481 
482  if( !SCIPisInfinity(scip, -consdata->lhs) )
483  {
484  /* compute maximal activity only if there is a finite left hand side */
485  SCIPintervalSetRoundingModeUpwards();
486 
487  for( i = 0; i < consdata->nlinvars; ++i )
488  {
489  assert(consdata->lineventdata[i] != NULL);
490  assert(SCIPvarGetLbLocal(consdata->linvars[i]) <= SCIPvarGetUbLocal(consdata->linvars[i]));
491  if( consdata->lincoefs[i] >= 0.0 )
492  {
493  bnd = SCIPvarGetUbLocal(consdata->linvars[i]);
494  if( SCIPisInfinity(scip, bnd) )
495  {
496  ++consdata->maxlinactivityinf;
497  continue;
498  }
499  assert(!SCIPisInfinity(scip, -bnd)); /* do not like variables that are fixed at -infinity */
500  }
501  else
502  {
503  bnd = SCIPvarGetLbLocal(consdata->linvars[i]);
504  if( SCIPisInfinity(scip, -bnd) )
505  {
506  ++consdata->maxlinactivityinf;
507  continue;
508  }
509  assert(!SCIPisInfinity(scip, bnd)); /* do not like variables that are fixed at +infinity */
510  }
511  consdata->maxlinactivity += consdata->lincoefs[i] * bnd;
512  }
513  }
514  assert(consdata->minlinactivity <= consdata->maxlinactivity || consdata->minlinactivityinf > 0 || consdata->maxlinactivityinf > 0);
515 
516  SCIPintervalSetRoundingMode(prevroundmode);
517 }
518 
519 /** update the linear activities after a change in the lower bound of a variable */
520 static
521 void consdataUpdateLinearActivityLbChange(
522  SCIP* scip, /**< SCIP data structure */
523  SCIP_CONSDATA* consdata, /**< constraint data */
524  SCIP_Real coef, /**< coefficient of variable in constraint */
525  SCIP_Real oldbnd, /**< previous lower bound of variable */
526  SCIP_Real newbnd /**< new lower bound of variable */
527  )
528 {
529  SCIP_ROUNDMODE prevroundmode;
530 
531  assert(scip != NULL);
532  assert(consdata != NULL);
533  /* we can't deal with lower bounds at infinity */
534  assert(!SCIPisInfinity(scip, oldbnd));
535  assert(!SCIPisInfinity(scip, newbnd));
536 
537  /* assume lhs <= a*x + y <= rhs, then the following boundchanges can be deduced:
538  * a > 0: y <= rhs - a*lb(x), y >= lhs - a*ub(x)
539  * a < 0: y <= rhs - a*ub(x), y >= lhs - a*lb(x)
540  */
541 
542  if( coef > 0.0 )
543  {
544  /* we should only be called if rhs is finite */
545  assert(!SCIPisInfinity(scip, consdata->rhs));
546 
547  /* we have no min activities computed so far, so cannot update */
548  if( consdata->minlinactivity == SCIP_INVALID ) /*lint !e777*/
549  return;
550 
551  assert(consdata->minlinactivity > -INTERVALINFTY);
552 
553  prevroundmode = SCIPintervalGetRoundingMode();
554  SCIPintervalSetRoundingModeDownwards();
555 
556  /* update min activity */
557  if( SCIPisInfinity(scip, -oldbnd) )
558  {
559  --consdata->minlinactivityinf;
560  assert(consdata->minlinactivityinf >= 0);
561  }
562  else
563  {
564  consdata->minlinactivity += SCIPintervalNegateReal(coef) * oldbnd;
565  }
566 
567  if( SCIPisInfinity(scip, -newbnd) )
568  {
569  ++consdata->minlinactivityinf;
570  }
571  else
572  {
573  consdata->minlinactivity += coef * newbnd;
574  }
575 
576  SCIPintervalSetRoundingMode(prevroundmode);
577  }
578  else
579  {
580  /* we should only be called if lhs is finite */
581  assert(!SCIPisInfinity(scip, -consdata->lhs));
582 
583  /* we have no max activities computed so far, so cannot update */
584  if( consdata->maxlinactivity == SCIP_INVALID ) /*lint !e777*/
585  return;
586 
587  assert(consdata->maxlinactivity < INTERVALINFTY);
588 
589  prevroundmode = SCIPintervalGetRoundingMode();
590  SCIPintervalSetRoundingModeUpwards();
591 
592  /* update max activity */
593  if( SCIPisInfinity(scip, -oldbnd) )
594  {
595  --consdata->maxlinactivityinf;
596  assert(consdata->maxlinactivityinf >= 0);
597  }
598  else
599  {
600  consdata->maxlinactivity += SCIPintervalNegateReal(coef) * oldbnd;
601  }
602 
603  if( SCIPisInfinity(scip, -newbnd) )
604  {
605  ++consdata->maxlinactivityinf;
606  }
607  else
608  {
609  consdata->maxlinactivity += coef * newbnd;
610  }
611 
612  SCIPintervalSetRoundingMode(prevroundmode);
613  }
614 }
615 
616 /** update the linear activities after a change in the upper bound of a variable */
617 static
618 void consdataUpdateLinearActivityUbChange(
619  SCIP* scip, /**< SCIP data structure */
620  SCIP_CONSDATA* consdata, /**< constraint data */
621  SCIP_Real coef, /**< coefficient of variable in constraint */
622  SCIP_Real oldbnd, /**< previous lower bound of variable */
623  SCIP_Real newbnd /**< new lower bound of variable */
624  )
625 {
626  SCIP_ROUNDMODE prevroundmode;
627 
628  assert(scip != NULL);
629  assert(consdata != NULL);
630  /* we can't deal with upper bounds at -infinity */
631  assert(!SCIPisInfinity(scip, -oldbnd));
632  assert(!SCIPisInfinity(scip, -newbnd));
633 
634  /* assume lhs <= a*x + y <= rhs, then the following boundchanges can be deduced:
635  * a > 0: y <= rhs - a*lb(x), y >= lhs - a*ub(x)
636  * a < 0: y <= rhs - a*ub(x), y >= lhs - a*lb(x)
637  */
638  if( coef > 0.0 )
639  {
640  /* we should only be called if lhs is finite */
641  assert(!SCIPisInfinity(scip, -consdata->lhs));
642 
643  /* we have no max activities computed so far, so cannot update */
644  if( consdata->maxlinactivity == SCIP_INVALID ) /*lint !e777*/
645  return;
646 
647  assert(consdata->maxlinactivity < INTERVALINFTY);
648 
649  prevroundmode = SCIPintervalGetRoundingMode();
650  SCIPintervalSetRoundingModeUpwards();
651 
652  /* update max activity */
653  if( SCIPisInfinity(scip, oldbnd) )
654  {
655  --consdata->maxlinactivityinf;
656  assert(consdata->maxlinactivityinf >= 0);
657  }
658  else
659  {
660  consdata->maxlinactivity += SCIPintervalNegateReal(coef) * oldbnd;
661  }
662 
663  if( SCIPisInfinity(scip, newbnd) )
664  {
665  ++consdata->maxlinactivityinf;
666  }
667  else
668  {
669  consdata->maxlinactivity += coef * newbnd;
670  }
671 
672  SCIPintervalSetRoundingMode(prevroundmode);
673  }
674  else
675  {
676  /* we should only be called if rhs is finite */
677  assert(!SCIPisInfinity(scip, consdata->rhs));
678 
679  /* we have no min activities computed so far, so cannot update */
680  if( consdata->minlinactivity == SCIP_INVALID ) /*lint !e777*/
681  return;
682 
683  assert(consdata->minlinactivity > -INTERVALINFTY);
684 
685  prevroundmode = SCIPintervalGetRoundingMode();
686  SCIPintervalSetRoundingModeDownwards();
687 
688  /* update min activity */
689  if( SCIPisInfinity(scip, oldbnd) )
690  {
691  --consdata->minlinactivityinf;
692  assert(consdata->minlinactivityinf >= 0);
693  }
694  else
695  {
696  consdata->minlinactivity += SCIPintervalNegateReal(coef) * oldbnd;
697  }
698 
699  if( SCIPisInfinity(scip, newbnd) )
700  {
701  ++consdata->minlinactivityinf;
702  }
703  else
704  {
705  consdata->minlinactivity += coef * newbnd;
706  }
707 
708  SCIPintervalSetRoundingMode(prevroundmode);
709  }
710 }
711 
712 /** processes variable fixing or bound change event */
713 static
714 SCIP_DECL_EVENTEXEC(processLinearVarEvent)
715 {
716  SCIP_CONS* cons;
717  SCIP_CONSDATA* consdata;
718  SCIP_EVENTTYPE eventtype;
719  int varidx;
720 
721  assert(scip != NULL);
722  assert(event != NULL);
723  assert(eventdata != NULL);
724  assert(eventhdlr != NULL);
725 
726  cons = ((LINVAREVENTDATA*)eventdata)->cons;
727  assert(cons != NULL);
728 
729  consdata = SCIPconsGetData(cons);
730  assert(consdata != NULL);
731 
732  varidx = ((LINVAREVENTDATA*)eventdata)->varidx;
733  assert(varidx >= 0);
734  assert(varidx < consdata->nlinvars);
735 
736  eventtype = SCIPeventGetType(event);
737 
738  if( eventtype & SCIP_EVENTTYPE_VARFIXED )
739  {
740  consdata->isremovedfixingslin = FALSE;
741  }
742 
743  if( eventtype & SCIP_EVENTTYPE_BOUNDCHANGED )
744  {
745  /* update activity bounds for linear terms */
746  if( eventtype & SCIP_EVENTTYPE_LBCHANGED )
747  consdataUpdateLinearActivityLbChange(scip, consdata, consdata->lincoefs[varidx], SCIPeventGetOldbound(event), SCIPeventGetNewbound(event));
748  else
749  consdataUpdateLinearActivityUbChange(scip, consdata, consdata->lincoefs[varidx], SCIPeventGetOldbound(event), SCIPeventGetNewbound(event));
750 
751  if( eventtype & SCIP_EVENTTYPE_BOUNDTIGHTENED )
752  {
753  assert(((LINVAREVENTDATA*)eventdata)->conshdlrdata != NULL);
754  ((LINVAREVENTDATA*)eventdata)->conshdlrdata->ispropagated = FALSE;
755 
756  /* mark constraint for propagation */
757  SCIP_CALL( SCIPmarkConsPropagate(scip, cons) );
758  }
759  }
760 
761  return SCIP_OKAY;
762 }
763 
764 /** processes bound change events for variables in expression graph */
765 static
766 SCIP_DECL_EVENTEXEC(processNonlinearVarEvent)
767 {
768  SCIP_CONSHDLRDATA* conshdlrdata;
769  SCIP_EVENTTYPE eventtype;
770 
771  assert(scip != NULL);
772  assert(event != NULL);
773  assert(eventdata != NULL);
774  assert(eventhdlr != NULL);
775 
776  conshdlrdata = (SCIP_CONSHDLRDATA*)SCIPeventhdlrGetData(eventhdlr);
777  assert(conshdlrdata != NULL);
778  assert(conshdlrdata->exprgraph != NULL);
779 
780  eventtype = SCIPeventGetType(event);
781  assert( eventtype & (SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED) );
782 
783  if( eventtype & SCIP_EVENTTYPE_BOUNDCHANGED )
784  {
785  SCIP_Real newbd;
786 
787  SCIPdebugMsg(scip, "changed %s bound on expression graph variable <%s> from %g to %g\n",
788  (eventtype & SCIP_EVENTTYPE_LBCHANGED) ? "lower" : "upper",
789  SCIPvarGetName(SCIPeventGetVar(event)), SCIPeventGetOldbound(event), SCIPeventGetNewbound(event));
790 
791  if( eventtype & SCIP_EVENTTYPE_BOUNDTIGHTENED )
792  conshdlrdata->ispropagated = FALSE;
793  /* @todo a global bound tightening may yield in convex/concave curvatures, maybe the iscurvcheck flag should be reset? */
794 
795  /* update variable bound in expression graph, with epsilon added */
796  if( eventtype & SCIP_EVENTTYPE_LBCHANGED )
797  {
798  newbd = -infty2infty(SCIPinfinity(scip), INTERVALINFTY, -SCIPeventGetNewbound(event)); /*lint !e666*/
799  /* if newbd in [0,eps], then relax to 0.0, otherwise relax by -epsilon
800  * this is to ensure that an original positive variable does not get negative by this, which may have an adverse effect on convexity recoginition, for example */
801  if( newbd >= 0.0 && newbd <= SCIPepsilon(scip) )
802  newbd = 0.0;
803  else
804  newbd -= SCIPepsilon(scip);
805  SCIPexprgraphSetVarNodeLb(conshdlrdata->exprgraph, (SCIP_EXPRGRAPHNODE*)eventdata, newbd);
806  }
807  else
808  {
809  newbd = +infty2infty(SCIPinfinity(scip), INTERVALINFTY, SCIPeventGetNewbound(event)); /*lint !e666*/
810  /* if newbd in [-eps,0], then relax to 0.0, otherwise relax by +epsilon */
811  if( newbd >= -SCIPepsilon(scip) && newbd <= 0.0 )
812  newbd = 0.0;
813  else
814  newbd += SCIPepsilon(scip);
815  SCIPexprgraphSetVarNodeUb(conshdlrdata->exprgraph, (SCIP_EXPRGRAPHNODE*)eventdata, newbd);
816  }
817  }
818  else
819  {
820  assert(eventtype & SCIP_EVENTTYPE_VARFIXED);
821  conshdlrdata->isremovedfixings = FALSE;
822  }
823 
824  return SCIP_OKAY;
825 }
826 
827 /** callback method for variable addition in expression graph */
828 static
829 SCIP_DECL_EXPRGRAPHVARADDED( exprgraphVarAdded )
830 {
831  SCIP_CONSHDLRDATA* conshdlrdata;
832  SCIP_INTERVAL varbounds;
833  SCIP_VAR* var_;
834 
835  assert(exprgraph != NULL);
836  assert(var != NULL);
837  assert(varnode != NULL);
838 
839  var_ = (SCIP_VAR*)var;
840 
841  conshdlrdata = (SCIP_CONSHDLRDATA*)userdata;
842  assert(conshdlrdata != NULL);
843  assert(conshdlrdata->exprgraph == exprgraph);
844 
845  /* catch variable bound change events */
846  SCIP_CALL( SCIPcatchVarEvent(conshdlrdata->scip, (SCIP_VAR*)var, SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED, conshdlrdata->nonlinvareventhdlr, (SCIP_EVENTDATA*)varnode, NULL) );
847  SCIPdebugMessage("catch boundchange events on new expression graph variable <%s>\n", SCIPvarGetName(var_));
848 
849  /* set current bounds in expression graph */
850  SCIPintervalSetBounds(&varbounds,
851  -infty2infty(SCIPinfinity(conshdlrdata->scip), INTERVALINFTY, -MIN(SCIPvarGetLbLocal(var_), SCIPvarGetUbLocal(var_))), /*lint !e666*/
852  +infty2infty(SCIPinfinity(conshdlrdata->scip), INTERVALINFTY, MAX(SCIPvarGetLbLocal(var_), SCIPvarGetUbLocal(var_))) /*lint !e666*/
853  );
854  SCIPexprgraphSetVarNodeBounds(exprgraph, varnode, varbounds);
855 
856  SCIP_CALL( SCIPaddVarLocksType(conshdlrdata->scip, var_, SCIP_LOCKTYPE_MODEL, 1, 1) );
857  SCIPdebugMessage("increased up- and downlocks of variable <%s>\n", SCIPvarGetName(var_));
858 
859  SCIP_CALL( SCIPcaptureVar(conshdlrdata->scip, var_) );
860  SCIPdebugMessage("capture variable <%s>\n", SCIPvarGetName(var_));
861 
862  conshdlrdata->isremovedfixings &= SCIPvarIsActive(var_);
863  conshdlrdata->ispropagated = FALSE;
864 
865  return SCIP_OKAY;
866 }
867 
868 /** callback method for variable removal in expression graph */
869 static
870 SCIP_DECL_EXPRGRAPHVARREMOVE( exprgraphVarRemove )
871 {
872  SCIP_CONSHDLRDATA* conshdlrdata;
873  SCIP_VAR* var_;
874 
875  assert(exprgraph != NULL);
876  assert(var != NULL);
877  assert(varnode != NULL);
878 
879  var_ = (SCIP_VAR*)var;
880 
881  conshdlrdata = (SCIP_CONSHDLRDATA*)userdata;
882  assert(conshdlrdata != NULL);
883  assert(conshdlrdata->exprgraph == exprgraph);
884 
885  SCIP_CALL( SCIPdropVarEvent(conshdlrdata->scip, var_, SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_VARFIXED, conshdlrdata->nonlinvareventhdlr, (SCIP_EVENTDATA*)varnode, -1) );
886  SCIPdebugMessage("drop boundchange events on expression graph variable <%s>\n", SCIPvarGetName(var_));
887 
888  SCIP_CALL( SCIPaddVarLocksType(conshdlrdata->scip, var_, SCIP_LOCKTYPE_MODEL, -1, -1) );
889  SCIPdebugMessage("decreased up- and downlocks of variable <%s>\n", SCIPvarGetName(var_));
890 
891  SCIPdebugMessage("release variable <%s>\n", SCIPvarGetName(var_));
892  SCIP_CALL( SCIPreleaseVar(conshdlrdata->scip, &var_) );
893 
894  return SCIP_OKAY;
895 }
896 
897 /* adds expression trees to constraint */
898 static
899 SCIP_RETCODE consdataAddExprtrees(
900  SCIP* scip, /**< SCIP data structure */
901  SCIP_CONSDATA* consdata, /**< nonlinear constraint data */
902  int nexprtrees, /**< number of expression trees */
903  SCIP_EXPRTREE** exprtrees, /**< expression trees */
904  SCIP_Real* coefs, /**< coefficients of expression trees, or NULL if all 1.0 */
905  SCIP_Bool copytrees /**< whether trees should be copied or ownership should be assumed */
906  )
907 {
908  int i;
909 
910  assert(scip != NULL);
911  assert(consdata != NULL);
912  assert(consdata->exprtrees != NULL || consdata->nexprtrees == 0);
913 
914  if( nexprtrees == 0 )
915  return SCIP_OKAY;
916 
917  /* invalidate activity information */
918  consdata->activity = SCIP_INVALID;
919 
920  /* invalidate nonlinear row */
921  if( consdata->nlrow != NULL )
922  {
923  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
924  }
925 
926  consdata->ispresolved = FALSE;
927  consdata->curvature = SCIP_EXPRCURV_UNKNOWN;
928  consdata->iscurvchecked = FALSE;
929 
930  if( consdata->nexprtrees == 0 )
931  {
932  assert(consdata->exprtrees == NULL);
933  assert(consdata->nonlincoefs == NULL);
934  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->exprtrees, nexprtrees) );
935  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->nonlincoefs, nexprtrees) );
936  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->curvatures, nexprtrees) );
937  }
938  else
939  {
940  assert(consdata->exprtrees != NULL);
941  assert(consdata->nonlincoefs != NULL);
942  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->exprtrees, consdata->nexprtrees, consdata->nexprtrees + nexprtrees) );
943  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->nonlincoefs, consdata->nexprtrees, consdata->nexprtrees + nexprtrees) );
944  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->curvatures, consdata->nexprtrees, consdata->nexprtrees + nexprtrees) );
945  }
946 
947  for( i = 0; i < nexprtrees; ++i )
948  {
949  assert(exprtrees[i] != NULL);
950  /* the expression tree need to have SCIP_VAR*'s stored */
951  assert(SCIPexprtreeGetNVars(exprtrees[i]) == 0 || SCIPexprtreeGetVars(exprtrees[i]) != NULL);
952 
953  if( copytrees )
954  {
955  SCIP_CALL( SCIPexprtreeCopy(SCIPblkmem(scip), &consdata->exprtrees[consdata->nexprtrees + i], exprtrees[i]) );
956  }
957  else
958  {
959  consdata->exprtrees[consdata->nexprtrees + i] = exprtrees[i];
960  }
961 
962  consdata->nonlincoefs[consdata->nexprtrees + i] = (coefs != NULL ? coefs[i] : 1.0);
963  consdata->curvatures[consdata->nexprtrees + i] = SCIP_EXPRCURV_UNKNOWN;
964  }
965  consdata->nexprtrees += nexprtrees;
966 
967  return SCIP_OKAY;
968 }
969 
970 /* sets expression trees of constraints, clears previously ones */
971 static
972 SCIP_RETCODE consdataSetExprtrees(
973  SCIP* scip, /**< SCIP data structure */
974  SCIP_CONSDATA* consdata, /**< nonlinear constraint data */
975  int nexprtrees, /**< number of expression trees */
976  SCIP_EXPRTREE** exprtrees, /**< expression trees */
977  SCIP_Real* coefs, /**< coefficients of expression trees, or NULL if all 1.0 */
978  SCIP_Bool copytrees /**< whether trees should be copied or ownership should be assumed */
979  )
980 {
981  int i;
982 
983  assert(scip != NULL);
984  assert(consdata != NULL);
985  assert(consdata->exprtrees != NULL || consdata->nexprtrees == 0);
986 
987  /* clear existing expression trees */
988  if( consdata->nexprtrees > 0 )
989  {
990  for( i = 0; i < consdata->nexprtrees; ++i )
991  {
992  assert(consdata->exprtrees[i] != NULL);
993  SCIP_CALL( SCIPexprtreeFree(&consdata->exprtrees[i]) );
994  }
995 
996  /* invalidate activity information */
997  consdata->activity = SCIP_INVALID;
998 
999  /* invalidate nonlinear row */
1000  if( consdata->nlrow != NULL )
1001  {
1002  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1003  }
1004 
1005  consdata->ispresolved = FALSE;
1006  consdata->curvature = SCIP_EXPRCURV_LINEAR;
1007  consdata->iscurvchecked = TRUE;
1008 
1009  SCIPfreeBlockMemoryArray(scip, &consdata->exprtrees, consdata->nexprtrees);
1010  SCIPfreeBlockMemoryArray(scip, &consdata->nonlincoefs, consdata->nexprtrees);
1011  SCIPfreeBlockMemoryArray(scip, &consdata->curvatures, consdata->nexprtrees);
1012  consdata->nexprtrees = 0;
1013  }
1014 
1015  SCIP_CALL( consdataAddExprtrees(scip, consdata, nexprtrees, exprtrees, coefs, copytrees) );
1016 
1017  return SCIP_OKAY;
1018 }
1019 
1020 /** ensures, that linear vars and coefs arrays can store at least num entries */
1021 static
1022 SCIP_RETCODE consdataEnsureLinearVarsSize(
1023  SCIP* scip, /**< SCIP data structure */
1024  SCIP_CONSDATA* consdata, /**< nonlinear constraint data */
1025  int num /**< minimum number of entries to store */
1026  )
1027 {
1028  assert(scip != NULL);
1029  assert(consdata != NULL);
1030  assert(consdata->nlinvars <= consdata->linvarssize);
1031 
1032  if( num > consdata->linvarssize )
1033  {
1034  int newsize;
1035 
1036  newsize = SCIPcalcMemGrowSize(scip, num);
1037  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->linvars, consdata->linvarssize, newsize) );
1038  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->lincoefs, consdata->linvarssize, newsize) );
1039  if( consdata->lineventdata != NULL )
1040  {
1041  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->lineventdata, consdata->linvarssize, newsize) );
1042  }
1043  consdata->linvarssize = newsize;
1044  }
1045  assert(num <= consdata->linvarssize);
1046 
1047  return SCIP_OKAY;
1048 }
1049 
1050 /** creates constraint data structure */
1051 static
1052 SCIP_RETCODE consdataCreate(
1053  SCIP* scip, /**< SCIP data structure */
1054  SCIP_CONSDATA** consdata, /**< a buffer to store pointer to new constraint data */
1055  SCIP_Real lhs, /**< left hand side of constraint */
1056  SCIP_Real rhs, /**< right hand side of constraint */
1057  int nlinvars, /**< number of linear variables */
1058  SCIP_VAR** linvars, /**< array of linear variables */
1059  SCIP_Real* lincoefs, /**< array of coefficients of linear variables */
1060  int nexprtrees, /**< number of expression trees */
1061  SCIP_EXPRTREE** exprtrees, /**< expression trees */
1062  SCIP_Real* nonlincoefs, /**< coefficients of expression trees */
1063  SCIP_Bool capturevars /**< whether we should capture variables */
1064  )
1065 {
1066  int i;
1067 
1068  assert(scip != NULL);
1069  assert(consdata != NULL);
1070 
1071  assert(nlinvars == 0 || linvars != NULL);
1072  assert(nlinvars == 0 || lincoefs != NULL);
1073  assert(nexprtrees == 0 || exprtrees != NULL);
1074  assert(nexprtrees == 0 || nonlincoefs != NULL);
1075 
1076  SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );
1077  BMSclearMemory(*consdata);
1078 
1079  (*consdata)->minlinactivity = SCIP_INVALID;
1080  (*consdata)->maxlinactivity = SCIP_INVALID;
1081  (*consdata)->minlinactivityinf = -1;
1082  (*consdata)->maxlinactivityinf = -1;
1083 
1084  (*consdata)->lhs = lhs;
1085  (*consdata)->rhs = rhs;
1086 
1087  if( nlinvars > 0 )
1088  {
1089  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->linvars, linvars, nlinvars) );
1090  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->lincoefs, lincoefs, nlinvars) );
1091  (*consdata)->nlinvars = nlinvars;
1092  (*consdata)->linvarssize = nlinvars;
1093 
1094  if( capturevars )
1095  for( i = 0; i < nlinvars; ++i )
1096  {
1097  SCIP_CALL( SCIPcaptureVar(scip, linvars[i]) );
1098  }
1099  }
1100  else
1101  {
1102  (*consdata)->linvarssorted = TRUE;
1103  (*consdata)->linvarsmerged = TRUE;
1104  }
1105 
1106  SCIP_CALL( consdataSetExprtrees(scip, *consdata, nexprtrees, exprtrees, nonlincoefs, TRUE) );
1107 
1108  (*consdata)->linvar_maydecrease = -1;
1109  (*consdata)->linvar_mayincrease = -1;
1110 
1111  (*consdata)->activity = SCIP_INVALID;
1112  (*consdata)->lhsviol = SCIPisInfinity(scip, -lhs) ? 0.0 : SCIP_INVALID;
1113  (*consdata)->rhsviol = SCIPisInfinity(scip, rhs) ? 0.0 : SCIP_INVALID;
1114 
1115  return SCIP_OKAY;
1116 }
1117 
1118 /** creates empty constraint data structure */
1119 static
1120 SCIP_RETCODE consdataCreateEmpty(
1121  SCIP* scip, /**< SCIP data structure */
1122  SCIP_CONSDATA** consdata /**< a buffer to store pointer to new constraint data */
1123  )
1124 {
1125  assert(scip != NULL);
1126  assert(consdata != NULL);
1127 
1128  SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );
1129  BMSclearMemory(*consdata);
1130 
1131  (*consdata)->lhs = -SCIPinfinity(scip);
1132  (*consdata)->rhs = SCIPinfinity(scip);
1133 
1134  (*consdata)->linvarssorted = TRUE;
1135  (*consdata)->linvarsmerged = TRUE;
1136 
1137  (*consdata)->isremovedfixingslin = TRUE;
1138 
1139  (*consdata)->linvar_maydecrease = -1;
1140  (*consdata)->linvar_mayincrease = -1;
1141 
1142  (*consdata)->minlinactivity = SCIP_INVALID;
1143  (*consdata)->maxlinactivity = SCIP_INVALID;
1144  (*consdata)->minlinactivityinf = -1;
1145  (*consdata)->maxlinactivityinf = -1;
1146 
1147  (*consdata)->curvature = SCIP_EXPRCURV_LINEAR;
1148  (*consdata)->iscurvchecked = TRUE;
1149 
1150  (*consdata)->ncuts = 0;
1151 
1152  return SCIP_OKAY;
1153 }
1154 
1155 /** frees constraint data structure */
1156 static
1157 SCIP_RETCODE consdataFree(
1158  SCIP* scip, /**< SCIP data structure */
1159  SCIP_CONSDATA** consdata /**< pointer to constraint data to free */
1160  )
1161 {
1162  int i;
1163 
1164  assert(scip != NULL);
1165  assert(consdata != NULL);
1166  assert(*consdata != NULL);
1167 
1168  /* release linear variables and free linear part */
1169  if( (*consdata)->linvarssize > 0 )
1170  {
1171  for( i = 0; i < (*consdata)->nlinvars; ++i )
1172  {
1173  assert((*consdata)->lineventdata == NULL || (*consdata)->lineventdata[i] == NULL); /* variable events should have been dropped earlier */
1174  SCIP_CALL( SCIPreleaseVar(scip, &(*consdata)->linvars[i]) );
1175  }
1176  SCIPfreeBlockMemoryArray(scip, &(*consdata)->linvars, (*consdata)->linvarssize);
1177  SCIPfreeBlockMemoryArray(scip, &(*consdata)->lincoefs, (*consdata)->linvarssize);
1178  SCIPfreeBlockMemoryArrayNull(scip, &(*consdata)->lineventdata, (*consdata)->linvarssize);
1179  }
1180  assert((*consdata)->linvars == NULL);
1181  assert((*consdata)->lincoefs == NULL);
1182  assert((*consdata)->lineventdata == NULL);
1183 
1184  if( (*consdata)->nexprtrees > 0 )
1185  {
1186  assert((*consdata)->exprtrees != NULL);
1187  assert((*consdata)->nonlincoefs != NULL);
1188  assert((*consdata)->curvatures != NULL);
1189  for( i = 0; i < (*consdata)->nexprtrees; ++i )
1190  {
1191  assert((*consdata)->exprtrees[i] != NULL);
1192  SCIP_CALL( SCIPexprtreeFree(&(*consdata)->exprtrees[i]) );
1193  assert((*consdata)->exprtrees[i] == NULL);
1194  }
1195  SCIPfreeBlockMemoryArray(scip, &(*consdata)->exprtrees, (*consdata)->nexprtrees);
1196  SCIPfreeBlockMemoryArray(scip, &(*consdata)->nonlincoefs, (*consdata)->nexprtrees);
1197  SCIPfreeBlockMemoryArray(scip, &(*consdata)->curvatures, (*consdata)->nexprtrees);
1198  }
1199  assert((*consdata)->exprtrees == NULL);
1200  assert((*consdata)->nonlincoefs == NULL);
1201  assert((*consdata)->curvatures == NULL);
1202 
1203  /* free nonlinear row representation */
1204  if( (*consdata)->nlrow != NULL )
1205  {
1206  SCIP_CALL( SCIPreleaseNlRow(scip, &(*consdata)->nlrow) );
1207  }
1208 
1209  SCIPfreeBlockMemory(scip, consdata);
1210  *consdata = NULL;
1211 
1212  return SCIP_OKAY;
1213 }
1214 
1215 /** sorts linear part of constraint data */
1216 static
1217 void consdataSortLinearVars(
1218  SCIP_CONSDATA* consdata /**< nonlinear constraint data */
1219  )
1220 {
1221  assert(consdata != NULL);
1222 
1223  if( consdata->linvarssorted )
1224  return;
1225 
1226  if( consdata->nlinvars <= 1 )
1227  {
1228  consdata->linvarssorted = TRUE;
1229  return;
1230  }
1231 
1232  if( consdata->lineventdata == NULL )
1233  {
1234  SCIPsortPtrReal((void**)consdata->linvars, consdata->lincoefs, SCIPvarComp, consdata->nlinvars);
1235  }
1236  else
1237  {
1238  int i;
1239 
1240  SCIPsortPtrPtrReal((void**)consdata->linvars, (void**)consdata->lineventdata, consdata->lincoefs, SCIPvarComp, consdata->nlinvars);
1241 
1242  /* update variable indices in event data */
1243  for( i = 0; i < consdata->nlinvars; ++i )
1244  if( consdata->lineventdata[i] != NULL )
1245  consdata->lineventdata[i]->varidx = i;
1246  }
1247 
1248  consdata->linvarssorted = TRUE;
1249 }
1250 
1251 /* this function is currently not needed, but also to nice to be deleted, so it is only deactivated */
1252 #ifdef SCIP_DISABLED_CODE
1253 /** returns the position of variable in the linear coefficients array of a constraint, or -1 if not found */
1254 static
1255 int consdataFindLinearVar(
1256  SCIP_CONSDATA* consdata, /**< nonlinear constraint data */
1257  SCIP_VAR* var /**< variable to search for */
1258  )
1259 {
1260  int pos;
1261 
1262  assert(consdata != NULL);
1263  assert(var != NULL);
1264 
1265  if( consdata->nlinvars == 0 )
1266  return -1;
1267 
1268  consdataSortLinearVars(consdata);
1269 
1270  if( !SCIPsortedvecFindPtr((void**)consdata->linvars, SCIPvarComp, (void*)var, consdata->nlinvars, &pos) )
1271  pos = -1;
1272 
1273  return pos;
1274 }
1275 #endif
1276 
1277 /** moves a linear variable from one position to another */
1278 static
1279 void consdataMoveLinearVar(
1280  SCIP_CONSDATA* consdata, /**< constraint data */
1281  int oldpos, /**< position of variable that shall be moved */
1282  int newpos /**< new position of variable */
1283  )
1284 {
1285  assert(consdata != NULL);
1286  assert(oldpos >= 0);
1287  assert(oldpos < consdata->nlinvars);
1288  assert(newpos >= 0);
1289  assert(newpos < consdata->linvarssize);
1290 
1291  if( newpos == oldpos )
1292  return;
1293 
1294  consdata->linvars [newpos] = consdata->linvars [oldpos];
1295  consdata->lincoefs[newpos] = consdata->lincoefs[oldpos];
1296 
1297  if( consdata->lineventdata != NULL )
1298  {
1299  assert(newpos >= consdata->nlinvars || consdata->lineventdata[newpos] == NULL);
1300 
1301  consdata->lineventdata[newpos] = consdata->lineventdata[oldpos];
1302  consdata->lineventdata[newpos]->varidx = newpos;
1303 
1304  consdata->lineventdata[oldpos] = NULL;
1305  }
1306 
1307  consdata->linvarssorted = FALSE;
1308 }
1309 
1310 /** adds linear coefficient in nonlinear constraint */
1311 static
1312 SCIP_RETCODE addLinearCoef(
1313  SCIP* scip, /**< SCIP data structure */
1314  SCIP_CONS* cons, /**< nonlinear constraint */
1315  SCIP_VAR* var, /**< variable of constraint entry */
1316  SCIP_Real coef /**< coefficient of constraint entry */
1317  )
1318 {
1319  SCIP_CONSDATA* consdata;
1320  SCIP_Bool transformed;
1321 
1322  assert(scip != NULL);
1323  assert(cons != NULL);
1324  assert(var != NULL);
1325 
1326  /* ignore coefficient if it is nearly zero */
1327  if( SCIPisZero(scip, coef) )
1328  return SCIP_OKAY;
1329 
1330  consdata = SCIPconsGetData(cons);
1331  assert(consdata != NULL);
1332 
1333  /* are we in the transformed problem? */
1334  transformed = SCIPconsIsTransformed(cons);
1335 
1336  /* always use transformed variables in transformed constraints */
1337  if( transformed )
1338  {
1339  SCIP_CALL( SCIPgetTransformedVar(scip, var, &var) );
1340  }
1341  assert(var != NULL);
1342  assert(transformed == SCIPvarIsTransformed(var));
1343 
1344  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars+1) );
1345  consdata->linvars [consdata->nlinvars] = var;
1346  consdata->lincoefs[consdata->nlinvars] = coef;
1347 
1348  ++consdata->nlinvars;
1349 
1350  /* catch variable events */
1351  if( SCIPconsIsEnabled(cons) )
1352  {
1353  /* catch bound change events of variable */
1354  SCIP_CALL( catchLinearVarEvents(scip, cons, consdata->nlinvars-1) );
1355  }
1356 
1357  /* invalidate activity information */
1358  consdata->activity = SCIP_INVALID;
1359  consdata->minlinactivity = SCIP_INVALID;
1360  consdata->maxlinactivity = SCIP_INVALID;
1361  consdata->minlinactivityinf = -1;
1362  consdata->maxlinactivityinf = -1;
1363 
1364  /* invalidate nonlinear row */
1365  if( consdata->nlrow != NULL )
1366  {
1367  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1368  }
1369 
1370  /* install rounding locks for new variable */
1371  SCIP_CALL( lockLinearVariable(scip, cons, var, coef) );
1372 
1373  /* capture new variable */
1374  SCIP_CALL( SCIPcaptureVar(scip, var) );
1375 
1376  consdata->ispresolved = FALSE;
1377  consdata->isremovedfixingslin = consdata->isremovedfixingslin && SCIPvarIsActive(var);
1378  if( consdata->nlinvars == 1 )
1379  consdata->linvarssorted = TRUE;
1380  else
1381  consdata->linvarssorted = consdata->linvarssorted &&
1382  (SCIPvarCompare(consdata->linvars[consdata->nlinvars-2], consdata->linvars[consdata->nlinvars-1]) == -1);
1383  /* always set too FALSE since the new linear variable should be checked if already existing as quad var term */
1384  consdata->linvarsmerged = FALSE;
1385 
1386  return SCIP_OKAY;
1387 }
1388 
1389 /** deletes linear coefficient at given position from nonlinear constraint data */
1390 static
1391 SCIP_RETCODE delLinearCoefPos(
1392  SCIP* scip, /**< SCIP data structure */
1393  SCIP_CONS* cons, /**< nonlinear constraint */
1394  int pos /**< position of coefficient to delete */
1395  )
1396 {
1397  SCIP_CONSDATA* consdata;
1398  SCIP_VAR* var;
1399  SCIP_Real coef;
1400 
1401  assert(scip != NULL);
1402  assert(cons != NULL);
1403 
1404  consdata = SCIPconsGetData(cons);
1405  assert(consdata != NULL);
1406  assert(0 <= pos && pos < consdata->nlinvars);
1407 
1408  var = consdata->linvars[pos];
1409  coef = consdata->lincoefs[pos];
1410  assert(var != NULL);
1411 
1412  /* remove rounding locks for deleted variable */
1413  SCIP_CALL( unlockLinearVariable(scip, cons, var, coef) );
1414 
1415  /* if constraint is enabled, drop the events on the variable */
1416  if( SCIPconsIsEnabled(cons) )
1417  {
1418  /* drop bound change events of variable */
1419  SCIP_CALL( dropLinearVarEvents(scip, cons, pos) );
1420  }
1421 
1422  /* release variable */
1423  SCIP_CALL( SCIPreleaseVar(scip, &consdata->linvars[pos]) );
1424 
1425  /* move the last variable to the free slot */
1426  consdataMoveLinearVar(consdata, consdata->nlinvars-1, pos);
1427 
1428  --consdata->nlinvars;
1429 
1430  /* invalidate activity */
1431  consdata->activity = SCIP_INVALID;
1432  consdata->minlinactivity = SCIP_INVALID;
1433  consdata->maxlinactivity = SCIP_INVALID;
1434  consdata->minlinactivityinf = -1;
1435  consdata->maxlinactivityinf = -1;
1436 
1437  /* invalidate nonlinear row */
1438  if( consdata->nlrow != NULL )
1439  {
1440  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1441  }
1442 
1443  consdata->ispresolved = FALSE;
1444 
1445  return SCIP_OKAY;
1446 }
1447 
1448 /** changes linear coefficient value at given position of nonlinear constraint */
1449 static
1450 SCIP_RETCODE chgLinearCoefPos(
1451  SCIP* scip, /**< SCIP data structure */
1452  SCIP_CONS* cons, /**< nonlinear constraint */
1453  int pos, /**< position of linear coefficient to change */
1454  SCIP_Real newcoef /**< new value of linear coefficient */
1455  )
1456 {
1457  SCIP_CONSDATA* consdata;
1458  SCIP_VAR* var;
1459  SCIP_Real coef;
1460 
1461  assert(scip != NULL);
1462  assert(cons != NULL);
1463  assert(!SCIPisZero(scip, newcoef));
1464 
1465  consdata = SCIPconsGetData(cons);
1466  assert(consdata != NULL);
1467  assert(0 <= pos);
1468  assert(pos < consdata->nlinvars);
1469  assert(!SCIPisZero(scip, newcoef));
1470 
1471  var = consdata->linvars[pos];
1472  coef = consdata->lincoefs[pos];
1473  assert(var != NULL);
1474  assert(SCIPconsIsTransformed(cons) == SCIPvarIsTransformed(var));
1475 
1476  /* invalidate activity */
1477  consdata->activity = SCIP_INVALID;
1478  consdata->minlinactivity = SCIP_INVALID;
1479  consdata->maxlinactivity = SCIP_INVALID;
1480  consdata->minlinactivityinf = -1;
1481  consdata->maxlinactivityinf = -1;
1482 
1483  /* invalidate nonlinear row */
1484  if( consdata->nlrow != NULL )
1485  {
1486  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1487  }
1488 
1489  /* if necessary, remove the rounding locks and event catching of the variable */
1490  if( newcoef * coef < 0.0 )
1491  {
1492  if( SCIPconsIsLocked(cons) )
1493  {
1494  assert(SCIPconsIsTransformed(cons));
1495 
1496  /* remove rounding locks for variable with old coefficient */
1497  SCIP_CALL( unlockLinearVariable(scip, cons, var, coef) );
1498  }
1499 
1500  if( SCIPconsIsEnabled(cons) )
1501  {
1502  /* drop bound change events of variable */
1503  SCIP_CALL( dropLinearVarEvents(scip, cons, pos) );
1504  }
1505  }
1506 
1507  /* change the coefficient */
1508  consdata->lincoefs[pos] = newcoef;
1509 
1510  /* if necessary, install the rounding locks and event catching of the variable again */
1511  if( newcoef * coef < 0.0 )
1512  {
1513  if( SCIPconsIsLocked(cons) )
1514  {
1515  /* install rounding locks for variable with new coefficient */
1516  SCIP_CALL( lockLinearVariable(scip, cons, var, newcoef) );
1517  }
1518 
1519  if( SCIPconsIsEnabled(cons) )
1520  {
1521  /* catch bound change events of variable */
1522  SCIP_CALL( catchLinearVarEvents(scip, cons, pos) );
1523  }
1524  }
1525 
1526  consdata->ispresolved = FALSE;
1527 
1528  return SCIP_OKAY;
1529 }
1530 
1531 
1532 /* merges entries with same linear variable into one entry and cleans up entries with coefficient 0.0 */
1533 static
1534 SCIP_RETCODE mergeAndCleanLinearVars(
1535  SCIP* scip, /**< SCIP data structure */
1536  SCIP_CONS* cons /**< nonlinear constraint */
1537  )
1538 {
1539  SCIP_CONSDATA* consdata;
1540  SCIP_Real newcoef;
1541  int i;
1542  int j;
1543 
1544  assert(scip != NULL);
1545  assert(cons != NULL);
1546 
1547  consdata = SCIPconsGetData(cons);
1548 
1549  if( consdata->linvarsmerged )
1550  return SCIP_OKAY;
1551 
1552  if( consdata->nlinvars == 0 )
1553  {
1554  consdata->linvarsmerged = TRUE;
1555  return SCIP_OKAY;
1556  }
1557 
1558  i = 0;
1559  while( i < consdata->nlinvars )
1560  {
1561  /* make sure linear variables are sorted (do this in every round, since we may move variables around) */
1562  consdataSortLinearVars(consdata);
1563 
1564  /* sum up coefficients that correspond to variable i */
1565  newcoef = consdata->lincoefs[i];
1566  for( j = i+1; j < consdata->nlinvars && consdata->linvars[i] == consdata->linvars[j]; ++j )
1567  newcoef += consdata->lincoefs[j];
1568  /* delete the additional variables in backward order */
1569  for( j = j-1; j > i; --j )
1570  {
1571  SCIP_CALL( delLinearCoefPos(scip, cons, j) );
1572  }
1573 
1574  /* delete also entry at position i, if it became zero (or was zero before) */
1575  if( SCIPisZero(scip, newcoef) )
1576  {
1577  SCIP_CALL( delLinearCoefPos(scip, cons, i) );
1578  }
1579  else
1580  {
1581  SCIP_CALL( chgLinearCoefPos(scip, cons, i, newcoef) );
1582  ++i;
1583  }
1584  }
1585 
1586  consdata->linvarsmerged = TRUE;
1587 
1588  return SCIP_OKAY;
1589 }
1590 
1591 /** removes fixes (or aggregated) linear variables from a nonlinear constraint */
1592 static
1593 SCIP_RETCODE removeFixedLinearVariables(
1594  SCIP* scip, /**< SCIP data structure */
1595  SCIP_CONS* cons /**< nonlinearconstraint */
1596  )
1597 {
1598  SCIP_CONSDATA* consdata;
1599  SCIP_Real coef;
1600  SCIP_Real offset;
1601  SCIP_VAR* var;
1602  int i;
1603  int j;
1604 
1605  assert(scip != NULL);
1606  assert(cons != NULL);
1607 
1608  consdata = SCIPconsGetData(cons);
1609 
1610  if( !consdata->isremovedfixingslin )
1611  {
1612  i = 0;
1613  while( i < consdata->nlinvars )
1614  {
1615  var = consdata->linvars[i];
1616 
1617  if( SCIPvarIsActive(var) )
1618  {
1619  ++i;
1620  continue;
1621  }
1622 
1623  coef = consdata->lincoefs[i];
1624  offset = 0.0;
1625 
1626  SCIP_CALL( SCIPgetProbvarSum(scip, &var, &coef, &offset) );
1627 
1628  SCIPdebugMsg(scip, " linear term %g*<%s> is replaced by %g * <%s> + %g\n", consdata->lincoefs[i], SCIPvarGetName(consdata->linvars[i]), coef, SCIPvarGetName(var), offset);
1629 
1630  /* delete previous variable (this will move another variable to position i) */
1631  SCIP_CALL( delLinearCoefPos(scip, cons, i) );
1632 
1633  /* put constant part into bounds */
1634  if( offset != 0.0 )
1635  {
1636  if( !SCIPisInfinity(scip, -consdata->lhs) )
1637  {
1638  consdata->lhs -= offset;
1639  assert(!SCIPisInfinity(scip, REALABS(consdata->lhs)));
1640  }
1641  if( !SCIPisInfinity(scip, consdata->rhs) )
1642  {
1643  consdata->rhs -= offset;
1644  assert(!SCIPisInfinity(scip, REALABS(consdata->rhs)));
1645  }
1646  }
1647 
1648  /* nothing left to do if variable had been fixed */
1649  if( coef == 0.0 )
1650  continue;
1651 
1652  /* if GetProbvar gave a linear variable, just add it
1653  * if it's a multilinear variable, add it's disaggregated variables */
1654  if( SCIPvarIsActive(var) )
1655  {
1656  SCIP_CALL( addLinearCoef(scip, cons, var, coef) );
1657  }
1658  else
1659  {
1660  int naggrs;
1661  SCIP_VAR** aggrvars;
1662  SCIP_Real* aggrscalars;
1663  SCIP_Real aggrconstant;
1664 
1665  assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_MULTAGGR);
1666 
1667  naggrs = SCIPvarGetMultaggrNVars(var);
1668  aggrvars = SCIPvarGetMultaggrVars(var);
1669  aggrscalars = SCIPvarGetMultaggrScalars(var);
1670  aggrconstant = SCIPvarGetMultaggrConstant(var);
1671 
1672  SCIP_CALL( consdataEnsureLinearVarsSize(scip, consdata, consdata->nlinvars + naggrs) );
1673 
1674  for( j = 0; j < naggrs; ++j )
1675  {
1676  SCIP_CALL( addLinearCoef(scip, cons, aggrvars[j], coef * aggrscalars[j]) );
1677  }
1678 
1679  if( aggrconstant != 0.0 )
1680  {
1681  if( !SCIPisInfinity(scip, -consdata->lhs) )
1682  {
1683  consdata->lhs -= coef * aggrconstant;
1684  assert(!SCIPisInfinity(scip, REALABS(consdata->lhs)));
1685  }
1686  if( !SCIPisInfinity(scip, consdata->rhs) )
1687  {
1688  consdata->rhs -= coef * aggrconstant;
1689  assert(!SCIPisInfinity(scip, REALABS(consdata->rhs)));
1690  }
1691  }
1692  }
1693  }
1694 
1695  SCIP_CALL( mergeAndCleanLinearVars(scip, cons) );
1696 
1697  consdata->isremovedfixingslin = TRUE;
1698  }
1699 
1700  SCIPdebugMsg(scip, "removed fixations of linear variables from <%s>\n -> ", SCIPconsGetName(cons));
1701  SCIPdebugPrintCons(scip, cons, NULL);
1702 
1703 #ifndef NDEBUG
1704  for( i = 0; i < consdata->nlinvars; ++i )
1705  assert(SCIPvarIsActive(consdata->linvars[i]));
1706 #endif
1707 
1708  return SCIP_OKAY;
1709 }
1710 
1711 /** removes fixed variables from expression graph */
1712 static
1713 SCIP_RETCODE removeFixedNonlinearVariables(
1714  SCIP* scip, /**< SCIP data structure */
1715  SCIP_CONSHDLR* conshdlr /**< constraint handler */
1716  )
1717 {
1718  SCIP_CONSHDLRDATA* conshdlrdata;
1719  SCIP_VAR* var;
1720  SCIP_VAR** vars;
1721  SCIP_Real* coefs;
1722  int nvars;
1723  int varssize;
1724  SCIP_Real constant;
1725  int i;
1726  int requsize;
1727  SCIPdebug( int j );
1728 
1729  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1730  assert(conshdlrdata != NULL);
1731  assert(conshdlrdata->exprgraph != NULL);
1732 
1733  if( conshdlrdata->isremovedfixings )
1734  return SCIP_OKAY;
1735 
1736  varssize = 5;
1737  SCIP_CALL( SCIPallocBufferArray(scip, &vars, varssize) );
1738  SCIP_CALL( SCIPallocBufferArray(scip, &coefs, varssize) );
1739 
1740  i = 0;
1741  while( i < SCIPexprgraphGetNVars(conshdlrdata->exprgraph) )
1742  {
1743  var = (SCIP_VAR*)SCIPexprgraphGetVars(conshdlrdata->exprgraph)[i];
1744  if( SCIPvarIsActive(var) )
1745  {
1746  ++i;
1747  continue;
1748  }
1749 
1750  vars[0] = var;
1751  coefs[0] = 1.0;
1752  constant = 0.0;
1753  nvars = 1;
1754  SCIP_CALL( SCIPgetProbvarLinearSum(scip, vars, coefs, &nvars, varssize, &constant, &requsize, TRUE) );
1755 
1756  if( requsize > varssize )
1757  {
1758  SCIP_CALL( SCIPreallocBufferArray(scip, &vars, requsize) );
1759  SCIP_CALL( SCIPreallocBufferArray(scip, &coefs, requsize) );
1760  varssize = requsize;
1761 
1762  SCIP_CALL( SCIPgetProbvarLinearSum(scip, vars, coefs, &nvars, varssize, &constant, &requsize, TRUE) );
1763  }
1764 
1765 #ifdef SCIP_DEBUG
1766  SCIPdebugMsg(scip, "replace fixed variable <%s> by %g", SCIPvarGetName(var), constant);
1767  for( j = 0; j < nvars; ++j )
1768  {
1769  SCIPdebugMsgPrint(scip, " %+g <%s>", coefs[j], SCIPvarGetName(vars[j]));
1770  }
1771  SCIPdebugMsgPrint(scip, "\n");
1772 #endif
1773 
1774  SCIP_CALL( SCIPexprgraphReplaceVarByLinearSum(conshdlrdata->exprgraph, var, nvars, coefs, (void**)vars, constant) );
1775 
1776  i = 0;
1777  }
1778 
1779  SCIPfreeBufferArray(scip, &coefs);
1780  SCIPfreeBufferArray(scip, &vars);
1781 
1782  conshdlrdata->isremovedfixings = TRUE;
1783 
1784  return SCIP_OKAY;
1785 }
1786 
1787 /** moves constant and linear part from expression graph node into constraint sides and linear part
1788  * frees expression graph node if expression is constant or linear */
1789 static
1790 SCIP_RETCODE splitOffLinearPart(
1791  SCIP* scip, /**< SCIP data structure */
1792  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1793  SCIP_CONS* cons, /**< nonlinear constraint */
1794  SCIP_Bool* infeasible /**< pointer to store whether the problem is infeasible or not */
1795  )
1796 {
1797  SCIP_CONSHDLRDATA* conshdlrdata;
1798  SCIP_CONSDATA* consdata;
1799  SCIP_VAR** linvars;
1800  SCIP_Real* lincoefs;
1801  SCIP_Real constant;
1802  int linvarssize;
1803  int nlinvars;
1804  int i;
1805 
1806  assert(scip != NULL);
1807  assert(conshdlr != NULL);
1808  assert(cons != NULL);
1809 
1810  consdata = SCIPconsGetData(cons);
1811  assert(consdata != NULL);
1812 
1813  *infeasible = FALSE;
1814 
1815  if( consdata->exprgraphnode == NULL )
1816  return SCIP_OKAY;
1817 
1818  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1819  assert(conshdlrdata != NULL);
1820  assert(conshdlrdata->exprgraph != NULL);
1821 
1822  /* number of children of expression graph node is a good upper estimate on number of linear variables */
1823  linvarssize = MAX(SCIPexprgraphGetNodeNChildren(consdata->exprgraphnode), 1); /*lint !e666*/
1824  SCIP_CALL( SCIPallocBufferArray(scip, &linvars, linvarssize) );
1825  SCIP_CALL( SCIPallocBufferArray(scip, &lincoefs, linvarssize) );
1826 
1827  /* get linear and constant part from expression graph node
1828  * releases expression graph node if not uses otherwise */
1829  SCIP_CALL( SCIPexprgraphNodeSplitOffLinear(conshdlrdata->exprgraph, &consdata->exprgraphnode, linvarssize, &nlinvars, (void**)linvars, lincoefs, &constant) );
1830 
1831  if( SCIPisInfinity(scip, constant) )
1832  {
1833  if( !SCIPisInfinity(scip, -consdata->lhs) )
1834  {
1835  /* setting constraint lhs to -infinity; this may change linear variable locks and events */
1836  for( i = 0; i < consdata->nlinvars; ++i )
1837  {
1838  if( SCIPconsIsLocked(cons) )
1839  {
1840  SCIP_CALL( unlockLinearVariable(scip, cons, consdata->linvars[i], consdata->lincoefs[i]) );
1841  }
1842  if( SCIPconsIsEnabled(cons) )
1843  {
1844  SCIP_CALL( dropLinearVarEvents(scip, cons, i) );
1845  }
1846  }
1847 
1848  consdata->lhs = -SCIPinfinity(scip);
1849 
1850  for( i = 0; i < consdata->nlinvars; ++i )
1851  {
1852  if( SCIPconsIsEnabled(cons) )
1853  {
1854  SCIP_CALL( catchLinearVarEvents(scip, cons, i) );
1855  }
1856  if( SCIPconsIsLocked(cons) )
1857  {
1858  SCIP_CALL( lockLinearVariable(scip, cons, consdata->linvars[i], consdata->lincoefs[i]) );
1859  }
1860  }
1861  }
1862 
1863  if( !SCIPisInfinity(scip, consdata->rhs) )
1864  {
1865  *infeasible = TRUE;
1866  goto TERMINATE;
1867  }
1868  }
1869  else if( SCIPisInfinity(scip, -constant) )
1870  {
1871  if( !SCIPisInfinity(scip, consdata->rhs) )
1872  {
1873  /* setting constraint rhs to infinity; this may change linear variable locks and events */
1874  for( i = 0; i < consdata->nlinvars; ++i )
1875  {
1876  if( SCIPconsIsLocked(cons) )
1877  {
1878  SCIP_CALL( unlockLinearVariable(scip, cons, consdata->linvars[i], consdata->lincoefs[i]) );
1879  }
1880  if( SCIPconsIsEnabled(cons) )
1881  {
1882  SCIP_CALL( dropLinearVarEvents(scip, cons, i) );
1883  }
1884  }
1885 
1886  consdata->rhs = SCIPinfinity(scip);
1887 
1888  for( i = 0; i < consdata->nlinvars; ++i )
1889  {
1890  if( SCIPconsIsEnabled(cons) )
1891  {
1892  SCIP_CALL( catchLinearVarEvents(scip, cons, i) );
1893  }
1894  if( SCIPconsIsLocked(cons) )
1895  {
1896  SCIP_CALL( lockLinearVariable(scip, cons, consdata->linvars[i], consdata->lincoefs[i]) );
1897  }
1898  }
1899  }
1900  if( !SCIPisInfinity(scip, -consdata->lhs) )
1901  {
1902  *infeasible = TRUE;
1903  goto TERMINATE;
1904  }
1905  }
1906  else if( constant != 0.0 )
1907  {
1908  if( !SCIPisInfinity(scip, -consdata->lhs) )
1909  {
1910  consdata->lhs -= constant;
1911  assert(!SCIPisInfinity(scip, REALABS(consdata->lhs)));
1912  }
1913  if( !SCIPisInfinity(scip, consdata->rhs) )
1914  {
1915  consdata->rhs -= constant;
1916  assert(!SCIPisInfinity(scip, REALABS(consdata->rhs)));
1917  }
1918  }
1919 
1920 TERMINATE:
1921  for( i = 0; i < nlinvars; ++i )
1922  {
1923  SCIP_CALL( addLinearCoef(scip, cons, linvars[i], lincoefs[i]) );
1924  }
1925 
1926  SCIPfreeBufferArray(scip, &lincoefs);
1927  SCIPfreeBufferArray(scip, &linvars);
1928 
1929  /* @todo linear variables that are also children of exprgraphnode could be moved into the expression graph for certain nonlinear operators (quadratic, polynomial), since that may allow better bound tightening */
1930 
1931  return SCIP_OKAY;
1932 }
1933 
1934 /** create a nonlinear row representation of the constraint and stores them in consdata */
1935 static
1936 SCIP_RETCODE createNlRow(
1937  SCIP* scip, /**< SCIP data structure */
1938  SCIP_CONS* cons /**< nonlinear constraint */
1939  )
1940 {
1941  SCIP_CONSDATA* consdata;
1942 
1943  assert(scip != NULL);
1944  assert(cons != NULL);
1945 
1946  consdata = SCIPconsGetData(cons);
1947  assert(consdata != NULL);
1948 
1949  if( consdata->nlrow != NULL )
1950  {
1951  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
1952  }
1953 
1954  if( consdata->nexprtrees == 0 )
1955  {
1956  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
1957  consdata->nlinvars, consdata->linvars, consdata->lincoefs,
1958  0, NULL, 0, NULL,
1959  NULL, consdata->lhs, consdata->rhs,
1960  consdata->curvature) );
1961  }
1962  else if( consdata->nexprtrees == 1 && consdata->nonlincoefs[0] == 1.0 )
1963  {
1964  assert(consdata->exprtrees[0] != NULL);
1965  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
1966  consdata->nlinvars, consdata->linvars, consdata->lincoefs,
1967  0, NULL, 0, NULL,
1968  consdata->exprtrees[0], consdata->lhs, consdata->rhs,
1969  consdata->curvature) );
1970  }
1971  else
1972  {
1973  /* since expression trees may share variable, we cannot easily sum them up,
1974  * but we can request a single expression tree from the expression graph
1975  */
1976  SCIP_CONSHDLRDATA* conshdlrdata;
1977  SCIP_EXPRTREE* exprtree;
1978 
1979  assert(consdata->exprgraphnode != NULL); /* since nexprtrees > 0 */
1980  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
1981  assert(conshdlrdata != NULL);
1982 
1983  SCIP_CALL( SCIPexprgraphGetTree(conshdlrdata->exprgraph, consdata->exprgraphnode, &exprtree) );
1984  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
1985  consdata->nlinvars, consdata->linvars, consdata->lincoefs,
1986  0, NULL, 0, NULL,
1987  exprtree, consdata->lhs, consdata->rhs,
1988  consdata->curvature) );
1989  SCIP_CALL( SCIPexprtreeFree(&exprtree) );
1990  }
1991 
1992  return SCIP_OKAY;
1993 }
1994 
1995 /** tries to automatically convert a nonlinear constraint (or a part of it) into a more specific and more specialized constraint */
1996 static
1997 SCIP_RETCODE presolveUpgrade(
1998  SCIP* scip, /**< SCIP data structure */
1999  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
2000  SCIP_CONS* cons, /**< source constraint to try to convert */
2001  SCIP_Bool* upgraded, /**< buffer to store whether constraint was upgraded */
2002  int* nupgdconss, /**< buffer to increase if constraint was upgraded */
2003  int* naddconss /**< buffer to increase with number of additional constraints created during upgrade */
2004  )
2005 {
2006  SCIP_CONSHDLRDATA* conshdlrdata;
2007  SCIP_CONS** upgdconss;
2008  int upgdconsssize;
2009  int nupgdconss_;
2010  int i;
2011 
2012  assert(scip != NULL);
2013  assert(conshdlr != NULL);
2014  assert(cons != NULL);
2015  assert(!SCIPconsIsModifiable(cons));
2016  assert(upgraded != NULL);
2017  assert(nupgdconss != NULL);
2018  assert(naddconss != NULL);
2019 
2020  *upgraded = FALSE;
2021 
2022  nupgdconss_ = 0;
2023 
2024  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2025  assert(conshdlrdata != NULL);
2026 
2027  /* if there are no upgrade methods, we can stop */
2028  if( conshdlrdata->nnlconsupgrades == 0 )
2029  return SCIP_OKAY;
2030 
2031  /* set debug solution in expression graph and evaluate nodes, so reformulation methods can compute debug solution values for new auxiliary variables */
2032 #ifdef WITH_DEBUG_SOLUTION
2033  if( SCIPdebugIsMainscip(scip) )
2034  {
2035  SCIP_Real* varvals;
2036 
2037  SCIP_CALL( SCIPallocBufferArray(scip, &varvals, SCIPexprgraphGetNVars(conshdlrdata->exprgraph)) );
2038 
2039  for( i = 0; i < SCIPexprgraphGetNVars(conshdlrdata->exprgraph); ++i )
2040  SCIP_CALL( SCIPdebugGetSolVal(scip, (SCIP_VAR*)SCIPexprgraphGetVars(conshdlrdata->exprgraph)[i], &varvals[i]) );
2041 
2042  SCIP_CALL( SCIPexprgraphEval(conshdlrdata->exprgraph, varvals) );
2043 
2044  SCIPfreeBufferArray(scip, &varvals);
2045  }
2046 #endif
2047 
2048  upgdconsssize = 2;
2049  SCIP_CALL( SCIPallocBufferArray(scip, &upgdconss, upgdconsssize) );
2050 
2051  /* call the upgrading methods */
2052 
2053  SCIPdebugMsg(scip, "upgrading nonlinear constraint <%s> (up to %d upgrade methods):\n",
2054  SCIPconsGetName(cons), conshdlrdata->nnlconsupgrades);
2055  SCIPdebugPrintCons(scip, cons, NULL);
2056 
2057  /* try all upgrading methods in priority order in case the upgrading step is enable */
2058  for( i = 0; i < conshdlrdata->nnlconsupgrades; ++i )
2059  {
2060  if( !conshdlrdata->nlconsupgrades[i]->active )
2061  continue;
2062  if( conshdlrdata->nlconsupgrades[i]->nlconsupgd == NULL )
2063  continue;
2064 
2065  SCIP_CALL( conshdlrdata->nlconsupgrades[i]->nlconsupgd(scip, cons, &nupgdconss_, upgdconss, upgdconsssize) );
2066 
2067  while( nupgdconss_ < 0 )
2068  {
2069  /* upgrade function requires more memory: resize upgdconss and call again */
2070  assert(-nupgdconss_ > upgdconsssize);
2071  upgdconsssize = -nupgdconss_;
2072  SCIP_CALL( SCIPreallocBufferArray(scip, &upgdconss, -nupgdconss_) );
2073 
2074  SCIP_CALL( conshdlrdata->nlconsupgrades[i]->nlconsupgd(scip, cons, &nupgdconss_, upgdconss, upgdconsssize) );
2075 
2076  assert(nupgdconss_ != 0);
2077  }
2078 
2079  if( nupgdconss_ > 0 )
2080  {
2081  /* got upgrade */
2082  int j;
2083 
2084  SCIPdebugMsg(scip, " -> upgraded to %d constraints:\n", nupgdconss_);
2085 
2086  /* add the upgraded constraints to the problem and forget them */
2087  for( j = 0; j < nupgdconss_; ++j )
2088  {
2089  SCIPdebugMsgPrint(scip, "\t");
2090  SCIPdebugPrintCons(scip, upgdconss[j], NULL);
2091 
2092  SCIP_CALL( SCIPaddCons(scip, upgdconss[j]) ); /*lint !e613*/
2093  SCIP_CALL( SCIPreleaseCons(scip, &upgdconss[j]) ); /*lint !e613*/
2094  }
2095 
2096  /* count the first upgrade constraint as constraint upgrade and the remaining ones as added constraints */
2097  *nupgdconss += 1;
2098  *naddconss += nupgdconss_ - 1;
2099  *upgraded = TRUE;
2100 
2101  /* delete upgraded constraint */
2102  SCIPdebugMsg(scip, "delete constraint <%s> after upgrade\n", SCIPconsGetName(cons));
2103  SCIP_CALL( SCIPdelCons(scip, cons) );
2104 
2105  break;
2106  }
2107  }
2108 
2109  SCIPfreeBufferArray(scip, &upgdconss);
2110 
2111  return SCIP_OKAY;
2112 }
2113 
2114 /** checks a nonlinear constraint for convexity and/or concavity */
2115 static
2116 SCIP_RETCODE checkCurvature(
2117  SCIP* scip, /**< SCIP data structure */
2118  SCIP_CONS* cons, /**< nonlinear constraint */
2119  SCIP_Bool expensivechecks, /**< whether also expensive checks should be executed */
2120  SCIP_Bool assumeconvex /**< whether to assume convexity in inequalities */
2121  )
2122 {
2123  SCIP_CONSDATA* consdata;
2124  SCIP_INTERVAL* varbounds;
2125  int varboundssize;
2126  SCIP_VAR* var;
2127  int i;
2128  int j;
2129 
2130  assert(scip != NULL);
2131  assert(cons != NULL);
2132 
2133  consdata = SCIPconsGetData(cons);
2134  assert(consdata != NULL);
2135 
2136  if( consdata->iscurvchecked )
2137  return SCIP_OKAY;
2138 
2139  SCIPdebugMsg(scip, "Checking curvature of constraint <%s>\n", SCIPconsGetName(cons));
2140 
2141  consdata->curvature = SCIP_EXPRCURV_LINEAR;
2142  consdata->iscurvchecked = TRUE;
2143 
2144  varbounds = NULL;
2145  varboundssize = 0;
2146 
2147  for( i = 0; i < consdata->nexprtrees; ++i )
2148  {
2149  assert(consdata->exprtrees[i] != NULL);
2150  assert(SCIPexprtreeGetNVars(consdata->exprtrees[i]) > 0 );
2151 
2152  if( assumeconvex )
2153  {
2154  /* for constraints a*f(x) <= rhs, we assume that it is convex */
2155  if( SCIPisInfinity(scip, -consdata->lhs) )
2156  consdata->curvatures[i] = SCIP_EXPRCURV_CONVEX;
2157 
2158  /* for constraints lhs <= a*f(x), we assume that it is concave */
2159  if( SCIPisInfinity(scip, consdata->rhs) )
2160  consdata->curvatures[i] = SCIP_EXPRCURV_CONCAVE;
2161  }
2162  else
2163  {
2164  if( varboundssize == 0 )
2165  {
2166  SCIP_CALL( SCIPallocBufferArray(scip, &varbounds, SCIPexprtreeGetNVars(consdata->exprtrees[i])) );
2167  varboundssize = SCIPexprtreeGetNVars(consdata->exprtrees[i]);
2168  }
2169  else if( varboundssize < SCIPexprtreeGetNVars(consdata->exprtrees[i]) )
2170  {
2171  SCIP_CALL( SCIPreallocBufferArray(scip, &varbounds, SCIPexprtreeGetNVars(consdata->exprtrees[i])) );
2172  varboundssize = SCIPexprtreeGetNVars(consdata->exprtrees[i]);
2173  }
2174  assert(varbounds != NULL);
2175 
2176  for( j = 0; j < SCIPexprtreeGetNVars(consdata->exprtrees[i]); ++j )
2177  {
2178  var = SCIPexprtreeGetVars(consdata->exprtrees[i])[j];
2179  SCIPintervalSetBounds(&varbounds[j],
2180  -infty2infty(SCIPinfinity(scip), INTERVALINFTY, -MIN(SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var))), /*lint !e666*/
2181  +infty2infty(SCIPinfinity(scip), INTERVALINFTY, MAX(SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var))) ); /*lint !e666*/
2182  }
2183 
2184  SCIP_CALL( SCIPexprtreeCheckCurvature(consdata->exprtrees[i], INTERVALINFTY, varbounds, &consdata->curvatures[i], NULL) );
2185  consdata->curvatures[i] = SCIPexprcurvMultiply(consdata->nonlincoefs[i], consdata->curvatures[i]);
2186 
2187  if( consdata->curvatures[i] == SCIP_EXPRCURV_UNKNOWN && SCIPconshdlrGetData(SCIPconsGetHdlr(cons))->isreformulated && SCIPexprGetOperator(SCIPexprtreeGetRoot(consdata->exprtrees[i])) != SCIP_EXPR_USER )
2188  {
2189  SCIPverbMessage(scip, SCIP_VERBLEVEL_NORMAL, NULL, "indefinite expression tree in constraint <%s>\n", SCIPconsGetName(cons));
2190  SCIPdebug( SCIP_CALL( SCIPexprtreePrintWithNames(consdata->exprtrees[i], SCIPgetMessagehdlr(scip), NULL) ) );
2191  SCIPdebugMsgPrint(scip, "\n");
2192  }
2193  }
2194 
2195  /* @todo implement some more expensive checks */
2196 
2197  consdata->curvature = SCIPexprcurvAdd(consdata->curvature, consdata->curvatures[i]);
2198 
2199  SCIPdebugMsg(scip, "-> tree %d with coef %g is %s -> nonlinear part is %s\n", i, consdata->nonlincoefs[i], SCIPexprcurvGetName(consdata->curvatures[i]), SCIPexprcurvGetName(consdata->curvature));
2200  }
2201 
2202  SCIPfreeBufferArrayNull(scip, &varbounds);
2203 
2204  return SCIP_OKAY;
2205 } /*lint !e715*/
2206 
2207 /* replaces a node by another node in expression graph
2208  * moves all parents of node to replacement
2209  * replaces all exprgraphnode's in constraints that are node by replacement
2210  * node may be freed, if captured only by given constraints
2211  */
2212 static
2213 SCIP_RETCODE reformReplaceNode(
2214  SCIP_EXPRGRAPH* exprgraph, /**< expression graph */
2215  SCIP_EXPRGRAPHNODE** node, /**< pointer to node to be replaced in expression graph */
2216  SCIP_EXPRGRAPHNODE* replacement, /**< node which takes node's place */
2217  SCIP_CONS** conss, /**< constraints */
2218  int nconss /**< number of constraints */
2219  )
2220 {
2221  SCIP_CONSDATA* consdata;
2222  int c;
2223 
2224  assert(exprgraph != NULL);
2225  assert(node != NULL);
2226  assert(*node != NULL);
2227  assert(replacement != NULL);
2228  assert(conss != NULL || nconss == 0);
2229 
2230  SCIP_CALL( SCIPexprgraphMoveNodeParents(exprgraph, node, replacement) );
2231 
2232  /* node was not captured by any constraint */
2233  if( *node == NULL )
2234  return SCIP_OKAY;
2235 
2236  /* if node still exists, then because it is captured by some constraint (it should not have parents anymore)
2237  * thus, look into the given constraints and replace their exprgraphnode by replacement
2238  * @todo may be expensive when this is done more often, esp. when we know that node will not be freed due to an added auxiliary constraint
2239  */
2240  assert(*node == NULL || SCIPexprgraphGetNodeNParents(*node) == 0);
2241  for( c = 0; c < nconss; ++c )
2242  {
2243  assert(conss[c] != NULL); /*lint !e613*/
2244 
2245  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
2246  assert(consdata != NULL);
2247 
2248  if( consdata->exprgraphnode == *node )
2249  {
2250  SCIP_CALL( SCIPexprgraphReleaseNode(exprgraph, &consdata->exprgraphnode) );
2251  consdata->exprgraphnode = replacement;
2252  SCIPexprgraphCaptureNode(replacement);
2253 
2254  /* since we change the node, also the constraint changes, so ensure that it is presolved again */
2255  consdata->ispresolved = FALSE;
2256  }
2257  }
2258  *node = NULL;
2259 
2260  return SCIP_OKAY;
2261 }
2262 
2263 /** creates a new auxiliary variable and a new auxiliary nonlinear constraint connecting the var and a given node
2264  * node is replaced by new new auxiliary variables node in all parents of node in expression graph and in all constraints that use node
2265  */
2266 static
2267 SCIP_RETCODE reformNode2Var(
2268  SCIP* scip, /**< SCIP data structure */
2269  SCIP_EXPRGRAPH* exprgraph, /**< expression graph */
2270  SCIP_EXPRGRAPHNODE* node, /**< expression graph node */
2271  SCIP_CONS** conss, /**< constraints where to update exprgraphnode */
2272  int nconss, /**< number of constraints */
2273  int* naddcons, /**< counter to increase when constraint is added */
2274  SCIP_Bool donotmultaggr /**< whether to mark auxiliary variable as not to multiaggregate */
2275  )
2276 {
2277  char name[SCIP_MAXSTRLEN];
2278  SCIP_VAR* auxvar;
2279  SCIP_CONS* auxcons;
2280  SCIP_EXPRGRAPHNODE* auxvarnode;
2281  SCIP_INTERVAL bounds;
2282  SCIP_Real minusone;
2283  SCIP_Bool cutoff;
2284 
2285  assert(scip != NULL);
2286  assert(exprgraph != NULL);
2287  assert(node != NULL);
2288  assert(naddcons != NULL);
2289  assert(SCIPexprgraphGetNodeDepth(node) >= 1); /* do not turn vars or consts into new vars */
2290 
2291  bounds = SCIPexprgraphGetNodeBounds(node);
2292  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%d", *naddcons);
2293 
2294  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s for node %p(%d,%d)\n", name, (void*)node, SCIPexprgraphGetNodeDepth(node), SCIPexprgraphGetNodePosition(node));
2295 
2296  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, SCIPintervalGetInf(bounds), SCIPintervalGetSup(bounds), 0.0,
2297  SCIP_VARTYPE_CONTINUOUS, TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
2298  SCIP_CALL( SCIPaddVar(scip, auxvar) );
2299  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&auxvar, &auxvarnode) );
2300 #ifdef WITH_DEBUG_SOLUTION
2301  if( SCIPdebugIsMainscip(scip) )
2302  {
2303  /* store debug sol value of node as value for auxvar in debug solution and as value for auxvarnode */
2304  SCIPexprgraphSetVarNodeValue(auxvarnode, SCIPexprgraphGetNodeVal(node));
2305  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SCIPexprgraphGetNodeVal(node)) );
2306  }
2307 #endif
2308 
2309  if( donotmultaggr )
2310  {
2311  SCIP_CALL( SCIPmarkDoNotMultaggrVar(scip, auxvar) );
2312  }
2313 
2314  /* set also bounds of auxvarnode to bounds, so it is available for new parent nodes (currently node->parents)
2315  * when updating their curvature information; avoid having to run domain propagation through exprgraph
2316  */
2317  SCIPexprgraphTightenNodeBounds(exprgraph, auxvarnode, bounds, BOUNDTIGHTENING_MINSTRENGTH, INTERVALINFTY, &cutoff);
2318  assert(!cutoff); /* we tightenend bounds from [-inf,+inf] to bounds, this should not be infeasible */
2319 
2320  /* add new constraint auxvar == node */
2321  minusone = -1.0;
2322  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &auxcons, name, 1, &auxvar, &minusone, node, 0.0, 0.0, TRUE, TRUE, TRUE, TRUE, TRUE,
2323  FALSE, FALSE, FALSE, FALSE, FALSE) );
2324  SCIP_CALL( SCIPaddCons(scip, auxcons) );
2325 
2326  /* move parents of node in expression graph to auxvarnode
2327  * replace node by auxvarnode in constraints that use node */
2328  SCIP_CALL( reformReplaceNode(exprgraph, &node, auxvarnode, conss, nconss) );
2329 
2330  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
2331  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
2332 
2333  ++*naddcons;
2334 
2335  return SCIP_OKAY;
2336 }
2337 
2338 /** ensures that all children of a node have at least a given curvature by adding auxiliary variables */
2339 static
2340 SCIP_RETCODE reformEnsureChildrenMinCurvature(
2341  SCIP* scip, /**< SCIP data structure */
2342  SCIP_EXPRGRAPH* exprgraph, /**< expression graph */
2343  SCIP_EXPRGRAPHNODE* node, /**< expression graph node */
2344  SCIP_EXPRCURV mincurv, /**< minimal desired curvature */
2345  SCIP_CONS** conss, /**< constraints to check whether they use one of the replaced nodes */
2346  int nconss, /**< number of constraints to check */
2347  int* naddcons /**< counter to increase when constraint is added */
2348  )
2349 {
2350  SCIP_EXPRGRAPHNODE* child;
2351  SCIP_Bool needupdate;
2352 
2353  int i;
2354  assert(scip != NULL);
2355  assert(exprgraph != NULL);
2356  assert(node != NULL);
2357  assert(naddcons != NULL);
2358  assert(SCIPexprgraphGetNodeDepth(node) >= 1); /* do not turn vars or consts into new vars */
2359  assert(mincurv != SCIP_EXPRCURV_UNKNOWN); /* this is trivial and makes no sense */
2360 
2361  needupdate = FALSE; /* whether we need to update curvature of node */
2362 
2363  for( i = 0; i < SCIPexprgraphGetNodeNChildren(node); ++i )
2364  {
2365  child = SCIPexprgraphGetNodeChildren(node)[i];
2366  assert(child != NULL);
2367 
2368  if( (SCIPexprgraphGetNodeCurvature(child) & mincurv) != mincurv )
2369  {
2370  SCIPdebugMsg(scip, "add auxiliary variable for child %p(%d,%d) with curvature %s\n",
2371  (void*)child, SCIPexprgraphGetNodeDepth(child), SCIPexprgraphGetNodePosition(child), SCIPexprcurvGetName(SCIPexprgraphGetNodeCurvature(child)) );
2372 
2373  SCIP_CALL( reformNode2Var(scip, exprgraph, child, conss, nconss, naddcons, FALSE) );
2374  needupdate = TRUE;
2375 
2376  /* i'th child of node should now be a variable */
2377  assert(SCIPexprgraphGetNodeChildren(node)[i] != child);
2378  assert(SCIPexprgraphGetNodeOperator(SCIPexprgraphGetNodeChildren(node)[i]) == SCIP_EXPR_VARIDX);
2379  }
2380 
2381  assert(SCIPexprgraphGetNodeCurvature(SCIPexprgraphGetNodeChildren(node)[i]) & mincurv);
2382  }
2383 
2384  if( needupdate )
2385  {
2386  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(node, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
2387  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(node)));
2388  }
2389 
2390  return SCIP_OKAY;
2391 }
2392 
2393 /** reformulates a monomial by adding auxiliary variables and constraints for bilinear terms */
2394 static
2395 SCIP_RETCODE reformMonomial(
2396  SCIP* scip, /**< SCIP data structure */
2397  SCIP_EXPRGRAPH* exprgraph, /**< expression graph */
2398  int nfactors, /**< number of factors */
2399  SCIP_EXPRGRAPHNODE** factors, /**< factors */
2400  SCIP_Real* exponents, /**< exponents, or NULL if all 1.0 */
2401  SCIP_EXPRGRAPHNODE** resultnode, /**< buffer to store node which represents the reformulated monomial */
2402  SCIP_Bool createauxcons, /**< whether to create auxiliary var/cons */
2403  int mindepth, /**< minimal depth of new nodes in expression graph, or -1 */
2404  int* naddcons /**< buffer to increase by number of added cons */
2405  )
2406 {
2407  char name[SCIP_MAXSTRLEN];
2408  SCIP_VAR* auxvar;
2409  SCIP_CONS* auxcons;
2410  SCIP_Real minusone;
2411 
2412  assert(scip != NULL);
2413  assert(exprgraph != NULL);
2414  assert(nfactors > 0);
2415  assert(factors != NULL);
2416  assert(resultnode != NULL);
2417  assert(naddcons != NULL);
2418 
2419  /* factors are just one node */
2420  if( nfactors == 1 && (exponents == NULL || exponents[0] == 1.0) )
2421  {
2422  *resultnode = factors[0];
2423  return SCIP_OKAY;
2424  }
2425 
2426  /* only one factor, but with exponent < 0.0 and factor has mixed sign, e.g., x^(-3)
2427  * reformulate as auxvar * factor^(-exponent) = 1 and return the node for auxvar in resultnode
2428  */
2429  if( nfactors == 1 && exponents[0] < 0.0 && SCIPexprgraphGetNodeBounds(factors[0]).inf < 0.0 && SCIPexprgraphGetNodeBounds(factors[0]).sup > 0.0 ) /*lint !e613*/
2430  {
2431  SCIP_EXPRGRAPHNODE* auxnode;
2432  SCIP_EXPRGRAPHNODE* reformfactors[2];
2433  SCIP_Real reformexp[2];
2434 
2435  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%d", *naddcons);
2436  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s\n", name);
2437 
2438  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2439  SCIP_VARTYPE_CONTINUOUS, TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
2440  SCIP_CALL( SCIPaddVar(scip, auxvar) );
2441  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&auxvar, resultnode) );
2442 
2443 #ifdef WITH_DEBUG_SOLUTION
2444  /* store debug sol value of node as value for auxvar in debug solution and as value for resultnode */
2445  if( SCIPdebugIsMainscip(scip) )
2446  {
2447  SCIP_Real debugval;
2448  debugval = pow(SCIPexprgraphGetNodeVal(factors[0]), exponents[0]);
2449  SCIPexprgraphSetVarNodeValue(*resultnode, debugval);
2450  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
2451  }
2452 #endif
2453 
2454  /* increase naddcons before next call to reformMonomial, to avoid duplicate variable names, which is bad for debugging */
2455  ++*naddcons;
2456 
2457  /* add reformulation for resultnode(=auxvar) * factor^(-exponent) = 1.0
2458  * if exponent != -1.0, then factor^(-exponent) should be moved into extra variable
2459  * finally one should get an EXPR_MUL node */
2460  reformfactors[0] = *resultnode;
2461  reformfactors[1] = factors[0];
2462  reformexp[0] = 1.0;
2463  reformexp[1] = -exponents[0]; /*lint !e613*/
2464  SCIP_CALL( reformMonomial(scip, exprgraph, 2, reformfactors, reformexp, &auxnode, FALSE, mindepth, naddcons) );
2465 
2466  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &auxcons, name, 0, NULL, NULL, auxnode, 1.0, 1.0,
2467  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
2468  SCIP_CALL( SCIPaddCons(scip, auxcons) );
2469 
2470  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
2471  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
2472 
2473  return SCIP_OKAY;
2474  }
2475 
2476  /* only one factor, but with exponent != 1.0 */
2477  if( nfactors == 1 )
2478  {
2479  /* create some power expression node, if not existing already */
2480  SCIP_EXPRGRAPHNODE* expnode;
2481  SCIP_EXPRGRAPHNODE* parent;
2482  int p;
2483 
2484  assert(exponents != NULL);
2485 
2486  /* check if there is already a node for factors[0]^exponents[0] */
2487  expnode = NULL;
2488  for( p = 0; p < SCIPexprgraphGetNodeNParents(factors[0]); ++p)
2489  {
2490  parent = SCIPexprgraphGetNodeParents(factors[0])[p];
2491  if( SCIPisIntegral(scip, exponents[0]) &&
2492  SCIPexprgraphGetNodeOperator(parent) == SCIP_EXPR_INTPOWER &&
2493  SCIPexprgraphGetNodeIntPowerExponent(parent) == (int)SCIPround(scip, exponents[0]) )
2494  {
2495  expnode = parent;
2496  break;
2497  }
2498  if( SCIPexprgraphGetNodeOperator(parent) == SCIP_EXPR_REALPOWER &&
2499  SCIPisEQ(scip, SCIPexprgraphGetNodeRealPowerExponent(parent), exponents[0]) )
2500  {
2501  expnode = parent;
2502  }
2503  }
2504  if( expnode == NULL )
2505  {
2506  if( SCIPisIntegral(scip, exponents[0]) )
2507  SCIP_CALL( SCIPexprgraphCreateNode(SCIPblkmem(scip), &expnode, SCIP_EXPR_INTPOWER, (int)SCIPround(scip, exponents[0])) );
2508  else
2509  SCIP_CALL( SCIPexprgraphCreateNode(SCIPblkmem(scip), &expnode, SCIP_EXPR_REALPOWER, exponents[0]) );
2510 
2511  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, expnode, mindepth, 1, &factors[0]) );
2512  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(expnode, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
2513  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(expnode)));
2514  }
2515 
2516  if( createauxcons )
2517  {
2518  /* @todo if there was already a node for factors[0]^exponents[0], then there may have been also been already an auxiliary variable and constraint (-> ex7_3_4) */
2519  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%d", *naddcons);
2520  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s\n", name);
2521 
2522  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2523  SCIP_VARTYPE_CONTINUOUS, TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
2524  SCIP_CALL( SCIPaddVar(scip, auxvar) );
2525  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&auxvar, resultnode) );
2526 
2527 #ifdef WITH_DEBUG_SOLUTION
2528  if( SCIPdebugIsMainscip(scip) )
2529  {
2530  SCIPexprgraphSetVarNodeValue(*resultnode, SCIPexprgraphGetNodeVal(expnode));
2531  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SCIPexprgraphGetNodeVal(expnode)) );
2532  }
2533 #endif
2534 
2535  /* add new constraint resultnode(=auxvar) = expnode */
2536  minusone = -1.0;
2537  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &auxcons, name, 1, &auxvar, &minusone, expnode, 0.0, 0.0, TRUE, TRUE, TRUE, TRUE, TRUE,
2538  FALSE, FALSE, FALSE, FALSE, FALSE) );
2539  SCIP_CALL( SCIPaddCons(scip, auxcons) );
2540 
2541  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
2542  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
2543 
2544  ++*naddcons;
2545  }
2546  else
2547  {
2548  *resultnode = expnode;
2549  }
2550 
2551  return SCIP_OKAY;
2552  }
2553 
2554  if( nfactors == 2 && exponents != NULL && exponents[0] != 1.0 && exponents[0] == exponents[1] ) /*lint !e777*/
2555  {
2556  /* factor0^exponent * factor1^exponent with exponent != 1.0, reform as (factor0*factor1)^exponent */
2557  SCIP_EXPRGRAPHNODE* productnode;
2558 
2559  /* create node for factor0*factor1 */
2560  SCIP_CALL( reformMonomial(scip, exprgraph, 2, factors, NULL, &productnode, TRUE, mindepth, naddcons) );
2561 
2562  /* create node for productnode^exponents[0] by just calling this method again */
2563  SCIP_CALL( reformMonomial(scip, exprgraph, 1, &productnode, &exponents[0], resultnode, createauxcons, mindepth, naddcons) );
2564 
2565  return SCIP_OKAY;
2566  }
2567 
2568  if( nfactors == 2 && exponents != NULL && exponents[0] == -exponents[1] ) /*lint !e777*/
2569  {
2570  /* factor0^exponent * factor1^(-exponent), reform as (factor0/factor1)^exponent or (factor1/factor0)^(-exponent) */
2571  SCIP_EXPRGRAPHNODE* auxvarnode;
2572  SCIP_EXPRGRAPHNODE* auxconsnode;
2573  SCIP_EXPRGRAPHNODE* leftright[2];
2574  SCIP_Real absexp;
2575 
2576  /* create variable and constraint for factor0 = auxvar * factor1 (if exponent > 0) or factor1 = auxvar * factor0 (if exponent < 0) */
2577 
2578  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%d", *naddcons);
2579  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s\n", name);
2580 
2581  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2582  SCIP_VARTYPE_CONTINUOUS, TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
2583  SCIP_CALL( SCIPaddVar(scip, auxvar) );
2584  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&auxvar, &auxvarnode) );
2585 
2586 #ifdef WITH_DEBUG_SOLUTION
2587  /* store debug sol value of node as value for auxvar in debug solution and as value for resultnode */
2588  if( SCIPdebugIsMainscip(scip) )
2589  {
2590  SCIP_Real debugval;
2591  if( exponents[0] > 0.0 )
2592  debugval = SCIPexprgraphGetNodeVal(factors[0]) / SCIPexprgraphGetNodeVal(factors[1]);
2593  else
2594  debugval = SCIPexprgraphGetNodeVal(factors[1]) / SCIPexprgraphGetNodeVal(factors[0]);
2595  SCIPexprgraphSetVarNodeValue(auxvarnode, debugval);
2596  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
2597  }
2598 #endif
2599 
2600  /* add new constraint resultnode(= auxvar) * factor1 - factor0 == 0 (exponent > 0) or auxvar * factor0 - factor1 == 0 (exponent < 0) */
2601  leftright[0] = auxvarnode;
2602  leftright[1] = exponents[0] > 0.0 ? factors[1] : factors[0];
2603 
2604  SCIP_CALL( SCIPexprgraphCreateNode(SCIPblkmem(scip), &auxconsnode, SCIP_EXPR_MUL, NULL) );
2605  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, auxconsnode, -1, 2, leftright) );
2606 
2607  leftright[0] = auxconsnode;
2608  leftright[1] = exponents[0] > 0.0 ? factors[0] : factors[1];
2609 
2610  SCIP_CALL( SCIPexprgraphCreateNode(SCIPblkmem(scip), &auxconsnode, SCIP_EXPR_MINUS, NULL) );
2611  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, auxconsnode, -1, 2, leftright) );
2612 
2613  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &auxcons, name, 0, NULL, NULL, auxconsnode, 0.0, 0.0,
2614  TRUE, TRUE, TRUE, TRUE, TRUE,
2615  FALSE, FALSE, FALSE, FALSE, FALSE) );
2616  SCIP_CALL( SCIPaddCons(scip, auxcons) );
2617 
2618  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
2619  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
2620 
2621  ++*naddcons;
2622 
2623  /* create node for auxvarnode^abs(exponents[0]) by just calling this method again */
2624  absexp = fabs(exponents[0]);
2625  SCIP_CALL( reformMonomial(scip, exprgraph, 1, &auxvarnode, &absexp, resultnode, createauxcons, mindepth, naddcons) );
2626 
2627  return SCIP_OKAY;
2628  }
2629 
2630  /* @todo if nfactors > 2, assemble groups of factors with same exponent and replace these by a single variable first */
2631 
2632  {
2633  /* at least two factors */
2634  /* create auxvar for left half (recursively) and auxvar for right half (recursively) and maybe new auxvar for product */
2635  /* @todo it may be enough to replace single factors in a monomial to get it convex or concave, see Westerlund et.al. */
2636 
2637  SCIP_EXPRGRAPHNODE* productnode;
2638  SCIP_EXPRGRAPHNODE* leftright[2]; /* {left, right} */
2639  SCIP_EXPRGRAPHNODE* parent;
2640  int half;
2641  int p;
2642 
2643  half = nfactors / 2;
2644  assert(half > 0);
2645  assert(half < nfactors);
2646 
2647  SCIP_CALL( reformMonomial(scip, exprgraph, half, factors, exponents, &leftright[0], TRUE, mindepth, naddcons) );
2648  SCIP_CALL( reformMonomial(scip, exprgraph, nfactors-half, &factors[half], exponents != NULL ? &exponents[half] : NULL, &leftright[1], TRUE, mindepth, naddcons) ); /*lint !e826*/
2649 
2650  /* check if there is already a node for left * right */
2651  productnode = NULL;
2652  for( p = 0; p < SCIPexprgraphGetNodeNParents(leftright[0]); ++p)
2653  {
2654  parent = SCIPexprgraphGetNodeParents(leftright[0])[p];
2655  if( SCIPexprgraphGetNodeOperator(parent) != SCIP_EXPR_MUL )
2656  continue;
2657 
2658  assert(SCIPexprgraphGetNodeNChildren(parent) == 2);
2659  if( (SCIPexprgraphGetNodeChildren(parent)[0] == leftright[0] && SCIPexprgraphGetNodeChildren(parent)[1] == leftright[1]) ||
2660  ( SCIPexprgraphGetNodeChildren(parent)[0] == leftright[1] && SCIPexprgraphGetNodeChildren(parent)[1] == leftright[0]) )
2661  {
2662  productnode = parent;
2663  break;
2664  }
2665  }
2666  if( productnode == NULL )
2667  {
2668  /* create node for left * right */
2669  SCIP_CALL( SCIPexprgraphCreateNode(SCIPblkmem(scip), &productnode, SCIP_EXPR_MUL, NULL) );
2670  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, productnode, mindepth, 2, leftright) );
2671  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(productnode, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
2672  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(productnode)));
2673  }
2674 
2675  if( createauxcons )
2676  {
2677  /* @todo if there was already a node for factors[0]^exponents[0], then there may have been also been already an auxiliary variable and constraint (-> ex7_3_4) */
2678  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%d", *naddcons);
2679  SCIPdebugMsg(scip, "add auxiliary variable and constraint %s\n", name);
2680 
2681  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0,
2682  SCIP_VARTYPE_CONTINUOUS, TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
2683  SCIP_CALL( SCIPaddVar(scip, auxvar) );
2684  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&auxvar, resultnode) );
2685 
2686 #ifdef WITH_DEBUG_SOLUTION
2687  if( SCIPdebugIsMainscip(scip) )
2688  {
2689  SCIPexprgraphSetVarNodeValue(*resultnode, SCIPexprgraphGetNodeVal(productnode));
2690  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, SCIPexprgraphGetNodeVal(productnode)) );
2691  }
2692 #endif
2693 
2694  /* add new constraint resultnode(= auxvar) == left * right */
2695  minusone = -1.0;
2696  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &auxcons, name, 1, &auxvar, &minusone, productnode, 0.0, 0.0, TRUE, TRUE, TRUE, TRUE, TRUE,
2697  FALSE, FALSE, FALSE, FALSE, FALSE) );
2698  SCIP_CALL( SCIPaddCons(scip, auxcons) );
2699 
2700  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
2701  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
2702 
2703  ++*naddcons;
2704  }
2705  else
2706  {
2707  *resultnode = productnode;
2708  }
2709  }
2710 
2711  return SCIP_OKAY;
2712 }
2713 
2714 /** reformulates expression graph into a form so that for each node under- and overestimators could be computed
2715  * similar to factorable reformulation in other global solvers, but sometimes does not split up complex operands (like quadratic)
2716  */
2717 static
2718 SCIP_RETCODE reformulate(
2719  SCIP* scip, /**< SCIP data structure */
2720  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2721  SCIP_CONS** conss, /**< constraints */
2722  int nconss, /**< number of constraints */
2723  int* naddcons /**< buffer to increase by the number of added constraints */
2724  )
2725 {
2726  SCIP_CONSHDLRDATA* conshdlrdata;
2727  SCIP_CONSDATA* consdata;
2728  SCIP_EXPRGRAPH* exprgraph;
2729  SCIP_EXPRGRAPHNODE* node;
2730  SCIP_EXPRGRAPHNODE** children;
2731  SCIP_EXPRGRAPHNODE* reformnode;
2732  SCIP_Bool havenonlinparent;
2733  SCIP_Bool domainerror;
2734  int nchildren;
2735  int c;
2736  int d;
2737  int i;
2738  int u;
2739 #ifndef NDEBUG
2740  int j;
2741 #endif
2742 
2743  assert(scip != NULL);
2744  assert(conshdlr != NULL);
2745  assert(conss != NULL || nconss == 0);
2746  assert(naddcons != NULL);
2747  assert(SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING);
2748  assert(!SCIPinProbing(scip));
2749 
2750  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2751  assert(conshdlrdata != NULL);
2752 
2753  if( conshdlrdata->isreformulated )
2754  {
2755  SCIPdebugMsg(scip, "skip reformulation, already done\n");
2756  return SCIP_OKAY;
2757  }
2758 
2759  exprgraph = conshdlrdata->exprgraph;
2760 
2761  /* make sure current variable bounds are variable nodes of exprgraph */
2762  SCIP_CALL( SCIPexprgraphPropagateVarBounds(exprgraph, INTERVALINFTY, FALSE, &domainerror) );
2763  assert(!domainerror); /* should have been found by domain propagation */
2764 
2765  /* set debug solution in expression graph and evaluate nodes, so we can compute debug solution values for auxiliary variables */
2766 #ifdef WITH_DEBUG_SOLUTION
2767  if( SCIPdebugIsMainscip(scip) )
2768  {
2769  SCIP_Real* varvals;
2770 
2771  SCIP_CALL( SCIPallocBufferArray(scip, &varvals, SCIPexprgraphGetNVars(exprgraph)) );
2772 
2773  for( i = 0; i < SCIPexprgraphGetNVars(exprgraph); ++i )
2774  SCIP_CALL( SCIPdebugGetSolVal(scip, (SCIP_VAR*)SCIPexprgraphGetVars(exprgraph)[i], &varvals[i]) );
2775 
2776  SCIP_CALL( SCIPexprgraphEval(exprgraph, varvals) );
2777 
2778  SCIPfreeBufferArray(scip, &varvals);
2779  }
2780 #endif
2781 
2782  for( d = 1; d < SCIPexprgraphGetDepth(exprgraph); ++d )
2783  {
2784  i = 0;
2785  while( i < SCIPexprgraphGetNNodes(exprgraph)[d] )
2786  {
2787  node = SCIPexprgraphGetNodes(exprgraph)[d][i];
2788  assert(node != NULL);
2789 
2790  /* skip disabled nodes, they should be removed soon */
2791  if( !SCIPexprgraphIsNodeEnabled(node) )
2792  {
2793  ++i;
2794  continue;
2795  }
2796 
2797  /* make sure bounds and curvature of node are uptodate */
2798  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(node, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
2799  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(node)));
2800 
2801  /* try external reformulation methods */
2802  for( u = 0; u < conshdlrdata->nnlconsupgrades; ++u )
2803  {
2804  if( conshdlrdata->nlconsupgrades[u]->nodereform != NULL && conshdlrdata->nlconsupgrades[u]->active )
2805  {
2806  SCIP_CALL( conshdlrdata->nlconsupgrades[u]->nodereform(scip, exprgraph, node, naddcons, &reformnode) );
2807  if( reformnode == NULL )
2808  continue;
2809 
2810  SCIPdebugMsg(scip, "external nodereform reformulated node %p(%d,%d), replace by %p\n",
2811  (void*)node, SCIPexprgraphGetNodeDepth(node), SCIPexprgraphGetNodePosition(node), (void*)reformnode);
2812 
2813  SCIP_CALL( reformReplaceNode(exprgraph, &node, reformnode, conss, nconss) );
2814  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(reformnode, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
2815  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(reformnode)));
2816 
2817  break;
2818  }
2819  }
2820  /* if node has been reformulated, continue with next node without increasing i */
2821  if( u < conshdlrdata->nnlconsupgrades )
2822  continue;
2823 
2824  /* leave nodes that are known to be convex/concave/linear as they are */
2825  if( SCIPexprgraphGetNodeCurvature(node) != SCIP_EXPRCURV_UNKNOWN )
2826  {
2827  SCIPdebugMsg(scip, "skip reformulating node %p(%d,%d) = ", (void*)node, SCIPexprgraphGetNodeDepth(node), SCIPexprgraphGetNodePosition(node));
2828  SCIPdebug( SCIPexprgraphPrintNode(node, SCIPgetMessagehdlr(scip), NULL) );
2829  SCIPdebugMsgPrint(scip, ", curv = %s\n", SCIPexprcurvGetName(SCIPexprgraphGetNodeCurvature(node)));
2830  ++i;
2831  continue;
2832  }
2833 
2834  /* get flag whether node has a nonlinear parent
2835  * we want to know whether the current node will be at the top of the tree after the next simplification run
2836  * due to the tricky reformulation of polynomials below, this may not be the case yet
2837  */
2838  havenonlinparent = SCIPexprgraphHasNodeNonlinearAncestor(node);
2839 
2840  /* take action */
2841  assert(SCIPexprgraphGetNodeCurvature(node) == SCIP_EXPRCURV_UNKNOWN);
2842  SCIPdebugMsg(scip, "think about reformulating %s node %p(%d,%d) = ", SCIPexpropGetName(SCIPexprgraphGetNodeOperator(node)), (void*)node, SCIPexprgraphGetNodeDepth(node), SCIPexprgraphGetNodePosition(node));
2843  SCIPdebug( SCIPexprgraphPrintNode(node, SCIPgetMessagehdlr(scip), NULL) );
2844  SCIPdebugMsgPrint(scip, "\n");
2845 
2846  children = SCIPexprgraphGetNodeChildren(node);
2847  nchildren = SCIPexprgraphGetNodeNChildren(node);
2848  assert(children != NULL || nchildren == 0);
2849 
2850 #ifndef NDEBUG
2851  /* at this place, all children nodes should have a known curvature, except if they only appear only linearly in constraints
2852  * the latter for cases where constraints with unknown curvature are upgraded to other constraint handler that can handle these (quadratic, signpower,...)
2853  */
2854  for( j = 0; j < nchildren; ++j )
2855  {
2856  assert(children[j] != NULL); /*lint !e613*/
2857  if( havenonlinparent ||
2858  ( SCIPexprgraphGetNodeOperator(node) != SCIP_EXPR_PLUS &&
2859  SCIPexprgraphGetNodeOperator(node) != SCIP_EXPR_MINUS &&
2860  SCIPexprgraphGetNodeOperator(node) != SCIP_EXPR_SUM &&
2861  SCIPexprgraphGetNodeOperator(node) != SCIP_EXPR_LINEAR) )
2862  assert(SCIPexprgraphGetNodeCurvature(children[j]) != SCIP_EXPRCURV_UNKNOWN || SCIPexprgraphGetNodeOperator(children[j]) == SCIP_EXPR_USER); /*lint !e613*/
2863  }
2864 #endif
2865 
2866  switch( SCIPexprgraphGetNodeOperator(node) )
2867  {
2868  case SCIP_EXPR_VARIDX:
2869  case SCIP_EXPR_PARAM:
2870  case SCIP_EXPR_CONST:
2871  SCIPerrorMessage("node with operator %d cannot have unknown curvature\n", SCIPexprgraphGetNodeOperator(node));
2872  SCIPABORT();
2873  break; /*lint !e527*/
2874 
2875  /* linear operands */
2876  case SCIP_EXPR_PLUS:
2877  case SCIP_EXPR_MINUS:
2878  case SCIP_EXPR_SUM:
2879  case SCIP_EXPR_LINEAR:
2880  /* children have conflicting curvature, we can handle such sums in cons_nonlinear
2881  * thus, turn node into variable, if it has nonlinear parents */
2882  if( havenonlinparent )
2883  {
2884  SCIP_CALL( reformNode2Var(scip, exprgraph, node, conss, nconss, naddcons, FALSE) );
2885  assert(node != NULL);
2886  assert(SCIPexprgraphGetNodeNParents(node) == 0); /* node should now be at top of graph */
2887  }
2888  ++i;
2889  break;
2890 
2891  /* quadratic operands */
2892  case SCIP_EXPR_MUL:
2893  case SCIP_EXPR_QUADRATIC:
2894  {
2895  SCIP_EXPRGRAPHNODE* child;
2896  SCIP_Bool needupdate;
2897 
2898  /* ensure all children are linear, so next simplifier run makes sure all children will be variables (by distributing the product)
2899  * however, that will not work for user-expressions, so we should also ensure that they are none (@todo as they are linear, they could actually be replaced by a regular linear expression)
2900  */
2901  SCIPdebugMessage("ensure children are linear\n");
2902  SCIP_CALL( reformEnsureChildrenMinCurvature(scip, exprgraph, node, SCIP_EXPRCURV_LINEAR, conss, nconss, naddcons) );
2903 
2904  needupdate = FALSE; /* whether we need to update curvature of node */
2905  for( c = 0; c < SCIPexprgraphGetNodeNChildren(node); ++c )
2906  {
2907  child = SCIPexprgraphGetNodeChildren(node)[c];
2908  assert(child != NULL);
2909 
2910  if( SCIPexprgraphGetNodeCurvature(child) != SCIP_EXPRCURV_LINEAR || SCIPexprgraphGetNodeOperator(child) == SCIP_EXPR_USER || SCIPexprgraphGetNodeOperator(child) == SCIP_EXPR_ABS )
2911  {
2912  SCIPdebugMessage("add auxiliary variable for child %p(%d,%d) with curvature %s operator %s\n",
2913  (void*)child, SCIPexprgraphGetNodeDepth(child), SCIPexprgraphGetNodePosition(child), SCIPexprcurvGetName(SCIPexprgraphGetNodeCurvature(child)), SCIPexpropGetName(SCIPexprgraphGetNodeOperator(child)) );
2914 
2915  SCIP_CALL( reformNode2Var(scip, exprgraph, child, conss, nconss, naddcons, FALSE) );
2916  needupdate = TRUE;
2917 
2918  /* c'th child of node should now be a variable */
2919  assert(SCIPexprgraphGetNodeChildren(node)[c] != child);
2920  assert(SCIPexprgraphGetNodeOperator(SCIPexprgraphGetNodeChildren(node)[c]) == SCIP_EXPR_VARIDX);
2921  }
2922  }
2923  if( needupdate )
2924  {
2925  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(node, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
2926  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(node)));
2927  }
2928 
2929  if( SCIPexprgraphGetNodeCurvature(node) != SCIP_EXPRCURV_UNKNOWN )
2930  {
2931  /* if curvature is now known then we are done */
2932  ++i;
2933  break;
2934  }
2935 
2936  /* if we have nonlinear parents or a sibling, then add auxiliary variable for this node, so an upgrade to cons_quadratic should take place
2937  * we assume that siblings are non-linear and non-quadratic, which should be the case if simplifier was run, and also if this node was created during reformulating a polynomial
2938  * @todo we could also add auxvars for the sibling nodes, e.g., if there is only one
2939  * @todo if sibling nodes are quadratic (or even linear) due to reformulation, then we do not need to reform here... (-> nvs16)
2940  * maybe this step should not be done here at all if havenonlinparent is FALSE? e.g., move into upgrade from quadratic?
2941  */
2942  if( havenonlinparent || SCIPexprgraphHasNodeSibling(node) )
2943  {
2944  SCIP_CALL( reformNode2Var(scip, exprgraph, node, conss, nconss, naddcons, FALSE) );
2945  assert(node != NULL);
2946  assert(SCIPexprgraphGetNodeNParents(node) == 0); /* node should now be at top of graph, so it can be upgraded by cons_quadratic */
2947  break;
2948  }
2949 
2950  ++i;
2951  break;
2952  }
2953 
2954  case SCIP_EXPR_DIV:
2955  {
2956  /* reformulate as bilinear term
2957  * note that in the reformulation, a zero in the denominator forces the nominator to be zero too, but the auxiliary can be arbitrary
2958  */
2959  SCIP_EXPRGRAPHNODE* auxvarnode;
2960  SCIP_EXPRGRAPHNODE* auxnode;
2961  SCIP_EXPRGRAPHNODE* auxchildren[3];
2962  SCIP_Real lincoefs[3];
2963  SCIP_QUADELEM quadelem;
2964  SCIP_VAR* auxvar;
2965  SCIP_CONS* auxcons;
2966  char name[SCIP_MAXSTRLEN];
2967  SCIP_INTERVAL bounds;
2968 
2969  assert(children != NULL);
2970  assert(nchildren == 2);
2971 
2972  bounds = SCIPexprgraphGetNodeBounds(node);
2973  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nlreform%d", *naddcons);
2974 
2975  SCIPdebugMsg(scip, "add auxiliary variable %s for division in node %p(%d,%d)\n", name, (void*)node, SCIPexprgraphGetNodeDepth(node), SCIPexprgraphGetNodePosition(node));
2976 
2977  SCIP_CALL( SCIPcreateVar(scip, &auxvar, name, SCIPintervalGetInf(bounds), SCIPintervalGetSup(bounds), 0.0,
2978  SCIP_VARTYPE_CONTINUOUS, TRUE, TRUE, NULL, NULL, NULL, NULL, NULL) );
2979  SCIP_CALL( SCIPaddVar(scip, auxvar) );
2980  SCIP_CALL( SCIPexprgraphAddVars(exprgraph, 1, (void**)&auxvar, &auxvarnode) );
2981 
2982 #ifdef WITH_DEBUG_SOLUTION
2983  if( SCIPdebugIsMainscip(scip) )
2984  {
2985  SCIP_Real debugval;
2986  debugval = SCIPexprgraphGetNodeVal(children[0]) / SCIPexprgraphGetNodeVal(children[1]);
2987  SCIPexprgraphSetVarNodeValue(auxvarnode, debugval);
2988  SCIP_CALL( SCIPdebugAddSolVal(scip, auxvar, debugval) );
2989  }
2990 #endif
2991 
2992  /* add new constraint auxvar * child[1] - child[0] == 0 */
2993  auxchildren[0] = children[0]; /*lint !e613*/
2994  auxchildren[1] = children[1]; /*lint !e613*/
2995  auxchildren[2] = auxvarnode;
2996 
2997  lincoefs[0] = -1.0;
2998  lincoefs[1] = 0.0;
2999  lincoefs[2] = 0.0;
3000 
3001  quadelem.idx1 = 1;
3002  quadelem.idx2 = 2;
3003  quadelem.coef = 1.0;
3004 
3005  SCIP_CALL( SCIPexprgraphCreateNodeQuadratic(SCIPblkmem(scip), &auxnode, 3, lincoefs, 1, &quadelem, 0.0) );
3006  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, auxnode, -1, 3, auxchildren) );
3007 
3008  SCIP_CALL( SCIPcreateConsNonlinear2(scip, &auxcons, name, 0, NULL, NULL, auxnode, 0.0, 0.0,
3009  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
3010  SCIP_CALL( SCIPaddCons(scip, auxcons) );
3011 
3012  /* replace node by auxvarnode in graph and constraints that use it */
3013  SCIP_CALL( reformReplaceNode(exprgraph, &node, auxvarnode, conss, nconss) );
3014 
3015  SCIP_CALL( SCIPreleaseCons(scip, &auxcons) );
3016  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
3017 
3018  ++*naddcons;
3019 
3020  /* do not increase i, since node was removed and not necessarily replaced here */
3021  break;
3022  }
3023 
3024  case SCIP_EXPR_MIN:
3025  {
3026  /* make sure that both children are concave, because min of concave functions is concave */
3027  SCIP_CALL( reformEnsureChildrenMinCurvature(scip, exprgraph, node, SCIP_EXPRCURV_CONCAVE, conss, nconss, naddcons) );
3028  assert(SCIPexprgraphGetNodeCurvature(node) & SCIP_EXPRCURV_CONCAVE);
3029  ++i;
3030  break;
3031  }
3032 
3033  case SCIP_EXPR_MAX:
3034  {
3035  /* make sure that both children are convex, because max of convex functions is convex */
3036  SCIP_CALL( reformEnsureChildrenMinCurvature(scip, exprgraph, node, SCIP_EXPRCURV_CONVEX, conss, nconss, naddcons) );
3037  assert(SCIPexprgraphGetNodeCurvature(node) & SCIP_EXPRCURV_CONVEX);
3038  ++i;
3039  break;
3040  }
3041 
3042  case SCIP_EXPR_INTPOWER:
3043  {
3044  assert(nchildren == 1);
3045 
3046  /* for an intpower with mixed sign in the base and negative exponent, we reformulate similar as for EXPR_DIV */
3047  if( SCIPexprgraphGetNodeIntPowerExponent(node) < 0 && SCIPintervalGetInf(SCIPexprgraphGetNodeBounds(children[0])) < 0.0 && SCIPintervalGetSup(SCIPexprgraphGetNodeBounds(children[0])) > 0.0 ) /*lint !e613*/
3048  {
3049  SCIP_EXPRGRAPHNODE* auxvarnode;
3050  SCIP_Real exponent;
3051 
3052  /* if we have something like x^(-3) with mixed sign for x, then add auxvar and reform as auxvar*x^3 = 1 via reformMonomial */
3053  exponent = (SCIP_Real)SCIPexprgraphGetNodeIntPowerExponent(node);
3054  SCIP_CALL( reformMonomial(scip, exprgraph, 1, children, &exponent, &auxvarnode, TRUE, SCIPexprgraphGetNodeDepth(node), naddcons) );
3055  /* replace node by auxvarnode */
3056  SCIP_CALL( reformReplaceNode(exprgraph, &node, auxvarnode, conss, nconss) );
3057  break;
3058  }
3059 
3060  /* otherwise, we continue as for other univariate operands */
3061  } /*lint -fallthrough*/
3062 
3063  /* univariate operands where the child does not have bounds and curvature from which we can deduce curvature of this node,
3064  * but where we can do more if the child is linear
3065  * thus, turn child into auxiliary variable
3066  */
3067  case SCIP_EXPR_SQUARE:
3068  case SCIP_EXPR_SQRT:
3069  case SCIP_EXPR_EXP:
3070  case SCIP_EXPR_LOG:
3071  case SCIP_EXPR_ABS:
3072  case SCIP_EXPR_REALPOWER:
3073  case SCIP_EXPR_SIGNPOWER:
3074  {
3075  assert(nchildren == 1);
3076 
3077  SCIP_CALL( reformEnsureChildrenMinCurvature(scip, exprgraph, node, SCIP_EXPRCURV_LINEAR, conss, nconss, naddcons) );
3078 
3079  if( SCIPexprgraphGetNodeCurvature(node) == SCIP_EXPRCURV_UNKNOWN )
3080  {
3081  /* the only case where f(x) for a linear term x is indefinite here is if f is intpower or signpower and x has mixed sign */
3082  assert(SCIPexprgraphGetNodeOperator(node) == SCIP_EXPR_INTPOWER || SCIPexprgraphGetNodeOperator(node) == SCIP_EXPR_SIGNPOWER);
3083  assert(SCIPintervalGetInf(SCIPexprgraphGetNodeBounds(children[0])) < 0.0); /*lint !e613*/
3084  assert(SCIPintervalGetSup(SCIPexprgraphGetNodeBounds(children[0])) > 0.0); /*lint !e613*/
3085  }
3086 
3087  /* update curvature of node */
3088  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(node, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
3089  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(node)));
3090 
3091  if( SCIPexprgraphGetNodeCurvature(node) == SCIP_EXPRCURV_UNKNOWN )
3092  {
3093  /* if intpower and signpower with positive exponent and a mixed sign in the child bounds still does not give a curvature,
3094  * we can do more if we make this node the root of a nonlinear constraints expression node, so it can be upgraded by cons_signpower
3095  * of course, this is only required if the node is still intermediate
3096  *
3097  * an intpower with negative exponent should have been handled above
3098  * for signpower, we assume the exponent is > 1.0
3099  */
3100  assert(SCIPexprgraphGetNodeOperator(node) == SCIP_EXPR_INTPOWER || SCIPexprgraphGetNodeOperator(node) == SCIP_EXPR_SIGNPOWER);
3101  assert(SCIPexprgraphGetNodeOperator(node) != SCIP_EXPR_INTPOWER || SCIPexprgraphGetNodeIntPowerExponent(node) > 1);
3102  assert(SCIPexprgraphGetNodeOperator(node) != SCIP_EXPR_SIGNPOWER || SCIPexprgraphGetNodeSignPowerExponent(node) > 1.0);
3103  if( havenonlinparent )
3104  {
3105  SCIP_CALL( reformNode2Var(scip, exprgraph, node, conss, nconss, naddcons, FALSE) );
3106  assert(node != NULL); /* it should be used by some auxiliary constraint now */
3107  assert(SCIPexprgraphGetNodeNParents(node) == 0); /* node should now be at top of graph (and used by new auxiliary constraint) */
3108  }
3109  }
3110  ++i;
3111 
3112  break;
3113  }
3114 
3115  case SCIP_EXPR_SIN:
3116  case SCIP_EXPR_COS:
3117  case SCIP_EXPR_TAN:
3118  case SCIP_EXPR_SIGN:
3119  /* case SCIP_EXPR_ERF : */
3120  /* case SCIP_EXPR_ERFI : */
3121  {
3122  SCIPerrorMessage("no support for trigonometric or sign operands yet\n");
3123  return SCIP_ERROR;
3124  }
3125 
3126  case SCIP_EXPR_PRODUCT:
3127  {
3128  /* ensure all children are linear */
3129  SCIP_CALL( reformEnsureChildrenMinCurvature(scip, exprgraph, node, SCIP_EXPRCURV_LINEAR, conss, nconss, naddcons) );
3130  if( SCIPexprgraphGetNodeCurvature(node) != SCIP_EXPRCURV_UNKNOWN )
3131  {
3132  ++i;
3133  break;
3134  }
3135 
3136  /* if curvature is still unknown (quite likely), then turn into a cascade of bilinear terms
3137  * if node has parents, then ensure that it has a known curvature, otherwise we are also fine with a node that is a product of two (aux)variables */
3138  SCIP_CALL( reformMonomial(scip, exprgraph, nchildren, children, NULL, &reformnode, havenonlinparent, SCIPexprgraphGetNodeDepth(node), naddcons) );
3139 
3140  /* replace node by reformnode in graph and in all constraints that use it */
3141  SCIP_CALL( reformReplaceNode(exprgraph, &node, reformnode, conss, nconss) );
3142 
3143  /* do not increase i, since node was removed and not necessarily replaced here */
3144  break;
3145  }
3146 
3147  case SCIP_EXPR_POLYNOMIAL:
3148  {
3149  /* if polynomial has several monomials, replace by a sum of nodes each having a single monomial and one that has all linear and quadratic monomials
3150  * if polynomial has only a single monomial, then reformulate that one
3151  */
3152  SCIP_EXPRDATA_MONOMIAL** monomials;
3153  SCIP_EXPRDATA_MONOMIAL* monomial;
3154  int nmonomials;
3155  SCIP_Real* exponents;
3156  SCIP_Real coef;
3157  int* childidxs;
3158  int nfactors;
3159  int f;
3160  SCIP_INTERVAL childbounds;
3161  SCIP_EXPRCURV childcurv;
3162  SCIP_Bool modified;
3163 
3164  monomials = SCIPexprgraphGetNodePolynomialMonomials(node);
3165  nmonomials = SCIPexprgraphGetNodePolynomialNMonomials(node);
3166  assert(nmonomials >= 1); /* constant polynomials should have been simplified away */
3167 
3168  if( nmonomials > 1 )
3169  {
3170  SCIP_EXPRGRAPHNODE* sumnode;
3171  SCIP_Real constant;
3172  int nquadelems;
3173  SCIP_QUADELEM* quadelems;
3174  SCIP_Real* lincoefs;
3175  int nmonomialnodes;
3176  SCIP_EXPRGRAPHNODE** childrennew;
3177  SCIP_EXPRGRAPHNODE** monomialnodes;
3178  SCIP_Bool foundlincoefs;
3179  int m;
3180 
3181  /* @todo if a monomial is a factor of another monomial, then we could (and should?) replace it there by the node we create for it here -> ex7_2_1
3182  * @todo factorizing the polynomial could be beneficial
3183  */
3184 
3185  /* constant part of polynomials, to add to first monomialnode, if any, or quadratic or linear part */
3186  constant = SCIPexprgraphGetNodePolynomialConstant(node);
3187 
3188  /* coefficients from linear monomials */
3189  foundlincoefs = FALSE;
3190 
3191  /* quadratic elements */
3192  nquadelems = 0;
3193 
3194  /* expression graph nodes representing single higher-degree monomials, and single node with linear and/or quadratic monomials */
3195  nmonomialnodes = 0;
3196  SCIP_CALL( SCIPallocBufferArray(scip, &monomialnodes, nmonomials) );
3197 
3198  /* allocate memory */
3199  SCIP_CALL( SCIPallocClearBufferArray(scip, &lincoefs, nchildren) );
3200  SCIP_CALL( SCIPallocBufferArray(scip, &quadelems, nmonomials) );
3201  SCIP_CALL( SCIPallocBufferArray(scip, &childrennew, nchildren) );
3202 
3203  for( m = 0; m < nmonomials; ++m )
3204  {
3205  monomial = monomials[m];
3206  assert(monomial != NULL);
3207 
3208  coef = SCIPexprGetMonomialCoef(monomial);
3209  exponents = SCIPexprGetMonomialExponents(monomial);
3210  childidxs = SCIPexprGetMonomialChildIndices(monomial);
3211  nfactors = SCIPexprGetMonomialNFactors(monomial);
3212  assert(nfactors >= 1); /* constant monomials should have been simplified away */
3213  assert(coef != 0.0); /* zero-monomials should have been simplified away */
3214 
3215  if( nfactors == 1 && exponents[0] == 1.0 )
3216  {
3217  /* linear monomial */
3218  foundlincoefs = TRUE;
3219  assert(0 <= childidxs[0] && childidxs[0] < nchildren);
3220  assert(lincoefs[childidxs[0]] == 0.0); /* monomials should have been merged */
3221  lincoefs[childidxs[0]] = coef;
3222  }
3223  else if( nfactors == 1 && exponents[0] == 2.0 )
3224  {
3225  /* square monomial */
3226  quadelems[nquadelems].idx1 = childidxs[0];
3227  quadelems[nquadelems].idx2 = childidxs[0];
3228  quadelems[nquadelems].coef = coef;
3229  ++nquadelems;
3230  }
3231  else if( nfactors == 2 && exponents[0] == 1.0 && exponents[1] == 1.0 )
3232  {
3233  /* bilinear monomial */
3234  if( childidxs[0] < childidxs[1] )
3235  {
3236  quadelems[nquadelems].idx1 = childidxs[0];
3237  quadelems[nquadelems].idx2 = childidxs[1];
3238  }
3239  else
3240  {
3241  quadelems[nquadelems].idx1 = childidxs[1];
3242  quadelems[nquadelems].idx2 = childidxs[0];
3243  }
3244  quadelems[nquadelems].coef = coef;
3245  ++nquadelems;
3246  }
3247  else
3248  {
3249  /* general monomial -> pass into separate expression graph node */
3250  SCIP_EXPRDATA_MONOMIAL* monomialnew;
3251 
3252  /* create new node for this monomial, children will be those associated with factors */
3253  SCIP_CALL( SCIPexprCreateMonomial(SCIPblkmem(scip), &monomialnew, coef, nfactors, NULL, exponents) );
3254  SCIP_CALL( SCIPexprgraphCreateNodePolynomial(SCIPblkmem(scip), &monomialnodes[nmonomialnodes], 1, &monomialnew, constant, FALSE) );
3255  constant = 0.0;
3256 
3257  assert(nfactors <= nchildren);
3258  for( f = 0; f < nfactors; ++f )
3259  childrennew[f] = children[childidxs[f]]; /*lint !e613*/
3260 
3261  /* add new node to same depth as this node, so we will reformulate it during this run
3262  * no need to refresh bounds/curvature here, since that will be done when we reach this node next */
3263  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, monomialnodes[nmonomialnodes], SCIPexprgraphGetNodeDepth(node), nfactors, childrennew) );
3264 
3265  ++nmonomialnodes;
3266  }
3267  }
3268  /* should have had at least one linear, quadratic, or general monomial */
3269  assert(foundlincoefs || nquadelems > 0 || nmonomialnodes > 0);
3270 
3271  if( nquadelems > 0 )
3272  {
3273  /* create and add additional node for quadratic and linear part, simplifier should take care of removing unused children later */
3274  SCIP_CALL( SCIPexprgraphCreateNodeQuadratic(SCIPblkmem(scip), &monomialnodes[nmonomialnodes], nchildren, foundlincoefs ? lincoefs : NULL, nquadelems, quadelems, constant) );
3275  constant = 0.0;
3276  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, monomialnodes[nmonomialnodes], SCIPexprgraphGetNodeDepth(node), nchildren, children) );
3277  ++nmonomialnodes;
3278  }
3279  else if( foundlincoefs )
3280  {
3281  /* create additional node for linear part, simplifier should take care of removing unused children later */
3282  SCIP_CALL( SCIPexprgraphCreateNodeLinear(SCIPblkmem(scip), &monomialnodes[nmonomialnodes], nchildren, lincoefs, constant) );
3283  constant = 0.0;
3284  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, monomialnodes[nmonomialnodes], SCIPexprgraphGetNodeDepth(node), nchildren, children) );
3285  ++nmonomialnodes;
3286  }
3287  assert(constant == 0.0); /* the constant should have been used somewhere */
3288 
3289  /* release memory */
3290  SCIPfreeBufferArray(scip, &childrennew);
3291  SCIPfreeBufferArray(scip, &quadelems);
3292  SCIPfreeBufferArray(scip, &lincoefs);
3293 
3294  assert(nmonomialnodes > 0);
3295  if( nmonomialnodes > 1 )
3296  {
3297  /* add node for sum of monomials to expression graph */
3298  SCIP_CALL( SCIPexprgraphCreateNode(SCIPblkmem(scip), &sumnode, nmonomialnodes == 2 ? SCIP_EXPR_PLUS : SCIP_EXPR_SUM) );
3299  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, sumnode, -1, nmonomialnodes, monomialnodes) );
3300  }
3301  else
3302  {
3303  /* if only one monomial, then because polynomial was linear or quadratic... */
3304  assert(SCIPexprgraphGetNodeOperator(monomialnodes[0]) == SCIP_EXPR_LINEAR || SCIPexprgraphGetNodeOperator(monomialnodes[0]) == SCIP_EXPR_QUADRATIC);
3305  sumnode = monomialnodes[0];
3306  }
3307  SCIPfreeBufferArray(scip, &monomialnodes);
3308 
3309  /* replace node by sumnode, and we are done */
3310  SCIP_CALL( reformReplaceNode(exprgraph, &node, sumnode, conss, nconss) );
3311 
3312  SCIPdebugMsg(scip, "splitup polynomial into sum of %d nodes\n", nmonomialnodes);
3313 
3314  break;
3315  }
3316 
3317  /* reformulate a monomial such that it becomes convex or concave, if necessary */
3318 
3319  monomial = monomials[0];
3320  assert(monomial != NULL);
3321 
3322  coef = SCIPexprGetMonomialCoef(monomial);
3323  exponents = SCIPexprGetMonomialExponents(monomial);
3324  childidxs = SCIPexprGetMonomialChildIndices(monomial);
3325  nfactors = SCIPexprGetMonomialNFactors(monomial);
3326  assert(nfactors >= 1); /* constant monomials should have been simplified away */
3327  assert(coef != 0.0); /* zero-monomials should have been simplified away */
3328  assert(children != NULL);
3329 
3330  /* check if we make monomial convex or concave by making a child linear */
3331  modified = FALSE;
3332  if( nfactors == 1 )
3333  {
3334  /* ensure that the child of an univariate monomial is linear if its current (bounds,curvature) yields an unknown curvature for the monomial
3335  * and with linear child it had a known curvature (rules out x^a, a negative, x not linear) */
3336  childcurv = SCIPexprgraphGetNodeCurvature(children[childidxs[0]]); /*lint !e613*/
3337  childbounds = SCIPexprgraphGetNodeBounds(children[childidxs[0]]); /*lint !e613*/
3338  assert(SCIPexprcurvPower(childbounds, childcurv, exponents[0]) == SCIP_EXPRCURV_UNKNOWN); /* this is exactly the curvature of the node, which is unknown */
3339 
3340  /* if monomial were convex or concave if child were linear, then make child linear */
3341  if( SCIPexprcurvPower(childbounds, SCIP_EXPRCURV_LINEAR, exponents[0]) != SCIP_EXPRCURV_UNKNOWN )
3342  {
3343  assert(childcurv != SCIP_EXPRCURV_LINEAR);
3344  SCIPdebugMsg(scip, "reform child %d (univar. monomial) with curv %s into var\n", childidxs[0], SCIPexprcurvGetName(childcurv));
3345  SCIP_CALL( reformNode2Var(scip, exprgraph, children[childidxs[0]], conss, nconss, naddcons, FALSE) ); /*lint !e613*/
3346  modified = TRUE;
3347  }
3348  }
3349  else
3350  {
3351  /* check if the conditions on the exponents allow for a convex or concave monomial, assuming that the children are linear
3352  * if one of these conditions is fulfilled but a children curvature does not fit, then make these children linear
3353  */
3354  int nnegative;
3355  int npositive;
3356  SCIP_Real sum;
3357  SCIP_Bool expcurvpos;
3358  SCIP_Bool expcurvneg;
3359  SCIP_EXPRCURV desiredcurv;
3360 
3361  nnegative = 0; /* number of negative exponents */
3362  npositive = 0; /* number of positive exponents */
3363  sum = 0.0; /* sum of exponents */
3364  expcurvpos = TRUE; /* whether exp_j * f_j''(x) >= 0 for all factors (assuming f_j >= 0) */
3365  expcurvneg = TRUE; /* whether exp_j * f_j''(x) <= 0 for all factors (assuming f_j >= 0) */
3366 
3367  /* ensure that none of the children have unknown curvature */
3368  for( c = 0; c < SCIPexprgraphGetNodeNChildren(node); ++c )
3369  {
3370  childcurv = SCIPexprgraphGetNodeCurvature(children[c]); /*lint !e613*/
3371  if( childcurv == SCIP_EXPRCURV_UNKNOWN )
3372  {
3373  SCIPdebugMessage("reform child %d with unknown curvature into var\n", c);
3374  SCIP_CALL( reformNode2Var(scip, exprgraph, children[c], conss, nconss, naddcons, FALSE) ); /*lint !e613*/
3375  modified = TRUE;
3376  }
3377  }
3378  if( modified )
3379  {
3380  /* refresh curvature information in node, since we changed children */
3381  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(node, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
3382  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(node)));
3383 
3384  modified = FALSE;
3385  }
3386 
3387  for( f = 0; f < nfactors; ++f )
3388  {
3389  childcurv = SCIPexprgraphGetNodeCurvature(children[childidxs[f]]); /*lint !e613*/
3390  assert(childcurv != SCIP_EXPRCURV_UNKNOWN);
3391  childbounds = SCIPexprgraphGetNodeBounds(children[childidxs[f]]); /*lint !e613*/
3392  if( childbounds.inf < 0.0 && childbounds.sup > 0.0 )
3393  break;
3394 
3395  if( exponents[f] < 0.0 )
3396  ++nnegative;
3397  else
3398  ++npositive;
3399  sum += exponents[f];
3400 
3401  /* negate curvature if factor is negative */
3402  if( childbounds.inf < 0.0 )
3403  childcurv = SCIPexprcurvNegate(childcurv);
3404 
3405  /* check if exp_j * checkcurv is convex (>= 0) and/or concave */
3406  childcurv = SCIPexprcurvMultiply(exponents[f], childcurv);
3407  if( !(childcurv & SCIP_EXPRCURV_CONVEX) )
3408  expcurvpos = FALSE;
3409  if( !(childcurv & SCIP_EXPRCURV_CONCAVE) )
3410  expcurvneg = FALSE;
3411  }
3412 
3413  /* if some child can be both positive and negative, then nothing we can do here to get the monomial convex or concave
3414  * otherwise (i.e., f == nfactors), look further */
3415  desiredcurv = SCIP_EXPRCURV_UNKNOWN;
3416  if( f == nfactors )
3417  {
3418  /* if all factors are linear, then a product f_j^exp_j with f_j >= 0 is convex if
3419  * - all exponents are negative, or
3420  * - all except one exponent j* are negative and exp_j* >= 1 - sum_{j!=j*}exp_j, but the latter is equivalent to sum_j exp_j >= 1
3421  * further, the product is concave if
3422  * - all exponents are positive and the sum of exponents is <= 1.0
3423  *
3424  * if factors are nonlinear, then we require additionally, that for convexity
3425  * - each factor is convex if exp_j >= 0, or concave if exp_j <= 0, i.e., exp_j*f_j'' >= 0
3426  * and for concavity, we require that
3427  * - all factors are concave, i.e., exp_j*f_j'' <= 0
3428  */
3429 
3430  if( nnegative == nfactors || (nnegative == nfactors-1 && SCIPisGE(scip, sum, 1.0)) )
3431  {
3432  /* if exponents are such that we can be convex, but children curvature does not fit, make some children linear */
3433  SCIPdebugMsg(scip, "%d-variate monomial is convex (modulo sign), child curv fits = %u\n", nfactors, expcurvpos);
3434  /* since current node curvature is set to unknown, there must be such a child, since otherwise the node curvature had to be convex */
3435  assert(!expcurvpos);
3436  desiredcurv = SCIP_EXPRCURV_CONVEX;
3437  }
3438  else if( npositive == nfactors && SCIPisLE(scip, sum, 1.0) )
3439  {
3440  /* if exponents are such that we can be concave, but children curvature does not fit, make some children linear */
3441  SCIPdebugMsg(scip, "%d-variate monomial is concave (modulo sign), child curv fits = %u\n", nfactors, expcurvneg);
3442  /* since current node curvature is set to unknown, there must be such a child, since otherwise the node curvature had to be concave */
3443  assert(!expcurvneg);
3444  desiredcurv = SCIP_EXPRCURV_CONCAVE;
3445  }
3446  else
3447  {
3448  /* exponents are such that monomial is neither convex nor concave even if children were linear
3449  * thus, reformulate monomial below
3450  */
3451  }
3452  }
3453 
3454  if( desiredcurv != SCIP_EXPRCURV_UNKNOWN )
3455  {
3456  for( f = 0; f < nfactors; ++f )
3457  {
3458  childcurv = SCIPexprgraphGetNodeCurvature(children[childidxs[f]]); /*lint !e613*/
3459  assert(childcurv != SCIP_EXPRCURV_UNKNOWN);
3460  childbounds = SCIPexprgraphGetNodeBounds(children[childidxs[f]]); /*lint !e613*/
3461  assert(childbounds.inf >= 0.0 || childbounds.sup <= 0.0);
3462 
3463  /* negate curvature if factor is negative */
3464  if( childbounds.inf < 0.0 )
3465  childcurv = SCIPexprcurvNegate(childcurv);
3466 
3467  /* check if exp_j * checkcurv is convex (>= 0) and/or concave */
3468  childcurv = SCIPexprcurvMultiply(SCIPexprGetMonomialExponents(monomial)[f], childcurv);
3469  if( (desiredcurv == SCIP_EXPRCURV_CONVEX && !(childcurv & SCIP_EXPRCURV_CONVEX )) ||
3470  (desiredcurv == SCIP_EXPRCURV_CONCAVE && !(childcurv & SCIP_EXPRCURV_CONCAVE)) )
3471  {
3472  SCIPdebugMsg(scip, "reform child %d (factor %d) with curv %s into var\n",
3473  childidxs[f], f, SCIPexprcurvGetName(SCIPexprgraphGetNodeCurvature(children[childidxs[f]]))); /*lint !e613*/
3474  SCIP_CALL( reformNode2Var(scip, exprgraph, children[childidxs[f]], conss, nconss, naddcons, FALSE) ); /*lint !e613*/
3475  modified = TRUE;
3476  }
3477  }
3478  }
3479  }
3480 
3481  if( modified )
3482  {
3483  /* refresh curvature information in node, since we changed children, it should be convex or concave now */
3484  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(node, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
3485  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(node)));
3486  assert(SCIPexprgraphGetNodeCurvature(node) != SCIP_EXPRCURV_UNKNOWN);
3487 
3488  /* we are done and can proceed with the next node */
3489  ++i;
3490  break;
3491  }
3492 
3493  /* monomial can only have unknown curvature here, if it has several factors
3494  * or is of form x^a with x both negative and positive and a an odd or negative integer (-> INTPOWER expression)
3495  */
3496  assert(nfactors > 1 ||
3497  (SCIPexprgraphGetNodeBounds(children[childidxs[0]]).inf < 0.0 && SCIPexprgraphGetNodeBounds(children[childidxs[0]]).sup > 0.0 &&
3498  SCIPisIntegral(scip, exponents[0]) && (exponents[0] < 0.0 || ((int)SCIPround(scip, exponents[0]) % 2 != 0)))
3499  ); /*lint !e613*/
3500 
3501  /* bilinear monomials should not come up here, since simplifier should have turned them into quadratic expression nodes */
3502  assert(!(nfactors == 2 && exponents[0] == 1.0 && exponents[1] == 1.0));
3503 
3504  /* reform monomial if it is a product, or we need it to be on the top of the graph, or if it of the form x^a with a < 0.0 (and thus x having mixed sign, see assert above)
3505  * thus, in the case x^a with a an odd positive integer we assume that cons_signpower will do something */
3506  if( nfactors > 1 || havenonlinparent || exponents[0] < 0.0 )
3507  {
3508  SCIP_EXPRGRAPHNODE* auxnode;
3509  SCIP_EXPRGRAPHNODE** factors;
3510 
3511  if( nfactors > 1 )
3512  {
3513  SCIP_CALL( SCIPallocBufferArray(scip, &factors, nfactors) );
3514  for( f = 0; f < nfactors; ++f )
3515  factors[f] = children[childidxs[f]]; /*lint !e613*/
3516  }
3517  else
3518  factors = &children[childidxs[0]]; /*lint !e613*/
3519 
3520  SCIPdebugMsg(scip, "reform monomial node, create auxvar = %u\n", havenonlinparent);
3521  /* get new auxnode for monomial
3522  * if node has parents and monomial is of indefinite form x^a, then also create auxvar for it, since otherwise we create a auxnode with unknown curvature
3523  * note, that the case x^a with positive and odd a will still give an indefinite node (without parents), where we assume that signpower will pick it up at some point
3524  */
3525  SCIP_CALL( reformMonomial(scip, exprgraph, nfactors, factors, exponents, &auxnode, havenonlinparent, SCIPexprgraphGetNodeDepth(node), naddcons) );
3526 
3527  if( nfactors > 1 )
3528  {
3529  SCIPfreeBufferArray(scip, &factors);
3530  }
3531 
3532  /* create node for monomialcoef * auxnode + monomialconstant, if not identical to auxnode */
3533  if( SCIPexprgraphGetNodePolynomialConstant(node) != 0.0 || coef != 1.0 )
3534  {
3535  SCIP_EXPRGRAPHNODE* replnode;
3536 
3537  SCIP_CALL( SCIPexprgraphCreateNodeLinear(SCIPblkmem(scip), &replnode, 1, &coef, SCIPexprgraphGetNodePolynomialConstant(node)) );
3538  SCIP_CALL( SCIPexprgraphAddNode(exprgraph, replnode, -1, 1, &auxnode) );
3539  auxnode = replnode;
3540  }
3541 
3542  /* replace node by auxnode and refresh its curvature */
3543  SCIP_CALL( reformReplaceNode(exprgraph, &node, auxnode, conss, nconss) );
3544  SCIP_CALL( SCIPexprgraphUpdateNodeBoundsCurvature(auxnode, INTERVALINFTY, BOUNDTIGHTENING_MINSTRENGTH, TRUE) );
3545  assert(!SCIPintervalIsEmpty(INTERVALINFTY, SCIPexprgraphGetNodeBounds(auxnode)));
3546 
3547  break;
3548  }
3549  else
3550  {
3551  SCIPdebugMsg(scip, "no reformulation of monomial node, assume signpower will take care of it\n");
3552  }
3553 
3554  ++i;
3555  break;
3556  }
3557 
3558  case SCIP_EXPR_USER:
3559  {
3560  /* ensure all children are linear */
3561  SCIP_CALL( reformEnsureChildrenMinCurvature( scip, exprgraph, node, SCIP_EXPRCURV_LINEAR, conss, nconss, naddcons ) );
3562 
3563  /* unknown curvature can be handled by user estimator callback or interval gradient */
3564  /*
3565  if( SCIPexprgraphGetNodeCurvature( node ) == SCIP_EXPRCURV_UNKNOWN )
3566  {
3567  SCIPerrorMessage("user expression with unknown curvature not supported\n");
3568  return SCIP_ERROR;
3569  }
3570  */
3571 
3572  ++i;
3573  break;
3574  }
3575 
3576  case SCIP_EXPR_LAST:
3577  SCIPABORT();
3578  break;
3579  }
3580  }
3581  }
3582 
3583  /* for constraints with concave f(g(x)) with linear g:R^n -> R, n>1, reformulate to get a univariate concave function, since this is easier to underestimate
3584  * @todo this does not work yet for sums of functions other than polynomials
3585  */
3586  for( c = 0; c < nconss; ++c )
3587  {
3588  SCIP_EXPRGRAPHNODE* multivarnode;
3589  SCIP_EXPRCURV curv;
3590 
3591  assert(conss[c] != NULL); /*lint !e613*/
3592 
3593  /* skip constraints that are to be deleted */
3594  if( SCIPconsIsDeleted(conss[c]) ) /*lint !e613*/
3595  continue;
3596 
3597  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
3598  assert(consdata != NULL);
3599 
3600  if( consdata->exprgraphnode == NULL )
3601  continue;
3602 
3603  /* after reformulation, force a round of backpropagation in expression graph for all constraints,
3604  * since new variables (nlreform*) may now be used in existing constraints and we want domain restrictions
3605  * of operators propagated for these variables
3606  */
3607  consdata->forcebackprop = TRUE;
3608 
3609  if( SCIPexprgraphGetNodeOperator(consdata->exprgraphnode) == SCIP_EXPR_POLYNOMIAL )
3610  {
3611  SCIP_EXPRDATA_MONOMIAL* monomial;
3612  int m;
3613  int f;
3614 
3615  for( m = 0; m < SCIPexprgraphGetNodePolynomialNMonomials(consdata->exprgraphnode); ++m )
3616  {
3617  SCIP_CALL( SCIPexprgraphGetNodePolynomialMonomialCurvature(consdata->exprgraphnode, m, INTERVALINFTY, &curv) );
3618 
3619  monomial = SCIPexprgraphGetNodePolynomialMonomials(consdata->exprgraphnode)[m];
3620  assert(monomial != NULL);
3621 
3622  /* if nothing concave, then continue */
3623  if( (SCIPisInfinity(scip, consdata->rhs) || curv != SCIP_EXPRCURV_CONCAVE) &&
3624  ( SCIPisInfinity(scip, -consdata->lhs) || curv != SCIP_EXPRCURV_CONVEX) )
3625  continue;
3626 
3627  for( f = 0; f < SCIPexprGetMonomialNFactors(monomial); ++f )
3628  {
3629  multivarnode = SCIPexprgraphGetNodeChildren(consdata->exprgraphnode)[SCIPexprGetMonomialChildIndices(monomial)[f]];
3630 
3631  /* search for a descendant of node that has > 1 children
3632  * after simplifier run, there should be no constant expressions left
3633  */
3634  while( SCIPexprgraphGetNodeNChildren(multivarnode) == 1 )
3635  multivarnode = SCIPexprgraphGetNodeChildren(multivarnode)[0];
3636 
3637  /* if node expression is obviously univariate, then continue */
3638  if( SCIPexprgraphGetNodeNChildren(multivarnode) == 0 )
3639  {
3640  assert(SCIPexprgraphGetNodeOperator(multivarnode) == SCIP_EXPR_CONST || SCIPexprgraphGetNodeOperator(multivarnode) == SCIP_EXPR_VARIDX);
3641  continue;
3642  }
3643 
3644  /* if multivarnode is a linear expression, then replace this by an auxiliary variable/node
3645  * mark auxiliary variable as not to multiaggregate, so SCIP cannot undo what we just did
3646  */
3647  if( SCIPexprgraphGetNodeCurvature(multivarnode) == SCIP_EXPRCURV_LINEAR )
3648  {
3649  SCIPdebugMsg(scip, "replace linear multivariate node %p(%d,%d) in expression of cons <%s> by auxvar\n",
3650  (void*)multivarnode, SCIPexprgraphGetNodeDepth(multivarnode), SCIPexprgraphGetNodePosition(multivarnode), SCIPconsGetName(conss[c])); /*lint !e613*/
3651  SCIPdebugPrintCons(scip, conss[c], NULL); /*lint !e613*/
3652  SCIP_CALL( reformNode2Var(scip, exprgraph, multivarnode, conss, nconss, naddcons, TRUE) );
3653  }
3654  }
3655  }
3656  }
3657  else
3658  {
3659  curv = SCIPexprgraphGetNodeCurvature(consdata->exprgraphnode);
3660 
3661  /* if nothing concave, then continue */
3662  if( (SCIPisInfinity(scip, consdata->rhs) || curv != SCIP_EXPRCURV_CONCAVE) &&
3663  ( SCIPisInfinity(scip, -consdata->lhs) || curv != SCIP_EXPRCURV_CONVEX) )
3664  continue;
3665 
3666  /* search for a descendant of node that has > 1 children
3667  * after simplifier run, there should be no constant expressions left
3668  */
3669  multivarnode = consdata->exprgraphnode;
3670  while( SCIPexprgraphGetNodeNChildren(multivarnode) == 1 )
3671  multivarnode = SCIPexprgraphGetNodeChildren(multivarnode)[0];
3672 
3673  /* if node expression is obviously univariate, then continue */
3674  if( SCIPexprgraphGetNodeNChildren(multivarnode) == 0 )
3675  {
3676  assert(SCIPexprgraphGetNodeOperator(multivarnode) == SCIP_EXPR_CONST || SCIPexprgraphGetNodeOperator(multivarnode) == SCIP_EXPR_VARIDX);
3677  continue;
3678  }
3679 
3680  /* if node itself is multivariate, then continue */
3681  if( multivarnode == consdata->exprgraphnode )
3682  continue;
3683 
3684  /* if multivarnode is a linear expression, then replace this by an auxiliary variable/node
3685  * mark auxiliary variable as not to multiaggregate, so SCIP cannot undo what we just did
3686  */
3687  if( SCIPexprgraphGetNodeCurvature(multivarnode) == SCIP_EXPRCURV_LINEAR )
3688  {
3689  SCIPdebugMsg(scip, "replace linear multivariate node %p(%d,%d) in expression of cons <%s> by auxvar\n",
3690  (void*)multivarnode, SCIPexprgraphGetNodeDepth(multivarnode), SCIPexprgraphGetNodePosition(multivarnode), SCIPconsGetName(conss[c])); /*lint !e613*/
3691  SCIPdebugPrintCons(scip, conss[c], NULL); /*lint !e613*/
3692  SCIP_CALL( reformNode2Var(scip, exprgraph, multivarnode, conss, nconss, naddcons, TRUE) );
3693  }
3694  }
3695  }
3696 
3697  conshdlrdata->isreformulated = TRUE;
3698 
3699  return SCIP_OKAY;
3700 }
3701 
3702 /** computes activity and violation of a constraint
3703  *
3704  * During presolving and if the constraint is active, it is assumes that SCIPexprgraphEval has been called for sol before.
3705  *
3706  * If a solution is found to violate the variable bounds, then violation calculation is stopped and solviolbounds is set to TRUE.
3707  */
3708 static
3709 SCIP_RETCODE computeViolation(
3710  SCIP* scip, /**< SCIP data structure */
3711  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3712  SCIP_CONS* cons, /**< nonlinear constraint */
3713  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
3714  SCIP_Bool* solviolbounds /**< buffer to indicate whether solution is found to violate variable bounds by more than feastol */
3715  )
3716 { /*lint --e{666}*/
3717  SCIP_CONSHDLRDATA* conshdlrdata;
3718  SCIP_CONSDATA* consdata;
3719  SCIP_VAR* var;
3720  SCIP_Real varval;
3721  int i;
3722 
3723  assert(scip != NULL);
3724  assert(conshdlr != NULL);
3725  assert(cons != NULL);
3726  assert(solviolbounds != NULL);
3727 
3728  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3729  assert(conshdlrdata != NULL);
3730  assert(conshdlrdata->exprinterpreter != NULL);
3731 
3732  consdata = SCIPconsGetData(cons);
3733  assert(consdata != NULL);
3734 
3735  consdata->activity = 0.0;
3736  consdata->lhsviol = 0.0;
3737  consdata->rhsviol = 0.0;
3738  varval = 0.0;
3739  *solviolbounds = FALSE;
3740 
3741  for( i = 0; i < consdata->nlinvars; ++i )
3742  {
3743  SCIP_Real activity;
3744 
3745  var = consdata->linvars[i];
3746  varval = SCIPgetSolVal(scip, sol, var);
3747 
3748  /* project onto local box, in case the LP solution is slightly outside the bounds (which is not our job to enforce) */
3749  if( sol == NULL )
3750  {
3751  /* with non-initial columns, this might fail because variables can shortly be a column variable before entering the LP and have value 0.0 in this case */
3752  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) && !SCIPisFeasGE(scip, varval, SCIPvarGetLbLocal(var))) ||
3753  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) && !SCIPisFeasLE(scip, varval, SCIPvarGetUbLocal(var))) )
3754  {
3755  *solviolbounds = TRUE;
3756  return SCIP_OKAY;
3757  }
3758  varval = MAX(SCIPvarGetLbLocal(var), MIN(SCIPvarGetUbLocal(var), varval));
3759  }
3760  activity = consdata->lincoefs[i] * varval;
3761 
3762  /* the contribution of a variable with |varval| = +inf is +inf when activity > 0.0, -inf when activity < 0.0, and
3763  * 0.0 otherwise
3764  */
3765  if( SCIPisInfinity(scip, REALABS(varval)) )
3766  {
3767  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
3768  {
3769  consdata->activity = SCIPinfinity(scip);
3770  consdata->rhsviol = SCIPinfinity(scip);
3771  return SCIP_OKAY;
3772  }
3773 
3774  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
3775  {
3776  consdata->activity = -SCIPinfinity(scip);
3777  consdata->lhsviol = SCIPinfinity(scip);
3778  return SCIP_OKAY;
3779  }
3780  }
3781 
3782  consdata->activity += activity;
3783  }
3784 
3785  for( i = 0; i < consdata->nexprtrees; ++i )
3786  {
3787  SCIP_Real activity;
3788  SCIP_Real val;
3789  int nvars;
3790 
3791  /* compile expression tree, if not done before */
3792  if( SCIPexprtreeGetInterpreterData(consdata->exprtrees[i]) == NULL )
3793  {
3794  SCIP_CALL( SCIPexprintCompile(conshdlrdata->exprinterpreter, consdata->exprtrees[i]) );
3795  }
3796 
3797  nvars = SCIPexprtreeGetNVars(consdata->exprtrees[i]);
3798 
3799  if( nvars == 1 )
3800  {
3801  /* in the not so unusual case that an expression has only one variable, we do not need to extra allocate memory */
3802  var = SCIPexprtreeGetVars(consdata->exprtrees[i])[0];
3803  varval = SCIPgetSolVal(scip, sol, var);
3804 
3805  /* project onto local box, in case the LP solution is slightly outside the bounds (and then cannot be evaluated) */
3806  if( sol == NULL )
3807  {
3808  /* with non-initial columns, this might fail because variables can shortly be a column variable before entering the LP and have value 0.0 in this case */
3809  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) && !SCIPisFeasGE(scip, varval, SCIPvarGetLbLocal(var))) ||
3810  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) && !SCIPisFeasLE(scip, varval, SCIPvarGetUbLocal(var))) )
3811  {
3812  *solviolbounds = TRUE;
3813  return SCIP_OKAY;
3814  }
3815  varval = MAX(SCIPvarGetLbLocal(var), MIN(SCIPvarGetUbLocal(var), varval));
3816  }
3817 
3818  /* coverity[callee_ptr_arith] */
3819  SCIP_CALL( SCIPexprintEval(conshdlrdata->exprinterpreter, consdata->exprtrees[i], &varval, &val) );
3820  }
3821  else
3822  {
3823  SCIP_Real* x;
3824  int j;
3825 
3826  SCIP_CALL( SCIPallocBufferArray(scip, &x, nvars) );
3827 
3828  for( j = 0; j < nvars; ++j )
3829  {
3830  var = SCIPexprtreeGetVars(consdata->exprtrees[i])[j];
3831  varval = SCIPgetSolVal(scip, sol, var);
3832 
3833  /* project onto local box, in case the LP solution is slightly outside the bounds (and then cannot be evaluated) */
3834  if( sol == NULL )
3835  {
3836  /* with non-initial columns, this might fail because variables can shortly be a column variable before entering the LP and have value 0.0 in this case */
3837  if( (!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) && !SCIPisFeasGE(scip, varval, SCIPvarGetLbLocal(var))) ||
3838  (!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) && !SCIPisFeasLE(scip, varval, SCIPvarGetUbLocal(var))) )
3839  {
3840  *solviolbounds = TRUE;
3841  SCIPfreeBufferArray(scip, &x);
3842  return SCIP_OKAY;
3843  }
3844  varval = MAX(SCIPvarGetLbLocal(var), MIN(SCIPvarGetUbLocal(var), varval));
3845  }
3846 
3847  x[j] = varval;
3848  }
3849 
3850  SCIP_CALL( SCIPexprintEval(conshdlrdata->exprinterpreter, consdata->exprtrees[i], x, &val) );
3851 
3852  SCIPfreeBufferArray(scip, &x);
3853  }
3854 
3855  /* set the activity to infinity if a function evaluation was not valid (e.g., sqrt(-1) ) */
3856  if( !SCIPisFinite(val) )
3857  {
3858  consdata->activity = SCIPinfinity(scip);
3859  if( !SCIPisInfinity(scip, -consdata->lhs) )
3860  consdata->lhsviol = SCIPinfinity(scip);
3861  if( !SCIPisInfinity(scip, consdata->rhs) )
3862  consdata->rhsviol = SCIPinfinity(scip);
3863  return SCIP_OKAY;
3864  }
3865 
3866  /* the contribution of an expression with |val| = +inf is +inf when its coefficient is > 0.0, -inf when its coefficient is < 0.0, and
3867  * 0.0 otherwise
3868  */
3869  activity = consdata->nonlincoefs[i] * val;
3870  if( SCIPisInfinity(scip, REALABS(val)) )
3871  {
3872  if( activity > 0.0 && !SCIPisInfinity(scip, consdata->rhs) )
3873  {
3874  consdata->activity = SCIPinfinity(scip);
3875  consdata->rhsviol = SCIPinfinity(scip);
3876  return SCIP_OKAY;
3877  }
3878 
3879  if( activity < 0.0 && !SCIPisInfinity(scip, -consdata->lhs) )
3880  {
3881  consdata->activity = -SCIPinfinity(scip);
3882  consdata->lhsviol = SCIPinfinity(scip);
3883  return SCIP_OKAY;
3884  }
3885  }
3886 
3887  consdata->activity += activity;
3888  }
3889 
3890  if( consdata->nexprtrees == 0 && consdata->exprgraphnode != NULL )
3891  {
3892  SCIP_Real val;
3893 
3894  assert(SCIPgetStage(scip) >= SCIP_STAGE_INITPRESOLVE && SCIPgetStage(scip) <= SCIP_STAGE_EXITPRESOLVE);
3895 
3896  val = SCIPexprgraphGetNodeVal(consdata->exprgraphnode);
3897  assert(val != SCIP_INVALID); /*lint !e777*/
3898 
3899  /* set the activity to infinity if a function evaluation was not valid (e.g., sqrt(-1) ) */
3900  if( !SCIPisFinite(val) )
3901  {
3902  consdata->activity = SCIPinfinity(scip);
3903  if( !SCIPisInfinity(scip, -consdata->lhs) )
3904  consdata->lhsviol = SCIPinfinity(scip);
3905  if( !SCIPisInfinity(scip, consdata->rhs) )
3906  consdata->rhsviol = SCIPinfinity(scip);
3907  return SCIP_OKAY;
3908  }
3909 
3910  if( SCIPisInfinity(scip, val) && !SCIPisInfinity(scip, consdata->rhs) )
3911  {
3912  consdata->activity = SCIPinfinity(scip);
3913  consdata->rhsviol = SCIPinfinity(scip);
3914  return SCIP_OKAY;
3915  }
3916  else if( SCIPisInfinity(scip, -val) && !SCIPisInfinity(scip, -consdata->lhs) )
3917  {
3918  consdata->activity = -SCIPinfinity(scip);
3919  consdata->lhsviol = SCIPinfinity(scip);
3920  return SCIP_OKAY;
3921  }
3922 
3923  consdata->activity += val;
3924  }
3925 
3926  if( !SCIPisInfinity(scip, -consdata->lhs) && SCIPisGT(scip, consdata->lhs - consdata->activity, SCIPfeastol(scip)) )
3927  consdata->lhsviol = consdata->lhs - consdata->activity;
3928  else
3929  consdata->lhsviol = 0.0;
3930 
3931  if( !SCIPisInfinity(scip, consdata->rhs) && SCIPisGT(scip, consdata->activity - consdata->rhs, SCIPfeastol(scip)) )
3932  consdata->rhsviol = consdata->activity - consdata->rhs;
3933  else
3934  consdata->rhsviol = 0.0;
3935 
3936  /* update absolute and relative violation of the solution */
3937  if( sol != NULL )
3938  {
3939  SCIP_Real absviol;
3940  SCIP_Real relviol;
3941  SCIP_Real lhsrelviol;
3942  SCIP_Real rhsrelviol;
3943 
3944  absviol = MAX(consdata->lhsviol, consdata->rhsviol);
3945 
3946  lhsrelviol = SCIPrelDiff(consdata->lhs, consdata->activity);
3947  rhsrelviol = SCIPrelDiff(consdata->activity, consdata->rhs);
3948  relviol = MAX(lhsrelviol, rhsrelviol);
3949 
3950  SCIPupdateSolConsViolation(scip, sol, absviol, relviol);
3951  }
3952 
3953  return SCIP_OKAY;
3954 }
3955 
3956 /** computes violation of a set of constraints
3957  *
3958  * If the solution is found to violate bounds of some variable in some constraint, then violation computation is stopped and solviolbounds is set to TRUE.
3959  */
3960 static
3961 SCIP_RETCODE computeViolations(
3962  SCIP* scip, /**< SCIP data structure */
3963  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3964  SCIP_CONS** conss, /**< constraints */
3965  int nconss, /**< number of constraints */
3966  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
3967  SCIP_Bool* solviolbounds, /**< buffer to indicate whether solution violates bounds of some variable by more than feastol */
3968  SCIP_CONS** maxviolcon /**< buffer to store constraint with largest violation, or NULL if solution is feasible */
3969  )
3970 {
3971  SCIP_CONSDATA* consdata;
3972  SCIP_Real viol;
3973  SCIP_Real maxviol;
3974  int c;
3975 
3976  assert(scip != NULL);
3977  assert(conshdlr != NULL);
3978  assert(conss != NULL || nconss == 0);
3979  assert(maxviolcon != NULL);
3980  assert(solviolbounds != NULL);
3981 
3982  if( SCIPgetStage(scip) >= SCIP_STAGE_INITPRESOLVE && SCIPgetStage(scip) <= SCIP_STAGE_EXITPRESOLVE )
3983  {
3984  SCIP_CONSHDLRDATA* conshdlrdata;
3985  SCIP_Real* varvals;
3986 
3987  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3988  assert(conshdlrdata != NULL);
3989  assert(conshdlrdata->exprgraph != NULL);
3990 
3991  SCIP_CALL( SCIPallocBufferArray(scip, &varvals, SCIPexprgraphGetNVars(conshdlrdata->exprgraph)) );
3992  SCIP_CALL( SCIPgetSolVals(scip, sol, SCIPexprgraphGetNVars(conshdlrdata->exprgraph), (SCIP_VAR**)SCIPexprgraphGetVars(conshdlrdata->exprgraph), varvals) );
3993 
3994  SCIP_CALL( SCIPexprgraphEval(conshdlrdata->exprgraph, varvals) );
3995 
3996  SCIPfreeBufferArray(scip, &varvals);
3997  }
3998 
3999  *maxviolcon = NULL;
4000 
4001  maxviol = 0.0;
4002 
4003  for( c = 0; c < nconss; ++c )
4004  {
4005  assert(conss != NULL);
4006  assert(conss[c] != NULL);
4007 
4008  SCIP_CALL( computeViolation(scip, conshdlr, conss[c], sol, solviolbounds) );
4009 
4010  /* stop if solution violates bounds */
4011  if( *solviolbounds )
4012  break;
4013 
4014  consdata = SCIPconsGetData(conss[c]);
4015  assert(consdata != NULL);
4016 
4017  viol = MAX(consdata->lhsviol, consdata->rhsviol);
4018  if( viol > maxviol && SCIPisGT(scip, viol, SCIPfeastol(scip)) )
4019  {
4020  maxviol = viol;
4021  *maxviolcon = conss[c];
4022  }
4023 
4024  /* SCIPdebugMsg(scip, "constraint <%s> violated by (%g, %g), activity = %g\n", SCIPconsGetName(conss[c]), consdata->lhsviol, consdata->rhsviol, consdata->activity); */
4025  }
4026 
4027  return SCIP_OKAY;
4028 }
4029 
4030 /** adds linearization of a constraints expression tree in reference point to a row */
4031 static
4032 SCIP_RETCODE addLinearization(
4033  SCIP* scip, /**< SCIP data structure */
4034  SCIP_EXPRINT* exprint, /**< expression interpreter */
4035  SCIP_CONS* cons, /**< constraint */
4036  int exprtreeidx, /**< for which tree a linearization should be added */
4037  SCIP_Real* x, /**< value of expression tree variables where to generate cut */
4038  SCIP_Bool newx, /**< whether the last evaluation of the expression with the expression interpreter was not at x */
4039  SCIP_ROWPREP* rowprep, /**< rowprep where to add linearization */
4040  SCIP_Bool* success /**< buffer to store whether a linearization was succefully added to the row */
4041  )
4042 {
4043  SCIP_CONSDATA* consdata;
4044  SCIP_EXPRTREE* exprtree;
4045  SCIP_Real treecoef;
4046  SCIP_Real val;
4047  SCIP_Real* grad;
4048  SCIP_Real constant = 0.0;
4049  SCIP_Bool perturbedx;
4050  int nvars;
4051  int i;
4052 
4053  assert(scip != NULL);
4054  assert(cons != NULL);
4055  assert(x != NULL);
4056  assert(rowprep != NULL);
4057  assert(success != NULL);
4058 
4059  consdata = SCIPconsGetData(cons);
4060  assert(consdata != NULL);
4061  assert(exprtreeidx >= 0);
4062  assert(exprtreeidx < consdata->nexprtrees);
4063  assert(consdata->exprtrees != NULL);
4064 
4065  exprtree = consdata->exprtrees[exprtreeidx];
4066  assert(exprtree != NULL);
4067  assert(newx || SCIPexprtreeGetInterpreterData(exprtree) != NULL);
4068 
4069  treecoef = consdata->nonlincoefs[exprtreeidx];
4070 
4071  *success = FALSE;
4072 
4073  /* compile expression if evaluated the first time; can only happen if newx is FALSE */
4074  if( newx && SCIPexprtreeGetInterpreterData(exprtree) == NULL )
4075  {
4076  SCIP_CALL( SCIPexprintCompile(exprint, exprtree) );
4077  }
4078 
4079  nvars = SCIPexprtreeGetNVars(exprtree);
4080  SCIP_CALL( SCIPallocBufferArray(scip, &grad, nvars) );
4081 
4082  perturbedx = FALSE;
4083  do
4084  {
4085  /* get value and gradient */
4086  SCIP_CALL( SCIPexprintGrad(exprint, exprtree, x, newx, &val, grad) );
4087  if( SCIPisFinite(val) && !SCIPisInfinity(scip, REALABS(val)) )
4088  {
4089  val *= treecoef;
4090  /* check gradient entries and compute constant f(refx) - grad * refx */
4091  constant = val;
4092  for( i = 0; i < nvars; ++i )
4093  {
4094  if( !SCIPisFinite(grad[i]) || SCIPisInfinity(scip, grad[i]) || SCIPisInfinity(scip, -grad[i]) )
4095  break;
4096 
4097  grad[i] *= treecoef;
4098  constant -= grad[i] * x[i];
4099 
4100  /* try to perturb x if the constant is too large */
4101  if( SCIPisInfinity(scip, REALABS(constant)) )
4102  break;
4103 
4104  /* coefficients smaller than epsilon are rounded to 0.0 when added to row, this can be wrong if variable value is very large (bad numerics)
4105  * in this case, set gradient to 0.0 here, but modify constant so that cut is still valid (if possible)
4106  * i.e., estimate grad[i]*x >= grad[i] * bound(x) or grad[i]*x <= grad[i] * bound(x), depending on whether we compute an underestimator (convex) or an overestimator (concave)
4107  * if required bound of x is not finite, then give up
4108  */
4109  if( grad[i] != 0.0 && SCIPisZero(scip, grad[i]) )
4110  {
4111  SCIP_VAR* var;
4112  SCIP_Real xbnd;
4113 
4114  var = SCIPexprtreeGetVars(exprtree)[i];
4115  if( consdata->curvatures[exprtreeidx] & SCIP_EXPRCURV_CONVEX )
4116  {
4117  xbnd = grad[i] > 0.0 ? SCIPvarGetLbGlobal(var) : SCIPvarGetUbGlobal(var);
4118  }
4119  else
4120  {
4121  assert(consdata->curvatures[exprtreeidx] & SCIP_EXPRCURV_CONCAVE);
4122  xbnd = grad[i] > 0.0 ? SCIPvarGetUbGlobal(var) : SCIPvarGetLbGlobal(var);
4123  }
4124  if( !SCIPisInfinity(scip, REALABS(xbnd)) )
4125  {
4126  SCIPdebugMsg(scip, "var <%s> [%g,%g] has tiny gradient %g, replace coefficient by constant %g\n",
4127  SCIPvarGetName(var), SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var), grad[i], grad[i] * xbnd);
4128  constant += grad[i] * xbnd;
4129  grad[i] = 0.0;
4130  }
4131  else
4132  {
4133  *success = FALSE;
4134  SCIPdebugMsg(scip, "skipping linearization, var <%s> [%g,%g] has tiny gradient %g but no finite bound in this direction\n",
4135  SCIPvarGetName(var), SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var), grad[i]);
4136  SCIPfreeBufferArray(scip, &grad);
4137  return SCIP_OKAY;
4138  }
4139  }
4140  }
4141 
4142  if( i == nvars )
4143  break;
4144  }
4145 
4146  SCIPdebugMsg(scip, "got nonfinite value in evaluation or gradient of <%s>: ", SCIPconsGetName(cons));
4147  if( !perturbedx )
4148  {
4149  SCIP_Real lb;
4150  SCIP_Real ub;
4151 
4152  SCIPdebugMsgPrint(scip, "perturbing reference point and trying again\n");
4153  for( i = 0; i < nvars; ++i )
4154  {
4155  lb = SCIPvarGetLbGlobal(SCIPexprtreeGetVars(exprtree)[i]);
4156  ub = SCIPvarGetUbGlobal(SCIPexprtreeGetVars(exprtree)[i]);
4157  if( SCIPisEQ(scip, x[i], lb) )
4158  x[i] += MIN(0.9*(ub-lb), i*SCIPfeastol(scip)); /*lint !e666*/
4159  else if( SCIPisEQ(scip, x[i], ub) )
4160  x[i] -= MIN(0.9*(ub-lb), i*SCIPfeastol(scip)); /*lint !e666*/
4161  else
4162  x[i] += MIN3(0.9*(ub-x[i]), 0.9*(x[i]-lb), i*SCIPfeastol(scip)) * (i%2 != 0 ? -1.0 : 1.0); /*lint !e666*/
4163  }
4164  newx = TRUE;
4165  perturbedx = TRUE;
4166  }
4167  else
4168  {
4169  SCIPdebugMsgPrint(scip, "skipping linearization\n");
4170  SCIPfreeBufferArray(scip, &grad);
4171  return SCIP_OKAY;
4172  }
4173  }
4174  while( TRUE ); /*lint !e506*/
4175 
4176  /* add linearization to SCIP row */
4177  SCIPaddRowprepConstant(rowprep, constant);
4178  SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, nvars, SCIPexprtreeGetVars(exprtree), grad) );
4179 
4180  *success = TRUE;
4181 
4182  SCIPfreeBufferArray(scip, &grad);
4183 
4184  SCIPdebugMsg(scip, "added linearization for tree %d of constraint <%s>\n", exprtreeidx, SCIPconsGetName(cons));
4185  SCIPdebug( SCIPprintRowprep(scip, rowprep, NULL) );
4186 
4187  return SCIP_OKAY;
4188 }
4189 
4190 /** adds secant of a constraints univariate expression tree in reference point to a row */
4191 static
4192 SCIP_RETCODE addConcaveEstimatorUnivariate(
4193  SCIP* scip, /**< SCIP data structure */
4194  SCIP_CONS* cons, /**< constraint */
4195  int exprtreeidx, /**< for which tree a secant should be added */
4196  SCIP_ROWPREP* rowprep, /**< rowprep where to add estimator */
4197  SCIP_Bool* success /**< buffer to store whether a secant was succefully added to the row */
4198  )
4199 {
4200  SCIP_CONSDATA* consdata;
4201  SCIP_EXPRTREE* exprtree;
4202  SCIP_Real treecoef;
4203  SCIP_VAR* var;
4204  SCIP_Real xlb;
4205  SCIP_Real xub;
4206  SCIP_Real vallb;
4207  SCIP_Real valub;
4208  SCIP_Real slope;
4209  SCIP_Real constant;
4210 
4211  assert(scip != NULL);
4212  assert(cons != NULL);
4213  assert(rowprep != NULL);
4214  assert(success != NULL);
4215 
4216  consdata = SCIPconsGetData(cons);
4217  assert(consdata != NULL);
4218  assert(exprtreeidx >= 0);
4219  assert(exprtreeidx < consdata->nexprtrees);
4220  assert(consdata->exprtrees != NULL);
4221 
4222  exprtree = consdata->exprtrees[exprtreeidx];
4223  assert(exprtree != NULL);
4224  assert(SCIPexprtreeGetNVars(exprtree) == 1);
4225 
4226  treecoef = consdata->nonlincoefs[exprtreeidx];
4227 
4228  *success = FALSE;
4229 
4230  var = SCIPexprtreeGetVars(exprtree)[0];
4231  xlb = SCIPvarGetLbLocal(var);
4232  xub = SCIPvarGetUbLocal(var);
4233 
4234  /* if variable is unbounded, then cannot really compute secant */
4235  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) )
4236  {
4237  SCIPdebugMsg(scip, "skip secant for tree %d of constraint <%s> since variable is unbounded\n", exprtreeidx, SCIPconsGetName(cons));
4238  return SCIP_OKAY;
4239  }
4240  assert(SCIPisLE(scip, xlb, xub));
4241 
4242  /* coverity[callee_ptr_arith] */
4243  SCIP_CALL( SCIPexprtreeEval(exprtree, &xlb, &vallb) );
4244  if( !SCIPisFinite(vallb) || SCIPisInfinity(scip, REALABS(vallb)) )
4245  {
4246  SCIPdebugMsg(scip, "skip secant for tree %d of constraint <%s> since function cannot be evaluated in lower bound\n", exprtreeidx, SCIPconsGetName(cons));
4247  return SCIP_OKAY;
4248  }
4249  vallb *= treecoef;
4250 
4251  /* coverity[callee_ptr_arith] */
4252  SCIP_CALL( SCIPexprtreeEval(exprtree, &xub, &valub) );
4253  if( !SCIPisFinite(valub) || SCIPisInfinity(scip, REALABS(valub)) )
4254  {
4255  SCIPdebugMsg(scip, "skip secant for tree %d of constraint <%s> since function cannot be evaluated in upper bound\n", exprtreeidx, SCIPconsGetName(cons));
4256  return SCIP_OKAY;
4257  }
4258  valub *= treecoef;
4259 
4260  if( SCIPisEQ(scip, xlb, xub) )
4261  {
4262  slope = 0.0;
4263  /* choose most conservative value for the cut */
4264  if( rowprep->sidetype == SCIP_SIDETYPE_LEFT )
4265  constant = MAX(vallb, valub);
4266  else
4267  constant = MIN(vallb, valub);
4268  }
4269  else
4270  {
4271  slope = (valub - vallb) / (xub - xlb);
4272  constant = vallb - slope * xlb;
4273  }
4274 
4275  /* add secant to SCIP row */
4276  SCIPaddRowprepConstant(rowprep, constant);
4277  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, slope) );
4278 
4279  *success = TRUE;
4280 
4281  SCIPdebugMsg(scip, "added secant for tree %d of constraint <%s>, slope = %g\n", exprtreeidx, SCIPconsGetName(cons), slope);
4282  SCIPdebug( SCIPprintRowprep(scip, rowprep, NULL) );
4283 
4284  return SCIP_OKAY;
4285 }
4286 
4287 /** adds estimator of a constraints bivariate expression tree to a row
4288  * a reference point is given to decide which hyperplane to choose
4289  */
4290 static
4291 SCIP_RETCODE addConcaveEstimatorBivariate(
4292  SCIP* scip, /**< SCIP data structure */
4293  SCIP_CONS* cons, /**< constraint */
4294  int exprtreeidx, /**< for which tree a secant should be added */
4295  SCIP_Real* ref, /**< reference values of expression tree variables where to generate cut */
4296  SCIP_ROWPREP* rowprep, /**< rowprep where to add estimator */
4297  SCIP_Bool* success /**< buffer to store whether a secant was succefully added to the row */
4298  )
4299 {
4300  SCIP_CONSDATA* consdata;
4301  SCIP_EXPRTREE* exprtree;
4302  SCIP_Real treecoef;
4303  SCIP_VAR* x;
4304  SCIP_VAR* y;
4305  SCIP_Real xlb;
4306  SCIP_Real xub;
4307  SCIP_Real ylb;
4308  SCIP_Real yub;
4309 
4310  SCIP_Real coefx;
4311  SCIP_Real coefy;
4312  SCIP_Real constant;
4313 
4314  SCIP_Real p1[2];
4315  SCIP_Real p2[2];
4316  SCIP_Real p3[2];
4317  SCIP_Real p4[2];
4318  SCIP_Real p1val, p2val, p3val, p4val;
4319 
4320  assert(scip != NULL);
4321  assert(cons != NULL);
4322  assert(ref != NULL);
4323  assert(rowprep != NULL);
4324  assert(success != NULL);
4325 
4326  consdata = SCIPconsGetData(cons);
4327  assert(consdata != NULL);
4328  assert(exprtreeidx >= 0);
4329  assert(exprtreeidx < consdata->nexprtrees);
4330  assert(consdata->exprtrees != NULL);
4331 
4332  exprtree = consdata->exprtrees[exprtreeidx];
4333  assert(exprtree != NULL);
4334  assert(SCIPexprtreeGetNVars(exprtree) == 2);
4335 
4336  treecoef = consdata->nonlincoefs[exprtreeidx];
4337 
4338  *success = FALSE;
4339 
4340  x = SCIPexprtreeGetVars(exprtree)[0];
4341  y = SCIPexprtreeGetVars(exprtree)[1];
4342  xlb = SCIPvarGetLbLocal(x);
4343  xub = SCIPvarGetUbLocal(x);
4344  ylb = SCIPvarGetLbLocal(y);
4345  yub = SCIPvarGetUbLocal(y);
4346 
4347  if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) || SCIPisInfinity(scip, -ylb) || SCIPisInfinity(scip, yub) )
4348  {
4349  SCIPdebugMsg(scip, "skip bivariate secant since <%s> or <%s> is unbounded\n", SCIPvarGetName(x), SCIPvarGetName(y));
4350  return SCIP_OKAY;
4351  }
4352 
4353  /* reference point should not be outside of bounds */
4354  assert(SCIPisFeasLE(scip, xlb, ref[0]));
4355  assert(SCIPisFeasGE(scip, xub, ref[0]));
4356  ref[0] = MIN(xub, MAX(xlb, ref[0]));
4357  assert(SCIPisFeasLE(scip, ylb, ref[1]));
4358  assert(SCIPisFeasGE(scip, yub, ref[1]));
4359  ref[1] = MIN(yub, MAX(ylb, ref[1]));
4360 
4361  /* lower left */
4362  p1[0] = xlb;
4363  p1[1] = ylb;
4364 
4365  /* lower right */
4366  p2[0] = xub;
4367  p2[1] = ylb;
4368 
4369  /* upper right */
4370  p3[0] = xub;
4371  p3[1] = yub;
4372 
4373  /* upper left */
4374  p4[0] = xlb;
4375  p4[1] = yub;
4376 
4377  if( SCIPisEQ(scip, xlb, xub) && SCIPisEQ(scip, ylb, yub) )
4378  {
4379  SCIP_CALL( SCIPexprtreeEval(exprtree, p1, &p1val) );
4380 
4381  if( !SCIPisFinite(p1val) || SCIPisInfinity(scip, REALABS(p1val)) )
4382  {
4383  SCIPdebugMsg(scip, "skip secant for tree %d of constraint <%s> since function cannot be evaluated\n", exprtreeidx, SCIPconsGetName(cons));
4384  return SCIP_OKAY;
4385  }
4386 
4387  p1val *= treecoef;
4388 
4389  coefx = 0.0;
4390  coefy = 0.0;
4391  constant = p1val;
4392  }
4393  else if( SCIPisEQ(scip, xlb, xub) )
4394  {
4395  /* secant between p1 and p4: p1val + [(p4val - p1val) / (yub - ylb)] * (y - ylb) */
4396  assert(!SCIPisEQ(scip, ylb, yub));
4397 
4398  SCIP_CALL( SCIPexprtreeEval(exprtree, p1, &p1val) );
4399  SCIP_CALL( SCIPexprtreeEval(exprtree, p4, &p4val) );
4400  if( !SCIPisFinite(p1val) || SCIPisInfinity(scip, REALABS(p1val)) || !SCIPisFinite(p4val) || SCIPisInfinity(scip, REALABS(p4val)) )
4401  {
4402  SCIPdebugMsg(scip, "skip secant for tree %d of constraint <%s> since function cannot be evaluated\n", exprtreeidx, SCIPconsGetName(cons));
4403  return SCIP_OKAY;
4404  }
4405  p1val *= treecoef;
4406  p4val *= treecoef;
4407 
4408  coefx = 0.0;
4409  coefy = (p4val - p1val) / (yub - ylb);
4410  constant = p1val - coefy * ylb;
4411  }
4412  else if( SCIPisEQ(scip, ylb, yub) )
4413  {
4414  /* secant between p1 and p2: p1val + [(p2val - p1val) / (xub - xlb)] * (x - xlb) */
4415  assert(!SCIPisEQ(scip, xlb, xub));
4416 
4417  SCIP_CALL( SCIPexprtreeEval(exprtree, p1, &p1val) );
4418  SCIP_CALL( SCIPexprtreeEval(exprtree, p2, &p2val) );
4419  if( !SCIPisFinite(p1val) || SCIPisInfinity(scip, REALABS(p1val)) || !SCIPisFinite(p2val) || SCIPisInfinity(scip, REALABS(p2val)) )
4420  {
4421  SCIPdebugMsg(scip, "skip secant for tree %d of constraint <%s> since function cannot be evaluated\n", exprtreeidx, SCIPconsGetName(cons));
4422  return SCIP_OKAY;
4423  }
4424 
4425  p1val *= treecoef;
4426  p2val *= treecoef;
4427 
4428  coefx = (p2val - p1val) / (xub - xlb);
4429  coefy = 0.0;
4430  constant = p1val - coefx * xlb;
4431  }
4432  else
4433  {
4434  SCIP_Real alpha, beta, gamma_, delta;
4435  SCIP_Bool tryother;
4436  SCIP_Bool doover;
4437 
4438  /* if function is convex, then we want an overestimator, otherwise we want an underestimator */
4439  assert(consdata->curvatures[exprtreeidx] == SCIP_EXPRCURV_CONVEX || consdata->curvatures[exprtreeidx] == SCIP_EXPRCURV_CONCAVE);
4440  doover = (consdata->curvatures[exprtreeidx] & SCIP_EXPRCURV_CONVEX); /*lint !e641*/
4441 
4442  SCIP_CALL( SCIPexprtreeEval(exprtree, p1, &p1val) );
4443  SCIP_CALL( SCIPexprtreeEval(exprtree, p2, &p2val) );
4444  SCIP_CALL( SCIPexprtreeEval(exprtree, p3, &p3val) );
4445  SCIP_CALL( SCIPexprtreeEval(exprtree, p4, &p4val) );
4446  if( !SCIPisFinite(p1val) || SCIPisInfinity(scip, REALABS(p1val)) || !SCIPisFinite(p2val) || SCIPisInfinity(scip, REALABS(p2val)) ||
4447  ! SCIPisFinite(p3val) || SCIPisInfinity(scip, REALABS(p3val)) || !SCIPisFinite(p4val) || SCIPisInfinity(scip, REALABS(p4val)) )
4448  {
4449  SCIPdebugMsg(scip, "skip secant for tree %d of constraint <%s> since function cannot be evaluated\n", exprtreeidx, SCIPconsGetName(cons));
4450  return SCIP_OKAY;
4451  }
4452  p1val *= treecoef;
4453  p2val *= treecoef;
4454  p3val *= treecoef;
4455  p4val *= treecoef;
4456 
4457  /* if we want an underestimator, flip f(x,y), i.e., do as if we compute an overestimator for -f(x,y) */
4458  if( !doover )
4459  {
4460  p1val = -p1val;
4461  p2val = -p2val;
4462  p3val = -p3val;
4463  p4val = -p4val;
4464  }
4465 
4466  SCIPdebugMsg(scip, "p1 = (%g, %g), f(p1) = %g\n", p1[0], p1[1], p1val);
4467  SCIPdebugMsg(scip, "p2 = (%g, %g), f(p2) = %g\n", p2[0], p2[1], p2val);
4468  SCIPdebugMsg(scip, "p3 = (%g, %g), f(p3) = %g\n", p3[0], p3[1], p3val);
4469  SCIPdebugMsg(scip, "p4 = (%g, %g), f(p4) = %g\n", p4[0], p4[1], p4val);
4470 
4471  /* Compute coefficients alpha, beta, gamma (>0), delta such that
4472  * alpha*x + beta*y + gamma*z = delta
4473  * is satisfied by at least three of the corner points (p1,f(p1)), ..., (p4,f(p4)) and
4474  * the fourth corner point lies below this hyperplane.
4475  * Since we assume that f is convex, we then know that all points (x,y,f(x,y)) are below this hyperplane, i.e.,
4476  * alpha*x + beta*y - delta <= -gamma * f(x,y),
4477  * or, equivalently,
4478  * -alpha/gamma*x - beta/gamma*y + delta/gamma >= f(x,y).
4479  */
4480 
4481  tryother = FALSE;
4482  if( ref[1] <= ylb + (yub - ylb)/(xub - xlb) * (ref[0] - xlb) )
4483  {
4484  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p2[0], p2[1], p2val, p3[0], p3[1], p3val,
4485  &alpha, &beta, &gamma_, &delta) );
4486 
4487  assert(SCIPisRelEQ(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4488  assert(SCIPisRelEQ(scip, alpha * p2[0] + beta * p2[1] - delta, -gamma_ * p2val));
4489  assert(SCIPisRelEQ(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4490 
4491  /* if hyperplane through p1,p2,p3 does not overestimate f(p4), then it must be the other variant */
4492  if( alpha * p4[0] + beta * p4[1] + gamma_ * p4val > delta )
4493  tryother = TRUE;
4494  else if( (!SCIPisZero(scip, alpha) && SCIPisZero(scip, alpha/gamma_)) ||
4495  ( !SCIPisZero(scip, beta) && SCIPisZero(scip, beta /gamma_)) )
4496  {
4497  /* if numerically bad, take alternative hyperplane */
4498  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p3[0], p3[1], p3val, p4[0], p4[1],
4499  p4val, &alpha, &beta, &gamma_, &delta) );
4500 
4501  assert(SCIPisRelEQ(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4502  assert(SCIPisRelEQ(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4503  assert(SCIPisRelEQ(scip, alpha * p4[0] + beta * p4[1] - delta, -gamma_ * p4val));
4504 
4505  /* if hyperplane through p1,p3,p4 does not overestimate f(p2), then it must be the other variant */
4506  if( alpha * p2[0] + beta * p2[1] + gamma_ * p2val > delta )
4507  tryother = TRUE;
4508  }
4509  }
4510  else
4511  {
4512  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p3[0], p3[1], p3val, p4[0], p4[1], p4val,
4513  &alpha, &beta, &gamma_, &delta) );
4514 
4515  assert(SCIPisRelEQ(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4516  assert(SCIPisRelEQ(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4517  assert(SCIPisRelEQ(scip, alpha * p4[0] + beta * p4[1] - delta, -gamma_ * p4val));
4518 
4519  /* if hyperplane through p1,p3,p4 does not overestimate f(p2), then it must be the other variant */
4520  if( alpha * p2[0] + beta * p2[1] + gamma_ * p2val > delta )
4521  tryother = TRUE;
4522  else if( (!SCIPisZero(scip, alpha) && SCIPisZero(scip, alpha/gamma_)) ||
4523  ( !SCIPisZero(scip, beta) && SCIPisZero(scip, beta /gamma_)) )
4524  {
4525  /* if numerically bad, take alternative */
4526  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p2[0], p2[1], p2val, p3[0], p3[1],
4527  p3val, &alpha, &beta, &gamma_, &delta) );
4528 
4529  assert(SCIPisRelEQ(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4530  assert(SCIPisRelEQ(scip, alpha * p2[0] + beta * p2[1] - delta, -gamma_ * p2val));
4531  assert(SCIPisRelEQ(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4532 
4533  /* if hyperplane through p1,p2,p3 does not overestimate f(p4), then it must be the other variant */
4534  if( alpha * p4[0] + beta * p4[1] + gamma_ * p4val > delta )
4535  tryother = TRUE;
4536  }
4537  }
4538 
4539  if( tryother )
4540  {
4541  if( ref[1] <= yub + (ylb - yub)/(xub - xlb) * (ref[0] - xlb) )
4542  {
4543  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p2[0], p2[1], p2val, p4[0], p4[1],
4544  p4val, &alpha, &beta, &gamma_, &delta) );
4545 
4546  /* hyperplane should be above (p3,f(p3)) and other points should lie on hyperplane */
4547  assert(SCIPisRelEQ(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4548  assert(SCIPisRelEQ(scip, alpha * p2[0] + beta * p2[1] - delta, -gamma_ * p2val));
4549  assert(SCIPisRelLE(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4550  assert(SCIPisRelEQ(scip, alpha * p4[0] + beta * p4[1] - delta, -gamma_ * p4val));
4551 
4552  if( (!SCIPisZero(scip, alpha) && SCIPisZero(scip, alpha/gamma_)) ||
4553  ( !SCIPisZero(scip, beta) && SCIPisZero(scip, beta /gamma_)) )
4554  {
4555  /* if numerically bad, take alternative */
4556  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p2[0], p2[1], p2val, p3[0], p3[1], p3val, p4[0], p4[1],
4557  p4val, &alpha, &beta, &gamma_, &delta) );
4558 
4559  /* hyperplane should be above (p1,f(p1)) and other points should lie on hyperplane */
4560  assert(SCIPisRelLE(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4561  assert(SCIPisRelEQ(scip, alpha * p2[0] + beta * p2[1] - delta, -gamma_ * p2val));
4562  assert(SCIPisRelEQ(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4563  assert(SCIPisRelEQ(scip, alpha * p4[0] + beta * p4[1] - delta, -gamma_ * p4val));
4564  }
4565  }
4566  else
4567  {
4568  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p2[0], p2[1], p2val, p3[0], p3[1], p3val, p4[0], p4[1],
4569  p4val, &alpha, &beta, &gamma_, &delta) );
4570 
4571  /* hyperplane should be above (p1,f(p1)) and other points should lie on hyperplane */
4572  assert(SCIPisRelLE(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4573  assert(SCIPisRelEQ(scip, alpha * p2[0] + beta * p2[1] - delta, -gamma_ * p2val));
4574  assert(SCIPisRelEQ(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4575  assert(SCIPisRelEQ(scip, alpha * p4[0] + beta * p4[1] - delta, -gamma_ * p4val));
4576 
4577  if( (!SCIPisZero(scip, alpha) && SCIPisZero(scip, alpha/gamma_)) ||
4578  ( !SCIPisZero(scip, beta) && SCIPisZero(scip, beta /gamma_)) )
4579  {
4580  /* if numerically bad, take alternative */
4581  SCIP_CALL( SCIPcomputeHyperplaneThreePoints(scip, p1[0], p1[1], p1val, p2[0], p2[1], p2val, p4[0], p4[1],
4582  p4val, &alpha, &beta, &gamma_, &delta) );
4583 
4584  /* hyperplane should be above (p3,f(p3)) and other points should lie on hyperplane */
4585  assert(SCIPisRelEQ(scip, alpha * p1[0] + beta * p1[1] - delta, -gamma_ * p1val));
4586  assert(SCIPisRelEQ(scip, alpha * p2[0] + beta * p2[1] - delta, -gamma_ * p2val));
4587  assert(SCIPisRelLE(scip, alpha * p3[0] + beta * p3[1] - delta, -gamma_ * p3val));
4588  assert(SCIPisRelEQ(scip, alpha * p4[0] + beta * p4[1] - delta, -gamma_ * p4val));
4589  }
4590  }
4591  }
4592 
4593  SCIPdebugMsg(scip, "alpha = %g, beta = %g, gamma = %g, delta = %g\n", alpha, beta, gamma_, delta);
4594 
4595  /* check if bad luck: should not happen if xlb != xub and ylb != yub and numerics are fine */
4596  if( SCIPisZero(scip, gamma_) )
4597  return SCIP_OKAY;
4598  assert(!SCIPisNegative(scip, gamma_));
4599 
4600  /* flip hyperplane */
4601  if( !doover )
4602  gamma_ = -gamma_;
4603 
4604  coefx = -alpha / gamma_;
4605  coefy = -beta / gamma_;
4606  constant = delta / gamma_;
4607 
4608  /* if we loose coefficients because division by gamma makes them < SCIPepsilon(scip), then better not generate a cut here */
4609  if( (!SCIPisZero(scip, alpha) && SCIPisZero(scip, coefx)) ||
4610  ( !SCIPisZero(scip, beta) && SCIPisZero(scip, coefy)) )
4611  {
4612  SCIPdebugMsg(scip, "skip bivar secant for <%s> tree %d due to bad numerics\n", SCIPconsGetName(cons), exprtreeidx);
4613  return SCIP_OKAY;
4614  }
4615  }
4616 
4617  /* add hyperplane coefs to SCIP row */
4618  SCIPaddRowprepConstant(rowprep, constant);
4619  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, x, coefx) );
4620  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, y, coefy) );
4621 
4622  *success = TRUE;
4623 
4624  SCIPdebugMsg(scip, "added bivariate secant for tree %d of constraint <%s>\n", exprtreeidx, SCIPconsGetName(cons));
4625  SCIPdebug( SCIPprintRowprep(scip, rowprep, NULL) );
4626 
4627  return SCIP_OKAY;
4628 }
4629 
4630 /** internal method using an auxiliary LPI, see addConcaveEstimatorMultivariate() */
4631 static
4632 SCIP_RETCODE _addConcaveEstimatorMultivariate(
4633  SCIP* scip, /**< SCIP data structure */
4634  SCIP_LPI* lpi, /**< auxiliary LPI */
4635  SCIP_CONS* cons, /**< constraint */
4636  int exprtreeidx, /**< for which tree a secant should be added */
4637  SCIP_Real* ref, /**< reference values of expression tree variables where to generate cut */
4638  SCIP_ROWPREP* rowprep, /**< rowprep where to add estimator */
4639  SCIP_VAR** vars, /**< variables of the constraint */
4640  SCIP_EXPRTREE* exprtree, /**< expression tree of constraint */
4641  int nvars, /**< number of variables */
4642  SCIP_Bool doupper, /**< should an upper estimator be computed */
4643  SCIP_Bool* success /**< buffer to store whether a secant was succefully added to the row */
4644  )
4645 {
4646  SCIP_CONSDATA* consdata;
4647  SCIP_Real* val;
4648  SCIP_Real* obj;
4649  SCIP_Real* lb;
4650  SCIP_Real* ub;
4651  SCIP_Real* corner;
4652  SCIP_Real* lhs;
4653  SCIP_Real* rhs;
4654  int* beg;
4655  int* ind;
4656  SCIP_Real lpobj;
4657  int ncols;
4658  int nrows;
4659  int nnonz;
4660  SCIP_Real funcval;
4661  SCIP_Real treecoef;
4662 
4663  int i;
4664  int j;
4665  int idx;
4666 
4667  SCIP_RETCODE lpret;
4668 
4669  assert(lpi != NULL);
4670  assert(nvars <= 10);
4671 
4672  consdata = SCIPconsGetData(cons);
4673  treecoef = consdata->nonlincoefs[exprtreeidx];
4674 
4675  /* columns are cut coefficients plus constant */
4676  ncols = nvars + 1;
4677  SCIP_CALL( SCIPallocBufferArray(scip, &obj, ncols) );
4678  SCIP_CALL( SCIPallocBufferArray(scip, &lb, ncols) );
4679  SCIP_CALL( SCIPallocBufferArray(scip, &ub, ncols) );
4680  corner = lb; /* will not use lb and corner simultaneously, so can share memory */
4681 
4682  /* one row for each corner of domain, i.e., 2^nvars many */
4683  nrows = (int)(1u << nvars);
4684  SCIP_CALL( SCIPallocBufferArray(scip, &lhs, nrows) );
4685  SCIP_CALL( SCIPallocBufferArray(scip, &rhs, nrows) );
4686 
4687  /* the coefficients matrix will have at most ncols * nrows many nonzeros */
4688  nnonz = nrows * ncols;
4689  SCIP_CALL( SCIPallocBufferArray(scip, &beg, nrows+1) );
4690  SCIP_CALL( SCIPallocBufferArray(scip, &ind, nnonz) );
4691  SCIP_CALL( SCIPallocBufferArray(scip, &val, nnonz) );
4692 
4693  /* setup LP data */
4694  idx = 0;
4695  for( i = 0; i < nrows; ++i )
4696  {
4697  /* assemble corner point */
4698  SCIPdebugMsg(scip, "f(");
4699  for( j = 0; j < nvars; ++j )
4700  {
4701  /* if j'th bit of row index i is set, then take upper bound on var j, otherwise lower bound var j
4702  * we check this by shifting i for j positions to the right and checking whether the j'th bit is set */
4703  if( ((unsigned int)i >> j) & 0x1 )
4704  corner[j] = SCIPvarGetUbLocal(vars[j]);
4705  else
4706  corner[j] = SCIPvarGetLbLocal(vars[j]);
4707  SCIPdebugMsgPrint(scip, "%g, ", corner[j]);
4708  assert(!SCIPisInfinity(scip, REALABS(corner[j])));
4709  }
4710 
4711  /* evaluate function in current corner */
4712  SCIP_CALL( SCIPexprtreeEval(exprtree, corner, &funcval) );
4713  SCIPdebugMsgPrint(scip, ") = %g\n", funcval);
4714 
4715  if( !SCIPisFinite(funcval) || SCIPisInfinity(scip, REALABS(funcval)) )
4716  {
4717  SCIPdebugMsg(scip, "cannot compute underestimator for concave because constaint <%s> cannot be evaluated\n", SCIPconsGetName(cons));
4718  goto TERMINATE;
4719  }
4720 
4721  funcval *= treecoef;
4722 
4723  if( !doupper )
4724  {
4725  lhs[i] = -SCIPlpiInfinity(lpi);
4726  rhs[i] = funcval;
4727  }
4728  else
4729  {
4730  lhs[i] = funcval;
4731  rhs[i] = SCIPlpiInfinity(lpi);
4732  }
4733 
4734  /* add nonzeros of corner to matrix */
4735  beg[i] = idx;
4736  for( j = 0; j < nvars; ++j )
4737  {
4738  if( corner[j] != 0.0 )
4739  {
4740  ind[idx] = j;
4741  val[idx] = corner[j];
4742  ++idx;
4743  }
4744  }
4745 
4746  /* coefficient for constant is 1.0 */
4747  val[idx] = 1.0;
4748  ind[idx] = nvars;
4749  ++idx;
4750  }
4751  nnonz = idx;
4752  beg[nrows] = nnonz;
4753 
4754  for( j = 0; j < ncols; ++j )
4755  {
4756  lb[j] = -SCIPlpiInfinity(lpi);
4757  ub[j] = SCIPlpiInfinity(lpi);
4758  }
4759 
4760  /* objective coefficients are reference points, and an additional 1.0 for the constant */
4761  BMScopyMemoryArray(obj, ref, nvars);
4762  obj[nvars] = 1.0;
4763 
4764  /* get function value in reference point, so we can use this as a cutoff */
4765  SCIP_CALL( SCIPexprtreeEval(exprtree, ref, &funcval) );
4766  funcval *= treecoef;
4767 
4768  SCIP_CALL( SCIPlpiAddCols(lpi, ncols, obj, lb, ub, NULL, 0, NULL, NULL, NULL) );
4769  SCIP_CALL( SCIPlpiAddRows(lpi, nrows, lhs, rhs, NULL, nnonz, beg, ind, val) );
4770 
4771  /* make use of this convenient features, since for us nrows >> ncols */
4772  /*SCIP_CALL( SCIPlpiSetRealpar(lpi, SCIP_LPPAR_ROWREPSWITCH, 5.0) ); */
4773  /* get accurate coefficients */
4774  SCIP_CALL( SCIPlpiSetRealpar(lpi, SCIP_LPPAR_FEASTOL, SCIPfeastol(scip)/100.0) );
4775  SCIP_CALL( SCIPlpiSetRealpar(lpi, SCIP_LPPAR_OBJLIM, funcval) );
4776  SCIP_CALL( SCIPlpiSetIntpar(lpi, SCIP_LPPAR_LPITLIM, 10 * nvars) );
4777  SCIP_CALL( SCIPlpiSetIntpar(lpi, SCIP_LPPAR_SCALING, 1) );
4778  SCIP_CALL( SCIPlpiSetIntpar(lpi, SCIP_LPPAR_FROMSCRATCH, 1) );
4779 
4780  /* SCIPdebug( SCIP_CALL( SCIPlpiSetIntpar(lpi, SCIP_LPPAR_LPINFO, 1) ) ); */
4781 
4782  lpret = SCIPlpiSolveDual(lpi);
4783  if( lpret != SCIP_OKAY )
4784  {
4785  SCIPwarningMessage(scip, "solving auxiliary LP for underestimator of concave function returned %d\n", lpret);
4786  goto TERMINATE;
4787  }
4788 
4789  if( !SCIPlpiIsPrimalFeasible(lpi) )
4790  {
4791  SCIPdebugMsg(scip, "failed to find feasible solution for auxiliary LP for underestimator of concave function, iterlimexc = %u, cutoff = %u, unbounded = %u\n", SCIPlpiIsIterlimExc(lpi), SCIPlpiIsObjlimExc(lpi), SCIPlpiIsPrimalUnbounded(lpi));
4792  goto TERMINATE;
4793  }
4794  /* should be either solved to optimality, or the objective or iteration limit be hit */
4795  assert(SCIPlpiIsOptimal(lpi) || SCIPlpiIsObjlimExc(lpi) || SCIPlpiIsIterlimExc(lpi));
4796 
4797  /* setup row coefficient, reuse obj array to store LP sol values */
4798  SCIP_CALL( SCIPlpiGetSol(lpi, &lpobj, obj, NULL, NULL, NULL) );
4799 
4800  /* check that computed hyperplane is on right side of function in refpoint
4801  * if numerics is very bad (e.g., st_e32), then even this can happen */
4802  if( (!doupper && SCIPisFeasGT(scip, lpobj, funcval)) || (doupper && SCIPisFeasGT(scip, funcval, lpobj)) )
4803  {
4804  SCIPwarningMessage(scip, "computed cut does not underestimate concave function in refpoint\n");
4805  goto TERMINATE;
4806  }
4807  assert( doupper || SCIPisFeasLE(scip, lpobj, funcval) );
4808  assert(!doupper || SCIPisFeasLE(scip, funcval, lpobj) );
4809 
4810  /* add estimator to rowprep */
4811  SCIPaddRowprepConstant(rowprep, obj[nvars]);
4812  SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, nvars, vars, obj) );
4813 
4814  *success = TRUE;
4815 
4816 TERMINATE:
4817  SCIPfreeBufferArray(scip, &val);
4818  SCIPfreeBufferArray(scip, &ind);
4819  SCIPfreeBufferArray(scip, &beg);
4820  SCIPfreeBufferArray(scip, &rhs);
4821  SCIPfreeBufferArray(scip, &lhs);
4822  SCIPfreeBufferArray(scip, &ub);
4823  SCIPfreeBufferArray(scip, &lb);
4824  SCIPfreeBufferArray(scip, &obj);
4825 
4826  return SCIP_OKAY;
4827 }
4828 
4829 
4830 
4831 /** adds estimator of a constraints multivariate expression tree to a row
4832  * Given concave function f(x) and reference point ref.
4833  * Let (v_i: i=1,...,n) be corner points of current domain of x.
4834  * Find (coef,constant) such that <coef,v_i> + constant <= f(v_i) (cut validity) and
4835  * such that <coef, ref> + constant is maximized (cut efficacy).
4836  * Then <coef, x> + constant <= f(x) for all x in current domain.
4837  *
4838  * Similar to compute an overestimator for a convex function f(x).
4839  * Find (coef,constant) such that <coef,v_i> + constant >= f(v_i) and
4840  * such that <coef, ref> + constant is minimized.
4841  * Then <coef, x> + constant >= f(x) for all x in current domain.
4842  */
4843 static
4844 SCIP_RETCODE addConcaveEstimatorMultivariate(
4845  SCIP* scip, /**< SCIP data structure */
4846  SCIP_CONS* cons, /**< constraint */
4847  int exprtreeidx, /**< for which tree a secant should be added */
4848  SCIP_Real* ref, /**< reference values of expression tree variables where to generate cut */
4849  SCIP_ROWPREP* rowprep, /**< rowprep where to add estimator */
4850  SCIP_Bool* success /**< buffer to store whether a secant was succefully added to the row */
4851  )
4852 {
4853  SCIP_VAR** vars;
4854  SCIP_CONSDATA* consdata;
4855  SCIP_EXPRTREE* exprtree;
4856  SCIP_LPI* lpi;
4857  int nvars;
4858  int j;
4859  SCIP_Bool doupper;
4860 
4861  SCIP_RETCODE retcode;
4862 
4863  static SCIP_Bool warned_highdim_concave = FALSE;
4864 
4865  assert(scip != NULL);
4866  assert(cons != NULL);
4867  assert(ref != NULL);
4868  assert(rowprep != NULL);
4869  assert(success != NULL);
4870 
4871  consdata = SCIPconsGetData(cons);
4872  assert(consdata != NULL);
4873  assert(exprtreeidx >= 0);
4874  assert(exprtreeidx < consdata->nexprtrees);
4875  assert(consdata->exprtrees != NULL);
4876 
4877  exprtree = consdata->exprtrees[exprtreeidx];
4878  assert(exprtree != NULL);
4879 
4880  nvars = SCIPexprtreeGetNVars(exprtree);
4881  assert(nvars >= 2);
4882 
4883  *success = FALSE;
4884 
4885  /* size of LP is exponential in number of variables of tree, so do only for small trees */
4886  if( nvars > 10 )
4887  {
4888  if( !warned_highdim_concave )
4889  {
4890  SCIPwarningMessage(scip, "concave function in constraint <%s> too high-dimensional to compute underestimator\n", SCIPconsGetName(cons));
4891  warned_highdim_concave = TRUE;
4892  }
4893  return SCIP_OKAY;
4894  }
4895 
4896  vars = SCIPexprtreeGetVars(exprtree);
4897 
4898  /* check whether bounds are finite
4899  * make sure reference point is strictly within bounds
4900  * otherwise we can easily get an unbounded LP below, e.g., with instances like ex6_2_* from GlobalLib
4901  */
4902  for( j = 0; j < nvars; ++j )
4903  {
4904  if( SCIPisInfinity(scip, -SCIPvarGetLbLocal(vars[j])) || SCIPisInfinity(scip, SCIPvarGetUbLocal(vars[j])) )
4905  {
4906  SCIPdebugMsg(scip, "cannot compute underestimator for concave because variable <%s> is unbounded\n", SCIPvarGetName(vars[j]));
4907  return SCIP_OKAY;
4908  }
4909  assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(vars[j]), ref[j]));
4910  assert(SCIPisFeasGE(scip, SCIPvarGetUbLocal(vars[j]), ref[j]));
4911  ref[j] = MIN(SCIPvarGetUbLocal(vars[j]), MAX(SCIPvarGetLbLocal(vars[j]), ref[j])); /*lint !e666*/
4912  }
4913 
4914  /* create empty auxiliary LP and decide its objective sense */
4915  assert(consdata->curvatures[exprtreeidx] == SCIP_EXPRCURV_CONVEX || consdata->curvatures[exprtreeidx] == SCIP_EXPRCURV_CONCAVE);
4916  doupper = (consdata->curvatures[exprtreeidx] & SCIP_EXPRCURV_CONVEX); /*lint !e641*/
4917  SCIP_CALL( SCIPlpiCreate(&lpi, SCIPgetMessagehdlr(scip), "concaveunderest", doupper ? SCIP_OBJSEN_MINIMIZE : SCIP_OBJSEN_MAXIMIZE) );
4918  if( lpi == NULL )
4919  {
4920  SCIPerrorMessage("failed to create auxiliary LP\n");
4921  return SCIP_ERROR;
4922  }
4923 
4924  /* capture the retcode, free the LPI afterwards */
4925  retcode = _addConcaveEstimatorMultivariate(scip, lpi, cons, exprtreeidx, ref, rowprep, vars, exprtree, nvars, doupper, success);
4926 
4927  assert(lpi != NULL);
4928  SCIP_CALL( SCIPlpiFree(&lpi) );
4929 
4930  return retcode;
4931 }
4932 
4933 /** Computes the linear coeffs and the constant in a linear expression
4934  * both scaled by a given scalar value.
4935  * The coeffs of the variables will be stored in the given array at
4936  * their variable index.
4937  * The constant of the given linear expression will be added to the given
4938  * buffer.
4939  */
4940 static
4941 SCIP_RETCODE getCoeffsAndConstantFromLinearExpr(
4942  SCIP_EXPR* expr, /**< the linear expression */
4943  SCIP_Real scalar, /**< the scalar value, i.e. the coeff of the given expression */
4944  SCIP_Real* varcoeffs, /**< buffer array to store the computed coefficients */
4945  SCIP_Real* constant /**< buffer to hold the constant value of the given expression */
4946  )
4947 {
4948  switch( SCIPexprGetOperator( expr ) )
4949  {
4950  case SCIP_EXPR_VARIDX: /* set coeff for this variable to current scalar */
4951  {
4952  /* TODO: can a linear expression contain the same variable twice?
4953  * if yes varcoeffs need to be initialized to zero before calling this function
4954  * and coeff must not be overridden but summed up instead. */
4955  varcoeffs[SCIPexprGetOpIndex( expr )] = scalar;
4956  return SCIP_OKAY;
4957  }
4958 
4959  case SCIP_EXPR_CONST:
4960  {
4961  /* constant value increases */
4962  *constant += scalar * SCIPexprGetOpReal( expr );
4963  return SCIP_OKAY;
4964  }
4965 
4966  case SCIP_EXPR_MUL: /* need to find the constant part of the muliplication and then recurse */
4967  {
4968  SCIP_EXPR** children;
4969  children = SCIPexprGetChildren( expr );
4970 
4971  /* first part is constant */
4972  if( SCIPexprGetOperator( children[0] ) == SCIP_EXPR_CONST )
4973  {
4974  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[1], scalar * SCIPexprGetOpReal( children[0] ), varcoeffs, constant ) );
4975  return SCIP_OKAY;
4976  }
4977 
4978  /* second part is constant */
4979  if( SCIPexprGetOperator( children[1] ) == SCIP_EXPR_CONST )
4980  {
4981  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[0], scalar * SCIPexprGetOpReal( children[1] ), varcoeffs, constant ) );
4982  return SCIP_OKAY;
4983  }
4984 
4985  /* nonlinear -> break out to error case */
4986  break;
4987  }
4988 
4989  case SCIP_EXPR_PLUS: /* just recurse */
4990  {
4991  SCIP_EXPR** children;
4992  children = SCIPexprGetChildren( expr );
4993  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[0], scalar, varcoeffs, constant ) );
4994  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[1], scalar, varcoeffs, constant ) );
4995  return SCIP_OKAY;
4996  }
4997 
4998  case SCIP_EXPR_MINUS: /* recursion on second child is called with negated scalar */
4999  {
5000  SCIP_EXPR** children;
5001  children = SCIPexprGetChildren( expr );
5002  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[0], scalar, varcoeffs, constant ) );
5003  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[1], -scalar, varcoeffs, constant ) );
5004  return SCIP_OKAY;
5005  }
5006 
5007  case SCIP_EXPR_SUM: /* just recurse */
5008  {
5009  SCIP_EXPR** children;
5010  int nchildren;
5011  int c;
5012 
5013  children = SCIPexprGetChildren(expr);
5014  nchildren = SCIPexprGetNChildren(expr);
5015 
5016  for( c = 0; c < nchildren; ++c )
5017  {
5018  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[c], scalar, varcoeffs, constant ) );
5019  }
5020 
5021  return SCIP_OKAY;
5022  }
5023 
5024  case SCIP_EXPR_LINEAR: /* add scaled constant and recurse on children with their coeff multiplied into scalar */
5025  {
5026  SCIP_Real* childCoeffs;
5027  SCIP_EXPR** children;
5028  int i;
5029 
5030  *constant += scalar * SCIPexprGetLinearConstant( expr );
5031 
5032  children = SCIPexprGetChildren( expr );
5033  childCoeffs = SCIPexprGetLinearCoefs( expr );
5034 
5035  for( i = 0; i < SCIPexprGetNChildren( expr ); ++i )
5036  {
5037  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[i], scalar * childCoeffs[i], varcoeffs, constant ) );
5038  }
5039 
5040  return SCIP_OKAY;
5041  }
5042 
5043  default:
5044  break;
5045  } /*lint !e788*/
5046 
5047  SCIPerrorMessage( "Cannot extract linear coefficients from expressions with operator %d %s\n", SCIPexprGetOperator( expr ), SCIPexpropGetName(SCIPexprGetOperator( expr )));
5048  SCIPABORT();
5049  return SCIP_ERROR; /*lint !e527*/
5050 }
5051 
5052 /** adds estimator from user callback of a constraints user expression tree to a row
5053  */
5054 static
5055 SCIP_RETCODE addUserEstimator(
5056  SCIP* scip, /**< SCIP data structure */
5057  SCIP_CONS* cons, /**< constraint */
5058  int exprtreeidx, /**< for which tree an estimator should be added */
5059  SCIP_Real* x, /**< value of expression tree variables where to generate cut */
5060  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
5061  SCIP_ROWPREP* rowprep, /**< rowprep where to add estimator */
5062  SCIP_Bool* success /**< buffer to store whether an estimator was succefully added to the rowprep */
5063  )
5064 {
5065  SCIP_CONSDATA* consdata;
5066  SCIP_EXPRTREE* exprtree;
5067  SCIP_EXPR** children;
5068  SCIP_VAR** vars;
5069  SCIP_Real* params;
5070  SCIP_INTERVAL* varbounds;
5071 
5072  SCIP_INTERVAL* childbounds;
5073  SCIP_Real* childvals;
5074  SCIP_Real* childcoeffs;
5075 
5076  SCIP_Real constant;
5077  SCIP_Real treecoef;
5078  int nvars;
5079  int nchildren;
5080  int i;
5081 
5082  consdata = SCIPconsGetData( cons );
5083  assert( consdata != NULL );
5084  assert( exprtreeidx >= 0 );
5085  assert( exprtreeidx < consdata->nexprtrees );
5086  assert( consdata->exprtrees != NULL );
5087  assert( rowprep != NULL );
5088  assert( success != NULL );
5089 
5090  exprtree = consdata->exprtrees[exprtreeidx];
5091  assert( exprtree != NULL );
5092  assert( SCIPexprGetOperator(SCIPexprtreeGetRoot(exprtree)) == SCIP_EXPR_USER );
5093 
5094  /* if user did not implement estimator callback, then we cannot do anything */
5095  if( !SCIPexprHasUserEstimator(SCIPexprtreeGetRoot(exprtree)) )
5096  {
5097  *success = FALSE;
5098  return SCIP_OKAY;
5099  }
5100 
5101  params = SCIPexprtreeGetParamVals( exprtree );
5102  nvars = SCIPexprtreeGetNVars( exprtree );
5103  vars = SCIPexprtreeGetVars( exprtree );
5104  nchildren = SCIPexprGetNChildren( SCIPexprtreeGetRoot( exprtree ) );
5105  children = SCIPexprGetChildren( SCIPexprtreeGetRoot( exprtree ) );
5106 
5107  /* Get bounds of variables */
5108  SCIP_CALL( SCIPallocBufferArray( scip, &varbounds, nchildren ) );
5109 
5110  for( i = 0; i < nvars; ++i )
5111  {
5112  double lb = SCIPvarGetLbLocal( vars[i] );
5113  double ub = SCIPvarGetUbLocal( vars[i] );
5114  SCIPintervalSetBounds( &varbounds[i],
5115  -infty2infty( SCIPinfinity( scip ), INTERVALINFTY, -MIN( lb, ub ) ),
5116  +infty2infty( SCIPinfinity( scip ), INTERVALINFTY, MAX( lb, ub ) ) );
5117  }
5118 
5119  /* Compute bounds and solution value for the user expressions children */
5120  SCIP_CALL( SCIPallocBufferArray( scip, &childcoeffs, nchildren ) );
5121  SCIP_CALL( SCIPallocBufferArray( scip, &childbounds, nchildren ) );
5122  SCIP_CALL( SCIPallocBufferArray( scip, &childvals, nchildren ) );
5123 
5124  for( i = 0; i < nchildren; ++i )
5125  {
5126  SCIP_CALL( SCIPexprEval( children[i], x, params, &childvals[i] ) );
5127  SCIP_CALL( SCIPexprEvalInt( children[i], INTERVALINFTY, varbounds, params, &childbounds[i] ) );
5128  }
5129 
5130  /* varbounds not needed any longer */
5131  SCIPfreeBufferArray( scip, &varbounds );
5132 
5133  /* call estimator for user expressions to compute coeffs and constant for the user expressions children */
5134  SCIP_CALL( SCIPexprEstimateUser( SCIPexprtreeGetRoot( exprtree ), INTERVALINFTY, childvals, childbounds, overestimate, childcoeffs, &constant, success ) );
5135 
5136  if( *success )
5137  {
5138  SCIP_Real* varcoeffs;
5139  SCIP_CALL( SCIPallocBufferArray( scip, &varcoeffs, nvars ) );
5140 
5141  treecoef = consdata->nonlincoefs[exprtreeidx];
5142  constant *= treecoef;
5143 
5144  for( i = 0; i < nchildren; ++i )
5145  {
5146  SCIP_CALL( getCoeffsAndConstantFromLinearExpr( children[i], childcoeffs[i]*treecoef, varcoeffs, &constant ) );
5147  }
5148 
5149  SCIPaddRowprepConstant(rowprep, constant);
5150  SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, nvars, vars, varcoeffs) );
5151 
5152  SCIPfreeBufferArray( scip, &varcoeffs );
5153  }
5154 
5155  SCIPfreeBufferArray( scip, &childcoeffs );
5156  SCIPfreeBufferArray( scip, &childbounds );
5157  SCIPfreeBufferArray( scip, &childvals );
5158 
5159  return SCIP_OKAY;
5160 }
5161 
5162 /** adds estimator from interval gradient of a constraints univariate expression tree to a row
5163  * a reference point is used to decide in which corner to generate the cut
5164  */
5165 static
5166 SCIP_RETCODE addIntervalGradientEstimator(
5167  SCIP* scip, /**< SCIP data structure */
5168  SCIP_EXPRINT* exprint, /**< expression interpreter */
5169  SCIP_CONS* cons, /**< constraint */
5170  int exprtreeidx, /**< for which tree a secant should be added */
5171  SCIP_Real* x, /**< value of expression tree variables where to generate cut */
5172  SCIP_Bool newx, /**< whether the last evaluation of the expression with the expression interpreter was not at x */
5173  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
5174  SCIP_ROWPREP* rowprep, /**< rowprep where to add estimator */
5175  SCIP_Bool* success /**< buffer to store whether an estimator was succefully added to the rowprep */
5176  )
5177 {
5178  SCIP_CONSDATA* consdata;
5179  SCIP_EXPRTREE* exprtree;
5180  SCIP_Real treecoef;
5181  SCIP_Real* coefs;
5182  SCIP_Real constant;
5183  SCIP_Real val;
5184  SCIP_Real lb;
5185  SCIP_Real ub;
5186  SCIP_INTERVAL* box;
5187  SCIP_INTERVAL* intgrad;
5188  SCIP_INTERVAL intval;
5189  SCIP_VAR** vars;
5190  int nvars;
5191  int i;
5192 
5193  assert(scip != NULL);
5194  assert(cons != NULL);
5195  assert(x != NULL);
5196  assert(rowprep != NULL);
5197  assert(success != NULL);
5198 
5199  consdata = SCIPconsGetData(cons);
5200  assert(consdata != NULL);
5201  assert(exprtreeidx >= 0);
5202  assert(exprtreeidx < consdata->nexprtrees);
5203  assert(consdata->exprtrees != NULL);
5204 
5205  exprtree = consdata->exprtrees[exprtreeidx];
5206  assert(exprtree != NULL);
5207  assert(newx || SCIPexprtreeGetInterpreterData(exprtree) != NULL);
5208 
5209  *success = FALSE;
5210 
5211  /* skip interval gradient if expression interpreter cannot compute interval gradients */
5212  if( !(SCIPexprintGetCapability() & SCIP_EXPRINTCAPABILITY_INTGRADIENT) )
5213  return SCIP_OKAY;
5214 
5215  nvars = SCIPexprtreeGetNVars(exprtree);
5216  vars = SCIPexprtreeGetVars(exprtree);
5217 
5218  box = NULL;
5219  intgrad = NULL;
5220  coefs = NULL;
5221 
5222  SCIP_CALL( SCIPallocBufferArray(scip, &box, nvars) );
5223 
5224  /* move reference point to bounds, setup box */
5225  for( i = 0; i < nvars; ++i )
5226  {
5227  lb = SCIPvarGetLbLocal(vars[i]);
5228  ub = SCIPvarGetUbLocal(vars[i]);
5229  if( SCIPisInfinity(scip, -lb) )
5230  {
5231  if( SCIPisInfinity(scip, ub) )
5232  {
5233  SCIPdebugMsg(scip, "skip interval gradient estimator for constraint <%s> because variable <%s> is still unbounded.\n", SCIPconsGetName(cons), SCIPvarGetName(vars[i]));
5234  goto INTGRADESTIMATOR_CLEANUP;
5235  }
5236  x[i] = ub;
5237  }
5238  else
5239  {
5240  if( SCIPisInfinity(scip, ub) )
5241  x[i] = lb;
5242  else
5243  x[i] = (2.0*x[i] < lb+ub) ? lb : ub;
5244  }
5245  SCIPintervalSetBounds(&box[i],
5246  -infty2infty(SCIPinfinity(scip), INTERVALINFTY, -MIN(lb, ub)),
5247  +infty2infty(SCIPinfinity(scip), INTERVALINFTY, MAX(lb, ub)));
5248  }
5249 
5250  /* compile expression if evaluated the first time; can only happen if newx is FALSE */
5251  if( newx && SCIPexprtreeGetInterpreterData(exprtree) == NULL )
5252  {
5253  SCIP_CALL( SCIPexprintCompile(exprint, exprtree) );
5254  }
5255 
5256  /* evaluate in reference point */
5257  SCIP_CALL( SCIPexprintEval(exprint, exprtree, x, &val) );
5258  if( !SCIPisFinite(val) )
5259  {
5260  SCIPdebugMsg(scip, "Got nonfinite function value from evaluation of constraint %s tree %d. skipping interval gradient estimator.\n", SCIPconsGetName(cons), exprtreeidx);
5261  goto INTGRADESTIMATOR_CLEANUP;
5262  }
5263 
5264  treecoef = consdata->nonlincoefs[exprtreeidx];
5265  val *= treecoef;
5266  constant = val;
5267 
5268  /* compute interval gradient */
5269  SCIP_CALL( SCIPallocBufferArray(scip, &intgrad, nvars) );
5270  SCIP_CALL( SCIPexprintGradInt(exprint, exprtree, INTERVALINFTY, box, TRUE, &intval, intgrad) );
5271  SCIPintervalMulScalar(INTERVALINFTY, &intval, intval, treecoef);
5272 
5273  /* printf("nvars %d side %d xref = %g x = [%g,%g] intval = [%g,%g] intgrad = [%g,%g]\n", nvars, side, x[0],
5274  box[0].inf, box[0].sup, intval.inf, intval.sup, intgrad[0].inf, intgrad[0].sup); */
5275 
5276  /* compute coefficients and constant */
5277  SCIP_CALL( SCIPallocBufferArray(scip, &coefs, nvars) );
5278  for( i = 0; i < nvars; ++i )
5279  {
5280  val = x[i];
5281  lb = SCIPintervalGetInf(box[i]);
5282  ub = SCIPintervalGetSup(box[i]);
5283 
5284  SCIPintervalMulScalar(INTERVALINFTY, &intgrad[i], intgrad[i], treecoef);
5285 
5286  if( SCIPisEQ(scip, lb, ub) )
5287  coefs[i] = 0.0;
5288  else if( (overestimate && val == ub) || /*lint !e777*/
5289  (!overestimate && val == lb) ) /*lint !e777*/
5290  coefs[i] = SCIPintervalGetInf(intgrad[i]);
5291  else
5292  coefs[i] = SCIPintervalGetSup(intgrad[i]);
5293 
5294  if( SCIPisZero(scip, coefs[i]) )
5295  continue;
5296 
5297  if( SCIPisInfinity(scip, -coefs[i]) || SCIPisInfinity(scip, coefs[i]) )
5298  {
5299  SCIPdebugMsg(scip, "skip intgrad estimator because of infinite interval bound\n");
5300  goto INTGRADESTIMATOR_CLEANUP;
5301  }
5302 
5303  constant -= coefs[i] * val;
5304  }
5305 
5306  /* add interval gradient estimator to row */
5307  SCIPaddRowprepConstant(rowprep, constant);
5308  SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, nvars, vars, coefs) );
5309 
5310  INTGRADESTIMATOR_CLEANUP:
5311  SCIPfreeBufferArrayNull(scip, &coefs);
5312  SCIPfreeBufferArrayNull(scip, &intgrad);
5313  SCIPfreeBufferArrayNull(scip, &box);
5314 
5315  return SCIP_OKAY;
5316 }
5317 
5318 /** generates a cut based on linearization (if convex), secant (if concave), or intervalgradient (if indefinite)
5319  */
5320 static
5321 SCIP_RETCODE generateCut(
5322  SCIP* scip, /**< SCIP data structure */
5323  SCIP_EXPRINT* exprint, /**< expression interpreter */
5324  SCIP_CONS* cons, /**< constraint */
5325  SCIP_Real** ref, /**< reference point for each exprtree, or NULL if sol should be used */
5326  SCIP_SOL* sol, /**< reference solution where cut should be generated, or NULL if LP solution should be used */
5327  SCIP_Bool newsol, /**< whether the last evaluation of the expression with the expression interpreter was not at sol */
5328  SCIP_SIDETYPE side, /**< for which side a cut should be generated */
5329  SCIP_ROW** row, /**< storage for cut */
5330  SCIP_Real minviol, /**< minimal absolute violation we try to achieve */
5331  SCIP_Real maxrange, /**< maximal range allowed */
5332  SCIP_Bool expensivecurvchecks,/**< whether also expensive checks should be executed */
5333  SCIP_Bool assumeconvex /**< whether to assume convexity in inequalities */
5334  )
5335 {
5336  SCIP_ROWPREP* rowprep;
5337  SCIP_CONSDATA* consdata;
5338  SCIP_Bool success;
5339  SCIP_Real* x;
5340  int i;
5341 
5342  assert(scip != NULL);
5343  assert(cons != NULL);
5344  assert(row != NULL);
5345 
5346  SCIPdebugMsg(scip, "constructing cut for %s hand side of constraint <%s>\n", side == SCIP_SIDETYPE_LEFT ? "left" : "right", SCIPconsGetName(cons));
5347 
5348  SCIP_CALL( checkCurvature(scip, cons, expensivecurvchecks, assumeconvex) );
5349 
5350  consdata = SCIPconsGetData(cons);
5351  assert(consdata != NULL);
5352 
5353  if( consdata->nexprtrees == 0 )
5354  {
5355  char rowname[SCIP_MAXSTRLEN];
5356  (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "%s_%u", SCIPconsGetName(cons), ++(consdata->ncuts));
5357 
5358  /* if we are actually linear, add the constraint as row to the LP */
5359  SCIP_CALL( SCIPcreateEmptyRowCons(scip, row, SCIPconsGetHdlr(cons), rowname, consdata->lhs, consdata->rhs, SCIPconsIsLocal(cons), FALSE , TRUE) );
5360  SCIP_CALL( SCIPaddVarsToRow(scip, *row, consdata->nlinvars, consdata->linvars, consdata->lincoefs) );
5361  return SCIP_OKAY;
5362  }
5363 
5364  SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, side,
5365  !(side == SCIP_SIDETYPE_LEFT && (consdata->curvature & SCIP_EXPRCURV_CONCAVE)) &&
5366  !(side == SCIP_SIDETYPE_RIGHT && (consdata->curvature & SCIP_EXPRCURV_CONVEX ))) );
5367  SCIPaddRowprepSide(rowprep, side == SCIP_SIDETYPE_LEFT ? consdata->lhs : consdata->rhs);
5368  (void) SCIPsnprintf(rowprep->name, SCIP_MAXSTRLEN, "%s_%u", SCIPconsGetName(cons), ++(consdata->ncuts));
5369 
5370  if( ref == NULL )
5371  {
5372  SCIP_CALL( SCIPallocBufferArray(scip, &x, SCIPexprtreeGetNVars(consdata->exprtrees[0])) );
5373  }
5374 
5375  success = TRUE;
5376  for( i = 0; i < consdata->nexprtrees; ++i )
5377  {
5378  if( ref == NULL )
5379  {
5380  SCIP_CALL( SCIPreallocBufferArray(scip, &x, SCIPexprtreeGetNVars(consdata->exprtrees[i])) ); /*lint !e644*/
5381  SCIP_CALL( SCIPgetSolVals(scip, sol, SCIPexprtreeGetNVars(consdata->exprtrees[i]), SCIPexprtreeGetVars(consdata->exprtrees[i]), x) );
5382  }
5383  else
5384  {
5385  x = ref[i];
5386  }
5387 
5388  if( (side == SCIP_SIDETYPE_LEFT && (consdata->curvatures[i] & SCIP_EXPRCURV_CONCAVE)) ||
5389  (side == SCIP_SIDETYPE_RIGHT && (consdata->curvatures[i] & SCIP_EXPRCURV_CONVEX )) )
5390  {
5391  SCIP_CALL( addLinearization(scip, exprint, cons, i, x, newsol, rowprep, &success) );
5392  }
5393  else if( (side == SCIP_SIDETYPE_LEFT && (consdata->curvatures[i] & SCIP_EXPRCURV_CONVEX)) ||
5394  ( side == SCIP_SIDETYPE_RIGHT && (consdata->curvatures[i] & SCIP_EXPRCURV_CONCAVE)) )
5395  {
5396  switch( SCIPexprtreeGetNVars(consdata->exprtrees[i]) )
5397  {
5398  case 1:
5399  SCIP_CALL( addConcaveEstimatorUnivariate(scip, cons, i, rowprep, &success) );
5400  break;
5401 
5402  case 2:
5403  SCIP_CALL( addConcaveEstimatorBivariate(scip, cons, i, x, rowprep, &success) );
5404  break;
5405 
5406  default:
5407  SCIP_CALL( addConcaveEstimatorMultivariate(scip, cons, i, x, rowprep, &success) );
5408  break;
5409  }
5410  if( !success )
5411  {
5412  SCIPdebugMsg(scip, "failed to generate polyhedral estimator for %d-dim concave function in exprtree %d, fall back to intervalgradient cut\n", SCIPexprtreeGetNVars(consdata->exprtrees[i]), i);
5413  SCIP_CALL( addIntervalGradientEstimator(scip, exprint, cons, i, x, newsol, side == SCIP_SIDETYPE_LEFT, rowprep, &success) );
5414  }
5415  }
5416  else if( SCIPexprGetOperator( SCIPexprtreeGetRoot( consdata->exprtrees[i] ) ) == SCIP_EXPR_USER )
5417  {
5418  SCIP_CALL( addUserEstimator( scip, cons, i, x, side == SCIP_SIDETYPE_LEFT, rowprep, &success ) );
5419  if( !success ) /* the user estimation may not be implemented -> try interval estimator */
5420  {
5421  SCIP_CALL( addIntervalGradientEstimator(scip, exprint, cons, i, x, newsol, side == SCIP_SIDETYPE_LEFT, rowprep, &success) );
5422  }
5423  }
5424  else
5425  {
5426  SCIP_CALL( addIntervalGradientEstimator(scip, exprint, cons, i, x, newsol, side == SCIP_SIDETYPE_LEFT, rowprep, &success) );
5427  }
5428 
5429  if( !success )
5430  break;
5431  }
5432 
5433  if( ref == NULL )
5434  {
5435  SCIPfreeBufferArray(scip, &x);
5436  }
5437 
5438  if( success )
5439  {
5440  SCIP_Real coefrange;
5441 
5442  /* add coefficients for linear variables */
5443  SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, consdata->nlinvars, consdata->linvars, consdata->lincoefs) );
5444 
5445  /* merge terms in same variable */
5446  SCIPmergeRowprepTerms(scip, rowprep);
5447 
5448  /* cleanup row */
5449  SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, sol, maxrange, minviol, &coefrange, NULL) );
5450 
5451  /* check that coefficient range is ok */
5452  success = coefrange <= maxrange;
5453 
5454  /* check that side is finite */ /*lint --e{514} */
5455  success &= !SCIPisInfinity(scip, REALABS(rowprep->side));
5456 
5457  /* check whether maximal coef is finite, if any */
5458  success &= (rowprep->nvars == 0) || !SCIPisInfinity(scip, REALABS(rowprep->coefs[0]));
5459  }
5460 
5461  if( success )
5462  {
5463  SCIP_CALL( SCIPgetRowprepRowCons(scip, row, rowprep, SCIPconsGetHdlr(cons)) );
5464  }
5465  else
5466  *row = NULL;
5467 
5468  SCIPfreeRowprep(scip, &rowprep);
5469 
5470  return SCIP_OKAY;
5471 }
5472 
5473 /** tries to separate solution or LP solution by a linear cut
5474  *
5475  * assumes that constraint violations have been computed
5476  */
5477 static
5478 SCIP_RETCODE separatePoint(
5479  SCIP* scip, /**< SCIP data structure */
5480  SCIP_CONSHDLR* conshdlr, /**< nonlinear constraints handler */
5481  SCIP_CONS** conss, /**< constraints */
5482  int nconss, /**< number of constraints */
5483  int nusefulconss, /**< number of constraints that seem to be useful */
5484  SCIP_SOL* sol, /**< solution to separate, or NULL if LP solution should be used */
5485  SCIP_Bool newsol, /**< have the constraints just been evaluated at this point? */
5486  SCIP_Real minefficacy, /**< minimal efficacy of a cut if it should be added to the LP */
5487  SCIP_Bool inenforcement, /**< whether we are in constraint enforcement */
5488  SCIP_RESULT* result, /**< result of separation */
5489  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 */
5490  )
5491 {
5492  SCIP_CONSHDLRDATA* conshdlrdata;
5493  SCIP_CONSDATA* consdata;
5494  SCIP_Real efficacy;
5495  SCIP_Real feasibility;
5496  SCIP_SIDETYPE violside;
5497  int c;
5498  SCIP_ROW* row;
5499 
5500  assert(scip != NULL);
5501  assert(conshdlr != NULL);
5502  assert(conss != NULL || nconss == 0);
5503  assert(nusefulconss <= nconss);
5504  assert(result != NULL);
5505 
5506  *result = SCIP_FEASIBLE;
5507 
5508  conshdlrdata = SCIPconshdlrGetData(conshdlr);
5509  assert(conshdlrdata != NULL);
5510 
5511  if( bestefficacy != NULL )
5512  *bestefficacy = 0.0;
5513 
5514  for( c = 0, violside = SCIP_SIDETYPE_LEFT; c < nconss; c = (violside == SCIP_SIDETYPE_RIGHT ? c+1 : c), violside = (violside == SCIP_SIDETYPE_LEFT ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT) )
5515  {
5516  assert(conss != NULL);
5517 
5518  /* skip constraints that are not enabled */
5519  if( !SCIPconsIsEnabled(conss[c]) || SCIPconsIsDeleted(conss[c]) )
5520  continue;
5521  assert(SCIPconsIsActive(conss[c]));
5522 
5523  consdata = SCIPconsGetData(conss[c]);
5524  assert(consdata != NULL);
5525