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_convex.c
26   	 * @ingroup DEFPLUGINS_NLHDLR
27   	 * @brief  nonlinear handlers for convex and concave expressions
28   	 * @author Benjamin Mueller
29   	 * @author Stefan Vigerske
30   	 */
31   	
32   	#include <string.h>
33   	
34   	#include "scip/nlhdlr_convex.h"
35   	#include "scip/pub_nlhdlr.h"
36   	#include "scip/scip_expr.h"
37   	#include "scip/cons_nonlinear.h"
38   	#include "scip/expr_var.h"
39   	#include "scip/expr_abs.h"
40   	#include "scip/pub_misc_rowprep.h"
41   	#include "scip/dbldblarith.h"
42   	
43   	/* fundamental nonlinear handler properties */
44   	#define CONVEX_NLHDLR_NAME             "convex"
45   	#define CONVEX_NLHDLR_DESC             "handler that identifies and estimates convex expressions"
46   	#define CONVEX_NLHDLR_DETECTPRIORITY   50
47   	#define CONVEX_NLHDLR_ENFOPRIORITY     50
48   	
49   	#define CONCAVE_NLHDLR_NAME            "concave"
50   	#define CONCAVE_NLHDLR_DESC            "handler that identifies and estimates concave expressions"
51   	#define CONCAVE_NLHDLR_DETECTPRIORITY  40
52   	#define CONCAVE_NLHDLR_ENFOPRIORITY    40
53   	
54   	#define DEFAULT_DETECTSUM             FALSE
55   	#define DEFAULT_EXTENDEDFORM          TRUE
56   	#define DEFAULT_CVXQUADRATIC_CONVEX   TRUE
57   	#define DEFAULT_CVXQUADRATIC_CONCAVE  FALSE
58   	#define DEFAULT_CVXSIGNOMIAL          TRUE
59   	#define DEFAULT_CVXPRODCOMP           TRUE
60   	#define DEFAULT_HANDLETRIVIAL         FALSE
61   	#define DEFAULT_MAXPERTURB            0.01
62   	
63   	#define INITLPMAXVARVAL          1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */
64   	#define RANDNUMINITSEED          220802 /**< initial seed for random number generator for point perturbation */
65   	
66   	/*lint -e440*/
67   	/*lint -e441*/
68   	/*lint -e666*/
69   	/*lint -e777*/
70   	
71   	/*
72   	 * Data structures
73   	 */
74   	
75   	/** nonlinear handler expression data */
76   	struct SCIP_NlhdlrExprData
77   	{
78   	   SCIP_EXPR*            nlexpr;             /**< expression (copy) for which this nlhdlr estimates */
79   	   SCIP_HASHMAP*         nlexpr2origexpr;    /**< mapping of our copied expression to original expression */
80   	
81   	   int                   nleafs;             /**< number of distinct leafs of nlexpr, i.e., number of distinct (auxiliary) variables handled */
82   	   SCIP_EXPR**           leafexprs;          /**< distinct leaf expressions (excluding value-expressions), thus variables */
83   	};
84   	
85   	/** nonlinear handler data */
86   	struct SCIP_NlhdlrData
87   	{
88   	   SCIP_Bool             isnlhdlrconvex;     /**< whether this data is used for the convex nlhdlr (TRUE) or the concave one (FALSE) */
89   	   SCIP_SOL*             evalsol;            /**< solution used for evaluating expression in a different point,
90   	                                                  e.g., for facet computation of vertex-polyhedral function */
91   	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator used to perturb reference point in estimateGradient() */
92   	
93   	   /* parameters */
94   	   SCIP_Bool             detectsum;          /**< whether to run detection when the root of an expression is a non-quadratic sum */
95   	   SCIP_Bool             extendedform;       /**< whether to create extended formulations instead of looking for maximal possible subexpression */
96   	   SCIP_Real             maxperturb;         /**< maximal perturbation of non-differentiable reference point */
97   	
98   	   /* advanced parameters (maybe remove some day) */
99   	   SCIP_Bool             cvxquadratic;       /**< whether to use convexity check on quadratics */
100  	   SCIP_Bool             cvxsignomial;       /**< whether to use convexity check on signomials */
101  	   SCIP_Bool             cvxprodcomp;        /**< whether to use convexity check on product composition f(h)*h */
102  	   SCIP_Bool             handletrivial;      /**< whether to handle trivial expressions, i.e., those where all children are variables */
103  	};
104  	
105  	/** data struct to be be passed on to vertexpoly-evalfunction (see SCIPcomputeFacetVertexPolyhedralNonlinear) */
106  	typedef struct
107  	{
108  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata;
109  	   SCIP_SOL*             evalsol;
110  	   SCIP*                 scip;
111  	} VERTEXPOLYFUN_EVALDATA;
112  	
113  	/** stack used in constructExpr to store expressions that need to be investigated ("to do list") */
114  	typedef struct
115  	{
116  	   SCIP_EXPR**           stack;              /**< stack elements */
117  	   int                   stacksize;          /**< allocated space (in number of pointers) */
118  	   int                   stackpos;           /**< position of top element of stack */
119  	} EXPRSTACK;
120  	
121  	#define DECL_CURVCHECK(x) SCIP_RETCODE x( \
122  	   SCIP*                 scip,               /**< SCIP data structure */ \
123  	   SCIP_EXPR*            nlexpr,             /**< nlhdlr-expr to check */ \
124  	   SCIP_Bool             isrootexpr,         /**< whether nlexpr is the root from where detection has been started */ \
125  	   EXPRSTACK*            stack,              /**< stack where to add generated leafs */ \
126  	   SCIP_HASHMAP*         nlexpr2origexpr,    /**< mapping from our expression copy to original expression */ \
127  	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< data of nlhdlr */ \
128  	   SCIP_HASHMAP*         assumevarfixed,     /**< hashmap containing variables that should be assumed to be fixed, or NULL */ \
129  	   SCIP_Bool*            success             /**< whether we found something */ \
130  	   )
131  	
132  	/*
133  	 * static methods
134  	 */
135  	
136  	/** create nlhdlr-expression
137  	 *
138  	 * does not create children, i.e., assumes that this will be a leaf
139  	 */
140  	static
141  	SCIP_RETCODE nlhdlrExprCreate(
142  	   SCIP*                 scip,               /**< SCIP data structure */
143  	   SCIP_HASHMAP*         nlexpr2origexpr,    /**< mapping from copied to original expression */
144  	   SCIP_EXPR**           nlhdlrexpr,         /**< buffer to store created expr */
145  	   SCIP_EXPR*            origexpr,           /**< original expression to be copied */
146  	   SCIP_EXPRCURV         curv                /**< curvature to achieve */
147  	   )
148  	{
149  	   assert(scip != NULL);
150  	   assert(nlexpr2origexpr != NULL);
151  	   assert(nlhdlrexpr != NULL);
152  	   assert(origexpr != NULL);
153  	
154  	   if( SCIPexprGetNChildren(origexpr) == 0 )
155  	   {
156  	      /* for leaves, do not copy */
157  	      *nlhdlrexpr = origexpr;
158  	      SCIPcaptureExpr(*nlhdlrexpr);
159  	      if( !SCIPhashmapExists(nlexpr2origexpr, (void*)*nlhdlrexpr) )
160  	      {
161  	         SCIP_CALL( SCIPhashmapInsert(nlexpr2origexpr, (void*)*nlhdlrexpr, (void*)origexpr) );
162  	      }
163  	      return SCIP_OKAY;
164  	   }
165  	
166  	   /* create copy of expression, but without children */
167  	   SCIP_CALL( SCIPduplicateExprShallow(scip, origexpr, nlhdlrexpr, NULL, NULL) );
168  	   assert(*nlhdlrexpr != NULL);  /* copies within the same SCIP must always work */
169  	
170  	   /* store the curvature we want to get in the curvature flag of the copied expression
171  	    * it's a bit of a misuse, but once we are done with everything, this is actually correct
172  	    */
173  	   SCIPexprSetCurvature(*nlhdlrexpr, curv);
174  	
175  	   /* remember which the original expression was */
176  	   SCIP_CALL( SCIPhashmapInsert(nlexpr2origexpr, (void*)*nlhdlrexpr, (void*)origexpr) );
177  	
178  	   return SCIP_OKAY;
179  	}
180  	
181  	/** expand nlhdlr-expression by adding children according to original expression */
182  	static
183  	SCIP_RETCODE nlhdlrExprGrowChildren(
184  	   SCIP*                 scip,               /**< SCIP data structure */
185  	   SCIP_HASHMAP*         nlexpr2origexpr,    /**< mapping from copied to original expression */
186  	   SCIP_EXPR*            nlhdlrexpr,         /**< expression for which to create children */
187  	   SCIP_EXPRCURV*        childrencurv        /**< curvature required for children, or NULL if to set to UNKNOWN */
188  	   )
189  	{
190  	   SCIP_EXPR* origexpr;
191  	   SCIP_EXPR* child;
192  	   int nchildren;
193  	   int i;
194  	
195  	   assert(scip != NULL);
196  	   assert(nlhdlrexpr != NULL);
197  	   assert(SCIPexprGetNChildren(nlhdlrexpr) == 0);
198  	
199  	   origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlhdlrexpr);
200  	
201  	   nchildren = SCIPexprGetNChildren(origexpr);
202  	   if( nchildren == 0 )
203  	      return SCIP_OKAY;
204  	
205  	   for( i = 0; i < nchildren; ++i )
206  	   {
207  	      SCIP_CALL( nlhdlrExprCreate(scip, nlexpr2origexpr, &child, SCIPexprGetChildren(origexpr)[i],
208  	         childrencurv != NULL ? childrencurv[i] : SCIP_EXPRCURV_UNKNOWN) );
209  	      SCIP_CALL( SCIPappendExprChild(scip, nlhdlrexpr, child) );
210  	      /* append captures child, so we can release the capture from nlhdlrExprCreate */
211  	      SCIP_CALL( SCIPreleaseExpr(scip, &child) );
212  	   }
213  	
214  	   assert(SCIPexprGetNChildren(nlhdlrexpr) == SCIPexprGetNChildren(origexpr));
215  	
216  	   return SCIP_OKAY;
217  	}
218  	
219  	/** evaluate expression at solution w.r.t. auxiliary variables */
220  	static
221  	SCIP_DECL_VERTEXPOLYFUN(nlhdlrExprEvalConcave)
222  	{
223  	   VERTEXPOLYFUN_EVALDATA* evaldata = (VERTEXPOLYFUN_EVALDATA*)funcdata;
224  	   int i;
225  	
226  	   assert(args != NULL);
227  	   assert(nargs == evaldata->nlhdlrexprdata->nleafs);
228  	   assert(evaldata != NULL);
229  	
230  	#ifdef SCIP_MORE_DEBUG
231  	   SCIPdebugMsg(evaldata->scip, "eval vertexpolyfun at\n");
232  	#endif
233  	   for( i = 0; i < nargs; ++i )
234  	   {
235  	#ifdef SCIP_MORE_DEBUG
236  	      SCIPdebugMsg(evaldata->scip, "  <%s> = %g\n",
237  	            SCIPvarGetName(SCIPgetVarExprVar(evaldata->nlhdlrexprdata->leafexprs[i])), args[i]);
238  	#endif
239  	      SCIP_CALL_ABORT( SCIPsetSolVal(evaldata->scip, evaldata->evalsol,
240  	            SCIPgetVarExprVar(evaldata->nlhdlrexprdata->leafexprs[i]), args[i]) );
241  	   }
242  	
243  	   SCIP_CALL_ABORT( SCIPevalExpr(evaldata->scip, evaldata->nlhdlrexprdata->nlexpr, evaldata->evalsol, 0L) );
244  	
245  	   return SCIPexprGetEvalValue(evaldata->nlhdlrexprdata->nlexpr);
246  	}
247  	
248  	/** initialize expression stack */
249  	static
250  	SCIP_RETCODE exprstackInit(
251  	   SCIP*                 scip,               /**< SCIP data structure */
252  	   EXPRSTACK*            exprstack,          /**< stack to initialize */
253  	   int                   initsize            /**< initial size */
254  	   )
255  	{
256  	   assert(scip != NULL);
257  	   assert(exprstack != NULL);
258  	   assert(initsize > 0);
259  	
260  	   SCIP_CALL( SCIPallocBufferArray(scip, &exprstack->stack, initsize) );
261  	   exprstack->stacksize = initsize;
262  	   exprstack->stackpos = -1;
263  	
264  	   return SCIP_OKAY;
265  	}
266  	
267  	/** free expression stack */
268  	static
269  	void exprstackFree(
270  	   SCIP*                 scip,               /**< SCIP data structure */
271  	   EXPRSTACK*            exprstack           /**< free expression stack */
272  	   )
273  	{
274  	   assert(scip != NULL);
275  	   assert(exprstack != NULL);
276  	
277  	   SCIPfreeBufferArray(scip, &exprstack->stack);
278  	}
279  	
280  	/** add expressions to expression stack */
281  	static
282  	SCIP_RETCODE exprstackPush(
283  	   SCIP*                 scip,               /**< SCIP data structure */
284  	   EXPRSTACK*            exprstack,          /**< expression stack */
285  	   int                   nexprs,             /**< number of expressions to push */
286  	   SCIP_EXPR**           exprs               /**< expressions to push */
287  	   )
288  	{
289  	   assert(scip != NULL);
290  	   assert(exprstack != NULL);
291  	
292  	   if( nexprs == 0 )
293  	      return SCIP_OKAY;
294  	
295  	   assert(exprs != NULL);
296  	
297  	   if( exprstack->stackpos+1 + nexprs > exprstack->stacksize )  /*lint !e644*/
298  	   {
299  	      exprstack->stacksize = SCIPcalcMemGrowSize(scip, exprstack->stackpos+1 + nexprs);    /*lint !e644*/
300  	      SCIP_CALL( SCIPreallocBufferArray(scip, &exprstack->stack, exprstack->stacksize) );
301  	   }
302  	
303  	   memcpy(exprstack->stack + (exprstack->stackpos+1), exprs, nexprs * sizeof(SCIP_EXPR*));  /*lint !e679*/ /*lint !e737*/
304  	   exprstack->stackpos += nexprs;
305  	
306  	   return SCIP_OKAY;
307  	}
308  	
309  	/** gives expression from top of expression stack and removes it from stack */
310  	static
311  	SCIP_EXPR* exprstackPop(
312  	   EXPRSTACK*            exprstack           /**< expression stack */
313  	   )
314  	{
315  	   assert(exprstack != NULL);
316  	   assert(exprstack->stackpos >= 0);
317  	
318  	   return exprstack->stack[exprstack->stackpos--];
319  	}
320  	
321  	/** indicate whether expression stack is empty */
322  	static
323  	SCIP_Bool exprstackIsEmpty(
324  	   EXPRSTACK*            exprstack           /**< expression stack */
325  	   )
326  	{
327  	   assert(exprstack != NULL);
328  	
329  	   return exprstack->stackpos < 0;
330  	}
331  	
332  	/** looks whether given expression is (proper) quadratic and has a given curvature
333  	 *
334  	 * If having a given curvature, currently require all arguments of quadratic to be linear.
335  	 * Hence, not using this for a simple square term, as curvCheckExprhdlr may provide a better condition on argument curvature then.
336  	 * Also we wouldn't do anything useful for a single bilinear term.
337  	 * Thus, run on sum's only.
338  	 */
339  	static
340  	DECL_CURVCHECK(curvCheckQuadratic)
341  	{  /*lint --e{715}*/
342  	   SCIP_EXPR* expr;
343  	   SCIP_EXPRCURV presentcurv;
344  	   SCIP_EXPRCURV wantedcurv;
345  	   SCIP_HASHSET* lonelysquares = NULL;
346  	   SCIP_Bool isquadratic;
347  	   int nbilinexprs;
348  	   int nquadexprs;
349  	   int i;
350  	
351  	   assert(nlexpr != NULL);
352  	   assert(stack != NULL);
353  	   assert(nlexpr2origexpr != NULL);
354  	   assert(success != NULL);
355  	
356  	   *success = FALSE;
357  	
358  	   if( !nlhdlrdata->cvxquadratic )
359  	      return SCIP_OKAY;
360  	
361  	   if( !SCIPisExprSum(scip, nlexpr) )
362  	      return SCIP_OKAY;
363  	
364  	   wantedcurv = SCIPexprGetCurvature(nlexpr);
365  	   if( wantedcurv == SCIP_EXPRCURV_LINEAR )
366  	      return SCIP_OKAY;
367  	   assert(wantedcurv == SCIP_EXPRCURV_CONVEX || wantedcurv == SCIP_EXPRCURV_CONCAVE);
368  	
369  	   expr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr);
370  	   assert(expr != NULL);
371  	
372  	   /* check whether quadratic */
373  	   SCIP_CALL( SCIPcheckExprQuadratic(scip, expr, &isquadratic) );
374  	
375  	   /* if not quadratic, then give up here */
376  	   if( !isquadratic )
377  	      return SCIP_OKAY;
378  	
379  	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, &nbilinexprs, NULL, NULL);
380  	
381  	   /* if only single square term (+linear), then give up here (let curvCheckExprhdlr handle this) */
382  	   if( nquadexprs <= 1 )
383  	      return SCIP_OKAY;
384  	
385  	   /* if root expression is only sum of squares (+linear) and detectsum is disabled, then give up here, too */
386  	   if( isrootexpr && !nlhdlrdata->detectsum && nbilinexprs == 0 )
387  	      return SCIP_OKAY;
388  	
389  	   /* get curvature of quadratic
390  	    * TODO as we know what curvature we want, we could first do some simple checks like computing xQx for a random x
391  	    */
392  	   SCIP_CALL( SCIPcomputeExprQuadraticCurvature(scip, expr, &presentcurv, assumevarfixed, FALSE) );
393  	
394  	   /* if not having desired curvature, return */
395  	   if( presentcurv != wantedcurv )
396  	      return SCIP_OKAY;
397  	
398  	   *success = TRUE;
399  	
400  	   if( !nlhdlrdata->detectsum )
401  	   {
402  	      /* first step towards block-decomposition of quadratic term:
403  	       * collect all square-expressions (in original expr) which have no adjacent bilinear term
404  	       * we will treat these x^2 as linear, i.e., add an auxvar for them, so x^2 maybe linearized
405  	       * more efficiently (in particular if x is discrete)
406  	       */
407  	      SCIP_CALL( SCIPhashsetCreate(&lonelysquares, SCIPblkmem(scip), nquadexprs) );
408  	      for( i = 0; i < nquadexprs; ++i )
409  	      {
410  	         int nadjbilin;
411  	         SCIP_EXPR* sqrexpr;
412  	
413  	         SCIPexprGetQuadraticQuadTerm(expr, i, NULL, NULL, NULL, &nadjbilin, NULL, &sqrexpr);
414  	         if( nadjbilin == 0 )
415  	         {
416  	            assert(sqrexpr != NULL);
417  	            SCIP_CALL( SCIPhashsetInsert(lonelysquares, SCIPblkmem(scip), (void*)sqrexpr) );
418  	         }
419  	      }
420  	   }
421  	
422  	   /* add immediate children to nlexpr */
423  	   SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, NULL) );
424  	   assert(SCIPexprGetNChildren(nlexpr) == SCIPexprGetNChildren(expr));
425  	
426  	   /* put children that are not square or product on stack
427  	    * grow child for children that are square or product and put this child on stack
428  	    * require all children to be linear
429  	    */
430  	   for( i = 0; i < SCIPexprGetNChildren(nlexpr); ++i )
431  	   {
432  	      SCIP_EXPR* child;
433  	      SCIP_EXPRCURV curvlinear[2] = { SCIP_EXPRCURV_LINEAR, SCIP_EXPRCURV_LINEAR };
434  	
435  	      child = SCIPexprGetChildren(nlexpr)[i];
436  	      assert(child != NULL);
437  	
438  	      assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)child) == SCIPexprGetChildren(expr)[i]);
439  	
440  	      if( SCIPisExprPower(scip, child) && SCIPgetExponentExprPow(child) == 2.0 &&
441  	         (lonelysquares == NULL || !SCIPhashsetExists(lonelysquares, SCIPexprGetChildren(expr)[i])) )
442  	      {
443  	         /* square term that isn't lonely, i.e., orig-version of child is a square-expr and nadjbilin>0 */
444  	         SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, child, curvlinear) );
445  	         assert(SCIPexprGetNChildren(child) == 1);
446  	         SCIP_CALL( exprstackPush(scip, stack, 1, SCIPexprGetChildren(child)) );
447  	      }
448  	      else if( SCIPisExprProduct(scip, child) && SCIPexprGetNChildren(SCIPexprGetChildren(expr)[i]) == 2 )
449  	         /* using original version of child here as NChildren(child)==0 atm */
450  	      {
451  	         /* bilinear term */
452  	         SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, child, curvlinear) );
453  	         assert(SCIPexprGetNChildren(child) == 2);
454  	         SCIP_CALL( exprstackPush(scip, stack, 2, SCIPexprGetChildren(child)) );
455  	      }
456  	      else
457  	      {
458  	         /* linear term (or term to be considered as linear) or lonely square term
459  	          * if we want extended formulations, then require linearity, so an auxvar will be introduced if it is nonlinear
460  	          * if we do not want extended formulations, then the term needs to have curvature "wantedcurv"
461  	          *   thus, if the coef is negative, then the child needs to have the curvature opposite to "wantedcurv"
462  	          */
463  	         if( nlhdlrdata->extendedform )
464  	            SCIPexprSetCurvature(child, SCIP_EXPRCURV_LINEAR);
465  	         else
466  	            SCIPexprSetCurvature(child, SCIPexprcurvMultiply(SCIPgetCoefsExprSum(nlexpr)[i], wantedcurv));
467  	         SCIP_CALL( exprstackPush(scip, stack, 1, &child) );
468  	      }
469  	   }
470  	
471  	   if( lonelysquares != NULL )
472  	      SCIPhashsetFree(&lonelysquares, SCIPblkmem(scip));
473  	
474  	   return SCIP_OKAY;
475  	}
476  	
477  	/** looks whether top of given expression looks like a signomial that can have a given curvature
478  	 *
479  	 * e.g., sqrt(x)*sqrt(y) is convex if x,y >= 0 and x and y are convex
480  	 *
481  	 * unfortunately, doesn't work for tls, because i) it's originally sqrt(x*y), and ii) it is expanded into some sqrt(z*y+y);
482  	 * but works for cvxnonsep_nsig
483  	 */
484  	static
485  	DECL_CURVCHECK(curvCheckSignomial)
486  	{  /*lint --e{715}*/
487  	   SCIP_EXPR* expr;
488  	   SCIP_EXPR* child;
489  	   SCIP_Real* exponents;
490  	   SCIP_INTERVAL* bounds;
491  	   SCIP_EXPRCURV* curv;
492  	   int nfactors;
493  	   int i;
494  	
495  	   assert(nlexpr != NULL);
496  	   assert(stack != NULL);
497  	   assert(nlexpr2origexpr != NULL);
498  	   assert(success != NULL);
499  	
500  	   *success = FALSE;
501  	
502  	   if( !nlhdlrdata->cvxsignomial )
503  	      return SCIP_OKAY;
504  	
505  	   if( !SCIPisExprProduct(scip, nlexpr) )
506  	      return SCIP_OKAY;
507  	
508  	   expr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr);
509  	   assert(expr != NULL);
510  	
511  	   nfactors = SCIPexprGetNChildren(expr);
512  	   if( nfactors <= 1 )  /* boooring */
513  	      return SCIP_OKAY;
514  	
515  	   SCIP_CALL( SCIPallocBufferArray(scip, &exponents, nfactors) );
516  	   SCIP_CALL( SCIPallocBufferArray(scip, &bounds, nfactors) );
517  	   SCIP_CALL( SCIPallocBufferArray(scip, &curv, nfactors) );
518  	
519  	   for( i = 0; i < nfactors; ++i )
520  	   {
521  	      child = SCIPexprGetChildren(expr)[i];
522  	      assert(child != NULL);
523  	
524  	      if( !SCIPisExprPower(scip, child) )
525  	      {
526  	         exponents[i] = 1.0;
527  	         SCIP_CALL( SCIPevalExprActivity(scip, child) );
528  	         bounds[i] = SCIPexprGetActivity(child);
529  	      }
530  	      else
531  	      {
532  	         exponents[i] = SCIPgetExponentExprPow(child);
533  	         SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(child)[0]) );
534  	         bounds[i] = SCIPexprGetActivity(SCIPexprGetChildren(child)[0]);
535  	      }
536  	   }
537  	
538  	   if( !SCIPexprcurvMonomialInv(SCIPexprcurvMultiply(SCIPgetCoefExprProduct(expr), SCIPexprGetCurvature(nlexpr)),
539  	       nfactors, exponents, bounds, curv) )
540  	      goto TERMINATE;
541  	
542  	   /* add immediate children to nlexpr
543  	    * some entries in curv actually apply to arguments of pow's, will correct this next
544  	    */
545  	   SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, curv) );
546  	   assert(SCIPexprGetNChildren(nlexpr) == nfactors);
547  	
548  	   /* put children that are not power on stack
549  	    * grow child for children that are power and put this child on stack
550  	    * if extendedform, then require children to be linear
551  	    * unless they are linear, an auxvar will be introduced for them and thus they will be handled as var here
552  	    */
553  	   for( i = 0; i < nfactors; ++i )
554  	   {
555  	      child = SCIPexprGetChildren(nlexpr)[i];
556  	      assert(child != NULL);
557  	
558  	      if( SCIPisExprPower(scip, child) )
559  	      {
560  	         SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, child, &curv[i]) );
561  	         assert(SCIPexprGetNChildren(child) == 1);
562  	         child = SCIPexprGetChildren(child)[0];
563  	      }
564  	      assert(SCIPexprGetNChildren(child) == 0);
565  	
566  	      if( nlhdlrdata->extendedform )
567  	      {
568  	         SCIPexprSetCurvature(child, SCIP_EXPRCURV_LINEAR);
569  	#ifdef SCIP_DEBUG
570  	         SCIPinfoMessage(scip, NULL, "Extendedform: Require linearity for ");
571  	         SCIPprintExpr(scip, child, NULL);
572  	         SCIPinfoMessage(scip, NULL, "\n");
573  	#endif
574  	      }
575  	
576  	      SCIP_CALL( exprstackPush(scip, stack, 1, &child) );
577  	   }
578  	
579  	   *success = TRUE;
580  	
581  	TERMINATE:
582  	   SCIPfreeBufferArray(scip, &curv);
583  	   SCIPfreeBufferArray(scip, &bounds);
584  	   SCIPfreeBufferArray(scip, &exponents);
585  	
586  	   return SCIP_OKAY;
587  	}
588  	
589  	/** looks for \f$f(c h(x)+d) h(x) \cdot \text{constant}\f$ and tries to conclude conditions on curvature
590  	 *
591  	 * Assume \f$h\f$ is univariate:
592  	 * - First derivative is \f$f'(c h + d) c h' h + f(c h + d) h'\f$.
593  	 * - Second derivative is \f{align}{&f''(c h + d) c h' c h' h + f'(c h + d) (c h'' h + c h' h') + f'(c h + d) c h' h' + f(c h + d) h'' \\
594  	 *   =& f''(c h + d) c^2 h'^2 h + f'(c h + d) c h'' h + 2 f'(c h + d) c h'^2 + f(c h + d) h''.\f}
595  	 *   Remove always positive factors leaves \f[f''(c h + d) h,\quad f'(c h + d) c h'' h,\quad f'(c h + d) c,\quad f(c h + d) h''.\f]
596  	 *   For convexity we want all these terms to be nonnegative. For concavity we want all of them to be nonpositive.
597  	 *   Note, that in each term either both \f$f'(c h + d)\f$ and \f$c\f$ occur, or none of them.
598  	 * - Thus, \f$f(c h(x) + d)h(x)\f$ is convex if \f$cf\f$ is monotonically increasing \f$(c f' \geq 0)\f$ and either
599  	 *   - \f$f\f$ is convex \f$(f'' \geq 0)\f$ and \f$h\f$ is nonnegative \f$(h \geq 0)\f$ and \f$h\f$ is convex \f$(h'' \geq 0)\f$ and [\f$f\f$ is nonnegative \f$(f \geq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$], or
600  	 *   - \f$f\f$ is concave \f$(f'' \leq 0)\f$ and \f$h\f$ is nonpositive \f$(h \leq 0)\f$ and \f$h\f$ is concave \f$(h'' \leq 0)\f$ and [\f$f\f$ is nonpositive \f$(f \leq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$].
601  	 * - Further, \f$f(c h(x) + d)h(x)\f$ is concave if \f$cf\f$ is monotonically decreasing \f$(c f' \leq 0)\f$ and either
602  	 *   - f is convex \f$(f'' \geq 0)\f$ and \f$h\f$ is nonpositive \f$(h \leq 0)\f$ and \f$h\f$ is concave \f$(h'' \leq 0)\f$ and [\f$f\f$ is nonnegative \f$(f \geq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$], or
603  	 *   - f is concave \f$(f'' \leq 0)\f$ and \f$h\f$ is nonnegative \f$(h >= 0)\f$ and \f$h\f$ is convex \f$(h'' \geq 0)\f$ and [\f$f\f$ is nonpositive \f$(f \leq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$].
604  	 *
605  	 * This should hold also for multivariate and linear \f$h\f$, as things are invariant under linear transformations.
606  	 * Similar to signomial, I'll assume that this will also hold for other multivariate \f$h\f$ (someone has a formal proof?).
607  	 */
608  	static
609  	DECL_CURVCHECK(curvCheckProductComposite)
610  	{  /*lint --e{715}*/
611  	   SCIP_EXPR* expr;
612  	   SCIP_EXPR* f;
613  	   SCIP_EXPR* h = NULL;
614  	   SCIP_Real c = 0.0;
615  	   SCIP_EXPR* ch = NULL; /* c * h */
616  	   SCIP_Real d;
617  	   SCIP_INTERVAL fbounds;
618  	   SCIP_INTERVAL hbounds;
619  	   SCIP_MONOTONE fmonotonicity;
620  	   SCIP_EXPRCURV desiredcurv;
621  	   SCIP_EXPRCURV hcurv;
622  	   SCIP_EXPRCURV dummy;
623  	   int fidx;
624  	
625  	   assert(nlexpr != NULL);
626  	   assert(stack != NULL);
627  	   assert(nlexpr2origexpr != NULL);
628  	   assert(success != NULL);
629  	
630  	   *success = FALSE;
631  	
632  	   if( !nlhdlrdata->cvxprodcomp )
633  	      return SCIP_OKAY;
634  	
635  	   if( !SCIPisExprProduct(scip, nlexpr) )
636  	      return SCIP_OKAY;
637  	
638  	   expr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr);
639  	   assert(expr != NULL);
640  	
641  	   if( SCIPexprGetNChildren(expr) != 2 )
642  	      return SCIP_OKAY;
643  	
644  	   /* check whether we have f(c * h(x)) * h(x) or h(x) * f(c * h(x)) */
645  	   for( fidx = 0; fidx <= 1; ++fidx )
646  	   {
647  	      f = SCIPexprGetChildren(expr)[fidx];
648  	
649  	      if( SCIPexprGetNChildren(f) != 1 )
650  	         continue;
651  	
652  	      ch = SCIPexprGetChildren(f)[0];
653  	      c = 1.0;
654  	      h = ch;
655  	
656  	      /* check whether ch is of the form c*h(x), then switch h to child ch */
657  	      if( SCIPisExprSum(scip, ch) && SCIPexprGetNChildren(ch) == 1 )
658  	      {
659  	         c = SCIPgetCoefsExprSum(ch)[0];
660  	         h = SCIPexprGetChildren(ch)[0];
661  	         assert(c != 1.0 || SCIPgetConstantExprSum(ch) != 0.0);  /* we could handle this, but it should have been simplified away */
662  	      }
663  	
664  	#ifndef NLHDLR_CONVEX_UNITTEST
665  	      /* can assume that duplicate subexpressions have been identified and comparing pointer is sufficient */
666  	      if( SCIPexprGetChildren(expr)[1-fidx] == h )
667  	#else
668  	      /* called from unittest -> duplicate subexpressions were not identified -> compare more expensively */
669  	      if( SCIPcompareExpr(scip, SCIPexprGetChildren(expr)[1-fidx], h) == 0 )
670  	#endif
671  	         break;
672  	   }
673  	   if( fidx == 2 )
674  	      return SCIP_OKAY;
675  	
676  	   /* constant of c*h(x)+d */
677  	   d = h != ch ? SCIPgetConstantExprSum(ch) : 0.0;
678  	
679  	#ifdef SCIP_MORE_DEBUG
680  	   SCIPinfoMessage(scip, NULL, "f(c*h+d)*h with f = %s, c = %g, d = %g, h = ", SCIPexprhdlrGetName(SCIPexprGetHdlr(f)), c, d);
681  	   SCIPprintExpr(scip, h, NULL);
682  	   SCIPinfoMessage(scip, NULL, "\n");
683  	#endif
684  	
685  	   assert(c != 0.0);
686  	
687  	   SCIP_CALL( SCIPevalExprActivity(scip, f) );
688  	   SCIP_CALL( SCIPevalExprActivity(scip, h) );
689  	   fbounds = SCIPexprGetActivity(f);
690  	   hbounds = SCIPexprGetActivity(h);
691  	
692  	   /* if h has mixed sign, then cannot conclude anything */
693  	   if( hbounds.inf < 0.0 && hbounds.sup > 0.0 )
694  	      return SCIP_OKAY;
695  	
696  	   /* If we have some convex or concave x*abs(c*x+d), then gradients at x=-d/c may be very wrong due to
697  	    * rounding errors and non-differentiability of abs() at zero (#3411). Therefore, we skip handling
698  	    * such expression in this nonlinear handler when one of the bounds of c*x+d is very close to zero.
699  	    * (If zero is in between the bounds of c*x+d, then the composition wouldn't be regarded as convex/concave anyway.)
700  	    */
701  	   if( SCIPisExprAbs(scip, f) && (SCIPisZero(scip, c*hbounds.inf+d) || SCIPisZero(scip, c*hbounds.sup+d)) )
702  	      return SCIP_OKAY;
703  	
704  	   SCIP_CALL( SCIPcallExprMonotonicity(scip, f, 0, &fmonotonicity) );
705  	
706  	   /* if f is not monotone, then cannot conclude anything */
707  	   if( fmonotonicity == SCIP_MONOTONE_UNKNOWN )
708  	      return SCIP_OKAY;
709  	
710  	   /* curvature we want to achieve (negate if product has negative coef) */
711  	   desiredcurv = SCIPexprcurvMultiply(SCIPgetCoefExprProduct(nlexpr), SCIPexprGetCurvature(nlexpr));
712  	
713  	   /* now check the conditions as stated above */
714  	   if( desiredcurv == SCIP_EXPRCURV_CONVEX )
715  	   {
716  	      /* f(c h(x)+d)h(x) is convex if c*f is monotonically increasing (c f' >= 0) and either
717  	      *   - f is convex (f'' >= 0) and h is nonnegative (h >= 0) and h is convex (h'' >= 0) and [f is nonnegative (f >= 0) or h is linear (h''=0)], or
718  	      *   - f is concave (f'' <= 0) and h is nonpositive (h <= 0) and h is concave (h'' <= 0) and [f is nonpositive (f <= 0) or h is linear (h''=0)]
719  	      *  as the curvature requirements on f are on f only and not the composition f(h), we can ignore the requirements returned by SCIPcallExprCurvature (last arg)
720  	      */
721  	      if( (c > 0.0 && fmonotonicity != SCIP_MONOTONE_INC) || (c < 0.0 && fmonotonicity != SCIP_MONOTONE_DEC) )
722  	         return SCIP_OKAY;
723  	
724  	      /* check whether f can be convex (h>=0) or concave (h<=0), resp., and derive requirements for h */
725  	      if( hbounds.inf >= 0 )
726  	      {
727  	         SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONVEX, success, &dummy) );
728  	
729  	         /* now h also needs to be convex; and if f < 0, then h actually needs to be linear */
730  	         if( fbounds.inf < 0.0 )
731  	            hcurv = SCIP_EXPRCURV_LINEAR;
732  	         else
733  	            hcurv = SCIP_EXPRCURV_CONVEX;
734  	      }
735  	      else
736  	      {
737  	         SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONCAVE, success, &dummy) );
738  	
739  	         /* now h also needs to be concave; and if f > 0, then h actually needs to be linear */
740  	         if( fbounds.sup > 0.0 )
741  	            hcurv = SCIP_EXPRCURV_LINEAR;
742  	         else
743  	            hcurv = SCIP_EXPRCURV_CONCAVE;
744  	      }
745  	   }
746  	   else
747  	   {
748  	      /* f(c h(x)+d)*h(x) is concave if c*f is monotonically decreasing (c f' <= 0) and either
749  	      *   - f is convex (f'' >= 0) and h is nonpositive (h <= 0) and h is concave (h'' <= 0) and [f is nonnegative (f >= 0) or h is linear (h''=0)], or
750  	      *   - f is concave (f'' <= 0) and h is nonnegative (h >= 0) and h is convex (h'' >= 0) and [f is nonpositive (f <= 0) or h is linear (h''=0)]
751  	      *  as the curvature requirements on f are on f only and not the composition f(h), we can ignore the requirements returned by SCIPcallExprCurvature (last arg)
752  	      */
753  	      if( (c > 0.0 && fmonotonicity != SCIP_MONOTONE_DEC) || (c < 0.0 && fmonotonicity != SCIP_MONOTONE_INC) )
754  	         return SCIP_OKAY;
755  	
756  	      /* check whether f can be convex (h<=0) or concave (h>=0), resp., and derive requirements for h */
757  	      if( hbounds.sup <= 0 )
758  	      {
759  	         SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONVEX, success, &dummy) );
760  	
761  	         /* now h also needs to be concave; and if f < 0, then h actually needs to be linear */
762  	         if( fbounds.inf < 0.0 )
763  	            hcurv = SCIP_EXPRCURV_LINEAR;
764  	         else
765  	            hcurv = SCIP_EXPRCURV_CONCAVE;
766  	      }
767  	      else
768  	      {
769  	         SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONCAVE, success, &dummy) );
770  	
771  	         /* now h also needs to be convex; and if f > 0, then h actually needs to be linear */
772  	         if( fbounds.sup > 0.0 )
773  	            hcurv = SCIP_EXPRCURV_LINEAR;
774  	         else
775  	            hcurv = SCIP_EXPRCURV_CONVEX;
776  	      }
777  	   }
778  	
779  	   if( !*success )
780  	      return SCIP_OKAY;
781  	
782  	   /* add immediate children (f and ch) to nlexpr; we set required curvature for h further below */
783  	   SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, NULL) );
784  	   assert(SCIPexprGetNChildren(nlexpr) == 2);
785  	
786  	   /* copy of f (and h) should have same child position in nlexpr as f (and h) has on expr (resp) */
787  	   assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)SCIPexprGetChildren(nlexpr)[fidx]) == (void*)f);
788  	#ifndef NLHDLR_CONVEX_UNITTEST
789  	   assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)SCIPexprGetChildren(nlexpr)[1-fidx]) == (void*)h);
790  	#endif
791  	   /* push this h onto stack for further checking */
792  	   SCIP_CALL( exprstackPush(scip, stack, 1, &(SCIPexprGetChildren(nlexpr)[1-fidx])) );
793  	
794  	   /* if we prefer extended formulations, then we always want h() to be linear */
795  	   if( nlhdlrdata->extendedform )
796  	      hcurv = SCIP_EXPRCURV_LINEAR;
797  	
798  	   /* h-child of product should have curvature hcurv */
799  	   SCIPexprSetCurvature(SCIPexprGetChildren(nlexpr)[1-fidx], hcurv);
800  	
801  	   if( h != ch )
802  	   {
803  	      /* add copy of ch as child to copy of f */
804  	      SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, SCIPexprGetChildren(nlexpr)[fidx], NULL) );
805  	      assert(SCIPexprGetNChildren(SCIPexprGetChildren(nlexpr)[fidx]) == 1);
806  	      assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)SCIPexprGetChildren(SCIPexprGetChildren(nlexpr)[fidx])[0]) == (void*)ch);
807  	
808  	      /* add copy of h (created above as child of product) as child in copy of ch */
809  	      SCIP_CALL( SCIPappendExprChild(scip,
810  	         SCIPexprGetChildren(SCIPexprGetChildren(nlexpr)[fidx])[0] /* copy of ch */,
811  	         SCIPexprGetChildren(nlexpr)[1-fidx] /* copy of h */) );
812  	   }
813  	   else
814  	   {
815  	      /* add copy of h (created above as child of product) as child in copy of f */
816  	      SCIP_CALL( SCIPappendExprChild(scip,
817  	         SCIPexprGetChildren(nlexpr)[fidx] /* copy of f */,
818  	         SCIPexprGetChildren(nlexpr)[1-fidx] /* copy of h */) );
819  	   }
820  	
821  	   return SCIP_OKAY;
822  	}
823  	
824  	/** use expression handlers curvature callback to check whether given curvature can be achieved */
825  	static
826  	DECL_CURVCHECK(curvCheckExprhdlr)
827  	{  /*lint --e{715}*/
828  	   SCIP_EXPR* origexpr;
829  	   int nchildren;
830  	   SCIP_EXPRCURV* childcurv;
831  	
832  	   assert(nlexpr != NULL);
833  	   assert(stack != NULL);
834  	   assert(nlexpr2origexpr != NULL);
835  	   assert(success != NULL);
836  	
837  	   origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, nlexpr);
838  	   assert(origexpr != NULL);
839  	   nchildren = SCIPexprGetNChildren(origexpr);
840  	
841  	   if( nchildren == 0 )
842  	   {
843  	      /* if originally no children, then should be var or value, which should have every curvature,
844  	       * so should always be success
845  	       */
846  	      SCIP_CALL( SCIPcallExprCurvature(scip, origexpr, SCIPexprGetCurvature(nlexpr), success, NULL) );
847  	      assert(*success);
848  	
849  	      return SCIP_OKAY;
850  	   }
851  	
852  	   /* ignore sums if > 1 children
853  	    * NOTE: this means that for something like 1+f(x), even if f is a trivial convex expression, we would handle 1+f(x)
854  	    * with this nlhdlr, instead of formulating this as 1+z and handling z=f(x) with the default nlhdlr, i.e., the exprhdlr
855  	    * today, I prefer handling this here, as it avoids introducing an extra auxiliary variable
856  	    */
857  	   if( isrootexpr && !nlhdlrdata->detectsum && SCIPisExprSum(scip, nlexpr) && nchildren > 1 )
858  	      return SCIP_OKAY;
859  	
860  	   SCIP_CALL( SCIPallocBufferArray(scip, &childcurv, nchildren) );
861  	
862  	   /* check whether and under which conditions origexpr can have desired curvature */
863  	   SCIP_CALL( SCIPcallExprCurvature(scip, origexpr, SCIPexprGetCurvature(nlexpr), success, childcurv) );
864  	#ifdef SCIP_MORE_DEBUG
865  	   SCIPprintExpr(scip, origexpr, NULL);
866  	   SCIPinfoMessage(scip, NULL, " is %s? %d\n", SCIPexprcurvGetName(SCIPexprGetCurvature(nlexpr)), *success);
867  	#endif
868  	   if( !*success )
869  	      goto TERMINATE;
870  	
871  	   /* if origexpr can have curvature curv, then don't treat it as leaf, but include its children */
872  	   SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, childcurv) );
873  	   assert(SCIPexprGetChildren(nlexpr) != NULL);
874  	   assert(SCIPexprGetNChildren(nlexpr) == nchildren);
875  	
876  	   /* If we prefer extended formulations, then require all children to be linear.
877  	    * Unless they are, auxvars will be introduced and they will be handles as variables, which can be an
878  	    * advantage in the context of extended formulations.
879  	    */
880  	   if( nlhdlrdata->extendedform )
881  	   {
882  	      int i;
883  	      for( i = 0; i < nchildren; ++i )
884  	         SCIPexprSetCurvature(SCIPexprGetChildren(nlexpr)[i], SCIP_EXPRCURV_LINEAR);
885  	#ifdef SCIP_DEBUG
886  	      SCIPinfoMessage(scip, NULL, "require linearity for children of ");
887  	      SCIPprintExpr(scip, origexpr, NULL);
888  	      SCIPinfoMessage(scip, NULL, "\n");
889  	#endif
890  	   }
891  	
892  	   /* add children expressions to to-do list (stack) */
893  	   SCIP_CALL( exprstackPush(scip, stack, nchildren, SCIPexprGetChildren(nlexpr)) );
894  	
895  	TERMINATE:
896  	   SCIPfreeBufferArray(scip, &childcurv);
897  	
898  	   return SCIP_OKAY;
899  	}
900  	
901  	/** curvature check and expression-growing methods
902  	 *
903  	 * some day this could be plugins added by users at runtime, but for now we have a fixed list here
904  	 * @note curvCheckExprhdlr should be last
905  	 */
906  	static DECL_CURVCHECK((*CURVCHECKS[])) = { curvCheckProductComposite, curvCheckSignomial, curvCheckQuadratic, curvCheckExprhdlr };
907  	/** number of curvcheck methods */
908  	static const int NCURVCHECKS = sizeof(CURVCHECKS) / sizeof(void*);
909  	
910  	/** checks whether expression is a sum with more than one child and each child being a variable or going to be a variable if `expr` is a nlhdlr-specific copy
911  	 *
912  	 * Within constructExpr(), we can have an expression of any type which is a copy of an original expression,
913  	 * but without children. At the end of constructExpr() (after the loop with the stack), these expressions
914  	 * will remain as leafs and will eventually be turned into variables in collectLeafs(). Thus, we treat
915  	 * every child that has no children as if it were a variable. Theoretically, there is still the possibility
916  	 * that it could be a constant (value-expression), but simplify should have removed these.
917  	 */
918  	static
919  	SCIP_Bool exprIsMultivarLinear(
920  	   SCIP*                 scip,               /**< SCIP data structure */
921  	   SCIP_EXPR*            expr                /**< expression to check */
922  	   )
923  	{
924  	   int nchildren;
925  	   int c;
926  	
927  	   assert(expr != NULL);
928  	
929  	   if( !SCIPisExprSum(scip, expr) )
930  	      return FALSE;
931  	
932  	   nchildren = SCIPexprGetNChildren(expr);
933  	   if( nchildren <= 1 )
934  	      return FALSE;
935  	
936  	   for( c = 0; c < nchildren; ++c )
937  	      /*if( !SCIPisExprVar(scip, SCIPexprGetChildren(expr)[c]) ) */
938  	      if( SCIPexprGetNChildren(SCIPexprGetChildren(expr)[c]) > 0 )
939  	         return FALSE;
940  	
941  	   return TRUE;
942  	}
943  	
944  	/** constructs a subexpression (as nlhdlr-expression) of maximal size that has a given curvature
945  	 *
946  	 * If the curvature cannot be achieved for an expression in the original expression graph,
947  	 * then this expression becomes a leaf in the nlhdlr-expression.
948  	 *
949  	 * Sets `*rootnlexpr` to NULL if failed.
950  	 */
951  	static
952  	SCIP_RETCODE constructExpr(
953  	   SCIP*                 scip,               /**< SCIP data structure */
954  	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nonlinear handler data */
955  	   SCIP_EXPR**           rootnlexpr,         /**< buffer to store created expression */
956  	   SCIP_HASHMAP*         nlexpr2origexpr,    /**< mapping from our expression copy to original expression */
957  	   int*                  nleafs,             /**< number of leafs in constructed expression */
958  	   SCIP_EXPR*            rootexpr,           /**< expression */
959  	   SCIP_EXPRCURV         curv,               /**< curvature to achieve */
960  	   SCIP_HASHMAP*         assumevarfixed,     /**< hashmap containing variables that should be assumed to be fixed, or NULL */
961  	   SCIP_Bool             assumecurvature,    /**< whether to assume that desired curvature is given (skips curvature checks) */
962  	   SCIP_Bool*            curvsuccess         /**< pointer to store whether the curvature could be achieved
963  	                                              *   w.r.t. the original variables (might be NULL) */
964  	   )
965  	{
966  	   SCIP_EXPR* nlexpr;
967  	   EXPRSTACK stack; /* to do list: expressions where to check whether they can have the desired curvature when taking their children into account */
968  	   int oldstackpos;
969  	   SCIP_Bool isrootexpr = TRUE;
970  	
971  	   assert(scip != NULL);
972  	   assert(nlhdlrdata != NULL);
973  	   assert(rootnlexpr != NULL);
974  	   assert(nlexpr2origexpr != NULL);
975  	   assert(nleafs != NULL);
976  	   assert(rootexpr != NULL);
977  	   assert(curv == SCIP_EXPRCURV_CONVEX || curv == SCIP_EXPRCURV_CONCAVE);
978  	
979  	   /* create root expression */
980  	   SCIP_CALL( nlhdlrExprCreate(scip, nlexpr2origexpr, rootnlexpr, rootexpr, curv) );
981  	
982  	   *nleafs = 0;
983  	   if( curvsuccess != NULL )
984  	      *curvsuccess = TRUE;
985  	
986  	   SCIP_CALL( exprstackInit(scip, &stack, 20) );
987  	   SCIP_CALL( exprstackPush(scip, &stack, 1, rootnlexpr) );
988  	   while( !exprstackIsEmpty(&stack) )
989  	   {
990  	      /* take expression from stack */
991  	      nlexpr = exprstackPop(&stack);
992  	      assert(nlexpr != NULL);
993  	      assert(SCIPexprGetNChildren(nlexpr) == 0);
994  	
995  	      oldstackpos = stack.stackpos;
996  	      if( nlhdlrdata->isnlhdlrconvex && !SCIPexprhdlrHasBwdiff(SCIPexprGetHdlr(nlexpr)) )
997  	      {
998  	         /* if bwdiff is not implemented, then we could not generate cuts in the convex nlhdlr, so "stop" (treat nlexpr as variable) */
999  	      }
1000 	      else if( !nlhdlrdata->isnlhdlrconvex && exprIsMultivarLinear(scip, (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr)) )
1001 	      {
1002 	         /* if we are in the concave handler, we would like to treat linear multivariate subexpressions by a new auxvar always,
1003 	          * e.g., handle log(x+y) as log(z), z=x+y, because the estimation problem will be smaller then without making the estimator worse
1004 	          * (cons_nonlinear does this, too)
1005 	          * this check takes care of this when x and y are original variables
1006 	          * however, it isn't unlikely that we will have sums that become linear after we add auxvars for some children
1007 	          * this will be handled in a postprocessing below
1008 	          * for now, the check is performed on the original expression since there is not enough information in nlexpr yet
1009 	          */
1010 	#ifdef SCIP_MORE_DEBUG
1011 	         SCIPprintExpr(scip, SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr), NULL);
1012 	         SCIPinfoMessage(scip, NULL, "... is a multivariate linear sum that we'll treat as auxvar\n");
1013 	#endif
1014 	      }
1015 	      else if( SCIPexprGetCurvature(nlexpr) != SCIP_EXPRCURV_UNKNOWN && !assumecurvature )
1016 	      {
1017 	         /* if we are here, either convexity or concavity is required; try to check for this curvature */
1018 	         SCIP_Bool success;
1019 	         int method;
1020 	
1021 	         /* try through curvature check methods until one succeeds */
1022 	         for( method = 0; method < NCURVCHECKS; ++method )
1023 	         {
1024 	            SCIP_CALL( CURVCHECKS[method](scip, nlexpr, isrootexpr, &stack, nlexpr2origexpr, nlhdlrdata, assumevarfixed, &success) );
1025 	            if( success )
1026 	               break;
1027 	         }
1028 	      }
1029 	      else
1030 	      {
1031 	         /* if we don't care about curvature in this subtree anymore (very unlikely),
1032 	          * or we are told to assume that the desired curvature is present (assumecurvature==TRUE),
1033 	          * then only continue iterating this subtree to assemble leaf expressions
1034 	          */
1035 	         SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, NULL) );
1036 	
1037 	         /* add children expressions, if any, to to-do list (stack) */
1038 	         SCIP_CALL( exprstackPush(scip, &stack, SCIPexprGetNChildren(nlexpr), SCIPexprGetChildren(nlexpr)) );
1039 	      }
1040 	      assert(stack.stackpos >= oldstackpos);  /* none of the methods above should have removed something from the stack */
1041 	
1042 	      isrootexpr = FALSE;
1043 	
1044 	      /* if nothing was added, then none of the successors of nlexpr were added to the stack
1045 	       * this is either because nlexpr was already a variable or value expressions, thus a leaf,
1046 	       * or because the desired curvature could not be achieved, so it will be handled as variables, thus a leaf
1047 	       */
1048 	      if( stack.stackpos == oldstackpos )
1049 	      {
1050 	         ++*nleafs;
1051 	
1052 	         /* check whether the new leaf is not an original variable (or constant) */
1053 	         if( curvsuccess != NULL && !SCIPisExprVar(scip, nlexpr) && !SCIPisExprValue(scip, nlexpr) )
1054 	            *curvsuccess = FALSE;
1055 	      }
1056 	   }
1057 	
1058 	   exprstackFree(scip, &stack);
1059 	
1060 	   if( !nlhdlrdata->isnlhdlrconvex && *rootnlexpr != NULL )
1061 	   {
1062 	      /* remove multivariate linear subexpressions, that is, change some f(z1+z2) into f(z3) (z3=z1+z2 will be done by nlhdlr_default)
1063 	       * this handles the case that was not covered by the above check, which could recognize f(x+y) for x, y original variables
1064 	       */
1065 	      SCIP_EXPRITER* it;
1066 	
1067 	      SCIP_CALL( SCIPcreateExpriter(scip, &it) );
1068 	      SCIP_CALL( SCIPexpriterInit(it, *rootnlexpr, SCIP_EXPRITER_DFS, FALSE) );
1069 	      SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD);
1070 	
1071 	      while( !SCIPexpriterIsEnd(it) )
1072 	      {
1073 	         SCIP_EXPR* child;
1074 	
1075 	         child = SCIPexpriterGetChildExprDFS(it);
1076 	         assert(child != NULL);
1077 	
1078 	         /* We want to change some f(x+y+z) into just f(), where f is the expression the iterator points to
1079 	          * and x+y+z is child. A child of a child, e.g., z, may not be a variable yet (these are added in collectLeafs later),
1080 	          * but an expression of some nonlinear type without children.
1081 	          */
1082 	         if( exprIsMultivarLinear(scip, child) )
1083 	         {
1084 	            /* turn child (x+y+z) into a sum without children
1085 	             * collectLeafs() should then replace this by an auxvar
1086 	             */
1087 	#ifdef SCIP_MORE_DEBUG
1088 	            SCIPprintExpr(scip, child, NULL);
1089 	            SCIPinfoMessage(scip, NULL, "... is a multivariate linear sum that we'll treat as auxvar instead (postprocess)\n");
1090 	#endif
1091 	
1092 	            /* TODO remove children from nlexpr2origexpr ?
1093 	             * should also do this if they are not used somewhere else; we could check nuses for this
1094 	             * however, it shouldn't matter to have some stray entries in the hashmap either
1095 	             */
1096 	            SCIP_CALL( SCIPremoveExprChildren(scip, child) );
1097 	            assert(SCIPexprGetNChildren(child) == 0);
1098 	
1099 	            (void) SCIPexpriterSkipDFS(it);
1100 	         }
1101 	         else
1102 	         {
1103 	            (void) SCIPexpriterGetNext(it);
1104 	         }
1105 	      }
1106 	
1107 	      SCIPfreeExpriter(&it);
1108 	   }
1109 	
1110 	   if( *rootnlexpr != NULL )
1111 	   {
1112 	      SCIP_Bool istrivial = TRUE;
1113 	
1114 	      /* if handletrivial is enabled, then only require that rootnlexpr itself has required curvature (so has children; see below) and
1115 	       * that we are not a trivial sum  (because the previous implementation of this nlhdlr didn't allow this, either)
1116 	       */
1117 	      if( !nlhdlrdata->handletrivial || SCIPisExprSum(scip, *rootnlexpr) )
1118 	      {
1119 	         /* if all children do not have children, i.e., are variables, or will be replaced by auxvars, then free
1120 	          * also if rootnlexpr has no children, then free
1121 	          */
1122 	         int i;
1123 	         for( i = 0; i < SCIPexprGetNChildren(*rootnlexpr); ++i )
1124 	         {
1125 	            if( SCIPexprGetNChildren(SCIPexprGetChildren(*rootnlexpr)[i]) > 0 )
1126 	            {
1127 	               istrivial = FALSE;
1128 	               break;
1129 	            }
1130 	         }
1131 	      }
1132 	      else if( SCIPexprGetNChildren(*rootnlexpr) > 0 )  /* if handletrivial, then just require children */
1133 	            istrivial = FALSE;
1134 	
1135 	      if( istrivial )
1136 	      {
1137 	         SCIP_CALL( SCIPreleaseExpr(scip, rootnlexpr) );
1138 	      }
1139 	   }
1140 	
1141 	   return SCIP_OKAY;
1142 	}
1143 	
1144 	/** collects (non-value) leaf expressions and ensure that they correspond to a variable (original or auxiliary)
1145 	 *
1146 	 * For children where we could not achieve the desired curvature, get the auxvar and replace the child by a
1147 	 * var-expression that points to this auxvar.
1148 	 * Collect all leaf expressions (if not a value-expression) and index them.
1149 	 */
1150 	static
1151 	SCIP_RETCODE collectLeafs(
1152 	   SCIP*                 scip,               /**< SCIP data structure */
1153 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata      /**< nlhdlr expression data */
1154 	   )
1155 	{
1156 	   SCIP_EXPRITER* it;
1157 	   SCIP_EXPR* nlexpr;
1158 	   SCIP_HASHMAP* leaf2index;
1159 	   int i;
1160 	
1161 	   assert(nlhdlrexprdata != NULL);
1162 	   assert(nlhdlrexprdata->nlexpr != NULL);
1163 	   assert(nlhdlrexprdata->nlexpr2origexpr != NULL);
1164 	   /* nleafs should be the upper bound on the number of variables given by constructExpr
1165 	    * leafexprs should be NULL, as this is what we want to setup here
1166 	    */
1167 	   assert(nlhdlrexprdata->nleafs > 0);
1168 	   assert(nlhdlrexprdata->leafexprs == NULL);
1169 	
1170 	   /* collect all auxvars and collect all variables */
1171 	   SCIP_CALL( SCIPhashmapCreate(&leaf2index, SCIPblkmem(scip), nlhdlrexprdata->nleafs) );
1172 	   nlhdlrexprdata->nleafs = 0;  /* we start a new count, this time skipping value-expressions */
1173 	
1174 	   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
1175 	   SCIP_CALL( SCIPexpriterInit(it, nlhdlrexprdata->nlexpr, SCIP_EXPRITER_DFS, FALSE) );
1176 	   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD);
1177 	
1178 	   for( nlexpr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); nlexpr = SCIPexpriterGetNext(it) )
1179 	   {
1180 	      SCIP_EXPR* child;
1181 	      SCIP_EXPR* origexpr;
1182 	
1183 	      assert(nlexpr != NULL);
1184 	
1185 	      child = SCIPexpriterGetChildExprDFS(it);
1186 	
1187 	      /* if the to-be-visited child has children, then it doesn't need to be replaced by a new expression (representing the auxvar) */
1188 	      if( SCIPexprGetNChildren(child) > 0 )
1189 	         continue;
1190 	
1191 	      origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)child);
1192 	      assert(origexpr != NULL);
1193 	
1194 	      if( SCIPexprGetNChildren(origexpr) > 0 )
1195 	      {
1196 	         SCIP_EXPR* newchild;
1197 	         int childidx;
1198 	         SCIP_VAR* var;
1199 	
1200 	         /* having a child that had children in original but not in copy means that we could not achieve the desired curvature
1201 	          * thus, replace by a new child that points to the auxvar of the original expression
1202 	          * we registered in createNlhdlrExprData that we need an auxvar, so it should exist now
1203 	          */
1204 	         var = SCIPgetExprAuxVarNonlinear(origexpr);
1205 	         assert(var != NULL);
1206 	
1207 	         SCIP_CALL( SCIPcreateExprVar(scip, &newchild, var, NULL, NULL) );  /* this captures newchild once */
1208 	
1209 	         childidx = SCIPexpriterGetChildIdxDFS(it);
1210 	         SCIP_CALL( SCIPreplaceExprChild(scip, nlexpr, childidx, newchild) );  /* this captures newchild again */
1211 	
1212 	         /* do not remove child->origexpr from hashmap, as child may appear again due to common subexprs
1213 	          * (created by curvCheckProductComposite, for example)
1214 	          * if it doesn't reappear, though, but the memory address is reused, we need to make sure it
1215 	          * points to the right origexpr
1216 	          */
1217 	         /* SCIP_CALL( SCIPhashmapRemove(nlexpr2origexpr, (void*)child) ); */
1218 	         SCIP_CALL( SCIPhashmapSetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)newchild, (void*)origexpr) );
1219 	
1220 	         if( !SCIPhashmapExists(leaf2index, (void*)newchild) )
1221 	         {
1222 	            /* new leaf -> new index and remember in hashmap */
1223 	            SCIP_CALL( SCIPhashmapInsertInt(leaf2index, (void*)newchild, nlhdlrexprdata->nleafs++) );
1224 	         }
1225 	
1226 	         child = newchild;
1227 	         SCIP_CALL( SCIPreleaseExpr(scip, &newchild) );  /* because it was captured by both create and replace */
1228 	      }
1229 	      else if( SCIPisExprVar(scip, child) )
1230 	      {
1231 	         /* if variable, then add to hashmap, if not already there */
1232 	         if( !SCIPhashmapExists(leaf2index, (void*)child) )
1233 	         {
1234 	            SCIP_CALL( SCIPhashmapInsertInt(leaf2index, (void*)child, nlhdlrexprdata->nleafs++) );
1235 	         }
1236 	      }
1237 	      /* else: it's probably a value-expression, nothing to do */
1238 	
1239 	      /* update integrality flag for future leaf expressions: convex nlhdlr may use this information */
1240 	      SCIP_CALL( SCIPcomputeExprIntegrality(scip, child) );
1241 	   }
1242 	   assert(nlhdlrexprdata->nleafs > 0);
1243 	
1244 	   SCIPfreeExpriter(&it);
1245 	
1246 	   /* assemble auxvars array */
1247 	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(nlhdlrexprdata->leafexprs), nlhdlrexprdata->nleafs) );
1248 	   for( i = 0; i < SCIPhashmapGetNEntries(leaf2index); ++i )
1249 	   {
1250 	      SCIP_HASHMAPENTRY* entry;
1251 	      SCIP_EXPR* leaf;
1252 	      int idx;
1253 	
1254 	      entry = SCIPhashmapGetEntry(leaf2index, i);
1255 	      if( entry == NULL )
1256 	         continue;
1257 	
1258 	      leaf = (SCIP_EXPR*) SCIPhashmapEntryGetOrigin(entry);
1259 	      assert(leaf != NULL);
1260 	      assert(SCIPisExprVar(scip, leaf));
1261 	
1262 	      idx = SCIPhashmapEntryGetImageInt(entry);
1263 	      assert(idx >= 0);
1264 	      assert(idx < nlhdlrexprdata->nleafs);
1265 	
1266 	      nlhdlrexprdata->leafexprs[idx] = leaf;
1267 	
1268 	      SCIPdebugMsg(scip, "leaf %d: <%s>\n", idx, SCIPvarGetName(SCIPgetVarExprVar(leaf)));
1269 	   }
1270 	
1271 	   SCIPhashmapFree(&leaf2index);
1272 	
1273 	   return SCIP_OKAY;
1274 	}
1275 	
1276 	/** creates nonlinear handler expression data structure and registers expr usage */
1277 	static
1278 	SCIP_RETCODE createNlhdlrExprData(
1279 	   SCIP*                 scip,               /**< SCIP data structure */
1280 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
1281 	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< pointer to store nlhdlr expression data */
1282 	   SCIP_EXPR*            expr,               /**< original expression */
1283 	   SCIP_EXPR*            nlexpr,             /**< our copy of expression */
1284 	   SCIP_HASHMAP*         nlexpr2origexpr,    /**< mapping of expression copy to original */
1285 	   int                   nleafs,             /**< number of leafs as counted by constructExpr */
1286 	   SCIP_NLHDLR_METHOD    participating       /**< the enfo methods in which we plan to participate */
1287 	   )
1288 	{
1289 	   SCIP_EXPRITER* it;
1290 	   SCIP_Bool usingaux;
1291 	
1292 	   assert(scip != NULL);
1293 	   assert(expr != NULL);
1294 	   assert(nlhdlrexprdata != NULL);
1295 	   assert(*nlhdlrexprdata == NULL);
1296 	   assert(nlexpr != NULL);
1297 	   assert(nlexpr2origexpr != NULL);
1298 	
1299 	   assert(SCIPexprGetNChildren(nlexpr) > 0);
1300 	   assert(SCIPexprGetChildren(nlexpr) != NULL);
1301 	
1302 	   SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) );
1303 	   (*nlhdlrexprdata)->nlexpr = nlexpr;
1304 	   (*nlhdlrexprdata)->nlexpr2origexpr = nlexpr2origexpr;
1305 	   (*nlhdlrexprdata)->nleafs = nleafs;
1306 	
1307 	   usingaux = FALSE;
1308 	
1309 	   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
1310 	   SCIP_CALL( SCIPexpriterInit(it, nlexpr, SCIP_EXPRITER_DFS, FALSE) );
1311 	   SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD);
1312 	
1313 	   for( ; !SCIPexpriterIsEnd(it); (void) SCIPexpriterGetNext(it) )
1314 	   {
1315 	      SCIP_EXPR* child;
1316 	      SCIP_EXPR* origexpr;
1317 	
1318 	      /* check whether to-be-visited child needs to be replaced by a new expression (representing the auxvar)
1319 	       * if child has children, then that is not the case
1320 	       * if child has no children, but also corresponding origexpr has no chilren, then this is also not the case
1321 	       */
1322 	      child = SCIPexpriterGetChildExprDFS(it);
1323 	      if( SCIPexprGetNChildren(child) > 0 )
1324 	         continue;
1325 	
1326 	      origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)child);
1327 	      assert(origexpr != NULL);
1328 	
1329 	      /* if child had children in original but not in copy means that we could not achieve the desired curvature
1330 	       * thus, we will later replace by a new child that points to the auxvar of the original expression
1331 	       * as we do not have the auxvar now, we will only register that we will need the auxvar later (if origexpr isn't a variable or constant)
1332 	       * if we are working for the concave nlhdlr, then we also indicate interest on the exprs activity for estimate (distinguish below or above)
1333 	       */
1334 	      SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, origexpr,
1335 	         SCIPexprGetNChildren(origexpr) > 0, FALSE,
1336 	         !nlhdlrdata->isnlhdlrconvex && (participating & SCIP_NLHDLR_METHOD_SEPABELOW),
1337 	         !nlhdlrdata->isnlhdlrconvex && (participating & SCIP_NLHDLR_METHOD_SEPAABOVE)) );
1338 	
1339 	      /* remember that we use an auxvar */
1340 	      if( SCIPexprGetNChildren(origexpr) > 0 )
1341 	         usingaux = TRUE;
1342 	   }
1343 	
1344 	   SCIPfreeExpriter(&it);
1345 	
1346 	#ifdef SCIP_DEBUG
1347 	   SCIPprintExpr(scip, nlexpr, NULL);
1348 	   SCIPinfoMessage(scip, NULL, " (%p) is handled as %s\n", SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr),
1349 	         SCIPexprcurvGetName(SCIPexprGetCurvature(nlexpr)));
1350 	#endif
1351 	
1352 	   /* If we don't work on the extended formulation, then set curvature also in original expression
1353 	    * (in case someone wants to pick this up; this might be removed again).
1354 	    * This doesn't ensure that every convex or concave original expression is actually marked here.
1355 	    * Not only because our tests are incomprehensive, but also because we may not detect on sums,
1356 	    * prefer extended formulations (in nlhdlr_convex), or introduce auxvars for linear subexpressions
1357 	    * on purpose (in nlhdlr_concave).
1358 	    */
1359 	   if( !usingaux )
1360 	      SCIPexprSetCurvature(expr, SCIPexprGetCurvature(nlexpr));
1361 	
1362 	   return SCIP_OKAY;
1363 	}
1364 	
1365 	/** adds an estimator for a vertex-polyhedral (e.g., concave) function to a given rowprep
1366 	 *
1367 	 * Calls \ref SCIPcomputeFacetVertexPolyhedralNonlinear() for given function and
1368 	 * box set to local bounds of auxiliary variables.
1369 	 */
1370 	static
1371 	SCIP_RETCODE estimateVertexPolyhedral(
1372 	   SCIP*                 scip,               /**< SCIP data structure */
1373 	   SCIP_CONSHDLR*        conshdlr,           /**< nonlinear constraint handler */
1374 	   SCIP_NLHDLR*          nlhdlr,             /**< nonlinear handler */
1375 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
1376 	   SCIP_SOL*             sol,                /**< solution to use, unless usemidpoint is TRUE */
1377 	   SCIP_Bool             usemidpoint,        /**< whether to use the midpoint of the domain instead of sol */
1378 	   SCIP_Bool             overestimate,       /**< whether over- or underestimating */
1379 	   SCIP_Real             targetvalue,        /**< a target value to achieve; if not reachable, then can give up early */
1380 	   SCIP_ROWPREP*         rowprep,            /**< rowprep where to store estimator */
1381 	   SCIP_Bool*            success             /**< buffer to store whether successful */
1382 	   )
1383 	{
1384 	   SCIP_NLHDLRDATA* nlhdlrdata;
1385 	   VERTEXPOLYFUN_EVALDATA evaldata;
1386 	   SCIP_Real* xstar;
1387 	   SCIP_Real* box;
1388 	   SCIP_Real facetconstant;
1389 	   SCIP_VAR* var;
1390 	   int i;
1391 	   SCIP_Bool allfixed;
1392 	
1393 	   assert(scip != NULL);
1394 	   assert(nlhdlr != NULL);
1395 	   assert(nlhdlrexprdata != NULL);
1396 	   assert(rowprep != NULL);
1397 	   assert(success != NULL);
1398 	
1399 	   *success = FALSE;
1400 	
1401 	   /* caller is responsible to have checked whether we can estimate, i.e., expression curvature and overestimate flag match */
1402 	   assert( overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE);  /* if underestimate, then must be concave */
1403 	   assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX);   /* if overestimate, then must be convex */
1404 	
1405 	#ifdef SCIP_DEBUG
1406 	   SCIPinfoMessage(scip, NULL, "%sestimate expression ", overestimate ? "over" : "under");
1407 	   SCIPprintExpr(scip, nlhdlrexprdata->nlexpr, NULL);
1408 	   SCIPinfoMessage(scip, NULL, " at point\n");
1409 	   for( i = 0; i < nlhdlrexprdata->nleafs; ++i )
1410 	   {
1411 	      var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]);
1412 	      assert(var != NULL);
1413 	
1414 	      SCIPinfoMessage(scip, NULL, "  <%s> = %g [%g,%g]\n", SCIPvarGetName(var),
1415 	         usemidpoint ? 0.5 * (SCIPvarGetLbLocal(var) + SCIPvarGetUbLocal(var)) : SCIPgetSolVal(scip, sol, var),
1416 	        SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var));
1417 	   }
1418 	#endif
1419 	
1420 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1421 	   assert(nlhdlrdata != NULL);
1422 	
1423 	   if( nlhdlrdata->evalsol == NULL )
1424 	   {
1425 	      SCIP_CALL( SCIPcreateSol(scip, &nlhdlrdata->evalsol, NULL) );
1426 	   }
1427 	
1428 	   evaldata.nlhdlrexprdata = nlhdlrexprdata;
1429 	   evaldata.evalsol = nlhdlrdata->evalsol;
1430 	   evaldata.scip = scip;
1431 	
1432 	   SCIP_CALL( SCIPallocBufferArray(scip, &xstar, nlhdlrexprdata->nleafs) );
1433 	   SCIP_CALL( SCIPallocBufferArray(scip, &box, 2*nlhdlrexprdata->nleafs) );
1434 	
1435 	   allfixed = TRUE;
1436 	   for( i = 0; i < nlhdlrexprdata->nleafs; ++i )
1437 	   {
1438 	      var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]);
1439 	      assert(var != NULL);
1440 	
1441 	      box[2*i] = SCIPvarGetLbLocal(var);
1442 	      if( SCIPisInfinity(scip, -box[2*i]) )
1443 	      {
1444 	         SCIPdebugMsg(scip, "lower bound at -infinity, no estimate possible\n");
1445 	         goto TERMINATE;
1446 	      }
1447 	
1448 	      box[2*i+1] = SCIPvarGetUbLocal(var);
1449 	      if( SCIPisInfinity(scip, box[2*i+1]) )
1450 	      {
1451 	         SCIPdebugMsg(scip, "upper bound at +infinity, no estimate possible\n");
1452 	         goto TERMINATE;
1453 	      }
1454 	
1455 	      if( !SCIPisRelEQ(scip, box[2*i], box[2*i+1]) )
1456 	         allfixed = FALSE;
1457 	
1458 	      if( usemidpoint )
1459 	         xstar[i] = 0.5 * (box[2*i] + box[2*i+1]);
1460 	      else
1461 	         xstar[i] = SCIPgetSolVal(scip, sol, var);
1462 	      assert(xstar[i] != SCIP_INVALID);
1463 	   }
1464 	
1465 	   if( allfixed )
1466 	   {
1467 	      /* SCIPcomputeFacetVertexPolyhedralNonlinear prints a warning and does not succeed if all is fixed */
1468 	      SCIPdebugMsg(scip, "all variables fixed, skip estimate\n");
1469 	      goto TERMINATE;
1470 	   }
1471 	
1472 	   SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, nlhdlrexprdata->nleafs + 1) );
1473 	
1474 	   SCIP_CALL( SCIPcomputeFacetVertexPolyhedralNonlinear(scip, conshdlr, overestimate, nlhdlrExprEvalConcave, (void*)&evaldata,
1475 	      xstar, box, nlhdlrexprdata->nleafs, targetvalue, success, SCIProwprepGetCoefs(rowprep), &facetconstant) );
1476 	
1477 	   if( !*success )
1478 	   {
1479 	      SCIPdebugMsg(scip, "failed to compute facet of convex hull\n");
1480 	      goto TERMINATE;
1481 	   }
1482 	
1483 	   SCIProwprepSetLocal(rowprep, TRUE);
1484 	   SCIProwprepAddConstant(rowprep, facetconstant);
1485 	   for( i = 0; i < nlhdlrexprdata->nleafs; ++i )
1486 	   {
1487 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]), SCIProwprepGetCoefs(rowprep)[i]) );
1488 	   }
1489 	
1490 	#ifdef SCIP_DEBUG
1491 	   SCIPinfoMessage(scip, NULL, "computed estimator: ");
1492 	   SCIPprintRowprep(scip, rowprep, NULL);
1493 	#endif
1494 	
1495 	 TERMINATE:
1496 	   SCIPfreeBufferArray(scip, &box);
1497 	   SCIPfreeBufferArray(scip, &xstar);
1498 	
1499 	   return SCIP_OKAY;
1500 	}
1501 	
1502 	/** adds an estimator computed via a gradient to a given rowprep */
1503 	static
1504 	SCIP_RETCODE estimateGradientInner(
1505 	   SCIP*                 scip,               /**< SCIP data structure */
1506 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
1507 	   SCIP_SOL*             sol,                /**< solution to use */
1508 	   SCIP_ROWPREP*         rowprep,            /**< rowprep where to store estimator */
1509 	   SCIP_Bool*            success             /**< buffer to store whether successful */
1510 	   )
1511 	{
1512 	   SCIP_EXPR* nlexpr;
1513 	   SCIP_Real QUAD(constant);
1514 	   int i;
1515 	
1516 	   assert(scip != NULL);
1517 	   assert(nlhdlrexprdata != NULL);
1518 	   assert(rowprep != NULL);
1519 	   assert(success != NULL);
1520 	
1521 	   nlexpr = nlhdlrexprdata->nlexpr;
1522 	   assert(nlexpr != NULL);
1523 	
1524 	   /* compute gradient (TODO: this also re-evaluates (soltag=0), which shouldn't be necessary unless we tried ConvexSecant before or are called from Sollinearize callback) */
1525 	   SCIP_CALL( SCIPevalExprGradient(scip, nlexpr, sol, 0L) );
1526 	
1527 	   /* if gradient evaluation error, then return */
1528 	   if( SCIPexprGetDerivative(nlexpr) == SCIP_INVALID )
1529 	   {
1530 	      SCIPdebugMsg(scip, "gradient evaluation error for %p\n", (void*)nlexpr);
1531 	      return SCIP_OKAY;
1532 	   }
1533 	
1534 	   /* add gradient underestimator to rowprep: f(sol) + (x - sol) \nabla f(sol)
1535 	    * constant will store f(sol) - sol * \nabla f(sol)
1536 	    * to avoid some cancellation errors when linear variables take huge values (like 1e20),
1537 	    * we use double-double arithmetic here
1538 	    */
1539 	   QUAD_ASSIGN(constant, SCIPexprGetEvalValue(nlexpr)); /* f(sol) */
1540 	   for( i = 0; i < nlhdlrexprdata->nleafs; ++i )
1541 	   {
1542 	      SCIP_VAR* var;
1543 	      SCIP_Real deriv;
1544 	      SCIP_Real varval;
1545 	
1546 	      assert(SCIPexprGetDiffTag(nlhdlrexprdata->leafexprs[i]) == SCIPexprGetDiffTag(nlexpr));
1547 	      deriv = SCIPexprGetDerivative(nlhdlrexprdata->leafexprs[i]);
1548 	      if( deriv == SCIP_INVALID )
1549 	      {
1550 	         SCIPdebugMsg(scip, "gradient evaluation error for component %d of %p\n", i, (void*)nlexpr);
1551 	         return SCIP_OKAY;
1552 	      }
1553 	
1554 	      var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]);
1555 	      assert(var != NULL);
1556 	
1557 	      varval = SCIPgetSolVal(scip, sol, var);
1558 	
1559 	      SCIPdebugMsg(scip, "add %g * (<%s> - %g) to rowprep\n", deriv, SCIPvarGetName(var), varval);
1560 	
1561 	      /* add deriv * var to rowprep and deriv * (-varval) to constant */
1562 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, deriv) );
1563 	      SCIPquadprecSumQD(constant, constant, -deriv * varval);
1564 	   }
1565 	
1566 	   SCIProwprepAddConstant(rowprep, QUAD_TO_DBL(constant));
1567 	   SCIProwprepSetLocal(rowprep, FALSE);
1568 	
1569 	   *success = TRUE;
1570 	
1571 	   return SCIP_OKAY;
1572 	}
1573 	
1574 	/** adds an estimator computed via a gradient to a given rowprep, possibly perturbing solution */
1575 	static
1576 	SCIP_RETCODE estimateGradient(
1577 	   SCIP*                 scip,               /**< SCIP data structure */
1578 	   SCIP_NLHDLR*          nlhdlr,             /**< nonlinear handler */
1579 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
1580 	   SCIP_SOL*             sol,                /**< solution to use */
1581 	   SCIP_ROWPREP*         rowprep,            /**< rowprep where to store estimator */
1582 	   SCIP_Bool*            success             /**< buffer to store whether successful */
1583 	   )
1584 	{
1585 	   SCIP_NLHDLRDATA* nlhdlrdata;
1586 	   int i;
1587 	
1588 	   assert(nlhdlrexprdata != NULL);
1589 	   assert(rowprep != NULL);
1590 	   assert(success != NULL);
1591 	
1592 	#ifdef SCIP_DEBUG
1593 	   SCIPinfoMessage(scip, NULL, "estimate expression ");
1594 	   SCIPprintExpr(scip, nlhdlrexprdata->nlexpr, NULL);
1595 	   SCIPinfoMessage(scip, NULL, " by gradient\n");
1596 	#endif
1597 	
1598 	   *success = FALSE;
1599 	
1600 	   SCIP_CALL( estimateGradientInner(scip, nlhdlrexprdata, sol, rowprep, success) );
1601 	
1602 	   /* if succeeded, then there was no gradient evaluation error, so we are done */
1603 	   if( *success )
1604 	      return SCIP_OKAY;
1605 	
1606 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1607 	   assert(nlhdlrdata != NULL);
1608 	
1609 	   /* we take maxperturb == 0 as signal to not try perturbation */
1610 	   if( nlhdlrdata->maxperturb == 0.0 )
1611 	   {
1612 	      SCIPdebugMsg(scip, "gradient evaluation error, perturbation disabled\n");
1613 	      return SCIP_OKAY;
1614 	   }
1615 	
1616 	   SCIPdebugMsg(scip, "gradient evaluation error, try perturbed point\n");
1617 	
1618 	   if( nlhdlrdata->evalsol == NULL )
1619 	   {
1620 	      SCIP_CALL( SCIPcreateSol(scip, &nlhdlrdata->evalsol, NULL) );
1621 	   }
1622 	
1623 	   if( nlhdlrdata->randnumgen == NULL )
1624 	   {
1625 	      SCIP_CALL( SCIPcreateRandom(scip, &nlhdlrdata->randnumgen, RANDNUMINITSEED, TRUE) );
1626 	   }
1627 	
1628 	   for( i = 0; i < nlhdlrexprdata->nleafs; ++i )
1629 	   {
1630 	      SCIP_VAR* var;
1631 	      SCIP_Real lb;
1632 	      SCIP_Real ub;
1633 	      SCIP_Real val;
1634 	      SCIP_Real p;
1635 	
1636 	      var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]);
1637 	      assert(var != NULL);
1638 	
1639 	      lb = SCIPvarGetLbGlobal(var);
1640 	      ub = SCIPvarGetUbGlobal(var);
1641 	      val = SCIPgetSolVal(scip, sol, var);
1642 	      val = MIN(ub, MAX(lb, val));
1643 	
1644 	      p = SCIPrandomGetReal(nlhdlrdata->randnumgen, -nlhdlrdata->maxperturb, nlhdlrdata->maxperturb);
1645 	      if( !SCIPisZero(scip, val) )
1646 	         p *= REALABS(val);
1647 	
1648 	      /* if perturbation to left underruns lower bound, then try perturb to right
1649 	       * if perturbation to right exceeds lower bound, then try perturb to left
1650 	       */
1651 	      if( val + p <= lb )
1652 	         val -= p;
1653 	      else if( val + p >= ub )
1654 	         val -= p;
1655 	      else
1656 	         val += p;
1657 	
1658 	      /* if still exceeding bound, then |ub-lb| < 2*maxperturb and we pick a random value between bounds
1659 	       * (if ub-lb < 2eps, then SCIPrandomGetReal() still gives a reasonable number in [lb,ub])
1660 	       */
1661 	      if( val <= lb || val >= ub )
1662 	         val = SCIPrandomGetReal(nlhdlrdata->randnumgen, lb + SCIPepsilon(scip), ub - SCIPepsilon(scip));
1663 	
1664 	      SCIP_CALL( SCIPsetSolVal(scip, nlhdlrdata->evalsol, var, val) );
1665 	   }
1666 	
1667 	   /* try again with perturbed point */
1668 	   SCIP_CALL( estimateGradientInner(scip, nlhdlrexprdata, nlhdlrdata->evalsol, rowprep, success) );
1669 	
1670 	   return SCIP_OKAY;
1671 	}
1672 	
1673 	/** adds an estimator generated by putting a secant through the coordinates given by the two closest integer points */
1674 	static
1675 	SCIP_RETCODE estimateConvexSecant(
1676 	   SCIP*                 scip,               /**< SCIP data structure */
1677 	   SCIP_NLHDLR*          nlhdlr,             /**< nonlinear handler */
1678 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
1679 	   SCIP_SOL*             sol,                /**< solution to use, unless usemidpoint is TRUE */
1680 	   SCIP_ROWPREP*         rowprep,            /**< rowprep where to store estimator */
1681 	   SCIP_Bool*            success             /**< buffer to store whether successful */
1682 	   )
1683 	{
1684 	   SCIP_NLHDLRDATA* nlhdlrdata;
1685 	   SCIP_EXPR* nlexpr;
1686 	   SCIP_VAR* var;
1687 	   SCIP_Real x;
1688 	   SCIP_Real left, right;
1689 	   SCIP_Real fleft, fright;
1690 	
1691 	   assert(nlhdlrexprdata != NULL);
1692 	   assert(nlhdlrexprdata->nleafs == 1);
1693 	   assert(rowprep != NULL);
1694 	   assert(success != NULL);
1695 	
1696 	   nlexpr = nlhdlrexprdata->nlexpr;
1697 	   assert(nlexpr != NULL);
1698 	
1699 	   *success = FALSE;
1700 	
1701 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1702 	   assert(nlhdlrdata != NULL);
1703 	
1704 	   var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[0]);
1705 	   assert(var != NULL);
1706 	
1707 	   x = SCIPgetSolVal(scip, sol, var);
1708 	
1709 	#ifdef SCIP_DEBUG
1710 	   SCIPinfoMessage(scip, NULL, "estimate expression ");
1711 	   SCIPprintExpr(scip, nlexpr, NULL);
1712 	   SCIPinfoMessage(scip, NULL, " by secant\n");
1713 	   SCIPinfoMessage(scip, NULL, "integral variable <%s> = %g [%g,%g]\n", SCIPvarGetName(var),
1714 	         x, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
1715 	#endif
1716 	
1717 	   /* find out coordinates of var left and right to sol */
1718 	   if( SCIPisIntegral(scip, x) )
1719 	   {
1720 	      x = SCIPround(scip, x);
1721 	      if( SCIPisEQ(scip, x, SCIPvarGetLbGlobal(var)) )
1722 	      {
1723 	         left = x;
1724 	         right = left + 1.0;
1725 	      }
1726 	      else
1727 	      {
1728 	         right = x;
1729 	         left = right - 1.0;
1730 	      }
1731 	   }
1732 	   else
1733 	   {
1734 	      left = SCIPfloor(scip, x);
1735 	      right = SCIPceil(scip, x);
1736 	   }
1737 	   assert(left != right);
1738 	
1739 	   /* now evaluate at left and right */
1740 	   if( nlhdlrdata->evalsol == NULL )
1741 	   {
1742 	      SCIP_CALL( SCIPcreateSol(scip, &nlhdlrdata->evalsol, NULL) );
1743 	   }
1744 	
1745 	   SCIP_CALL( SCIPsetSolVal(scip, nlhdlrdata->evalsol, var, left) );
1746 	   SCIP_CALL( SCIPevalExpr(scip, nlexpr, nlhdlrdata->evalsol, 0L) );
1747 	
1748 	   /* evaluation error or a too large constant -> skip */
1749 	   fleft = SCIPexprGetEvalValue(nlexpr);
1750 	   if( SCIPisInfinity(scip, REALABS(fleft)) )
1751 	   {
1752 	      SCIPdebugMsg(scip, "evaluation error / too large value (%g) for %p\n", SCIPexprGetEvalValue(nlexpr), (void*)nlexpr);
1753 	      return SCIP_OKAY;
1754 	   }
1755 	
1756 	   SCIP_CALL( SCIPsetSolVal(scip, nlhdlrdata->evalsol, var, right) );
1757 	   SCIP_CALL( SCIPevalExpr(scip, nlexpr, nlhdlrdata->evalsol, 0L) );
1758 	
1759 	   /* evaluation error or a too large constant -> skip */
1760 	   fright = SCIPexprGetEvalValue(nlexpr);
1761 	   if( SCIPisInfinity(scip, REALABS(fright)) )
1762 	   {
1763 	      SCIPdebugMsg(scip, "evaluation error / too large value (%g) for %p\n", SCIPexprGetEvalValue(nlexpr), (void*)nlexpr);
1764 	      return SCIP_OKAY;
1765 	   }
1766 	
1767 	   SCIPdebugMsg(scip, "f(%g)=%g, f(%g)=%g\n", left, fleft, right, fright);
1768 	
1769 	   /* skip if too steep
1770 	    * for clay0204h, this resulted in a wrong cut from f(0)=1e12 f(1)=0.99998,
1771 	    * since due to limited precision, this was handled as if f(1)=1
1772 	    */
1773 	   if( (!SCIPisZero(scip, fleft)  && REALABS(fright/fleft)*SCIPepsilon(scip) > 1.0) ||
1774 	       (!SCIPisZero(scip, fright) && REALABS(fleft/fright)*SCIPepsilon(scip) > 1.0) )
1775 	   {
1776 	      SCIPdebugMsg(scip, "function is too steep, abandoning\n");
1777 	      return SCIP_OKAY;
1778 	   }
1779 	
1780 	   /* now add f(left) + (f(right) - f(left)) * (x - left) as estimator to rowprep */
1781 	   SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, fright - fleft) );
1782 	   SCIProwprepAddConstant(rowprep, fleft - (fright - fleft) * left);
1783 	   SCIProwprepSetLocal(rowprep, FALSE);
1784 	
1785 	   *success = TRUE;
1786 	
1787 	   return SCIP_OKAY;
1788 	}
1789 	
1790 	/*
1791 	 * Callback methods of convex nonlinear handler
1792 	 */
1793 	
1794 	/** free handler data of convex or concave nlhdlr */
1795 	static
1796 	SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrfreeHdlrDataConvexConcave)
1797 	{  /*lint --e{715}*/
1798 	   assert(scip != NULL);
1799 	   assert(nlhdlrdata != NULL);
1800 	   assert(*nlhdlrdata != NULL);
1801 	   assert((*nlhdlrdata)->evalsol == NULL);
1802 	   assert((*nlhdlrdata)->randnumgen == NULL);
1803 	
1804 	   SCIPfreeBlockMemory(scip, nlhdlrdata);
1805 	
1806 	   return SCIP_OKAY;
1807 	}
1808 	
1809 	/** callback to free expression specific data */
1810 	static
1811 	SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrfreeExprDataConvexConcave)
1812 	{  /*lint --e{715}*/
1813 	   assert(scip != NULL);
1814 	   assert(nlhdlrexprdata != NULL);
1815 	   assert(*nlhdlrexprdata != NULL);
1816 	
1817 	   SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrexprdata)->leafexprs, (*nlhdlrexprdata)->nleafs);
1818 	   SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->nlexpr) );
1819 	   SCIPhashmapFree(&(*nlhdlrexprdata)->nlexpr2origexpr);
1820 	
1821 	   SCIPfreeBlockMemory(scip, nlhdlrexprdata);
1822 	
1823 	   return SCIP_OKAY;
1824 	}
1825 	
1826 	/** deinitialization of problem-specific data */
1827 	static
1828 	SCIP_DECL_NLHDLREXIT(nlhdlrExitConvex)
1829 	{
1830 	   SCIP_NLHDLRDATA* nlhdlrdata;
1831 	
1832 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1833 	   assert(nlhdlrdata != NULL);
1834 	
1835 	   if( nlhdlrdata->evalsol != NULL )
1836 	   {
1837 	      SCIP_CALL( SCIPfreeSol(scip, &nlhdlrdata->evalsol) );
1838 	   }
1839 	
1840 	   if( nlhdlrdata->randnumgen != NULL )
1841 	      SCIPfreeRandom(scip, &nlhdlrdata->randnumgen);
1842 	
1843 	   return SCIP_OKAY;
1844 	}
1845 	
1846 	/** checks whether expression (or -expression) is convex, possibly after introducing auxiliary variables */
1847 	static
1848 	SCIP_DECL_NLHDLRDETECT(nlhdlrDetectConvex)
1849 	{ /*lint --e{715}*/
1850 	   SCIP_NLHDLRDATA* nlhdlrdata;
1851 	   SCIP_EXPR* nlexpr = NULL;
1852 	   SCIP_HASHMAP* nlexpr2origexpr;
1853 	   int nleafs = 0;
1854 	
1855 	   assert(scip != NULL);
1856 	   assert(nlhdlr != NULL);
1857 	   assert(expr != NULL);
1858 	   assert(enforcing != NULL);
1859 	   assert(participating != NULL);
1860 	   assert(nlhdlrexprdata != NULL);
1861 	
1862 	   /* we currently do not participate if only activity computation is required */
1863 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH )
1864 	      return SCIP_OKAY;
1865 	
1866 	   /* ignore pure constants and variables */
1867 	   if( SCIPexprGetNChildren(expr) == 0 )
1868 	      return SCIP_OKAY;
1869 	
1870 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1871 	   assert(nlhdlrdata != NULL);
1872 	   assert(nlhdlrdata->isnlhdlrconvex);
1873 	
1874 	   SCIPdebugMsg(scip, "nlhdlr_convex detect for expr %p\n", (void*)expr);
1875 	
1876 	   /* initialize mapping from copied expression to original one
1877 	    * 20 is not a bad estimate for the size of convex subexpressions that we can usually discover
1878 	    * when expressions will be allowed to store "user"data, we could get rid of this hashmap (TODO)
1879 	    */
1880 	   SCIP_CALL( SCIPhashmapCreate(&nlexpr2origexpr, SCIPblkmem(scip), 20) );
1881 	
1882 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABELOW) == 0 )  /* if no separation below yet */
1883 	   {
1884 	      SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr,
1885 	         SCIP_EXPRCURV_CONVEX, NULL, SCIPassumeConvexNonlinear(conshdlr), NULL) );
1886 	      if( nlexpr != NULL )
1887 	      {
1888 	         assert(SCIPexprGetNChildren(nlexpr) > 0);  /* should not be trivial */
1889 	
1890 	         *participating |= SCIP_NLHDLR_METHOD_SEPABELOW;
1891 	
1892 	         SCIPdebugMsg(scip, "detected expr %p to be convex -> can enforce expr <= auxvar\n", (void*)expr);
1893 	      }
1894 	      else
1895 	      {
1896 	         SCIP_CALL( SCIPhashmapRemoveAll(nlexpr2origexpr) );
1897 	      }
1898 	   }
1899 	
1900 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPAABOVE) == 0 && nlexpr == NULL )  /* if no separation above and not convex */
1901 	   {
1902 	      SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr,
1903 	         SCIP_EXPRCURV_CONCAVE, NULL, SCIPassumeConvexNonlinear(conshdlr), NULL) );
1904 	      if( nlexpr != NULL )
1905 	      {
1906 	         assert(SCIPexprGetNChildren(nlexpr) > 0);  /* should not be trivial */
1907 	
1908 	         *participating |= SCIP_NLHDLR_METHOD_SEPAABOVE;
1909 	
1910 	         SCIPdebugMsg(scip, "detected expr %p to be concave -> can enforce expr >= auxvar\n", (void*)expr);
1911 	      }
1912 	   }
1913 	
1914 	   /* everything we participate in we also enforce */
1915 	   *enforcing |= *participating;
1916 	
1917 	   assert(*participating || nlexpr == NULL);
1918 	   if( !*participating )
1919 	   {
1920 	      SCIPhashmapFree(&nlexpr2origexpr);
1921 	      return SCIP_OKAY;
1922 	   }
1923 	
1924 	   /* create the expression data of the nonlinear handler
1925 	    * notify conshdlr about expr for which we will require auxiliary variables
1926 	    */
1927 	   SCIP_CALL( createNlhdlrExprData(scip, nlhdlrdata, nlhdlrexprdata, expr, nlexpr, nlexpr2origexpr, nleafs, *participating) );
1928 	
1929 	   return SCIP_OKAY;
1930 	}
1931 	
1932 	/** auxiliary evaluation callback */
1933 	static
1934 	SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalAuxConvexConcave)
1935 	{ /*lint --e{715}*/
1936 	   assert(nlhdlrexprdata != NULL);
1937 	   assert(nlhdlrexprdata->nlexpr != NULL);
1938 	   assert(auxvalue != NULL);
1939 	
1940 	   SCIP_CALL( SCIPevalExpr(scip, nlhdlrexprdata->nlexpr, sol, 0L) );
1941 	   *auxvalue = SCIPexprGetEvalValue(nlhdlrexprdata->nlexpr);
1942 	
1943 	   return SCIP_OKAY;
1944 	}
1945 	
1946 	/** init sepa callback that initializes LP */
1947 	static
1948 	SCIP_DECL_NLHDLRINITSEPA(nlhdlrInitSepaConvex)
1949 	{  /*lint --e{715}*/
1950 	   SCIP_EXPR* nlexpr;
1951 	   SCIP_EXPRCURV curvature;
1952 	   SCIP_Bool success;
1953 	   SCIP_ROWPREP* rowprep = NULL;
1954 	   SCIP_ROW* row;
1955 	   SCIP_Real lb;
1956 	   SCIP_Real ub;
1957 	   SCIP_Real lambda;
1958 	   SCIP_SOL* sol;
1959 	   int k;
1960 	
1961 	   assert(scip != NULL);
1962 	   assert(expr != NULL);
1963 	   assert(nlhdlrexprdata != NULL);
1964 	
1965 	   /* setup nlhdlrexprdata->leafexprs */
1966 	   SCIP_CALL( collectLeafs(scip, nlhdlrexprdata) );
1967 	
1968 	   nlexpr = nlhdlrexprdata->nlexpr;
1969 	   assert(nlexpr != NULL);
1970 	   assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlexpr) == expr);
1971 	
1972 	   curvature = SCIPexprGetCurvature(nlexpr);
1973 	   assert(curvature == SCIP_EXPRCURV_CONVEX || curvature == SCIP_EXPRCURV_CONCAVE);
1974 	
1975 	   /* we can only be estimating on the convex side */
1976 	   if( curvature == SCIP_EXPRCURV_CONVEX )
1977 	      overestimate = FALSE;
1978 	   else if( curvature == SCIP_EXPRCURV_CONCAVE )
1979 	      underestimate = FALSE;
1980 	   if( !overestimate && !underestimate )
1981 	      return SCIP_OKAY;
1982 	
1983 	   /* linearizes at 5 different points obtained as convex combination of the lower and upper bound of the variables
1984 	    * present in the convex expression; whether more weight is given to the lower or upper bound of a variable depends
1985 	    * on whether the fixing of the variable to that value is better for the objective function
1986 	    */
1987 	   SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) );
1988 	
1989 	   *infeasible = FALSE;
1990 	
1991 	   for( k = 0; k < 5; ++k )
1992 	   {
1993 	      int i;
1994 	      lambda = 0.1 * (k+1); /* lambda = 0.1, 0.2, 0.3, 0.4, 0.5 */
1995 	
1996 	      for( i = 0; i < nlhdlrexprdata->nleafs; ++i )
1997 	      {
1998 	         SCIP_VAR* var;
1999 	
2000 	         var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]);
2001 	
2002 	         lb = SCIPvarGetLbGlobal(var);
2003 	         ub = SCIPvarGetUbGlobal(var);
2004 	
2005 	         /* make sure the absolute values of bounds are not too large */
2006 	         if( ub > -INITLPMAXVARVAL )
2007 	            lb = MAX(lb, -INITLPMAXVARVAL);
2008 	         if( lb <  INITLPMAXVARVAL )
2009 	            ub = MIN(ub,  INITLPMAXVARVAL);
2010 	
2011 	         /* in the case when ub < -maxabsbnd or lb > maxabsbnd, we still want to at least make bounds finite */
2012 	         if( SCIPisInfinity(scip, -lb) )
2013 	            lb = MIN(-10.0, ub - 0.1*REALABS(ub));
2014 	         if( SCIPisInfinity(scip,  ub) )
2015 	            ub = MAX( 10.0, lb + 0.1*REALABS(lb));
2016 	
2017 	         if( SCIPvarGetBestBoundType(var) == SCIP_BOUNDTYPE_LOWER )
2018 	            SCIP_CALL( SCIPsetSolVal(scip, sol, var, lambda * ub + (1.0 - lambda) * lb) );
2019 	         else
2020 	            SCIP_CALL( SCIPsetSolVal(scip, sol, var, lambda * lb + (1.0 - lambda) * ub) );
2021 	      }
2022 	
2023 	      SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
2024 	      SCIP_CALL( estimateGradient(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, &success) );
2025 	      if( !success )
2026 	      {
2027 	         SCIPdebugMsg(scip, "failed to linearize for k = %d\n", k);
2028 	         SCIPfreeRowprep(scip, &rowprep);
2029 	         continue;
2030 	      }
2031 	
2032 	      /* add auxiliary variable */
2033 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) );
2034 	
2035 	      /* straighten out numerics */
2036 	      SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) );
2037 	      if( !success )
2038 	      {
2039 	         SCIPdebugMsg(scip, "failed to cleanup rowprep numerics for k = %d\n", k);
2040 	         SCIPfreeRowprep(scip, &rowprep);
2041 	         continue;
2042 	      }
2043 	
2044 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_gradient%p_initsepa_%d",
2045 	            overestimate ? "over" : "under", (void*)expr, k);
2046 	      SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );
2047 	      SCIPfreeRowprep(scip, &rowprep);
2048 	
2049 	#ifdef SCIP_DEBUG
2050 	      SCIPinfoMessage(scip, NULL, "initsepa computed row: ");
2051 	      SCIPprintRow(scip, row, NULL);
2052 	#endif
2053 	
2054 	      SCIP_CALL( SCIPaddRow(scip, row, FALSE, infeasible) );
2055 	      SCIP_CALL( SCIPreleaseRow(scip, &row) );
2056 	
2057 	      if( *infeasible )
2058 	         break;
2059 	   }
2060 	
2061 	   SCIP_CALL( SCIPfreeSol(scip, &sol) );
2062 	
2063 	   return SCIP_OKAY;
2064 	}
2065 	
2066 	/** estimator callback */
2067 	static
2068 	SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateConvex)
2069 	{ /*lint --e{715}*/
2070 	   SCIP_ROWPREP* rowprep;
2071 	
2072 	   assert(scip != NULL);
2073 	   assert(expr != NULL);
2074 	   assert(nlhdlrexprdata != NULL);
2075 	   assert(nlhdlrexprdata->nlexpr != NULL);
2076 	   assert(rowpreps != NULL);
2077 	   assert(success != NULL);
2078 	
2079 	   assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlhdlrexprdata->nlexpr) == expr);
2080 	
2081 	   /* we must be called only for the side that we indicated to participate in during DETECT */
2082 	   assert(SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX
2083 	          || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE);
2084 	   assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE);
2085 	   assert( overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX);
2086 	
2087 	   *success = FALSE;
2088 	   *addedbranchscores = FALSE;
2089 	
2090 	   /* we can skip eval as nlhdlrEvalAux should have been called for same solution before */
2091 	   /* SCIP_CALL( nlhdlrExprEval(scip, nlexpr, sol) ); */
2092 	   assert(auxvalue == SCIPexprGetEvalValue(nlhdlrexprdata->nlexpr)); /* given value (originally from
2093 	         nlhdlrEvalAuxConvexConcave) should coincide with the one stored in nlexpr */
2094 	
2095 	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
2096 	
2097 	   if( nlhdlrexprdata->nleafs == 1 && SCIPexprIsIntegral(nlhdlrexprdata->leafexprs[0]) )
2098 	   {
2099 	      SCIP_CALL( estimateConvexSecant(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, success) );
2100 	
2101 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexsecant%p_%s%" SCIP_LONGINT_FORMAT,
2102 	         overestimate ? "over" : "under",
2103 	         (void*)expr,
2104 	         sol != NULL ? "sol" : "lp",
2105 	         sol != NULL ? (SCIP_Longint) SCIPsolGetIndex(sol) : SCIPgetNLPs(scip));
2106 	   }
2107 	
2108 	   /* if secant method was not used or failed, then try with gradient (unless we had an evaluation error in sol before) */
2109 	   if( !*success && auxvalue != SCIP_INVALID )
2110 	   {
2111 	      SCIP_CALL( estimateGradient(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, success) );
2112 	
2113 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexgradient%p_%s%" SCIP_LONGINT_FORMAT,
2114 	         overestimate ? "over" : "under",
2115 	         (void*)expr,
2116 	         sol != NULL ? "sol" : "lp",
2117 	         sol != NULL ? (SCIP_Longint) SCIPsolGetIndex(sol) : SCIPgetNLPs(scip));
2118 	   }
2119 	
2120 	   if( *success )
2121 	   {
2122 	      SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
2123 	   }
2124 	   else
2125 	   {
2126 	      SCIPfreeRowprep(scip, &rowprep);
2127 	   }
2128 	
2129 	   return SCIP_OKAY;
2130 	}
2131 	
2132 	/** solution notification callback */
2133 	static
2134 	SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeConvex)
2135 	{ /*lint --e{715}*/
2136 	   SCIP_ROWPREP* rowprep;
2137 	   SCIP_Bool success = FALSE;
2138 	
2139 	   assert(scip != NULL);
2140 	   assert(expr != NULL);
2141 	   assert(nlhdlrexprdata != NULL);
2142 	   assert(nlhdlrexprdata->nlexpr != NULL);
2143 	   assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlhdlrexprdata->nlexpr) == expr);
2144 	   assert(sol != NULL);
2145 	
2146 	   /* we must be called only for the side that we indicated to participate in during DETECT */
2147 	   assert(SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX
2148 	          || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE);
2149 	   assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE);
2150 	   assert(!underestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX);
2151 	   assert(underestimate == !overestimate);  /* should be exactly one of under- and overestimate */
2152 	
2153 	   /* evaluate nlexpr in solution */
2154 	   SCIP_CALL( SCIPevalExpr(scip, nlhdlrexprdata->nlexpr, sol, 0L) );
2155 	
2156 	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
2157 	
2158 	   if( nlhdlrexprdata->nleafs == 1 && SCIPexprIsIntegral(nlhdlrexprdata->leafexprs[0]) )
2159 	   {
2160 	      SCIP_CALL( estimateConvexSecant(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, &success) );
2161 	
2162 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexsecant%p_sol%dnotify",
2163 	         overestimate ? "over" : "under", (void*)expr, SCIPsolGetIndex(sol));
2164 	   }
2165 	
2166 	   /* if secant method was not used or failed, then try with gradient */
2167 	   if( !success )
2168 	   {
2169 	      SCIP_CALL( estimateGradient(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, &success) );
2170 	
2171 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexgradient%p_sol%dnotify",
2172 	         overestimate ? "over" : "under", (void*)expr, SCIPsolGetIndex(sol));
2173 	   }
2174 	
2175 	   if( success )
2176 	   {
2177 	      /* complete estimator to cut and clean it up */
2178 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) );
2179 	      SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) );
2180 	   }
2181 	
2182 	   /* if cleanup succeeded and rowprep is still global, add to cutpool */
2183 	   if( success && !SCIProwprepIsLocal(rowprep) )
2184 	   {
2185 	      SCIP_ROW* row;
2186 	
2187 	      SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );
2188 	      SCIP_CALL( SCIPaddPoolCut(scip, row) );
2189 	      SCIP_CALL( SCIPreleaseRow(scip, &row) );
2190 	   }
2191 	
2192 	   SCIPfreeRowprep(scip, &rowprep);
2193 	
2194 	   return SCIP_OKAY;
2195 	}
2196 	
2197 	/** include nlhdlr in another scip instance */
2198 	static
2199 	SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrConvex)
2200 	{ /*lint --e{715}*/
2201 	   assert(targetscip != NULL);
2202 	   assert(sourcenlhdlr != NULL);
2203 	   assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), CONVEX_NLHDLR_NAME) == 0);
2204 	
2205 	   SCIP_CALL( SCIPincludeNlhdlrConvex(targetscip) );
2206 	
2207 	   return SCIP_OKAY;
2208 	}
2209 	
2210 	/** includes convex nonlinear handler in nonlinear constraint handler */
2211 	SCIP_RETCODE SCIPincludeNlhdlrConvex(
2212 	   SCIP*                 scip                /**< SCIP data structure */
2213 	   )
2214 	{
2215 	   SCIP_NLHDLR* nlhdlr;
2216 	   SCIP_NLHDLRDATA* nlhdlrdata;
2217 	
2218 	   assert(scip != NULL);
2219 	
2220 	   SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) );
2221 	   nlhdlrdata->isnlhdlrconvex = TRUE;
2222 	   nlhdlrdata->evalsol = NULL;
2223 	   nlhdlrdata->randnumgen = NULL;
2224 	
2225 	   SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, CONVEX_NLHDLR_NAME, CONVEX_NLHDLR_DESC,
2226 	      CONVEX_NLHDLR_DETECTPRIORITY, CONVEX_NLHDLR_ENFOPRIORITY, nlhdlrDetectConvex, nlhdlrEvalAuxConvexConcave, nlhdlrdata) );
2227 	   assert(nlhdlr != NULL);
2228 	
2229 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/detectsum",
2230 	      "whether to run convexity detection when the root of an expression is a non-quadratic sum",
2231 	      &nlhdlrdata->detectsum, FALSE, DEFAULT_DETECTSUM, NULL, NULL) );
2232 	
2233 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/extendedform",
2234 	      "whether to create extended formulations instead of looking for maximal convex expressions",
2235 	      &nlhdlrdata->extendedform, FALSE, DEFAULT_EXTENDEDFORM, NULL, NULL) );
2236 	
2237 	   SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/maxperturb",
2238 	      "maximal relative perturbation of non-differentiable reference point",
2239 	      &nlhdlrdata->maxperturb, FALSE, DEFAULT_MAXPERTURB, 0.0, 1.0, NULL, NULL) );
2240 	
2241 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/cvxquadratic",
2242 	      "whether to use convexity check on quadratics",
2243 	      &nlhdlrdata->cvxquadratic, TRUE, DEFAULT_CVXQUADRATIC_CONVEX, NULL, NULL) );
2244 	
2245 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/cvxsignomial",
2246 	      "whether to use convexity check on signomials",
2247 	      &nlhdlrdata->cvxsignomial, TRUE, DEFAULT_CVXSIGNOMIAL, NULL, NULL) );
2248 	
2249 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/cvxprodcomp",
2250 	      "whether to use convexity check on product composition f(h)*h",
2251 	      &nlhdlrdata->cvxprodcomp, TRUE, DEFAULT_CVXPRODCOMP, NULL, NULL) );
2252 	
2253 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/handletrivial",
2254 	      "whether to also handle trivial convex expressions",
2255 	      &nlhdlrdata->handletrivial, TRUE, DEFAULT_HANDLETRIVIAL, NULL, NULL) );
2256 	
2257 	   SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrfreeHdlrDataConvexConcave);
2258 	   SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrConvex);
2259 	   SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrfreeExprDataConvexConcave);
2260 	   SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaConvex, NULL, nlhdlrEstimateConvex, NULL);
2261 	   SCIPnlhdlrSetSollinearize(nlhdlr, nlhdlrSollinearizeConvex);
2262 	   SCIPnlhdlrSetInitExit(nlhdlr, NULL, nlhdlrExitConvex);
2263 	
2264 	   return SCIP_OKAY;
2265 	}
2266 	
2267 	/*
2268 	 * Callback methods of concave nonlinear handler
2269 	 */
2270 	
2271 	/** deinitialization of problem-specific data */
2272 	static
2273 	SCIP_DECL_NLHDLREXIT(nlhdlrExitConcave)
2274 	{
2275 	   SCIP_NLHDLRDATA* nlhdlrdata;
2276 	
2277 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
2278 	   assert(nlhdlrdata != NULL);
2279 	   assert(nlhdlrdata->randnumgen == NULL); /* not used for concave nlhdlr so far */
2280 	
2281 	   if( nlhdlrdata->evalsol != NULL )
2282 	   {
2283 	      SCIP_CALL( SCIPfreeSol(scip, &nlhdlrdata->evalsol) );
2284 	   }
2285 	
2286 	   return SCIP_OKAY;
2287 	}
2288 	
2289 	/** checks whether expression (or -expression) is concave, possibly after introducing auxiliary variables */
2290 	static
2291 	SCIP_DECL_NLHDLRDETECT(nlhdlrDetectConcave)
2292 	{ /*lint --e{715}*/
2293 	   SCIP_NLHDLRDATA* nlhdlrdata;
2294 	   SCIP_EXPR* nlexpr = NULL;
2295 	   SCIP_HASHMAP* nlexpr2origexpr;
2296 	   int nleafs = 0;
2297 	
2298 	   assert(scip != NULL);
2299 	   assert(nlhdlr != NULL);
2300 	   assert(expr != NULL);
2301 	   assert(enforcing != NULL);
2302 	   assert(participating != NULL);
2303 	   assert(nlhdlrexprdata != NULL);
2304 	
2305 	   /* we currently do not participate if only activity computation is required */
2306 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH )
2307 	      return SCIP_OKAY;
2308 	
2309 	   /* ignore pure constants and variables */
2310 	   if( SCIPexprGetNChildren(expr) == 0 )
2311 	      return SCIP_OKAY;
2312 	
2313 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
2314 	   assert(nlhdlrdata != NULL);
2315 	   assert(!nlhdlrdata->isnlhdlrconvex);
2316 	
2317 	   SCIPdebugMsg(scip, "nlhdlr_concave detect for expr %p\n", (void*)expr);
2318 	
2319 	   /* initialize mapping from copied expression to original one
2320 	    * 20 is not a bad estimate for the size of concave subexpressions that we can usually discover
2321 	    * when expressions will be allowed to store "user"data, we could get rid of this hashmap (TODO)
2322 	    */
2323 	   SCIP_CALL( SCIPhashmapCreate(&nlexpr2origexpr, SCIPblkmem(scip), 20) );
2324 	
2325 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABELOW) == 0 )  /* if no separation below yet */
2326 	   {
2327 	      SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr,
2328 	         SCIP_EXPRCURV_CONCAVE, NULL, FALSE, NULL) );
2329 	
2330 	      if( nlexpr != NULL && nleafs > SCIP_MAXVERTEXPOLYDIM )
2331 	      {
2332 	         SCIPdebugMsg(scip, "Too many variables (%d) in constructed expression. Will not be able to estimate. Rejecting.\n", nleafs);
2333 	         SCIP_CALL( SCIPreleaseExpr(scip, &nlexpr) );
2334 	      }
2335 	
2336 	      if( nlexpr != NULL )
2337 	      {
2338 	         assert(SCIPexprGetNChildren(nlexpr) > 0);  /* should not be trivial */
2339 	
2340 	         *participating |= SCIP_NLHDLR_METHOD_SEPABELOW;
2341 	
2342 	         SCIPdebugMsg(scip, "detected expr %p to be concave -> can enforce expr <= auxvar\n", (void*)expr);
2343 	      }
2344 	      else
2345 	      {
2346 	         SCIP_CALL( SCIPhashmapRemoveAll(nlexpr2origexpr) );
2347 	      }
2348 	   }
2349 	
2350 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPAABOVE) == 0 && nlexpr == NULL )  /* if no separation above and not concave */
2351 	   {
2352 	      SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr,
2353 	         SCIP_EXPRCURV_CONVEX, NULL, FALSE, NULL) );
2354 	
2355 	      if( nlexpr != NULL && nleafs > SCIP_MAXVERTEXPOLYDIM )
2356 	      {
2357 	         SCIPdebugMsg(scip, "Too many variables (%d) in constructed expression. Will not be able to estimate. Rejecting.\n", nleafs);
2358 	         SCIP_CALL( SCIPreleaseExpr(scip, &nlexpr) );
2359 	      }
2360 	
2361 	      if( nlexpr != NULL )
2362 	      {
2363 	         assert(SCIPexprGetNChildren(nlexpr) > 0);  /* should not be trivial */
2364 	
2365 	         *participating |= SCIP_NLHDLR_METHOD_SEPAABOVE;
2366 	
2367 	         SCIPdebugMsg(scip, "detected expr %p to be convex -> can enforce expr >= auxvar\n", (void*)expr);
2368 	      }
2369 	   }
2370 	
2371 	   /* everything we participate in we also enforce (at the moment) */
2372 	   *enforcing |= *participating;
2373 	
2374 	   assert(*participating || nlexpr == NULL);
2375 	   if( !*participating )
2376 	   {
2377 	      SCIPhashmapFree(&nlexpr2origexpr);
2378 	      return SCIP_OKAY;
2379 	   }
2380 	
2381 	   /* create the expression data of the nonlinear handler
2382 	    * notify conshdlr about expr for which we will require auxiliary variables and use activity
2383 	    */
2384 	   SCIP_CALL( createNlhdlrExprData(scip, nlhdlrdata, nlhdlrexprdata, expr, nlexpr, nlexpr2origexpr, nleafs, *participating) );
2385 	
2386 	   return SCIP_OKAY;
2387 	}
2388 	
2389 	/** init sepa callback that initializes LP */
2390 	static
2391 	SCIP_DECL_NLHDLRINITSEPA(nlhdlrInitSepaConcave)
2392 	{
2393 	   SCIP_EXPR* nlexpr;
2394 	   SCIP_EXPRCURV curvature;
2395 	   SCIP_Bool success;
2396 	   SCIP_ROWPREP* rowprep = NULL;
2397 	   SCIP_ROW* row;
2398 	
2399 	   assert(scip != NULL);
2400 	   assert(expr != NULL);
2401 	   assert(nlhdlrexprdata != NULL);
2402 	
2403 	   nlexpr = nlhdlrexprdata->nlexpr;
2404 	   assert(nlexpr != NULL);
2405 	   assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlexpr) == expr);
2406 	
2407 	   /* setup nlhdlrexprdata->leafexprs */
2408 	   SCIP_CALL( collectLeafs(scip, nlhdlrexprdata) );
2409 	
2410 	   curvature = SCIPexprGetCurvature(nlexpr);
2411 	   assert(curvature == SCIP_EXPRCURV_CONVEX || curvature == SCIP_EXPRCURV_CONCAVE);
2412 	   /* we can only be estimating on non-convex side */
2413 	   if( curvature == SCIP_EXPRCURV_CONCAVE )
2414 	      overestimate = FALSE;
2415 	   else if( curvature == SCIP_EXPRCURV_CONVEX )
2416 	      underestimate = FALSE;
2417 	   if( !overestimate && !underestimate )
2418 	      return SCIP_OKAY;
2419 	
2420 	   /* compute estimator and store in rowprep */
2421 	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
2422 	   SCIP_CALL( estimateVertexPolyhedral(scip, conshdlr, nlhdlr, nlhdlrexprdata, NULL, TRUE, overestimate,
2423 	         overestimate ? SCIPinfinity(scip) : -SCIPinfinity(scip), rowprep, &success) );
2424 	   if( !success )
2425 	   {
2426 	      SCIPdebugMsg(scip, "failed to compute facet of convex hull\n");
2427 	      goto TERMINATE;
2428 	   }
2429 	
2430 	   /* add auxiliary variable */
2431 	   SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) );
2432 	
2433 	   /* straighten out numerics */
2434 	   SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) );
2435 	   if( !success )
2436 	   {
2437 	      SCIPdebugMsg(scip, "failed to cleanup rowprep numerics\n");
2438 	      goto TERMINATE;
2439 	   }
2440 	
2441 	   (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_concave%p_initsepa",
2442 	         overestimate ? "over" : "under", (void*)expr);
2443 	   SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );
2444 	
2445 	#ifdef SCIP_DEBUG
2446 	   SCIPinfoMessage(scip, NULL, "initsepa computed row: ");
2447 	   SCIPprintRow(scip, row, NULL);
2448 	#endif
2449 	
2450 	   SCIP_CALL( SCIPaddRow(scip, row, FALSE, infeasible) );
2451 	   SCIP_CALL( SCIPreleaseRow(scip, &row) );
2452 	
2453 	 TERMINATE:
2454 	   if( rowprep != NULL )
2455 	      SCIPfreeRowprep(scip, &rowprep);
2456 	
2457 	   return SCIP_OKAY;
2458 	}
2459 	
2460 	/** estimator callback */
2461 	static
2462 	SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateConcave)
2463 	{ /*lint --e{715}*/
2464 	   SCIP_ROWPREP* rowprep;
2465 	
2466 	   assert(scip != NULL);
2467 	   assert(expr != NULL);
2468 	   assert(nlhdlrexprdata != NULL);
2469 	   assert(nlhdlrexprdata->nlexpr != NULL);
2470 	   assert(rowpreps != NULL);
2471 	   assert(success != NULL);
2472 	
2473 	   assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlhdlrexprdata->nlexpr) == expr);
2474 	
2475 	   /* we must be called only for the side that we indicated to participate in during DETECT */
2476 	   assert(SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX
2477 	         || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE);
2478 	   assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX);
2479 	   assert( overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE);
2480 	
2481 	   *success = FALSE;
2482 	   *addedbranchscores = FALSE;
2483 	
2484 	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
2485 	
2486 	   SCIP_CALL( estimateVertexPolyhedral(scip, conshdlr, nlhdlr, nlhdlrexprdata, sol, FALSE, overestimate, targetvalue, rowprep, success) );
2487 	
2488 	   if( *success )
2489 	   {
2490 	      SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
2491 	
2492 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_concave%p_%s%" SCIP_LONGINT_FORMAT,
2493 	         overestimate ? "over" : "under",
2494 	         (void*)expr,
2495 	         sol != NULL ? "sol" : "lp",
2496 	         sol != NULL ? (SCIP_Longint) SCIPsolGetIndex(sol) : SCIPgetNLPs(scip));
2497 	   }
2498 	   else
2499 	   {
2500 	      SCIPfreeRowprep(scip, &rowprep);
2501 	   }
2502 	
2503 	   if( addbranchscores )
2504 	   {
2505 	      SCIP_Real violation;
2506 	
2507 	      /* check how much is the violation on the side that we estimate */
2508 	      if( auxvalue == SCIP_INVALID )
2509 	      {
2510 	         /* if cannot evaluate, then always branch */
2511 	         violation = SCIPinfinity(scip);
2512 	      }
2513 	      else
2514 	      {
2515 	         SCIP_Real auxval;
2516 	
2517 	         /* get value of auxiliary variable of this expression */
2518 	         assert(SCIPgetExprAuxVarNonlinear(expr) != NULL);
2519 	         auxval = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr));
2520 	
2521 	         /* compute the violation
2522 	          * if we underestimate, then we enforce expr <= auxval, so violation is (positive part of) auxvalue - auxval
2523 	          * if we overestimate,  then we enforce expr >= auxval, so violation is (positive part of) auxval - auxvalue
2524 	          */
2525 	         if( !overestimate )
2526 	            violation = MAX(0.0, auxvalue - auxval);
2527 	         else
2528 	            violation = MAX(0.0, auxval - auxvalue);
2529 	      }
2530 	      assert(violation >= 0.0);
2531 	
2532 	      /* add violation as branching-score to expressions; the core will take care distributing this onto variables */
2533 	      if( nlhdlrexprdata->nleafs == 1 )
2534 	      {
2535 	         SCIP_EXPR* e;
2536 	         e = (SCIP_EXPR*)SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, nlhdlrexprdata->leafexprs[0]);
2537 	         SCIP_CALL( SCIPaddExprsViolScoreNonlinear(scip, &e, 1, violation, sol, addedbranchscores) );
2538 	      }
2539 	      else
2540 	      {
2541 	         SCIP_EXPR** exprs;
2542 	         int c;
2543 	
2544 	         /* map leaf expressions back to original expressions
2545 	          * TODO do this once at end of detect and store in nlhdlrexprdata
2546 	          */
2547 	         SCIP_CALL( SCIPallocBufferArray(scip, &exprs, nlhdlrexprdata->nleafs) );
2548 	         for( c = 0; c < nlhdlrexprdata->nleafs; ++c )
2549 	               exprs[c] = (SCIP_EXPR*)SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, nlhdlrexprdata->leafexprs[c]);
2550 	
2551 	         SCIP_CALL( SCIPaddExprsViolScoreNonlinear(scip, exprs, nlhdlrexprdata->nleafs, violation, sol, addedbranchscores) );
2552 	
2553 	         SCIPfreeBufferArray(scip, &exprs);
2554 	      }
2555 	   }
2556 	
2557 	   return SCIP_OKAY;
2558 	}
2559 	
2560 	/** includes nonlinear handler in another scip instance */
2561 	static
2562 	SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrConcave)
2563 	{ /*lint --e{715}*/
2564 	   assert(targetscip != NULL);
2565 	   assert(sourcenlhdlr != NULL);
2566 	   assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), CONCAVE_NLHDLR_NAME) == 0);
2567 	
2568 	   SCIP_CALL( SCIPincludeNlhdlrConcave(targetscip) );
2569 	
2570 	   return SCIP_OKAY;
2571 	}
2572 	
2573 	/** includes concave nonlinear handler in nonlinear constraint handler */
2574 	SCIP_RETCODE SCIPincludeNlhdlrConcave(
2575 	   SCIP*                 scip                /**< SCIP data structure */
2576 	   )
2577 	{
2578 	   SCIP_NLHDLR* nlhdlr;
2579 	   SCIP_NLHDLRDATA* nlhdlrdata;
2580 	
2581 	   assert(scip != NULL);
2582 	
2583 	   SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) );
2584 	   nlhdlrdata->isnlhdlrconvex = FALSE;
2585 	   nlhdlrdata->evalsol = NULL;
2586 	   nlhdlrdata->randnumgen = NULL;
2587 	
2588 	   SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, CONCAVE_NLHDLR_NAME, CONCAVE_NLHDLR_DESC,
2589 	      CONCAVE_NLHDLR_DETECTPRIORITY, CONCAVE_NLHDLR_ENFOPRIORITY, nlhdlrDetectConcave, nlhdlrEvalAuxConvexConcave, nlhdlrdata) );
2590 	   assert(nlhdlr != NULL);
2591 	
2592 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/detectsum",
2593 	      "whether to run convexity detection when the root of an expression is a sum",
2594 	      &nlhdlrdata->detectsum, FALSE, DEFAULT_DETECTSUM, NULL, NULL) );
2595 	
2596 	   /* "extended" formulations of a concave expressions can give worse estimators */
2597 	   nlhdlrdata->extendedform = FALSE;
2598 	
2599 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/cvxquadratic",
2600 	      "whether to use convexity check on quadratics",
2601 	      &nlhdlrdata->cvxquadratic, TRUE, DEFAULT_CVXQUADRATIC_CONCAVE, NULL, NULL) );
2602 	
2603 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/cvxsignomial",
2604 	      "whether to use convexity check on signomials",
2605 	      &nlhdlrdata->cvxsignomial, TRUE, DEFAULT_CVXSIGNOMIAL, NULL, NULL) );
2606 	
2607 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/cvxprodcomp",
2608 	      "whether to use convexity check on product composition f(h)*h",
2609 	      &nlhdlrdata->cvxprodcomp, TRUE, DEFAULT_CVXPRODCOMP, NULL, NULL) );
2610 	
2611 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/handletrivial",
2612 	      "whether to also handle trivial convex expressions",
2613 	      &nlhdlrdata->handletrivial, TRUE, DEFAULT_HANDLETRIVIAL, NULL, NULL) );
2614 	
2615 	   SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrfreeHdlrDataConvexConcave);
2616 	   SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrConcave);
2617 	   SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrfreeExprDataConvexConcave);
2618 	   SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaConcave, NULL, nlhdlrEstimateConcave, NULL);
2619 	   SCIPnlhdlrSetInitExit(nlhdlr, NULL, nlhdlrExitConcave);
2620 	
2621 	   return SCIP_OKAY;
2622 	}
2623 	
2624 	/** checks whether a given expression is convex or concave w.r.t. the original variables
2625 	 *
2626 	 * This function uses the methods that are used in the detection algorithm of the convex nonlinear handler.
2627 	 */
2628 	SCIP_RETCODE SCIPhasExprCurvature(
2629 	   SCIP*                 scip,               /**< SCIP data structure */
2630 	   SCIP_EXPR*            expr,               /**< expression */
2631 	   SCIP_EXPRCURV         curv,               /**< curvature to check for */
2632 	   SCIP_Bool*            success,            /**< buffer to store whether expression has curvature curv (w.r.t. original variables) */
2633 	   SCIP_HASHMAP*         assumevarfixed      /**< hashmap containing variables that should be assumed to be fixed, or NULL */
2634 	   )
2635 	{
2636 	   SCIP_NLHDLRDATA nlhdlrdata;
2637 	   SCIP_EXPR* rootnlexpr;
2638 	   SCIP_HASHMAP* nlexpr2origexpr;
2639 	   int nleafs;
2640 	
2641 	   assert(expr != NULL);
2642 	   assert(curv != SCIP_EXPRCURV_UNKNOWN);
2643 	   assert(success != NULL);
2644 	
2645 	   /* create temporary hashmap */
2646 	   SCIP_CALL( SCIPhashmapCreate(&nlexpr2origexpr, SCIPblkmem(scip), 20) );
2647 	
2648 	   /* prepare nonlinear handler data */
2649 	   nlhdlrdata.isnlhdlrconvex = TRUE;
2650 	   nlhdlrdata.evalsol = NULL;
2651 	   nlhdlrdata.detectsum = TRUE;
2652 	   nlhdlrdata.extendedform = FALSE;
2653 	   nlhdlrdata.cvxquadratic = TRUE;
2654 	   nlhdlrdata.cvxsignomial = TRUE;
2655 	   nlhdlrdata.cvxprodcomp = TRUE;
2656 	   nlhdlrdata.handletrivial = TRUE;
2657 	
2658 	   SCIP_CALL( constructExpr(scip, &nlhdlrdata, &rootnlexpr, nlexpr2origexpr, &nleafs, expr, curv, assumevarfixed, FALSE, success) );
2659 	
2660 	   /* free created expression */
2661 	   if( rootnlexpr != NULL )
2662 	   {
2663 	      SCIP_CALL( SCIPreleaseExpr(scip, &rootnlexpr) );
2664 	   }
2665 	
2666 	   /* free hashmap */
2667 	   SCIPhashmapFree(&nlexpr2origexpr);
2668 	
2669 	   return SCIP_OKAY;
2670 	}
2671