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