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