Scippy

SCIP

Solving Constraint Integer Programs

scip_nonlinear.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2020 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file scip_nonlinear.c
17  * @ingroup OTHER_CFILES
18  * @brief public methods for nonlinear functions
19  * @author Tobias Achterberg
20  * @author Timo Berthold
21  * @author Gerald Gamrath
22  * @author Leona Gottwald
23  * @author Stefan Heinz
24  * @author Gregor Hendel
25  * @author Thorsten Koch
26  * @author Alexander Martin
27  * @author Marc Pfetsch
28  * @author Michael Winkler
29  * @author Kati Wolter
30  *
31  * @todo check all SCIP_STAGE_* switches, and include the new stages TRANSFORMED and INITSOLVE
32  */
33 
34 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35 
36 #include "blockmemshell/memory.h"
37 #include "nlpi/nlpi.h"
38 #include "nlpi/pub_expr.h"
39 #include "nlpi/type_expr.h"
40 #include "scip/dbldblarith.h"
41 #include "scip/pub_lp.h"
42 #include "scip/pub_message.h"
43 #include "scip/pub_misc.h"
44 #include "scip/pub_nlp.h"
45 #include "scip/pub_var.h"
46 #include "scip/scip_mem.h"
47 #include "scip/scip_message.h"
48 #include "scip/scip_nonlinear.h"
49 #include "scip/scip_numerics.h"
50 #include "scip/scip_prob.h"
51 
52 /** computes coefficients of linearization of a square term in a reference point */
54  SCIP* scip, /**< SCIP data structure */
55  SCIP_Real sqrcoef, /**< coefficient of square term */
56  SCIP_Real refpoint, /**< point where to linearize */
57  SCIP_Bool isint, /**< whether corresponding variable is a discrete variable, and thus linearization could be moved */
58  SCIP_Real* lincoef, /**< buffer to add coefficient of linearization */
59  SCIP_Real* linconstant, /**< buffer to add constant of linearization */
60  SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
61  )
62 {
63  assert(scip != NULL);
64  assert(lincoef != NULL);
65  assert(linconstant != NULL);
66  assert(success != NULL);
67 
68  if( sqrcoef == 0.0 )
69  return;
70 
71  if( SCIPisInfinity(scip, REALABS(refpoint)) )
72  {
73  *success = FALSE;
74  return;
75  }
76 
77  if( !isint || SCIPisIntegral(scip, refpoint) )
78  {
79  SCIP_Real tmp;
80 
81  /* sqrcoef * x^2 -> tangent in refpoint = sqrcoef * 2 * refpoint * (x - refpoint) */
82 
83  tmp = sqrcoef * refpoint;
84 
85  if( SCIPisInfinity(scip, 2.0 * REALABS(tmp)) )
86  {
87  *success = FALSE;
88  return;
89  }
90 
91  *lincoef += 2.0 * tmp;
92  tmp *= refpoint;
93  *linconstant -= tmp;
94  }
95  else
96  {
97  /* sqrcoef * x^2 -> secant between f=floor(refpoint) and f+1 = sqrcoef * (f^2 + ((f+1)^2 - f^2) * (x-f))
98  * = sqrcoef * (-f*(f+1) + (2*f+1)*x)
99  */
100  SCIP_Real f;
101  SCIP_Real coef;
102  SCIP_Real constant;
103 
104  f = SCIPfloor(scip, refpoint);
105 
106  coef = sqrcoef * (2.0 * f + 1.0);
107  constant = -sqrcoef * f * (f + 1.0);
108 
109  if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
110  {
111  *success = FALSE;
112  return;
113  }
114 
115  *lincoef += coef;
116  *linconstant += constant;
117  }
118 }
119 
120 /** computes coefficients of secant of a square term */
122  SCIP* scip, /**< SCIP data structure */
123  SCIP_Real sqrcoef, /**< coefficient of square term */
124  SCIP_Real lb, /**< lower bound on variable */
125  SCIP_Real ub, /**< upper bound on variable */
126  SCIP_Real refpoint, /**< point for which to compute value of linearization */
127  SCIP_Real* lincoef, /**< buffer to add coefficient of secant */
128  SCIP_Real* linconstant, /**< buffer to add constant of secant */
129  SCIP_Bool* success /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */
130  )
131 {
132  SCIP_Real coef;
133  SCIP_Real constant;
134 
135  assert(scip != NULL);
136  assert(!SCIPisInfinity(scip, lb));
137  assert(!SCIPisInfinity(scip, -ub));
138  assert(SCIPisLE(scip, lb, ub));
139  assert(SCIPisLE(scip, lb, refpoint));
140  assert(SCIPisGE(scip, ub, refpoint));
141  assert(lincoef != NULL);
142  assert(linconstant != NULL);
143  assert(success != NULL);
144 
145  if( sqrcoef == 0.0 )
146  return;
147 
148  if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) )
149  {
150  /* unboundedness */
151  *success = FALSE;
152  return;
153  }
154 
155  /* sqrcoef * x^2 -> sqrcoef * (lb * lb + (ub*ub - lb*lb)/(ub-lb) * (x-lb)) = sqrcoef * (lb*lb + (ub+lb)*(x-lb))
156  * = sqrcoef * ((lb+ub)*x - lb*ub)
157  */
158  coef = sqrcoef * (lb + ub);
159  constant = -sqrcoef * lb * ub;
160  if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
161  {
162  *success = FALSE;
163  return;
164  }
165 
166  *lincoef += coef;
167  *linconstant += constant;
168 }
169 
170 /** computes coefficients of linearization of a bilinear term in a reference point */
172  SCIP* scip, /**< SCIP data structure */
173  SCIP_Real bilincoef, /**< coefficient of bilinear term */
174  SCIP_Real refpointx, /**< point where to linearize first variable */
175  SCIP_Real refpointy, /**< point where to linearize second variable */
176  SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
177  SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
178  SCIP_Real* linconstant, /**< buffer to add constant of linearization */
179  SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
180  )
181 {
182  SCIP_Real constant;
183 
184  assert(scip != NULL);
185  assert(lincoefx != NULL);
186  assert(lincoefy != NULL);
187  assert(linconstant != NULL);
188  assert(success != NULL);
189 
190  if( bilincoef == 0.0 )
191  return;
192 
193  if( SCIPisInfinity(scip, REALABS(refpointx)) || SCIPisInfinity(scip, REALABS(refpointy)) )
194  {
195  *success = FALSE;
196  return;
197  }
198 
199  /* bilincoef * x * y -> bilincoef * (refpointx * refpointy + refpointy * (x - refpointx) + refpointx * (y - refpointy))
200  * = -bilincoef * refpointx * refpointy + bilincoef * refpointy * x + bilincoef * refpointx * y
201  */
202 
203  constant = -bilincoef * refpointx * refpointy;
204 
205  if( SCIPisInfinity(scip, REALABS(bilincoef * refpointx)) || SCIPisInfinity(scip, REALABS(bilincoef * refpointy))
206  || SCIPisInfinity(scip, REALABS(constant)) )
207  {
208  *success = FALSE;
209  return;
210  }
211 
212  *lincoefx += bilincoef * refpointy;
213  *lincoefy += bilincoef * refpointx;
214  *linconstant += constant;
215 }
216 
217 /** computes coefficients of McCormick under- or overestimation of a bilinear term */
219  SCIP* scip, /**< SCIP data structure */
220  SCIP_Real bilincoef, /**< coefficient of bilinear term */
221  SCIP_Real lbx, /**< lower bound on first variable */
222  SCIP_Real ubx, /**< upper bound on first variable */
223  SCIP_Real refpointx, /**< reference point for first variable */
224  SCIP_Real lby, /**< lower bound on second variable */
225  SCIP_Real uby, /**< upper bound on second variable */
226  SCIP_Real refpointy, /**< reference point for second variable */
227  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
228  SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
229  SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
230  SCIP_Real* linconstant, /**< buffer to add constant of linearization */
231  SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
232  )
233 {
234  SCIP_Real constant;
235  SCIP_Real coefx;
236  SCIP_Real coefy;
237 
238  assert(scip != NULL);
239  assert(!SCIPisInfinity(scip, lbx));
240  assert(!SCIPisInfinity(scip, -ubx));
241  assert(!SCIPisInfinity(scip, lby));
242  assert(!SCIPisInfinity(scip, -uby));
243  assert(SCIPisLE(scip, lbx, ubx));
244  assert(SCIPisLE(scip, lby, uby));
245  assert(SCIPisLE(scip, lbx, refpointx));
246  assert(SCIPisGE(scip, ubx, refpointx));
247  assert(SCIPisLE(scip, lby, refpointy));
248  assert(SCIPisGE(scip, uby, refpointy));
249  assert(lincoefx != NULL);
250  assert(lincoefy != NULL);
251  assert(linconstant != NULL);
252  assert(success != NULL);
253 
254  if( bilincoef == 0.0 )
255  return;
256 
257  if( overestimate )
258  bilincoef = -bilincoef;
259 
260  if( SCIPisRelEQ(scip, lbx, ubx) && SCIPisRelEQ(scip, lby, uby) )
261  {
262  /* both x and y are mostly fixed */
263  SCIP_Real cand1;
264  SCIP_Real cand2;
265  SCIP_Real cand3;
266  SCIP_Real cand4;
267 
268  coefx = 0.0;
269  coefy = 0.0;
270 
271  /* estimate x * y by constant */
272  cand1 = lbx * lby;
273  cand2 = lbx * uby;
274  cand3 = ubx * lby;
275  cand4 = ubx * uby;
276 
277  /* take most conservative value for underestimator */
278  if( bilincoef < 0.0 )
279  constant = bilincoef * MAX( MAX(cand1, cand2), MAX(cand3, cand4) );
280  else
281  constant = bilincoef * MIN( MIN(cand1, cand2), MIN(cand3, cand4) );
282  }
283  else if( bilincoef > 0.0 )
284  {
285  /* either x or y is not fixed and coef > 0.0 */
286  if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, -lby) &&
287  (SCIPisInfinity(scip, ubx) || SCIPisInfinity(scip, uby)
288  || (uby - refpointy) * (ubx - refpointx) >= (refpointy - lby) * (refpointx - lbx)) )
289  {
290  if( SCIPisRelEQ(scip, lbx, ubx) )
291  {
292  /* x*y = lbx * y + (x-lbx) * y >= lbx * y + (x-lbx) * lby >= lbx * y + min{(ubx-lbx) * lby, 0 * lby} */
293  coefx = 0.0;
294  coefy = bilincoef * lbx;
295  constant = bilincoef * (lby < 0.0 ? (ubx-lbx) * lby : 0.0);
296  }
297  else if( SCIPisRelEQ(scip, lby, uby) )
298  {
299  /* x*y = lby * x + (y-lby) * x >= lby * x + (y-lby) * lbx >= lby * x + min{(uby-lby) * lbx, 0 * lbx} */
300  coefx = bilincoef * lby;
301  coefy = 0.0;
302  constant = bilincoef * (lbx < 0.0 ? (uby-lby) * lbx : 0.0);
303  }
304  else
305  {
306  coefx = bilincoef * lby;
307  coefy = bilincoef * lbx;
308  constant = -bilincoef * lbx * lby;
309  }
310  }
311  else if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, uby) )
312  {
313  if( SCIPisRelEQ(scip, lbx, ubx) )
314  {
315  /* x*y = ubx * y + (x-ubx) * y >= ubx * y + (x-ubx) * uby >= ubx * y + min{(lbx-ubx) * uby, 0 * uby} */
316  coefx = 0.0;
317  coefy = bilincoef * ubx;
318  constant = bilincoef * (uby > 0.0 ? (lbx - ubx) * uby : 0.0);
319  }
320  else if( SCIPisRelEQ(scip, lby, uby) )
321  {
322  /* x*y = uby * x + (y-uby) * x >= uby * x + (y-uby) * ubx >= uby * x + min{(lby-uby) * ubx, 0 * ubx} */
323  coefx = bilincoef * uby;
324  coefy = 0.0;
325  constant = bilincoef * (ubx > 0.0 ? (lby - uby) * ubx : 0.0);
326  }
327  else
328  {
329  coefx = bilincoef * uby;
330  coefy = bilincoef * ubx;
331  constant = -bilincoef * ubx * uby;
332  }
333  }
334  else
335  {
336  *success = FALSE;
337  return;
338  }
339  }
340  else
341  {
342  /* either x or y is not fixed and coef < 0.0 */
343  if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, -lby) &&
344  (SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip, uby)
345  || (ubx - lbx) * (refpointy - lby) <= (uby - lby) * (refpointx - lbx)) )
346  {
347  if( SCIPisRelEQ(scip, lbx, ubx) )
348  {
349  /* x*y = ubx * y + (x-ubx) * y <= ubx * y + (x-ubx) * lby <= ubx * y + max{(lbx-ubx) * lby, 0 * lby} */
350  coefx = 0.0;
351  coefy = bilincoef * ubx;
352  constant = bilincoef * (lby < 0.0 ? (lbx - ubx) * lby : 0.0);
353  }
354  else if( SCIPisRelEQ(scip, lby, uby) )
355  {
356  /* x*y = lby * x + (y-lby) * x <= lby * x + (y-lby) * ubx <= lby * x + max{(uby-lby) * ubx, 0 * ubx} */
357  coefx = bilincoef * lby;
358  coefy = 0.0;
359  constant = bilincoef * (ubx > 0.0 ? (uby - lby) * ubx : 0.0);
360  }
361  else
362  {
363  coefx = bilincoef * lby;
364  coefy = bilincoef * ubx;
365  constant = -bilincoef * ubx * lby;
366  }
367  }
368  else if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, uby) )
369  {
370  if( SCIPisRelEQ(scip, lbx, ubx) )
371  {
372  /* x*y = lbx * y + (x-lbx) * y <= lbx * y + (x-lbx) * uby <= lbx * y + max{(ubx-lbx) * uby, 0 * uby} */
373  coefx = 0.0;
374  coefy = bilincoef * lbx;
375  constant = bilincoef * (uby > 0.0 ? (ubx - lbx) * uby : 0.0);
376  }
377  else if( SCIPisRelEQ(scip, lby, uby) )
378  {
379  /* x*y = uby * x + (y-uby) * x <= uby * x + (y-uby) * lbx <= uby * x + max{(lby-uby) * lbx, 0 * lbx} */
380  coefx = bilincoef * uby;
381  coefy = 0.0;
382  constant = bilincoef * (lbx < 0.0 ? (lby - uby) * lbx : 0.0);
383  }
384  else
385  {
386  coefx = bilincoef * uby;
387  coefy = bilincoef * lbx;
388  constant = -bilincoef * lbx * uby;
389  }
390  }
391  else
392  {
393  *success = FALSE;
394  return;
395  }
396  }
397 
398  if( SCIPisInfinity(scip, REALABS(coefx)) || SCIPisInfinity(scip, REALABS(coefy))
399  || SCIPisInfinity(scip, REALABS(constant)) )
400  {
401  *success = FALSE;
402  return;
403  }
404 
405  if( overestimate )
406  {
407  coefx = -coefx;
408  coefy = -coefy;
409  constant = -constant;
410  }
411 
412  SCIPdebugMsg(scip, "%.15g * x[%.15g,%.15g] * y[%.15g,%.15g] %c= %.15g * x %+.15g * y %+.15g\n", bilincoef, lbx, ubx,
413  lby, uby, overestimate ? '<' : '>', coefx, coefy, constant);
414 
415  *lincoefx += coefx;
416  *lincoefy += coefy;
417  *linconstant += constant;
418 }
419 
420 
421 /** computes coefficients of linearization of a bilinear term in a reference point when given a linear inequality
422  * involving only the variables of the bilinear term
423  *
424  * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
425  * by Marco Locatelli
426  */
428  SCIP* scip, /**< SCIP data structure */
429  SCIP_Real bilincoef, /**< coefficient of bilinear term */
430  SCIP_Real lbx, /**< lower bound on first variable */
431  SCIP_Real ubx, /**< upper bound on first variable */
432  SCIP_Real refpointx, /**< reference point for first variable */
433  SCIP_Real lby, /**< lower bound on second variable */
434  SCIP_Real uby, /**< upper bound on second variable */
435  SCIP_Real refpointy, /**< reference point for second variable */
436  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
437  SCIP_Real xcoef, /**< x coefficient of linear inequality; must be in {-1,0,1} */
438  SCIP_Real ycoef, /**< y coefficient of linear inequality */
439  SCIP_Real constant, /**< constant of linear inequality */
440  SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
441  SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
442  SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
443  SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
444  )
445 {
446  SCIP_Real xs[2] = {lbx, ubx};
447  SCIP_Real ys[2] = {lby, uby};
448  SCIP_Real minx;
449  SCIP_Real maxx;
450  SCIP_Real miny;
451  SCIP_Real maxy;
452  SCIP_Real QUAD(lincoefyq);
453  SCIP_Real QUAD(lincoefxq);
454  SCIP_Real QUAD(linconstantq);
455  SCIP_Real QUAD(denomq);
456  SCIP_Real QUAD(mjq);
457  SCIP_Real QUAD(qjq);
458  SCIP_Real QUAD(xjq);
459  SCIP_Real QUAD(yjq);
460  SCIP_Real QUAD(tmpq);
461  SCIP_Real vx;
462  SCIP_Real vy;
463  int n;
464  int i;
465 
466  assert(scip != NULL);
467  assert(!SCIPisInfinity(scip, lbx));
468  assert(!SCIPisInfinity(scip, -ubx));
469  assert(!SCIPisInfinity(scip, lby));
470  assert(!SCIPisInfinity(scip, -uby));
471  assert(SCIPisLE(scip, lbx, ubx));
472  assert(SCIPisLE(scip, lby, uby));
473  assert(SCIPisLE(scip, lbx, refpointx));
474  assert(SCIPisGE(scip, ubx, refpointx));
475  assert(SCIPisLE(scip, lby, refpointy));
476  assert(SCIPisGE(scip, uby, refpointy));
477  assert(lincoefx != NULL);
478  assert(lincoefy != NULL);
479  assert(linconstant != NULL);
480  assert(success != NULL);
481  assert(xcoef == 0.0 || xcoef == -1.0 || xcoef == 1.0); /*lint !e777*/
482  assert(ycoef != SCIP_INVALID && ycoef != 0.0); /*lint !e777*/
483  assert(constant != SCIP_INVALID); /*lint !e777*/
484 
485  *success = FALSE;
486  *lincoefx = SCIP_INVALID;
487  *lincoefy = SCIP_INVALID;
488  *linconstant = SCIP_INVALID;
489 
490  /* reference point does not satisfy linear inequality */
491  if( SCIPisFeasGT(scip, xcoef * refpointx - ycoef * refpointy - constant, 0.0) )
492  return;
493 
494  /* compute minimal and maximal bounds on x and y for accepting the reference point */
495  minx = lbx + 0.01 * (ubx-lbx);
496  maxx = ubx - 0.01 * (ubx-lbx);
497  miny = lby + 0.01 * (uby-lby);
498  maxy = uby - 0.01 * (uby-lby);
499 
500  /* check whether the reference point is in [minx,maxx]x[miny,maxy] */
501  if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
502  || SCIPisLE(scip, refpointy, miny) || SCIPisGE(scip, refpointy, maxy) )
503  return;
504 
505  /* always consider xy without the bilinear coefficient */
506  if( bilincoef < 0.0 )
507  overestimate = !overestimate;
508 
509  /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
510  /* mj = xcoef / ycoef */
511  SCIPquadprecDivDD(mjq, xcoef, ycoef);
512 
513  /* qj = -constant / ycoef */
514  SCIPquadprecDivDD(qjq, -constant, ycoef);
515 
516  /* mj > 0 => underestimate; mj < 0 => overestimate */
517  if( SCIPisNegative(scip, QUAD_TO_DBL(mjq)) != overestimate )
518  return;
519 
520  /* get the corner point that satisfies the linear inequality xcoef*x <= ycoef*y + constant */
521  if( !overestimate )
522  {
523  ys[0] = uby;
524  ys[1] = lby;
525  }
526 
527  vx = SCIP_INVALID;
528  vy = SCIP_INVALID;
529  n = 0;
530  for( i = 0; i < 2; ++i )
531  {
532  SCIP_Real activity = xcoef * xs[i] - ycoef * ys[i] - constant;
533  if( SCIPisLE(scip, activity, 0.0) )
534  {
535  /* corner point is satisfies inequality */
536  vx = xs[i];
537  vy = ys[i];
538  }
539  else if( SCIPisFeasGT(scip, activity, 0.0) )
540  /* corner point is clearly cut off */
541  ++n;
542  }
543 
544  /* skip if no corner point satisfies the inequality or if no corner point is cut off (that is, all corner points satisfy the inequality almost [1e-9..1e-6]) */
545  if( n != 1 || vx == SCIP_INVALID || vy == SCIP_INVALID ) /*lint !e777*/
546  return;
547 
548  /* denom = mj*(refpointx - vx) + vy - refpointy */
549  SCIPquadprecSumDD(denomq, refpointx, -vx); /* refpoint - vx */
550  SCIPquadprecProdQQ(denomq, denomq, mjq); /* mj * (refpoint - vx) */
551  SCIPquadprecSumQD(denomq, denomq, vy); /* mj * (refpoint - vx) + vy */
552  SCIPquadprecSumQD(denomq, denomq, -refpointy); /* mj * (refpoint - vx) + vy - refpointy */
553 
554  if( SCIPisZero(scip, QUAD_TO_DBL(denomq)) )
555  return;
556 
557  /* (xj,yj) is the projection onto the line xcoef*x = ycoef*y + constant */
558  /* xj = (refpointx*(vy - qj) - vx*(refpointy - qj)) / denom */
559  SCIPquadprecProdQD(xjq, qjq, -1.0); /* - qj */
560  SCIPquadprecSumQD(xjq, xjq, vy); /* vy - qj */
561  SCIPquadprecProdQD(xjq, xjq, refpointx); /* refpointx * (vy - qj) */
562  SCIPquadprecProdQD(tmpq, qjq, -1.0); /* - qj */
563  SCIPquadprecSumQD(tmpq, tmpq, refpointy); /* refpointy - qj */
564  SCIPquadprecProdQD(tmpq, tmpq, -vx); /* - vx * (refpointy - qj) */
565  SCIPquadprecSumQQ(xjq, xjq, tmpq); /* refpointx * (vy - qj) - vx * (refpointy - qj) */
566  SCIPquadprecDivQQ(xjq, xjq, denomq); /* (refpointx * (vy - qj) - vx * (refpointy - qj)) / denom */
567 
568  /* yj = mj * xj + qj */
569  SCIPquadprecProdQQ(yjq, mjq, xjq);
570  SCIPquadprecSumQQ(yjq, yjq, qjq);
571 
572  assert(SCIPisFeasEQ(scip, xcoef*QUAD_TO_DBL(xjq) - ycoef*QUAD_TO_DBL(yjq) - constant, 0.0));
573 
574  /* check whether the projection is in [minx,maxx] x [miny,maxy]; this avoids numerical difficulties when the
575  * projection is close to the variable bounds
576  */
577  if( SCIPisLE(scip, QUAD_TO_DBL(xjq), minx) || SCIPisGE(scip, QUAD_TO_DBL(xjq), maxx)
578  || SCIPisLE(scip, QUAD_TO_DBL(yjq), miny) || SCIPisGE(scip, QUAD_TO_DBL(yjq), maxy) )
579  return;
580 
581  assert(vy - QUAD_TO_DBL(mjq)*vx - QUAD_TO_DBL(qjq) != 0.0);
582 
583  /* lincoefy = (mj*SQR(xj) - 2.0*mj*vx*xj - qj*vx + vx*vy) / (vy - mj*vx - qj) */
584  SCIPquadprecSquareQ(lincoefyq, xjq); /* xj^2 */
585  SCIPquadprecProdQQ(lincoefyq, lincoefyq, mjq); /* mj * xj^2 */
586  SCIPquadprecProdQQ(tmpq, mjq, xjq); /* mj * xj */
587  SCIPquadprecProdQD(tmpq, tmpq, -2.0 * vx); /* -2 * vx * mj * xj */
588  SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj */
589  SCIPquadprecProdQD(tmpq, qjq, -vx); /* -qj * vx */
590  SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx */
591  SCIPquadprecProdDD(tmpq, vx, vy); /* vx * vy */
592  SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy */
593  SCIPquadprecProdQD(tmpq, mjq, vx); /* mj * vx */
594  SCIPquadprecSumQD(tmpq, tmpq, -vy); /* -vy + mj * vx */
595  SCIPquadprecSumQQ(tmpq, tmpq, qjq); /* -vy + mj * vx + qj */
596  QUAD_SCALE(tmpq, -1.0); /* vy - mj * vx - qj */
597  SCIPquadprecDivQQ(lincoefyq, lincoefyq, tmpq); /* (mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy) / (vy - mj * vx - qj) */
598 
599  /* lincoefx = 2.0*mj*xj + qj - mj*(*lincoefy) */
600  SCIPquadprecProdQQ(lincoefxq, mjq, xjq); /* mj * xj */
601  QUAD_SCALE(lincoefxq, 2.0); /* 2 * mj * xj */
602  SCIPquadprecSumQQ(lincoefxq, lincoefxq, qjq); /* 2 * mj * xj + qj */
603  SCIPquadprecProdQQ(tmpq, mjq, lincoefyq); /* mj * lincoefy */
604  QUAD_SCALE(tmpq, -1.0); /* - mj * lincoefy */
605  SCIPquadprecSumQQ(lincoefxq, lincoefxq, tmpq); /* 2 * mj * xj + qj - mj * lincoefy */
606 
607  /* linconstant = -mj*SQR(xj) - (*lincoefy)*qj */
608  SCIPquadprecSquareQ(linconstantq, xjq); /* xj^2 */
609  SCIPquadprecProdQQ(linconstantq, linconstantq, mjq); /* mj * xj^2 */
610  QUAD_SCALE(linconstantq, -1.0); /* - mj * xj^2 */
611  SCIPquadprecProdQQ(tmpq, lincoefyq, qjq); /* lincoefy * qj */
612  QUAD_SCALE(tmpq, -1.0); /* - lincoefy * qj */
613  SCIPquadprecSumQQ(linconstantq, linconstantq, tmpq); /* - mj * xj^2 - lincoefy * qj */
614 
615  /* consider the bilinear coefficient */
616  SCIPquadprecProdQD(lincoefxq, lincoefxq, bilincoef);
617  SCIPquadprecProdQD(lincoefyq, lincoefyq, bilincoef);
618  SCIPquadprecProdQD(linconstantq, linconstantq, bilincoef);
619  *lincoefx = QUAD_TO_DBL(lincoefxq);
620  *lincoefy = QUAD_TO_DBL(lincoefyq);
621  *linconstant = QUAD_TO_DBL(linconstantq);
622 
623  /* cut needs to be tight at (vx,vy) and (xj,yj); otherwise we consider the cut to be numerically bad */
624  *success = SCIPisFeasEQ(scip, (*lincoefx)*vx + (*lincoefy)*vy + (*linconstant), bilincoef*vx*vy)
625  && SCIPisFeasEQ(scip, (*lincoefx)*QUAD_TO_DBL(xjq) + (*lincoefy)*QUAD_TO_DBL(yjq) + (*linconstant), bilincoef*QUAD_TO_DBL(xjq)*QUAD_TO_DBL(yjq));
626 
627 #ifndef NDEBUG
628  {
629  SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
630 
631  /* cut needs to under- or overestimate the bilinear term at the reference point */
632  if( bilincoef < 0.0 )
633  overestimate = !overestimate;
634 
635  if( overestimate )
636  assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
637  else
638  assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
639  }
640 #endif
641 }
642 
643 /** helper function to compute the convex envelope of a bilinear term when two linear inequalities are given; we
644  * use the same notation and formulas as in Locatelli 2016
645  */
646 static
648  SCIP* scip, /**< SCIP data structure */
649  SCIP_Real x, /**< reference point for x */
650  SCIP_Real y, /**< reference point for y */
651  SCIP_Real mi, /**< coefficient of x in the first linear inequality */
652  SCIP_Real qi, /**< constant in the first linear inequality */
653  SCIP_Real mj, /**< coefficient of x in the second linear inequality */
654  SCIP_Real qj, /**< constant in the second linear inequality */
655  SCIP_Real* RESTRICT xi, /**< buffer to store x coordinate of the first point */
656  SCIP_Real* RESTRICT yi, /**< buffer to store y coordinate of the first point */
657  SCIP_Real* RESTRICT xj, /**< buffer to store x coordinate of the second point */
658  SCIP_Real* RESTRICT yj, /**< buffer to store y coordinate of the second point */
659  SCIP_Real* RESTRICT xcoef, /**< buffer to store the x coefficient of the envelope */
660  SCIP_Real* RESTRICT ycoef, /**< buffer to store the y coefficient of the envelope */
661  SCIP_Real* RESTRICT constant /**< buffer to store the constant of the envelope */
662  )
663 {
664  SCIP_Real QUAD(xiq);
665  SCIP_Real QUAD(yiq);
666  SCIP_Real QUAD(xjq);
667  SCIP_Real QUAD(yjq);
668  SCIP_Real QUAD(xcoefq);
669  SCIP_Real QUAD(ycoefq);
670  SCIP_Real QUAD(constantq);
671  SCIP_Real QUAD(tmpq);
672 
673  assert(xi != NULL);
674  assert(yi != NULL);
675  assert(xj != NULL);
676  assert(yj != NULL);
677  assert(xcoef != NULL);
678  assert(ycoef != NULL);
679  assert(constant != NULL);
680 
681  if( SCIPisEQ(scip, mi, mj) )
682  {
683  /* xi = (x + mi * y - qi) / (2.0*mi) */
684  SCIPquadprecProdDD(xiq, mi, y);
685  SCIPquadprecSumQD(xiq, xiq, x);
686  SCIPquadprecSumQD(xiq, xiq, -qi);
687  SCIPquadprecDivQD(xiq, xiq, 2.0 * mi);
688  assert(EPSEQ((x + mi * y - qi) / (2.0*mi), QUAD_TO_DBL(xiq), 1e-3));
689 
690  /* yi = mi*(*xi) + qi */
691  SCIPquadprecProdQD(yiq, xiq, mi);
692  SCIPquadprecSumQD(yiq, yiq, qi);
693  assert(EPSEQ(mi*QUAD_TO_DBL(xiq) + qi, QUAD_TO_DBL(yiq), 1e-3));
694 
695  /* xj = (*xi) + (qi - qj)/ (2.0*mi) */
696  SCIPquadprecSumDD(xjq, qi, -qj);
697  SCIPquadprecDivQD(xjq, xjq, 2.0 * mi);
698  SCIPquadprecSumQQ(xjq, xjq, xiq);
699  assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj)/ (2.0*mi), QUAD_TO_DBL(xjq), 1e-3));
700 
701  /* yj = mj * (*xj) + qj */
702  SCIPquadprecProdQD(yjq, xjq, mj);
703  SCIPquadprecSumQD(yjq, yjq, qj);
704  assert(EPSEQ(mj * QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
705 
706  /* ycoef = (*xi) + (qi - qj) / (4.0*mi) note that this is wrong in Locatelli 2016 */
707  SCIPquadprecSumDD(ycoefq, qi, -qj);
708  SCIPquadprecDivQD(ycoefq, ycoefq, 4.0 * mi);
709  SCIPquadprecSumQQ(ycoefq, ycoefq, xiq);
710  assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj) / (4.0*mi), QUAD_TO_DBL(ycoefq), 1e-3));
711 
712  /* xcoef = 2.0*mi*(*xi) - mi * (*ycoef) + qi */
713  SCIPquadprecProdQD(xcoefq, xiq, 2.0 * mi);
714  SCIPquadprecProdQD(tmpq, ycoefq, -mi);
715  SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
716  SCIPquadprecSumQD(xcoefq, xcoefq, qi);
717  assert(EPSEQ(2.0*mi*QUAD_TO_DBL(xiq) - mi * QUAD_TO_DBL(ycoefq) + qi, QUAD_TO_DBL(xcoefq), 1e-3));
718 
719  /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
720  SCIPquadprecSquareQ(constantq, xjq);
721  SCIPquadprecProdQD(constantq, constantq, -mj);
722  SCIPquadprecProdQD(tmpq, ycoefq, -qj);
723  SCIPquadprecSumQQ(constantq, constantq, tmpq);
724  /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
725 
726  *xi = QUAD_TO_DBL(xiq);
727  *yi = QUAD_TO_DBL(yiq);
728  *xj = QUAD_TO_DBL(xjq);
729  *yj = QUAD_TO_DBL(yjq);
730  *ycoef = QUAD_TO_DBL(ycoefq);
731  *xcoef = QUAD_TO_DBL(xcoefq);
732  *constant = QUAD_TO_DBL(constantq);
733  }
734  else if( mi > 0.0 )
735  {
736  assert(mj > 0.0);
737 
738  /* xi = (y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)) */
739  SCIPquadprecProdDD(xiq, mi, mj);
740  SCIPquadprecSqrtQ(xiq, xiq);
741  SCIPquadprecProdQD(xiq, xiq, x);
742  SCIPquadprecSumQD(xiq, xiq, y);
743  SCIPquadprecSumQD(xiq, xiq, -qi); /* (y + SQRT(mi*mj)*x - qi) */
744  SCIPquadprecProdDD(tmpq, mi, mj);
745  SCIPquadprecSqrtQ(tmpq, tmpq);
746  SCIPquadprecSumQD(tmpq, tmpq, REALABS(mi)); /* REALABS(mi) + SQRT(mi*mj) */
747  SCIPquadprecDivQQ(xiq, xiq, tmpq);
748  assert(EPSEQ((y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)), QUAD_TO_DBL(xiq), 1e-3));
749 
750  /* yi = mi*(*xi) + qi */
751  SCIPquadprecProdQD(yiq, xiq, mi);
752  SCIPquadprecSumQD(yiq, yiq, qi);
753  assert(EPSEQ(mi*(QUAD_TO_DBL(xiq)) + qi, QUAD_TO_DBL(yiq), 1e-3));
754 
755  /* xj = (y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)) */
756  SCIPquadprecProdDD(xjq, mi, mj);
757  SCIPquadprecSqrtQ(xjq, xjq);
758  SCIPquadprecProdQD(xjq, xjq, x);
759  SCIPquadprecSumQD(xjq, xjq, y);
760  SCIPquadprecSumQD(xjq, xjq, -qj); /* (y + SQRT(mi*mj)*x - qj) */
761  SCIPquadprecProdDD(tmpq, mi, mj);
762  SCIPquadprecSqrtQ(tmpq, tmpq);
763  SCIPquadprecSumQD(tmpq, tmpq, REALABS(mj)); /* REALABS(mj) + SQRT(mi*mj) */
764  SCIPquadprecDivQQ(xjq, xjq, tmpq);
765  assert(EPSEQ((y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)), QUAD_TO_DBL(xjq), 1e-3));
766 
767  /* yj = mj*(*xj) + qj */
768  SCIPquadprecProdQD(yjq, xjq, mj);
769  SCIPquadprecSumQD(yjq, yjq, qj);
770  assert(EPSEQ(mj*QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
771 
772  /* ycoef = (2.0*mj*(*xj) + qj - 2.0*mi*(*xi) - qi) / (mj - mi) */
773  SCIPquadprecProdQD(ycoefq, xjq, 2.0 * mj);
774  SCIPquadprecSumQD(ycoefq, ycoefq, qj);
775  SCIPquadprecProdQD(tmpq, xiq, -2.0 * mi);
776  SCIPquadprecSumQQ(ycoefq, ycoefq, tmpq);
777  SCIPquadprecSumQD(ycoefq, ycoefq, -qi);
778  SCIPquadprecSumDD(tmpq, mj, -mi);
779  SCIPquadprecDivQQ(ycoefq, ycoefq, tmpq);
780  assert(EPSEQ((2.0*mj*QUAD_TO_DBL(xjq) + qj - 2.0*mi*QUAD_TO_DBL(xiq) - qi) / (mj - mi), QUAD_TO_DBL(ycoefq), 1e-3));
781 
782  /* xcoef = 2.0*mj*(*xj) + qj - mj*(*ycoef) */
783  SCIPquadprecProdQD(xcoefq, xjq, 2.0 * mj);
784  SCIPquadprecSumQD(xcoefq, xcoefq, qj);
785  SCIPquadprecProdQD(tmpq, ycoefq, -mj);
786  SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
787  assert(EPSEQ(2.0*mj*QUAD_TO_DBL(xjq) + qj - mj*QUAD_TO_DBL(ycoefq), QUAD_TO_DBL(xcoefq), 1e-3));
788 
789  /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
790  SCIPquadprecSquareQ(constantq, xjq);
791  SCIPquadprecProdQD(constantq, constantq, -mj);
792  SCIPquadprecProdQD(tmpq, ycoefq, -qj);
793  SCIPquadprecSumQQ(constantq, constantq, tmpq);
794  /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
795 
796  *xi = QUAD_TO_DBL(xiq);
797  *yi = QUAD_TO_DBL(yiq);
798  *xj = QUAD_TO_DBL(xjq);
799  *yj = QUAD_TO_DBL(yjq);
800  *ycoef = QUAD_TO_DBL(ycoefq);
801  *xcoef = QUAD_TO_DBL(xcoefq);
802  *constant = QUAD_TO_DBL(constantq);
803  }
804  else
805  {
806  assert(mi < 0.0 && mj < 0.0);
807 
808  /* apply variable transformation x = -x in case for overestimation */
809  computeBilinEnvelope2(scip, -x, y, -mi, qi, -mj, qj, xi, yi, xj, yj, xcoef, ycoef, constant);
810 
811  /* revert transformation; multiply cut by -1 and change -x by x */
812  *xi = -(*xi);
813  *xj = -(*xj);
814  *ycoef = -(*ycoef);
815  *constant = -(*constant);
816  }
817 }
818 
819 /** computes coefficients of linearization of a bilinear term in a reference point when given two linear inequality
820  * involving only the variables of the bilinear term
821  *
822  * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
823  * by Marco Locatelli
824  *
825  */
827  SCIP* scip, /**< SCIP data structure */
828  SCIP_Real bilincoef, /**< coefficient of bilinear term */
829  SCIP_Real lbx, /**< lower bound on first variable */
830  SCIP_Real ubx, /**< upper bound on first variable */
831  SCIP_Real refpointx, /**< reference point for first variable */
832  SCIP_Real lby, /**< lower bound on second variable */
833  SCIP_Real uby, /**< upper bound on second variable */
834  SCIP_Real refpointy, /**< reference point for second variable */
835  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
836  SCIP_Real xcoef1, /**< x coefficient of linear inequality; must be in {-1,0,1} */
837  SCIP_Real ycoef1, /**< y coefficient of linear inequality */
838  SCIP_Real constant1, /**< constant of linear inequality */
839  SCIP_Real xcoef2, /**< x coefficient of linear inequality; must be in {-1,0,1} */
840  SCIP_Real ycoef2, /**< y coefficient of linear inequality */
841  SCIP_Real constant2, /**< constant of linear inequality */
842  SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
843  SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
844  SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
845  SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
846  )
847 {
848  SCIP_Real mi, mj, qi, qj, xi, xj, yi, yj;
849  SCIP_Real xcoef, ycoef, constant;
850  SCIP_Real minx, maxx, miny, maxy;
851 
852  assert(scip != NULL);
853  assert(!SCIPisInfinity(scip, lbx));
854  assert(!SCIPisInfinity(scip, -ubx));
855  assert(!SCIPisInfinity(scip, lby));
856  assert(!SCIPisInfinity(scip, -uby));
857  assert(SCIPisLE(scip, lbx, ubx));
858  assert(SCIPisLE(scip, lby, uby));
859  assert(SCIPisLE(scip, lbx, refpointx));
860  assert(SCIPisGE(scip, ubx, refpointx));
861  assert(SCIPisLE(scip, lby, refpointy));
862  assert(SCIPisGE(scip, uby, refpointy));
863  assert(lincoefx != NULL);
864  assert(lincoefy != NULL);
865  assert(linconstant != NULL);
866  assert(success != NULL);
867  assert(xcoef1 != 0.0 && xcoef1 != SCIP_INVALID); /*lint !e777*/
868  assert(ycoef1 != SCIP_INVALID && ycoef1 != 0.0); /*lint !e777*/
869  assert(constant1 != SCIP_INVALID); /*lint !e777*/
870  assert(xcoef2 != 0.0 && xcoef2 != SCIP_INVALID); /*lint !e777*/
871  assert(ycoef2 != SCIP_INVALID && ycoef2 != 0.0); /*lint !e777*/
872  assert(constant2 != SCIP_INVALID); /*lint !e777*/
873 
874  *success = FALSE;
875  *lincoefx = SCIP_INVALID;
876  *lincoefy = SCIP_INVALID;
877  *linconstant = SCIP_INVALID;
878 
879  /* reference point does not satisfy linear inequalities */
880  if( SCIPisFeasGT(scip, xcoef1 * refpointx - ycoef1 * refpointy - constant1, 0.0)
881  || SCIPisFeasGT(scip, xcoef2 * refpointx - ycoef2 * refpointy - constant2, 0.0) )
882  return;
883 
884  /* compute minimal and maximal bounds on x and y for accepting the reference point */
885  minx = lbx + 0.01 * (ubx-lbx);
886  maxx = ubx - 0.01 * (ubx-lbx);
887  miny = lby + 0.01 * (uby-lby);
888  maxy = uby - 0.01 * (uby-lby);
889 
890  /* check the reference point is in the interior of the domain */
891  if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
892  || SCIPisLE(scip, refpointy, miny) || SCIPisFeasGE(scip, refpointy, maxy) )
893  return;
894 
895  /* the sign of the x-coefficients of the two inequalities must be different; otherwise the convex or concave
896  * envelope can be computed via SCIPcomputeBilinEnvelope1 for each inequality separately
897  */
898  if( (xcoef1 > 0) == (xcoef2 > 0) )
899  return;
900 
901  /* always consider xy without the bilinear coefficient */
902  if( bilincoef < 0.0 )
903  overestimate = !overestimate;
904 
905  /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
906  mi = xcoef1 / ycoef1;
907  qi = -constant1 / ycoef1;
908  mj = xcoef2 / ycoef2;
909  qj = -constant2 / ycoef2;
910 
911  /* mi, mj > 0 => underestimate; mi, mj < 0 => overestimate */
912  if( SCIPisNegative(scip, mi) != overestimate || SCIPisNegative(scip, mj) != overestimate )
913  return;
914 
915  /* compute cut according to Locatelli 2016 */
916  computeBilinEnvelope2(scip, refpointx, refpointy, mi, qi, mj, qj, &xi, &yi, &xj, &yj, &xcoef, &ycoef, &constant);
917  assert(SCIPisRelEQ(scip, mi*xi + qi, yi));
918  assert(SCIPisRelEQ(scip, mj*xj + qj, yj));
919 
920  /* it might happen that (xi,yi) = (xj,yj) if the two lines intersect */
921  if( SCIPisEQ(scip, xi, xj) && SCIPisEQ(scip, yi, yj) )
922  return;
923 
924  /* check whether projected points are in the interior */
925  if( SCIPisLE(scip, xi, minx) || SCIPisGE(scip, xi, maxx) || SCIPisLE(scip, yi, miny) || SCIPisGE(scip, yi, maxy) )
926  return;
927  if( SCIPisLE(scip, xj, minx) || SCIPisGE(scip, xj, maxx) || SCIPisLE(scip, yj, miny) || SCIPisGE(scip, yj, maxy) )
928  return;
929 
930  *lincoefx = bilincoef * xcoef;
931  *lincoefy = bilincoef * ycoef;
932  *linconstant = bilincoef * constant;
933 
934  /* cut needs to be tight at (vx,vy) and (xj,yj) */
935  *success = SCIPisFeasEQ(scip, (*lincoefx)*xi + (*lincoefy)*yi + (*linconstant), bilincoef*xi*yi)
936  && SCIPisFeasEQ(scip, (*lincoefx)*xj + (*lincoefy)*yj + (*linconstant), bilincoef*xj*yj);
937 
938 #ifndef NDEBUG
939  {
940  SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
941 
942  /* cut needs to under- or overestimate the bilinear term at the reference point */
943  if( bilincoef < 0.0 )
944  overestimate = !overestimate;
945 
946  if( overestimate )
947  assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
948  else
949  assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
950  }
951 #endif
952 }
953 
954 /** creates an NLP relaxation and stores it in a given NLPI problem; the function computes for each variable which the
955  * number of non-linearly occurrences and stores it in the nlscore array
956  *
957  * @note the first row corresponds always to the cutoff row (even if cutoffbound is SCIPinfinity(scip))
958  **/
960  SCIP* scip, /**< SCIP data structure */
961  SCIP_NLPI* nlpi, /**< interface to NLP solver */
962  SCIP_NLROW** nlrows, /**< nonlinear rows */
963  int nnlrows, /**< total number of nonlinear rows */
964  SCIP_NLPIPROBLEM* nlpiprob, /**< empty nlpi problem */
965  SCIP_HASHMAP* var2idx, /**< empty hash map to store mapping between variables and indices in nlpi
966  * problem */
967  SCIP_HASHMAP* nlrow2idx, /**< empty hash map to store mapping between variables and indices in nlpi
968  * problem, can be NULL */
969  SCIP_Real* nlscore, /**< array to store the score of each nonlinear variable (NULL if not
970  * needed) */
971  SCIP_Real cutoffbound, /**< cutoff bound */
972  SCIP_Bool setobj, /**< should the objective function be set? */
973  SCIP_Bool onlyconvex /**< filter only for convex constraints */
974  )
975 {
976  SCIP_EXPRTREE** exprtrees;
977  int** exprvaridxs;
978  SCIP_QUADELEM** quadelems;
979  int* nquadelems;
980  SCIP_Real** linvals;
981  int** lininds;
982  int* nlininds;
983  SCIP_Real* lhss;
984  SCIP_Real* rhss;
985  const char** names;
986  SCIP_VAR** vars;
987  int nvars;
988  SCIP_Real* lbs;
989  SCIP_Real* ubs;
990  SCIP_Real* objvals = NULL;
991  int* objinds = NULL;
992  const char** varnames;
993  int nobjinds;
994  int nconss;
995  int i;
996 
997  assert(nlpiprob != NULL);
998  assert(var2idx != NULL);
999  assert(nlrows != NULL);
1000  assert(nnlrows > 0);
1001  assert(nlpi != NULL);
1002 
1003  SCIPdebugMsg(scip, "call SCIPcreateConvexNlpNlobbt() with cutoffbound %g\n", cutoffbound);
1004 
1005  if( nlscore != NULL )
1006  {
1007  BMSclearMemoryArray(nlscore, SCIPgetNVars(scip));
1008  }
1009  vars = SCIPgetVars(scip);
1010  nvars = SCIPgetNVars(scip);
1011  nconss = 0;
1012 
1013  SCIP_CALL( SCIPallocBufferArray(scip, &exprtrees, nnlrows + 1) );
1014  SCIP_CALL( SCIPallocBufferArray(scip, &exprvaridxs, nnlrows + 1) );
1015  SCIP_CALL( SCIPallocBufferArray(scip, &quadelems, nnlrows + 1) );
1016  SCIP_CALL( SCIPallocBufferArray(scip, &nquadelems, nnlrows + 1) );
1017  SCIP_CALL( SCIPallocBufferArray(scip, &linvals, nnlrows + 1) );
1018  SCIP_CALL( SCIPallocBufferArray(scip, &lininds, nnlrows + 1) );
1019  SCIP_CALL( SCIPallocBufferArray(scip, &nlininds, nnlrows + 1) );
1020  SCIP_CALL( SCIPallocBufferArray(scip, &names, nnlrows + 1) );
1021  SCIP_CALL( SCIPallocBufferArray(scip, &lhss, nnlrows + 1) );
1022  SCIP_CALL( SCIPallocBufferArray(scip, &rhss, nnlrows + 1) );
1023 
1024  if( setobj )
1025  {
1026  SCIP_CALL( SCIPallocBufferArray(scip, &objvals, nvars) );
1027  SCIP_CALL( SCIPallocBufferArray(scip, &objinds, nvars) );
1028  }
1029 
1030  SCIP_CALL( SCIPallocBufferArray(scip, &lbs, nvars) );
1031  SCIP_CALL( SCIPallocBufferArray(scip, &ubs, nvars) );
1032  SCIP_CALL( SCIPallocBufferArray(scip, &varnames, nvars) );
1033 
1034  /* create a unique mapping between variables and {0,..,nvars-1} */
1035  nobjinds = 0;
1036  for( i = 0; i < nvars; ++i )
1037  {
1038  assert(vars[i] != NULL);
1039  SCIP_CALL( SCIPhashmapInsertInt(var2idx, (void*)vars[i], i) );
1040 
1041  lbs[i] = SCIPvarGetLbLocal(vars[i]);
1042  ubs[i] = SCIPvarGetUbLocal(vars[i]);
1043  varnames[i] = SCIPvarGetName(vars[i]);
1044 
1045  /* collect non-zero objective coefficients */
1046  if( setobj && !SCIPisZero(scip, SCIPvarGetObj(vars[i])) )
1047  {
1048  assert(objvals != NULL);
1049  assert(objinds != NULL);
1050 
1051  objvals[nobjinds] = SCIPvarGetObj(vars[i]);
1052  objinds[nobjinds] = i;
1053  ++nobjinds;
1054  }
1055  }
1056 
1057  /* add variables */
1058  SCIP_CALL( SCIPnlpiAddVars(nlpi, nlpiprob, nvars, lbs, ubs, varnames) );
1059  SCIPfreeBufferArray(scip, &varnames);
1060  SCIPfreeBufferArray(scip, &ubs);
1061  SCIPfreeBufferArray(scip, &lbs);
1062 
1063  /* set the objective function */
1064  if( setobj )
1065  {
1066  if( nobjinds > 0 )
1067  {
1068  SCIP_CALL( SCIPnlpiSetObjective(nlpi, nlpiprob, nobjinds, objinds, objvals, 0, NULL, NULL, NULL, 0.0) );
1069  }
1070 
1071  SCIPfreeBufferArray(scip, &objinds);
1072  SCIPfreeBufferArray(scip, &objvals);
1073  }
1074 
1075  /* add row for cutoff bound even if cutoffbound == SCIPinfinity() */
1076  lhss[nconss] = -SCIPinfinity(scip);
1077  rhss[nconss] = cutoffbound;
1078  names[nconss] = "objcutoff";
1079  lininds[nconss] = NULL;
1080  linvals[nconss] = NULL;
1081  nlininds[nconss] = 0;
1082  nquadelems[nconss] = 0;
1083  quadelems[nconss] = NULL;
1084  exprtrees[nconss] = NULL;
1085  exprvaridxs[nconss] = NULL;
1086 
1087  SCIP_CALL( SCIPallocBufferArray(scip, &lininds[nconss], nvars) ); /*lint !e866*/
1088  SCIP_CALL( SCIPallocBufferArray(scip, &linvals[nconss], nvars) ); /*lint !e866*/
1089 
1090  for( i = 0; i < nvars; ++i )
1091  {
1092  if( !SCIPisZero(scip, SCIPvarGetObj(vars[i])) )
1093  {
1094  linvals[nconss][nlininds[nconss]] = SCIPvarGetObj(vars[i]);
1095  lininds[nconss][nlininds[nconss]] = i;
1096  ++nlininds[nconss];
1097  }
1098  }
1099  ++nconss;
1100 
1101  /* add convex nonlinear rows to NLPI problem */
1102  for( i = 0; i < nnlrows; ++i )
1103  {
1104  SCIP_Bool userhs;
1105  SCIP_Bool uselhs;
1106  int k;
1107 
1108  assert(nlrows[i] != NULL);
1109 
1110  uselhs = FALSE;
1111  userhs = FALSE;
1112 
1113  /* check curvature together with constraint sides of a nonlinear row */
1114  if( SCIPnlrowGetNQuadElems(nlrows[i]) == 0 && SCIPnlrowGetExprtree(nlrows[i]) == NULL )
1115  {
1116  uselhs = TRUE;
1117  userhs = TRUE;
1118  }
1119  else
1120  {
1121  if( (!onlyconvex || SCIPnlrowGetCurvature(nlrows[i]) == SCIP_EXPRCURV_CONVEX)
1122  && !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrows[i])) )
1123  userhs = TRUE;
1124  if( (!onlyconvex || SCIPnlrowGetCurvature(nlrows[i]) == SCIP_EXPRCURV_CONCAVE)
1125  && !SCIPisInfinity(scip, SCIPnlrowGetLhs(nlrows[i])) )
1126  uselhs = TRUE;
1127  }
1128 
1129  if( !uselhs && !userhs )
1130  continue;
1131 
1132  lhss[nconss] = uselhs ? SCIPnlrowGetLhs(nlrows[i]) - SCIPnlrowGetConstant(nlrows[i]) : -SCIPinfinity(scip);
1133  rhss[nconss] = userhs ? SCIPnlrowGetRhs(nlrows[i]) - SCIPnlrowGetConstant(nlrows[i]) : SCIPinfinity(scip);
1134  names[nconss] = SCIPnlrowGetName(nlrows[i]);
1135  nlininds[nconss] = 0;
1136  lininds[nconss] = NULL;
1137  linvals[nconss] = NULL;
1138  nquadelems[nconss] = 0;
1139  quadelems[nconss] = NULL;
1140  exprtrees[nconss] = NULL;
1141  exprvaridxs[nconss] = NULL;
1142 
1143  /* copy linear part */
1144  if( SCIPnlrowGetNLinearVars(nlrows[i]) > 0 )
1145  {
1146  SCIP_VAR* var;
1147 
1148  nlininds[nconss] = SCIPnlrowGetNLinearVars(nlrows[i]);
1149 
1150  SCIP_CALL( SCIPallocBufferArray(scip, &lininds[nconss], nlininds[nconss]) ); /*lint !e866*/
1151  SCIP_CALL( SCIPallocBufferArray(scip, &linvals[nconss], nlininds[nconss]) ); /*lint !e866*/
1152 
1153  for( k = 0; k < nlininds[nconss]; ++k )
1154  {
1155  var = SCIPnlrowGetLinearVars(nlrows[i])[k];
1156  assert(var != NULL);
1157  assert(SCIPhashmapExists(var2idx, (void*)var));
1158 
1159  lininds[nconss][k] = SCIPhashmapGetImageInt(var2idx, (void*)var);
1160  assert(var == vars[lininds[nconss][k]]);
1161  linvals[nconss][k] = SCIPnlrowGetLinearCoefs(nlrows[i])[k];
1162  }
1163  }
1164 
1165  /* copy quadratic part */
1166  if( SCIPnlrowGetNQuadElems(nlrows[i]) > 0 )
1167  {
1168  SCIP_QUADELEM quadelem;
1169  SCIP_VAR* var1;
1170  SCIP_VAR* var2;
1171 
1172  nquadelems[nconss] = SCIPnlrowGetNQuadElems(nlrows[i]);
1173  SCIP_CALL( SCIPallocBufferArray(scip, &quadelems[nconss], nquadelems[nconss]) ); /*lint !e866*/
1174 
1175  for( k = 0; k < nquadelems[nconss]; ++k )
1176  {
1177  quadelem = SCIPnlrowGetQuadElems(nlrows[i])[k];
1178 
1179  var1 = SCIPnlrowGetQuadVars(nlrows[i])[quadelem.idx1];
1180  assert(var1 != NULL);
1181  assert(SCIPhashmapExists(var2idx, (void*)var1));
1182 
1183  var2 = SCIPnlrowGetQuadVars(nlrows[i])[quadelem.idx2];
1184  assert(var2 != NULL);
1185  assert(SCIPhashmapExists(var2idx, (void*)var2));
1186 
1187  quadelems[nconss][k].coef = quadelem.coef;
1188  quadelems[nconss][k].idx1 = SCIPhashmapGetImageInt(var2idx, (void*)var1);
1189  quadelems[nconss][k].idx2 = SCIPhashmapGetImageInt(var2idx, (void*)var2);
1190 
1191  /* expr.c assumes that the indices are ordered */
1192  if( quadelems[nconss][k].idx1 > quadelems[nconss][k].idx2 )
1193  {
1194  SCIPswapInts(&quadelems[nconss][k].idx1, &quadelems[nconss][k].idx2);
1195  }
1196  assert(quadelems[nconss][k].idx1 <= quadelems[nconss][k].idx2);
1197 
1198  /* update nlscore */
1199  if( nlscore != NULL )
1200  {
1201  ++nlscore[quadelems[nconss][k].idx1];
1202  if( quadelems[nconss][k].idx1 != quadelems[nconss][k].idx2 )
1203  ++nlscore[quadelems[nconss][k].idx2];
1204  }
1205  }
1206  }
1207 
1208  /* copy expression tree */
1209  if( SCIPnlrowGetExprtree(nlrows[i]) != NULL )
1210  {
1211  SCIP_VAR* var;
1212 
1213  /* note that we don't need to copy the expression tree here since only the mapping between variables in the
1214  * tree and the corresponding indices change; this mapping is stored in the exprvaridxs array
1215  */
1216  exprtrees[nconss] = SCIPnlrowGetExprtree(nlrows[i]);
1217 
1218  SCIP_CALL( SCIPallocBufferArray(scip, &exprvaridxs[nconss], SCIPexprtreeGetNVars(exprtrees[nconss])) ); /*lint !e866*/
1219 
1220  for( k = 0; k < SCIPexprtreeGetNVars(exprtrees[nconss]); ++k )
1221  {
1222  var = SCIPexprtreeGetVars(exprtrees[nconss])[k];
1223  assert(var != NULL);
1224  assert(SCIPhashmapExists(var2idx, (void*)var));
1225 
1226  exprvaridxs[nconss][k] = SCIPhashmapGetImageInt(var2idx, (void*)var);
1227 
1228  /* update nlscore */
1229  if( nlscore != NULL )
1230  ++nlscore[exprvaridxs[nconss][k]];
1231  }
1232  }
1233 
1234  /* if the row to index hash map is provided, we need to store the row index */
1235  if( nlrow2idx != NULL )
1236  {
1237  SCIP_CALL( SCIPhashmapInsertInt(nlrow2idx, nlrows[i], nconss) );
1238  }
1239 
1240  ++nconss;
1241  }
1242  assert(nconss > 0);
1243 
1244  /* pass all constraint information to nlpi */
1245  SCIP_CALL( SCIPnlpiAddConstraints(nlpi, nlpiprob, nconss, lhss, rhss, nlininds, lininds, linvals, nquadelems,
1246  quadelems, exprvaridxs, exprtrees, names) );
1247 
1248  /* free memory */
1249  for( i = nconss - 1; i > 0; --i )
1250  {
1251  if( exprtrees[i] != NULL )
1252  {
1253  assert(exprvaridxs[i] != NULL);
1254  SCIPfreeBufferArray(scip, &exprvaridxs[i]);
1255  }
1256 
1257  if( nquadelems[i] > 0 )
1258  {
1259  assert(quadelems[i] != NULL);
1260  SCIPfreeBufferArray(scip, &quadelems[i]);
1261  }
1262 
1263  if( nlininds[i] > 0 )
1264  {
1265  assert(linvals[i] != NULL);
1266  assert(lininds[i] != NULL);
1267  SCIPfreeBufferArray(scip, &linvals[i]);
1268  SCIPfreeBufferArray(scip, &lininds[i]);
1269  }
1270  }
1271  /* free row for cutoff bound even if objective is 0 */
1272  SCIPfreeBufferArray(scip, &linvals[i]);
1273  SCIPfreeBufferArray(scip, &lininds[i]);
1274 
1275  SCIPfreeBufferArray(scip, &rhss);
1276  SCIPfreeBufferArray(scip, &lhss);
1277  SCIPfreeBufferArray(scip, &names);
1278  SCIPfreeBufferArray(scip, &nlininds);
1279  SCIPfreeBufferArray(scip, &lininds);
1280  SCIPfreeBufferArray(scip, &linvals);
1281  SCIPfreeBufferArray(scip, &nquadelems);
1282  SCIPfreeBufferArray(scip, &quadelems);
1283  SCIPfreeBufferArray(scip, &exprvaridxs);
1284  SCIPfreeBufferArray(scip, &exprtrees);
1285 
1286  return SCIP_OKAY;
1287 }
1288 
1289 /** updates bounds of each variable and the cutoff row in the nlpiproblem */
1291  SCIP* scip, /**< SCIP data structure */
1292  SCIP_NLPI* nlpi, /**< interface to NLP solver */
1293  SCIP_NLPIPROBLEM* nlpiprob, /**< nlpi problem representing the convex NLP relaxation */
1294  SCIP_HASHMAP* var2nlpiidx, /**< mapping between variables and nlpi indices */
1295  SCIP_VAR** nlpivars, /**< array containing all variables of the nlpi */
1296  int nlpinvars, /**< total number of nlpi variables */
1297  SCIP_Real cutoffbound /**< new cutoff bound */
1298  )
1299 {
1300  SCIP_Real* lbs;
1301  SCIP_Real* ubs;
1302  SCIP_Real lhs;
1303  SCIP_Real rhs;
1304  int* inds;
1305  int i;
1306 
1307  SCIPdebugMsg(scip, "call SCIPupdateConvexNlpNlobbt()\n");
1308 
1309  /* update variable bounds */
1310  SCIP_CALL( SCIPallocBufferArray(scip, &lbs, nlpinvars) );
1311  SCIP_CALL( SCIPallocBufferArray(scip, &ubs, nlpinvars) );
1312  SCIP_CALL( SCIPallocBufferArray(scip, &inds, nlpinvars) );
1313 
1314  for( i = 0; i < nlpinvars; ++i )
1315  {
1316  assert(nlpivars[i] != NULL);
1317  assert(SCIPhashmapExists(var2nlpiidx, (void*)nlpivars[i]));
1318 
1319  lbs[i] = SCIPvarGetLbLocal(nlpivars[i]);
1320  ubs[i] = SCIPvarGetUbLocal(nlpivars[i]);
1321  inds[i] = SCIPhashmapGetImageInt(var2nlpiidx, (void*)nlpivars[i]);
1322  assert(inds[i] >= 0 && inds[i] < nlpinvars);
1323  }
1324 
1325  SCIP_CALL( SCIPnlpiChgVarBounds(nlpi, nlpiprob, nlpinvars, inds, lbs, ubs) );
1326 
1327  SCIPfreeBufferArray(scip, &inds);
1328  SCIPfreeBufferArray(scip, &ubs);
1329  SCIPfreeBufferArray(scip, &lbs);
1330 
1331  /* update cutoff row */
1332  lhs = -SCIPinfinity(scip);
1333  rhs = cutoffbound;
1334  i = 0;
1335 
1336  SCIP_CALL( SCIPnlpiChgConsSides(nlpi, nlpiprob, 1, &i, &lhs, &rhs) );
1337 
1338  return SCIP_OKAY;
1339 }
1340 
1341 /** adds linear rows to the NLP relaxation */
1343  SCIP* scip, /**< SCIP data structure */
1344  SCIP_NLPI* nlpi, /**< interface to NLP solver */
1345  SCIP_NLPIPROBLEM* nlpiprob, /**< nlpi problem */
1346  SCIP_HASHMAP* var2idx, /**< empty hash map to store mapping between variables and indices in nlpi
1347  * problem */
1348  SCIP_ROW** rows, /**< rows to add */
1349  int nrows /**< total number of rows to add */
1350  )
1351 {
1352  const char** names;
1353  SCIP_Real* lhss;
1354  SCIP_Real* rhss;
1355  SCIP_Real** linvals;
1356  int** lininds;
1357  int* nlininds;
1358  int i;
1359 
1360  assert(nlpi != NULL);
1361  assert(nlpiprob != NULL);
1362  assert(var2idx != NULL);
1363  assert(nrows == 0 || rows != NULL);
1364 
1365  SCIPdebugMsg(scip, "call SCIPaddConvexNlpRowsNlobbt() with %d rows\n", nrows);
1366 
1367  if( nrows <= 0 )
1368  return SCIP_OKAY;
1369 
1370  SCIP_CALL( SCIPallocBufferArray(scip, &names, nrows) );
1371  SCIP_CALL( SCIPallocBufferArray(scip, &lhss, nrows) );
1372  SCIP_CALL( SCIPallocBufferArray(scip, &rhss, nrows) );
1373  SCIP_CALL( SCIPallocBufferArray(scip, &linvals, nrows) );
1374  SCIP_CALL( SCIPallocBufferArray(scip, &lininds, nrows) );
1375  SCIP_CALL( SCIPallocBufferArray(scip, &nlininds, nrows) );
1376 
1377  for( i = 0; i < nrows; ++i )
1378  {
1379  int k;
1380 
1381  assert(rows[i] != NULL);
1382  assert(SCIProwGetNNonz(rows[i]) <= SCIPgetNVars(scip));
1383 
1384  names[i] = SCIProwGetName(rows[i]);
1385  lhss[i] = SCIProwGetLhs(rows[i]) - SCIProwGetConstant(rows[i]);
1386  rhss[i] = SCIProwGetRhs(rows[i]) - SCIProwGetConstant(rows[i]);
1387  nlininds[i] = SCIProwGetNNonz(rows[i]);
1388  linvals[i] = SCIProwGetVals(rows[i]);
1389  lininds[i] = NULL;
1390 
1391  SCIP_CALL( SCIPallocBufferArray(scip, &lininds[i], SCIProwGetNNonz(rows[i])) ); /*lint !e866*/
1392 
1393  for( k = 0; k < SCIProwGetNNonz(rows[i]); ++k )
1394  {
1395  SCIP_VAR* var;
1396 
1397  var = SCIPcolGetVar(SCIProwGetCols(rows[i])[k]);
1398  assert(var != NULL);
1399  assert(SCIPhashmapExists(var2idx, (void*)var));
1400 
1401  lininds[i][k] = SCIPhashmapGetImageInt(var2idx, (void*)var);
1402  assert(lininds[i][k] >= 0 && lininds[i][k] < SCIPgetNVars(scip));
1403  }
1404  }
1405 
1406  /* pass all linear rows to the nlpi */
1407  SCIP_CALL( SCIPnlpiAddConstraints(nlpi, nlpiprob, nrows, lhss, rhss, nlininds, lininds, linvals, NULL,
1408  NULL, NULL, NULL, names) );
1409 
1410  /* free memory */
1411  for( i = nrows - 1; i >= 0; --i )
1412  {
1413  SCIPfreeBufferArray(scip, &lininds[i]);
1414  }
1415  SCIPfreeBufferArray(scip, &nlininds);
1416  SCIPfreeBufferArray(scip, &lininds);
1417  SCIPfreeBufferArray(scip, &linvals);
1418  SCIPfreeBufferArray(scip, &rhss);
1419  SCIPfreeBufferArray(scip, &lhss);
1420  SCIPfreeBufferArray(scip, &names);
1421 
1422  return SCIP_OKAY;
1423 }
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisRelEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPnlpiAddVars(SCIP_NLPI *nlpi, SCIP_NLPIPROBLEM *problem, int nvars, const SCIP_Real *lbs, const SCIP_Real *ubs, const char **varnames)
Definition: nlpi.c:251
SCIP_EXPRTREE * SCIPnlrowGetExprtree(SCIP_NLROW *nlrow)
Definition: nlp.c:3370
SCIP_Bool SCIPisIntegral(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPquadprecSumQD(r, a, b)
Definition: dbldblarith.h:53
public methods for memory management
int SCIProwGetNNonz(SCIP_ROW *row)
Definition: lp.c:17062
internal methods for NLPI solver interfaces
SCIP_RETCODE SCIPhashmapInsertInt(SCIP_HASHMAP *hashmap, void *origin, int image)
Definition: misc.c:3131
#define SCIPquadprecSquareQ(r, a)
Definition: dbldblarith.h:59
SCIP_Real SCIProwGetConstant(SCIP_ROW *row)
Definition: lp.c:17107
#define SCIPquadprecProdQQ(r, a, b)
Definition: dbldblarith.h:57
SCIP_RETCODE SCIPnlpiChgVarBounds(SCIP_NLPI *nlpi, SCIP_NLPIPROBLEM *problem, int nvars, const int *indices, const SCIP_Real *lbs, const SCIP_Real *ubs)
Definition: nlpi.c:326
int SCIPgetNVars(SCIP *scip)
Definition: scip_prob.c:1986
#define FALSE
Definition: def.h:73
#define EPSEQ(x, y, eps)
Definition: def.h:188
SCIP_EXPORT SCIP_Real SCIPvarGetObj(SCIP_VAR *var)
Definition: var.c:17510
#define TRUE
Definition: def.h:72
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_Bool SCIPisFeasLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
public methods for problem variables
static void computeBilinEnvelope2(SCIP *scip, SCIP_Real x, SCIP_Real y, SCIP_Real mi, SCIP_Real qi, SCIP_Real mj, SCIP_Real qj, SCIP_Real *RESTRICT xi, SCIP_Real *RESTRICT yi, SCIP_Real *RESTRICT xj, SCIP_Real *RESTRICT yj, SCIP_Real *RESTRICT xcoef, SCIP_Real *RESTRICT ycoef, SCIP_Real *RESTRICT constant)
#define QUAD_SCALE(x, a)
Definition: dbldblarith.h:41
void SCIPcomputeBilinEnvelope1(SCIP *scip, SCIP_Real bilincoef, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real refpointx, SCIP_Real lby, SCIP_Real uby, SCIP_Real refpointy, SCIP_Bool overestimate, SCIP_Real xcoef, SCIP_Real ycoef, SCIP_Real constant, SCIP_Real *RESTRICT lincoefx, SCIP_Real *RESTRICT lincoefy, SCIP_Real *RESTRICT linconstant, SCIP_Bool *RESTRICT success)
void SCIPaddSquareLinearization(SCIP *scip, SCIP_Real sqrcoef, SCIP_Real refpoint, SCIP_Bool isint, SCIP_Real *lincoef, SCIP_Real *linconstant, SCIP_Bool *success)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:123
defines macros for basic operations in double-double arithmetic giving roughly twice the precision of...
#define SCIPdebugMsg
Definition: scip_message.h:69
#define SCIPquadprecDivQD(r, a, b)
Definition: dbldblarith.h:56
SCIP_VAR ** x
Definition: circlepacking.c:54
#define QUAD_TO_DBL(x)
Definition: dbldblarith.h:40
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
public methods for numerical tolerances
SCIP_RETCODE SCIPnlpiSetObjective(SCIP_NLPI *nlpi, SCIP_NLPIPROBLEM *problem, int nlins, const int *lininds, const SCIP_Real *linvals, int nquadelems, const SCIP_QUADELEM *quadelems, const int *exprvaridxs, const SCIP_EXPRTREE *exprtree, const SCIP_Real constant)
Definition: nlpi.c:301
public methods for expressions, expression trees, expression graphs, and related stuff ...
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
int SCIPhashmapGetImageInt(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3220
SCIP_Bool SCIPhashmapExists(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3362
SCIP_Real coef
Definition: type_expr.h:104
void SCIPaddBilinLinearization(SCIP *scip, SCIP_Real bilincoef, SCIP_Real refpointx, SCIP_Real refpointy, SCIP_Real *lincoefx, SCIP_Real *lincoefy, SCIP_Real *linconstant, SCIP_Bool *success)
SCIP_RETCODE SCIPnlpiAddConstraints(SCIP_NLPI *nlpi, SCIP_NLPIPROBLEM *problem, int nconss, const SCIP_Real *lhss, const SCIP_Real *rhss, const int *nlininds, int *const *lininds, SCIP_Real *const *linvals, const int *nquadelems, SCIP_QUADELEM *const *quadelems, int *const *exprvaridxs, SCIP_EXPRTREE *const *exprtrees, const char **names)
Definition: nlpi.c:269
SCIP_EXPORT const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17012
SCIP_COL ** SCIProwGetCols(SCIP_ROW *row)
Definition: lp.c:17087
SCIP_Bool SCIPisFeasEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
public methods for nonlinear functions
void SCIPcomputeBilinEnvelope2(SCIP *scip, SCIP_Real bilincoef, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real refpointx, SCIP_Real lby, SCIP_Real uby, SCIP_Real refpointy, SCIP_Bool overestimate, SCIP_Real xcoef1, SCIP_Real ycoef1, SCIP_Real constant1, SCIP_Real xcoef2, SCIP_Real ycoef2, SCIP_Real constant2, SCIP_Real *RESTRICT lincoefx, SCIP_Real *RESTRICT lincoefy, SCIP_Real *RESTRICT linconstant, SCIP_Bool *RESTRICT success)
SCIP_VAR ** SCIPgetVars(SCIP *scip)
Definition: scip_prob.c:1941
SCIP_Real * SCIProwGetVals(SCIP_ROW *row)
Definition: lp.c:17097
#define QUAD(x)
Definition: dbldblarith.h:38
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
int SCIPnlrowGetNQuadElems(SCIP_NLROW *nlrow)
Definition: nlp.c:3329
#define NULL
Definition: lpi_spx1.cpp:155
#define REALABS(x)
Definition: def.h:187
SCIP_RETCODE SCIPcreateNlpiProb(SCIP *scip, SCIP_NLPI *nlpi, SCIP_NLROW **nlrows, int nnlrows, SCIP_NLPIPROBLEM *nlpiprob, SCIP_HASHMAP *var2idx, SCIP_HASHMAP *nlrow2idx, SCIP_Real *nlscore, SCIP_Real cutoffbound, SCIP_Bool setobj, SCIP_Bool onlyconvex)
int SCIPexprtreeGetNVars(SCIP_EXPRTREE *tree)
Definition: expr.c:8614
#define SCIP_CALL(x)
Definition: def.h:364
void SCIPswapInts(int *value1, int *value2)
Definition: misc.c:10192
int SCIPnlrowGetNLinearVars(SCIP_NLROW *nlrow)
Definition: nlp.c:3252
SCIP_QUADELEM * SCIPnlrowGetQuadElems(SCIP_NLROW *nlrow)
Definition: nlp.c:3339
#define SCIPquadprecProdDD(r, a, b)
Definition: dbldblarith.h:49
public methods for NLP management
#define SCIPquadprecSqrtQ(r, a)
Definition: dbldblarith.h:62
SCIP_EXPRCURV SCIPnlrowGetCurvature(SCIP_NLROW *nlrow)
Definition: nlp.c:3400
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:111
SCIP_Real SCIPinfinity(SCIP *scip)
public data structures and miscellaneous methods
SCIP_Real SCIPnlrowGetRhs(SCIP_NLROW *nlrow)
Definition: nlp.c:3390
#define SCIP_Bool
Definition: def.h:70
#define MAX(x, y)
Definition: tclique_def.h:83
SCIP_Real SCIProwGetLhs(SCIP_ROW *row)
Definition: lp.c:17141
public methods for LP management
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
SCIP_VAR ** SCIPexprtreeGetVars(SCIP_EXPRTREE *tree)
Definition: nlp.c:103
#define SCIPquadprecProdQD(r, a, b)
Definition: dbldblarith.h:54
SCIP_Real SCIPnlrowGetLhs(SCIP_NLROW *nlrow)
Definition: nlp.c:3380
SCIP_Bool SCIPisFeasGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Real * SCIPnlrowGetLinearCoefs(SCIP_NLROW *nlrow)
Definition: nlp.c:3272
SCIP_EXPORT SCIP_Real SCIPvarGetLbLocal(SCIP_VAR *var)
Definition: var.c:17718
#define SCIPquadprecSumQQ(r, a, b)
Definition: dbldblarith.h:58
SCIP_VAR * SCIPcolGetVar(SCIP_COL *col)
Definition: lp.c:16901
SCIP_Real SCIPfloor(SCIP *scip, SCIP_Real val)
SCIP_EXPORT SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17728
#define SCIPquadprecDivDD(r, a, b)
Definition: dbldblarith.h:52
SCIP_RETCODE SCIPnlpiChgConsSides(SCIP_NLPI *nlpi, SCIP_NLPIPROBLEM *problem, int nconss, const int *indices, const SCIP_Real *lhss, const SCIP_Real *rhss)
Definition: nlpi.c:344
SCIP_Bool SCIPisFeasGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPupdateNlpiProb(SCIP *scip, SCIP_NLPI *nlpi, SCIP_NLPIPROBLEM *nlpiprob, SCIP_HASHMAP *var2nlpiidx, SCIP_VAR **nlpivars, int nlpinvars, SCIP_Real cutoffbound)
type definitions for expressions and expression trees
public methods for message output
const char * SCIPnlrowGetName(SCIP_NLROW *nlrow)
Definition: nlp.c:3419
#define SCIP_Real
Definition: def.h:163
SCIP_VAR ** y
Definition: circlepacking.c:55
const char * SCIProwGetName(SCIP_ROW *row)
Definition: lp.c:17200
public methods for message handling
#define SCIP_INVALID
Definition: def.h:183
#define SCIPquadprecSumDD(r, a, b)
Definition: dbldblarith.h:51
SCIP_VAR ** SCIPnlrowGetQuadVars(SCIP_NLROW *nlrow)
Definition: nlp.c:3292
SCIP_Real SCIPnlrowGetConstant(SCIP_NLROW *nlrow)
Definition: nlp.c:3242
SCIP_Real SCIProwGetRhs(SCIP_ROW *row)
Definition: lp.c:17151
SCIP_VAR ** SCIPnlrowGetLinearVars(SCIP_NLROW *nlrow)
Definition: nlp.c:3262
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:122
void SCIPaddBilinMcCormick(SCIP *scip, SCIP_Real bilincoef, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real refpointx, SCIP_Real lby, SCIP_Real uby, SCIP_Real refpointy, SCIP_Bool overestimate, SCIP_Real *lincoefx, SCIP_Real *lincoefy, SCIP_Real *linconstant, SCIP_Bool *success)
public methods for global and local (sub)problems
#define SCIPquadprecDivQQ(r, a, b)
Definition: dbldblarith.h:60
void SCIPaddSquareSecant(SCIP *scip, SCIP_Real sqrcoef, SCIP_Real lb, SCIP_Real ub, SCIP_Real refpoint, SCIP_Real *lincoef, SCIP_Real *linconstant, SCIP_Bool *success)
SCIP_RETCODE SCIPaddNlpiProbRows(SCIP *scip, SCIP_NLPI *nlpi, SCIP_NLPIPROBLEM *nlpiprob, SCIP_HASHMAP *var2idx, SCIP_ROW **rows, int nrows)
memory allocation routines