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