1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
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   nlhdlr_quotient.c
26   	 * @ingroup DEFPLUGINS_NLHDLR
27   	 * @brief  quotient nonlinear handler
28   	 * @author Benjamin Mueller
29   	 * @author Fabian Wegscheider
30   	 *
31   	 * @todo implement INITSEPA
32   	 * @todo use the convex envelope for x/y described in Tawarmalani and Sahinidis (2002) if y has a finite upper bound
33   	 */
34   	
35   	#include <string.h>
36   	
37   	#include "scip/nlhdlr_quotient.h"
38   	#include "scip/cons_nonlinear.h"
39   	#include "scip/pub_misc_rowprep.h"
40   	#include "scip/nlhdlr.h"
41   	#include "scip/nlhdlr_bilinear.h"
42   	
43   	/* fundamental nonlinear handler properties */
44   	#define NLHDLR_NAME         "quotient"
45   	#define NLHDLR_DESC         "nonlinear handler for quotient expressions"
46   	#define NLHDLR_DETECTPRIORITY   20
47   	#define NLHDLR_ENFOPRIORITY     20
48   	
49   	/** translate from one value of infinity to another
50   	 *
51   	 *  if val is &ge; infty1, then give infty2, else give val
52   	 */
53   	#define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))
54   	
55   	/*lint -e666*/
56   	
57   	/*
58   	 * Data structures
59   	 */
60   	
61   	/** nonlinear handler expression data */
62   	struct SCIP_NlhdlrExprData
63   	{
64   	   SCIP_EXPR*            numexpr;            /**< expression of the numerator */
65   	   SCIP_Real             numcoef;            /**< coefficient of the numerator */
66   	   SCIP_Real             numconst;           /**< constant of the numerator */
67   	   SCIP_EXPR*            denomexpr;          /**< expression of the denominator */
68   	   SCIP_Real             denomcoef;          /**< coefficient of the denominator */
69   	   SCIP_Real             denomconst;         /**< constant of the denominator */
70   	   SCIP_Real             constant;           /**< constant */
71   	};
72   	
73   	/*
74   	 * Local methods
75   	 */
76   	
77   	/** helper method to create nonlinear handler expression data */
78   	static
79   	SCIP_RETCODE exprdataCreate(
80   	   SCIP*                 scip,               /**< SCIP data structure */
81   	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< nonlinear handler expression data */
82   	   SCIP_EXPR*            numexpr,            /**< expression of the numerator */
83   	   SCIP_Real             numcoef,            /**< coefficient of the numerator */
84   	   SCIP_Real             numconst,           /**< constant of the numerator */
85   	   SCIP_EXPR*            denomexpr,          /**< expression of the denominator */
86   	   SCIP_Real             denomcoef,          /**< coefficient of the denominator */
87   	   SCIP_Real             denomconst,         /**< constant of the denominator */
88   	   SCIP_Real             constant            /**< constant */
89   	   )
90   	{
91   	   assert(nlhdlrexprdata != NULL);
92   	   assert(numexpr != NULL);
93   	   assert(denomexpr != NULL);
94   	   assert(!SCIPisZero(scip, numcoef));
95   	   assert(!SCIPisZero(scip, denomcoef));
96   	
97   	   /* allocate memory */
98   	   SCIP_CALL( SCIPallocBlockMemory(scip, nlhdlrexprdata) );
99   	
100  	   /* store values */
101  	   (*nlhdlrexprdata)->numexpr = numexpr;
102  	   (*nlhdlrexprdata)->numcoef = numcoef;
103  	   (*nlhdlrexprdata)->numconst = numconst;
104  	   (*nlhdlrexprdata)->denomexpr = denomexpr;
105  	   (*nlhdlrexprdata)->denomcoef = denomcoef;
106  	   (*nlhdlrexprdata)->denomconst = denomconst;
107  	   (*nlhdlrexprdata)->constant = constant;
108  	
109  	   /* capture expressions */
110  	   SCIPcaptureExpr(numexpr);
111  	   SCIPcaptureExpr(denomexpr);
112  	
113  	   return SCIP_OKAY;
114  	}
115  	
116  	/** helper method to free nonlinear handler expression data */
117  	static
118  	SCIP_RETCODE exprdataFree(
119  	   SCIP*                 scip,               /**< SCIP data structure */
120  	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata      /**< nonlinear handler expression data */
121  	   )
122  	{
123  	   assert(nlhdlrexprdata != NULL);
124  	   assert(*nlhdlrexprdata != NULL);
125  	   assert((*nlhdlrexprdata)->numexpr != NULL);
126  	   assert((*nlhdlrexprdata)->denomexpr != NULL);
127  	
128  	   /* release expressions */
129  	   SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->denomexpr) );
130  	   SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->numexpr) );
131  	
132  	   /* free expression data of nonlinear handler */
133  	   SCIPfreeBlockMemory(scip, nlhdlrexprdata);
134  	
135  	   return SCIP_OKAY;
136  	}
137  	
138  	/** helper method to transform an expression g(x) into a*f(x) + b */
139  	static
140  	void transformExpr(
141  	   SCIP*                 scip,               /**< SCIP data structure */
142  	   SCIP_EXPR*            expr,               /**< expression */
143  	   SCIP_EXPR**           target,             /**< pointer to store the expression f(x) */
144  	   SCIP_Real*            coef,               /**< pointer to store the coefficient */
145  	   SCIP_Real*            constant            /**< pointer to store the constant */
146  	   )
147  	{
148  	   assert(expr != NULL);
149  	   assert(target != NULL);
150  	   assert(coef != NULL);
151  	   assert(constant != NULL);
152  	
153  	   /* expression is a sum with one child */
154  	   if( SCIPisExprSum(scip, expr) && SCIPexprGetNChildren(expr) == 1 )
155  	   {
156  	      *target = SCIPexprGetChildren(expr)[0];
157  	      *coef = SCIPgetCoefsExprSum(expr)[0];
158  	      *constant = SCIPgetConstantExprSum(expr);
159  	   }
160  	   else /* otherwise return 1 * f(x) + 0 */
161  	   {
162  	      *target = expr;
163  	      *coef = 1.0;
164  	      *constant = 0.0;
165  	   }
166  	}
167  	
168  	/** helper method to detect an expression of the form (a*x + b) / (c*y + d) + e
169  	 *
170  	 * Due to the expansion of products, there are two types of expressions that can be detected:
171  	 *
172  	 * 1. prod(f(x), pow(g(y),-1))
173  	 * 2. sum(prod(f(x),pow(g(y),-1)), pow(g(y),-1))
174  	 *
175  	 * @todo At the moment quotients like xy / z are not detected, because they are turned into a product expression
176  	 * with three children, i.e., x * y * (1 / z).
177  	 */
178  	static
179  	SCIP_RETCODE detectExpr(
180  	   SCIP*                 scip,               /**< SCIP data structure */
181  	   SCIP_EXPR*            expr,               /**< expression */
182  	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< pointer to store nonlinear handler expression data */
183  	   SCIP_Bool*            success             /**< pointer to store whether nonlinear handler should be called for this expression */
184  	   )
185  	{
186  	   SCIP_EXPR** children;
187  	   SCIP_EXPR* denomexpr = NULL;
188  	   SCIP_EXPR* numexpr = NULL;
189  	   SCIP_EXPR* xexpr = NULL;
190  	   SCIP_EXPR* yexpr = NULL;
191  	   SCIP_Real a, b, c, d, e;
192  	   SCIP_Real nomfac = 1.0;
193  	   SCIP_Real numconst = 0.0;
194  	
195  	   assert(scip != NULL);
196  	   assert(expr != NULL);
197  	
198  	   *success = FALSE;
199  	   a = 0.0;
200  	   b = 0.0;
201  	   c = 0.0;
202  	   d = 0.0;
203  	   e = 0.0;
204  	
205  	   /* possible structures only have two children */
206  	   if( SCIPexprGetNChildren(expr) != 2 )
207  	      return SCIP_OKAY;
208  	
209  	   /* expression must be either a product or a sum */
210  	   if( !SCIPisExprProduct(scip, expr) && !SCIPisExprSum(scip, expr) )
211  	      return SCIP_OKAY;
212  	
213  	   children = SCIPexprGetChildren(expr);
214  	   assert(children != NULL);
215  	
216  	   /* case: prod(f(x), pow(g(y),-1)) */
217  	   if( SCIPisExprProduct(scip, expr) )
218  	   {
219  	      if( SCIPisExprPower(scip, children[0]) && SCIPgetExponentExprPow(children[0]) == -1.0 )  /*lint !e777*/
220  	      {
221  	         denomexpr = SCIPexprGetChildren(children[0])[0];
222  	         numexpr = children[1];
223  	      }
224  	      else if( SCIPisExprPower(scip, children[1]) && SCIPgetExponentExprPow(children[1]) == -1.0 )  /*lint !e777*/
225  	      {
226  	         denomexpr = SCIPexprGetChildren(children[1])[0];
227  	         numexpr = children[0];
228  	      }
229  	
230  	      /* remember to scale the numerator by the coefficient stored in the product expression */
231  	      nomfac = SCIPgetCoefExprProduct(expr);
232  	   }
233  	   /* case: sum(prod(f(x),pow(g(y),-1)), pow(g(y),-1)) */
234  	   else
235  	   {
236  	      SCIP_Real* sumcoefs;
237  	
238  	      assert(SCIPisExprSum(scip, expr));
239  	      sumcoefs = SCIPgetCoefsExprSum(expr);
240  	
241  	      /* children[0] is 1/g(y) and children[1] is a product of f(x) and 1/g(y) */
242  	      if( SCIPisExprPower(scip, children[0]) && SCIPgetExponentExprPow(children[0]) == -1.0
243  	         && SCIPisExprProduct(scip, children[1]) && SCIPexprGetNChildren(children[1]) == 2 )  /* lint !e777 */
244  	      {
245  	         SCIP_Real prodcoef = SCIPgetCoefExprProduct(children[1]);
246  	
247  	         if( children[0] == SCIPexprGetChildren(children[1])[0] )
248  	         {
249  	            denomexpr = SCIPexprGetChildren(children[0])[0];
250  	            numexpr = SCIPexprGetChildren(children[1])[1];
251  	         }
252  	         else if( children[0] == SCIPexprGetChildren(children[1])[1] )
253  	         {
254  	            denomexpr = SCIPexprGetChildren(children[0])[0];
255  	            numexpr = SCIPexprGetChildren(children[1])[0];
256  	         }
257  	
258  	         /* remember scalar and constant for numerator */
259  	         nomfac = sumcoefs[1] * prodcoef;
260  	         numconst = sumcoefs[0];
261  	      }
262  	      /* children[1] is 1/g(y) and children[0] is a product of f(x) and 1/g(y) */
263  	      else if( SCIPisExprPower(scip, children[1]) && SCIPgetExponentExprPow(children[1]) == -1.0
264  	         && SCIPisExprProduct(scip, children[0]) && SCIPexprGetNChildren(children[0]) == 2 )  /* lint !e777 */
265  	      {
266  	         SCIP_Real prodcoef = SCIPgetCoefExprProduct(children[0]);
267  	
268  	         if( children[1] == SCIPexprGetChildren(children[0])[0] )
269  	         {
270  	            denomexpr = SCIPexprGetChildren(children[1])[0];
271  	            numexpr = SCIPexprGetChildren(children[0])[1];
272  	         }
273  	         else if( children[1] == SCIPexprGetChildren(children[0])[1] )
274  	         {
275  	            denomexpr = SCIPexprGetChildren(children[1])[0];
276  	            numexpr = SCIPexprGetChildren(children[0])[0];
277  	         }
278  	
279  	         /* remember scalar and constant for numerator */
280  	         nomfac = sumcoefs[0] * prodcoef;
281  	         numconst = sumcoefs[1];
282  	      }
283  	
284  	      /* remember the constant of the sum expression */
285  	      e = SCIPgetConstantExprSum(expr);
286  	   }
287  	
288  	   if( denomexpr != NULL && numexpr != NULL )
289  	   {
290  	      /* transform numerator and denominator to detect structures like (a * f(x) + b) / (c * f(x) + d) */
291  	      transformExpr(scip, numexpr, &xexpr, &a, &b);
292  	      transformExpr(scip, denomexpr, &yexpr, &c, &d);
293  	
294  	      SCIPdebugMsg(scip, "detected numerator (%g * %p + %g) and denominator (%g * %p + %g)\n", a, (void*)xexpr, b,
295  	         c, (void*)yexpr, d);
296  	
297  	      /* detection is only be successful if the expression of the numerator an denominator are the same
298  	       * (so boundtightening can be stronger than default) or we are going to provide estimators (there will be an auxvar)
299  	       */
300  	      *success = (xexpr == yexpr) || (SCIPgetExprNAuxvarUsesNonlinear(expr) > 0);
301  	
302  	#ifdef SCIP_DEBUG
303  	      SCIPinfoMessage(scip, NULL, "Expression for numerator: ");
304  	      SCIP_CALL( SCIPprintExpr(scip, xexpr, NULL) );
305  	      SCIPinfoMessage(scip, NULL, "\nExpression for denominator: ");
306  	      SCIP_CALL( SCIPprintExpr(scip, yexpr, NULL) );
307  	      SCIPinfoMessage(scip, NULL, "\n");
308  	#endif
309  	   }
310  	
311  	   /* register usage of xexpr and yexpr
312  	    * create nonlinear handler expression data
313  	    */
314  	   if( *success )
315  	   {
316  	      assert(xexpr != NULL);
317  	      assert(xexpr != NULL);
318  	      assert(a != 0.0);
319  	      assert(c != 0.0);
320  	
321  	      assert(SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 || xexpr == yexpr);
322  	
323  	      /* request auxiliary variables for xexpr and yexpr if we will estimate
324  	       * mark that the bounds of the expression are important to construct the estimators
325  	       *   (TODO check the curvature of the univariate quotient, as bounds may actually not be used)
326  	       * if univariate, then we also do inteval and reverseprop, so mark that the activities will be used for inteval
327  	       */
328  	      SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, xexpr,
329  	         SCIPgetExprNAuxvarUsesNonlinear(expr) > 0,
330  	         xexpr == yexpr,
331  	         SCIPgetExprNAuxvarUsesNonlinear(expr) > 0,
332  	         SCIPgetExprNAuxvarUsesNonlinear(expr) > 0) );
333  	
334  	      if( xexpr != yexpr && SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 )
335  	      {
336  	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, yexpr, TRUE, FALSE, TRUE, TRUE) );
337  	      }
338  	
339  	      a = nomfac * a;
340  	      b = nomfac * b + numconst;
341  	
342  	      SCIPdebug( SCIP_CALL( SCIPprintExpr(scip, expr, NULL) ); )
343  	      SCIPdebug( SCIPinfoMessage(scip, NULL, "\n") );
344  	      SCIPdebugMsg(scip, "detected quotient expression (%g * %p + %g) / (%g * %p + %g) + %g\n", a, (void*)xexpr,
345  	         b, c, (void*)yexpr, d, e);
346  	      SCIP_CALL( exprdataCreate(scip, nlhdlrexprdata, xexpr, a, b, yexpr, c, d, e) );
347  	   }
348  	
349  	   return SCIP_OKAY;
350  	}
351  	
352  	/** helper method to compute interval for (a x + b) / (c x + d) + e */
353  	static
354  	SCIP_INTERVAL intEvalQuotient(
355  	   SCIP*                 scip,               /**< SCIP data structure */
356  	   SCIP_INTERVAL         bnds,               /**< bounds on x */
357  	   SCIP_Real             a,                  /**< coefficient in numerator */
358  	   SCIP_Real             b,                  /**< constant in numerator */
359  	   SCIP_Real             c,                  /**< coefficient in denominator */
360  	   SCIP_Real             d,                  /**< constant in denominator */
361  	   SCIP_Real             e                   /**< constant */
362  	   )
363  	{
364  	   SCIP_INTERVAL result;
365  	   SCIP_INTERVAL denominterval;
366  	   SCIP_INTERVAL numinterval;
367  	   int i;
368  	
369  	   assert(scip != NULL);
370  	
371  	   /* return empty interval if the domain of x is empty */
372  	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bnds) )
373  	   {
374  	      SCIPintervalSetEmpty(&result);
375  	      return result;
376  	   }
377  	
378  	   /* compute bounds for denominator */
379  	   SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominterval, bnds, c);
380  	   SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, d);
381  	
382  	   /* there is no useful interval if 0 is in the interior of the interval of the denominator */
383  	   if( SCIPintervalGetInf(denominterval) < 0.0 && SCIPintervalGetSup(denominterval) > 0.0 )
384  	   {
385  	      SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &result);
386  	      return result;
387  	   }
388  	
389  	   /* a d = b c implies that f(x) = b / d + e, i.e., f is constant */
390  	   if( a*d - b*c == 0.0 )
391  	   {
392  	      SCIPintervalSet(&result, b / d + e);
393  	      return result;
394  	   }
395  	
396  	   /*
397  	    * evaluate for [x.inf,x.inf] and [x.sup,x.sup] independently
398  	    */
399  	   SCIPintervalSetEmpty(&result);
400  	
401  	   for( i = 0; i < 2; ++i )
402  	   {
403  	      SCIP_INTERVAL quotinterval;
404  	      SCIP_Real val = (i == 0) ? bnds.inf : bnds.sup;
405  	
406  	      /* set the resulting interval to a / c if the bounds is infinite */
407  	      if( SCIPisInfinity(scip, REALABS(val)) )
408  	      {
409  	         SCIPintervalSet(&quotinterval, a);
410  	         SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, &quotinterval, quotinterval, c);
411  	      }
412  	      else
413  	      {
414  	         /* a x' + b */
415  	         SCIPintervalSet(&numinterval, val);
416  	         SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &numinterval, numinterval, a);
417  	         SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &numinterval, numinterval, b);
418  	
419  	         /* c x' + d */
420  	         SCIPintervalSet(&denominterval, val);
421  	         SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, c);
422  	         SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, d);
423  	
424  	         /* (a x' + b) / (c x' + d) + e */
425  	         SCIPintervalDiv(SCIP_INTERVAL_INFINITY, &quotinterval, numinterval, denominterval);
426  	         SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &quotinterval, quotinterval, e);
427  	      }
428  	
429  	      /* unify with the resulting interval */
430  	      SCIPintervalUnify(&result, result, quotinterval);
431  	   }
432  	
433  	   return result;
434  	}
435  	
436  	/** helper method to compute reverse propagation for (a x + b) / (c x + d) + e */
437  	static
438  	SCIP_INTERVAL reversepropQuotient(
439  	   SCIP_INTERVAL         bnds,               /**< bounds on (a x + b) / (c x + d) + e */
440  	   SCIP_Real             a,                  /**< coefficient in numerator */
441  	   SCIP_Real             b,                  /**< constant in numerator */
442  	   SCIP_Real             c,                  /**< coefficient in denominator */
443  	   SCIP_Real             d,                  /**< constant in denominator */
444  	   SCIP_Real             e                   /**< constant */
445  	   )
446  	{
447  	   SCIP_INTERVAL result;
448  	   int i;
449  	
450  	   SCIPintervalSetEmpty(&result);
451  	
452  	   /* return empty interval if the domain of the expression is empty */
453  	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bnds) )
454  	      return result;
455  	
456  	   /* substract constant from bounds of the expression */
457  	   SCIPintervalSubScalar(SCIP_INTERVAL_INFINITY, &bnds, bnds, e);
458  	
459  	   /* if the expression is constant or the limit lies inside the domain, nothing can be propagated */
460  	   if( a*d - b*c == 0.0 || (bnds.inf < a / c && bnds.sup > a / c) )
461  	   {
462  	      SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &result);
463  	      return result;
464  	   }
465  	
466  	   /* compute bounds for [x.inf,x.inf] and [x.sup,x.sup] independently */
467  	   for( i = 0; i < 2; ++i )
468  	   {
469  	      SCIP_INTERVAL denominator;
470  	      SCIP_INTERVAL numerator;
471  	      SCIP_INTERVAL quotient;
472  	      SCIP_Real val = (i == 0) ? bnds.inf : bnds.sup;
473  	
474  	      /* (d * x' - b) */
475  	      SCIPintervalSet(&numerator, d);
476  	      SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &numerator, numerator, val);
477  	      SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &numerator, numerator, -b);
478  	
479  	      /* (a - c * x') */
480  	      SCIPintervalSet(&denominator, -c);
481  	      SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominator, denominator, val);
482  	      SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominator, denominator, a);
483  	
484  	      /* (d * x' - b) / (a - c * x') */
485  	      SCIPintervalDiv(SCIP_INTERVAL_INFINITY, &quotient, numerator, denominator);
486  	
487  	      /* unify with the resulting interval */
488  	      SCIPintervalUnify(&result, result, quotient);
489  	   }
490  	
491  	   return result;
492  	}
493  	
494  	/** adds data to given rowprep; the generated estimator is always locally valid
495  	 *
496  	 *  @note the constant is moved to the left- or right-hand side
497  	 *  @note other than the name of this function may indicate, it does not create a rowprep
498  	 */
499  	static
500  	SCIP_RETCODE createRowprep(
501  	   SCIP*                 scip,               /**< SCIP data structure */
502  	   SCIP_ROWPREP*         rowprep,            /**< a rowprep where to store the estimator */
503  	   SCIP_VAR**            vars,               /**< variables */
504  	   SCIP_Real*            coefs,              /**< coefficients */
505  	   SCIP_Real             constant,           /**< constant */
506  	   int                   nlinvars            /**< total number of variables */
507  	   )
508  	{
509  	   assert(scip != NULL);
510  	   assert(rowprep != NULL);
511  	   assert(coefs != NULL);
512  	   assert(vars != NULL);
513  	
514  	   /* create rowprep */
515  	   SCIProwprepAddSide(rowprep, -constant);
516  	   SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, nlinvars + 1) );
517  	
518  	   /* add coefficients */
519  	   SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, nlinvars, vars, coefs) );
520  	
521  	   return SCIP_OKAY;
522  	}
523  	
524  	/** computes an estimator at a given point for the univariate case (ax + b) / (cx + d) + e
525  	 *
526  	 *  Depending on the reference point, the estimator is a tangent or a secant on the graph.
527  	 *  It depends on whether we are under- or overestimating, whether we are on the left or
528  	 *  on the right side of the singularity at -d/c, and whether it is the monotone increasing
529  	 *  (ad - bc > 0) or decreasing part (ad - bc < 0). Together, there are 8 cases:
530  	 *
531  	 *  - mon. incr. + overestimate + left hand side  -->  secant
532  	 *  - mon. incr. + overestimate + right hand side -->  tangent
533  	 *  - mon. incr. + understimate + left hand side  -->  tangent
534  	 *  - mon. incr. + understimate + right hand side -->  secant
535  	 *  - mon. decr. + overestimate + left hand side  -->  tangent
536  	 *  - mon. decr. + overestimate + right hand side -->  secant
537  	 *  - mon. decr. + understimate + left hand side  -->  secant
538  	 *  - mon. decr. + understimate + right hand side -->  tangent
539  	 */
540  	static
541  	SCIP_RETCODE estimateUnivariate(
542  	   SCIP*                 scip,               /**< SCIP data structure */
543  	   SCIP_Real             lbx,                /**< local lower bound of x */
544  	   SCIP_Real             ubx,                /**< local upper bound of x */
545  	   SCIP_Real             gllbx,              /**< global lower bound of x */
546  	   SCIP_Real             glubx,              /**< global upper bound of x */
547  	   SCIP_Real             solx,               /**< solution value of x */
548  	   SCIP_Real             a,                  /**< coefficient in numerator */
549  	   SCIP_Real             b,                  /**< constant in numerator */
550  	   SCIP_Real             c,                  /**< coefficient in denominator */
551  	   SCIP_Real             d,                  /**< constant in denominator */
552  	   SCIP_Real             e,                  /**< constant */
553  	   SCIP_Real*            coef,               /**< pointer to store the coefficient */
554  	   SCIP_Real*            constant,           /**< pointer to store the constant */
555  	   SCIP_Bool             overestimate,       /**< whether the expression should be overestimated */
556  	   SCIP_Bool*            local,              /**< pointer to store whether the estimate is locally valid */
557  	   SCIP_Bool*            branchinguseful,    /**< pointer to store whether branching on the expression would improve the estimator */
558  	   SCIP_Bool*            success             /**< buffer to store whether separation was successful */
559  	   )
560  	{
561  	   SCIP_Real singularity;
562  	   SCIP_Bool isinleftpart;
563  	   SCIP_Bool monincreasing;
564  	
565  	   assert(lbx <= solx && solx <= ubx);
566  	   assert(coef != NULL);
567  	   assert(constant != NULL);
568  	   assert(local != NULL);
569  	   assert(branchinguseful != NULL);
570  	   assert(success != NULL);
571  	
572  	   *branchinguseful = TRUE;
573  	   *success = FALSE;
574  	   *coef = 0.0;
575  	   *constant = 0.0;
576  	   singularity = -d / c;
577  	
578  	   /* estimate is globally valid if local and global bounds are equal */
579  	   *local = gllbx != lbx || glubx != ubx; /*lint !e777*/
580  	
581  	   /* if 0 is in the denom interval, estimation is not possible */
582  	   if( SCIPisLE(scip, lbx, singularity) && SCIPisGE(scip, ubx, singularity) )
583  	      return SCIP_OKAY;
584  	
585  	   isinleftpart = (ubx < singularity);
586  	   monincreasing = (a * d - b * c > 0.0);
587  	
588  	   /* this encodes the 8 cases explained above */
589  	   if( monincreasing == (overestimate == isinleftpart) )
590  	   {
591  	      SCIP_Real lbeval;
592  	      SCIP_Real ubeval;
593  	
594  	      /* if one of the bounds is infinite, secant cannot be computed */
595  	      if( SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip, ubx) )
596  	         return SCIP_OKAY;
597  	
598  	      lbeval = (a * lbx + b) / (c * lbx + d) + e;
599  	      ubeval = (a * ubx + b) / (c * ubx + d) + e;
600  	
601  	      /* compute coefficient and constant of linear estimator */
602  	      *coef = (ubeval - lbeval) / (ubx - lbx);
603  	      *constant = ubeval - (*coef) * ubx;
604  	   }
605  	   else
606  	   {
607  	      SCIP_Real soleval;
608  	
609  	      soleval = (a * solx + b) / (c * solx + d) + e;
610  	
611  	      /* compute coefficient and constant of linear estimator */
612  	      *coef = (a * d - b * c) / SQR(d + c * solx);
613  	      *constant = soleval - (*coef) * solx;
614  	
615  	      /* gradient cuts are globally valid if the singularity is not in [gllbx,glubx] */
616  	      *local = SCIPisLE(scip, gllbx, singularity) && SCIPisGE(scip, glubx, singularity);
617  	
618  	      /* branching will not improve the convexification via tangent cuts */
619  	      *branchinguseful = FALSE;
620  	   }
621  	
622  	   /* avoid huge values in the cut */
623  	   if( SCIPisHugeValue(scip, REALABS(*coef)) || SCIPisHugeValue(scip, REALABS(*constant)) )
624  	      return SCIP_OKAY;
625  	
626  	   *success = TRUE;
627  	
628  	   return SCIP_OKAY;
629  	}
630  	
631  	/** helper method to compute estimator for the univariate case; the estimator is stored in a given rowprep */
632  	static
633  	SCIP_RETCODE estimateUnivariateQuotient(
634  	   SCIP*                 scip,               /**< SCIP data structure */
635  	   SCIP_SOL*             sol,                /**< solution point (or NULL for the LP solution) */
636  	   SCIP_EXPR*            xexpr,              /**< argument expression */
637  	   SCIP_Real             a,                  /**< coefficient in numerator */
638  	   SCIP_Real             b,                  /**< constant in numerator */
639  	   SCIP_Real             c,                  /**< coefficient in denominator */
640  	   SCIP_Real             d,                  /**< constant in denominator */
641  	   SCIP_Real             e,                  /**< constant */
642  	   SCIP_Bool             overestimate,       /**< whether the expression should be overestimated */
643  	   SCIP_ROWPREP*         rowprep,            /**< a rowprep where to store the estimator */
644  	   SCIP_Bool*            branchinguseful,    /**< pointer to store whether branching on the expression would improve the estimator */
645  	   SCIP_Bool*            success             /**< buffer to store whether separation was successful */
646  	   )
647  	{
648  	   SCIP_VAR* x;
649  	   SCIP_Real constant;
650  	   SCIP_Real coef;
651  	   SCIP_Real gllbx;
652  	   SCIP_Real glubx;
653  	   SCIP_Real lbx;
654  	   SCIP_Real ubx;
655  	   SCIP_Real solx;
656  	   SCIP_Bool local;
657  	   SCIP_INTERVAL bnd;
658  	
659  	   assert(rowprep != NULL);
660  	   assert(branchinguseful != NULL);
661  	   assert(success != NULL);
662  	
663  	   x = SCIPgetExprAuxVarNonlinear(xexpr);
664  	
665  	   /* get local bounds on xexpr */
666  	   SCIPintervalSetBounds(&bnd,
667  	      -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(x)),
668  	       infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY,  SCIPvarGetUbGlobal(x)));
669  	   SCIP_CALL( SCIPevalExprActivity(scip, xexpr) );
670  	   SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(xexpr));
671  	   lbx = bnd.inf;
672  	   ubx = bnd.sup;
673  	
674  	   /* check whether variable has been fixed or has empty interval */
675  	   if( SCIPisEQ(scip, lbx, ubx) || ubx < lbx )
676  	   {
677  	      *success = FALSE;
678  	      return SCIP_OKAY;
679  	   }
680  	
681  	   /* get global variable bounds */
682  	   gllbx = SCIPvarGetLbGlobal(x);
683  	   glubx = SCIPvarGetUbGlobal(x);
684  	
685  	   /* get and adjust solution value */
686  	   solx = SCIPgetSolVal(scip, sol, x);
687  	   solx = MIN(MAX(solx, lbx), ubx);
688  	
689  	   /* compute an estimator */
690  	   SCIP_CALL( estimateUnivariate(scip, lbx, ubx, gllbx, glubx, solx, a, b, c, d, e, &coef, &constant, overestimate, &local, branchinguseful, success) );
691  	
692  	   /* add estimator to rowprep, if successful */
693  	   if( *success )
694  	   {
695  	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_%lld", SCIPvarGetName(x), SCIPgetNLPs(scip));
696  	      SCIP_CALL( createRowprep(scip, rowprep, &x, &coef, constant, 1) );
697  	      SCIProwprepSetLocal(rowprep, local);
698  	   }
699  	
700  	   return SCIP_OKAY;
701  	}
702  	
703  	/** helper method to compute a gradient cut for
704  	 *  \f[
705  	 *     h^c(x,y) := \frac{1}{y} \left(\frac{x + \sqrt{\text{lbx}\cdot\text{ubx}}}{\sqrt{\text{lbx}} + \sqrt{\text{ubx}}}\right)^2
706  	 *  \f]
707  	 *  at a given reference point
708  	 *
709  	 *  See Zamora and Grossmann (1988) for more details.
710  	 */
711  	static
712  	void hcGradCut(
713  	   SCIP_Real             lbx,                /**< lower bound of x */
714  	   SCIP_Real             ubx,                /**< upper bound of x */
715  	   SCIP_Real             solx,               /**< solution value of x */
716  	   SCIP_Real             soly,               /**< solution value of y */
717  	   SCIP_Real*            coefx,              /**< pointer to store the coefficient of x */
718  	   SCIP_Real*            coefy,              /**< pointer to store the coefficient of y */
719  	   SCIP_Real*            constant            /**< pointer to store the constant */
720  	   )
721  	{
722  	   SCIP_Real tmp1;
723  	   SCIP_Real tmp2;
724  	
725  	   assert(lbx >= 0.0);
726  	   assert(lbx <= ubx);
727  	   assert(soly > 0.0);
728  	   assert(coefx != NULL);
729  	   assert(coefy != NULL);
730  	   assert(constant != NULL);
731  	
732  	   tmp1 = SQRT(lbx * ubx) + solx;
733  	   tmp2 = SQR(SQRT(lbx) + SQRT(ubx)) * soly; /*lint !e666*/
734  	   assert(tmp2 > 0.0);
735  	
736  	   *coefx = 2.0 * tmp1 / tmp2;
737  	   *coefy = -SQR(tmp1) / (tmp2 * soly);
738  	   *constant = 2.0 * SQRT(lbx * ubx) * tmp1 / tmp2;
739  	}
740  	
741  	/** computes an over- or underestimator at a given point for the bivariate case x/y &le;/&ge; z
742  	 *
743  	 *  There are the following cases for y > 0:
744  	 *
745  	 *    1. lbx < 0 < ubx:
746  	 *          Rewrite x / y = z as x = y * z and use McCormick to compute a valid inequality of the form
747  	 *          x = y * z &le; a * y +  b * z + c. Note that b > 0 because of y > 0. The inequality is then transformed
748  	 *          to x / b - a/b * y - c/b &le; z, which results in a valid underestimator for x / y over the set
749  	 *          {(x,y) | lbz &le; x / y &le; ubz}. Note that overestimating/underestimating the bilinear term with McCormick
750  	 *          results in an underestimator/overestimator for x / y.
751  	 *
752  	 *    2. lbx &ge; 0 or ubx &le; 0:
753  	 *       - overestimation:  use \f$z \leq \frac{1}{\text{lby}\cdot\text{uby}} \min(\text{uby}\cdot x - \text{lbx}\cdot y + \text{lbx}\cdot\text{lby}, \text{lby}\cdot x - \text{ubx}\cdot y + \text{ubx}\cdot\text{uby})\f$
754  	 *       - underestimation: use \f$z \geq x/y \geq \frac{1}{y} \frac{x + \sqrt{\text{lbx}\cdot\text{ubx}}}{\sqrt{\text{lbx} + \sqrt{\text{ubx}}}}\f$ and build gradient cut
755  	 *
756  	 *    If y < 0, swap and negate its bounds and compute the respective opposite estimator (and negate it).
757  	 *
758  	 *    If 0 is in the interval of y, nothing is possible.
759  	 */
760  	static
761  	SCIP_RETCODE estimateBivariate(
762  	   SCIP*                 scip,               /**< SCIP data structure */
763  	   SCIP_Real             lbx,                /**< lower bound of x */
764  	   SCIP_Real             ubx,                /**< upper bound of x */
765  	   SCIP_Real             lby,                /**< lower bound of y */
766  	   SCIP_Real             uby,                /**< upper bound of y */
767  	   SCIP_Real             lbz,                /**< lower bound of z */
768  	   SCIP_Real             ubz,                /**< lower bound of z */
769  	   SCIP_Real             solx,               /**< reference point for x */
770  	   SCIP_Real             soly,               /**< reference point for y */
771  	   SCIP_Real             solz,               /**< reference point for z */
772  	   SCIP_Bool             overestimate,       /**< whether the expression should be overestimated */
773  	   SCIP_Real*            coefx,              /**< pointer to store the x coefficient */
774  	   SCIP_Real*            coefy,              /**< pointer to store the y coefficient */
775  	   SCIP_Real*            constant,           /**< pointer to store the constant */
776  	   SCIP_Bool*            branchingusefulx,   /**< pointer to store whether branching on x would improve the estimator */
777  	   SCIP_Bool*            branchingusefuly,   /**< pointer to store whether branching on y would improve the estimator */
778  	   SCIP_Bool*            success             /**< buffer to store whether computing the estimator was successful */
779  	   )
780  	{
781  	   SCIP_Bool negatedx = FALSE;
782  	   SCIP_Bool negatedy = FALSE;
783  	
784  	   assert(lbx <= solx && solx <= ubx);
785  	   assert(lby <= soly && soly <= uby);
786  	   assert(lbz <= solz && solz <= ubz);
787  	   assert(coefx != NULL);
788  	   assert(coefy != NULL);
789  	   assert(constant != NULL);
790  	   assert(branchingusefulx != NULL);
791  	   assert(branchingusefuly != NULL);
792  	   assert(success != NULL);
793  	
794  	   *branchingusefulx = TRUE;
795  	   *branchingusefuly = TRUE;
796  	   *success = TRUE;
797  	   *coefx = 0.0;
798  	   *coefy = 0.0;
799  	   *constant = 0.0;
800  	
801  	   /* if 0 is in [lby,uby], then it is not possible to compute an estimator */
802  	   if( SCIPisLE(scip, lby, 0.0) && SCIPisGE(scip, uby, 0.0) )
803  	   {
804  	      *success = FALSE;
805  	      return SCIP_OKAY;
806  	   }
807  	
808  	   /* negate bounds of y if it is not positive */
809  	   if( uby < 0.0 )
810  	   {
811  	      SCIP_Real tmp = uby;
812  	
813  	      uby = -lby;
814  	      lby = -tmp;
815  	      soly = -soly;
816  	      negatedy = TRUE;
817  	      overestimate = !overestimate;
818  	   }
819  	
820  	   /* case 1: 0 is in the interior of [lbx,ubx] */
821  	   if( lbx < 0.0 && 0.0 < ubx )
822  	   {
823  	      SCIP_Real mccoefy = 0.0;
824  	      SCIP_Real mccoefaux = 0.0;
825  	      SCIP_Real mcconst = 0.0;
826  	
827  	      /* as explained in the description of this method, overestimating/underestimating the bilinear term results in an
828  	       * underestimator/overestimator for x / y
829  	       */
830  	      SCIPaddBilinMcCormick(scip, 1.0, lbz, ubz, solz, lby, uby, soly, !overestimate, &mccoefaux, &mccoefy, &mcconst,
831  	         success);
832  	      assert(mccoefaux >= 0.0);
833  	
834  	      if( !(*success) )
835  	         return SCIP_OKAY;
836  	
837  	      /* resulting estimator is x/b - a/b * y - c/b, where a*y +  b*z + c is the estimator for y*z */
838  	      *coefx = 1.0 / mccoefaux;
839  	      *coefy = -mccoefy / mccoefaux;
840  	      *constant = -mcconst / mccoefaux;
841  	   }
842  	   /* case 2: 0 is not in the interior of [lbx,ubx] */
843  	   else
844  	   {
845  	      /* negate bounds of x if it is negative */
846  	      if( ubx <= 0.0 )
847  	      {
848  	         SCIP_Real tmp = ubx;
849  	
850  	         ubx = -lbx;
851  	         lbx = -tmp;
852  	         solx = -solx;
853  	         negatedx = TRUE;
854  	         overestimate = !overestimate;
855  	      }
856  	
857  	      /* case 2a */
858  	      if( overestimate )
859  	      {
860  	         /* check where the minimum is attained */
861  	         if( uby * solx - lbx * soly + lbx * lby <= lby * solx - ubx * soly + ubx * uby )
862  	         {
863  	            *coefx = 1.0 / lby;
864  	            *coefy = -lbx / (lby * uby);
865  	            *constant = lbx / uby;
866  	         }
867  	         else
868  	         {
869  	            *coefx = 1.0 / uby;
870  	            *coefy = -ubx / (lby * uby);
871  	            *constant = ubx / lby;
872  	         }
873  	      }
874  	      /* case 2b */
875  	      else
876  	      {
877  	         /* compute gradient cut for h^c(x,y) at (solx,soly) */
878  	         hcGradCut(lbx, ubx, solx, soly, coefx, coefy, constant);
879  	
880  	         /* estimator is independent of the bounds of y */
881  	         *branchingusefuly = FALSE;
882  	      }
883  	   }
884  	
885  	   /* reverse negations of x and y in the resulting estimator */
886  	   if( negatedx )
887  	      *coefx = -(*coefx);
888  	   if( negatedy )
889  	      *coefy = -(*coefy);
890  	
891  	   /* if exactly one variable has been negated, then we have computed an underestimate/overestimate for the negated
892  	    * expression, which results in an overestimate/underestimate for the original expression
893  	    */
894  	   if( negatedx != negatedy )
895  	   {
896  	      *coefx = -(*coefx);
897  	      *coefy = -(*coefy);
898  	      *constant = -(*constant);
899  	   }
900  	
901  	   /* avoid huge values in the estimator */
902  	   if( SCIPisHugeValue(scip, REALABS(*coefx)) || SCIPisHugeValue(scip, REALABS(*coefy))
903  	      || SCIPisHugeValue(scip, REALABS(*constant)) )
904  	   {
905  	      *success = FALSE;
906  	      return SCIP_OKAY;
907  	   }
908  	
909  	   return SCIP_OKAY;
910  	}
911  	
912  	/** construct an estimator for a quotient expression of the form (ax + b) / (cy + d) + e
913  	 *
914  	 *  The resulting estimator is stored in a rowprep.
915  	 *
916  	 *  The method first computes an estimator for x' / y' with x := ax + b and y := cy + d
917  	 *  and then transforms this estimator to one for the quotient (ax + b) / (cy + d) + e.
918  	 */
919  	static
920  	SCIP_RETCODE estimateBivariateQuotient(
921  	   SCIP*                 scip,               /**< SCIP data structure */
922  	   SCIP_EXPR*            xexpr,              /**< numerator expression */
923  	   SCIP_EXPR*            yexpr,              /**< denominator expression */
924  	   SCIP_VAR*             auxvar,             /**< auxiliary variable */
925  	   SCIP_SOL*             sol,                /**< solution point (or NULL for the LP solution) */
926  	   SCIP_Real             a,                  /**< coefficient of numerator */
927  	   SCIP_Real             b,                  /**< constant of numerator */
928  	   SCIP_Real             c,                  /**< coefficient of denominator */
929  	   SCIP_Real             d,                  /**< constant of denominator */
930  	   SCIP_Real             e,                  /**< constant term */
931  	   SCIP_Bool             overestimate,       /**< whether the expression should be overestimated */
932  	   SCIP_ROWPREP*         rowprep,            /**< a rowprep where to store the estimator */
933  	   SCIP_Bool*            branchingusefulx,   /**< pointer to store whether branching on x would improve the estimator */
934  	   SCIP_Bool*            branchingusefuly,   /**< pointer to store whether branching on y would improve the estimator */
935  	   SCIP_Bool*            success             /**< buffer to store whether separation was successful */
936  	   )
937  	{
938  	   SCIP_VAR* vars[2];
939  	   SCIP_Real coefs[2] = {0.0, 0.0};
940  	   SCIP_Real constant = 0.0;
941  	   SCIP_Real solx;
942  	   SCIP_Real soly;
943  	   SCIP_Real solz;
944  	   SCIP_Real lbx;
945  	   SCIP_Real ubx;
946  	   SCIP_Real lby;
947  	   SCIP_Real uby;
948  	   SCIP_Real lbz;
949  	   SCIP_Real ubz;
950  	   SCIP_INTERVAL bnd;
951  	
952  	   assert(xexpr != NULL);
953  	   assert(yexpr != NULL);
954  	   assert(xexpr != yexpr);
955  	   assert(auxvar != NULL);
956  	   assert(rowprep != NULL);
957  	   assert(branchingusefulx != NULL);
958  	   assert(branchingusefuly != NULL);
959  	   assert(success != NULL);
960  	
961  	   vars[0] = SCIPgetExprAuxVarNonlinear(xexpr);
962  	   vars[1] = SCIPgetExprAuxVarNonlinear(yexpr);
963  	
964  	   /* get bounds for x, y, and z */
965  	   SCIPintervalSetBounds(&bnd,
966  	      -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(vars[0])),
967  	       infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY,  SCIPvarGetUbGlobal(vars[0])));
968  	   SCIP_CALL( SCIPevalExprActivity(scip, xexpr) );
969  	   SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(xexpr));
970  	   lbx = bnd.inf;
971  	   ubx = bnd.sup;
972  	
973  	   SCIPintervalSetBounds(&bnd,
974  	      -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(vars[1])),
975  	       infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY,  SCIPvarGetUbGlobal(vars[1])));
976  	   SCIP_CALL( SCIPevalExprActivity(scip, yexpr) );
977  	   SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(yexpr));
978  	   lby = bnd.inf;
979  	   uby = bnd.sup;
980  	
981  	   lbz = SCIPvarGetLbLocal(auxvar);
982  	   ubz = SCIPvarGetUbLocal(auxvar);
983  	
984  	   /* check whether one of the variables has been fixed or has empty domain */
985  	   if( SCIPisEQ(scip, lbx, ubx) || SCIPisEQ(scip, lby, uby) || ubx < lbx || uby < lby )
986  	   {
987  	      *success = FALSE;
988  	      return SCIP_OKAY;
989  	   }
990  	
991  	   /* get and adjust solution values */
992  	   solx = SCIPgetSolVal(scip, sol, vars[0]);
993  	   soly = SCIPgetSolVal(scip, sol, vars[1]);
994  	   solz = SCIPgetSolVal(scip, sol, auxvar);
995  	   solx = MIN(MAX(solx, lbx), ubx);
996  	   soly = MIN(MAX(soly, lby), uby);
997  	   solz = MIN(MAX(solz, lbz), ubz);
998  	
999  	   /* compute an estimator */
1000 	   SCIP_CALL( estimateBivariate(scip,
1001 	      MIN(a * lbx, a * ubx) + b, MAX(a * lbx, a * ubx) + b, /* bounds of x' */
1002 	      MIN(c * lby, c * uby) + d, MAX(c * lby, c * uby) + d, /* bounds of y' */
1003 	      lbz, ubz, a * solx + b, c * soly + d, solz, overestimate, &coefs[0], &coefs[1], &constant,
1004 	      branchingusefulx, branchingusefuly, success) );
1005 	
1006 	   /* add estimator to rowprep, if successful */
1007 	   if( *success )
1008 	   {
1009 	      /* transform estimator Ax' + By'+ C = A(ax + b) + B (cy + d) + C = (Aa) x + (Bc) y + (C + Ab + Bd);
1010 	       * add the constant e separately
1011 	       */
1012 	      constant += coefs[0] * b + coefs[1] * d + e;
1013 	      coefs[0] *= a;
1014 	      coefs[1] *= c;
1015 	
1016 	      /* prepare rowprep */
1017 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_%s_%lld", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]),
1018 	         SCIPgetNLPs(scip));
1019 	      SCIP_CALL( createRowprep(scip, rowprep, vars, coefs, constant, 2) );
1020 	   }
1021 	
1022 	   return SCIP_OKAY;
1023 	}
1024 	
1025 	/*
1026 	 * Callback methods of nonlinear handler
1027 	 */
1028 	
1029 	/** nonlinear handler copy callback */
1030 	static
1031 	SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrQuotient)
1032 	{ /*lint --e{715}*/
1033 	   assert(targetscip != NULL);
1034 	   assert(sourcenlhdlr != NULL);
1035 	   assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
1036 	
1037 	   SCIP_CALL( SCIPincludeNlhdlrQuotient(targetscip) );
1038 	
1039 	   return SCIP_OKAY;
1040 	}
1041 	
1042 	
1043 	/** callback to free expression specific data */
1044 	static
1045 	SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataQuotient)
1046 	{  /*lint --e{715}*/
1047 	   assert(nlhdlrexprdata != NULL);
1048 	   assert(*nlhdlrexprdata != NULL);
1049 	
1050 	   /* free expression data of nonlinear handler */
1051 	   SCIP_CALL( exprdataFree(scip, nlhdlrexprdata) );
1052 	
1053 	   return SCIP_OKAY;
1054 	}
1055 	
1056 	
1057 	/** callback to detect structure in expression tree */
1058 	static
1059 	SCIP_DECL_NLHDLRDETECT(nlhdlrDetectQuotient)
1060 	{ /*lint --e{715}*/
1061 	   SCIP_Bool success;
1062 	
1063 	   assert(nlhdlrexprdata != NULL);
1064 	
1065 	   /* call detection routine */
(1) Event cond_false: Condition "(_restat_ = detectExpr(scip, expr, nlhdlrexprdata, &success)) != SCIP_OKAY", taking false branch.
(2) Event if_end: End of if statement.
1066 	   SCIP_CALL( detectExpr(scip, expr, nlhdlrexprdata, &success) );
1067 	
(3) Event cond_true: Condition "success", taking true branch.
1068 	   if( success )
1069 	   {
(4) Event cond_true: Condition "SCIPgetExprNAuxvarUsesNonlinear(expr) > 0", taking true branch.
1070 	      if( SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 )
1071 	         *participating = SCIP_NLHDLR_METHOD_SEPABOTH;
1072 	
(5) Event deref_parm: Directly dereferencing parameter "*nlhdlrexprdata".
1073 	      if( (*nlhdlrexprdata)->numexpr == (*nlhdlrexprdata)->denomexpr )
1074 	      {
1075 	         /* if univariate, then we also do inteval and reverseprop */
1076 	         *participating |= SCIP_NLHDLR_METHOD_ACTIVITY;
1077 	
1078 	         /* if univariate, then all our methods are enforcing */
1079 	         *enforcing |= *participating;
1080 	      }
1081 	   }
1082 	
1083 	   return SCIP_OKAY;
1084 	}
1085 	
1086 	
1087 	/** auxiliary evaluation callback of nonlinear handler */
1088 	static
1089 	SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxQuotient)
1090 	{ /*lint --e{715}*/
1091 	   SCIP_VAR* auxvarx;
1092 	   SCIP_VAR* auxvary;
1093 	   SCIP_Real solvalx;
1094 	   SCIP_Real solvaly;
1095 	   SCIP_Real nomval;
1096 	   SCIP_Real denomval;
1097 	
1098 	   assert(expr != NULL);
1099 	   assert(auxvalue != NULL);
1100 	
1101 	   /**! [SnippetNlhdlrEvalauxQuotient] */
1102 	   /* get auxiliary variables */
1103 	   auxvarx = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->numexpr);
1104 	   auxvary = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->denomexpr);
1105 	   assert(auxvarx != NULL);
1106 	   assert(auxvary != NULL);
1107 	
1108 	   /* get solution values of the auxiliary variables */
1109 	   solvalx = SCIPgetSolVal(scip, sol, auxvarx);
1110 	   solvaly = SCIPgetSolVal(scip, sol, auxvary);
1111 	
1112 	   /* evaluate expression w.r.t. the values of the auxiliary variables */
1113 	   nomval = nlhdlrexprdata->numcoef *  solvalx + nlhdlrexprdata->numconst;
1114 	   denomval = nlhdlrexprdata->denomcoef *  solvaly + nlhdlrexprdata->denomconst;
1115 	
1116 	   /* return SCIP_INVALID if the denominator evaluates to zero */
1117 	   *auxvalue = (denomval != 0.0) ? nlhdlrexprdata->constant + nomval / denomval : SCIP_INVALID;
1118 	   /**! [SnippetNlhdlrEvalauxQuotient] */
1119 	
1120 	   return SCIP_OKAY;
1121 	}
1122 	
1123 	
1124 	/** nonlinear handler under/overestimation callback
1125 	 *
1126 	 * @todo which of the paramters did I not use, but have to be taken into consideration?
1127 	*/
1128 	static
1129 	SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateQuotient)
1130 	{ /*lint --e{715}*/
1131 	   SCIP_Bool branchingusefulx = FALSE;
1132 	   SCIP_Bool branchingusefuly = FALSE;
1133 	   SCIP_ROWPREP* rowprep;
1134 	
1135 	   assert(nlhdlr != NULL);
1136 	   assert(expr != NULL);
1137 	   assert(nlhdlrexprdata != NULL);
1138 	   assert(rowpreps != NULL);
1139 	
1140 	   /** ![SnippetNlhdlrEstimateQuotient] */
1141 	   *addedbranchscores = FALSE;
1142 	   *success = FALSE;
1143 	
1144 	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
1145 	
1146 	   if( nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr )
1147 	   {
1148 	      /* univariate case */
1149 	      SCIP_CALL( estimateUnivariateQuotient(scip, sol, nlhdlrexprdata->numexpr, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst,
1150 	         nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant, overestimate, rowprep,
1151 	         &branchingusefulx, success) );
1152 	   }
1153 	   else
1154 	   {
1155 	      /* bivariate case */
1156 	      SCIP_CALL( estimateBivariateQuotient(scip, nlhdlrexprdata->numexpr, nlhdlrexprdata->denomexpr, SCIPgetExprAuxVarNonlinear(expr), sol,
1157 	         nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst,
1158 	         nlhdlrexprdata->constant, overestimate, rowprep,
1159 	         &branchingusefulx, &branchingusefuly, success) );
1160 	   }
1161 	
1162 	   if( *success )
1163 	   {
1164 	      SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
1165 	   }
1166 	   else
1167 	   {
1168 	      SCIPfreeRowprep(scip, &rowprep);
1169 	   }
1170 	
1171 	   /* add branching scores if requested */
1172 	   if( addbranchscores )
1173 	   {
1174 	      SCIP_EXPR* exprs[2];
1175 	      SCIP_Real violation;
1176 	      int nexprs = 0;
1177 	
1178 	      if( branchingusefulx )
1179 	         exprs[nexprs++] = nlhdlrexprdata->numexpr;
1180 	      if( branchingusefuly )
1181 	         exprs[nexprs++] = nlhdlrexprdata->denomexpr;
1182 	
1183 	      /* compute violation w.r.t. the auxiliary variable(s) */
1184 	#ifndef BRSCORE_ABSVIOL
1185 	      SCIP_CALL( SCIPgetExprRelAuxViolationNonlinear(scip, expr, auxvalue, sol, &violation, NULL, NULL) );
1186 	#else
1187 	      SCIP_CALL( SCIPgetExprAbsAuxViolationNonlinear(scip, expr, auxvalue, sol, &violation, NULL, NULL) );
1188 	#endif
1189 	      assert(violation > 0.0);  /* there should be a violation if we were called to enforce */
1190 	
1191 	      SCIP_CALL( SCIPaddExprsViolScoreNonlinear(scip, exprs, nexprs, violation, sol, addedbranchscores) );
1192 	   }
1193 	   /** ![SnippetNlhdlrEstimateQuotient] */
1194 	
1195 	   return SCIP_OKAY;
1196 	}
1197 	
1198 	/** nonlinear handler solution linearization callback */
1199 	static
1200 	SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeQuotient)
1201 	{ /*lint --e{715}*/
1202 	   SCIP_VAR* x;
1203 	   SCIP_Real lbx;
1204 	   SCIP_Real ubx;
1205 	   SCIP_Real solx;
1206 	   int c;
1207 	
1208 	   assert(nlhdlr != NULL);
1209 	   assert(expr != NULL);
1210 	   assert(nlhdlrexprdata != NULL);
1211 	
1212 	   /* in the bivariate case, we do not get globally valid estimators */
1213 	   if( nlhdlrexprdata->numexpr != nlhdlrexprdata->denomexpr )
1214 	      return SCIP_OKAY;
1215 	
1216 	   x = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->numexpr);
1217 	   lbx = SCIPvarGetLbGlobal(x);
1218 	   ubx = SCIPvarGetUbGlobal(x);
1219 	   solx = SCIPgetSolVal(scip, sol, x);
1220 	   solx = MAX(MIN(solx, ubx), lbx);
1221 	
1222 	   for( c = (overestimate ? 0 : 1); c < (underestimate ? 2 : 1); ++c )  /* c == 0: overestimate, c == 1: underestimate */
1223 	   {
1224 	      SCIP_ROWPREP* rowprep;
1225 	      SCIP_Bool success = FALSE;
1226 	      SCIP_Bool local = TRUE;
1227 	      SCIP_Bool branchinguseful;
1228 	      SCIP_Real coef;
1229 	      SCIP_Real constant;
1230 	
1231 	      /* compute estimator, will be secant or gradient */
1232 	      SCIP_CALL( estimateUnivariate(scip, lbx, ubx, lbx, ubx, solx,
1233 	         nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant,
1234 	         &coef, &constant, c == 0, &local, &branchinguseful, &success) );
1235 	
1236 	      /* skip if not successful, only locally valid, or just a secant (branchinguseful is TRUE) */
1237 	      if( !success || local || branchinguseful )
1238 	         continue;
1239 	
1240 	      SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, c == 0 ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, FALSE) );
1241 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, x, coef) );
1242 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) );
1243 	      SCIProwprepAddConstant(rowprep, constant);
1244 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_sol%d", SCIPvarGetName(x), SCIPsolGetIndex(sol));
1245 	
1246 	      SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) );
1247 	
1248 	      /* if cleanup succeeded and rowprep is still global, add to cutpool */
1249 	      if( success && !SCIProwprepIsLocal(rowprep) )
1250 	      {
1251 	         SCIP_ROW* row;
1252 	
1253 	         SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );
1254 	         SCIP_CALL( SCIPaddPoolCut(scip, row) );
1255 	         SCIP_CALL( SCIPreleaseRow(scip, &row) );
1256 	      }
1257 	
1258 	      SCIPfreeRowprep(scip, &rowprep);
1259 	   }
1260 	
1261 	   return SCIP_OKAY;
1262 	}
1263 	
1264 	/** nonlinear handler interval evaluation callback */
1265 	static
1266 	SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalQuotient)
1267 	{ /*lint --e{715}*/
1268 	   SCIP_INTERVAL bnds;
1269 	
1270 	   assert(nlhdlrexprdata != NULL);
1271 	   assert(nlhdlrexprdata->numexpr != NULL);
1272 	   assert(nlhdlrexprdata->denomexpr != NULL);
1273 	
1274 	   /* it is not possible to compute tighter intervals if both expressions are different
1275 	    * we should not be called in this case, as we haven't said we would participate in this activity in detect
1276 	    */
1277 	   assert(nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr);
1278 	
1279 	   /**! [SnippetNlhdlrIntevalQuotient] */
1280 	   /* get activity of the numerator (= denominator) expression */
1281 	   bnds = SCIPexprGetActivity(nlhdlrexprdata->numexpr);
1282 	
1283 	   /* call interval evaluation for the univariate quotient expression */
1284 	   *interval = intEvalQuotient(scip, bnds, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst,
1285 	      nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant);
1286 	   /**! [SnippetNlhdlrIntevalQuotient] */
1287 	
1288 	   return SCIP_OKAY;
1289 	}
1290 	
1291 	
1292 	/** nonlinear handler callback for reverse propagation */
1293 	static
1294 	SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropQuotient)
1295 	{ /*lint --e{715}*/
1296 	   SCIP_INTERVAL result;
1297 	
1298 	   assert(nlhdlrexprdata != NULL);
1299 	   assert(nlhdlrexprdata->numexpr != NULL);
1300 	   assert(nlhdlrexprdata->denomexpr != NULL);
1301 	
1302 	   /* it is not possible to compute tighter intervals if both expressions are different
1303 	    * we should not be called in this case, as we haven't said we would participate in this activity in detect
1304 	    */
1305 	   assert(nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr);
1306 	
1307 	   SCIPdebugMsg(scip, "call reverse propagation for expression (%g %p + %g) / (%g %p + %g) + %g bounds [%g,%g]\n",
1308 	      nlhdlrexprdata->numcoef, (void*)nlhdlrexprdata->numexpr, nlhdlrexprdata->numconst,
1309 	      nlhdlrexprdata->denomcoef, (void*)nlhdlrexprdata->denomexpr, nlhdlrexprdata->denomconst,
1310 	      nlhdlrexprdata->constant, bounds.inf, bounds.sup);
1311 	
1312 	   /* call reverse propagation */
1313 	   /**! [SnippetNlhdlrReversepropQuotient] */
1314 	   result = reversepropQuotient(bounds, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst,
1315 	      nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant);
1316 	
1317 	   SCIPdebugMsg(scip, "try to tighten bounds of %p: [%g,%g] -> [%g,%g]\n",
1318 	      (void*)nlhdlrexprdata->numexpr, SCIPgetExprBoundsNonlinear(scip, nlhdlrexprdata->numexpr).inf,
1319 	      SCIPgetExprBoundsNonlinear(scip, nlhdlrexprdata->numexpr).sup, result.inf, result.sup);
1320 	
1321 	   /* tighten bounds of the expression */
1322 	   SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, nlhdlrexprdata->numexpr, result, infeasible, nreductions) );
1323 	   /**! [SnippetNlhdlrReversepropQuotient] */
1324 	
1325 	   return SCIP_OKAY;
1326 	}
1327 	
1328 	
1329 	/*
1330 	 * nonlinear handler specific interface methods
1331 	 */
1332 	
1333 	/** includes quotient nonlinear handler in nonlinear constraint handler */
1334 	SCIP_RETCODE SCIPincludeNlhdlrQuotient(
1335 	   SCIP*                 scip                /**< SCIP data structure */
1336 	   )
1337 	{
1338 	   SCIP_NLHDLRDATA* nlhdlrdata;
1339 	   SCIP_NLHDLR* nlhdlr;
1340 	
1341 	   assert(scip != NULL);
1342 	
1343 	   /* create nonlinear handler data */
1344 	   nlhdlrdata = NULL;
1345 	
1346 	   SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY,
1347 	      NLHDLR_ENFOPRIORITY, nlhdlrDetectQuotient, nlhdlrEvalauxQuotient, nlhdlrdata) );
1348 	   assert(nlhdlr != NULL);
1349 	
1350 	   SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrQuotient);
1351 	   SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataQuotient);
1352 	   SCIPnlhdlrSetSepa(nlhdlr, NULL, NULL, nlhdlrEstimateQuotient, NULL);
1353 	   SCIPnlhdlrSetSollinearize(nlhdlr, nlhdlrSollinearizeQuotient);
1354 	   SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalQuotient, nlhdlrReversepropQuotient);
1355 	
1356 	   return SCIP_OKAY;
1357 	}
1358