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-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file cons_nonlinear.c
17  * @ingroup DEFPLUGINS_CONS
18  * @brief constraint handler for nonlinear constraints specified by algebraic expressions
19  * @author Ksenia Bestuzheva
20  * @author Benjamin Mueller
21  * @author Felipe Serrano
22  * @author Stefan Vigerske
23  */
24 
25 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
26 
27 #ifdef SCIP_DEBUG
28 #define ENFO_LOGGING
29 #endif
30 
31 /* enable to get log output for enforcement */
32 /* #define ENFO_LOGGING */
33 /* define to get enforcement logging into file */
34 /* #define ENFOLOGFILE "consexpr_enfo.log" */
35 
36 /* define to get more debug output from domain propagation */
37 /* #define DEBUG_PROP */
38 
39 /*lint -e440*/
40 /*lint -e441*/
41 /*lint -e528*/
42 /*lint -e666*/
43 /*lint -e777*/
44 /*lint -e866*/
45 
46 #include <assert.h>
47 #include <ctype.h>
48 
49 #include "scip/cons_nonlinear.h"
50 #include "scip/nlhdlr.h"
51 #include "scip/expr_var.h"
52 #include "scip/expr_sum.h"
53 #include "scip/expr_value.h"
54 #include "scip/expr_pow.h"
55 #include "scip/nlhdlr_convex.h"
56 #include "scip/cons_linear.h"
57 #include "scip/cons_varbound.h"
58 #include "scip/cons_and.h"
60 #include "scip/heur_subnlp.h"
61 #include "scip/heur_trysol.h"
62 #include "scip/nlpi_ipopt.h" /* for SCIPsolveLinearEquationsIpopt */
63 #include "scip/debug.h"
64 #include "scip/dialog_default.h"
65 
66 /* fundamental constraint handler properties */
67 #define CONSHDLR_NAME "nonlinear"
68 #define CONSHDLR_DESC "handler for nonlinear constraints specified by algebraic expressions"
69 #define CONSHDLR_ENFOPRIORITY -60 /**< priority of the constraint handler for constraint enforcing */
70 #define CONSHDLR_CHECKPRIORITY -4000010 /**< priority of the constraint handler for checking feasibility */
71 #define CONSHDLR_EAGERFREQ 100 /**< frequency for using all instead of only the useful constraints in separation,
72  * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
73 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
74 
75 /* optional constraint handler properties */
76 #define CONSHDLR_SEPAPRIORITY 10 /**< priority of the constraint handler for separation */
77 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
78 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
79 
80 #define CONSHDLR_PROPFREQ 1 /**< frequency for propagating domains; zero means only preprocessing propagation */
81 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
82 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP /**< propagation timing mask of the constraint handler*/
83 
84 #define CONSHDLR_PRESOLTIMING SCIP_PRESOLTIMING_ALWAYS /**< presolving timing of the constraint handler (fast, medium, or exhaustive) */
85 #define CONSHDLR_MAXPREROUNDS -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
86 
87 /* properties of the nonlinear constraint handler statistics table */
88 #define TABLE_NAME_NONLINEAR "cons_nonlinear"
89 #define TABLE_DESC_NONLINEAR "nonlinear constraint handler statistics"
90 #define TABLE_POSITION_NONLINEAR 14600 /**< the position of the statistics table */
91 #define TABLE_EARLIEST_STAGE_NONLINEAR SCIP_STAGE_TRANSFORMED /**< output of the statistics table is only printed from this stage onwards */
92 
93 /* properties of the nonlinear handler statistics table */
94 #define TABLE_NAME_NLHDLR "nlhdlr"
95 #define TABLE_DESC_NLHDLR "nonlinear handler statistics"
96 #define TABLE_POSITION_NLHDLR 14601 /**< the position of the statistics table */
97 #define TABLE_EARLIEST_STAGE_NLHDLR SCIP_STAGE_PRESOLVING /**< output of the statistics table is only printed from this stage onwards */
98 
99 #define DIALOG_NAME "nlhdlrs"
100 #define DIALOG_DESC "display nonlinear handlers"
101 #define DIALOG_ISSUBMENU FALSE
103 #define VERTEXPOLY_MAXPERTURBATION 1e-3 /**< maximum perturbation */
104 #define VERTEXPOLY_USEDUALSIMPLEX TRUE /**< use dual or primal simplex algorithm? */
105 #define VERTEXPOLY_RANDNUMINITSEED 20181029 /**< seed for random number generator, which is used to move points away from the boundary */
106 #define VERTEXPOLY_ADJUSTFACETFACTOR 1e1 /**< adjust resulting facets in checkRikun() up to a violation of this value times lpfeastol */
108 #define BRANCH_RANDNUMINITSEED 20191229 /**< seed for random number generator, which is used to select from several similar good branching candidates */
110 #define BILIN_MAXNAUXEXPRS 10 /**< maximal number of auxiliary expressions per bilinear term */
112 /** translate from one value of infinity to another
113  *
114  * if val is &ge; infty1, then give infty2, else give val
115  */
116 #define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))
118 /** translates x to 2^x for non-negative integer x */
119 #define POWEROFTWO(x) (0x1u << (x))
121 #ifdef ENFO_LOGGING
122 #define ENFOLOG(x) if( SCIPgetSubscipDepth(scip) == 0 && SCIPgetVerbLevel(scip) >= SCIP_VERBLEVEL_NORMAL ) { x }
123 FILE* enfologfile = NULL;
124 #else
125 #define ENFOLOG(x)
126 #endif
127 
128 /*
129  * Data structures
130  */
131 
132 /** enforcement data of an expression */
133 typedef struct
134 {
135  SCIP_NLHDLR* nlhdlr; /**< nonlinear handler */
136  SCIP_NLHDLREXPRDATA* nlhdlrexprdata; /**< data of nonlinear handler */
137  SCIP_NLHDLR_METHOD nlhdlrparticipation; /**< methods where nonlinear handler participates */
138  SCIP_Bool issepainit; /**< was the initsepa callback of nlhdlr called */
139  SCIP_Real auxvalue; /**< auxiliary value of expression w.r.t. currently enforced solution */
140  SCIP_Bool sepabelowusesactivity;/**< whether sepabelow uses activity of some expression */
141  SCIP_Bool sepaaboveusesactivity;/**< whether sepaabove uses activity of some expression */
143 
144 /** data stored by constraint handler in an expression that belongs to a nonlinear constraint */
145 struct SCIP_Expr_OwnerData
146 {
147  SCIP_CONSHDLR* conshdlr; /** nonlinear constraint handler */
148 
149  /* locks and monotonicity */
150  int nlockspos; /**< positive locks counter */
151  int nlocksneg; /**< negative locks counter */
152  SCIP_MONOTONE* monotonicity; /**< array containing monotonicity of expression w.r.t. each child */
153  int monotonicitysize; /**< length of monotonicity array */
154 
155  /* propagation (in addition to activity that is stored in expr) */
156  SCIP_INTERVAL propbounds; /**< bounds to propagate in reverse propagation */
157  unsigned int propboundstag; /**< tag to indicate whether propbounds are valid for the current propagation rounds */
158  SCIP_Bool inpropqueue; /**< whether expression is queued for propagation */
159 
160  /* enforcement of expr == auxvar (or expr <= auxvar, or expr >= auxvar) */
161  EXPRENFO** enfos; /**< enforcements */
162  int nenfos; /**< number of enforcements, or -1 if not initialized */
163  unsigned int lastenforced; /**< last enforcement round where expression was enforced successfully */
164  unsigned int nactivityusesprop; /**< number of nonlinear handlers whose activity computation (or domain propagation) depends on the activity of the expression */
165  unsigned int nactivityusessepa; /**< number of nonlinear handlers whose separation (estimate or enfo) depends on the activity of the expression */
166  unsigned int nauxvaruses; /**< number of nonlinear handlers whose separation uses an auxvar in the expression */
167  SCIP_VAR* auxvar; /**< auxiliary variable used for outer approximation cuts */
168 
169  /* branching */
170  SCIP_Real violscoresum; /**< sum of violation scores for branching stored for this expression */
171  SCIP_Real violscoremax; /**< max of violation scores for branching stored for this expression */
172  int nviolscores; /**< number of violation scores stored for this expression */
173  unsigned int violscoretag; /**< tag to decide whether a violation score of an expression needs to be initialized */
174 
175  /* additional data for variable expressions (TODO move into sub-struct?) */
176  SCIP_CONS** conss; /**< constraints in which this variable appears */
177  int nconss; /**< current number of constraints in conss */
178  int consssize; /**< length of conss array */
179  SCIP_Bool consssorted; /**< is the array of constraints sorted */
180 
181  int filterpos; /**< position of eventdata in SCIP's event filter, -1 if not catching events */
182 };
183 
184 /** constraint data for nonlinear constraints */
185 struct SCIP_ConsData
186 {
187  /* data that defines the constraint: expression and sides */
188  SCIP_EXPR* expr; /**< expression that represents this constraint */
189  SCIP_Real lhs; /**< left-hand side */
190  SCIP_Real rhs; /**< right-hand side */
191 
192  /* variables */
193  SCIP_EXPR** varexprs; /**< array containing all variable expressions */
194  int nvarexprs; /**< total number of variable expressions */
195  SCIP_Bool catchedevents; /**< do we catch events on variables? */
196 
197  /* constraint violation */
198  SCIP_Real lhsviol; /**< violation of left-hand side by current solution */
199  SCIP_Real rhsviol; /**< violation of right-hand side by current solution */
200  SCIP_Real gradnorm; /**< norm of gradient of constraint function in current solution (if evaluated) */
201  SCIP_Longint gradnormsoltag; /**< tag of solution used that gradnorm corresponds to */
202 
203  /* status flags */
204  unsigned int ispropagated:1; /**< did we propagate the current bounds already? */
205  unsigned int issimplified:1; /**< did we simplify the expression tree already? */
206 
207  /* locks */
208  int nlockspos; /**< number of positive locks */
209  int nlocksneg; /**< number of negative locks */
210 
211  /* repair infeasible solutions */
212  SCIP_VAR* linvardecr; /**< variable that may be decreased without making any other constraint infeasible, or NULL if none */
213  SCIP_VAR* linvarincr; /**< variable that may be increased without making any other constraint infeasible, or NULL if none */
214  SCIP_Real linvardecrcoef; /**< linear coefficient of linvardecr */
215  SCIP_Real linvarincrcoef; /**< linear coefficient of linvarincr */
216 
217  /* miscellaneous */
218  SCIP_EXPRCURV curv; /**< curvature of the root expression w.r.t. the original variables */
219  SCIP_NLROW* nlrow; /**< a nonlinear row representation of this constraint */
220  int consindex; /**< an index of the constraint that is unique among all expr-constraints in this SCIP instance and is constant */
221 };
222 
223 /** constraint upgrade method */
224 typedef struct
225 {
226  SCIP_DECL_NONLINCONSUPGD((*consupgd)); /**< method to call for upgrading nonlinear constraint */
227  int priority; /**< priority of upgrading method */
228  SCIP_Bool active; /**< is upgrading enabled */
230 
231 /** constraint handler data */
232 struct SCIP_ConshdlrData
233 {
234  /* nonlinear handler */
235  SCIP_NLHDLR** nlhdlrs; /**< nonlinear handlers */
236  int nnlhdlrs; /**< number of nonlinear handlers */
237  int nlhdlrssize; /**< size of nlhdlrs array */
238  SCIP_Bool indetect; /**< whether we are currently in detectNlhdlr */
239  SCIP_Bool registerusesactivitysepabelow; /**< a flag that is used only during \ref @detectNlhdlr() */
240  SCIP_Bool registerusesactivitysepaabove; /**< a flag that is used only during \ref @detectNlhdlr() */
241 
242  /* constraint upgrades */
243  CONSUPGRADE** consupgrades; /**< constraint upgrade methods for specializing nonlinear constraints */
244  int consupgradessize; /**< size of consupgrades array */
245  int nconsupgrades; /**< number of constraint upgrade methods */
246 
247  /* other plugins */
248  SCIP_EVENTHDLR* eventhdlr; /**< handler for variable bound change events */
249  SCIP_HEUR* subnlpheur; /**< a pointer to the subnlp heuristic, if available */
250  SCIP_HEUR* trysolheur; /**< a pointer to the trysol heuristic, if available */
251 
252  /* tags and counters */
253  int auxvarid; /**< unique id for the next auxiliary variable */
254  SCIP_Longint curboundstag; /**< tag indicating current variable bounds */
255  SCIP_Longint lastboundrelax; /**< tag when bounds where most recently relaxed */
256  SCIP_Longint lastvaractivitymethodchange; /**< tag when method used to evaluate activity of variables changed last */
257  unsigned int enforound; /**< total number of enforcement calls, including current one */
258  int lastconsindex; /**< last used consindex, plus one */
259 
260  /* activity intervals and domain propagation */
261  SCIP_DECL_EXPR_INTEVALVAR((*intevalvar)); /**< method currently used for activity calculation of variable expressions */
262  SCIP_Bool globalbounds; /**< whether global variable bounds should be used for activity calculation */
263  SCIP_QUEUE* reversepropqueue; /**< expression queue to be used in reverse propagation, filled by SCIPtightenExprIntervalNonlinear */
264  SCIP_Bool forceboundtightening; /**< whether bound change passed to SCIPtightenExprIntervalNonlinear should be forced */
265  unsigned int curpropboundstag; /**< tag indicating current propagation rounds, to match with expr->propboundstag */
266 
267  /* parameters */
268  int maxproprounds; /**< limit on number of propagation rounds for a set of constraints within one round of SCIP propagation */
269  SCIP_Bool propauxvars; /**< whether to check bounds of all auxiliary variable to seed reverse propagation */
270  char varboundrelax; /**< strategy on how to relax variable bounds during bound tightening */
271  SCIP_Real varboundrelaxamount; /**< by how much to relax variable bounds during bound tightening */
272  SCIP_Real conssiderelaxamount; /**< by how much to relax constraint sides during bound tightening */
273  SCIP_Real vp_maxperturb; /**< maximal relative perturbation of reference point */
274  SCIP_Real vp_adjfacetthreshold; /**< adjust computed facet up to a violation of this value times lpfeastol */
275  SCIP_Bool vp_dualsimplex; /**< whether to use dual simplex instead of primal simplex for facet computing LP */
276  SCIP_Bool reformbinprods; /**< whether to reformulate products of binary variables during presolving */
277  SCIP_Bool reformbinprodsand; /**< whether to use the AND constraint handler for reformulating binary products */
278  int reformbinprodsfac; /**< minimum number of terms to reformulate bilinear binary products by factorizing variables (<= 1: disabled) */
279  SCIP_Bool forbidmultaggrnlvar; /**< whether to forbid multiaggregation of variables that appear in a nonlinear term of a constraint */
280  SCIP_Bool tightenlpfeastol; /**< whether to tighten LP feasibility tolerance during enforcement, if it seems useful */
281  SCIP_Bool propinenforce; /**< whether to (re)run propagation in enforcement */
282  SCIP_Real weakcutthreshold; /**< threshold for when to regard a cut from an estimator as weak */
283  SCIP_Real strongcutmaxcoef; /**< "strong" cuts will be scaled to have their maximal coef in [1/strongcutmaxcoef,strongcutmaxcoef] */
284  SCIP_Bool strongcutefficacy; /**< consider efficacy requirement when deciding whether a cut is "strong" */
285  SCIP_Bool forcestrongcut; /**< whether to force "strong" cuts in enforcement */
286  SCIP_Real enfoauxviolfactor; /**< an expression will be enforced if the "auxiliary" violation is at least enfoauxviolfactor times the "original" violation */
287  SCIP_Real weakcutminviolfactor; /**< retry with weak cuts for constraints with violation at least this factor of maximal violated constraints */
288  char rownotremovable; /**< whether to make rows to be non-removable in the node where they are added (can prevent some cycling): 'o'ff, in 'e'nforcement only, 'a'lways */
289  char violscale; /**< method how to scale violations to make them comparable (not used for feasibility check) */
290  char checkvarlocks; /**< whether variables contained in a single constraint should be forced to be at their lower or upper bounds ('d'isable, change 't'ype, add 'b'ound disjunction) */
291  int branchauxmindepth; /**< from which depth on to allow branching on auxiliary variables */
292  SCIP_Bool branchexternal; /**< whether to use external branching candidates for branching */
293  SCIP_Real branchhighviolfactor; /**< consider a constraint highly violated if its violation is >= this factor * maximal violation among all constraints */
294  SCIP_Real branchhighscorefactor; /**< consider a variable branching score high if its branching score >= this factor * maximal branching score among all variables */
295  SCIP_Real branchviolweight; /**< weight by how much to consider the violation assigned to a variable for its branching score */
296  SCIP_Real branchdualweight; /**< weight by how much to consider the dual values of rows that contain a variable for its branching score */
297  SCIP_Real branchpscostweight; /**< weight by how much to consider the pseudo cost of a variable for its branching score */
298  SCIP_Real branchdomainweight; /**< weight by how much to consider the domain width in branching score */
299  SCIP_Real branchvartypeweight;/**< weight by how much to consider variable type in branching score */
300  char branchscoreagg; /**< how to aggregate several branching scores given for the same expression ('a'verage, 'm'aximum, or 's'um) */
301  char branchviolsplit; /**< method used to split violation in expression onto variables ('u'niform, 'm'idness of solution, 'd'omain width, 'l'ogarithmic domain width) */
302  SCIP_Real branchpscostreliable; /**< minimum pseudo-cost update count required to consider pseudo-costs reliable */
303  char linearizeheursol; /**< whether tight linearizations of nonlinear constraints should be added to cutpool when some heuristics finds a new solution ('o'ff, on new 'i'ncumbents, on 'e'very solution) */
304  SCIP_Bool assumeconvex; /**< whether to assume that any constraint is convex */
305 
306  /* statistics */
307  SCIP_Longint nweaksepa; /**< number of times we used "weak" cuts for enforcement */
308  SCIP_Longint ntightenlp; /**< number of times we requested solving the LP with a smaller feasibility tolerance when enforcing */
309  SCIP_Longint ndesperatetightenlp; /**< number of times we requested solving the LP with a smaller feasibility tolerance when enforcing because we didn't know anything better */
310  SCIP_Longint ndesperatebranch; /**< number of times we branched on some variable because normal enforcement was not successful */
311  SCIP_Longint ndesperatecutoff; /**< number of times we cut off a node in enforcement because no branching candidate could be found */
312  SCIP_Longint nforcelp; /**< number of times we forced solving the LP when enforcing a pseudo solution */
313  SCIP_CLOCK* canonicalizetime; /**< time spend for canonicalization */
314  SCIP_Longint ncanonicalizecalls; /**< number of times we called canonicalization */
315 
316  /* facets of envelops of vertex-polyhedral functions */
317  SCIP_RANDNUMGEN* vp_randnumgen; /**< random number generator used to perturb reference point */
318  SCIP_LPI* vp_lp[SCIP_MAXVERTEXPOLYDIM+1]; /**< LPs used to compute facets for functions of different dimension */
319 
320  /* hashing of bilinear terms */
321  SCIP_HASHTABLE* bilinhashtable; /**< hash table for bilinear terms */
322  SCIP_CONSNONLINEAR_BILINTERM* bilinterms; /**< bilinear terms */
323  int nbilinterms; /**< total number of bilinear terms */
324  int bilintermssize; /**< size of bilinterms array */
325  int bilinmaxnauxexprs; /**< maximal number of auxiliary expressions per bilinear term */
326 
327  /* branching */
328  SCIP_RANDNUMGEN* branchrandnumgen; /**< random number generated used in branching variable selection */
329  char branchpscostupdatestrategy; /**< value of parameter branching/lpgainnormalize */
330 
331  /* misc */
332  SCIP_Bool checkedvarlocks; /**< whether variables contained in a single constraint have been already considered */
333  SCIP_HASHMAP* var2expr; /**< hashmap to map SCIP variables to variable-expressions */
334  int newsoleventfilterpos; /**< filter position of new solution event handler, if caught */
335 };
336 
337 /** branching candidate with various scores */
338 typedef struct
339 {
340  SCIP_EXPR* expr; /**< expression that holds branching candidate */
341  SCIP_Real auxviol; /**< aux-violation score of candidate */
342  SCIP_Real domain; /**< domain score of candidate */
343  SCIP_Real dual; /**< dual score of candidate */
344  SCIP_Real pscost; /**< pseudo-cost score of candidate */
345  SCIP_Real vartype; /**< variable type score of candidate */
346  SCIP_Real weighted; /**< weighted sum of other scores, see scoreBranchingCandidates() */
348 
349 /*
350  * Local methods
351  */
352 
353 /* forward declaration */
354 static
356  SCIP* scip, /**< SCIP data structure */
357  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
358  SCIP_EXPR* rootexpr, /**< expression */
359  SCIP_Bool tightenauxvars, /**< should the bounds of auxiliary variables be tightened? */
360  SCIP_Bool* infeasible, /**< buffer to store whether the problem is infeasible (NULL if not needed) */
361  int* ntightenings /**< buffer to store the number of auxiliary variable tightenings (NULL if not needed) */
362  );
363 
364 /** frees auxiliary variables of expression, if any */
365 static
367  SCIP* scip, /**< SCIP data structure */
368  SCIP_EXPR* expr /**< expression which auxvar to free, if any */
369  )
370 {
371  SCIP_EXPR_OWNERDATA* mydata;
372 
373  assert(scip != NULL);
374  assert(expr != NULL);
375 
376  mydata = SCIPexprGetOwnerData(expr);
377  assert(mydata != NULL);
378 
379  if( mydata->auxvar == NULL )
380  return SCIP_OKAY;
381 
382  SCIPdebugMsg(scip, "remove auxiliary variable <%s> for expression %p\n", SCIPvarGetName(mydata->auxvar), (void*)expr);
383 
384  /* remove variable locks
385  * as this is a relaxation-only variable, no other plugin should use it for deducing any type of reductions or cutting planes
386  */
387  SCIP_CALL( SCIPaddVarLocks(scip, mydata->auxvar, -1, -1) );
388 
389  /* release auxiliary variable */
390  SCIP_CALL( SCIPreleaseVar(scip, &mydata->auxvar) );
391  assert(mydata->auxvar == NULL);
392 
393  return SCIP_OKAY;
394 }
395 
396 /** frees data used for enforcement of expression, that is, nonlinear handlers
397  *
398  * can also clear indicators whether expr needs enforcement methods, that is,
399  * free an associated auxiliary variable and reset the nactivityuses counts
400  */
401 static
403  SCIP* scip, /**< SCIP data structure */
404  SCIP_EXPR* expr, /**< expression whose enforcement data will be released */
405  SCIP_Bool freeauxvar /**< whether aux var should be released and activity usage counts be reset */
406  )
407 {
408  SCIP_EXPR_OWNERDATA* mydata;
409  int e;
410 
411  mydata = SCIPexprGetOwnerData(expr);
412  assert(mydata != NULL);
413 
414  if( freeauxvar )
415  {
416  /* free auxiliary variable */
417  SCIP_CALL( freeAuxVar(scip, expr) );
418  assert(mydata->auxvar == NULL);
419 
420  /* reset count on activity and auxvar usage */
421  mydata->nactivityusesprop = 0;
422  mydata->nactivityusessepa = 0;
423  mydata->nauxvaruses = 0;
424  }
425 
426  /* free data stored by nonlinear handlers */
427  for( e = 0; e < mydata->nenfos; ++e )
428  {
429  SCIP_NLHDLR* nlhdlr;
430 
431  assert(mydata->enfos[e] != NULL);
432 
433  nlhdlr = mydata->enfos[e]->nlhdlr;
434  assert(nlhdlr != NULL);
435 
436  if( mydata->enfos[e]->issepainit )
437  {
438  /* call the separation deinitialization callback of the nonlinear handler */
439  SCIP_CALL( SCIPnlhdlrExitsepa(scip, nlhdlr, expr, mydata->enfos[e]->nlhdlrexprdata) );
440  mydata->enfos[e]->issepainit = FALSE;
441  }
442 
443  /* free nlhdlr exprdata, if there is any and there is a method to free this data */
444  if( mydata->enfos[e]->nlhdlrexprdata != NULL )
445  {
446  SCIP_CALL( SCIPnlhdlrFreeexprdata(scip, nlhdlr, expr, &mydata->enfos[e]->nlhdlrexprdata) );
447  assert(mydata->enfos[e]->nlhdlrexprdata == NULL);
448  }
449 
450  /* free enfo data */
451  SCIPfreeBlockMemory(scip, &mydata->enfos[e]);
452  }
453 
454  /* free array with enfo data */
455  SCIPfreeBlockMemoryArrayNull(scip, &mydata->enfos, mydata->nenfos);
456 
457  /* we need to look at this expression in detect again */
458  mydata->nenfos = -1;
459 
460  return SCIP_OKAY;
461 }
462 
463 /** callback that frees data that this conshdlr stored in an expression */
464 static
465 SCIP_DECL_EXPR_OWNERFREE(exprownerFree)
466 {
467  assert(scip != NULL);
468  assert(expr != NULL);
469  assert(ownerdata != NULL);
470  assert(*ownerdata != NULL);
471 
472  /* expression should not be locked anymore */
473  assert((*ownerdata)->nlockspos == 0);
474  assert((*ownerdata)->nlocksneg == 0);
475 
476  SCIP_CALL( freeEnfoData(scip, expr, TRUE) );
477 
478  /* expression should not be enforced anymore */
479  assert((*ownerdata)->nenfos <= 0);
480  assert((*ownerdata)->auxvar == NULL);
481 
482  if( SCIPisExprVar(scip, expr) )
483  {
484  SCIP_CONSHDLRDATA* conshdlrdata;
485  SCIP_VAR* var;
486 
487  /* there should be no constraints left that still use this variable */
488  assert((*ownerdata)->nconss == 0);
489  /* thus, there should also be no variable event catched (via this exprhdlr) */
490  assert((*ownerdata)->filterpos == -1);
491 
492  SCIPfreeBlockMemoryArrayNull(scip, &(*ownerdata)->conss, (*ownerdata)->consssize);
493 
494  /* update var2expr hashmap in conshdlrdata */
495  conshdlrdata = SCIPconshdlrGetData((*ownerdata)->conshdlr);
496  assert(conshdlrdata != NULL);
497 
498  var = SCIPgetVarExprVar(expr);
499  assert(var != NULL);
500 
501  /* remove var -> expr map from hashmap if present
502  * (if no variable-expression stored for var hashmap, then the var hasn't been used in any constraint, so do nothing
503  * if variable-expression stored for var is different, then also do nothing)
504  */
505  if( SCIPhashmapGetImage(conshdlrdata->var2expr, var) == (void*)expr )
506  {
507  SCIP_CALL( SCIPhashmapRemove(conshdlrdata->var2expr, var) );
508  }
509  }
510 
511  SCIPfreeBlockMemory(scip, ownerdata);
512 
513  return SCIP_OKAY;
514 }
515 
516 static
517 SCIP_DECL_EXPR_OWNERPRINT(exprownerPrint)
518 { /*lint --e{715}*/
519  assert(ownerdata != NULL);
520 
521  /* print nl handlers associated to expr */
522  if( ownerdata->nenfos > 0 )
523  {
524  int i;
525  SCIPinfoMessage(scip, file, " {");
526 
527  for( i = 0; i < ownerdata->nenfos; ++i )
528  {
529  SCIPinfoMessage(scip, file, "%s:", SCIPnlhdlrGetName(ownerdata->enfos[i]->nlhdlr));
530  if( ownerdata->enfos[i]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY )
531  SCIPinfoMessage(scip, file, "a");
532  if( ownerdata->enfos[i]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPABELOW )
533  SCIPinfoMessage(scip, file, "u");
534  if( ownerdata->enfos[i]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPAABOVE )
535  SCIPinfoMessage(scip, file, "o");
536  if( i < ownerdata->nenfos-1 )
537  SCIPinfoMessage(scip, file, ", ");
538  }
539 
540  SCIPinfoMessage(scip, file, "}");
541  }
542 
543  /* print aux var associated to expr */
544  if( ownerdata->auxvar != NULL )
545  {
546  SCIPinfoMessage(scip, file, " (<%s> in [%g, %g])", SCIPvarGetName(ownerdata->auxvar), SCIPvarGetLbLocal(ownerdata->auxvar), SCIPvarGetUbLocal(ownerdata->auxvar));
547  }
548  SCIPinfoMessage(scip, file, "\n");
549 
550  return SCIP_OKAY;
551 }
552 
553 /** possibly reevaluates and then returns the activity of the expression
554  *
555  * Reevaluate activity if currently stored is not up to date (some bound was changed since last evaluation).
556  */
557 static
558 SCIP_DECL_EXPR_OWNEREVALACTIVITY(exprownerEvalactivity)
559 {
560  SCIP_CONSHDLRDATA* conshdlrdata;
561 
562  assert(scip != NULL);
563  assert(expr != NULL);
564  assert(ownerdata != NULL);
565 
566  conshdlrdata = SCIPconshdlrGetData(ownerdata->conshdlr);
567  assert(conshdlrdata != NULL);
568 
569  if( SCIPexprGetActivityTag(expr) < conshdlrdata->curboundstag )
570  {
571  /* update activity of expression */
572  SCIP_CALL( forwardPropExpr(scip, ownerdata->conshdlr, expr, FALSE, NULL, NULL) );
573 
574  assert(SCIPexprGetActivityTag(expr) == conshdlrdata->curboundstag);
575  }
576 
577  return SCIP_OKAY;
578 }
579 
580 /** callback that creates data that this conshdlr wants to store in an expression */
581 static
582 SCIP_DECL_EXPR_OWNERCREATE(exprownerCreate)
583 {
584  assert(scip != NULL);
585  assert(expr != NULL);
586  assert(ownerdata != NULL);
587 
588  SCIP_CALL( SCIPallocClearBlockMemory(scip, ownerdata) );
589  (*ownerdata)->nenfos = -1;
590  (*ownerdata)->conshdlr = (SCIP_CONSHDLR*)ownercreatedata;
591 
592  if( SCIPisExprVar(scip, expr) )
593  {
594  SCIP_CONSHDLRDATA* conshdlrdata;
595  SCIP_VAR* var;
596 
597  (*ownerdata)->filterpos = -1;
598 
599  /* add to var2expr hashmap if not having expr for var yet */
600 
601  conshdlrdata = SCIPconshdlrGetData((*ownerdata)->conshdlr);
602  assert(conshdlrdata != NULL);
603 
604  var = SCIPgetVarExprVar(expr);
605 
606  if( !SCIPhashmapExists(conshdlrdata->var2expr, (void*)var) )
607  {
608  /* store the variable expression in the hashmap */
609  SCIP_CALL( SCIPhashmapInsert(conshdlrdata->var2expr, (void*)var, (void*)expr) );
610  }
611  else
612  {
613  /* if expr was just created, then it shouldn't already be stored as image of var */
614  assert(SCIPhashmapGetImage(conshdlrdata->var2expr, (void*)var) != (void*)expr);
615  }
616  }
617  else
618  {
619  /* just so that we can use filterpos to recognize whether an expr is a varexpr if not having a SCIP pointer around */
620  (*ownerdata)->filterpos = -2;
621  }
622 
623  *ownerfree = exprownerFree;
624  *ownerprint = exprownerPrint;
625  *ownerevalactivity = exprownerEvalactivity;
626 
627  return SCIP_OKAY;
628 }
629 
630 /** creates a variable expression or retrieves from hashmap in conshdlr data */
631 static
633  SCIP* scip, /**< SCIP data structure */
634  SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */
635  SCIP_EXPR** expr, /**< pointer where to store expression */
636  SCIP_VAR* var /**< variable to be stored */
637  )
638 {
639  assert(conshdlr != NULL);
640  assert(expr != NULL);
641  assert(var != NULL);
642 
643  /* get variable expression representing the given variable if there is one already */
644  *expr = (SCIP_EXPR*) SCIPhashmapGetImage(SCIPconshdlrGetData(conshdlr)->var2expr, (void*) var);
645 
646  if( *expr == NULL )
647  {
648  /* create a new variable expression; this also captures the expression */
649  SCIP_CALL( SCIPcreateExprVar(scip, expr, var, exprownerCreate, (void*)conshdlr) );
650  assert(*expr != NULL);
651  /* exprownerCreate should have added var->expr to var2expr */
652  assert(SCIPhashmapGetImage(SCIPconshdlrGetData(conshdlr)->var2expr, (void*)var) == (void*)*expr);
653  }
654  else
655  {
656  /* only capture already existing expr to get a consistent uses-count */
657  SCIPcaptureExpr(*expr);
658  }
659 
660  return SCIP_OKAY;
661 }
662 
663 /* map var exprs to var-expr from var2expr hashmap */
664 static
665 SCIP_DECL_EXPR_MAPEXPR(mapexprvar)
666 { /*lint --e{715}*/
667  SCIP_CONSHDLR* conshdlr = (SCIP_CONSHDLR*)mapexprdata;
668 
669  assert(sourcescip != NULL);
670  assert(targetscip != NULL);
671  assert(sourceexpr != NULL);
672  assert(targetexpr != NULL);
673  assert(*targetexpr == NULL);
674  assert(mapexprdata != NULL);
675 
676  /* do not provide map if not variable */
677  if( !SCIPisExprVar(sourcescip, sourceexpr) )
678  return SCIP_OKAY;
679 
680  SCIP_CALL( createExprVar(targetscip, conshdlr, targetexpr, SCIPgetVarExprVar(sourceexpr)) );
681 
682  return SCIP_OKAY;
683 }
684 
685 /* map var exprs to var-expr from var2expr hashmap corresponding to transformed var */
686 static
687 SCIP_DECL_EXPR_MAPEXPR(mapexprtransvar)
688 { /*lint --e{715}*/
689  SCIP_CONSHDLR* conshdlr = (SCIP_CONSHDLR*)mapexprdata;
690  SCIP_VAR* var;
691 
692  assert(sourcescip != NULL);
693  assert(targetscip != NULL);
694  assert(sourceexpr != NULL);
695  assert(targetexpr != NULL);
696  assert(*targetexpr == NULL);
697  assert(mapexprdata != NULL);
698 
699  /* do not provide map if not variable */
700  if( !SCIPisExprVar(sourcescip, sourceexpr) )
701  return SCIP_OKAY;
702 
703  var = SCIPgetVarExprVar(sourceexpr);
704  assert(var != NULL);
705 
706  /* transform variable */
707  SCIP_CALL( SCIPgetTransformedVar(sourcescip, var, &var) );
708  assert(var != NULL);
709 
710  SCIP_CALL( createExprVar(targetscip, conshdlr, targetexpr, var) );
711 
712  return SCIP_OKAY;
713 }
714 
715 /** stores all variable expressions into a given constraint */
716 static
718  SCIP* scip, /**< SCIP data structure */
719  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
720  SCIP_CONSDATA* consdata /**< constraint data */
721  )
722 {
723  SCIP_CONSHDLRDATA* conshdlrdata;
724  int i;
725 
726  assert(consdata != NULL);
727 
728  /* skip if we have stored the variable expressions already */
729  if( consdata->varexprs != NULL )
730  return SCIP_OKAY;
731 
732  assert(consdata->varexprs == NULL);
733  assert(consdata->nvarexprs == 0);
734 
735  /* create array to store all variable expressions; the number of variable expressions is bounded by SCIPgetNTotalVars() */
736  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->varexprs, SCIPgetNTotalVars(scip)) );
737 
738  SCIP_CALL( SCIPgetExprVarExprs(scip, consdata->expr, consdata->varexprs, &(consdata->nvarexprs)) );
739  assert(SCIPgetNTotalVars(scip) >= consdata->nvarexprs);
740 
741  /* shrink array if there are less variables in the expression than in the problem */
742  if( SCIPgetNTotalVars(scip) > consdata->nvarexprs )
743  {
744  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->varexprs, SCIPgetNTotalVars(scip), consdata->nvarexprs) );
745  }
746 
747  conshdlrdata = SCIPconshdlrGetData(conshdlr);
748  assert(conshdlrdata != NULL);
749  assert(conshdlrdata->var2expr != NULL);
750 
751  /* ensure that for every variable an entry exists in the var2expr hashmap
752  * when removing duplicate subexpressions it can happen than a var->varexpr map was removed from the hashmap
753  */
754  for( i = 0; i < consdata->nvarexprs; ++i )
755  {
756  if( !SCIPhashmapExists(conshdlrdata->var2expr, SCIPgetVarExprVar(consdata->varexprs[i])) )
757  {
758  SCIP_CALL( SCIPhashmapInsert(conshdlrdata->var2expr, SCIPgetVarExprVar(consdata->varexprs[i]), consdata->varexprs[i]) );
759  }
760  }
761 
762  return SCIP_OKAY;
763 }
764 
765 /** frees all variable expression stored in storeVarExprs() */
766 static
768  SCIP* scip, /**< SCIP data structure */
769  SCIP_CONSDATA* consdata /**< constraint data */
770  )
771 {
772  int i;
773 
774  assert(consdata != NULL);
775 
776  /* skip if we have stored the variable expressions already*/
777  if( consdata->varexprs == NULL )
778  return SCIP_OKAY;
779 
780  assert(consdata->varexprs != NULL);
781  assert(consdata->nvarexprs >= 0);
782 
783  /* release variable expressions */
784  for( i = 0; i < consdata->nvarexprs; ++i )
785  {
786  assert(consdata->varexprs[i] != NULL);
787  SCIP_CALL( SCIPreleaseExpr(scip, &consdata->varexprs[i]) );
788  assert(consdata->varexprs[i] == NULL);
789  }
790 
791  /* free variable expressions */
792  SCIPfreeBlockMemoryArrayNull(scip, &consdata->varexprs, consdata->nvarexprs);
793  consdata->varexprs = NULL;
794  consdata->nvarexprs = 0;
795 
796  return SCIP_OKAY;
797 }
798 
799 /** interval evaluation of variables as used in bound tightening
800  *
801  * Returns slightly relaxed local variable bounds of a variable as interval.
802  * Does not relax beyond integer values, thus does not relax bounds on integer variables at all.
803  */
804 static
805 SCIP_DECL_EXPR_INTEVALVAR(intEvalVarBoundTightening)
806 {
807  SCIP_INTERVAL interval;
808  SCIP_CONSHDLRDATA* conshdlrdata;
809  SCIP_Real lb;
810  SCIP_Real ub;
811 
812  assert(scip != NULL);
813  assert(var != NULL);
814 
815  conshdlrdata = (SCIP_CONSHDLRDATA*)intevalvardata;
816  assert(conshdlrdata != NULL);
817 
818  if( conshdlrdata->globalbounds )
819  {
820  lb = SCIPvarGetLbGlobal(var);
821  ub = SCIPvarGetUbGlobal(var);
822  }
823  else
824  {
825  lb = SCIPvarGetLbLocal(var);
826  ub = SCIPvarGetUbLocal(var);
827  }
828  assert(lb <= ub); /* SCIP should ensure that variable bounds are not contradicting */
829 
830  /* implicit integer variables may have non-integer bounds, apparently (run space25a) */
832  {
833  lb = EPSROUND(lb, 0.0); /*lint !e835*/
834  ub = EPSROUND(ub, 0.0); /*lint !e835*/
835  }
836 
837  /* integer variables should always have integral bounds in SCIP */
838  assert(EPSFRAC(lb, 0.0) == 0.0 || !SCIPvarIsIntegral(var)); /*lint !e835*/
839  assert(EPSFRAC(ub, 0.0) == 0.0 || !SCIPvarIsIntegral(var)); /*lint !e835*/
840 
841  switch( conshdlrdata->varboundrelax )
842  {
843  case 'n' : /* no relaxation */
844  break;
845 
846  case 'a' : /* relax by absolute value */
847  {
848  /* do not look at integer variables, they already have integral bounds, so wouldn't be relaxed */
849  if( SCIPvarIsIntegral(var) )
850  break;
851 
852  if( !SCIPisInfinity(scip, -lb) )
853  {
854  /* reduce lb by epsilon, or to the next integer value, which ever is larger */
855  SCIP_Real bnd = floor(lb);
856  lb = MAX(bnd, lb - conshdlrdata->varboundrelaxamount);
857  }
858 
859  if( !SCIPisInfinity(scip, ub) )
860  {
861  /* increase ub by epsilon, or to the next integer value, which ever is smaller */
862  SCIP_Real bnd = ceil(ub);
863  ub = MIN(bnd, ub + conshdlrdata->varboundrelaxamount);
864  }
865 
866  break;
867  }
868 
869  case 'b' : /* relax always by absolute value */
870  {
871  /* do not look at integer variables, they already have integral bounds, so wouldn't be relaxed */
872  if( SCIPvarIsIntegral(var) )
873  break;
874 
875  if( !SCIPisInfinity(scip, -lb) )
876  lb -= conshdlrdata->varboundrelaxamount;
877 
878  if( !SCIPisInfinity(scip, ub) )
879  ub += conshdlrdata->varboundrelaxamount;
880 
881  break;
882  }
883 
884  case 'r' : /* relax by relative value */
885  {
886  /* do not look at integer variables, they already have integral bounds, so wouldn't be relaxed */
887  if( SCIPvarIsIntegral(var) )
888  break;
889 
890  /* relax bounds by epsilon*max(1,|bnd|), instead of just epsilon as in case 'a', thus we trust the first log(epsilon) digits
891  * however, when domains get small, relaxing can excessively weaken bound tightening, thus do only fraction of |ub-lb| if that is smaller
892  * further, do not relax beyond next integer value
893  */
894  if( !SCIPisInfinity(scip, -lb) )
895  {
896  SCIP_Real bnd = floor(lb);
897  lb = MAX(bnd, lb - MIN(conshdlrdata->varboundrelaxamount * MAX(1.0, REALABS(lb)), 0.001 * REALABS(ub-lb)));
898  }
899 
900  if( !SCIPisInfinity(scip, ub) )
901  {
902  SCIP_Real bnd = ceil(ub);
903  ub = MIN(bnd, ub + MIN(conshdlrdata->varboundrelaxamount * MAX(1.0, REALABS(ub)), 0.001 * REALABS(ub-lb)));
904  }
905 
906  break;
907  }
908 
909  default :
910  {
911  SCIPerrorMessage("Unsupported value '%c' for varboundrelax option.\n", conshdlrdata->varboundrelax);
912  SCIPABORT();
913  break;
914  }
915  }
916 
917  /* convert SCIPinfinity() to SCIP_INTERVAL_INFINITY */
920  assert(lb <= ub);
921 
922  SCIPintervalSetBounds(&interval, lb, ub);
923 
924  return interval;
925 }
926 
927 /** compares two nonlinear constraints by its index
928  *
929  * Usable as compare operator in array sort functions.
930  */
931 static
932 SCIP_DECL_SORTPTRCOMP(compIndexConsNonlinear)
933 {
934  SCIP_CONSDATA* consdata1 = SCIPconsGetData((SCIP_CONS*)elem1);
935  SCIP_CONSDATA* consdata2 = SCIPconsGetData((SCIP_CONS*)elem2);
936 
937  assert(consdata1 != NULL);
938  assert(consdata2 != NULL);
939 
940  return consdata1->consindex - consdata2->consindex;
941 }
942 
943 /** processes variable fixing or bound change event */
944 static
945 SCIP_DECL_EVENTEXEC(processVarEvent)
946 { /*lint --e{715}*/
947  SCIP_EVENTTYPE eventtype;
948  SCIP_EXPR* expr;
949  SCIP_EXPR_OWNERDATA* ownerdata;
950 
951  eventtype = SCIPeventGetType(event);
953 
954  assert(eventdata != NULL);
955  expr = (SCIP_EXPR*) eventdata;
956  assert(SCIPisExprVar(scip, expr));
957 
958  SCIPdebugMsg(scip, " exec event %" SCIP_EVENTTYPE_FORMAT " for variable <%s> (local [%g,%g], global [%g,%g])\n", eventtype,
962 
963  ownerdata = SCIPexprGetOwnerData(expr);
964  assert(ownerdata != NULL);
965  /* we only catch varevents for variables in constraints, so there should be constraints */
966  assert(ownerdata->nconss > 0);
967  assert(ownerdata->conss != NULL);
968 
969  /* notify constraints that use this variable expression (expr) to repropagate and possibly resimplify
970  * - propagation can only find something new if a bound was tightened
971  * - simplify can only find something new if a var is fixed (or maybe a bound is tightened)
972  * and we look at global changes (that is, we are not looking at boundchanges in probing)
973  */
975  {
976  SCIP_CONSDATA* consdata;
977  int c;
978 
979  for( c = 0; c < ownerdata->nconss; ++c )
980  {
981  assert(ownerdata->conss[c] != NULL);
982  consdata = SCIPconsGetData(ownerdata->conss[c]);
983 
984  /* if bound tightening, then mark constraints to be propagated again
985  * TODO we could try be more selective here and only trigger a propagation if a relevant bound has changed,
986  * that is, we don't need to repropagate x + ... <= rhs if only the upper bound of x has been tightened
987  * the locks don't help since they are not available separately for each constraint
988  */
989  if( eventtype & SCIP_EVENTTYPE_BOUNDTIGHTENED )
990  {
991  consdata->ispropagated = FALSE;
992  SCIPdebugMsg(scip, " marked <%s> for propagate\n", SCIPconsGetName(ownerdata->conss[c]));
993  }
994 
995  /* if still in presolve (but not probing), then mark constraints to be unsimplified */
996  if( SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING && !SCIPinProbing(scip) )
997  {
998  consdata->issimplified = FALSE;
999  SCIPdebugMsg(scip, " marked <%s> for simplify\n", SCIPconsGetName(ownerdata->conss[c]));
1000  }
1001  }
1002  }
1003 
1004  /* update curboundstag, lastboundrelax, and expr activity */
1005  if( eventtype & SCIP_EVENTTYPE_BOUNDCHANGED )
1006  {
1007  SCIP_CONSHDLRDATA* conshdlrdata;
1008  SCIP_INTERVAL activity;
1009 
1010  conshdlrdata = SCIPconshdlrGetData(ownerdata->conshdlr);
1011  assert(conshdlrdata != NULL);
1012 
1013  /* increase tag on bounds */
1014  ++conshdlrdata->curboundstag;
1015  assert(conshdlrdata->curboundstag > 0);
1016 
1017  /* remember also if we relaxed bounds now */
1018  if( eventtype & SCIP_EVENTTYPE_BOUNDRELAXED )
1019  conshdlrdata->lastboundrelax = conshdlrdata->curboundstag;
1020 
1021  /* update the activity of the var-expr here immediately
1022  * (we could call expr->activity = intevalvar(var, consdhlr) directly, but then the exprhdlr statistics are not updated)
1023  */
1024  SCIP_CALL( SCIPcallExprInteval(scip, expr, &activity, conshdlrdata->intevalvar, conshdlrdata) );
1025  /* activity = conshdlrdata->intevalvar(scip, SCIPgetVarExprVar(expr), conshdlrdata); */
1026 #ifdef DEBUG_PROP
1027  SCIPdebugMsg(scip, " var-exprhdlr::inteval = [%.20g, %.20g]\n", activity.inf, activity.sup);
1028 #endif
1029  SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);
1030  }
1031 
1032  return SCIP_OKAY;
1033 }
1034 
1035 /** registers event handler to catch variable events on variable
1036  *
1037  * Additionally, the given constraint is stored in the ownerdata of the variable-expression.
1038  * When an event occurs, all stored constraints are notified.
1039  */
1040 static
1042  SCIP* scip, /**< SCIP data structure */
1043  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
1044  SCIP_EXPR* expr, /**< variable expression */
1045  SCIP_CONS* cons /**< nonlinear constraint */
1046  )
1047 {
1048  SCIP_EXPR_OWNERDATA* ownerdata;
1049 
1050  assert(eventhdlr != NULL);
1051  assert(expr != NULL);
1052  assert(SCIPisExprVar(scip, expr));
1053  assert(cons != NULL);
1054 
1055  ownerdata = SCIPexprGetOwnerData(expr);
1056  assert(ownerdata != NULL);
1057 
1058 #ifndef NDEBUG
1059  /* assert that constraint does not double-catch variable */
1060  {
1061  int i;
1062  for( i = 0; i < ownerdata->nconss; ++i )
1063  assert(ownerdata->conss[i] != cons);
1064  }
1065 #endif
1066 
1067  /* append cons to ownerdata->conss */
1068  SCIP_CALL( SCIPensureBlockMemoryArray(scip, &ownerdata->conss, &ownerdata->consssize, ownerdata->nconss + 1) );
1069  ownerdata->conss[ownerdata->nconss++] = cons;
1070  /* we're not capturing the constraint here to avoid circular references */
1071 
1072  /* updated sorted flag */
1073  if( ownerdata->nconss <= 1 )
1074  ownerdata->consssorted = TRUE;
1075  else if( ownerdata->consssorted )
1076  ownerdata->consssorted = compIndexConsNonlinear(ownerdata->conss[ownerdata->nconss-2], ownerdata->conss[ownerdata->nconss-1]) > 0;
1077 
1078  /* catch variable events, if not done so yet (first constraint) */
1079  if( ownerdata->filterpos < 0 )
1080  {
1081  SCIP_EVENTTYPE eventtype;
1082 
1083  assert(ownerdata->nconss == 1);
1084 
1086 
1087  SCIP_CALL( SCIPcatchVarEvent(scip, SCIPgetVarExprVar(expr), eventtype, eventhdlr, (SCIP_EVENTDATA*)expr, &ownerdata->filterpos) );
1088  assert(ownerdata->filterpos >= 0);
1089  }
1090 
1091  return SCIP_OKAY;
1092 }
1093 
1094 /** catch variable events */
1095 static
1097  SCIP* scip, /**< SCIP data structure */
1098  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
1099  SCIP_CONS* cons /**< constraint for which to catch bound change events */
1100  )
1101 {
1102  SCIP_CONSHDLRDATA* conshdlrdata;
1103  SCIP_CONSDATA* consdata;
1104  SCIP_EXPR* expr;
1105  int i;
1106 
1107  assert(eventhdlr != NULL);
1108  assert(cons != NULL);
1109 
1110  consdata = SCIPconsGetData(cons);
1111  assert(consdata != NULL);
1112  assert(consdata->varexprs != NULL);
1113  assert(consdata->nvarexprs >= 0);
1114 
1115  /* check if we have catched variable events already */
1116  if( consdata->catchedevents )
1117  return SCIP_OKAY;
1118 
1119  conshdlrdata = SCIPconshdlrGetData(SCIPconsGetHdlr(cons));
1120  assert(conshdlrdata != NULL);
1121  assert(conshdlrdata->intevalvar == intEvalVarBoundTightening);
1122 
1123  SCIPdebugMsg(scip, "catchVarEvents for %s\n", SCIPconsGetName(cons));
1124 
1125  for( i = 0; i < consdata->nvarexprs; ++i )
1126  {
1127  expr = consdata->varexprs[i];
1128 
1129  assert(expr != NULL);
1130  assert(SCIPisExprVar(scip, expr));
1131 
1132  SCIP_CALL( catchVarEvent(scip, eventhdlr, expr, cons) );
1133 
1134  /* from now on, activity of var-expr will usually be updated in processVarEvent if variable bound is changing
1135  * since we just registered this eventhdlr, we should make sure that the activity is also up to date now
1136  */
1137  if( SCIPexprGetActivityTag(expr) < conshdlrdata->curboundstag )
1138  {
1139  SCIP_INTERVAL activity;
1140  SCIP_CALL( SCIPcallExprInteval(scip, expr, &activity, intEvalVarBoundTightening, conshdlrdata) );
1141  /* activity = intEvalVarBoundTightening(scip, SCIPgetVarExprVar(expr), conshdlrdata); */
1142  SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);
1143 #ifdef DEBUG_PROP
1144  SCIPdebugMsg(scip, "var-exprhdlr::inteval for var <%s> = [%.20g, %.20g]\n", SCIPvarGetName(SCIPgetVarExprVar(expr)), activity.inf, activity.sup);
1145 #endif
1146  }
1147  }
1148 
1149  consdata->catchedevents = TRUE;
1150 
1151  return SCIP_OKAY;
1152 }
1153 
1154 /** unregisters event handler to catch variable events on variable
1155  *
1156  * The given constraint is removed from the constraints array in the ownerdata of the variable-expression.
1157  * If this was the last constraint, then the event handler is unregistered for this variable.
1158  */
1159 static
1161  SCIP* scip, /**< SCIP data structure */
1162  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
1163  SCIP_EXPR* expr, /**< variable expression */
1164  SCIP_CONS* cons /**< expr constraint */
1165  )
1166 {
1167  SCIP_EXPR_OWNERDATA* ownerdata;
1168  int pos;
1169 
1170  assert(eventhdlr != NULL);
1171  assert(expr != NULL);
1172  assert(SCIPisExprVar(scip, expr));
1173  assert(cons != NULL);
1174 
1175  ownerdata = SCIPexprGetOwnerData(expr);
1176  assert(ownerdata != NULL);
1177  assert(ownerdata->nconss > 0);
1178 
1179  if( ownerdata->conss[ownerdata->nconss-1] == cons )
1180  {
1181  pos = ownerdata->nconss-1;
1182  }
1183  else
1184  {
1185  if( !ownerdata->consssorted )
1186  {
1187  SCIPsortPtr((void**)ownerdata->conss, compIndexConsNonlinear, ownerdata->nconss);
1188  ownerdata->consssorted = TRUE;
1189  }
1190 
1191  if( !SCIPsortedvecFindPtr((void**)ownerdata->conss, compIndexConsNonlinear, cons, ownerdata->nconss, &pos) )
1192  {
1193  SCIPerrorMessage("Constraint <%s> not in constraint array of expression for variable <%s>\n", SCIPconsGetName(cons), SCIPvarGetName(SCIPgetVarExprVar(expr)));
1194  return SCIP_ERROR;
1195  }
1196  assert(pos >= 0 && pos < ownerdata->nconss);
1197  }
1198  assert(ownerdata->conss[pos] == cons);
1199 
1200  /* move last constraint into position of removed constraint */
1201  if( pos < ownerdata->nconss-1 )
1202  {
1203  ownerdata->conss[pos] = ownerdata->conss[ownerdata->nconss-1];
1204  ownerdata->consssorted = FALSE;
1205  }
1206  --ownerdata->nconss;
1207 
1208  /* drop variable events if that was the last constraint */
1209  if( ownerdata->nconss == 0 )
1210  {
1211  SCIP_EVENTTYPE eventtype;
1212 
1213  assert(ownerdata->filterpos >= 0);
1214 
1216 
1217  SCIP_CALL( SCIPdropVarEvent(scip, SCIPgetVarExprVar(expr), eventtype, eventhdlr, (SCIP_EVENTDATA*)expr, ownerdata->filterpos) );
1218  ownerdata->filterpos = -1;
1219  }
1220 
1221  return SCIP_OKAY;
1222 }
1223 
1224 /** drop variable events */
1225 static
1227  SCIP* scip, /**< SCIP data structure */
1228  SCIP_EVENTHDLR* eventhdlr, /**< event handler */
1229  SCIP_CONS* cons /**< constraint for which to drop bound change events */
1230  )
1231 {
1232  SCIP_CONSDATA* consdata;
1233  int i;
1234 
1235  assert(eventhdlr != NULL);
1236  assert(cons != NULL);
1237 
1238  consdata = SCIPconsGetData(cons);
1239  assert(consdata != NULL);
1240 
1241  /* check if we have catched variable events already */
1242  if( !consdata->catchedevents )
1243  return SCIP_OKAY;
1244 
1245  assert(consdata->varexprs != NULL);
1246  assert(consdata->nvarexprs >= 0);
1247 
1248  SCIPdebugMsg(scip, "dropVarEvents for %s\n", SCIPconsGetName(cons));
1249 
1250  for( i = consdata->nvarexprs - 1; i >= 0; --i )
1251  {
1252  assert(consdata->varexprs[i] != NULL);
1253 
1254  SCIP_CALL( dropVarEvent(scip, eventhdlr, consdata->varexprs[i], cons) );
1255  }
1256 
1257  consdata->catchedevents = FALSE;
1258 
1259  return SCIP_OKAY;
1260 }
1261 
1262 /** creates and captures a nonlinear constraint
1263  *
1264  * @attention Use copyexpr=FALSE only if expr is already "owned" by conshdlr, that is, if expressions were created with exprownerCreate() and ownerdata passed in the last two arguments
1265  */
1266 static
1268  SCIP* scip, /**< SCIP data structure */
1269  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1270  SCIP_CONS** cons, /**< pointer to hold the created constraint */
1271  const char* name, /**< name of constraint */
1272  SCIP_EXPR* expr, /**< expression of constraint (must not be NULL) */
1273  SCIP_Real lhs, /**< left hand side of constraint */
1274  SCIP_Real rhs, /**< right hand side of constraint */
1275  SCIP_Bool copyexpr, /**< whether to copy the expression or reuse the given expr (capture it) */
1276  SCIP_Bool initial, /**< should the LP relaxation of constraint be in the initial LP?
1277  * Usually set to TRUE. Set to FALSE for 'lazy constraints'. */
1278  SCIP_Bool separate, /**< should the constraint be separated during LP processing?
1279  * Usually set to TRUE. */
1280  SCIP_Bool enforce, /**< should the constraint be enforced during node processing?
1281  * TRUE for model constraints, FALSE for additional, redundant constraints. */
1282  SCIP_Bool check, /**< should the constraint be checked for feasibility?
1283  * TRUE for model constraints, FALSE for additional, redundant constraints. */
1284  SCIP_Bool propagate, /**< should the constraint be propagated during node processing?
1285  * Usually set to TRUE. */
1286  SCIP_Bool local, /**< is constraint only valid locally?
1287  * Usually set to FALSE. Has to be set to TRUE, e.g., for branching constraints. */
1288  SCIP_Bool modifiable, /**< is constraint modifiable (subject to column generation)?
1289  * Usually set to FALSE. In column generation applications, set to TRUE if pricing
1290  * adds coefficients to this constraint. */
1291  SCIP_Bool dynamic, /**< is constraint subject to aging?
1292  * Usually set to FALSE. Set to TRUE for own cuts which
1293  * are separated as constraints. */
1294  SCIP_Bool removable /**< should the relaxation be removed from the LP due to aging or cleanup?
1295  * Usually set to FALSE. Set to TRUE for 'lazy constraints' and 'user cuts'. */
1296  )
1297 {
1298  SCIP_CONSHDLRDATA* conshdlrdata;
1299  SCIP_CONSDATA* consdata;
1300 
1301  assert(conshdlr != NULL);
1302  assert(expr != NULL);
1303 
1304  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1305  assert(conshdlrdata != NULL);
1306 
1307  if( local && SCIPgetDepth(scip) != 0 )
1308  {
1309  SCIPerrorMessage("Locally valid nonlinear constraints are not supported, yet.\n");
1310  return SCIP_INVALIDCALL;
1311  }
1312 
1313  /* TODO we should allow for non-initial nonlinear constraints */
1314  if( !initial )
1315  {
1316  SCIPerrorMessage("Non-initial nonlinear constraints are not supported, yet.\n");
1317  return SCIP_INVALIDCALL;
1318  }
1319 
1320  /* create constraint data */
1321  SCIP_CALL( SCIPallocClearBlockMemory(scip, &consdata) );
1322 
1323  if( copyexpr )
1324  {
1325  /* copy expression, thereby map variables expressions to already existing variables expressions in var2expr map, or augment var2expr map */
1326  SCIP_CALL( SCIPduplicateExpr(scip, expr, &consdata->expr, mapexprvar, conshdlr, exprownerCreate, (void*)conshdlr) );
1327  }
1328  else
1329  {
1330  consdata->expr = expr;
1331  SCIPcaptureExpr(consdata->expr);
1332  }
1333  consdata->lhs = lhs;
1334  consdata->rhs = rhs;
1335  consdata->consindex = conshdlrdata->lastconsindex++;
1336  consdata->curv = SCIP_EXPRCURV_UNKNOWN;
1337 
1338  /* create constraint */
1339  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate,
1340  local, modifiable, dynamic, removable, FALSE) );
1341 
1342  return SCIP_OKAY;
1343 }
1344 
1345 /** returns absolute violation for auxvar relation in an expression w.r.t. original variables
1346  *
1347  * Assume the expression is f(x), where x are original (i.e., not auxiliary) variables.
1348  * Assume that f(x) is associated with auxiliary variable z.
1349  *
1350  * If there are negative locks, then return the violation of z &le; f(x) and sets `violover` to TRUE.
1351  * If there are positive locks, then return the violation of z &ge; f(x) and sets `violunder` to TRUE.
1352  * Of course, if there both negative and positive locks, then return the violation of z = f(x).
1353  * If f could not be evaluated, then return SCIPinfinity() and set both `violover` and `violunder` to TRUE.
1354  *
1355  * @note This does not reevaluate the violation, but assumes that the expression has been evaluated
1356  */
1357 static
1359  SCIP* scip, /**< SCIP data structure */
1360  SCIP_EXPR* expr, /**< expression */
1361  SCIP_SOL* sol, /**< solution that has been evaluated */
1362  SCIP_Bool* violunder, /**< buffer to store whether z >= f(x) is violated, or NULL */
1363  SCIP_Bool* violover /**< buffer to store whether z <= f(x) is violated, or NULL */
1364  )
1365 {
1366  SCIP_EXPR_OWNERDATA* ownerdata;
1367  SCIP_Real auxvarvalue;
1368 
1369  assert(expr != NULL);
1370 
1371  ownerdata = SCIPexprGetOwnerData(expr);
1372  assert(ownerdata != NULL);
1373  assert(ownerdata->auxvar != NULL);
1374 
1375  if( SCIPexprGetEvalValue(expr) == SCIP_INVALID )
1376  {
1377  if( violunder != NULL )
1378  *violunder = TRUE;
1379  if( violover != NULL )
1380  *violover = TRUE;
1381  return SCIPinfinity(scip);
1382  }
1383 
1384  auxvarvalue = SCIPgetSolVal(scip, sol, ownerdata->auxvar);
1385 
1386  if( ownerdata->nlocksneg > 0 && auxvarvalue > SCIPexprGetEvalValue(expr) )
1387  {
1388  if( violunder != NULL )
1389  *violunder = FALSE;
1390  if( violover != NULL )
1391  *violover = TRUE;
1392  return auxvarvalue - SCIPexprGetEvalValue(expr);
1393  }
1394 
1395  if( ownerdata->nlockspos > 0 && SCIPexprGetEvalValue(expr) > auxvarvalue )
1396  {
1397  if( violunder != NULL )
1398  *violunder = TRUE;
1399  if( violover != NULL )
1400  *violover = FALSE;
1401  return SCIPexprGetEvalValue(expr) - auxvarvalue;
1402  }
1403 
1404  if( violunder != NULL )
1405  *violunder = FALSE;
1406  if( violover != NULL )
1407  *violover = FALSE;
1408  return 0.0;
1409 }
1410 
1411 /** returns absolute violation for auxvar relation in an expression w.r.t. auxiliary variables
1412  *
1413  * Assume the expression is f(w), where w are auxiliary variables that were introduced by some nlhdlr.
1414  * Assume that f(w) is associated with auxiliary variable z.
1415  *
1416  * If there are negative locks, then return the violation of z &le; f(w) and sets `violover` to TRUE.
1417  * If there are positive locks, then return the violation of z &ge; f(w) and sets `violunder` to TRUE.
1418  * Of course, if there both negative and positive locks, then return the violation of z = f(w).
1419  * If f could not be evaluated, then return SCIPinfinity() and set both `violover` and `violunder` to TRUE.
1420  *
1421  * @note This does not reevaluate the violation, but assumes that f(w) is passed in with auxvalue.
1422  */
1423 static
1425  SCIP* scip, /**< SCIP data structure */
1426  SCIP_EXPR* expr, /**< expression */
1427  SCIP_Real auxvalue, /**< value of f(w) */
1428  SCIP_SOL* sol, /**< solution that has been evaluated */
1429  SCIP_Bool* violunder, /**< buffer to store whether z >= f(w) is violated, or NULL */
1430  SCIP_Bool* violover /**< buffer to store whether z <= f(w) is violated, or NULL */
1431  )
1432 {
1433  SCIP_EXPR_OWNERDATA* ownerdata;
1434  SCIP_Real auxvarvalue;
1435 
1436  assert(expr != NULL);
1437 
1438  ownerdata = SCIPexprGetOwnerData(expr);
1439  assert(ownerdata != NULL);
1440  assert(ownerdata->auxvar != NULL);
1441 
1442  if( auxvalue == SCIP_INVALID )
1443  {
1444  if( violunder != NULL )
1445  *violunder = TRUE;
1446  if( violover != NULL )
1447  *violover = TRUE;
1448  return SCIPinfinity(scip);
1449  }
1450 
1451  auxvarvalue = SCIPgetSolVal(scip, sol, ownerdata->auxvar);
1452 
1453  if( ownerdata->nlocksneg > 0 && auxvarvalue > auxvalue )
1454  {
1455  if( violunder != NULL )
1456  *violunder = FALSE;
1457  if( violover != NULL )
1458  *violover = TRUE;
1459  return auxvarvalue - auxvalue;
1460  }
1461 
1462  if( ownerdata->nlockspos > 0 && auxvalue > auxvarvalue )
1463  {
1464  if( violunder != NULL )
1465  *violunder = TRUE;
1466  if( violover != NULL )
1467  *violover = FALSE;
1468  return auxvalue - auxvarvalue;
1469  }
1470 
1471  if( violunder != NULL )
1472  *violunder = FALSE;
1473  if( violover != NULL )
1474  *violover = FALSE;
1475 
1476  return 0.0;
1477 }
1478 
1479 /** computes violation of a constraint */
1480 static
1482  SCIP* scip, /**< SCIP data structure */
1483  SCIP_CONS* cons, /**< constraint */
1484  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
1485  SCIP_Longint soltag /**< tag that uniquely identifies the solution (with its values), or 0. */
1486  )
1487 {
1488  SCIP_CONSDATA* consdata;
1489  SCIP_Real activity;
1490 
1491  assert(scip != NULL);
1492  assert(cons != NULL);
1493 
1494  consdata = SCIPconsGetData(cons);
1495  assert(consdata != NULL);
1496 
1497  SCIP_CALL( SCIPevalExpr(scip, consdata->expr, sol, soltag) );
1498  activity = SCIPexprGetEvalValue(consdata->expr);
1499 
1500  /* consider constraint as violated if it is undefined in the current point */
1501  if( activity == SCIP_INVALID )
1502  {
1503  consdata->lhsviol = SCIPinfinity(scip);
1504  consdata->rhsviol = SCIPinfinity(scip);
1505  return SCIP_OKAY;
1506  }
1507 
1508  /* compute violations */
1509  consdata->lhsviol = SCIPisInfinity(scip, -consdata->lhs) ? -SCIPinfinity(scip) : consdata->lhs - activity;
1510  consdata->rhsviol = SCIPisInfinity(scip, consdata->rhs) ? -SCIPinfinity(scip) : activity - consdata->rhs;
1511 
1512  return SCIP_OKAY;
1513 }
1514 
1515 /** returns absolute violation of a constraint
1516  *
1517  * @note This does not reevaluate the violation, but assumes that computeViolation() has been called before.
1518  */
1519 static
1521  SCIP_CONS* cons /**< constraint */
1522  )
1523 {
1524  SCIP_CONSDATA* consdata;
1525 
1526  assert(cons != NULL);
1527 
1528  consdata = SCIPconsGetData(cons);
1529  assert(consdata != NULL);
1530 
1531  return MAX3(0.0, consdata->lhsviol, consdata->rhsviol);
1532 }
1533 
1534 /** computes relative violation of a constraint
1535  *
1536  * @note This does not reevaluate the violation, but assumes that computeViolation() has been called before.
1537  */
1538 static
1540  SCIP* scip, /**< SCIP data structure */
1541  SCIP_CONS* cons, /**< constraint */
1542  SCIP_Real* viol, /**< buffer to store violation */
1543  SCIP_SOL* sol, /**< solution or NULL if LP solution should be used */
1544  SCIP_Longint soltag /**< tag that uniquely identifies the solution (with its values), or 0 */
1545  )
1546 {
1547  SCIP_CONSHDLR* conshdlr;
1548  SCIP_CONSHDLRDATA* conshdlrdata;
1549  SCIP_CONSDATA* consdata;
1550  SCIP_Real scale;
1551 
1552  assert(cons != NULL);
1553  assert(viol != NULL);
1554 
1555  conshdlr = SCIPconsGetHdlr(cons);
1556  assert(conshdlr != NULL);
1557 
1558  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1559  assert(conshdlrdata != NULL);
1560 
1561  *viol = getConsAbsViolation(cons);
1562 
1563  if( conshdlrdata->violscale == 'n' )
1564  return SCIP_OKAY;
1565 
1566  if( SCIPisInfinity(scip, *viol) )
1567  return SCIP_OKAY;
1568 
1569  consdata = SCIPconsGetData(cons);
1570  assert(consdata != NULL);
1571 
1572  if( conshdlrdata->violscale == 'a' )
1573  {
1574  scale = MAX(1.0, REALABS(SCIPexprGetEvalValue(consdata->expr)));
1575 
1576  /* consider value of side that is violated for scaling, too */
1577  if( consdata->lhsviol > 0.0 && REALABS(consdata->lhs) > scale )
1578  {
1579  assert(!SCIPisInfinity(scip, -consdata->lhs));
1580  scale = REALABS(consdata->lhs);
1581  }
1582  else if( consdata->rhsviol > 0.0 && REALABS(consdata->rhs) > scale )
1583  {
1584  assert(!SCIPisInfinity(scip, consdata->rhs));
1585  scale = REALABS(consdata->rhs);
1586  }
1587 
1588  *viol /= scale;
1589  return SCIP_OKAY;
1590  }
1591 
1592  /* if not 'n' or 'a', then it has to be 'g' at the moment */
1593  assert(conshdlrdata->violscale == 'g');
1594  if( soltag == 0L || consdata->gradnormsoltag != soltag )
1595  {
1596  /* we need the varexprs to conveniently access the gradient */
1597  SCIP_CALL( storeVarExprs(scip, conshdlr, consdata) );
1598 
1599  /* update cached value of norm of gradient */
1600  consdata->gradnorm = 0.0;
1601 
1602  /* compute gradient */
1603  SCIP_CALL( SCIPevalExprGradient(scip, consdata->expr, sol, soltag) );
1604 
1605  /* gradient evaluation error -> no scaling */
1606  if( SCIPexprGetDerivative(consdata->expr) != SCIP_INVALID )
1607  {
1608  int i;
1609  for( i = 0; i < consdata->nvarexprs; ++i )
1610  {
1611  SCIP_Real deriv;
1612 
1613  assert(SCIPexprGetDiffTag(consdata->expr) == SCIPexprGetDiffTag(consdata->varexprs[i]));
1614  deriv = SCIPexprGetDerivative(consdata->varexprs[i]);
1615  if( deriv == SCIP_INVALID )
1616  {
1617  /* SCIPdebugMsg(scip, "gradient evaluation error for component %d\n", i); */
1618  consdata->gradnorm = 0.0;
1619  break;
1620  }
1621 
1622  consdata->gradnorm += deriv*deriv;
1623  }
1624  }
1625  consdata->gradnorm = sqrt(consdata->gradnorm);
1626  consdata->gradnormsoltag = soltag;
1627  }
1628 
1629  *viol /= MAX(1.0, consdata->gradnorm);
1630 
1631  return SCIP_OKAY;
1632 }
1633 
1634 /** returns whether constraint is currently violated
1635  *
1636  * @note This does not reevaluate the violation, but assumes that computeViolation() has been called before.
1637  */
1638 static
1640  SCIP* scip, /**< SCIP data structure */
1641  SCIP_CONS* cons /**< constraint */
1642  )
1643 {
1644  return getConsAbsViolation(cons) > SCIPfeastol(scip);
1645 }
1646 
1647 /** checks for a linear variable that can be increased or decreased without harming feasibility */
1648 static
1650  SCIP* scip, /**< SCIP data structure */
1651  SCIP_CONS* cons /**< constraint */
1652  )
1653 {
1654  SCIP_CONSDATA* consdata;
1655  int poslock;
1656  int neglock;
1657  int i;
1658 
1659  assert(cons != NULL);
1660 
1661  consdata = SCIPconsGetData(cons);
1662  assert(consdata != NULL);
1663 
1664  consdata->linvarincr = NULL;
1665  consdata->linvardecr = NULL;
1666  consdata->linvarincrcoef = 0.0;
1667  consdata->linvardecrcoef = 0.0;
1668 
1669  /* root expression is not a sum -> no unlocked linear variable available */
1670  if( !SCIPisExprSum(scip, consdata->expr) )
1671  return;
1672 
1673  for( i = 0; i < SCIPexprGetNChildren(consdata->expr); ++i )
1674  {
1675  SCIP_EXPR* child;
1676 
1677  child = SCIPexprGetChildren(consdata->expr)[i];
1678  assert(child != NULL);
1679 
1680  /* check whether the child is a variable expression */
1681  if( SCIPisExprVar(scip, child) )
1682  {
1683  SCIP_VAR* var = SCIPgetVarExprVar(child);
1684  SCIP_Real coef = SCIPgetCoefsExprSum(consdata->expr)[i];
1685 
1686  if( coef > 0.0 )
1687  {
1688  poslock = !SCIPisInfinity(scip, consdata->rhs) ? 1 : 0;
1689  neglock = !SCIPisInfinity(scip, -consdata->lhs) ? 1 : 0;
1690  }
1691  else
1692  {
1693  poslock = !SCIPisInfinity(scip, -consdata->lhs) ? 1 : 0;
1694  neglock = !SCIPisInfinity(scip, consdata->rhs) ? 1 : 0;
1695  }
1697 
1698  if( SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) - neglock == 0 )
1699  {
1700  /* for a*x + f(y) \in [lhs, rhs], we can decrease x without harming other constraints
1701  * if we have already one candidate, then take the one where the loss in the objective function is less
1702  */
1703  if( (consdata->linvardecr == NULL) ||
1704  (SCIPvarGetObj(consdata->linvardecr) / consdata->linvardecrcoef > SCIPvarGetObj(var) / coef) )
1705  {
1706  consdata->linvardecr = var;
1707  consdata->linvardecrcoef = coef;
1708  }
1709  }
1710 
1711  if( SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) - poslock == 0 )
1712  {
1713  /* for a*x + f(y) \in [lhs, rhs], we can increase x without harm
1714  * if we have already one candidate, then take the one where the loss in the objective function is less
1715  */
1716  if( (consdata->linvarincr == NULL) ||
1717  (SCIPvarGetObj(consdata->linvarincr) / consdata->linvarincrcoef > SCIPvarGetObj(var) / coef) )
1718  {
1719  consdata->linvarincr = var;
1720  consdata->linvarincrcoef = coef;
1721  }
1722  }
1723  }
1724  }
1725 
1726  assert(consdata->linvarincr == NULL || consdata->linvarincrcoef != 0.0);
1727  assert(consdata->linvardecr == NULL || consdata->linvardecrcoef != 0.0);
1728 
1729  if( consdata->linvarincr != NULL )
1730  {
1731  SCIPdebugMsg(scip, "may increase <%s> to become feasible\n", SCIPvarGetName(consdata->linvarincr));
1732  }
1733  if( consdata->linvardecr != NULL )
1734  {
1735  SCIPdebugMsg(scip, "may decrease <%s> to become feasible\n", SCIPvarGetName(consdata->linvardecr));
1736  }
1737 }
1738 
1739 /** Given a solution where every nonlinear constraint is either feasible or can be made feasible by
1740  * moving a linear variable, construct the corresponding feasible solution and pass it to the trysol heuristic.
1741  *
1742  * The method assumes that this is always possible and that not all constraints are feasible already.
1743  */
1744 static
1746  SCIP* scip, /**< SCIP data structure */
1747  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1748  SCIP_CONS** conss, /**< constraints to process */
1749  int nconss, /**< number of constraints */
1750  SCIP_SOL* sol, /**< solution to process */
1751  SCIP_Bool* success /**< buffer to store whether we succeeded to construct a solution that satisfies all provided constraints */
1752  )
1753 {
1754  SCIP_CONSHDLRDATA* conshdlrdata;
1755  SCIP_SOL* newsol;
1756  int c;
1757 
1758  assert(scip != NULL);
1759  assert(conshdlr != NULL);
1760  assert(conss != NULL || nconss == 0);
1761  assert(success != NULL);
1762 
1763  *success = FALSE;
1764 
1765  /* don't propose new solutions if not in presolve or solving */
1767  return SCIP_OKAY;
1768 
1769  conshdlrdata = SCIPconshdlrGetData(conshdlr);
1770  assert(conshdlrdata != NULL);
1771 
1772  if( sol != NULL )
1773  {
1774  SCIP_CALL( SCIPcreateSolCopy(scip, &newsol, sol) );
1775  }
1776  else
1777  {
1778  SCIP_CALL( SCIPcreateLPSol(scip, &newsol, NULL) );
1779  }
1780  SCIP_CALL( SCIPunlinkSol(scip, newsol) );
1781  SCIPdebugMsg(scip, "attempt to make solution from <%s> feasible by shifting linear variable\n",
1782  sol != NULL ? (SCIPsolGetHeur(sol) != NULL ? SCIPheurGetName(SCIPsolGetHeur(sol)) : "tree") : "LP");
1783 
1784  for( c = 0; c < nconss; ++c )
1785  {
1786  SCIP_CONSDATA* consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
1787  SCIP_Real viol = 0.0;
1788  SCIP_Real delta;
1789  SCIP_Real gap;
1790 
1791  assert(consdata != NULL);
1792 
1793  /* get absolute violation and sign */
1794  if( consdata->lhsviol > SCIPfeastol(scip) )
1795  viol = consdata->lhsviol; /* lhs - activity */
1796  else if( consdata->rhsviol > SCIPfeastol(scip) )
1797  viol = -consdata->rhsviol; /* rhs - activity */
1798  else
1799  continue; /* constraint is satisfied */
1800 
1801  if( consdata->linvarincr != NULL &&
1802  ((viol > 0.0 && consdata->linvarincrcoef > 0.0) || (viol < 0.0 && consdata->linvarincrcoef < 0.0)) )
1803  {
1804  SCIP_VAR* var = consdata->linvarincr;
1805 
1806  /* compute how much we would like to increase var */
1807  delta = viol / consdata->linvarincrcoef;
1808  assert(delta > 0.0);
1809 
1810  /* if var has an upper bound, may need to reduce delta */
1811  if( !SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) )
1812  {
1813  gap = SCIPvarGetUbGlobal(var) - SCIPgetSolVal(scip, newsol, var);
1814  delta = MIN(MAX(0.0, gap), delta);
1815  }
1816  if( SCIPisPositive(scip, delta) )
1817  {
1818  /* if variable is integral, round delta up so that it will still have an integer value */
1819  if( SCIPvarIsIntegral(var) )
1820  delta = SCIPceil(scip, delta);
1821 
1822  SCIP_CALL( SCIPincSolVal(scip, newsol, var, delta) );
1823  SCIPdebugMsg(scip, "increase <%s> by %g to %g to remedy lhs-violation %g of cons <%s>\n",
1824  SCIPvarGetName(var), delta, SCIPgetSolVal(scip, newsol, var), viol, SCIPconsGetName(conss[c])); /*lint !e613*/
1825 
1826  /* adjust constraint violation, if satisfied go on to next constraint */
1827  viol -= consdata->linvarincrcoef * delta;
1828  if( SCIPisZero(scip, viol) )
1829  continue;
1830  }
1831  }
1832 
1833  assert(viol != 0.0);
1834  if( consdata->linvardecr != NULL &&
1835  ((viol > 0.0 && consdata->linvardecrcoef < 0.0) || (viol < 0.0 && consdata->linvardecrcoef > 0.0)) )
1836  {
1837  SCIP_VAR* var = consdata->linvardecr;
1838 
1839  /* compute how much we would like to decrease var */
1840  delta = viol / consdata->linvardecrcoef;
1841  assert(delta < 0.0);
1842 
1843  /* if var has a lower bound, may need to reduce delta */
1844  if( !SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)) )
1845  {
1846  gap = SCIPgetSolVal(scip, newsol, var) - SCIPvarGetLbGlobal(var);
1847  delta = MAX(MIN(0.0, gap), delta);
1848  }
1849  if( SCIPisNegative(scip, delta) )
1850  {
1851  /* if variable is integral, round delta down so that it will still have an integer value */
1852  if( SCIPvarIsIntegral(var) )
1853  delta = SCIPfloor(scip, delta);
1854  SCIP_CALL( SCIPincSolVal(scip, newsol, consdata->linvardecr, delta) );
1855  /*lint --e{613} */
1856  SCIPdebugMsg(scip, "increase <%s> by %g to %g to remedy rhs-violation %g of cons <%s>\n",
1857  SCIPvarGetName(var), delta, SCIPgetSolVal(scip, newsol, var), viol, SCIPconsGetName(conss[c]));
1858 
1859  /* adjust constraint violation, if satisfied go on to next constraint */
1860  viol -= consdata->linvardecrcoef * delta;
1861  if( SCIPisZero(scip, viol) )
1862  continue;
1863  }
1864  }
1865 
1866  /* still here... so probably we could not make constraint feasible due to variable bounds, thus give up */
1867  break;
1868  }
1869 
1870  /* if we have a solution that should satisfy all quadratic constraints and has a better objective than the current upper bound,
1871  * then pass it to the trysol heuristic
1872  */
1873  if( c == nconss && (SCIPisInfinity(scip, SCIPgetUpperbound(scip)) || SCIPisSumLT(scip, SCIPgetSolTransObj(scip, newsol), SCIPgetUpperbound(scip))) )
1874  {
1875  SCIPdebugMsg(scip, "pass solution with objective val %g to trysol heuristic\n", SCIPgetSolTransObj(scip, newsol));
1876 
1877  assert(conshdlrdata->trysolheur != NULL);
1878  SCIP_CALL( SCIPheurPassSolTrySol(scip, conshdlrdata->trysolheur, newsol) );
1879 
1880  *success = TRUE;
1881  }
1882 
1883  SCIP_CALL( SCIPfreeSol(scip, &newsol) );
1884 
1885  return SCIP_OKAY;
1886 }
1887 
1888 /** adds globally valid tight estimators in a given solution as cut to cutpool
1889  *
1890  * Called by addTightEstimatorCuts() for a specific expression, nlhdlr, and estimate-direction (over or under).
1891  */
1892 static
1894  SCIP* scip, /**< SCIP data structure */
1895  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1896  SCIP_CONS* cons, /**< constraint */
1897  SCIP_EXPR* expr, /**< expression */
1898  EXPRENFO* exprenfo, /**< expression enfo data, e.g., nlhdlr to use */
1899  SCIP_SOL* sol, /**< reference point where to estimate */
1900  SCIP_Bool overestimate, /**< whether to overestimate */
1901  SCIP_PTRARRAY* rowpreps /**< array for rowpreps */
1902  )
1903 {
1904  SCIP_Bool estimatesuccess = FALSE;
1905  SCIP_Bool branchscoresuccess = FALSE;
1906  int minidx;
1907  int maxidx;
1908  int r;
1909 
1910  assert(scip != NULL);
1911  assert(expr != NULL);
1912  assert(exprenfo != NULL);
1913  assert(rowpreps != NULL);
1914 
1915  ENFOLOG( SCIPinfoMessage(scip, enfologfile, " %sestimate using nlhdlr <%s> for expr %p (%s)\n",
1916  overestimate ? "over" : "under", SCIPnlhdlrGetName(exprenfo->nlhdlr), (void*)expr, SCIPexprhdlrGetName(SCIPexprGetHdlr(expr))); )
1917 
1918  SCIP_CALL( SCIPnlhdlrEstimate(scip, conshdlr, exprenfo->nlhdlr, expr, exprenfo->nlhdlrexprdata, sol,
1919  exprenfo->auxvalue, overestimate, SCIPinfinity(scip), FALSE, rowpreps, &estimatesuccess, &branchscoresuccess) );
1920 
1921  minidx = SCIPgetPtrarrayMinIdx(scip, rowpreps);
1922  maxidx = SCIPgetPtrarrayMaxIdx(scip, rowpreps);
1923  assert(estimatesuccess == (minidx <= maxidx));
1924 
1925  if( !estimatesuccess )
1926  {
1927  ENFOLOG( SCIPinfoMessage(scip, enfologfile, " estimate of nlhdlr %s failed\n", SCIPnlhdlrGetName(exprenfo->nlhdlr)); )
1928  return SCIP_OKAY;
1929  }
1930 
1931  for( r = minidx; r <= maxidx; ++r )
1932  {
1933  SCIP_ROWPREP* rowprep;
1934  SCIP_ROW* row;
1935  SCIP_Real estimateval;
1936  int i;
1937 
1938  rowprep = (SCIP_ROWPREP*) SCIPgetPtrarrayVal(scip, rowpreps, r);
1939  assert(rowprep != NULL);
1940  assert(SCIProwprepGetSidetype(rowprep) == (overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT));
1941 
1942  /* if estimators is only local valid, then skip */
1943  if( SCIProwprepIsLocal(rowprep) )
1944  {
1945  ENFOLOG( SCIPinfoMessage(scip, enfologfile, " skip local estimator\n"); )
1946  SCIPfreeRowprep(scip, &rowprep);
1947  continue;
1948  }
1949 
1950  /* compute value of estimator */
1951  estimateval = -SCIProwprepGetSide(rowprep);
1952  for( i = 0; i < SCIProwprepGetNVars(rowprep); ++i )
1953  estimateval += SCIProwprepGetCoefs(rowprep)[i] * SCIPgetSolVal(scip, sol, SCIProwprepGetVars(rowprep)[i]);
1954 
1955  /* if estimator value is not tight (or even "more than tight", e.g., when estimating in integer vars), then skip */
1956  if( (overestimate && !SCIPisFeasLE(scip, estimateval, SCIPexprGetEvalValue(expr))) ||
1957  (!overestimate && !SCIPisFeasGE(scip, estimateval, SCIPexprGetEvalValue(expr))) )
1958  {
1959  ENFOLOG( SCIPinfoMessage(scip, enfologfile, " skip non-tight estimator with value %g, expr value %g\n", estimateval, SCIPexprGetEvalValue(expr)); )
1960  SCIPfreeRowprep(scip, &rowprep);
1961  continue;
1962  }
1963 
1964  /* complete estimator to cut and clean it up */
1965  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) );
1966  SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPinfinity(scip), &estimatesuccess) );
1967 
1968  /* if cleanup failed or rowprep is local now, then skip */
1969  if( !estimatesuccess || SCIProwprepIsLocal(rowprep) )
1970  {
1971  ENFOLOG( SCIPinfoMessage(scip, enfologfile, " skip after cleanup failed or made estimator locally valid\n"); )
1972  SCIPfreeRowprep(scip, &rowprep);
1973  continue;
1974  }
1975 
1976  /* generate row and add to cutpool */
1977  SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );
1978 
1979  ENFOLOG( SCIPinfoMessage(scip, enfologfile, " adding cut ");
1980  SCIP_CALL( SCIPprintRow(scip, row, enfologfile) ); )
1981 
1982  SCIP_CALL( SCIPaddPoolCut(scip, row) );
1983  /* SCIPnlhdlrIncrementNSeparated(nlhdlr); */
1984 
1985  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1986  SCIPfreeRowprep(scip, &rowprep);
1987  }
1988 
1989  SCIP_CALL( SCIPclearPtrarray(scip, rowpreps) );
1990 
1991  return SCIP_OKAY;
1992 }
1993 
1994 /** adds globally valid tight estimators in a given solution as cuts to cutpool
1995  *
1996  * Essentially we want to ensure that the LP relaxation is tight in the new solution, if possible.
1997  * For convex constraints, we would achieve this by linearizing.
1998  * To avoid checking explicitly for convexity, we compute estimators via any nlhdlr that didn't say it would
1999  * use bound information and check whether the estimator is tight.
2000  *
2001  * Since linearization may happen in auxiliary variables, we ensure that auxiliary variables are set
2002  * to the eval-value of its expression, i.e., we change sol so it is also feasible in the extended formulation.
2003  */
2004 static
2006  SCIP* scip, /**< SCIP data structure */
2007  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2008  SCIP_CONS** conss, /**< constraints */
2009  int nconss, /**< number of constraints */
2010  SCIP_SOL* sol /**< reference point where to estimate */
2011  )
2012 {
2013  SCIP_CONSDATA* consdata;
2014  SCIP_Longint soltag;
2015  SCIP_EXPRITER* it;
2016  SCIP_EXPR* expr;
2017  SCIP_PTRARRAY* rowpreps;
2018  int c, e;
2019 
2020  assert(scip != NULL);
2021  assert(conshdlr != NULL);
2022  assert(conss != NULL || nconss == 0);
2023 
2024  ENFOLOG( SCIPinfoMessage(scip, enfologfile, "add tight estimators in new solution from <%s> to cutpool\n", SCIPheurGetName(SCIPsolGetHeur(sol))); )
2025 
2026  /* TODO probably we just evaluated all expressions when checking the sol before it was added
2027  * would be nice to recognize this and skip reevaluating
2028  */
2029  soltag = SCIPgetExprNewSoltag(scip);
2030 
2031  SCIP_CALL( SCIPcreatePtrarray(scip, &rowpreps) );
2032 
2033  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
2036 
2037  for( c = 0; c < nconss; ++c )
2038  {
2039  /* skip constraints that are not enabled or deleted or have separation disabled */
2040  if( !SCIPconsIsEnabled(conss[c]) || SCIPconsIsDeleted(conss[c]) || !SCIPconsIsSeparationEnabled(conss[c]) )
2041  continue;
2042  assert(SCIPconsIsActive(conss[c]));
2043 
2044  consdata = SCIPconsGetData(conss[c]);
2045  assert(consdata != NULL);
2046 
2047  /* TODO we could remember for which constraints there is a chance that we would add anything,
2048  * i.e., there is some convex-like expression, and skip other constraints
2049  */
2050 
2051  ENFOLOG(
2052  {
2053  int i;
2054  SCIPinfoMessage(scip, enfologfile, " constraint ");
2055  SCIP_CALL( SCIPprintCons(scip, conss[c], enfologfile) );
2056  SCIPinfoMessage(scip, enfologfile, "\n and point\n");
2057  for( i = 0; i < consdata->nvarexprs; ++i )
2058  {
2059  SCIP_VAR* var;
2060  var = SCIPgetVarExprVar(consdata->varexprs[i]);
2061  SCIPinfoMessage(scip, enfologfile, " %-10s = %15g bounds: [%15g,%15g]\n", SCIPvarGetName(var),
2062  SCIPgetSolVal(scip, sol, var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
2063  }
2064  })
2065 
2066  SCIP_CALL( SCIPevalExpr(scip, consdata->expr, sol, soltag) );
2067  assert(SCIPexprGetEvalValue(consdata->expr) != SCIP_INVALID);
2068 
2069  for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
2070  {
2071  SCIP_EXPR_OWNERDATA* ownerdata;
2072 
2073  ownerdata = SCIPexprGetOwnerData(expr);
2074  assert(ownerdata != NULL);
2075 
2076  /* we can only generate a cut from an estimator if there is an auxvar */
2077  if( ownerdata->auxvar == NULL )
2078  continue;
2079 
2080  /* set value for auxvar in sol to value of expr, in case it is used to compute estimators higher up of this expression */
2081  assert(SCIPexprGetEvalTag(expr) == soltag);
2082  assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID);
2083  SCIP_CALL( SCIPsetSolVal(scip, sol, ownerdata->auxvar, SCIPexprGetEvalValue(expr)) );
2084 
2085  /* generate cuts from estimators of each nonlinear handler that provides estimates */
2086  for( e = 0; e < ownerdata->nenfos; ++e )
2087  {
2088  SCIP_NLHDLR* nlhdlr;
2089 
2090  nlhdlr = ownerdata->enfos[e]->nlhdlr;
2091  assert(nlhdlr != NULL);
2092 
2093  /* skip nlhdlr that does not implement estimate (so it does enfo) */
2094  if( !SCIPnlhdlrHasEstimate(nlhdlr) )
2095  continue;
2096 
2097  /* skip nlhdlr that does not participate in separation or looks like it would give only locally-valid estimators
2098  * (because it uses activities on vars/auxvars)
2099  */
2100  if( ((ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPAABOVE) == 0 || ownerdata->enfos[e]->sepaaboveusesactivity) &&
2101  ((ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPABELOW) == 0 || ownerdata->enfos[e]->sepabelowusesactivity) )
2102  continue;
2103 
2104  /* skip nlhdlr_default on sum, as the estimator doesn't depend on the reference point (expr is linear in auxvars) */
2105  if( SCIPisExprSum(scip, expr) && strcmp(SCIPnlhdlrGetName(nlhdlr), "default") == 0 )
2106  continue;
2107 
2108  /* evaluate the expression w.r.t. the nlhdlrs auxiliary variables, since some nlhdlr expect this before their estimate is called */
2109  SCIP_CALL( SCIPnlhdlrEvalaux(scip, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata, &ownerdata->enfos[e]->auxvalue, sol) );
2110  ENFOLOG(
2111  SCIPinfoMessage(scip, enfologfile, " expr ");
2112  SCIPprintExpr(scip, expr, enfologfile);
2113  SCIPinfoMessage(scip, enfologfile, " (%p): evalvalue %.15g auxvarvalue %.15g, nlhdlr <%s> auxvalue: %.15g\n",
2114  (void*)expr, SCIPexprGetEvalValue(expr), SCIPgetSolVal(scip, sol, ownerdata->auxvar), SCIPnlhdlrGetName(nlhdlr), ownerdata->enfos[e]->auxvalue);
2115  )
2116  /* due to setting values of auxvars to expr values in sol, the auxvalue should equal to expr evalvalue */
2117  assert(SCIPisEQ(scip, ownerdata->enfos[e]->auxvalue, SCIPexprGetEvalValue(expr)));
2118 
2119  /* if nlhdlr wants to be called for overestimate and does not use local bounds, then call estimate of nlhdlr */
2120  if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPAABOVE) && !ownerdata->enfos[e]->sepaaboveusesactivity )
2121  {
2122  SCIP_CALL( addTightEstimatorCut(scip, conshdlr, conss[c], expr, ownerdata->enfos[e], sol, TRUE, rowpreps) );
2123  }
2124 
2125  /* if nlhdlr wants to be called for underestimate and does not use local bounds, then call estimate of nlhdlr */
2126  if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_SEPABELOW) && !ownerdata->enfos[e]->sepabelowusesactivity )
2127  {
2128  SCIP_CALL( addTightEstimatorCut(scip, conshdlr, conss[c], expr, ownerdata->enfos[e], sol, FALSE, rowpreps) );
2129  }
2130  }
2131  }
2132  }
2133 
2134  SCIPfreeExpriter(&it);
2135  SCIP_CALL( SCIPfreePtrarray(scip, &rowpreps) );
2136 
2137  return SCIP_OKAY;
2138 }
2139 
2140 /** processes the event that a new primal solution has been found */
2141 static
2142 SCIP_DECL_EVENTEXEC(processNewSolutionEvent)
2144  SCIP_CONSHDLR* conshdlr;
2145  SCIP_CONSHDLRDATA* conshdlrdata;
2146  SCIP_SOL* sol;
2147 
2148  assert(scip != NULL);
2149  assert(event != NULL);
2150  assert(eventdata != NULL);
2151  assert(eventhdlr != NULL);
2152  assert(SCIPeventGetType(event) & SCIP_EVENTTYPE_SOLFOUND);
2153 
2154  conshdlr = (SCIP_CONSHDLR*)eventdata;
2155 
2156  if( SCIPconshdlrGetNConss(conshdlr) == 0 )
2157  return SCIP_OKAY;
2158 
2159  sol = SCIPeventGetSol(event);
2160  assert(sol != NULL);
2161 
2162  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2163  assert(conshdlrdata != NULL);
2164 
2165  /* we are only interested in solution coming from some heuristic other than trysol, but not from the tree
2166  * the reason for ignoring trysol solutions is that they may come ~~from an NLP solve in sepalp, where we already added linearizations, or are~~
2167  * from the tree, but postprocessed via proposeFeasibleSolution
2168  */
2169  if( SCIPsolGetHeur(sol) == NULL || SCIPsolGetHeur(sol) == conshdlrdata->trysolheur )
2170  return SCIP_OKAY;
2171 
2172  SCIPdebugMsg(scip, "caught new sol event %" SCIP_EVENTTYPE_FORMAT " from heur <%s>\n", SCIPeventGetType(event), SCIPheurGetName(SCIPsolGetHeur(sol)));
2173 
2174  SCIP_CALL( addTightEstimatorCuts(scip, conshdlr, SCIPconshdlrGetConss(conshdlr), SCIPconshdlrGetNConss(conshdlr), sol) );
2175 
2176  return SCIP_OKAY;
2177 }
2178 
2179 /** tightens the bounds of the auxiliary variable associated with an expression (or original variable if being a variable-expression) according to given bounds
2180  *
2181  * The given bounds may very well be the exprs activity (when called from forwardPropExpr()), but can also be some
2182  * tighter bounds (when called from SCIPtightenExprIntervalNonlinear()).
2183  *
2184  * Nothing will happen if SCIP is not in presolve or solve.
2185  */
2186 static
2188  SCIP* scip, /**< SCIP data structure */
2189  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2190  SCIP_EXPR* expr, /**< expression whose auxvar is to be tightened */
2191  SCIP_INTERVAL bounds, /**< bounds to be used for tightening (must not be empty) */
2192  SCIP_Bool* cutoff, /**< buffer to store whether a cutoff was detected */
2193  int* ntightenings /**< buffer to add the total number of tightenings, or NULL */
2194  )
2195 {
2196  SCIP_VAR* var;
2197  SCIP_Bool tightenedlb;
2198  SCIP_Bool tightenedub;
2199  SCIP_Bool force;
2200 
2201  assert(scip != NULL);
2202  assert(conshdlr != NULL);
2203  assert(expr != NULL);
2204  assert(cutoff != NULL);
2205 
2206  /* the given bounds must not be empty (we could cope, but we shouldn't be called in this situation) */
2207  assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bounds));
2208 
2209  *cutoff = FALSE;
2210 
2211  /* do not tighten variable in problem stage (important for unittests)
2212  * TODO put some kind of #ifdef UNITTEST around this
2213  */
2215  return SCIP_OKAY;
2216 
2217  var = SCIPgetExprAuxVarNonlinear(expr);
2218  if( var == NULL )
2219  return SCIP_OKAY;
2220 
2221  /* force tightening if conshdlrdata says so or it would mean fixing the variable */
2222  force = SCIPconshdlrGetData(conshdlr)->forceboundtightening || SCIPisEQ(scip, bounds.inf, bounds.sup);
2223 
2224  /* try to tighten lower bound of (auxiliary) variable */
2225  SCIP_CALL( SCIPtightenVarLb(scip, var, bounds.inf, force, cutoff, &tightenedlb) );
2226  if( tightenedlb )
2227  {
2228  if( ntightenings != NULL )
2229  ++*ntightenings;
2230  SCIPdebugMsg(scip, "tightened lb on auxvar <%s> to %.15g (forced:%u)\n", SCIPvarGetName(var), SCIPvarGetLbLocal(var), force);
2231  }
2232  if( *cutoff )
2233  {
2234  SCIPdebugMsg(scip, "cutoff when tightening lb on auxvar <%s> to %.15g\n", SCIPvarGetName(var), bounds.inf);
2235  return SCIP_OKAY;
2236  }
2237 
2238  /* try to tighten upper bound of (auxiliary) variable */
2239  SCIP_CALL( SCIPtightenVarUb(scip, var, bounds.sup, force, cutoff, &tightenedub) );
2240  if( tightenedub )
2241  {
2242  if( ntightenings != NULL )
2243  ++*ntightenings;
2244  SCIPdebugMsg(scip, "tightened ub on auxvar <%s> to %.15g (forced:%u)\n", SCIPvarGetName(var), SCIPvarGetUbLocal(var), force);
2245  }
2246  if( *cutoff )
2247  {
2248  SCIPdebugMsg(scip, "cutoff when tightening ub on auxvar <%s> to %.15g\n", SCIPvarGetName(var), bounds.sup);
2249  return SCIP_OKAY;
2250  }
2251 
2252  /* TODO expr->activity should have been reevaluated now due to boundchange-events, but it used to relax bounds
2253  * that seems unnecessary and we could easily undo this here, e.g.,
2254  * if( tightenedlb ) expr->activity.inf = bounds.inf
2255  */
2256 
2257  return SCIP_OKAY;
2258 }
2259 
2260 /** propagate bounds of the expressions in a given expression tree (that is, updates activity intervals)
2261  * and tries to tighten the bounds of the auxiliary variables accordingly
2262  */
2263 static
2265  SCIP* scip, /**< SCIP data structure */
2266  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2267  SCIP_EXPR* rootexpr, /**< expression */
2268  SCIP_Bool tightenauxvars, /**< should the bounds of auxiliary variables be tightened? */
2269  SCIP_Bool* infeasible, /**< buffer to store whether the problem is infeasible (NULL if not needed) */
2270  int* ntightenings /**< buffer to store the number of auxiliary variable tightenings (NULL if not needed) */
2271  )
2272 {
2273  SCIP_EXPRITER* it;
2274  SCIP_EXPR* expr;
2275  SCIP_EXPR_OWNERDATA* ownerdata;
2276  SCIP_CONSHDLRDATA* conshdlrdata;
2277 
2278  assert(scip != NULL);
2279  assert(rootexpr != NULL);
2280 
2281  if( infeasible != NULL )
2282  *infeasible = FALSE;
2283  if( ntightenings != NULL )
2284  *ntightenings = 0;
2285 
2286  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2287  assert(conshdlrdata != NULL);
2288 
2289  /* if value is valid and empty, then we cannot improve, so do nothing */
2290  if( SCIPexprGetActivityTag(rootexpr) >= conshdlrdata->lastboundrelax && SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(rootexpr)) )
2291  {
2292  SCIPdebugMsg(scip, "stored activity of root expr is empty and valid (activitytag >= lastboundrelax (%" SCIP_LONGINT_FORMAT ")), skip forwardPropExpr -> cutoff\n", conshdlrdata->lastboundrelax);
2293 
2294  if( infeasible != NULL )
2295  *infeasible = TRUE;
2296 
2297  /* just update tag to curboundstag */
2298  SCIPexprSetActivity(rootexpr, SCIPexprGetActivity(rootexpr), conshdlrdata->curboundstag);
2299 
2300  return SCIP_OKAY;
2301  }
2302 
2303  /* if value is up-to-date, then nothing to do */
2304  if( SCIPexprGetActivityTag(rootexpr) == conshdlrdata->curboundstag )
2305  {
2306  SCIPdebugMsg(scip, "activitytag of root expr equals curboundstag (%" SCIP_LONGINT_FORMAT "), skip forwardPropExpr\n", conshdlrdata->curboundstag);
2307 
2308  assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(rootexpr))); /* handled in previous if() */
2309 
2310  return SCIP_OKAY;
2311  }
2312 
2313  ownerdata = SCIPexprGetOwnerData(rootexpr);
2314  assert(ownerdata != NULL);
2315 
2316  /* if activity of rootexpr is not used, but expr participated in detect (nenfos >= 0), then we do nothing
2317  * it seems wrong to be called for such an expression (unless we are in detect at the moment), so I add a SCIPABORT()
2318  * during detect, we are in some in-between state where we may want to eval activity
2319  * on exprs that we did not notify about their activity usage
2320  */
2321  if( ownerdata->nenfos >= 0 && ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 && !conshdlrdata->indetect)
2322  {
2323 #ifdef DEBUG_PROP
2324  SCIPdebugMsg(scip, "root expr activity is not used but enfo initialized, skip inteval\n");
2325 #endif
2326  SCIPABORT();
2327  return SCIP_OKAY;
2328  }
2329 
2330  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
2331  SCIP_CALL( SCIPexpriterInit(it, rootexpr, SCIP_EXPRITER_DFS, TRUE) );
2333 
2334  for( expr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); )
2335  {
2336  switch( SCIPexpriterGetStageDFS(it) )
2337  {
2339  {
2340  /* skip child if it has been evaluated already */
2341  SCIP_EXPR* child;
2342 
2343  child = SCIPexpriterGetChildExprDFS(it);
2344  if( conshdlrdata->curboundstag == SCIPexprGetActivityTag(child) )
2345  {
2346  if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(child)) && infeasible != NULL )
2347  *infeasible = TRUE;
2348 
2349  expr = SCIPexpriterSkipDFS(it);
2350  continue;
2351  }
2352 
2353  break;
2354  }
2355 
2357  {
2358  SCIP_INTERVAL activity;
2359 
2360  /* we should not have entered this expression if its activity was already up to date */
2361  assert(SCIPexprGetActivityTag(expr) < conshdlrdata->curboundstag);
2362 
2363  ownerdata = SCIPexprGetOwnerData(expr);
2364  assert(ownerdata != NULL);
2365 
2366  /* for var exprs where varevents are catched, activity is updated immediately when the varbound has been changed
2367  * so we can assume that the activity is up to date for all these variables
2368  * UNLESS we changed the method used to evaluate activity of variable expressions
2369  * or we currently use global bounds (varevents are catched for local bound changes only)
2370  */
2371  if( SCIPisExprVar(scip, expr) && ownerdata->filterpos >= 0 &&
2372  SCIPexprGetActivityTag(expr) >= conshdlrdata->lastvaractivitymethodchange && !conshdlrdata->globalbounds )
2373  {
2374 #ifndef NDEBUG
2375  SCIP_INTERVAL exprhdlrinterval;
2376 
2377  SCIP_CALL( SCIPcallExprInteval(scip, expr, &exprhdlrinterval, conshdlrdata->intevalvar, conshdlrdata) );
2378  assert(SCIPisRelEQ(scip, exprhdlrinterval.inf, SCIPexprGetActivity(expr).inf));
2379  assert(SCIPisRelEQ(scip, exprhdlrinterval.sup, SCIPexprGetActivity(expr).sup));
2380 #endif
2381 #ifdef DEBUG_PROP
2382  SCIPdebugMsg(scip, "skip interval evaluation of expr for var <%s> [%g,%g]\n", SCIPvarGetName(SCIPgetVarExprVar(expr)), SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup);
2383 #endif
2384  SCIPexprSetActivity(expr, SCIPexprGetActivity(expr), conshdlrdata->curboundstag);
2385 
2386  break;
2387  }
2388 
2389  if( SCIPexprGetActivityTag(expr) < conshdlrdata->lastboundrelax )
2390  {
2391  /* start with entire activity if current one is invalid */
2393  }
2395  {
2396  /* If already empty, then don't try to compute even better activity.
2397  * If cons_nonlinear were alone, then we should have noted that we are infeasible
2398  * so an assert(infeasible == NULL || *infeasible) should work here.
2399  * However, after reporting a cutoff due to expr->activity being empty,
2400  * SCIP may wander to a different node and call propagation again.
2401  * If no bounds in a nonlinear constraint have been relaxed when switching nodes
2402  * (so expr->activitytag >= conshdlrdata->lastboundrelax), then
2403  * we will still have expr->activity being empty, but will have forgotten
2404  * that we found infeasibility here before (!2221#note_134120).
2405  * Therefore we just set *infeasibility=TRUE here and stop.
2406  */
2407  if( infeasible != NULL )
2408  *infeasible = TRUE;
2409  SCIPdebugMsg(scip, "expr %p already has empty activity -> cutoff\n", (void*)expr);
2410  break;
2411  }
2412  else
2413  {
2414  /* start with current activity, since it is valid */
2415  activity = SCIPexprGetActivity(expr);
2416  }
2417 
2418  /* if activity of expr is not used, but expr participated in detect (nenfos >= 0), then do nothing */
2419  if( ownerdata->nenfos >= 0 && ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 && !conshdlrdata->indetect )
2420  {
2421 #ifdef DEBUG_PROP
2422  SCIPdebugMsg(scip, "expr %p activity is not used but enfo initialized, skip inteval\n", (void*)expr);
2423 #endif
2424  break;
2425  }
2426 
2427 #ifdef DEBUG_PROP
2428  SCIPdebugMsg(scip, "interval evaluation of expr %p ", (void*)expr);
2429  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2430  SCIPdebugMsgPrint(scip, ", current activity = [%.20g, %.20g]\n", SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup);
2431 #endif
2432 
2433  /* run interval eval of nonlinear handlers or expression handler */
2434  if( ownerdata->nenfos > 0 )
2435  {
2436  SCIP_NLHDLR* nlhdlr;
2437  SCIP_INTERVAL nlhdlrinterval;
2438  int e;
2439 
2440  /* for expressions with enforcement, nlhdlrs take care of interval evaluation */
2441  for( e = 0; e < ownerdata->nenfos && !SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, activity); ++e )
2442  {
2443  /* skip nlhdlr if it does not want to participate in activity computation */
2444  if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY) == 0 )
2445  continue;
2446 
2447  nlhdlr = ownerdata->enfos[e]->nlhdlr;
2448  assert(nlhdlr != NULL);
2449 
2450  /* skip nlhdlr if it does not provide interval evaluation (so it may only provide reverse propagation) */
2451  if( !SCIPnlhdlrHasIntEval(nlhdlr) )
2452  continue;
2453 
2454  /* let nlhdlr evaluate current expression */
2455  nlhdlrinterval = activity;
2456  SCIP_CALL( SCIPnlhdlrInteval(scip, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata,
2457  &nlhdlrinterval, conshdlrdata->intevalvar, conshdlrdata) );
2458 #ifdef DEBUG_PROP
2459  SCIPdebugMsg(scip, " nlhdlr <%s>::inteval = [%.20g, %.20g]", SCIPnlhdlrGetName(nlhdlr), nlhdlrinterval.inf, nlhdlrinterval.sup);
2460 #endif
2461 
2462  /* update activity by intersecting with computed activity */
2463  SCIPintervalIntersectEps(&activity, SCIPepsilon(scip), activity, nlhdlrinterval);
2464 #ifdef DEBUG_PROP
2465  SCIPdebugMsgPrint(scip, " -> new activity: [%.20g, %.20g]\n", activity.inf, activity.sup);
2466 #endif
2467  }
2468  }
2469  else
2470  {
2471  /* for node without enforcement (before or during detect), call the callback of the exprhdlr directly */
2472  SCIP_INTERVAL exprhdlrinterval = activity;
2473  SCIP_CALL( SCIPcallExprInteval(scip, expr, &exprhdlrinterval, conshdlrdata->intevalvar, conshdlrdata) );
2474 #ifdef DEBUG_PROP
2475  SCIPdebugMsg(scip, " exprhdlr <%s>::inteval = [%.20g, %.20g]", SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), exprhdlrinterval.inf, exprhdlrinterval.sup);
2476 #endif
2477 
2478  /* update expr->activity by intersecting with computed activity */
2479  SCIPintervalIntersectEps(&activity, SCIPepsilon(scip), activity, exprhdlrinterval);
2480 #ifdef DEBUG_PROP
2481  SCIPdebugMsgPrint(scip, " -> new activity: [%.20g, %.20g]\n", activity.inf, activity.sup);
2482 #endif
2483  }
2484 
2485  /* if expression is integral, then we try to tighten the interval bounds a bit
2486  * this should undo the addition of some unnecessary safety added by use of nextafter() in interval arithmetics, e.g., when doing pow()
2487  * it would be ok to use ceil() and floor(), but for safety we use SCIPceil and SCIPfloor for now
2488  * do this only if using boundtightening-inteval and not in redundancy check (there we really want to relax all variables)
2489  * boundtightening-inteval does not relax integer variables, so can omit expressions without children
2490  * (constants should be ok, too)
2491  */
2492  if( SCIPexprIsIntegral(expr) && conshdlrdata->intevalvar == intEvalVarBoundTightening && SCIPexprGetNChildren(expr) > 0 )
2493  {
2494  if( activity.inf > -SCIP_INTERVAL_INFINITY )
2495  activity.inf = SCIPceil(scip, activity.inf);
2496  if( activity.sup < SCIP_INTERVAL_INFINITY )
2497  activity.sup = SCIPfloor(scip, activity.sup);
2498 #ifdef DEBUG_PROP
2499  SCIPdebugMsg(scip, " applying integrality: [%.20g, %.20g]\n", activity.inf, activity.sup);
2500 #endif
2501  }
2502 
2503  /* mark the current node to be infeasible if either the lower/upper bound is above/below +/- SCIPinfinity()
2504  * TODO this is a problem if dual-presolve fixed a variable to +/- infinity
2505  */
2506  if( SCIPisInfinity(scip, activity.inf) || SCIPisInfinity(scip, -activity.sup) )
2507  {
2508  SCIPdebugMsg(scip, "cut off due to activity [%g,%g] beyond infinity\n", activity.inf, activity.sup);
2509  SCIPintervalSetEmpty(&activity);
2510  }
2511 
2512  /* now finally store activity in expr */
2513  SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);
2514 
2516  {
2517  if( infeasible != NULL )
2518  *infeasible = TRUE;
2519  }
2520  else if( tightenauxvars && ownerdata->auxvar != NULL )
2521  {
2522  SCIP_Bool tighteninfeasible;
2523 
2524  SCIP_CALL( tightenAuxVarBounds(scip, conshdlr, expr, activity, &tighteninfeasible, ntightenings) );
2525  if( tighteninfeasible )
2526  {
2527  if( infeasible != NULL )
2528  *infeasible = TRUE;
2529  SCIPintervalSetEmpty(&activity);
2530  SCIPexprSetActivity(expr, activity, conshdlrdata->curboundstag);
2531  }
2532  }
2533 
2534  break;
2535  }
2536 
2537  default:
2538  /* you should never be here */
2539  SCIPerrorMessage("unexpected iterator stage\n");
2540  SCIPABORT();
2541  break;
2542  }
2543 
2544  expr = SCIPexpriterGetNext(it);
2545  }
2546 
2547  SCIPfreeExpriter(&it);
2548 
2549  return SCIP_OKAY;
2550 }
2551 
2552 /** returns whether intersecting `oldinterval` with `newinterval` would provide a properly smaller interval
2553  *
2554  * If `subsetsufficient` is TRUE, then the intersection being smaller than oldinterval is sufficient.
2555  *
2556  * If `subsetsufficient` is FALSE, then we require
2557  * - a change from an unbounded interval to a bounded one, or
2558  * - or a change from an unfixed (width > epsilon) to a fixed interval, or
2559  * - a minimal tightening of one of the interval bounds as defined by SCIPis{Lb,Ub}Better().
2560  */
2561 static
2563  SCIP* scip, /**< SCIP data structure */
2564  SCIP_Bool subsetsufficient, /**< whether the intersection being a proper subset of oldinterval is sufficient */
2565  SCIP_INTERVAL newinterval, /**< new interval */
2566  SCIP_INTERVAL oldinterval /**< old interval */
2567  )
2568 {
2569  assert(scip != NULL);
2570  assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, newinterval));
2571  assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, oldinterval));
2572 
2573  if( subsetsufficient )
2574  /* oldinterval \cap newinterval < oldinterval iff not oldinterval is subset of newinterval */
2575  return !SCIPintervalIsSubsetEQ(SCIP_INTERVAL_INFINITY, oldinterval, newinterval);
2576 
2577  /* check whether lower bound of interval becomes finite */
2578  if( oldinterval.inf <= -SCIP_INTERVAL_INFINITY && newinterval.inf > -SCIP_INTERVAL_INFINITY )
2579  return TRUE;
2580 
2581  /* check whether upper bound of interval becomes finite */
2582  if( oldinterval.sup >= SCIP_INTERVAL_INFINITY && newinterval.sup > SCIP_INTERVAL_INFINITY )
2583  return TRUE;
2584 
2585  /* check whether intersection will have width <= epsilon, if oldinterval doesn't have yet */
2586  if( !SCIPisEQ(scip, oldinterval.inf, oldinterval.sup) && SCIPisEQ(scip, MAX(oldinterval.inf, newinterval.inf), MIN(oldinterval.sup, newinterval.sup)) )
2587  return TRUE;
2588 
2589  /* check whether lower bound on interval will be better by SCIP's quality measures for boundchanges */
2590  if( SCIPisLbBetter(scip, newinterval.inf, oldinterval.inf, oldinterval.sup) )
2591  return TRUE;
2592 
2593  /* check whether upper bound on interval will be better by SCIP's quality measures for boundchanges */
2594  if( SCIPisUbBetter(scip, newinterval.sup, oldinterval.inf, oldinterval.sup) )
2595  return TRUE;
2596 
2597  return FALSE;
2598 }
2599 
2600 /** propagates bounds for each sub-expression in the `reversepropqueue` by starting from the root expressions
2601  *
2602  * The expression will be traversed in breadth first search by using this queue.
2603  *
2604  * @note Calling this function requires feasible intervals for each sub-expression; this is guaranteed by calling
2605  * forwardPropExpr() before calling this function.
2606  *
2607  * @note Calling this function with `*infeasible` = TRUE will only empty the queue.
2608  */
2609 static
2611  SCIP* scip, /**< SCIP data structure */
2612  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2613  SCIP_Bool* infeasible, /**< buffer to update whether an expression's bounds were propagated to an empty interval */
2614  int* ntightenings /**< buffer to store the number of (variable) tightenings */
2615  )
2616 {
2617  SCIP_CONSHDLRDATA* conshdlrdata;
2618  SCIP_EXPR* expr;
2619  SCIP_EXPR_OWNERDATA* ownerdata;
2620 
2621  assert(infeasible != NULL);
2622  assert(ntightenings != NULL);
2623 
2624  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2625  assert(conshdlrdata != NULL);
2626 
2627  *ntightenings = 0;
2628 
2629  /* main loop that calls reverse propagation for expressions on the queue
2630  * when reverseprop finds a tightening for an expression, then that expression is added to the queue (within the reverseprop call)
2631  */
2632  while( !SCIPqueueIsEmpty(conshdlrdata->reversepropqueue) && !(*infeasible) )
2633  {
2634  SCIP_INTERVAL propbounds;
2635  int e;
2636 
2637  expr = (SCIP_EXPR*) SCIPqueueRemove(conshdlrdata->reversepropqueue);
2638  assert(expr != NULL);
2639 
2640  ownerdata = SCIPexprGetOwnerData(expr);
2641  assert(ownerdata != NULL);
2642 
2643  assert(ownerdata->inpropqueue);
2644  /* mark that the expression is not in the queue anymore */
2645  ownerdata->inpropqueue = FALSE;
2646 
2647  /* since the expr was in the propagation queue, the propbounds should belong to current propagation and should not be empty
2648  * (propbounds being entire doesn't make much sense, so assert this for now, too, but that could be removed)
2649  */
2650  assert(ownerdata->propboundstag == conshdlrdata->curpropboundstag);
2651  assert(!SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, ownerdata->propbounds));
2652  assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, ownerdata->propbounds));
2653 
2654  /* this intersects propbounds with activity and auxvar bounds
2655  * I doubt this would be much helpful, since propbounds are already subset of activity and we also propagate
2656  * auxvar bounds separately, so disabling this for now
2657  */
2658 #ifdef SCIP_DISABLED_CODE
2659  propbounds = SCIPgetExprBoundsNonlinear(scip, expr);
2660  if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, propbounds) )
2661  {
2662  *infeasible = TRUE;
2663  break;
2664  }
2665 #else
2666  propbounds = ownerdata->propbounds;
2667 #endif
2668 
2669  if( ownerdata->nenfos > 0 )
2670  {
2671  /* for nodes with enforcement, call reverse propagation callbacks of nlhdlrs */
2672  for( e = 0; e < ownerdata->nenfos && !*infeasible; ++e )
2673  {
2674  SCIP_NLHDLR* nlhdlr;
2675  int nreds;
2676 
2677  /* skip nlhdlr if it does not want to participate in activity computation */
2678  if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY) == 0 )
2679  continue;
2680 
2681  nlhdlr = ownerdata->enfos[e]->nlhdlr;
2682  assert(nlhdlr != NULL);
2683 
2684  /* call the reverseprop of the nlhdlr */
2685 #ifdef SCIP_DEBUG
2686  SCIPdebugMsg(scip, "call reverse propagation for ");
2687  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2688  SCIPdebugMsgPrint(scip, " in [%g,%g] using nlhdlr <%s>\n", propbounds.inf, propbounds.sup, SCIPnlhdlrGetName(nlhdlr));
2689 #endif
2690 
2691  nreds = 0;
2692  SCIP_CALL( SCIPnlhdlrReverseprop(scip, conshdlr, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata, propbounds, infeasible, &nreds) );
2693  assert(nreds >= 0);
2694  *ntightenings += nreds;
2695  }
2696  }
2698  {
2699  /* if expr without enforcement (before detect), call reverse propagation callback of exprhdlr directly */
2700  SCIP_INTERVAL* childrenbounds;
2701  int c;
2702 
2703 #ifdef SCIP_DEBUG
2704  SCIPdebugMsg(scip, "call reverse propagation for ");
2705  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2706  SCIPdebugMsgPrint(scip, " in [%g,%g] using exprhdlr <%s>\n", SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup, SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)));
2707 #endif
2708 
2709  /* if someone added an expr without nlhdlr into the reversepropqueue, then this must be because its enfo hasn't
2710  * been initialized in detectNlhdlr yet (nenfos < 0)
2711  */
2712  assert(ownerdata->nenfos < 0);
2713 
2714  SCIP_CALL( SCIPallocBufferArray(scip, &childrenbounds, SCIPexprGetNChildren(expr)) );
2715  for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
2716  childrenbounds[c] = SCIPgetExprBoundsNonlinear(scip, SCIPexprGetChildren(expr)[c]);
2717 
2718  /* call the reverseprop of the exprhdlr */
2719  SCIP_CALL( SCIPcallExprReverseprop(scip, expr, propbounds, childrenbounds, infeasible) );
2720 
2721  if( !*infeasible )
2722  for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
2723  {
2724  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[c], childrenbounds[c], infeasible, ntightenings) );
2725  }
2726 
2727  SCIPfreeBufferArray(scip, &childrenbounds);
2728  }
2729  }
2730 
2731  /* reset inpropqueue for all remaining expr's in queue (can happen in case of early stop due to infeasibility) */
2732  while( !SCIPqueueIsEmpty(conshdlrdata->reversepropqueue) )
2733  {
2734  expr = (SCIP_EXPR*) SCIPqueueRemove(conshdlrdata->reversepropqueue);
2735  assert(expr != NULL);
2736 
2737  ownerdata = SCIPexprGetOwnerData(expr);
2738  assert(ownerdata != NULL);
2739 
2740  /* mark that the expression is not in the queue anymore */
2741  ownerdata->inpropqueue = FALSE;
2742  }
2743 
2744  return SCIP_OKAY;
2745 }
2746 
2747 /** calls domain propagation for a given set of constraints
2748  *
2749  * The algorithm alternates calls of forward and reverse propagation.
2750  * Forward propagation ensures that activity of expressions is up to date.
2751  * Reverse propagation tries to derive tighter variable bounds by reversing the activity computation, using the constraints
2752  * [lhs,rhs] interval as starting point.
2753  *
2754  * The propagation algorithm works as follows:
2755  * 1. apply forward propagation (update activities) for all constraints not marked as propagated
2756  * 2. if presolve or propauxvars is disabled: collect expressions for which the constraint sides provide tighter bounds
2757  * if solve and propauxvars is enabled: collect expressions for which auxvars (including those in root exprs)
2758  * provide tighter bounds
2759  * 3. apply reverse propagation to all collected expressions; don't explore
2760  * sub-expressions which have not changed since the beginning of the propagation loop
2761  * 4. if we have found enough tightenings go to 1, otherwise leave propagation loop
2762  *
2763  * @note After calling forward propagation for a constraint, we mark this constraint as propagated. This flag might be
2764  * reset during the reverse propagation when we find a bound tightening of a variable expression contained in the
2765  * constraint. Resetting this flag is done in the EVENTEXEC callback of the event handler
2766  *
2767  * TODO should we distinguish between expressions where activity information is used for separation and those where not,
2768  * e.g., try less to propagate on convex constraints?
2769  */
2770 static
2772  SCIP* scip, /**< SCIP data structure */
2773  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2774  SCIP_CONS** conss, /**< constraints to propagate */
2775  int nconss, /**< total number of constraints */
2776  SCIP_Bool force, /**< force tightening even if below bound strengthening tolerance */
2777  SCIP_RESULT* result, /**< pointer to store the result */
2778  int* nchgbds /**< buffer to add the number of changed bounds */
2779  )
2780 {
2781  SCIP_CONSHDLRDATA* conshdlrdata;
2782  SCIP_CONSDATA* consdata;
2783  SCIP_EXPR_OWNERDATA* ownerdata;
2784  SCIP_Bool cutoff = FALSE;
2785  SCIP_INTERVAL conssides;
2786  int ntightenings;
2787  int roundnr;
2788  SCIP_EXPRITER* revpropcollectit = NULL;
2789  int i;
2790 
2791  assert(scip != NULL);
2792  assert(conshdlr != NULL);
2793  assert(conss != NULL);
2794  assert(nconss >= 0);
2795  assert(result != NULL);
2796  assert(nchgbds != NULL);
2797  assert(*nchgbds >= 0);
2798 
2799  /* no constraints to propagate */
2800  if( nconss == 0 )
2801  {
2802  *result = SCIP_DIDNOTRUN;
2803  return SCIP_OKAY;
2804  }
2805 
2806  conshdlrdata = SCIPconshdlrGetData(conshdlr);
2807  assert(conshdlrdata != NULL);
2808  assert(conshdlrdata->intevalvar == intEvalVarBoundTightening);
2809  assert(!conshdlrdata->globalbounds);
2810 
2811  *result = SCIP_DIDNOTFIND;
2812  roundnr = 0;
2813 
2814  /* tightenAuxVarBounds() needs to know whether boundtightenings are to be forced */
2815  conshdlrdata->forceboundtightening = force;
2816 
2817  /* invalidate all propbounds (probably not needed) */
2818  ++conshdlrdata->curpropboundstag;
2819 
2820  /* create iterator that we will use if we need to look at all auxvars */
2821  if( conshdlrdata->propauxvars )
2822  {
2823  SCIP_CALL( SCIPcreateExpriter(scip, &revpropcollectit) );
2824  }
2825 
2826  /* main propagation loop */
2827  do
2828  {
2829  SCIPdebugMsg(scip, "start propagation round %d\n", roundnr);
2830 
2831  assert(SCIPqueueIsEmpty(conshdlrdata->reversepropqueue));
2832 
2833  /* apply forward propagation (update expression activities)
2834  * and add promising root expressions into queue for reversepropagation
2835  */
2836  for( i = 0; i < nconss; ++i )
2837  {
2838  consdata = SCIPconsGetData(conss[i]);
2839  assert(consdata != NULL);
2840 
2841  /* skip deleted, non-active, or propagation-disabled constraints */
2842  if( SCIPconsIsDeleted(conss[i]) || !SCIPconsIsActive(conss[i]) || !SCIPconsIsPropagationEnabled(conss[i]) )
2843  continue;
2844 
2845  /* skip already propagated constraints, i.e., constraints where no (original) variable has changed and thus
2846  * activity didn't change
2847  */
2848  if( consdata->ispropagated )
2849  continue;
2850 
2851  /* update activities in expression */
2852  SCIPdebugMsg(scip, "call forwardPropExpr() for constraint <%s> (round %d): ", SCIPconsGetName(conss[i]), roundnr);
2853  SCIPdebugPrintCons(scip, conss[i], NULL);
2854 
2855  ntightenings = 0;
2856  SCIP_CALL( forwardPropExpr(scip, conshdlr, consdata->expr, TRUE, &cutoff, &ntightenings) );
2857  assert(cutoff || !SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(consdata->expr)));
2858 
2859  if( cutoff )
2860  {
2861  SCIPdebugMsg(scip, " -> cutoff in forwardPropExpr (due to domain error or auxvar tightening) of constraint <%s>\n", SCIPconsGetName(conss[i]));
2862  *result = SCIP_CUTOFF;
2863  break;
2864  }
2865 
2866  ownerdata = SCIPexprGetOwnerData(consdata->expr);
2867 
2868  /* TODO for a constraint that only has an auxvar for consdata->expr (e.g., convex quadratic), we could also just do the if(TRUE)-branch */
2869  if( !conshdlrdata->propauxvars || ownerdata->auxvar == NULL )
2870  {
2871  /* check whether constraint sides (relaxed by epsilon) or auxvar bounds provide a tightening
2872  * (if we have auxvar (not in presolve), then bounds of the auxvar are initially set to constraint sides,
2873  * so taking auxvar bounds is enough)
2874  */
2875  if( ownerdata->auxvar == NULL )
2876  {
2877  /* relax sides by SCIPepsilon() and handle infinite sides */
2878  SCIP_Real lhs = SCIPisInfinity(scip, -consdata->lhs) ? -SCIP_INTERVAL_INFINITY : consdata->lhs - conshdlrdata->conssiderelaxamount;
2879  SCIP_Real rhs = SCIPisInfinity(scip, consdata->rhs) ? SCIP_INTERVAL_INFINITY : consdata->rhs + conshdlrdata->conssiderelaxamount;
2880  SCIPintervalSetBounds(&conssides, lhs, rhs);
2881  }
2882  else
2883  {
2884  conssides = intEvalVarBoundTightening(scip, ownerdata->auxvar, (void*)conshdlrdata);
2885  }
2886  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, consdata->expr, conssides, &cutoff, &ntightenings) );
2887  }
2888  else
2889  {
2890  /* check whether bounds of any auxvar used in constraint provides a tightening
2891  * (for the root expression, bounds of auxvar are initially set to constraint sides)
2892  * but skip exprs that have an auxvar, but do not participate in propagation
2893  */
2894  SCIP_EXPR* expr;
2895 
2896  assert(revpropcollectit != NULL);
2897  SCIP_CALL( SCIPexpriterInit(revpropcollectit, consdata->expr, SCIP_EXPRITER_BFS, FALSE) );
2898  for( expr = SCIPexpriterGetCurrent(revpropcollectit); !SCIPexpriterIsEnd(revpropcollectit) && !cutoff; expr = SCIPexpriterGetNext(revpropcollectit) )
2899  {
2900  ownerdata = SCIPexprGetOwnerData(expr);
2901  assert(ownerdata != NULL);
2902 
2903  if( ownerdata->auxvar == NULL )
2904  continue;
2905 
2906  if( ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 )
2907  continue;
2908 
2909  conssides = intEvalVarBoundTightening(scip, ownerdata->auxvar, (void*)conshdlrdata);
2910  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, expr, conssides, &cutoff, &ntightenings) );
2911  }
2912  }
2913 
2914  if( cutoff )
2915  {
2916  SCIPdebugMsg(scip, " -> cutoff after intersect with conssides of constraint <%s>\n", SCIPconsGetName(conss[i]));
2917  *result = SCIP_CUTOFF;
2918  break;
2919  }
2920 
2921  assert(ntightenings >= 0);
2922  if( ntightenings > 0 )
2923  {
2924  *nchgbds += ntightenings;
2925  *result = SCIP_REDUCEDDOM;
2926  }
2927 
2928  /* mark constraint as propagated; this will be reset via the event system when we find a variable tightening */
2929  consdata->ispropagated = TRUE;
2930  }
2931 
2932  /* apply backward propagation (if cutoff is TRUE, then this call empties the queue) */
2933  SCIP_CALL( reversePropQueue(scip, conshdlr, &cutoff, &ntightenings) );
2934  assert(ntightenings >= 0);
2935  assert(SCIPqueueIsEmpty(conshdlrdata->reversepropqueue));
2936 
2937  if( cutoff )
2938  {
2939  SCIPdebugMsg(scip, " -> cutoff\n");
2940  *result = SCIP_CUTOFF;
2941  break;
2942  }
2943 
2944  if( ntightenings > 0 )
2945  {
2946  *nchgbds += ntightenings;
2947  *result = SCIP_REDUCEDDOM;
2948  }
2949  }
2950  while( ntightenings > 0 && ++roundnr < conshdlrdata->maxproprounds );
2951 
2952  if( conshdlrdata->propauxvars )
2953  {
2954  SCIPfreeExpriter(&revpropcollectit);
2955  }
2956 
2957  conshdlrdata->forceboundtightening = FALSE;
2958 
2959  /* invalidate propbounds in all exprs, so noone accidentally uses them outside propagation */
2960  ++conshdlrdata->curpropboundstag;
2961 
2962  return SCIP_OKAY;
2963 }
2964 
2965 /** calls the reverseprop callbacks of all nlhdlrs in all expressions in all constraints using activity as bounds
2966  *
2967  * This is meant to propagate any domain restrictions on functions onto variable bounds, if possible.
2968  *
2969  * Assumes that activities are still valid and curpropboundstag does not need to be increased.
2970  * Therefore, a good place to call this function is immediately after propConss() or after forwardPropExpr() if outside propagation.
2971  */
2972 static
2974  SCIP* scip, /**< SCIP data structure */
2975  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2976  SCIP_CONS** conss, /**< constraints to propagate */
2977  int nconss, /**< total number of constraints */
2978  SCIP_RESULT* result, /**< pointer to store the result */
2979  int* nchgbds /**< buffer to add the number of changed bounds */
2980  )
2981 {
2982  SCIP_CONSDATA* consdata;
2983  SCIP_EXPRITER* it;
2984  SCIP_EXPR* expr;
2985  SCIP_EXPR_OWNERDATA* ownerdata;
2986  SCIP_Bool cutoff = FALSE;
2987  int ntightenings;
2988  int c;
2989  int e;
2990 
2991  assert(scip != NULL);
2992  assert(conshdlr != NULL);
2993  assert(conss != NULL);
2994  assert(nconss >= 0);
2995  assert(result != NULL);
2996  assert(nchgbds != NULL);
2997  assert(*nchgbds >= 0);
2998 
2999  assert(SCIPconshdlrGetData(conshdlr)->intevalvar == intEvalVarBoundTightening);
3000  assert(!SCIPconshdlrGetData(conshdlr)->globalbounds);
3001  assert(SCIPqueueIsEmpty(SCIPconshdlrGetData(conshdlr)->reversepropqueue));
3002 
3003  *result = SCIP_DIDNOTFIND;
3004 
3005  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
3007 
3008  for( c = 0; c < nconss && !cutoff; ++c )
3009  {
3010  /* skip deleted, non-active, or propagation-disabled constraints */
3011  if( SCIPconsIsDeleted(conss[c]) || !SCIPconsIsActive(conss[c]) || !SCIPconsIsPropagationEnabled(conss[c]) )
3012  continue;
3013 
3014  consdata = SCIPconsGetData(conss[c]);
3015  assert(consdata != NULL);
3016 
3017  for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it) && !cutoff; expr = SCIPexpriterGetNext(it) )
3018  {
3019  ownerdata = SCIPexprGetOwnerData(expr);
3020  assert(ownerdata != NULL);
3021 
3022  /* call reverseprop for those nlhdlr that participate in this expr's activity computation
3023  * this will propagate the current activity
3024  */
3025  for( e = 0; e < ownerdata->nenfos; ++e )
3026  {
3027  SCIP_NLHDLR* nlhdlr;
3028  assert(ownerdata->enfos[e] != NULL);
3029 
3030  nlhdlr = ownerdata->enfos[e]->nlhdlr;
3031  assert(nlhdlr != NULL);
3032  if( (ownerdata->enfos[e]->nlhdlrparticipation & SCIP_NLHDLR_METHOD_ACTIVITY) == 0 )
3033  continue;
3034 
3035  SCIPdebugMsg(scip, "propExprDomains calling reverseprop for expression %p [%g,%g]\n", (void*)expr,
3036  SCIPexprGetActivity(expr).inf, SCIPexprGetActivity(expr).sup);
3037  ntightenings = 0;
3038  SCIP_CALL( SCIPnlhdlrReverseprop(scip, conshdlr, nlhdlr, expr, ownerdata->enfos[e]->nlhdlrexprdata,
3039  SCIPexprGetActivity(expr), &cutoff, &ntightenings) );
3040 
3041  if( cutoff )
3042  {
3043  /* stop everything if we detected infeasibility */
3044  SCIPdebugMsg(scip, "detect infeasibility for constraint <%s> during reverseprop()\n", SCIPconsGetName(conss[c]));
3045  *result = SCIP_CUTOFF;
3046  break;
3047  }
3048 
3049  assert(ntightenings >= 0);
3050  if( ntightenings > 0 )
3051  {
3052  *nchgbds += ntightenings;
3053  *result = SCIP_REDUCEDDOM;
3054  }
3055  }
3056  }
3057  }
3058 
3059  /* apply backward propagation (if cutoff is TRUE, then this call empties the queue) */
3060  SCIP_CALL( reversePropQueue(scip, conshdlr, &cutoff, &ntightenings) );
3061  assert(ntightenings >= 0);
3062 
3063  if( cutoff )
3064  {
3065  SCIPdebugMsg(scip, " -> cutoff\n");
3066  *result = SCIP_CUTOFF;
3067  }
3068  else if( ntightenings > 0 )
3069  {
3070  *nchgbds += ntightenings;
3071  *result = SCIP_REDUCEDDOM;
3072  }
3073 
3074  SCIPfreeExpriter(&it);
3075 
3076  /* invalidate propbounds in all exprs, so noone accidentally uses them outside propagation */
3077  ++SCIPconshdlrGetData(conshdlr)->curpropboundstag;
3078 
3079  return SCIP_OKAY;
3080 }
3081 
3082 /** propagates variable locks through expression and adds locks to variables */
3083 static
3085  SCIP* scip, /**< SCIP data structure */
3086  SCIP_EXPR* expr, /**< expression */
3087  int nlockspos, /**< number of positive locks */
3088  int nlocksneg /**< number of negative locks */
3089  )
3090 {
3091  SCIP_EXPR_OWNERDATA* ownerdata;
3092  SCIP_EXPRITER* it;
3093  SCIP_EXPRITER_USERDATA ituserdata;
3094 
3095  assert(expr != NULL);
3096 
3097  /* if no locks, then nothing to propagate */
3098  if( nlockspos == 0 && nlocksneg == 0 )
3099  return SCIP_OKAY;
3100 
3101  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
3104  assert(SCIPexpriterGetCurrent(it) == expr); /* iterator should not have moved */
3105 
3106  /* store locks in root node */
3107  ituserdata.intvals[0] = nlockspos;
3108  ituserdata.intvals[1] = nlocksneg;
3109  SCIPexpriterSetCurrentUserData(it, ituserdata);
3110 
3111  while( !SCIPexpriterIsEnd(it) )
3112  {
3113  /* collect locks */
3114  ituserdata = SCIPexpriterGetCurrentUserData(it);
3115  nlockspos = ituserdata.intvals[0];
3116  nlocksneg = ituserdata.intvals[1];
3117 
3118  ownerdata = SCIPexprGetOwnerData(expr);
3119 
3120  switch( SCIPexpriterGetStageDFS(it) )
3121  {
3123  {
3124  if( SCIPisExprVar(scip, expr) )
3125  {
3126  /* if a variable, then also add nlocksneg/nlockspos via SCIPaddVarLocks() */
3127  SCIP_CALL( SCIPaddVarLocks(scip, SCIPgetVarExprVar(expr), nlocksneg, nlockspos) );
3128  }
3129 
3130  /* add locks to expression */
3131  ownerdata->nlockspos += nlockspos;
3132  ownerdata->nlocksneg += nlocksneg;
3133 
3134  /* add monotonicity information if expression has been locked for the first time */
3135  if( ownerdata->nlockspos == nlockspos && ownerdata->nlocksneg == nlocksneg && SCIPexprGetNChildren(expr) > 0
3137  {
3138  int i;
3139 
3140  assert(ownerdata->monotonicity == NULL);
3141  assert(ownerdata->monotonicitysize == 0);
3142 
3143  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &ownerdata->monotonicity, SCIPexprGetNChildren(expr)) );
3144  ownerdata->monotonicitysize = SCIPexprGetNChildren(expr);
3145 
3146  /* store the monotonicity for each child */
3147  for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
3148  {
3149  SCIP_CALL( SCIPcallExprMonotonicity(scip, expr, i, &ownerdata->monotonicity[i]) );
3150  }
3151  }
3152  break;
3153  }
3154 
3156  {
3157  /* remove monotonicity information if expression has been unlocked */
3158  if( ownerdata->nlockspos == 0 && ownerdata->nlocksneg == 0 && ownerdata->monotonicity != NULL )
3159  {
3160  assert(ownerdata->monotonicitysize > 0);
3161  /* keep this assert for checking whether someone changed an expression without updating locks properly */
3162  assert(ownerdata->monotonicitysize == SCIPexprGetNChildren(expr));
3163 
3164  SCIPfreeBlockMemoryArray(scip, &ownerdata->monotonicity, ownerdata->monotonicitysize);
3165  ownerdata->monotonicitysize = 0;
3166  }
3167  break;
3168  }
3169 
3171  {
3172  SCIP_MONOTONE monotonicity;
3173 
3174  /* get monotonicity of child */
3175  /* NOTE: the monotonicity stored in an expression might be different from the result obtained by
3176  * SCIPcallExprMonotonicity
3177  */
3178  monotonicity = ownerdata->monotonicity != NULL ? ownerdata->monotonicity[SCIPexpriterGetChildIdxDFS(it)] : SCIP_MONOTONE_UNKNOWN;
3179 
3180  /* compute resulting locks of the child expression */
3181  switch( monotonicity )
3182  {
3183  case SCIP_MONOTONE_INC:
3184  ituserdata.intvals[0] = nlockspos;
3185  ituserdata.intvals[1] = nlocksneg;
3186  break;
3187  case SCIP_MONOTONE_DEC:
3188  ituserdata.intvals[0] = nlocksneg;
3189  ituserdata.intvals[1] = nlockspos;
3190  break;
3191  case SCIP_MONOTONE_UNKNOWN:
3192  ituserdata.intvals[0] = nlockspos + nlocksneg;
3193  ituserdata.intvals[1] = nlockspos + nlocksneg;
3194  break;
3195  case SCIP_MONOTONE_CONST:
3196  ituserdata.intvals[0] = 0;
3197  ituserdata.intvals[1] = 0;
3198  break;
3199  }
3200  /* set locks in child expression */
3201  SCIPexpriterSetChildUserData(it, ituserdata);
3202 
3203  break;
3204  }
3205 
3206  default :
3207  /* you should never be here */
3208  SCIPABORT();
3209  break;
3210  }
3211 
3212  expr = SCIPexpriterGetNext(it);
3213  }
3214 
3215  SCIPfreeExpriter(&it);
3216 
3217  return SCIP_OKAY;
3218 }
3219 
3220 /** main function for adding locks to expressions and variables
3221  *
3222  * Locks for a nonlinear constraint are used to update locks for all sub-expressions and variables.
3223  * Locks of expressions depend on the monotonicity of expressions w.r.t. their children, e.g.,
3224  * consider the constraint \f$x^2 \leq 1\f$ with \f$x \in [-2,-1]\f$ implies an up-lock for the root
3225  * expression (pow) and a down-lock for its child \f$x\f$ because \f$x^2\f$ is decreasing on [-2,-1].
3226  * Since the monotonicity (and thus the locks) might also depend on variable bounds, the function remembers
3227  * the computed monotonicity information of each expression until all locks of an expression have been removed,
3228  * which implies that updating the monotonicity information during the next locking of this expression does not
3229  * break existing locks.
3230  *
3231  * @note When modifying the structure of an expression, e.g., during simplification, it is necessary to remove all
3232  * locks from an expression and repropagating them after the structural changes have been applied.
3233  * Because of existing common sub-expressions, it might be necessary to remove the locks of all constraints
3234  * to ensure that an expression is unlocked (see canonicalizeConstraints() for an example)
3235  */
3236 static
3238  SCIP* scip, /**< SCIP data structure */
3239  SCIP_CONS* cons, /**< nonlinear constraint */
3240  int nlockspos, /**< number of positive rounding locks */
3241  int nlocksneg /**< number of negative rounding locks */
3242  )
3243 {
3244  SCIP_CONSDATA* consdata;
3245 
3246  assert(cons != NULL);
3247 
3248  if( nlockspos == 0 && nlocksneg == 0 )
3249  return SCIP_OKAY;
3250 
3251  consdata = SCIPconsGetData(cons);
3252  assert(consdata != NULL);
3253 
3254  /* no constraint sides -> nothing to lock */
3255  if( SCIPisInfinity(scip, consdata->rhs) && SCIPisInfinity(scip, -consdata->lhs) )
3256  return SCIP_OKAY;
3257 
3258  /* remember locks */
3259  consdata->nlockspos += nlockspos;
3260  consdata->nlocksneg += nlocksneg;
3261 
3262  assert(consdata->nlockspos >= 0);
3263  assert(consdata->nlocksneg >= 0);
3264 
3265  /* compute locks for lock propagation */
3266  if( !SCIPisInfinity(scip, consdata->rhs) && !SCIPisInfinity(scip, -consdata->lhs) )
3267  {
3268  SCIP_CALL( propagateLocks(scip, consdata->expr, nlockspos + nlocksneg, nlockspos + nlocksneg));
3269  }
3270  else if( !SCIPisInfinity(scip, consdata->rhs) )
3271  {
3272  SCIP_CALL( propagateLocks(scip, consdata->expr, nlockspos, nlocksneg));
3273  }
3274  else
3275  {
3276  assert(!SCIPisInfinity(scip, -consdata->lhs));
3277  SCIP_CALL( propagateLocks(scip, consdata->expr, nlocksneg, nlockspos));
3278  }
3279 
3280  return SCIP_OKAY;
3281 }
3282 
3283 /** create a nonlinear row representation of a nonlinear constraint and stores them in consdata */
3284 static
3286  SCIP* scip, /**< SCIP data structure */
3287  SCIP_CONS* cons /**< nonlinear constraint */
3288  )
3289 {
3290  SCIP_CONSDATA* consdata;
3291 
3292  assert(scip != NULL);
3293  assert(cons != NULL);
3294 
3295  consdata = SCIPconsGetData(cons);
3296  assert(consdata != NULL);
3297  assert(consdata->expr != NULL);
3298 
3299  if( consdata->nlrow != NULL )
3300  {
3301  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
3302  }
3303 
3304  /* better curvature info will be set in initSolve() just before nlrow is added to NLP */
3305  SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
3306  0, NULL, NULL, NULL, consdata->lhs, consdata->rhs, SCIP_EXPRCURV_UNKNOWN) );
3307 
3308  if( SCIPisExprSum(scip, consdata->expr) )
3309  {
3310  /* if root is a sum, then split into linear and nonlinear terms */
3311  SCIP_EXPR* nonlinpart;
3312  SCIP_EXPR* child;
3313  SCIP_Real* coefs;
3314  int i;
3315 
3316  coefs = SCIPgetCoefsExprSum(consdata->expr);
3317 
3318  /* constant term of sum */
3319  SCIP_CALL( SCIPchgNlRowConstant(scip, consdata->nlrow, SCIPgetConstantExprSum(consdata->expr)) );
3320 
3321  /* a sum-expression that will hold the nonlinear terms and be passed to the nlrow eventually */
3322  SCIP_CALL( SCIPcreateExprSum(scip, &nonlinpart, 0, NULL, NULL, 0.0, exprownerCreate, (void*)SCIPconsGetHdlr(cons)) );
3323 
3324  for( i = 0; i < SCIPexprGetNChildren(consdata->expr); ++i )
3325  {
3326  child = SCIPexprGetChildren(consdata->expr)[i];
3327  if( SCIPisExprVar(scip, child) )
3328  {
3329  /* linear term */
3330  SCIP_CALL( SCIPaddLinearCoefToNlRow(scip, consdata->nlrow, SCIPgetVarExprVar(child), coefs[i]) );
3331  }
3332  else
3333  {
3334  /* nonlinear term */
3335  SCIP_CALL( SCIPappendExprSumExpr(scip, nonlinpart, child, coefs[i]) );
3336  }
3337  }
3338 
3339  if( SCIPexprGetNChildren(nonlinpart) > 0 )
3340  {
3341  /* add expression to nlrow (this will make a copy) */
3342  SCIP_CALL( SCIPsetNlRowExpr(scip, consdata->nlrow, nonlinpart) );
3343  }
3344  SCIP_CALL( SCIPreleaseExpr(scip, &nonlinpart) );
3345  }
3346  else
3347  {
3348  SCIP_CALL( SCIPsetNlRowExpr(scip, consdata->nlrow, consdata->expr) );
3349  }
3350 
3351  return SCIP_OKAY;
3352 }
3353 
3354 /** compares enfodata by enforcement priority of nonlinear handler
3355  *
3356  * If handlers have same enforcement priority, then compare by detection priority, then by name.
3357  */
3358 static
3359 SCIP_DECL_SORTPTRCOMP(enfodataCmp)
3361  SCIP_NLHDLR* h1;
3362  SCIP_NLHDLR* h2;
3363 
3364  assert(elem1 != NULL);
3365  assert(elem2 != NULL);
3366 
3367  h1 = ((EXPRENFO*)elem1)->nlhdlr;
3368  h2 = ((EXPRENFO*)elem2)->nlhdlr;
3369 
3370  assert(h1 != NULL);
3371  assert(h2 != NULL);
3372 
3375 
3378 
3379  return strcmp(SCIPnlhdlrGetName(h1), SCIPnlhdlrGetName(h2));
3380 }
3381 
3382 /** install nlhdlrs in one expression */
3383 static
3385  SCIP* scip, /**< SCIP data structure */
3386  SCIP_EXPR* expr, /**< expression for which to run detection routines */
3387  SCIP_CONS* cons /**< constraint for which expr == consdata->expr, otherwise NULL */
3388  )
3389 {
3390  SCIP_EXPR_OWNERDATA* ownerdata;
3391  SCIP_CONSHDLRDATA* conshdlrdata;
3392  SCIP_NLHDLR_METHOD enforcemethodsallowed;
3393  SCIP_NLHDLR_METHOD enforcemethods;
3394  SCIP_NLHDLR_METHOD enforcemethodsnew;
3395  SCIP_NLHDLR_METHOD nlhdlrenforcemethods;
3396  SCIP_NLHDLR_METHOD nlhdlrparticipating;
3397  SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
3398  int enfossize; /* allocated length of expr->enfos array */
3399  int h;
3400 
3401  assert(expr != NULL);
3402 
3403  ownerdata = SCIPexprGetOwnerData(expr);
3404  assert(ownerdata != NULL);
3405 
3406  conshdlrdata = SCIPconshdlrGetData(ownerdata->conshdlr);
3407  assert(conshdlrdata != NULL);
3408  assert(conshdlrdata->auxvarid >= 0);
3409  assert(!conshdlrdata->indetect);
3410 
3411  /* there should be no enforcer yet and detection should not even have considered expr yet */
3412  assert(ownerdata->nenfos < 0);
3413  assert(ownerdata->enfos == NULL);
3414 
3415  /* check which enforcement methods are required by setting flags in enforcemethods for those that are NOT required
3416  * - if no auxiliary variable is used, then do not need sepabelow or sepaabove
3417  * - if auxiliary variable is used, but nobody positively (up) locks expr -> only need to enforce expr >= auxvar -> no need for underestimation
3418  * - if auxiliary variable is used, but nobody negatively (down) locks expr -> only need to enforce expr <= auxvar -> no need for overestimation
3419  * - if no one uses activity, then do not need activity methods
3420  */
3421  enforcemethods = SCIP_NLHDLR_METHOD_NONE;
3422  if( ownerdata->nauxvaruses == 0 )
3423  enforcemethods |= SCIP_NLHDLR_METHOD_SEPABOTH;
3424  else
3425  {
3426  if( ownerdata->nlockspos == 0 ) /* no need for underestimation */
3427  enforcemethods |= SCIP_NLHDLR_METHOD_SEPABELOW;
3428  if( ownerdata->nlocksneg == 0 ) /* no need for overestimation */
3429  enforcemethods |= SCIP_NLHDLR_METHOD_SEPAABOVE;
3430  }
3431  if( ownerdata->nactivityusesprop == 0 && ownerdata->nactivityusessepa == 0 )
3432  enforcemethods |= SCIP_NLHDLR_METHOD_ACTIVITY;
3433 
3434  /* it doesn't make sense to have been called on detectNlhdlr, if the expr isn't used for anything */
3435  assert(enforcemethods != SCIP_NLHDLR_METHOD_ALL);
3436 
3437  /* all methods that have not been flagged above are the ones that we want to be handled by nlhdlrs */
3438  enforcemethodsallowed = ~enforcemethods & SCIP_NLHDLR_METHOD_ALL;
3439 
3440  ownerdata->nenfos = 0;
3441  enfossize = 2;
3442  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &ownerdata->enfos, enfossize) );
3443  conshdlrdata->indetect = TRUE;
3444 
3445  SCIPdebugMsg(scip, "detecting nlhdlrs for %s expression %p (%s); requiring%s%s%s\n",
3446  cons != NULL ? "root" : "non-root", (void*)expr, SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)),
3447  (enforcemethods & SCIP_NLHDLR_METHOD_SEPABELOW) != 0 ? "" : " sepabelow",
3448  (enforcemethods & SCIP_NLHDLR_METHOD_SEPAABOVE) != 0 ? "" : " sepaabove",
3449  (enforcemethods & SCIP_NLHDLR_METHOD_ACTIVITY) != 0 ? "" : " activity");
3450 
3451  for( h = 0; h < conshdlrdata->nnlhdlrs; ++h )
3452  {
3453  SCIP_NLHDLR* nlhdlr;
3454 
3455  nlhdlr = conshdlrdata->nlhdlrs[h];
3456  assert(nlhdlr != NULL);
3457 
3458  /* skip disabled nlhdlrs */
3459  if( !SCIPnlhdlrIsEnabled(nlhdlr) )
3460  continue;
3461 
3462  /* call detect routine of nlhdlr */
3463  nlhdlrexprdata = NULL;
3464  enforcemethodsnew = enforcemethods;
3465  nlhdlrparticipating = SCIP_NLHDLR_METHOD_NONE;
3466  conshdlrdata->registerusesactivitysepabelow = FALSE; /* SCIPregisterExprUsageNonlinear() as called by detect may set this to TRUE */
3467  conshdlrdata->registerusesactivitysepaabove = FALSE; /* SCIPregisterExprUsageNonlinear() as called by detect may set this to TRUE */
3468  SCIP_CALL( SCIPnlhdlrDetect(scip, ownerdata->conshdlr, nlhdlr, expr, cons, &enforcemethodsnew, &nlhdlrparticipating, &nlhdlrexprdata) );
3469 
3470  /* nlhdlr might have claimed more than needed: clean up sepa flags */
3471  nlhdlrparticipating &= enforcemethodsallowed;
3472 
3473  /* detection is only allowed to augment to nlhdlrenforcemethods, so previous enforcemethods must still be set */
3474  assert((enforcemethodsnew & enforcemethods) == enforcemethods);
3475 
3476  /* Because of the previous assert, nlhdlrenforcenew ^ enforcemethods are the methods enforced by this nlhdlr.
3477  * They are also cleaned up here to ensure that only the needed methods are claimed.
3478  */
3479  nlhdlrenforcemethods = (enforcemethodsnew ^ enforcemethods) & enforcemethodsallowed;
3480 
3481  /* nlhdlr needs to participate for the methods it is enforcing */
3482  assert((nlhdlrparticipating & nlhdlrenforcemethods) == nlhdlrenforcemethods);
3483 
3484  if( nlhdlrparticipating == SCIP_NLHDLR_METHOD_NONE )
3485  {
3486  /* nlhdlr might not have detected anything, or all set flags might have been removed by
3487  * clean up; in the latter case, we may need to free nlhdlrexprdata */
3488 
3489  /* free nlhdlr exprdata, if there is any and there is a method to free this data */
3490  if( nlhdlrexprdata != NULL )
3491  {
3492  SCIP_CALL( SCIPnlhdlrFreeexprdata(scip, nlhdlr, expr, &nlhdlrexprdata) );
3493  }
3494  /* nlhdlr cannot have added an enforcement method if it doesn't participate (actually redundant due to previous asserts) */
3495  assert(nlhdlrenforcemethods == SCIP_NLHDLR_METHOD_NONE);
3496 
3497  SCIPdebugMsg(scip, "nlhdlr <%s> detect unsuccessful\n", SCIPnlhdlrGetName(nlhdlr));
3498 
3499  continue;
3500  }
3501 
3502  SCIPdebugMsg(scip, "nlhdlr <%s> detect successful; sepabelow: %s, sepaabove: %s, activity: %s\n",
3503  SCIPnlhdlrGetName(nlhdlr),
3504  ((nlhdlrenforcemethods & SCIP_NLHDLR_METHOD_SEPABELOW) != 0) ? "enforcing" : ((nlhdlrparticipating & SCIP_NLHDLR_METHOD_SEPABELOW) != 0) ? "participating" : "no",
3505  ((nlhdlrenforcemethods & SCIP_NLHDLR_METHOD_SEPAABOVE) != 0) ? "enforcing" : ((nlhdlrparticipating & SCIP_NLHDLR_METHOD_SEPAABOVE) != 0) ? "participating" : "no",
3506  ((nlhdlrenforcemethods & SCIP_NLHDLR_METHOD_ACTIVITY) != 0) ? "enforcing" : ((nlhdlrparticipating & SCIP_NLHDLR_METHOD_ACTIVITY) != 0) ? "participating" : "no");
3507 
3508  /* store nlhdlr and its data */
3509  SCIP_CALL( SCIPensureBlockMemoryArray(scip, &ownerdata->enfos, &enfossize, ownerdata->nenfos+1) );
3510  SCIP_CALL( SCIPallocBlockMemory(scip, &ownerdata->enfos[ownerdata->nenfos]) );
3511  ownerdata->enfos[ownerdata->nenfos]->nlhdlr = nlhdlr;
3512  ownerdata->enfos[ownerdata->nenfos]->nlhdlrexprdata = nlhdlrexprdata;
3513  ownerdata->enfos[ownerdata->nenfos]->nlhdlrparticipation = nlhdlrparticipating;
3514  ownerdata->enfos[ownerdata->nenfos]->issepainit = FALSE;
3515  ownerdata->enfos[ownerdata->nenfos]->sepabelowusesactivity = conshdlrdata->registerusesactivitysepabelow;
3516  ownerdata->enfos[ownerdata->nenfos]->sepaaboveusesactivity = conshdlrdata->registerusesactivitysepaabove;
3517  ownerdata->nenfos++;
3518 
3519  /* update enforcement flags */
3520  enforcemethods = enforcemethodsnew;
3521  }
3522 
3523  conshdlrdata->indetect = FALSE;
3524 
3525  /* stop if an enforcement method is missing but we are already in solving stage
3526  * (as long as the expression provides its callbacks, the default nlhdlr should have provided all enforcement methods)
3527  */
3528  if( enforcemethods != SCIP_NLHDLR_METHOD_ALL && SCIPgetStage(scip) == SCIP_STAGE_SOLVING )
3529  {
3530  SCIPerrorMessage("no nonlinear handler provided some of the required enforcement methods\n");
3531  return SCIP_ERROR;
3532  }
3533 
3534  assert(ownerdata->nenfos > 0);
3535 
3536  /* sort nonlinear handlers by enforcement priority, in decreasing order */
3537  if( ownerdata->nenfos > 1 )
3538  SCIPsortDownPtr((void**)ownerdata->enfos, enfodataCmp, ownerdata->nenfos);
3539 
3540  /* resize enfos array to be nenfos long */
3541  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &ownerdata->enfos, enfossize, ownerdata->nenfos) );
3542 
3543  return SCIP_OKAY;
3544 }
3545 
3546 /** detect nlhdlrs that can handle the expressions */
3547 static
3549  SCIP* scip, /**< SCIP data structure */
3550  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3551  SCIP_CONS** conss, /**< constraints for which to run nlhdlr detect */
3552  int nconss /**< total number of constraints */
3553  )
3554 {
3555  SCIP_CONSHDLRDATA* conshdlrdata;
3556  SCIP_CONSDATA* consdata;
3557  SCIP_EXPR* expr;
3558  SCIP_EXPR_OWNERDATA* ownerdata;
3559  SCIP_EXPRITER* it;
3560  int i;
3561 
3562  assert(conss != NULL || nconss == 0);
3563  assert(nconss >= 0);
3564  assert(SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING || SCIPgetStage(scip) == SCIP_STAGE_INITSOLVE || SCIPgetStage(scip) == SCIP_STAGE_SOLVING); /* should only be called in presolve or initsolve or consactive */
3565 
3566  conshdlrdata = SCIPconshdlrGetData(conshdlr);
3567  assert(conshdlrdata != NULL);
3568 
3569  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
3571 
3572  if( SCIPgetStage(scip) == SCIP_STAGE_SOLVING && SCIPgetDepth(scip) != 0 )
3573  {
3574  /* ensure that activities are recomputed w.r.t. the global variable bounds if CONSACTIVE is called in a local node;
3575  * for example, this happens if globally valid nonlinear constraints are added during the tree search
3576  */
3578  conshdlrdata->globalbounds = TRUE;
3579  conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
3580  }
3581 
3582  for( i = 0; i < nconss; ++i )
3583  {
3584  assert(conss != NULL && conss[i] != NULL);
3585 
3586  consdata = SCIPconsGetData(conss[i]);
3587  assert(consdata != NULL);
3588  assert(consdata->expr != NULL);
3589 
3590  /* if a constraint is separated, we currently need it to be initial, too
3591  * this is because INITLP will create the auxiliary variables that are used for any separation
3592  * TODO we may relax this with a little more programming effort when required, see also TODO in INITLP
3593  */
3594  assert((!SCIPconsIsSeparated(conss[i]) && !SCIPconsIsEnforced(conss[i])) || SCIPconsIsInitial(conss[i]));
3595 
3596  ownerdata = SCIPexprGetOwnerData(consdata->expr);
3597  assert(ownerdata != NULL);
3598 
3599  /* because of common sub-expressions it might happen that we already detected a nonlinear handler and added it to the expr
3600  * then we would normally skip to run DETECT again
3601  * HOWEVER: most likely we have been running DETECT with cons == NULL, which may interest less nlhdlrs
3602  * thus, if expr is the root expression, we rerun DETECT
3603  */
3604  if( ownerdata->nenfos > 0 )
3605  {
3606  SCIP_CALL( freeEnfoData(scip, consdata->expr, FALSE) );
3607  assert(ownerdata->nenfos < 0);
3608  }
3609 
3610  /* if constraint will be enforced, and we are in solve, then ensure auxiliary variable for root expression
3611  * this way we can treat the root expression like any other expression when enforcing via separation
3612  * if constraint will be propagated, then register activity usage of root expression
3613  * this can trigger a call to forwardPropExpr, for which we better have the indetect flag set
3614  */
3615  conshdlrdata->indetect = TRUE;
3616  SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, consdata->expr,
3617  SCIPgetStage(scip) >= SCIP_STAGE_INITSOLVE && (SCIPconsIsSeparated(conss[i]) || SCIPconsIsEnforced(conss[i])),
3618  SCIPconsIsPropagated(conss[i]),
3619  FALSE, FALSE) );
3620  conshdlrdata->indetect = FALSE;
3621 
3622  /* compute integrality information for all subexpressions */
3623  SCIP_CALL( SCIPcomputeExprIntegrality(scip, consdata->expr) );
3624 
3625  /* run detectNlhdlr on all expr where required */
3626  for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
3627  {
3628  ownerdata = SCIPexprGetOwnerData(expr);
3629  assert(ownerdata != NULL);
3630 
3631  /* skip exprs that we already looked at */
3632  if( ownerdata->nenfos >= 0 )
3633  continue;
3634 
3635  /* if there is use of the auxvar, then someone requires that
3636  * auxvar == expr (or auxvar >= expr or auxvar <= expr) or we are at the root expression (expr==consdata->expr)
3637  * thus, we need to find nlhdlrs that separate or estimate
3638  * if there is use of the activity, then there is someone requiring that
3639  * activity of this expression is updated; this someone would also benefit from better bounds on the activity of this expression
3640  * thus, we need to find nlhdlrs that do interval-evaluation
3641  */
3642  if( ownerdata->nauxvaruses > 0 || ownerdata->nactivityusesprop > 0 || ownerdata->nactivityusessepa > 0 )
3643  {
3644  SCIP_CALL( detectNlhdlr(scip, expr, expr == consdata->expr ? conss[i] : NULL) );
3645 
3646  assert(ownerdata->nenfos >= 0);
3647  }
3648  else
3649  {
3650  /* remember that we looked at this expression during detectNlhdlrs
3651  * even though we have not actually run detectNlhdlr, because no nlhdlr showed interest in this expr,
3652  * in some situations (forwardPropExpr, to be specific) we will have to distinguish between exprs for which
3653  * we have not initialized enforcement yet (nenfos < 0) and expressions which are just not used in enforcement (nenfos == 0)
3654  */
3655  ownerdata->nenfos = 0;
3656  }
3657  }
3658 
3659  /* include this constraint into the next propagation round because the added nlhdlr may do find tighter bounds now */
3660  if( SCIPconsIsPropagated(conss[i]) )
3661  consdata->ispropagated = FALSE;
3662  }
3663 
3664  if( SCIPgetStage(scip) == SCIP_STAGE_SOLVING && SCIPgetDepth(scip) != 0 )
3665  {
3666  /* ensure that the local bounds are used again when reevaluating the expressions later;
3667  * this is only needed if CONSACTIVE is called in a local node (see begin of this function)
3668  */
3670  conshdlrdata->globalbounds = FALSE;
3671  conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
3672  }
3673  else
3674  {
3675  /* ensure that all activities (except for var-exprs) are reevaluated since better methods may be available now */
3677  }
3678 
3679  SCIPfreeExpriter(&it);
3680 
3681  return SCIP_OKAY;
3682 }
3683 
3684 /** initializes (pre)solving data of constraints
3685  *
3686  * This initializes data in a constraint that is used for separation, propagation, etc, and assumes that expressions will
3687  * not be modified.
3688  * In particular, this function
3689  * - runs the detection method of nlhldrs
3690  * - looks for unlocked linear variables
3691  * - checks curvature (if not in presolve)
3692  * - creates and add row to NLP (if not in presolve)
3693  *
3694  * This function can be called in presolve and solve and can be called several times with different sets of constraints,
3695  * e.g., it should be called in INITSOL and for constraints that are added during solve.
3696  */
3697 static
3699  SCIP* scip, /**< SCIP data structure */
3700  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3701  SCIP_CONS** conss, /**< constraints */
3702  int nconss /**< number of constraints */
3703  )
3704 {
3705  int c;
3706 
3707  for( c = 0; c < nconss; ++c )
3708  {
3709  /* check for a linear variable that can be increase or decreased without harming feasibility */
3710  findUnlockedLinearVar(scip, conss[c]);
3711 
3713  {
3714  SCIP_CONSDATA* consdata;
3715  SCIP_Bool success = FALSE;
3716 
3717  consdata = SCIPconsGetData(conss[c]); /*lint !e613*/
3718  assert(consdata != NULL);
3719  assert(consdata->expr != NULL);
3720 
3721  if( !SCIPconshdlrGetData(conshdlr)->assumeconvex )
3722  {
3723  /* call the curvature detection algorithm of the convex nonlinear handler
3724  * Check only for those curvature that may result in a convex inequality, i.e.,
3725  * whether f(x) is concave when f(x) >= lhs and/or f(x) is convex when f(x) <= rhs.
3726  * Also we can assume that we are nonlinear, so do not check for convex if already concave.
3727  */
3728  if( !SCIPisInfinity(scip, -consdata->lhs) )
3729  {
3730  SCIP_CALL( SCIPhasExprCurvature(scip, consdata->expr, SCIP_EXPRCURV_CONCAVE, &success, NULL) );
3731  if( success )
3732  consdata->curv = SCIP_EXPRCURV_CONCAVE;
3733  }
3734  if( !success && !SCIPisInfinity(scip, consdata->rhs) )
3735  {
3736  SCIP_CALL( SCIPhasExprCurvature(scip, consdata->expr, SCIP_EXPRCURV_CONVEX, &success, NULL) );
3737  if( success )
3738  consdata->curv = SCIP_EXPRCURV_CONVEX;
3739  }
3740  }
3741  else
3742  {
3743  if( !SCIPisInfinity(scip, -consdata->lhs) && !SCIPisInfinity(scip, consdata->rhs) )
3744  {
3745  SCIPwarningMessage(scip, "Nonlinear constraint <%s> has finite left- and right-hand side, but constraints/nonlinear/assumeconvex is enabled.\n", SCIPconsGetName(conss[c]));
3746  consdata->curv = SCIP_EXPRCURV_LINEAR;
3747  }
3748  else
3749  {
3750  consdata->curv = !SCIPisInfinity(scip, consdata->rhs) ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
3751  }
3752  }
3753  SCIPdebugMsg(scip, "root curvature of constraint %s = %d\n", SCIPconsGetName(conss[c]), consdata->curv);
3754 
3755  /* add nlrow representation to NLP, if NLP had been constructed */
3756  if( SCIPisNLPConstructed(scip) && SCIPconsIsActive(conss[c]) )
3757  {
3758  if( consdata->nlrow == NULL )
3759  {
3760  SCIP_CALL( createNlRow(scip, conss[c]) );
3761  assert(consdata->nlrow != NULL);
3762  }
3763  SCIPnlrowSetCurvature(consdata->nlrow, consdata->curv);
3764  SCIP_CALL( SCIPaddNlRow(scip, consdata->nlrow) );
3765  }
3766  }
3767  }
3768 
3769  /* register non linear handlers */
3770  SCIP_CALL( detectNlhdlrs(scip, conshdlr, conss, nconss) );
3771 
3772  return SCIP_OKAY;
3773 }
3774 
3775 /** deinitializes (pre)solving data of constraints
3776  *
3777  * This removes the initialization data created in initSolve().
3778  *
3779  * This function can be called in presolve and solve.
3780  *
3781  * TODO At the moment, it should not be called for a constraint if there are other constraints
3782  * that use the same expressions but still require their nlhdlr.
3783  * We should probably only decrement the auxvar and activity usage for the root expr and then
3784  * proceed as in detectNlhdlrs(), i.e., free enfo data only where none is used.
3785  */
3786 static
3788  SCIP* scip, /**< SCIP data structure */
3789  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3790  SCIP_CONS** conss, /**< constraints */
3791  int nconss /**< number of constraints */
3792  )
3793 {
3794  SCIP_EXPRITER* it;
3795  SCIP_EXPR* expr;
3796  SCIP_CONSDATA* consdata;
3797  SCIP_Bool rootactivityvalid;
3798  int c;
3799 
3800  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
3803 
3804  /* call deinitialization callbacks of expression and nonlinear handlers
3805  * free nonlinear handlers information from expressions
3806  * remove auxiliary variables and nactivityuses counts from expressions
3807  */
3808  for( c = 0; c < nconss; ++c )
3809  {
3810  assert(conss != NULL);
3811  assert(conss[c] != NULL);
3812 
3813  consdata = SCIPconsGetData(conss[c]);
3814  assert(consdata != NULL);
3815  assert(consdata->expr != NULL);
3816 
3817  /* check and remember whether activity in root is valid */
3818  rootactivityvalid = SCIPexprGetActivityTag(consdata->expr) >= SCIPconshdlrGetData(conshdlr)->lastboundrelax;
3819 
3820  for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
3821  {
3822  SCIPdebugMsg(scip, "exitsepa and free nonlinear handler data for expression %p\n", (void*)expr);
3823 
3824  /* remove nonlinear handlers in expression and their data and auxiliary variables; reset activityusage count */
3825  SCIP_CALL( freeEnfoData(scip, expr, TRUE) );
3826 
3827  /* remove quadratic info */
3828  SCIPfreeExprQuadratic(scip, expr);
3829 
3830  if( rootactivityvalid )
3831  {
3832  /* ensure activity is valid if consdata->expr activity is valid
3833  * this is mainly to ensure that we do not leave invalid activities in parts of the expression tree where activity was not used,
3834  * e.g., an expr's activity was kept up to date by a nlhdlr, but without using some childs activity
3835  * so this childs activity would be invalid, which can generate confusion
3836  */
3837  SCIP_CALL( SCIPevalExprActivity(scip, expr) );
3838  }
3839  }
3840 
3841  if( consdata->nlrow != NULL )
3842  {
3843  /* remove row from NLP, if still in solving
3844  * if we are in exitsolve, the whole NLP will be freed anyway
3845  */
3846  if( SCIPgetStage(scip) == SCIP_STAGE_SOLVING )
3847  {
3848  SCIP_CALL( SCIPdelNlRow(scip, consdata->nlrow) );
3849  }
3850 
3851  SCIP_CALL( SCIPreleaseNlRow(scip, &consdata->nlrow) );
3852  }
3853 
3854  /* forget about linear variables that can be increased or decreased without harming feasibility */
3855  consdata->linvardecr = NULL;
3856  consdata->linvarincr = NULL;
3857 
3858  /* forget about curvature */
3859  consdata->curv = SCIP_EXPRCURV_UNKNOWN;
3860  }
3861 
3862  SCIPfreeExpriter(&it);
3863 
3864  return SCIP_OKAY;
3865 }
3866 
3867 /** helper method to decide whether a given expression is product of at least two binary variables */
3868 static
3870  SCIP* scip, /**< SCIP data structure */
3871  SCIP_EXPR* expr /**< expression */
3872  )
3873 {
3874  int i;
3875 
3876  assert(expr != NULL);
3877 
3878  /* check whether the expression is a product */
3879  if( !SCIPisExprProduct(scip, expr) )
3880  return FALSE;
3881 
3882  /* don't consider products with a coefficient != 1 and products with a single child
3883  * simplification will take care of this expression later
3884  */
3885  if( SCIPexprGetNChildren(expr) <= 1 || SCIPgetCoefExprProduct(expr) != 1.0 )
3886  return FALSE;
3887 
3888  for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
3889  {
3890  SCIP_EXPR* child;
3891  SCIP_VAR* var;
3892  SCIP_Real ub;
3893  SCIP_Real lb;
3894 
3895  child = SCIPexprGetChildren(expr)[i];
3896  assert(child != NULL);
3897 
3898  if( !SCIPisExprVar(scip, child) )
3899  return FALSE;
3900 
3901  var = SCIPgetVarExprVar(child);
3902  lb = SCIPvarGetLbLocal(var);
3903  ub = SCIPvarGetUbLocal(var);
3904 
3905  /* check whether variable is integer and has [0,1] as variable bounds */
3906  if( !SCIPvarIsIntegral(var) || !SCIPisEQ(scip, lb, 0.0) || !SCIPisEQ(scip, ub, 1.0) )
3907  return FALSE;
3908  }
3909 
3910  return TRUE;
3911 }
3912 
3913 /** helper method to collect all bilinear binary product terms */
3914 static
3916  SCIP* scip, /**< SCIP data structure */
3917  SCIP_EXPR* sumexpr, /**< sum expression */
3918  SCIP_VAR** xs, /**< array to collect first variable of each bilinear binary product */
3919  SCIP_VAR** ys, /**< array to collect second variable of each bilinear binary product */
3920  int* childidxs, /**< array to store the index of the child of each stored bilinear binary product */
3921  int* nterms /**< pointer to store the total number of bilinear binary terms */
3922  )
3923 {
3924  int i;
3925 
3926  assert(sumexpr != NULL);
3927  assert(SCIPisExprSum(scip, sumexpr));
3928  assert(xs != NULL);
3929  assert(ys != NULL);
3930  assert(childidxs != NULL);
3931  assert(nterms != NULL);
3932 
3933  *nterms = 0;
3934 
3935  for( i = 0; i < SCIPexprGetNChildren(sumexpr); ++i )
3936  {
3937  SCIP_EXPR* child;
3938 
3939  child = SCIPexprGetChildren(sumexpr)[i];
3940  assert(child != NULL);
3941 
3942  if( SCIPexprGetNChildren(child) == 2 && isBinaryProduct(scip, child) )
3943  {
3946 
3947  assert(x != NULL);
3948  assert(y != NULL);
3949 
3950  if( x != y )
3951  {
3952  xs[*nterms] = x;
3953  ys[*nterms] = y;
3954  childidxs[*nterms] = i;
3955  ++(*nterms);
3956  }
3957  }
3958  }
3959 
3960  return SCIP_OKAY;
3961 }
3962 
3963 /** helper method to reformulate \f$x_i \sum_j c_{ij} x_j\f$ */
3964 static
3966  SCIP* scip, /**< SCIP data structure */
3967  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
3968  SCIP_CONS* cons, /**< constraint */
3969  SCIP_VAR* facvar, /**< variable that has been factorized */
3970  SCIP_VAR** vars, /**< variables of sum_j c_ij x_j */
3971  SCIP_Real* coefs, /**< coefficients of sum_j c_ij x_j */
3972  int nvars, /**< total number of variables in sum_j c_ij x_j */
3973  SCIP_EXPR** newexpr, /**< pointer to store the new expression */
3974  int* naddconss /**< pointer to update the total number of added constraints (might be NULL) */
3975  )
3976 {
3977  SCIP_VAR* auxvar;
3978  SCIP_CONS* newcons;
3979  SCIP_Real minact = 0.0;
3980  SCIP_Real maxact = 0.0;
3981  SCIP_Bool integral = TRUE;
3982  char name [SCIP_MAXSTRLEN];
3983  int i;
3984 
3985  assert(facvar != NULL);
3986  assert(vars != NULL);
3987  assert(nvars > 1);
3988  assert(newexpr != NULL);
3989 
3990  /* compute minimum and maximum activity of sum_j c_ij x_j */
3991  /* TODO could compute minact and maxact for facvar=0 and facvar=1 separately, taking implied bounds into account, allowing for possibly tighter big-M's below */
3992  for( i = 0; i < nvars; ++i )
3993  {
3994  minact += MIN(coefs[i], 0.0);
3995  maxact += MAX(coefs[i], 0.0);
3996  integral = integral && SCIPisIntegral(scip, coefs[i]);
3997  }
3998  assert(minact <= maxact);
3999 
4000  /* create and add auxiliary variable */
4001  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s", SCIPconsGetName(cons), SCIPvarGetName(facvar));
4002  SCIP_CALL( SCIPcreateVarBasic(scip, &auxvar, name, minact, maxact, 0.0, integral ? SCIP_VARTYPE_IMPLINT : SCIP_VARTYPE_CONTINUOUS) );
4003  SCIP_CALL( SCIPaddVar(scip, auxvar) );
4004 
4005  /* create and add z - maxact x <= 0 */
4006  if( !SCIPisZero(scip, maxact) )
4007  {
4008  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_1", SCIPconsGetName(cons), SCIPvarGetName(facvar));
4009  SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &newcons, name, auxvar, facvar, -maxact, -SCIPinfinity(scip), 0.0) );
4010  SCIP_CALL( SCIPaddCons(scip, newcons) );
4011  SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
4012  if( naddconss != NULL )
4013  ++(*naddconss);
4014  }
4015 
4016  /* create and add 0 <= z - minact x */
4017  if( !SCIPisZero(scip, minact) )
4018  {
4019  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_2", SCIPconsGetName(cons), SCIPvarGetName(facvar));
4020  SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &newcons, name, auxvar, facvar, -minact, 0.0, SCIPinfinity(scip)) );
4021  SCIP_CALL( SCIPaddCons(scip, newcons) );
4022  SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
4023  if( naddconss != NULL )
4024  ++(*naddconss);
4025  }
4026 
4027  /* create and add minact <= sum_j c_j x_j - z + minact x_i */
4028  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_3", SCIPconsGetName(cons), SCIPvarGetName(facvar));
4029  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &newcons, name, nvars, vars, coefs, minact, SCIPinfinity(scip)) );
4030  SCIP_CALL( SCIPaddCoefLinear(scip, newcons, auxvar, -1.0) );
4031  if( !SCIPisZero(scip, minact) )
4032  {
4033  SCIP_CALL( SCIPaddCoefLinear(scip, newcons, facvar, minact) );
4034  }
4035  SCIP_CALL( SCIPaddCons(scip, newcons) );
4036  SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
4037  if( naddconss != NULL )
4038  ++(*naddconss);
4039 
4040  /* create and add sum_j c_j x_j - z + maxact x_i <= maxact */
4041  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_4", SCIPconsGetName(cons), SCIPvarGetName(facvar));
4042  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &newcons, name, nvars, vars, coefs, -SCIPinfinity(scip), maxact) );
4043  SCIP_CALL( SCIPaddCoefLinear(scip, newcons, auxvar, -1.0) );
4044  if( !SCIPisZero(scip, maxact) )
4045  {
4046  SCIP_CALL( SCIPaddCoefLinear(scip, newcons, facvar, maxact) );
4047  }
4048  SCIP_CALL( SCIPaddCons(scip, newcons) );
4049  SCIP_CALL( SCIPreleaseCons(scip, &newcons) );
4050  if( naddconss != NULL )
4051  ++(*naddconss);
4052 
4053  /* create variable expression */
4054  SCIP_CALL( createExprVar(scip, conshdlr, newexpr, auxvar) );
4055 
4056  /* release auxvar */
4057  SCIP_CALL( SCIPreleaseVar(scip, &auxvar) );
4058 
4059  return SCIP_OKAY;
4060 }
4061 
4062 /** helper method to generate an expression for a sum of products of binary variables; note that the method captures the generated expression */
4063 static
4065  SCIP* scip, /**< SCIP data structure */
4066  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4067  SCIP_CONS* cons, /**< constraint */
4068  SCIP_EXPR* sumexpr, /**< expression */
4069  int minterms, /**< minimum number of terms in a the sum of x_i sum_j c_j x_j */
4070  SCIP_EXPR** newexpr, /**< pointer to store the expression that represents the binary quadratic */
4071  int* naddconss /**< pointer to update the total number of added constraints (might be NULL) */
4072  )
4073 {
4074  SCIP_EXPR** exprs = NULL;
4075  SCIP_VAR** tmpvars = NULL;
4076  SCIP_VAR** vars = NULL;
4077  SCIP_VAR** xs = NULL;
4078  SCIP_VAR** ys = NULL;
4079  SCIP_Real* exprcoefs = NULL;
4080  SCIP_Real* tmpcoefs = NULL;
4081  SCIP_Real* sumcoefs;
4082  SCIP_Bool* isused = NULL;
4083  int* childidxs = NULL;
4084  int* count = NULL;
4085  int nchildren;
4086  int nexprs = 0;
4087  int nterms;
4088  int nvars;
4089  int ntotalvars;
4090  int i;
4091 
4092  assert(sumexpr != NULL);
4093  assert(minterms > 1);
4094  assert(newexpr != NULL);
4095 
4096  *newexpr = NULL;
4097 
4098  /* check whether sumexpr is indeed a sum */
4099  if( !SCIPisExprSum(scip, sumexpr) )
4100  return SCIP_OKAY;
4101 
4102  nchildren = SCIPexprGetNChildren(sumexpr);
4103  sumcoefs = SCIPgetCoefsExprSum(sumexpr);
4104  nvars = SCIPgetNVars(scip);
4105  ntotalvars = SCIPgetNTotalVars(scip);
4106 
4107  /* check whether there are enough terms available */
4108  if( nchildren < minterms )
4109  return SCIP_OKAY;
4110 
4111  /* allocate memory */
4112  SCIP_CALL( SCIPallocBufferArray(scip, &xs, nchildren) );
4113  SCIP_CALL( SCIPallocBufferArray(scip, &ys, nchildren) );
4114  SCIP_CALL( SCIPallocBufferArray(scip, &childidxs, nchildren) );
4115 
4116  /* collect all bilinear binary product terms */
4117  SCIP_CALL( getBilinearBinaryTerms(scip, sumexpr, xs, ys, childidxs, &nterms) );
4118 
4119  /* check whether there are enough terms available */
4120  if( nterms < minterms )
4121  goto TERMINATE;
4122 
4123  /* store how often each variable appears in a bilinear binary product */
4124  SCIP_CALL( SCIPduplicateBufferArray(scip, &vars, SCIPgetVars(scip), nvars) );
4125  SCIP_CALL( SCIPallocClearBufferArray(scip, &count, ntotalvars) );
4126  SCIP_CALL( SCIPallocClearBufferArray(scip, &isused, nchildren) );
4127 
4128  SCIP_CALL( SCIPallocBufferArray(scip, &exprs, nchildren) );
4129  SCIP_CALL( SCIPallocBufferArray(scip, &exprcoefs, nchildren) );
4130  SCIP_CALL( SCIPallocBufferArray(scip, &tmpvars, MIN(nterms, nvars)) );
4131  SCIP_CALL( SCIPallocBufferArray(scip, &tmpcoefs, MIN(nterms, nvars)) );
4132 
4133  for( i = 0; i < nterms; ++i )
4134  {
4135  int xidx;
4136  int yidx;
4137 
4138  assert(xs[i] != NULL);
4139  assert(ys[i] != NULL);
4140 
4141  xidx = SCIPvarGetIndex(xs[i]);
4142  assert(xidx < ntotalvars);
4143  yidx = SCIPvarGetIndex(ys[i]);
4144  assert(yidx < ntotalvars);
4145 
4146  ++count[xidx];
4147  ++count[yidx];
4148 
4149  SCIPdebugMsg(scip, "increase counter for %s to %d\n", SCIPvarGetName(xs[i]), count[xidx]);
4150  SCIPdebugMsg(scip, "increase counter for %s to %d\n", SCIPvarGetName(ys[i]), count[yidx]);
4151  }
4152 
4153  /* sort variables; don't change order of count array because it depends on problem indices */
4154  {
4155  int* tmpcount;
4156 
4157  SCIP_CALL( SCIPduplicateBufferArray(scip, &tmpcount, count, nvars) );
4158  SCIPsortDownIntPtr(tmpcount, (void**)vars, nvars);
4159  SCIPfreeBufferArray(scip, &tmpcount);
4160  }
4161 
4162  for( i = 0; i < nvars; ++i )
4163  {
4164  SCIP_VAR* facvar = vars[i];
4165  int ntmpvars = 0;
4166  int j;
4167 
4168  /* skip candidate if there are not enough terms left */
4169  if( count[SCIPvarGetIndex(vars[i])] < minterms )
4170  continue;
4171 
4172  SCIPdebugMsg(scip, "consider facvar = %s with count = %d\n", SCIPvarGetName(facvar), count[SCIPvarGetIndex(vars[i])]);
4173 
4174  /* collect variables for x_i * sum_j c_ij x_j */
4175  for( j = 0; j < nterms; ++j )
4176  {
4177  int childidx = childidxs[j];
4178  assert(childidx >= 0 && childidx < nchildren);
4179 
4180  if( !isused[childidx] && (xs[j] == facvar || ys[j] == facvar) )
4181  {
4182  SCIP_Real coef;
4183  int xidx;
4184  int yidx;
4185 
4186  coef = sumcoefs[childidx];
4187  assert(coef != 0.0);
4188 
4189  /* collect corresponding variable */
4190  tmpvars[ntmpvars] = (xs[j] == facvar) ? ys[j] : xs[j];
4191  tmpcoefs[ntmpvars] = coef;
4192  ++ntmpvars;
4193 
4194  /* update counters */
4195  xidx = SCIPvarGetIndex(xs[j]);
4196  assert(xidx < ntotalvars);
4197  yidx = SCIPvarGetIndex(ys[j]);
4198  assert(yidx < ntotalvars);
4199  --count[xidx];
4200  --count[yidx];
4201  assert(count[xidx] >= 0);
4202  assert(count[yidx] >= 0);
4203 
4204  /* mark term to be used */
4205  isused[childidx] = TRUE;
4206  }
4207  }
4208  assert(ntmpvars >= minterms);
4209  assert(SCIPvarGetIndex(facvar) < ntotalvars);
4210  assert(count[SCIPvarGetIndex(facvar)] == 0); /* facvar should not appear in any other bilinear term */
4211 
4212  /* create required constraints and store the generated expression */
4213  SCIP_CALL( reformulateFactorizedBinaryQuadratic(scip, conshdlr, cons, facvar, tmpvars, tmpcoefs, ntmpvars, &exprs[nexprs], naddconss) );
4214  exprcoefs[nexprs] = 1.0;
4215  ++nexprs;
4216  }
4217 
4218  /* factorization was only successful if at least one expression has been generated */
4219  if( nexprs > 0 )
4220  {
4221  int nexprsold = nexprs;
4222 
4223  /* add all children of the sum that have not been used */
4224  for( i = 0; i < nchildren; ++i )
4225  {
4226  if( !isused[i] )
4227  {
4228  exprs[nexprs] = SCIPexprGetChildren(sumexpr)[i];
4229  exprcoefs[nexprs] = sumcoefs[i];
4230  ++nexprs;
4231  }
4232  }
4233 
4234  /* create a new sum expression */
4235  SCIP_CALL( SCIPcreateExprSum(scip, newexpr, nexprs, exprs, exprcoefs, SCIPgetConstantExprSum(sumexpr), exprownerCreate, (void*)conshdlr) );
4236 
4237  /* release all expressions that have been generated by reformulateFactorizedBinaryQuadratic() */
4238  for( i = 0; i < nexprsold; ++i )
4239  {
4240  SCIP_CALL( SCIPreleaseExpr(scip, &exprs[i]) );
4241  }
4242  }
4243 
4244 TERMINATE:
4245  /* free memory */
4246  SCIPfreeBufferArrayNull(scip, &tmpcoefs);
4247  SCIPfreeBufferArrayNull(scip, &tmpvars);
4248  SCIPfreeBufferArrayNull(scip, &exprcoefs);
4249  SCIPfreeBufferArrayNull(scip, &exprs);
4250  SCIPfreeBufferArrayNull(scip, &vars);
4251  SCIPfreeBufferArrayNull(scip, &isused);
4252  SCIPfreeBufferArrayNull(scip, &count);
4253  SCIPfreeBufferArray(scip, &childidxs);
4254  SCIPfreeBufferArray(scip, &ys);
4255  SCIPfreeBufferArray(scip, &xs);
4256 
4257  return SCIP_OKAY;
4258 }
4259 
4260 /** helper method to create an AND constraint or varbound constraints for a given binary product expression */
4261 static
4263  SCIP* scip, /**< SCIP data structure */
4264  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4265  SCIP_EXPR* prodexpr, /**< product expression */
4266  SCIP_EXPR** newexpr, /**< pointer to store the expression that represents the product */
4267  int* naddconss, /**< pointer to update the total number of added constraints (might be NULL) */
4268  SCIP_Bool empathy4and /**< whether to use an AND constraint, if possible */
4269  )
4270 {
4271  SCIP_VAR** vars;
4272  SCIP_CONS* cons;
4273  SCIP_Real* coefs;
4274  SCIP_VAR* w;
4275  char name[SCIP_MAXSTRLEN];
4276  int nchildren;
4277  int i;
4278 
4279  assert(conshdlr != NULL);
4280  assert(prodexpr != NULL);
4281  assert(SCIPisExprProduct(scip, prodexpr));
4282  assert(newexpr != NULL);
4283 
4284  nchildren = SCIPexprGetNChildren(prodexpr);
4285  assert(nchildren >= 2);
4286 
4287  /* memory to store the variables of the variable expressions (+1 for w) */
4288  SCIP_CALL( SCIPallocBufferArray(scip, &vars, nchildren + 1) );
4289  SCIP_CALL( SCIPallocBufferArray(scip, &coefs, nchildren + 1) );
4290 
4291  /* prepare the names of the variable and the constraints */
4292  (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform");
4293  for( i = 0; i < nchildren; ++i )
4294  {
4295  vars[i] = SCIPgetVarExprVar(SCIPexprGetChildren(prodexpr)[i]);
4296  coefs[i] = 1.0;
4297  assert(vars[i] != NULL);
4298  (void) strcat(name, "_");
4299  (void) strcat(name, SCIPvarGetName(vars[i]));
4300  }
4301 
4302  /* create and add variable */
4303  SCIP_CALL( SCIPcreateVarBasic(scip, &w, name, 0.0, 1.0, 0.0, SCIP_VARTYPE_IMPLINT) );
4304  SCIP_CALL( SCIPaddVar(scip, w) );
4305  SCIPdebugMsg(scip, " created auxiliary variable %s\n", name);
4306 
4307  /* use variable bound constraints if it is a bilinear product and there is no empathy for an AND constraint */
4308  if( nchildren == 2 && !empathy4and )
4309  {
4310  SCIP_VAR* x = vars[0];
4311  SCIP_VAR* y = vars[1];
4312 
4313  assert(x != NULL);
4314  assert(y != NULL);
4315  assert(x != y);
4316 
4317  /* create and add x - w >= 0 */
4318  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_1", SCIPvarGetName(x), SCIPvarGetName(y));
4319  SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &cons, name, x, w, -1.0, 0.0, SCIPinfinity(scip)) );
4320  SCIP_CALL( SCIPaddCons(scip, cons) );
4321  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
4322 
4323  /* create and add y - w >= 0 */
4324  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_2", SCIPvarGetName(x), SCIPvarGetName(y));
4325  SCIP_CALL( SCIPcreateConsBasicVarbound(scip, &cons, name, y, w, -1.0, 0.0, SCIPinfinity(scip)) );
4326  SCIP_CALL( SCIPaddCons(scip, cons) );
4327  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
4328 
4329  /* create and add x + y - w <= 1 */
4330  vars[2] = w;
4331  coefs[2] = -1.0;
4332  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "binreform_%s_%s_3", SCIPvarGetName(x), SCIPvarGetName(y));
4333  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 3, vars, coefs, -SCIPinfinity(scip), 1.0) );
4334  SCIP_CALL( SCIPaddCons(scip, cons) );
4335  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
4336 
4337  /* update number of added constraints */
4338  if( naddconss != NULL )
4339  *naddconss += 3;
4340  }
4341  else
4342  {
4343  /* create, add, and release AND constraint */
4344  SCIP_CALL( SCIPcreateConsBasicAnd(scip, &cons, name, w, nchildren, vars) );
4345  SCIP_CALL( SCIPaddCons(scip, cons) );
4346  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
4347  SCIPdebugMsg(scip, " create AND constraint\n");
4348 
4349  /* update number of added constraints */
4350  if( naddconss != NULL )
4351  *naddconss += 1;
4352  }
4353 
4354  /* create variable expression */
4355  SCIP_CALL( createExprVar(scip, conshdlr, newexpr, w) );
4356 
4357  /* release created variable */
4358  SCIP_CALL( SCIPreleaseVar(scip, &w) );
4359 
4360  /* free memory */
4361  SCIPfreeBufferArray(scip, &coefs);
4362  SCIPfreeBufferArray(scip, &vars);
4363 
4364  return SCIP_OKAY;
4365 }
4366 
4367 /** helper method to generate an expression for the product of binary variables; note that the method captures the generated expression */
4368 static
4370  SCIP* scip, /**< SCIP data structure */
4371  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4372  SCIP_HASHMAP* exprmap, /**< map to remember generated variables for visited product expressions */
4373  SCIP_EXPR* prodexpr, /**< product expression */
4374  SCIP_EXPR** newexpr, /**< pointer to store the expression that represents the product */
4375  int* naddconss, /**< pointer to update the total number of added constraints (might be NULL) */
4376  int* nchgcoefs /**< pointer to update the total number of changed coefficients (might be NULL) */
4377  )
4378 {
4379  SCIP_CONSHDLRDATA* conshdlrdata;
4380  int nchildren;
4381 
4382  assert(prodexpr != NULL);
4383  assert(newexpr != NULL);
4384 
4385  *newexpr = NULL;
4386 
4387  /* only consider products of binary variables */
4388  if( !isBinaryProduct(scip, prodexpr) )
4389  return SCIP_OKAY;
4390 
4391  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4392  assert(conshdlrdata != NULL);
4393  nchildren = SCIPexprGetNChildren(prodexpr);
4394  assert(nchildren >= 2);
4395 
4396  /* check whether there is already an expression that represents the product */
4397  if( SCIPhashmapExists(exprmap, (void*)prodexpr) )
4398  {
4399  *newexpr = (SCIP_EXPR*) SCIPhashmapGetImage(exprmap, (void*)prodexpr);
4400  assert(*newexpr != NULL);
4401 
4402  /* capture expression */
4403  SCIPcaptureExpr(*newexpr);
4404  }
4405  else
4406  {
4407  SCIPdebugMsg(scip, " product expression %p has been considered for the first time\n", (void*)prodexpr);
4408 
4409  if( nchildren == 2 )
4410  {
4411  SCIP_CLIQUE** xcliques;
4412  SCIP_VAR* x;
4413  SCIP_VAR* y;
4414  SCIP_Bool found_clique = FALSE;
4415  int c;
4416 
4417  /* get variables from the product expression */
4418  x = SCIPgetVarExprVar(SCIPexprGetChildren(prodexpr)[0]);
4419  assert(x != NULL);
4420  y = SCIPgetVarExprVar(SCIPexprGetChildren(prodexpr)[1]);
4421  assert(y != NULL);
4422  assert(x != y);
4423 
4424  /* first try to find a clique containing both variables */
4425  xcliques = SCIPvarGetCliques(x, TRUE);
4426 
4427  /* look in cliques containing x */
4428  for( c = 0; c < SCIPvarGetNCliques(x, TRUE); ++c )
4429  {
4430  if( SCIPcliqueHasVar(xcliques[c], y, TRUE) ) /* x + y <= 1 => x*y = 0 */
4431  {
4432  /* create zero value expression */
4433  SCIP_CALL( SCIPcreateExprValue(scip, newexpr, 0.0, exprownerCreate, (void*)conshdlr) );
4434 
4435  if( nchgcoefs != NULL )
4436  *nchgcoefs += 1;
4437 
4438  found_clique = TRUE;
4439  break;
4440  }
4441 
4442  if( SCIPcliqueHasVar(xcliques[c], y, FALSE) ) /* x + (1-y) <= 1 => x*y = x */
4443  {
4444  /* create variable expression for x */
4445  SCIP_CALL( createExprVar(scip, conshdlr, newexpr, x) );
4446 
4447  if( nchgcoefs != NULL )
4448  *nchgcoefs += 2;
4449 
4450  found_clique = TRUE;
4451  break;
4452  }
4453  }
4454 
4455  if( !found_clique )
4456  {
4457  xcliques = SCIPvarGetCliques(x, FALSE);
4458 
4459  /* look in cliques containing complement of x */
4460  for( c = 0; c < SCIPvarGetNCliques(x, FALSE); ++c )
4461  {
4462  if( SCIPcliqueHasVar(xcliques[c], y, TRUE) ) /* (1-x) + y <= 1 => x*y = y */
4463  {
4464  /* create variable expression for y */
4465  SCIP_CALL( createExprVar(scip, conshdlr, newexpr, y) );
4466 
4467  if( nchgcoefs != NULL )
4468  *nchgcoefs += 1;
4469 
4470  found_clique = TRUE;
4471  break;
4472  }
4473 
4474  if( SCIPcliqueHasVar(xcliques[c], y, FALSE) ) /* (1-x) + (1-y) <= 1 => x*y = x + y - 1 */
4475  {
4476  /* create sum expression */
4477  SCIP_EXPR* sum_children[2];
4478  SCIP_Real sum_coefs[2];
4479  SCIP_CALL( createExprVar(scip, conshdlr, &sum_children[0], x) );
4480  SCIP_CALL( createExprVar(scip, conshdlr, &sum_children[1], y) );
4481  sum_coefs[0] = 1.0;
4482  sum_coefs[1] = 1.0;
4483  SCIP_CALL( SCIPcreateExprSum(scip, newexpr, 2, sum_children, sum_coefs, -1.0, exprownerCreate, (void*)conshdlr) );
4484 
4485  SCIP_CALL( SCIPreleaseExpr(scip, &sum_children[0]) );
4486  SCIP_CALL( SCIPreleaseExpr(scip, &sum_children[1]) );
4487 
4488  if( nchgcoefs != NULL )
4489  *nchgcoefs += 3;
4490 
4491  found_clique = TRUE;
4492  break;
4493  }
4494  }
4495  }
4496 
4497  /* if the variables are not in a clique, do standard linearization */
4498  if( !found_clique )
4499  {
4500  SCIP_CALL( getBinaryProductExprDo(scip, conshdlr, prodexpr, newexpr, naddconss, conshdlrdata->reformbinprodsand) );
4501  }
4502  }
4503  else
4504  {
4505  /* linearize binary product using an AND constraint because nchildren > 2 */
4506  SCIP_CALL( getBinaryProductExprDo(scip, conshdlr, prodexpr, newexpr, naddconss, conshdlrdata->reformbinprodsand) );
4507  }
4508 
4509  /* hash variable expression */
4510  SCIP_CALL( SCIPhashmapInsert(exprmap, (void*)prodexpr, *newexpr) );
4511  }
4512 
4513  return SCIP_OKAY;
4514 }
4515 
4516 /** helper function to replace binary products in a given constraint */
4517 static
4519  SCIP* scip, /**< SCIP data structure */
4520  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4521  SCIP_CONS* cons, /**< constraint */
4522  SCIP_HASHMAP* exprmap, /**< map to remember generated variables for visited product expressions */
4523  SCIP_EXPRITER* it, /**< expression iterator */
4524  int* naddconss, /**< pointer to update the total number of added constraints (might be NULL) */
4525  int* nchgcoefs /**< pointer to update the total number of changed coefficients (might be NULL) */
4526  )
4527 {
4528  SCIP_CONSHDLRDATA* conshdlrdata;
4529  SCIP_CONSDATA* consdata;
4530  SCIP_EXPR* expr;
4531 
4532  assert(conshdlr != NULL);
4533  assert(cons != NULL);
4534  assert(exprmap != NULL);
4535  assert(it != NULL);
4536 
4537  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4538  assert(conshdlrdata != NULL);
4539 
4540  consdata = SCIPconsGetData(cons);
4541  assert(consdata != NULL);
4542  assert(consdata->expr != NULL);
4543 
4544  SCIPdebugMsg(scip, " check constraint %s\n", SCIPconsGetName(cons));
4545 
4546  for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
4547  {
4548  SCIP_EXPR* newexpr = NULL;
4549  SCIP_EXPR* childexpr;
4550  int childexpridx;
4551 
4552  childexpridx = SCIPexpriterGetChildIdxDFS(it);
4553  assert(childexpridx >= 0 && childexpridx < SCIPexprGetNChildren(expr));
4554  childexpr = SCIPexpriterGetChildExprDFS(it);
4555  assert(childexpr != NULL);
4556 
4557  /* try to factorize variables in a sum expression that contains several products of binary variables */
4558  if( conshdlrdata->reformbinprodsfac > 1 )
4559  {
4560  SCIP_CALL( getFactorizedBinaryQuadraticExpr(scip, conshdlr, cons, childexpr, conshdlrdata->reformbinprodsfac, &newexpr, naddconss) );
4561  }
4562 
4563  /* try to create an expression that represents a product of binary variables */
4564  if( newexpr == NULL )
4565  {
4566  SCIP_CALL( getBinaryProductExpr(scip, conshdlr, exprmap, childexpr, &newexpr, naddconss, nchgcoefs) );
4567  }
4568 
4569  if( newexpr != NULL )
4570  {
4571  assert(naddconss == NULL || *naddconss > 0 || nchgcoefs == NULL || *nchgcoefs > 0);
4572 
4573  /* replace product expression */
4574  SCIP_CALL( SCIPreplaceExprChild(scip, expr, childexpridx, newexpr) );
4575 
4576  /* note that the expression has been captured by getBinaryProductExpr and SCIPreplaceExprChild */
4577  SCIP_CALL( SCIPreleaseExpr(scip, &newexpr) );
4578 
4579  /* mark the constraint to not be simplified anymore */
4580  consdata->issimplified = FALSE;
4581  }
4582  }
4583 
4584  return SCIP_OKAY;
4585 }
4586 
4587 /** reformulates products of binary variables during presolving in the following way:
4588  *
4589  * Let \f$\sum_{i,j} Q_{ij} x_i x_j\f$ be a subexpression that only contains binary variables.
4590  * Each term \f$x_i x_j\f$ is reformulated with the help of an extra (implicit integer) variable \f$z_{ij}\f$ in {0,1}:
4591  * \f[
4592  * z_{ij} \leq x_i, \qquad z_{ij} \leq x_j, \qquad x_i + x_j - z_{ij} \leq 1.
4593  * \f]
4594  *
4595  * Before reformulating \f$x_i x_j\f$ in this way, it is checked whether there is a clique that contains \f$x_i\f$ and \f$x_j\f$.
4596  * These cliques allow for a better reformulation. There are four cases:
4597  *
4598  * 1. \f$x_i + x_j \leq 1\f$ implies that \f$x_i x_j = 0\f$
4599  * 2. \f$x_i + (1 - x_j) \leq 1\f$ implies \f$x_i x_j = x_i\f$
4600  * 3. \f$(1 - x_i) + x_j \leq 1\f$ implies \f$x_i x_j = x_j\f$
4601  * 4. \f$(1 - x_i) + (1 - x_j) \leq 1\f$ implies \f$x_i x_j = x_i + x_j - 1\f$
4602  *
4603  * The reformulation using \f$z_{ij}\f$ or the cliques is implemented in getBinaryProductExpr().
4604  *
4605  * Introducing too many extra variables and constraints can have a negative impact on the performance (e.g., due to
4606  * slow probing). For this reason, it is checked in getFactorizedBinaryQuadraticExpr() whether \f$\sum_{i,j} Q_{ij} x_i x_j\f$
4607  * contains large (&ge; `reformbinprodsfac` parameter) lower sums of the form \f$x_i \sum_j Q_{ij} x_j\f$.
4608  * Such a lower sum is reformulated with only one extra variable w_i:
4609  * \f{align}{
4610  * \text{maxact} & := \sum_j \max(0, Q_{ij}), \\
4611  * \text{minact} & := \sum_j \min(0, Q_{ij}), \\
4612  * \text{minact}\, x_i & \leq w_i, \\
4613  * w_i &\leq \text{maxact}\, x_i, \\
4614  * \text{minact} &\leq \sum_j Q_{ij} x_j - w_i + \text{minact}\, x_i \\
4615  * \text{maxact} &\geq \sum_j Q_{ij} x_j - w_i + \text{maxact}\, x_i
4616  * \f}
4617  * We mark \f$w_i\f$ to be implicit integer if all \f$Q_{ij}\f$ are integer. After each replacement of a lower sum, it
4618  * is checked whether there are enough terms left to factorize other binary variables. Lower sums with a larger number
4619  * of terms are prioritized.
4620  */
4621 static
4623  SCIP* scip, /**< SCIP data structure */
4624  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4625  SCIP_CONS** conss, /**< constraints */
4626  int nconss, /**< total number of constraints */
4627  int* naddconss, /**< pointer to store the total number of added constraints (might be NULL) */
4628  int* nchgcoefs /**< pointer to store the total number of changed coefficients (might be NULL) */
4629  )
4630 {
4631  SCIP_CONSHDLRDATA* conshdlrdata;
4632  SCIP_HASHMAP* exprmap;
4633  SCIP_EXPRITER* it;
4634  int c;
4635 
4636  assert(conshdlr != NULL);
4637 
4638  /* no nonlinear constraints or binary variables -> skip */
4639  if( nconss == 0 || SCIPgetNBinVars(scip) == 0 )
4640  return SCIP_OKAY;
4641  assert(conss != NULL);
4642 
4643  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4644  assert(conshdlrdata != NULL);
4645 
4646  /* create expression hash map */
4647  SCIP_CALL( SCIPhashmapCreate(&exprmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );
4648 
4649  /* create expression iterator */
4650  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
4653 
4654  SCIPdebugMsg(scip, "call presolveBinaryProducts()\n");
4655 
4656  for( c = 0; c < nconss; ++c )
4657  {
4658  SCIP_CONSDATA* consdata;
4659  SCIP_EXPR* newexpr = NULL;
4660 
4661  assert(conss[c] != NULL);
4662 
4663  consdata = SCIPconsGetData(conss[c]);
4664  assert(consdata != NULL);
4665 
4666  /* try to reformulate the root expression */
4667  if( conshdlrdata->reformbinprodsfac > 1 )
4668  {
4669  SCIP_CALL( getFactorizedBinaryQuadraticExpr(scip, conshdlr, conss[c], consdata->expr, conshdlrdata->reformbinprodsfac, &newexpr, naddconss) );
4670  }
4671 
4672  /* release the root node if another expression has been found */
4673  if( newexpr != NULL )
4674  {
4675  SCIP_CALL( SCIPreleaseExpr(scip, &consdata->expr) );
4676  consdata->expr = newexpr;
4677 
4678  /* mark constraint to be not simplified anymore */
4679  consdata->issimplified = FALSE;
4680  }
4681 
4682  /* replace each product of binary variables separately */
4683  SCIP_CALL( replaceBinaryProducts(scip, conshdlr, conss[c], exprmap, it, naddconss, nchgcoefs) );
4684  }
4685 
4686  /* free memory */
4687  SCIPhashmapFree(&exprmap);
4688  SCIPfreeExpriter(&it);
4689 
4690  return SCIP_OKAY;
4691 }
4692 
4693 /** scales the sides of the constraint \f$\ell \leq \sum_i c_i f_i(x) \leq r\f$.
4694  *
4695  * Let \f$n_+\f$ the number of positive coefficients \f$c_i\f$ and \f$n_-\f$ be the number of negative coefficients.
4696  * Then scale by -1 if
4697  * - \f$n_+ < n_-\f$, or
4698  * - \f$n_+ = n_-\f$ and \f$r = \infty\f$.
4699  */
4700 static
4702  SCIP* scip, /**< SCIP data structure */
4703  SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */
4704  SCIP_CONS* cons, /**< nonlinear constraint */
4705  SCIP_Bool* changed /**< buffer to store if the expression of cons changed */
4706  )
4707 {
4708  SCIP_CONSDATA* consdata;
4709  int i;
4710 
4711  assert(cons != NULL);
4712 
4713  consdata = SCIPconsGetData(cons);
4714  assert(consdata != NULL);
4715 
4716  if( SCIPisExprSum(scip, consdata->expr) )
4717  {
4718  SCIP_Real* coefs;
4719  SCIP_Real constant;
4720  int nchildren;
4721  int counter = 0;
4722 
4723  coefs = SCIPgetCoefsExprSum(consdata->expr);
4724  constant = SCIPgetConstantExprSum(consdata->expr);
4725  nchildren = SCIPexprGetNChildren(consdata->expr);
4726 
4727  /* handle special case when constraint is l <= -f(x) <= r and f(x) not a sum: simplfy ensures f is not a sum */
4728  if( nchildren == 1 && constant == 0.0 && coefs[0] == -1.0 )
4729  {
4730  SCIP_EXPR* expr;
4731  expr = consdata->expr;
4732 
4733  consdata->expr = SCIPexprGetChildren(expr)[0];
4734  assert(!SCIPisExprSum(scip, consdata->expr));
4735 
4736  SCIPcaptureExpr(consdata->expr);
4737 
4738  SCIPswapReals(&consdata->lhs, &consdata->rhs);
4739  consdata->lhs = -consdata->lhs;
4740  consdata->rhs = -consdata->rhs;
4741 
4742  SCIP_CALL( SCIPreleaseExpr(scip, &expr) );
4743  *changed = TRUE;
4744  return SCIP_OKAY;
4745  }
4746 
4747  /* compute n_+ - n_i */
4748  for( i = 0; i < nchildren; ++i )
4749  counter += coefs[i] > 0 ? 1 : -1;
4750 
4751  if( counter < 0 || (counter == 0 && SCIPisInfinity(scip, consdata->rhs)) )
4752  {
4753  SCIP_EXPR* expr;
4754  SCIP_Real* newcoefs;
4755 
4756  /* allocate memory */
4757  SCIP_CALL( SCIPallocBufferArray(scip, &newcoefs, nchildren) );
4758 
4759  for( i = 0; i < nchildren; ++i )
4760  newcoefs[i] = -coefs[i];
4761 
4762  /* create a new sum expression */
4763  SCIP_CALL( SCIPcreateExprSum(scip, &expr, nchildren, SCIPexprGetChildren(consdata->expr), newcoefs, -constant, exprownerCreate, (void*)conshdlr) );
4764 
4765  /* replace expression in constraint data and scale sides */
4766  SCIP_CALL( SCIPreleaseExpr(scip, &consdata->expr) );
4767  consdata->expr = expr;
4768  SCIPswapReals(&consdata->lhs, &consdata->rhs);
4769  consdata->lhs = -consdata->lhs;
4770  consdata->rhs = -consdata->rhs;
4771 
4772  /* free memory */
4773  SCIPfreeBufferArray(scip, &newcoefs);
4774 
4775  *changed = TRUE;
4776  }
4777  }
4778 
4779  return SCIP_OKAY;
4780 }
4781 
4782 /** forbid multiaggrations of variables that appear nonlinear in constraints */
4783 static
4785  SCIP* scip, /**< SCIP data structure */
4786  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4787  SCIP_CONS** conss, /**< constraints */
4788  int nconss /**< number of constraints */
4789  )
4790 {
4791  SCIP_EXPRITER* it;
4792  SCIP_CONSDATA* consdata;
4793  SCIP_EXPR* expr;
4794  int c;
4795 
4796  assert(scip != NULL);
4797  assert(conshdlr != NULL);
4798 
4799  if( !SCIPconshdlrGetData(conshdlr)->forbidmultaggrnlvar )
4800  return SCIP_OKAY;
4801 
4802  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
4804 
4805  for( c = 0; c < nconss; ++c )
4806  {
4807  consdata = SCIPconsGetData(conss[c]);
4808  assert(consdata != NULL);
4809 
4810  /* if root expression is sum, then forbid multiaggregation only for variables that are not in linear terms of sum,
4811  * i.e., skip children of sum that are variables
4812  */
4813  if( SCIPisExprSum(scip, consdata->expr) )
4814  {
4815  int i;
4816  SCIP_EXPR* child;
4817  for( i = 0; i < SCIPexprGetNChildren(consdata->expr); ++i )
4818  {
4819  child = SCIPexprGetChildren(consdata->expr)[i];
4820 
4821  /* skip variable expression, as they correspond to a linear term */
4822  if( SCIPisExprVar(scip, child) )
4823  continue;
4824 
4825  for( expr = SCIPexpriterRestartDFS(it, child); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
4826  if( SCIPisExprVar(scip, expr) )
4827  {
4829  }
4830  }
4831  }
4832  else
4833  {
4834  for( expr = SCIPexpriterRestartDFS(it, consdata->expr); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
4835  if( SCIPisExprVar(scip, expr) )
4836  {
4838  }
4839  }
4840  }
4841 
4842  SCIPfreeExpriter(&it);
4843 
4844  return SCIP_OKAY;
4845 }
4846 
4847 /** simplifies expressions and replaces common subexpressions for a set of constraints
4848  * @todo put the constant to the constraint sides
4849  */
4850 static
4852  SCIP* scip, /**< SCIP data structure */
4853  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
4854  SCIP_CONS** conss, /**< constraints */
4855  int nconss, /**< total number of constraints */
4856  SCIP_PRESOLTIMING presoltiming, /**< presolve timing (SCIP_PRESOLTIMING_ALWAYS if not in presolving) */
4857  SCIP_Bool* infeasible, /**< buffer to store whether infeasibility has been detected */
4858  int* ndelconss, /**< counter to add number of deleted constraints, or NULL */
4859  int* naddconss, /**< counter to add number of added constraints, or NULL */
4860  int* nchgcoefs /**< counter to add number of changed coefficients, or NULL */
4861  )
4862 {
4863  SCIP_CONSHDLRDATA* conshdlrdata;
4864  SCIP_CONSDATA* consdata;
4865  int* nlockspos;
4866  int* nlocksneg;
4867  SCIP_Bool havechange;
4868  int i;
4869 
4870  assert(scip != NULL);
4871  assert(conshdlr != NULL);
4872  assert(conss != NULL);
4873  assert(nconss > 0);
4874  assert(infeasible != NULL);
4875 
4876  conshdlrdata = SCIPconshdlrGetData(conshdlr);
4877  assert(conshdlrdata != NULL);
4878 
4879  /* update number of canonicalize calls */
4880  ++(conshdlrdata->ncanonicalizecalls);
4881 
4882  SCIP_CALL( SCIPstartClock(scip, conshdlrdata->canonicalizetime) );
4883 
4884  *infeasible = FALSE;
4885 
4886  /* set havechange to TRUE in the first call of canonicalize; otherwise we might not replace common subexpressions */
4887  havechange = conshdlrdata->ncanonicalizecalls == 1;
4888 
4889  /* free nonlinear handlers information from expressions */ /* TODO can skip this in first presolve round */
4890  SCIP_CALL( deinitSolve(scip, conshdlr, conss, nconss) );
4891 
4892  /* allocate memory for storing locks of each constraint */
4893  SCIP_CALL( SCIPallocBufferArray(scip, &nlockspos, nconss) );
4894  SCIP_CALL( SCIPallocBufferArray(scip, &nlocksneg, nconss) );
4895 
4896  /* unlock all constraints */
4897  for( i = 0; i < nconss; ++i )
4898  {
4899  assert(conss[i] != NULL);
4900 
4901  consdata = SCIPconsGetData(conss[i]);
4902  assert(consdata != NULL);
4903 
4904  /* remember locks */
4905  nlockspos[i] = consdata->nlockspos;
4906  nlocksneg[i] = consdata->nlocksneg;
4907 
4908  /* remove locks */
4909  SCIP_CALL( addLocks(scip, conss[i], -consdata->nlockspos, -consdata->nlocksneg) );
4910  assert(consdata->nlockspos == 0);
4911  assert(consdata->nlocksneg == 0);
4912  }
4913 
4914 #ifndef NDEBUG
4915  /* check whether all locks of each expression have been removed */
4916  for( i = 0; i < nconss; ++i )
4917  {
4918  SCIP_EXPR* expr;
4919  SCIP_EXPRITER* it;
4920 
4921  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
4922 
4923  consdata = SCIPconsGetData(conss[i]);
4924  assert(consdata != NULL);
4925 
4926  SCIP_CALL( SCIPexpriterInit(it, consdata->expr, SCIP_EXPRITER_RTOPOLOGIC, TRUE) );
4927  for( expr = consdata->expr; !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
4928  {
4929  assert(expr != NULL);
4930  assert(SCIPexprGetOwnerData(expr)->nlocksneg == 0);
4931  assert(SCIPexprGetOwnerData(expr)->nlockspos == 0);
4932  }
4933  SCIPfreeExpriter(&it);
4934  }
4935 #endif
4936 
4937  /* reformulate products of binary variables */
4938  if( conshdlrdata->reformbinprods && SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING
4939  && (presoltiming & SCIP_PRESOLTIMING_EXHAUSTIVE) )
4940  {
4941  int tmpnaddconss = 0;
4942  int tmpnchgcoefs = 0;
4943 
4944  /* call this function before simplification because expressions might not be simplified after reformulating
4945  * binary products; the detection of some nonlinear handlers might assume that expressions are simplified
4946  */
4947  SCIP_CALL( presolveBinaryProducts(scip, conshdlr, conss, nconss, &tmpnaddconss, &tmpnchgcoefs) );
4948 
4949  /* update counters */
4950  if( naddconss != NULL )
4951  *naddconss = tmpnaddconss;
4952  if( nchgcoefs != NULL )
4953  *nchgcoefs = tmpnchgcoefs;
4954 
4955  /* check whether at least one expression has changed */
4956  if( tmpnaddconss + tmpnchgcoefs > 0 )
4957  havechange = TRUE;
4958  }
4959 
4960  for( i = 0; i < nconss; ++i )
4961  {
4962  consdata = SCIPconsGetData(conss[i]);
4963  assert(consdata != NULL);
4964 
4965  /* call simplify for each expression */
4966  if( !consdata->issimplified && consdata->expr != NULL )
4967  {
4968  SCIP_EXPR* simplified;
4969  SCIP_Bool changed;
4970 
4971  changed = FALSE;
4972  SCIP_CALL( SCIPsimplifyExpr(scip, consdata->expr, &simplified, &changed, infeasible, exprownerCreate, (void*)conshdlr) );
4973  consdata->issimplified = TRUE;
4974 
4975  if( changed )
4976  havechange = TRUE;
4977 
4978  /* If root expression changed, then we need to take care updating the locks as well (the consdata is the one holding consdata->expr "as a child").
4979  * If root expression did not change, some subexpression may still have changed, but the locks were taking care of in the corresponding SCIPreplaceExprChild() call.
4980  */
4981  if( simplified != consdata->expr )
4982  {
4983  assert(changed);
4984 
4985  /* release old expression */
4986  SCIP_CALL( SCIPreleaseExpr(scip, &consdata->expr) );
4987 
4988  /* store simplified expression */
4989  consdata->expr = simplified;
4990  }
4991  else
4992  {
4993  /* The simplify captures simplified in any case, also if nothing has changed.
4994  * Therefore, we have to release it here.
4995  */
4996  SCIP_CALL( SCIPreleaseExpr(scip, &simplified) );
4997  }
4998 
4999  if( *infeasible )
5000  break;
5001 
5002  /* scale constraint sides */
5003  SCIP_CALL( scaleConsSides(scip, conshdlr, conss[i], &changed) );
5004 
5005  if( changed )
5006  havechange = TRUE;
5007 
5008  /* handle constant root expression; either the problem is infeasible or the constraint is redundant */
5009  if( SCIPisExprValue(scip, consdata->expr) )
5010  {
5011  SCIP_Real value = SCIPgetValueExprValue(consdata->expr);
5012  if( (!SCIPisInfinity(scip, -consdata->lhs) && SCIPisFeasNegative(scip, value - consdata->lhs)) ||
5013  (!SCIPisInfinity(scip, consdata->rhs) && SCIPisFeasPositive(scip, value - consdata->rhs)) )
5014  {
5015  SCIPdebugMsg(scip, "<%s> with constant expression found infeasible\n", SCIPconsGetName(conss[i]));
5016  SCIPdebugPrintCons(scip, conss[i], NULL);
5017  *infeasible = TRUE;
5018  break;
5019  }
5020  else
5021  {
5022  SCIP_CALL( addLocks(scip, conss[i], nlockspos[i], nlocksneg[i]) );
5023  SCIP_CALL( SCIPdelCons(scip, conss[i]) );
5024  if( ndelconss != NULL )
5025  ++*ndelconss;
5026  havechange = TRUE;
5027  }
5028  }
5029  }
5030  }
5031 
5032  /* replace common subexpressions */
5033  if( havechange && !*infeasible )
5034  {
5035  SCIP_CONS** consssorted;
5036  SCIP_EXPR** rootexprs;
5037  SCIP_Bool replacedroot;
5038 
5039  SCIP_CALL( SCIPallocBufferArray(scip, &rootexprs, nconss) );
5040  for( i = 0; i < nconss; ++i )
5041  rootexprs[i] = SCIPconsGetData(conss[i])->expr;
5042 
5043  SCIP_CALL( SCIPreplaceCommonSubexpressions(scip, rootexprs, nconss, &replacedroot) );
5044 
5045  /* update pointer to root expr in constraints, if any has changed
5046  * SCIPreplaceCommonSubexpressions will have released the old expr and captures the new one
5047  */
5048  if( replacedroot )
5049  for( i = 0; i < nconss; ++i )
5050  SCIPconsGetData(conss[i])->expr = rootexprs[i];
5051 
5052  SCIPfreeBufferArray(scip, &rootexprs);
5053 
5054  /* TODO this is a possibly expensive way to update the variable expressions stored inside an expression which might have
5055  * been changed after simplification; now we completely recollect all variable expression and variable events
5056  */
5057 
5058  /* Each variable stores the constraints for which it catched varbound events sorted by the constraint index.
5059  * Thus, for performance reasons, it is better to call dropVarEvents in descending order of constraint index.
5060  */
5061  SCIP_CALL( SCIPduplicateBufferArray(scip, &consssorted, conss, nconss) );
5062  SCIPsortPtr((void**)consssorted, compIndexConsNonlinear, nconss);
5063 
5064  for( i = nconss-1; i >= 0; --i )
5065  {
5066  assert(i == 0 || compIndexConsNonlinear((void*)consssorted[i-1], (void*)consssorted[i]) < 0);
5067  if( SCIPconsIsDeleted(consssorted[i]) )
5068  continue;
5069 
5070  SCIP_CALL( dropVarEvents(scip, conshdlrdata->eventhdlr, consssorted[i]) );
5071  SCIP_CALL( freeVarExprs(scip, SCIPconsGetData(consssorted[i])) );
5072  }
5073  for( i = 0; i < nconss; ++i )
5074  {
5075  if( SCIPconsIsDeleted(consssorted[i]) )
5076  continue;
5077 
5078  SCIP_CALL( storeVarExprs(scip, conshdlr, SCIPconsGetData(consssorted[i])) );
5079  SCIP_CALL( catchVarEvents(scip, conshdlrdata->eventhdlr, consssorted[i]) );
5080  }
5081 
5082  SCIPfreeBufferArray(scip, &consssorted);
5083 
5084  /* forbid multiaggregation for nonlinear variables again (in case new variables appeared now)
5085  * a multiaggregation of a nonlinear variable can yield to a large increase in expressions due to
5086  * expanding terms in simplify, e.g. ,(sum_i x_i)^2, so we just forbid these
5087  */
5088  SCIP_CALL( forbidNonlinearVariablesMultiaggration(scip, conshdlr, conss, nconss) );
5089  }
5090 
5091  /* restore locks */
5092  for( i = 0; i < nconss; ++i )
5093  {
5094  if( SCIPconsIsDeleted(conss[i]) )
5095  continue;
5096 
5097  SCIP_CALL( addLocks(scip, conss[i], nlockspos[i], nlocksneg[i]) );
5098  }
5099 
5100  /* run nlhdlr detect if in presolving stage (that is, not in exitpre)
5101  * TODO can we skip this in presoltiming fast?
5102  */
5103  if( SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING && !*infeasible )
5104  {
5105  /* reset one of the number of detections counter to count only current presolving round */
5106  for( i = 0; i < conshdlrdata->nnlhdlrs; ++i )
5107  SCIPnlhdlrResetNDetectionslast(conshdlrdata->nlhdlrs[i]);
5108 
5109  SCIP_CALL( initSolve(scip, conshdlr, conss, nconss) );
5110  }
5111 
5112  /* free allocated memory */
5113  SCIPfreeBufferArray(scip, &nlocksneg);
5114  SCIPfreeBufferArray(scip, &nlockspos);
5115 
5116  SCIP_CALL( SCIPstopClock(scip, conshdlrdata->canonicalizetime) );
5117 
5118  return SCIP_OKAY;
5119 }
5120 
5121 /** merges constraints that have the same root expression */
5122 static
5124  SCIP* scip, /**< SCIP data structure */
5125  SCIP_CONS** conss, /**< constraints to process */
5126  int nconss, /**< number of constraints */
5127  SCIP_Bool* success /**< pointer to store whether at least one constraint could be deleted */
5128  )
5129 {
5130  SCIP_HASHMAP* expr2cons;
5131  SCIP_Bool* updatelocks;
5132  int* nlockspos;
5133  int* nlocksneg;
5134  int c;
5135 
5136  assert(success != NULL);
5137 
5138  *success = FALSE;
5139 
5140  /* not enough constraints available */
5141  if( nconss <= 1 )
5142  return SCIP_OKAY;
5143 
5144  SCIP_CALL( SCIPhashmapCreate(&expr2cons, SCIPblkmem(scip), nconss) );
5145  SCIP_CALL( SCIPallocClearBufferArray(scip, &updatelocks, nconss) );
5146  SCIP_CALL( SCIPallocBufferArray(scip, &nlockspos, nconss) );
5147  SCIP_CALL( SCIPallocBufferArray(scip, &nlocksneg, nconss) );
5148 
5149  for( c = 0; c < nconss; ++c )
5150  {
5151  SCIP_CONSDATA* consdata;
5152 
5153  /* ignore deleted constraints */
5154  if( SCIPconsIsDeleted(conss[c]) )
5155  continue;
5156 
5157  consdata = SCIPconsGetData(conss[c]);
5158  assert(consdata != NULL);
5159 
5160  /* add expression to the hash map if not seen so far */
5161  if( !SCIPhashmapExists(expr2cons, (void*)consdata->expr) )
5162  {
5163  SCIP_CALL( SCIPhashmapInsertInt(expr2cons, (void*)consdata->expr, c) );
5164  }
5165  else
5166  {
5167  SCIP_CONSDATA* imgconsdata;
5168  int idx;
5169 
5170  idx = SCIPhashmapGetImageInt(expr2cons, (void*)consdata->expr);
5171  assert(idx >= 0 && idx < nconss);
5172 
5173  imgconsdata = SCIPconsGetData(conss[idx]);
5174  assert(imgconsdata != NULL);
5175  assert(imgconsdata->expr == consdata->expr);
5176 
5177  SCIPdebugMsg(scip, "merge constraint %g <= %s <= %g with %g <= %s <= %g\n", consdata->lhs,
5178  SCIPconsGetName(conss[c]), consdata->rhs, imgconsdata->lhs, SCIPconsGetName(conss[idx]), imgconsdata->rhs);
5179 
5180  /* check whether locks need to be updated */
5181  if( !updatelocks[idx] && ((SCIPisInfinity(scip, -imgconsdata->lhs) && !SCIPisInfinity(scip, -consdata->lhs))
5182  || (SCIPisInfinity(scip, imgconsdata->rhs) && !SCIPisInfinity(scip, consdata->rhs))) )
5183  {
5184  nlockspos[idx] = imgconsdata->nlockspos;
5185  nlocksneg[idx] = imgconsdata->nlocksneg;
5186  SCIP_CALL( addLocks(scip, conss[idx], -imgconsdata->nlockspos, -imgconsdata->nlocksneg) );
5187  updatelocks[idx] = TRUE;
5188  }
5189 
5190  /* update constraint sides */
5191  imgconsdata->lhs = MAX(imgconsdata->lhs, consdata->lhs);
5192  imgconsdata->rhs = MIN(imgconsdata->rhs, consdata->rhs);
5193 
5194  /* delete constraint */
5195  SCIP_CALL( SCIPdelCons(scip, conss[c]) );
5196  *success = TRUE;
5197  }
5198  }
5199 
5200  /* restore locks of updated constraints */
5201  if( *success )
5202  {
5203  for( c = 0; c < nconss; ++c )
5204  {
5205  if( updatelocks[c] )
5206  {
5207  SCIP_CALL( addLocks(scip, conss[c], nlockspos[c], nlocksneg[c]) );
5208  }
5209  }
5210  }
5211 
5212  /* free memory */
5213  SCIPfreeBufferArray(scip, &nlocksneg);
5214  SCIPfreeBufferArray(scip, &nlockspos);
5215  SCIPfreeBufferArray(scip, &updatelocks);
5216  SCIPhashmapFree(&expr2cons);
5217 
5218  return SCIP_OKAY;
5219 }
5220 
5221 /** interval evaluation of variables as used in redundancy check
5222  *
5223  * Returns local variable bounds of a variable, relaxed by feastol, as interval.
5224  */
5225 static
5226 SCIP_DECL_EXPR_INTEVALVAR(intEvalVarRedundancyCheck)
5227 { /*lint --e{715}*/
5228  SCIP_CONSHDLRDATA* conshdlrdata;
5229  SCIP_INTERVAL interval;
5230  SCIP_Real lb;
5231  SCIP_Real ub;
5232 
5233  assert(scip != NULL);
5234  assert(var != NULL);
5235 
5236  conshdlrdata = (SCIP_CONSHDLRDATA*)intevalvardata;
5237  assert(conshdlrdata != NULL);
5238 
5239  if( conshdlrdata->globalbounds )
5240  {
5241  lb = SCIPvarGetLbGlobal(var);
5242  ub = SCIPvarGetUbGlobal(var);
5243  }
5244  else
5245  {
5246  lb = SCIPvarGetLbLocal(var);
5247  ub = SCIPvarGetUbLocal(var);
5248  }
5249  assert(lb <= ub); /* can SCIP ensure by now that variable bounds are not contradicting? */
5250 
5251  /* relax variable bounds, if there are bounds and variable is not fixed
5252  * (actually some assert complains if trying SCIPisRelEQ if both bounds are at different infinity)
5253  */
5254  if( !(SCIPisInfinity(scip, -lb) && SCIPisInfinity(scip, ub)) && !SCIPisRelEQ(scip, lb, ub) )
5255  {
5256  if( !SCIPisInfinity(scip, -lb) )
5257  lb -= SCIPfeastol(scip);
5258 
5259  if( !SCIPisInfinity(scip, ub) )
5260  ub += SCIPfeastol(scip);
5261  }
5262 
5263  /* convert SCIPinfinity() to SCIP_INTERVAL_INFINITY */
5264  lb = -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -lb);
5266  assert(lb <= ub);
5267 
5268  SCIPintervalSetBounds(&interval, lb, ub);
5269 
5270  return interval;
5271 }
5272 
5273 /** removes constraints that are always feasible or very simple
5274  *
5275  * Checks whether the activity of constraint functions is a subset of the constraint sides (relaxed by feastol).
5276  * To compute the activity, we use forwardPropExpr(), but relax variable bounds by feastol, because solutions to be checked
5277  * might violate variable bounds by up to feastol, too.
5278  * This is the main reason why the redundancy check is not done in propConss(), which relaxes variable bounds by epsilon only.
5279  *
5280  * Also removes constraints of the form lhs &le; variable &le; rhs.
5281  *
5282  * @todo it would be sufficient to check constraints for which we know that they are not currently violated by a valid solution
5283  *
5284  * @note This could should not run during solving, because the forwardProp takes the bounds of auxiliary variables into account.
5285  * For the root expression, these bounds are already set to the constraint sides, so that the activity of every expression
5286  * would appear as if the constraint is redundant.
5287  */
5288 static
5290  SCIP* scip, /**< SCIP data structure */
5291  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
5292  SCIP_CONS** conss, /**< constraints to propagate */
5293  int nconss, /**< total number of constraints */
5294  SCIP_Bool* cutoff, /**< pointer to store whether infeasibility has been identified */
5295  int* ndelconss, /**< buffer to add the number of deleted constraints */
5296  int* nchgbds /**< buffer to add the number of variable bound tightenings */
5297  )
5298 {
5299  SCIP_CONSHDLRDATA* conshdlrdata;
5300  SCIP_CONSDATA* consdata;
5301  SCIP_INTERVAL activity;
5302  SCIP_INTERVAL sides;
5303  int i;
5304 
5305  assert(scip != NULL);
5306  assert(conshdlr != NULL);
5307  assert(conss != NULL);
5308  assert(nconss >= 0);
5309  assert(cutoff != NULL);
5310  assert(ndelconss != NULL);
5311  assert(nchgbds != NULL);
5312 
5313  /* no constraints to check */
5314  if( nconss == 0 )
5315  return SCIP_OKAY;
5316 
5317  conshdlrdata = SCIPconshdlrGetData(conshdlr);
5318  assert(conshdlrdata != NULL);
5319 
5320  /* increase curboundstag and set lastvaractivitymethodchange
5321  * we do this here to trigger a reevaluation of all variable bounds, since we will relax variable bounds
5322  * for the redundancy check differently than for domain propagation
5323  * we also update lastboundrelax to ensure activites of all expressions are indeed reevaluated
5324  */
5325  ++conshdlrdata->curboundstag;
5326  assert(conshdlrdata->curboundstag > 0);
5327  conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
5328  conshdlrdata->lastboundrelax = conshdlrdata->curboundstag;
5329  conshdlrdata->intevalvar = intEvalVarRedundancyCheck;
5330 
5331  SCIPdebugMsg(scip, "checking %d constraints for redundancy\n", nconss);
5332 
5333  *cutoff = FALSE;
5334  for( i = 0; i < nconss; ++i )
5335  {
5336  if( !SCIPconsIsActive(conss[i]) || SCIPconsIsDeleted(conss[i]) )
5337  continue;
5338 
5339  consdata = SCIPconsGetData(conss[i]);
5340  assert(consdata != NULL);
5341 
5342  /* handle constant expressions separately: either the problem is infeasible or the constraint is redundant */
5343  if( SCIPisExprValue(scip, consdata->expr) )
5344  {
5345  SCIP_Real value = SCIPgetValueExprValue(consdata->expr);
5346 
5347  if( (!SCIPisInfinity(scip, -consdata->lhs) && value < consdata->lhs - SCIPfeastol(scip)) ||
5348  (!SCIPisInfinity(scip, consdata->rhs) && value > consdata->rhs + SCIPfeastol(scip)) )
5349  {
5350  SCIPdebugMsg(scip, "constant constraint <%s> is infeasible: %g in [%g,%g] ", SCIPconsGetName(conss[i]), value, consdata->lhs, consdata->rhs);
5351  *cutoff = TRUE;
5352 
5353  goto TERMINATE;
5354  }
5355 
5356  SCIPdebugMsg(scip, "constant constraint <%s> is redundant: %g in [%g,%g] ", SCIPconsGetName(conss[i]), value, consdata->lhs, consdata->rhs);
5357 
5358  SCIP_CALL( SCIPdelConsLocal(scip, conss[i]) );
5359  ++*ndelconss;
5360 
5361  continue;
5362  }
5363 
5364  /* handle variable expressions separately: tighten variable bounds to constraint sides, then remove constraint (now redundant) */
5365  if( SCIPisExprVar(scip, consdata->expr) )
5366  {
5367  SCIP_VAR* var;
5368  SCIP_Bool tightened;
5369 
5370  var = SCIPgetVarExprVar(consdata->expr);
5371  assert(var != NULL);
5372 
5373  SCIPdebugMsg(scip, "variable constraint <%s> can be made redundant: <%s>[%g,%g] in [%g,%g]\n", SCIPconsGetName(conss[i]), SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), consdata->lhs, consdata->rhs);
5374 
5375  /* ensure that variable bounds are within constraint sides */
5376  if( !SCIPisInfinity(scip, -consdata->lhs) )
5377  {
5378  SCIP_CALL( SCIPtightenVarLb(scip, var, consdata->lhs, TRUE, cutoff, &tightened) );
5379 
5380  if( tightened )
5381  ++*nchgbds;
5382 
5383  if( *cutoff )
5384  goto TERMINATE;
5385  }
5386 
5387  if( !SCIPisInfinity(scip, consdata->rhs) )
5388  {
5389  SCIP_CALL( SCIPtightenVarUb(scip, var, consdata->rhs, TRUE, cutoff, &tightened) );
5390 
5391  if( tightened )
5392  ++*nchgbds;
5393 
5394  if( *cutoff )
5395  goto TERMINATE;
5396  }
5397 
5398  /* delete the (now) redundant constraint locally */
5399  SCIP_CALL( SCIPdelConsLocal(scip, conss[i]) );
5400  ++*ndelconss;
5401 
5402  continue;
5403  }
5404 
5405  /* reevaluate expression activity, now using intEvalVarRedundancyCheck
5406  * we relax variable bounds by feastol here, as solutions that are checked later can also violate
5407  * variable bounds by up to feastol
5408  * (relaxing fixed variables seems to be too much, but they would be removed by presolve soon anyway)
5409  */
5410  SCIPdebugMsg(scip, "call forwardPropExpr() for constraint <%s>: ", SCIPconsGetName(conss[i]));
5411  SCIPdebugPrintCons(scip, conss[i], NULL);
5412 
5413  SCIP_CALL( forwardPropExpr(scip, conshdlr, consdata->expr, FALSE, cutoff, NULL) );
5414  assert(*cutoff || !SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(consdata->expr)));
5415 
5416  /* it is unlikely that we detect infeasibility by doing forward propagation */
5417  if( *cutoff )
5418  {
5419  SCIPdebugMsg(scip, " -> cutoff\n");
5420  goto TERMINATE;
5421  }
5422 
5423  assert(SCIPexprGetActivityTag(consdata->expr) == conshdlrdata->curboundstag);
5424  activity = SCIPexprGetActivity(consdata->expr);
5425 
5426  /* relax sides by feastol
5427  * we could accept every solution that violates constraints up to feastol as redundant, so this is the most permissive we can be
5428  */
5429  SCIPintervalSetBounds(&sides,
5430  SCIPisInfinity(scip, -consdata->lhs) ? -SCIP_INTERVAL_INFINITY : consdata->lhs - SCIPfeastol(scip),
5431  SCIPisInfinity(scip, consdata->rhs) ? SCIP_INTERVAL_INFINITY : consdata->rhs + SCIPfeastol(scip));
5432 
5433  if( SCIPintervalIsSubsetEQ(SCIP_INTERVAL_INFINITY, activity, sides) )
5434  {
5435  SCIPdebugMsg(scip, " -> redundant: activity [%g,%g] within sides [%g,%g]\n", activity.inf, activity.sup, consdata->lhs, consdata->rhs);
5436 
5437  SCIP_CALL( SCIPdelConsLocal(scip, conss[i]) );
5438  ++*ndelconss;
5439 
5440  continue;
5441  }
5442 
5443  SCIPdebugMsg(scip, " -> not redundant: activity [%g,%g] not within sides [%g,%g]\n", activity.inf, activity.sup, consdata->lhs, consdata->rhs);
5444  }
5445 
5446 TERMINATE:
5447  /* make sure all activities are reevaluated again, since we relaxed bounds in a different way */
5448  ++conshdlrdata->curboundstag;
5449  conshdlrdata->lastvaractivitymethodchange = conshdlrdata->curboundstag;
5450  conshdlrdata->lastboundrelax = conshdlrdata->curboundstag;
5451  conshdlrdata->intevalvar = intEvalVarBoundTightening;
5452 
5453  return SCIP_OKAY;
5454 }
5455 
5456 /** tries to automatically convert a nonlinear constraint into a more specific and more specialized constraint */
5457 static
5459  SCIP* scip, /**< SCIP data structure */
5460  SCIP_CONSHDLR* conshdlr, /**< constraint handler data structure */
5461  SCIP_CONS* cons, /**< source constraint to try to convert */
5462  SCIP_Bool* upgraded, /**< buffer to store whether constraint was upgraded */
5463  int* nupgdconss, /**< buffer to increase if constraint was upgraded */
5464  int* naddconss /**< buffer to increase with number of additional constraints created during upgrade */
5465  )
5466 {
5467  SCIP_CONSHDLRDATA* conshdlrdata;
5468  SCIP_CONSDATA* consdata;
5469  SCIP_CONS** upgdconss;
5470  int upgdconsssize;
5471  int nupgdconss_;
5472  int i;
5473 
5474  assert(scip != NULL);
5475  assert(conshdlr != NULL);
5476  assert(cons != NULL);
5477  assert(!SCIPconsIsModifiable(cons));
5478  assert(upgraded != NULL);
5479  assert(nupgdconss != NULL);
5480  assert(naddconss != NULL);
5481 
5482  *upgraded = FALSE;
5483 
5484  nupgdconss_ = 0;
5485 
5486  conshdlrdata = SCIPconshdlrGetData(conshdlr);
5487  assert(conshdlrdata != NULL);
5488 
5489  /* if there are no upgrade methods, we can stop */
5490  if( conshdlrdata->nconsupgrades == 0 )
5491  return SCIP_OKAY;
5492 
5493  upgdconsssize = 2;
5494  SCIP_CALL( SCIPallocBufferArray(scip, &upgdconss, upgdconsssize) );
5495 
5496  /* call the upgrading methods */
5497  SCIPdebugMsg(scip, "upgrading nonlinear constraint <%s> (up to %d upgrade methods): ", SCIPconsGetName(cons), conshdlrdata->nconsupgrades);
5498  SCIPdebugPrintCons(scip, cons, NULL);
5499 
5500  consdata = SCIPconsGetData(cons);
5501  assert(consdata != NULL);
5502 
5503  /* try all upgrading methods in priority order in case the upgrading step is enable */
5504  for( i = 0; i < conshdlrdata->nconsupgrades; ++i )
5505  {
5506  if( !conshdlrdata->consupgrades[i]->active )
5507  continue;
5508 
5509  assert(conshdlrdata->consupgrades[i]->consupgd != NULL);
5510 
5511  SCIP_CALL( conshdlrdata->consupgrades[i]->consupgd(scip, cons, consdata->nvarexprs, &nupgdconss_, upgdconss, upgdconsssize) );
5512 
5513  while( nupgdconss_ < 0 )
5514  {
5515  /* upgrade function requires more memory: resize upgdconss and call again */
5516  assert(-nupgdconss_ > upgdconsssize);
5517  upgdconsssize = -nupgdconss_;
5518  SCIP_CALL( SCIPreallocBufferArray(scip, &upgdconss, -nupgdconss_) );
5519 
5520  SCIP_CALL( conshdlrdata->consupgrades[i]->consupgd(scip, cons, consdata->nvarexprs, &nupgdconss_, upgdconss, upgdconsssize) );
5521 
5522  assert(nupgdconss_ != 0);
5523  }
5524 
5525  if( nupgdconss_ > 0 )
5526  {
5527  /* got upgrade */
5528  int j;
5529 
5530  SCIPdebugMsg(scip, " -> upgraded to %d constraints:\n", nupgdconss_);
5531 
5532  /* add the upgraded constraints to the problem and forget them */
5533  for( j = 0; j < nupgdconss_; ++j )
5534  {
5535  SCIPdebugMsgPrint(scip, "\t");
5536  SCIPdebugPrintCons(scip, upgdconss[j], NULL);
5537 
5538  SCIP_CALL( SCIPaddCons(scip, upgdconss[j]) ); /*lint !e613*/
5539  SCIP_CALL( SCIPreleaseCons(scip, &upgdconss[j]) ); /*lint !e613*/
5540  }
5541 
5542  /* count the first upgrade constraint as constraint upgrade and the remaining ones as added constraints */
5543  *nupgdconss += 1;
5544  *naddconss += nupgdconss_ - 1;
5545  *upgraded = TRUE;
5546 
5547  /* delete upgraded constraint */
5548  SCIPdebugMsg(scip, "delete constraint <%s> after upgrade\n", SCIPconsGetName(cons));
5549  SCIP_CALL( SCIPdelCons(scip, cons) );
5550 
5551  break;
5552  }
5553  }
5554 
5555  SCIPfreeBufferArray(scip, &upgdconss);
5556 
5557  return SCIP_OKAY;
5558 }
5559 
5560 /** returns whether the variable of a given variable expression is a candidate for presolveSingleLockedVars(), i.e.,
5561  * the variable is only contained in a single nonlinear constraint, has no objective coefficient, has finite
5562  * variable bounds, and is not binary
5563  */
5564 static
5566  SCIP* scip, /**< SCIP data structure */
5567  SCIP_EXPR* expr /**< variable expression */
5568  )
5569 {
5570  SCIP_VAR* var;
5571  SCIP_EXPR_OWNERDATA* ownerdata;
5572 
5573  assert(SCIPisExprVar(scip, expr));
5574 
5575  var = SCIPgetVarExprVar(expr);
5576  assert(var != NULL);
5577 
5578  ownerdata = SCIPexprGetOwnerData(expr);
5579  assert(ownerdata != NULL);
5580 
5581  return SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == ownerdata->nlocksneg
5582  && SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == ownerdata->nlockspos
5583  && ownerdata->nconss == 1 && SCIPisZero(scip, SCIPvarGetObj(var))
5584  && !SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)) && !SCIPisInfinity(scip, SCIPvarGetUbGlobal(var))
5586  && !SCIPisEQ(scip, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
5587 }
5588 
5589 /** removes all variable expressions that are contained in a given expression from a hash map */
5590 static
5592  SCIP* scip, /**< SCIP data structure */
5593  SCIP_EXPR* expr, /**< expression */
5594  SCIP_EXPRITER* it, /**< expression iterator */
5595  SCIP_HASHMAP* exprcands /**< map to hash variable expressions */
5596  )
5597 {
5598  SCIP_EXPR* e;
5599 
5600  for( e = SCIPexpriterRestartDFS(it, expr); !SCIPexpriterIsEnd(it); e = SCIPexpriterGetNext(it) )
5601  {
5602  if( SCIPisExprVar(scip, e) && SCIPhashmapExists(exprcands, (void*)e) )
5603  {
5604  SCIP_CALL( SCIPhashmapRemove(exprcands, (void*)e) );
5605  }
5606  }
5607 
5608  return SCIP_OKAY;
5609 }
5610 
5611 /** presolving method to fix a variable \f$x_i\f$ to one of its bounds if the variable is only contained in a single
5612  * nonlinear constraint g(x) &le; rhs (&ge; lhs) if g() is concave (convex) in \f$x_i\f$
5613  *
5614  * If a continuous variable has bounds [0,1], then the variable type is changed to be binary.
5615  * Otherwise, a bound disjunction constraint is added.
5616  *
5617  * @todo the same reduction can be applied if g(x) is not concave, but monotone in \f$x_i\f$ for g(x) &le; rhs
5618  * @todo extend this to cases where a variable can appear in a monomial with an exponent, essentially relax
5619  * g(x) to \f$\sum_i [a_i,b_i] x^{p_i}\f$ for a single variable \f$x\f$ and try to conclude montonicity or convexity/concavity
5620  * on this (probably have one or two flags per variable and update this whenever another \f$x^{p_i}\f$ is found)
5621  */
5622 static
5624  SCIP* scip, /**< SCIP data structure */
5625  SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */
5626  SCIP_CONS* cons, /**< nonlinear c