1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright 2002-2022 Zuse Institute Berlin                                */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   nlhdlr_signomial.c
26   	 * @ingroup DEFPLUGINS_NLHDLR
27   	 * @brief  signomial nonlinear handler
28   	 * @author Liding Xu
29   	 */
30   	
31   	#ifdef SCIP_SIGCUT_DEBUG
32   	#ifndef SCIP_SIG_DEBUG
33   	#define SCIP_SIG_DEBUG
34   	#endif
35   	#endif
36   	
37   	#ifdef SCIP_SIG_DEBUG
38   	#define SCIP_DEBUG
39   	#endif
40   	
41   	#include "scip/nlhdlr_signomial.h"
42   	#include "scip/cons_nonlinear.h"
43   	#include "scip/pub_misc_rowprep.h"
44   	#include "scip/pub_nlhdlr.h"
45   	#include "scip/scip_expr.h"
46   	#include "scip/expr_var.h"
47   	#include "scip/expr_pow.h"
48   	
49   	/* fundamental nonlinear handler properties */
50   	#define NLHDLR_NAME                    "signomial"
51   	#define NLHDLR_DESC                    "handler for signomial expressions"
52   	#define NLHDLR_DETECTPRIORITY          30
53   	#define NLHDLR_ENFOPRIORITY            30
54   	
55   	/* handler specific parameters */
56   	#define NLHDLR_MAXNUNDERVARS           14    /**< maximum number of variables when underestimating a concave power function (maximum: 14) */
57   	#define NLHDLR_MINCUTSCALE             1e-5  /**< minimum scale factor when scaling a cut (minimum: 1e-6) */
58   	
59   	
60   	/** nonlinear handler expression data
61   	 *
62   	 * A signomial expression admits the form \f$ cx^a = y \f$, where \f$ y \f$ is an auxiliary variable representing this
63   	 * expression and all \f$ x \f$ are positive.
64   	 *
65   	 * The natural formulation of the expression is defined as \f$ x^a = t = y/c \f$, where \f$ t \f$ is a
66   	 * non-existant slack variable denoting \f$ y/c \f$. The variables in \f$x\f$ with positive exponents form positive
67   	 * variables \f$ u \f$, and the associated exponents form positive exponents \f$ f \f$. The variables in \f$ x \f$ with
68   	 * negative exponents and \f$ t \f$  form negative variables \f$ v \f$, and the associated exponents form negative
69   	 * exponents \f$ g \f$. Let \f$ s =  \max(|f|,|g|) \f$ be a normalization constant, where \f$ |\cdot| \f$ denotes the L1 norm. Apply a scaling step:
70   	 * Dividing the entries of \f$ f \f$  by \f$ s \f$, and dividing the entries of \f$ g \f$  by \f$ s \f$ as well. Then \f$ x^a = t \f$ has
71   	 * a reformulation \f$ u^f = v^g \f$, where \f$ u^f, v^g \$ are two concave power functions.
72   	 */
73   	struct SCIP_NlhdlrExprData
74   	{
75   	   SCIP_Real             coef;               /**< coefficient \f$c\f$ */
76   	   SCIP_EXPR**           factors;            /**< expression factors representing \f$x\f$ */
77   	   int                   nfactors;           /**< number of factors */
78   	   int                   nvars;              /**< number of variables \f$(x,y)\f$ */
79   	   SCIP_Real*            exponents;          /**< exponents \f$e\f$ */
80   	   int                   nposvars;           /**< number of positive variables \f$u\f$ */
81   	   int                   nnegvars;           /**< number of negative variables  \f$v\f$ */
82   	   SCIP_Bool*            signs;              /**< indicators for sign of variables after reformulation, TRUE for positive, FALSE for negative */
83   	   SCIP_Real*            refexponents;       /**< exponents of \f$(x,t)\f$ after reformulation (always positive) */
84   	   SCIP_Bool             isstorecapture;     /**< are all variables already got? */
85   	
86   	   /* working parameters will be modified after getting all variables */
87   	   SCIP_VAR**            vars;               /**< variables \f$(x,y)\f$ */
88   	   SCIP_INTERVAL*        intervals;          /**< intervals storing lower and upper bounds of variables \f$(x,y)\f$ */
89   	   SCIP_Real*            box;                /**< the upper/lower bounds of variables, used in SCIPcomputeFacetVertexPolyhedralNonlinear() */
90   	   SCIP_Real*            xstar;              /**< the values of variables, used in SCIPcomputeFacetVertexPolyhedralNonlinear() */
91   	   SCIP_Real*            facetcoefs;         /**< the coefficients of variables, returned by SCIPcomputeFacetVertexPolyhedralNonlinear() */
92   	};
93   	
94   	/** nonlinear handler data */
95   	struct SCIP_NlhdlrData
96   	{
97   	   /* parameters */
98   	   int                   maxnundervars;      /**< maximum number of variables in underestimating a concave power function */
99   	   SCIP_Real             mincutscale;        /**< minimum scale factor when scaling a cut */
100  	};
101  	
102  	/** data struct to be passed on to vertexpoly-evalfunction (see SCIPcomputeFacetVertexPolyhedralNonlinear) */
103  	typedef struct
104  	{
105  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata;     /**< expression data */
106  	   SCIP_Bool             sign;               /**< the sign of variables in the reformulated constraint for vertexpoly-evalfunction */
107  	   int                   nsignvars;          /**< the number of variables in the reformulated constraint for vertexpoly-evalfunction */
108  	   SCIP*                 scip;               /**< SCIP data structure */
109  	} VERTEXPOLYFUN_EVALDATA;
110  	
111  	
112  	/*
113  	 * Local methods
114  	 */
115  	
116  	#ifdef SCIP_SIGCUT_DEBUG
117  	
118  	/** print the information on a signomial term */
119  	static
120  	SCIP_RETCODE printSignomial(
121  	   SCIP*                 scip,               /**< SCIP data structure */
122  	   SCIP_EXPR*            expr,               /**< expression */
123  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata      /**< expression data */
124  	   )
125  	{
126  	   assert(expr != NULL);
127  	
128  	   /* print overall statistics and the expression */
129  	   SCIPdebugMsg(scip, " #all variables: %d, #positive exponent variables: %d, #negative exponent variables: %d, auxvar: %s \n expr: ",
130  	      nlhdlrexprdata->nvars, nlhdlrexprdata->nposvars, nlhdlrexprdata->nnegvars, SCIPvarGetName(SCIPgetExprAuxVarNonlinear(expr)) );
131  	   SCIPprintExpr(scip, expr, NULL);
132  	
133  	   /* if variables are detected, we can print more information */
134  	   if( !nlhdlrexprdata->isstorecapture )
135  	   {
136  	      SCIPinfoMessage(scip, NULL, "\n");
137  	      return SCIP_OKAY;
138  	   }
139  	
140  	   /* print the natural formulation of the expression */
141  	   SCIPinfoMessage(scip, NULL, ". natural formulation c x^a = y: %.2f", nlhdlrexprdata->coef);
142  	   for( int i = 0; i < nlhdlrexprdata->nvars - 1; i++ )
143  	   {
144  	      SCIPinfoMessage(scip, NULL, "%s^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->exponents[i]);
145  	   }
146  	   SCIPinfoMessage(scip, NULL, " = %s", SCIPvarGetName(nlhdlrexprdata->vars[nlhdlrexprdata->nvars - 1]));
147  	
148  	   /* print the reformulation of the expression */
149  	   SCIPinfoMessage(scip, NULL, ". Reformulation u^f = v^g: ");
150  	   if( nlhdlrexprdata->nposvars == 0 )
151  	   {
152  	      SCIPinfoMessage(scip, NULL, "1");
153  	   }
154  	   else
155  	   {
156  	      for( int i = 0; i < nlhdlrexprdata->nvars; i++ )
157  	         if( nlhdlrexprdata->signs[i] )
158  	         {
159  	            SCIPinfoMessage(scip, NULL, "%s^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->refexponents[i]);
160  	         }
161  	   }
162  	   SCIPinfoMessage(scip, NULL, " = ");
163  	   if( nlhdlrexprdata->nnegvars == 0 )
164  	   {
165  	      SCIPinfoMessage(scip, NULL, "1");
166  	   }
167  	   else
168  	   {
169  	      for( int i = 0; i < nlhdlrexprdata->nvars; i++ )
170  	      {
171  	         if( !nlhdlrexprdata->signs[i] )
172  	         {
173  	            if( i == nlhdlrexprdata->nvars - 1 )
174  	            {
175  	               SCIPinfoMessage(scip, NULL, "(%s/%.2f)^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->coef, nlhdlrexprdata->refexponents[i]);
176  	            }
177  	            else
178  	            {
179  	               SCIPinfoMessage(scip, NULL, "%s^%.2f", SCIPvarGetName(nlhdlrexprdata->vars[i]), nlhdlrexprdata->refexponents[i]);
180  	            }
181  	         }
182  	      }
183  	   }
184  	   SCIPinfoMessage(scip, NULL, "\n");
185  	
186  	   return SCIP_OKAY;
187  	}
188  	
189  	#endif
190  	
191  	/** free the memory of expression data */
192  	static
193  	void freeExprDataMem(
194  	   SCIP*                 scip,               /**< SCIP data structure */
195  	   SCIP_NLHDLREXPRDATA** nlhdlrexprdata,     /**< expression data */
196  	   SCIP_Bool             ispartial           /**< free the partially allocated memory or the fully allocated memory? */
197  	   )
198  	{
199  	
200  	   SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrexprdata)->factors, (*nlhdlrexprdata)->nfactors);
201  	   SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->exponents, (*nlhdlrexprdata)->nfactors);
202  	   if( !ispartial )
203  	   {
204  	      int nvars = (*nlhdlrexprdata)->nvars;
205  	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->signs, nvars);
206  	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->refexponents, nvars);
207  	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, nvars);
208  	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->intervals, nvars);
209  	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->box, 2*nvars);  /*lint !e647*/
210  	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->xstar, nvars);
211  	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->facetcoefs, nvars);
212  	   }
213  	   SCIPfreeBlockMemory(scip, nlhdlrexprdata);
214  	   *nlhdlrexprdata = NULL;
215  	}
216  	
217  	/** reforms a rowprep to a standard form for nonlinear handlers
218  	 *
219  	 * The input rowprep is of the form  \f$ a_u u + a_v v + b \f$.
220  	 * If in the overestimate mode, then we relax \f$ t \le x^a \f$, i.e., \f$ 0 \le u^f - v^g \f$. This implies that \f$ t \f$ is in \f$ v = (v',t) \f$.
221  	 * Therefore, the valid inequality is \f$ 0 \le a_u u + a_v v + b \f$. As overestimate mode requires that \f$ t \f$ is in the left hand side,
222  	 * the coefficients of \f$ t \f$ must be made negative while keeping the sign of the inequality, we can show that \f$ a_t \le 0 \f$, so it suffices
223  	 * to divide the both sides by \f$ -a_t \ge 0\f$, which yields  \f$ t \le (a_u u + a_v' v' + b) / -a_t \f$.
224  	 * A rowprep in standard form only contains an estimator of the expression and no auxvar.
225  	 * If in the underestimate mode, then we relax \f$ x^a \le t \f$, i.e., \f$ u^f - v^g \le 0 \f$. This implies that \f$ t \f$ is in \f$ v = (v',t) \f$.
226  	 * Therefore, the valid inequality is \f$ a_u u + a_v v + b \le 0 \f$. As overestimate mode requires that \f$ t \f$ is in the left hand side,
227  	 * the coefficients of \f$ t \f$ must be made negative while keeping the sign of the inequality, we can show that \f$ a_t \le 0 \f$, so it suffices
228  	 * to divide the both sides by \f$ -a_t \ge 0 \f$, which yields \f$ (a_u u + a_v' v' + b) / -a_t \le t \f$.
229  	 * A rowprep in standard form only contains an estimator of the expression and no auxvar.
230  	 */
231  	static
232  	void reformRowprep(
233  	   SCIP*                 scip,               /**< SCIP data structure */
234  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< expression data */
235  	   SCIP_ROWPREP*         rowprep,            /**< cut to be reformulated */
236  	   SCIP_Real             mincutscale,        /**< min scaling factor for the cut in rowprep */
237  	   SCIP_Bool*            success             /**< pointer to store whether the reformulating was successful */
238  	)
239  	{
240  	   int i;
241  	   int nvars;
242  	   SCIP_Real coefauxvar;
243  	   SCIP_Real scale;
244  	   SCIP_Real* coefs;
245  	   SCIP_VAR* auxvar;
246  	   SCIP_VAR** vars;
247  	
248  	   assert(rowprep != NULL);
249  	   assert(nlhdlrexprdata != NULL);
250  	
251  	   coefs = SCIProwprepGetCoefs(rowprep);
252  	   vars = SCIProwprepGetVars(rowprep);
253  	   nvars = SCIProwprepGetNVars(rowprep);
254  	
255  	   /* find auxvar's cut coefficient and set it to zero */
256  	   auxvar = nlhdlrexprdata->vars[nlhdlrexprdata->nfactors];
257  	   coefauxvar = 1.0;
258  	   for( i = 0; i < nvars; i++ )
259  	   {
260  	      if( vars[i] == auxvar )
261  	      {
262  	         coefauxvar = coefs[i];
263  	         coefs[i] = 0.0;
264  	         break;
265  	      }
266  	   }
267  	
268  	   if( SCIPisZero(scip, coefauxvar) || coefauxvar >= 0.0 )
269  	   {
270  	      *success = FALSE;
271  	   }
272  	   else
273  	   {
274  	      /* the reformation scales the cut so that coefficients and constant are divided by the absolute value of coefauxvar */
275  	      assert(coefauxvar < 0.0);
276  	      scale = -1.0 / coefauxvar;
277  	      for( i = 0; i < nvars; i++ )
278  	      {
279  	         if( vars[i] == auxvar )
280  	            continue;
281  	         coefs[i] *= scale;
282  	      }
283  	      /* set the side to scale*side by adding -side and adding scale*side */
284  	      SCIProwprepAddSide(rowprep, SCIProwprepGetSide(rowprep) * (-1.0 + scale));
285  	
286  	      *success = SCIPisGT(scip, scale, mincutscale);
287  	   }
288  	}
289  	
290  	/** store and capture variables associated with the expression and its subexpressions */
291  	static
292  	SCIP_RETCODE storeCaptureVars(
293  	   SCIP*                 scip,               /**< SCIP data structure */
294  	   SCIP_EXPR*            expr,               /**< expression */
295  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata      /**< expression data */
296  	   )
297  	{
298  	   int c;
299  	
300  	   assert(!nlhdlrexprdata->isstorecapture);
301  	
302  	   /* get and capture variables \f$x\f$ */
303  	   for( c = 0; c < nlhdlrexprdata->nfactors; ++c )
304  	   {
305  	      nlhdlrexprdata->vars[c] = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->factors[c]);
306  	      assert(nlhdlrexprdata->vars[c] != NULL);
307  	
308  	      SCIP_CALL( SCIPcaptureVar(scip, nlhdlrexprdata->vars[c]) );
309  	   }
310  	
311  	   /* get and capture variable \f$y\f$ */
312  	   nlhdlrexprdata->vars[c] = SCIPgetExprAuxVarNonlinear(expr);
313  	   assert(nlhdlrexprdata->vars[c] != NULL);
314  	   SCIP_CALL( SCIPcaptureVar(scip, nlhdlrexprdata->vars[c]) );
315  	
316  	   nlhdlrexprdata->isstorecapture = TRUE;
317  	
318  	   return SCIP_OKAY;
319  	}
320  	
321  	/** get bounds of variables x,t and check whether they are box constrained signomial variables */
322  	static
323  	void checkSignomialBounds(
324  	   SCIP*                 scip,               /**< SCIP data structure */
325  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< expression data */
326  	   SCIP_Bool*            isboxsignomial      /**< buffer to store whether variables are box constrained signomial variables */
327  	   )
328  	{
329  	   int c;
330  	   SCIP_Real powinf;
331  	   SCIP_Real powsup;
332  	   SCIP_Real productinf = 1;
333  	   SCIP_Real productsup = 1;
334  	
335  	   assert(nlhdlrexprdata->isstorecapture);
336  	
337  	   *isboxsignomial = FALSE;
338  	
339  	   /* get bounds of x */
340  	   for( c = 0; c < nlhdlrexprdata->nfactors; c++ )
341  	   {
342  	      SCIP_Real inf = SCIPvarGetLbLocal(nlhdlrexprdata->vars[c]);
343  	      SCIP_Real sup = SCIPvarGetUbLocal(nlhdlrexprdata->vars[c]);
344  	
345  	      /* if the bounds of the variable are not positive and finite, or (bounds equal) then the expression is not a signomial */
346  	      if( !SCIPisPositive(scip, inf) || !SCIPisPositive(scip, sup) || SCIPisInfinity(scip, sup) || SCIPisEQ(scip, inf, sup) )
347  	         return;
348  	
349  	      nlhdlrexprdata->intervals[c].inf = inf;
350  	      nlhdlrexprdata->intervals[c].sup = MAX(sup, inf + 0.1);
351  	      powinf = pow(inf, nlhdlrexprdata->exponents[c]);
352  	      powsup = pow(sup, nlhdlrexprdata->exponents[c]);
353  	      productinf *= MIN(powinf, powsup);
354  	      productsup *= MAX(powinf, powsup);
355  	   }
356  	
357  	   /* compute bounds of t by bounds of x */
358  	   nlhdlrexprdata->intervals[c].inf = productinf;
359  	   nlhdlrexprdata->intervals[c].sup = productsup;
360  	
361  	   *isboxsignomial = TRUE;
362  	}
363  	
364  	/** evaluate expression at solution w.r.t. auxiliary variables */
365  	static
366  	SCIP_DECL_VERTEXPOLYFUN(nlhdlrExprEvalPower)
367  	{
368  	   int i;
369  	   int j;
370  	   SCIP_Real val;
371  	   VERTEXPOLYFUN_EVALDATA* evaldata;
372  	   SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
373  	
374  	   assert(args != NULL);
375  	
376  	   evaldata = (VERTEXPOLYFUN_EVALDATA *)funcdata;
377  	
378  	   assert(evaldata != NULL);
379  	   assert(nargs == evaldata->nsignvars);
380  	
381  	   nlhdlrexprdata = evaldata->nlhdlrexprdata;
382  	
383  	#ifdef SCIP_MORE_DEBUG
384  	   SCIPdebugMsg(evaldata->scip, "eval vertexpolyfun\n");
385  	#endif
386  	   val = 1.0;
387  	   for( i = 0, j = 0; i < nlhdlrexprdata->nvars; ++i )
388  	   {
389  	      if( nlhdlrexprdata->signs[i] != evaldata->sign )
390  	         continue;
391  	      /* the reformulated exponent of args[j] is found */
392  	      val *= pow(args[j], nlhdlrexprdata->refexponents[i]);
393  	      j++;
394  	   }
395  	
396  	   return val;
397  	}
398  	
399  	/** determine whether a power function \f$ w^h \f$ is special and add an overunderestimator or underestimator to a given rowprep
400  	 *
401  	 * \f$ w^h \f$ is special, if all variables are fixed, or it is a constant to estimate, a univariate power to estimate,
402  	 * or a bivariate power to underestimate. The estimator is multiplied by the multiplier and stored in the rowprep.
403  	 */
404  	static
405  	SCIP_RETCODE estimateSpecialPower(
406  	   SCIP*                 scip,               /**< SCIP data structure */
407  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
408  	   SCIP_Bool             sign,               /**< the sign of variables of the power function */
409  	   SCIP_Real             multiplier,         /**< the multiplier of the estimator */
410  	   SCIP_Bool             overestimate,       /**< whether overestimate or underestimator the power function */
411  	   SCIP_SOL*             sol,                /**< solution \f$ w' \f$  to use in estimation */
412  	   SCIP_ROWPREP*         rowprep,            /**< rowprep where to store estimator */
413  	   SCIP_Bool*            isspecial,          /**< buffer to store whether this function is a special case */
414  	   SCIP_Bool*            success             /**< buffer to store whether successful */
415  	   )
416  	{
417  	   int i;
418  	   SCIP_Bool allfixed = TRUE;
419  	   int nsignvars;
420  	
421  	   assert(scip != NULL);
422  	   assert(nlhdlrexprdata != NULL);
423  	   assert(rowprep != NULL);
424  	   assert(isspecial != NULL);
425  	   assert(success != NULL);
426  	
427  	   *success = FALSE;
428  	   /* check whether all variables are fixed */
429  	   for( i = 0; i < nlhdlrexprdata->nvars; ++i )
430  	   {
431  	      if( nlhdlrexprdata->signs[i] != sign )
432  	         continue;
433  	
434  	      if( !SCIPisRelEQ(scip, nlhdlrexprdata->intervals[i].inf, nlhdlrexprdata->intervals[i].sup) )
435  	      {
436  	         allfixed = FALSE;
437  	         break;
438  	      }
439  	   }
440  	
441  	   /* allfixed is a special case */
442  	   if( allfixed )
443  	   {
444  	      /* return a constant */
445  	      SCIP_Real funcval = 1.0;
446  	      SCIP_Real scale;
447  	      SCIP_Real val;
448  	
449  	      for( i = 0; i < nlhdlrexprdata->nvars; ++i )
450  	      {
451  	         if( nlhdlrexprdata->signs[i] != sign )
452  	            continue;
453  	
454  	         scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0;
455  	         val = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale;
456  	         funcval *= pow(val, nlhdlrexprdata->refexponents[i]);
457  	      }
458  	      SCIProwprepAddConstant(rowprep, multiplier * funcval);
459  	      *isspecial = TRUE;
460  	
461  	      return SCIP_OKAY;
462  	   }
463  	
464  	   /* number of variables in the power function */
465  	   nsignvars = sign ? nlhdlrexprdata->nposvars : nlhdlrexprdata->nnegvars;
466  	
467  	   /* if the power function has no more than 2 variables, this a special case */
468  	   *isspecial = ( nsignvars <= 1 ) || ( nsignvars == 2  && !overestimate );
469  	   if( !*isspecial )
470  	      return SCIP_OKAY;
471  	
472  	   if( nsignvars == 0 )
473  	   {
474  	      /* constant case */
475  	      SCIProwprepAddConstant(rowprep, multiplier);
476  	      *success = TRUE;
477  	   }
478  	   else if( nsignvars == 1 )
479  	   {
480  	      /* univariate case, w^h */
481  	      for( i = 0; i < nlhdlrexprdata->nvars; ++i )
482  	      {
483  	         if( nlhdlrexprdata->signs[i] == sign )
484  	         {
485  	            /* the variable w is found */
486  	            SCIP_Real scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1;
487  	            SCIP_VAR* var = nlhdlrexprdata->vars[i];
488  	            SCIP_Real refexponent = nlhdlrexprdata->refexponents[i];
489  	            if( refexponent == 1.0 )
490  	            {
491  	               /* h = 1, a univariate linear function. Only rescale, no need for linearization */
492  	               SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, multiplier / scale) );
493  	            }
494  	            else
495  	            {
496  	               /* a univariate power function */
497  	               SCIP_Real facetconstant;
498  	               SCIP_Real facetcoef;
499  	               SCIP_Real val = SCIPgetSolVal(scip, sol, var) / scale;
500  	               /* local (using bounds) depends on whether we under- or overestimate */
501  	               SCIP_Bool islocal = !overestimate;
502  	               SCIPestimateRoot(scip, refexponent, overestimate, nlhdlrexprdata->intervals[i].inf, nlhdlrexprdata->intervals[i].sup,
503  	                  val, &facetconstant, &facetcoef, &islocal, success);
504  	               SCIProwprepAddConstant(rowprep,  multiplier * facetconstant);
505  	               SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, multiplier * facetcoef / scale) );
506  	            }
507  	         }
508  	      }
509  	      *success = TRUE;
510  	   }
511  	   else if( nsignvars == 2 && !overestimate )
512  	   {
513  	      /* bivariate case, f(w) = w^h = f_0(w_0) f_1(w_1). The convex under envelope is the maxmimum of
514  	       * two affine functions, each of which is determined by three extreme points the box bound.
515  	       * One affine function is supported by three lower left points; the other affine function is
516  	       * supported by three upper right points. For a point xstar in the box, its corresponding affine function can be determined by
517  	       * which region (upper right or lower left half space) the point is in. Thus, we can determine the region, and use the
518  	       * associated three extreme point to interpolate the affine function.
519  	       */
520  	      SCIP_Bool isupperright;
521  	      SCIP_Real dw0, dw1;
522  	      SCIP_Real f0w0l, f0w0u, f1w1l, f1w1u;
523  	      SCIP_Real fw0lw1u, fw0uw1l;
524  	      SCIP_Real facetconstant;
525  	      SCIP_Real facetcoefs[2] = {0.0, 0.0};
526  	      SCIP_VAR* vars[2] = {NULL, NULL};
527  	      SCIP_Real refexponents[2] = {0.0, 0.0};
528  	      SCIP_Real xstar[2] = {0.0, 0.0};;
529  	      SCIP_Real scale[2] = {0.0, 0.0};;
530  	      SCIP_Real box[4] = {0.0, 0.0, 0.0, 0.0};
531  	      int j = 0;
532  	
533  	      /* get refexponents, xstar, scale, and box */
534  	      for( i = 0; i < nlhdlrexprdata->nvars; ++i )
535  	      {
536  	         if( nlhdlrexprdata->signs[i] != sign )
537  	            continue;
538  	         box[2 * j] = nlhdlrexprdata->intervals[i].inf;
539  	         box[2 * j + 1] = nlhdlrexprdata->intervals[i].sup;
540  	         refexponents[j] = nlhdlrexprdata->refexponents[i];
541  	         scale[j] = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1;
542  	         vars[j] = nlhdlrexprdata->vars[i];
543  	         xstar[j] = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[j]) / scale[j];
544  	         j++;
545  	      }
546  	
547  	      /* compute the box length*/
548  	      dw0 = box[1] - box[0];
549  	      dw1 = box[3] - box[2];
550  	
551  	      /* determine the location (upper right or lower left half sapce) of the xstar.
552  	       * (dw1, dw0) is the direction vector to the upper right half space.
553  	       */
554  	      isupperright = ( (xstar[0] - box[0]) * dw1 + (xstar[1] - box[3]) * dw0 ) > 0.0;
555  	
556  	      /* compute function values of f_0, f_1 at vertices */
557  	      f0w0l = pow(box[0], refexponents[0]);
558  	      f0w0u = pow(box[1], refexponents[0]);
559  	      f1w1l = pow(box[2], refexponents[1]);
560  	      f1w1u = pow(box[3], refexponents[1]);
561  	      fw0lw1u = f0w0l * f1w1u;
562  	      fw0uw1l = f0w0u * f1w1l;
563  	      if( isupperright )
564  	      {
565  	         /* underestimator: fw0uw1u + (fw0uw1u - fw0lw1u) / (dw0) * (x0 - w0u) + (fw0uw1u - fw0uw1l) / (dw1) * (x1 - w1u) */
566  	         SCIP_Real fw0uw1u = f0w0u * f1w1u;
567  	         facetcoefs[0] = (fw0uw1u - fw0lw1u) / dw0;
568  	         facetcoefs[1] = (fw0uw1u - fw0uw1l) / dw1;
569  	         facetconstant = fw0uw1u  - facetcoefs[0] * box[1] - facetcoefs[1] * box[3];
570  	      }
571  	      else
572  	      {
573  	         /* underestimator: fw0lw1l + (fw0uw1l - fw0lw1l) / (dw0) * (x0 - w0l) + (fw0lw1u- fw0lw1l) / (dw1) * (x1 - w1l) */
574  	         SCIP_Real fw0lw1l = f0w0l * f1w1l;
575  	         facetcoefs[0] = (fw0uw1l - fw0lw1l) / dw0;
576  	         facetcoefs[1] = (fw0lw1u - fw0lw1l) / dw1;
577  	         facetconstant = fw0lw1l  - facetcoefs[0] * box[0] - facetcoefs[1] * box[2];
578  	      }
579  	      SCIProwprepAddConstant(rowprep,  multiplier * facetconstant);
580  	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, vars[0], multiplier * facetcoefs[0] / scale[0]) );
581  	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, vars[1], multiplier * facetcoefs[1] / scale[1]) );
582  	      *success = TRUE;
583  	   }
584  	
585  	   return SCIP_OKAY;
586  	}
587  	
588  	/** adds an underestimator for a multivariate concave power function \f$ w^h \f$ to a given rowprep
589  	 *
590  	 * Calls \ref SCIPcomputeFacetVertexPolyhedralNonlinear() for \f$ w^h \f$  and
591  	 * box set to local bounds of auxiliary variables. The estimator is multiplied
592  	 * by the multiplier and stored in the rowprep.
593  	 */
594  	static
595  	SCIP_RETCODE underEstimatePower(
596  	   SCIP*                 scip,               /**< SCIP data structure */
597  	   SCIP_CONSHDLR*        conshdlr,           /**< nonlinear constraint handler */
598  	   SCIP_NLHDLR*          nlhdlr,             /**< nonlinear handler */
599  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
600  	   SCIP_Bool             sign,               /**< the sign of variables of the power function */
601  	   SCIP_Real             multiplier,         /**< the multiplier of the estimator */
602  	   SCIP_SOL*             sol,                /**< solution \f$ w' \f$ to use*/
603  	   SCIP_Real             targetvalue,        /**< a target value to achieve; if not reachable, then can give up early */
604  	   SCIP_ROWPREP*         rowprep,            /**< rowprep where to store estimator */
605  	   SCIP_Bool*            success             /**< buffer to store whether successful */
606  	   )
607  	{
608  	   int i;
609  	   int j;
610  	   int nsignvars;
611  	   SCIP_Real facetconstant;
612  	   SCIP_Real scale;
613  	   SCIP_Real* box;
614  	   SCIP_Real* facetcoefs;
615  	   SCIP_Real* xstar;
616  	   VERTEXPOLYFUN_EVALDATA evaldata;
617  	
618  	   assert(scip != NULL);
619  	   assert(nlhdlr != NULL);
620  	   assert(nlhdlrexprdata != NULL);
621  	   assert(rowprep != NULL);
622  	   assert(success != NULL);
623  	
624  	   *success = FALSE;
625  	
626  	   /* number of variables of sign */
627  	   nsignvars = sign ? nlhdlrexprdata->nposvars : nlhdlrexprdata->nnegvars;
628  	
629  	   /* data structure to evaluate the power function */
630  	   evaldata.nlhdlrexprdata = nlhdlrexprdata;
631  	   evaldata.sign = sign;
632  	   evaldata.nsignvars = nsignvars;
633  	   evaldata.scip = scip;
634  	
635  	   /* compute box constraints and reference value of variables*/
636  	   xstar = nlhdlrexprdata->xstar;
637  	   box = nlhdlrexprdata->box;
638  	   j = 0;
639  	   for( i = 0; i < nlhdlrexprdata->nvars; ++i )
640  	   {
641  	      if( nlhdlrexprdata->signs[i] != sign )
642  	         continue;
643  	
644  	      box[2 * j] = nlhdlrexprdata->intervals[i].inf;
645  	      box[2 * j + 1] = nlhdlrexprdata->intervals[i].sup;
646  	      scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0;
647  	      xstar[j] = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale;
648  	      j++;
649  	   }
650  	
651  	   /* find a facet of the underestimator */
652  	   facetcoefs = nlhdlrexprdata->facetcoefs;
653  	   SCIP_CALL( SCIPcomputeFacetVertexPolyhedralNonlinear(scip, conshdlr, FALSE, nlhdlrExprEvalPower, (void*)&evaldata, xstar, box,
654  	      nsignvars, targetvalue, success, facetcoefs, &facetconstant) );
655  	
656  	   if( !*success )
657  	   {
658  	      SCIPdebugMsg(scip, "failed to compute facet of convex hull\n");
659  	      return SCIP_OKAY;
660  	   }
661  	
662  	   /* set coefficients in the rowprep */
663  	   SCIProwprepAddConstant(rowprep, multiplier * facetconstant);
664  	   j = 0;
665  	   for( i = 0; i < nlhdlrexprdata->nvars; ++i )
666  	   {
667  	      if( nlhdlrexprdata->signs[i] != sign )
668  	         continue;
669  	      scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0;
670  	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, nlhdlrexprdata->vars[i], multiplier * facetcoefs[j] / scale) );
671  	      j++;
672  	   }
673  	
674  	   return SCIP_OKAY;
675  	}
676  	
677  	/** adds an overestimator for a concave power function \f$ w^h \f$ to a given rowprep
678  	 *
679  	 * The estimator is multiplied by the multiplier and stored in the rowprep.
680  	 */
681  	static
682  	SCIP_RETCODE overEstimatePower(
683  	   SCIP*                 scip,               /**< SCIP data structure */
684  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nonlinear handler expression data */
685  	   SCIP_Bool             sign,               /**< the sign of variables of the power function */
686  	   SCIP_Real             multiplier,         /**< the multiplier of the estimator */
687  	   SCIP_SOL*             sol,                /**< solution to use */
688  	   SCIP_ROWPREP*         rowprep,            /**< rowprep where to store estimator */
689  	   SCIP_Bool*            success             /**< buffer to store whether successful */
690  	   )
691  	{
692  	   int i;
693  	   SCIP_Real facetcoef;
694  	   SCIP_Real facetconstant;
695  	   SCIP_Real funcval;
696  	   SCIP_Real refexponent;
697  	   SCIP_Real sumrefexponents;
698  	   SCIP_Real scale;
699  	   SCIP_Real val;
700  	   SCIP_VAR* var;
701  	   assert(scip != NULL);
702  	   assert(nlhdlrexprdata != NULL);
703  	   assert(rowprep != NULL);
704  	   assert(success != NULL);
705  	
706  	   *success = FALSE;
707  	
708  	   /* compute the value and the sum of reformulated exponents of the power function */
709  	   sumrefexponents = 0;
710  	   funcval = 1;
711  	   for( i = 0; i < nlhdlrexprdata->nvars; ++i )
712  	   {
713  	      if( nlhdlrexprdata->signs[i] != sign )
714  	         continue;
715  	      scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0;
716  	      val = SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale;
717  	      val = MAX(val, 0.1);
718  	      refexponent = nlhdlrexprdata->refexponents[i];
719  	      sumrefexponents += refexponent;
720  	      funcval *= pow(val, refexponent);
721  	   }
722  	
723  	   /* overestimate by gradient cut: w'^h + h w'^(h - 1)(w - w') */
724  	   facetconstant = (1 - sumrefexponents) * funcval;
725  	   SCIProwprepAddConstant(rowprep,  multiplier * facetconstant);
726  	   for( i = 0; i < nlhdlrexprdata->nvars; ++i )
727  	   {
728  	      if( nlhdlrexprdata->signs[i] != sign )
729  	         continue;
730  	      scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0;
731  	      var = nlhdlrexprdata->vars[i];
732  	      val = SCIPgetSolVal(scip, sol, var) / scale;
733  	      val = MAX(val, 0.1);
734  	      facetcoef = nlhdlrexprdata->refexponents[i] * funcval / val;
735  	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, multiplier * facetcoef / scale) );
736  	   }
737  	
738  	   *success = TRUE;
739  	
740  	   return SCIP_OKAY;
741  	}
742  	
743  	/*
744  	 * Callback methods of nonlinear handler
745  	 */
746  	
747  	/** nonlinear handler estimation callback */
748  	static
749  	SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateSignomial)
750  	{ /*lint --e{715}*/
751  	   SCIP_Bool isboxsignomial;
752  	   SCIP_Bool isspecial;
753  	   SCIP_Bool undersign;
754  	   SCIP_Bool oversign;
755  	   int i;
756  	   int nundervars;
757  	   SCIP_Real undermultiplier;
758  	   SCIP_Real overmultiplier;
759  	   SCIP_NLHDLRDATA* nlhdlrdata;
760  	   SCIP_ROWPREP* rowprep;
761  	   SCIP_Real targetunder;
762  	   SCIP_Real scale;
763  	
764  	   assert(conshdlr != NULL);
765  	   assert(nlhdlr != NULL);
766  	   assert(expr != NULL);
767  	   assert(nlhdlrexprdata != NULL);
768  	   assert(rowpreps != NULL);
769  	   *success = FALSE;
770  	   *addedbranchscores = FALSE;
771  	
772  	   /* store and capture the vars of an expression, if the vars are not stored and captured yet */
773  	   if( !nlhdlrexprdata->isstorecapture )
774  	   {
775  	      SCIP_CALL( storeCaptureVars(scip, expr, nlhdlrexprdata) );
776  	   }
777  	
778  	   /* check whether all variables have finite positive bounds, which is necessary for the expression to be a signomial term */
779  	   /* TODO consider allowing 0.0 lower bounds for u variables (and handle gradients at 0.0) */
780  	   isboxsignomial = FALSE;
781  	   checkSignomialBounds(scip, nlhdlrexprdata, &isboxsignomial);
782  	
783  	   if( !isboxsignomial )
784  	      return SCIP_OKAY;
785  	
786  	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
787  	
788  	   /* determine the sign of variables for over- and underestimators, and the multiplier for estimators in the rowprep */
789  	   if( overestimate )
790  	   {
791  	      /* relax t <= x^a, i.e., 0 <= u^f - v^g, overestimate u^f, underestimate v^g */
792  	      oversign = TRUE;
793  	      undersign = FALSE;
794  	      overmultiplier = 1.0;
795  	      undermultiplier = -1.0;
796  	      nundervars = nlhdlrexprdata->nnegvars;
797  	   }
798  	   else
799  	   {
800  	      /* relax x^a <= t, i.e., u^f - v^g <= 0, underestimate u^f, overestimate v^g */
801  	      oversign = FALSE;
802  	      undersign = TRUE;
803  	      overmultiplier = -1.0;
804  	      undermultiplier = 1.0;
805  	      nundervars = nlhdlrexprdata->nposvars;
806  	   }
807  	
808  	   /* compute the value of the overestimator, which is a target value for computing the underestimator efficiently */
809  	   targetunder = 1.0;
810  	   for( i = 0; i < nlhdlrexprdata->nvars; i++ )
811  	   {
812  	      if( nlhdlrexprdata->signs[i] == oversign )
813  	      {
814  	         scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0;
815  	         targetunder *= pow(SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale, nlhdlrexprdata->refexponents[i]);
816  	      }
817  	   }
818  	
819  	#ifdef SCIP_SIGCUT_DEBUG
820  	   SCIP_Real targetover = 1.0;
821  	
822  	   /* print information on estimators */
823  	   SCIP_CALL( printSignomial(scip, expr, nlhdlrexprdata));
824  	   SCIPinfoMessage(scip, NULL, " Auxvalue: %f, targetvalue: %f, %sestimate.", auxvalue, targetvalue, overestimate ? "over" : "under");
825  	
826  	   targetunder = 1.0;
827  	   for( i = 0; i < nlhdlrexprdata->nvars; i++ )
828  	   {
829  	      if( nlhdlrexprdata->signs[i] == undersign )
830  	      {
831  	         scale = i == (nlhdlrexprdata->nvars - 1) ? nlhdlrexprdata->coef : 1.0;
832  	         targetover *= pow(SCIPgetSolVal(scip, sol, nlhdlrexprdata->vars[i]) / scale, nlhdlrexprdata->refexponents[i]);
833  	      }
834  	   }
835  	   SCIPinfoMessage(scip, NULL, " Underestimate: %s, overestimate: %s.",
836  	      undersign ? "positive" : "negative", oversign ? "positive" : "negative");
837  	   SCIPinfoMessage(scip, NULL, " Undervalue (targetover): %f, overvalue (targetunder): %f.", targetover, targetunder);
838  	#endif
839  	
840  	   /* create a rowprep and allocate memory for it */
841  	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
842  	   SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, nlhdlrexprdata->nvars + 1) );
843  	
844  	   /* only underestimate a concave function, if the number of variables is below a given threshold */
845  	   if( nundervars <= nlhdlrdata->maxnundervars )
846  	   {
847  	      /* compute underestimator */
848  	      isspecial = FALSE;
849  	      /* check and compute the special case */
850  	      SCIP_CALL( estimateSpecialPower(scip, nlhdlrexprdata, undersign, undermultiplier, FALSE, sol, rowprep, &isspecial, success) );
851  	      if( !isspecial )
852  	      {
853  	         SCIP_CALL( underEstimatePower(scip, conshdlr, nlhdlr, nlhdlrexprdata, undersign, undermultiplier, sol, targetunder, rowprep, success) );
854  	      }
855  	   }
856  	
857  	   if( *success )
858  	   {
859  	      /* compute overestimator */
860  	      isspecial = FALSE;
861  	      /* check and compute the special case */
862  	      SCIP_CALL( estimateSpecialPower(scip, nlhdlrexprdata, oversign, overmultiplier, TRUE, sol, rowprep, &isspecial, success) );
863  	      if( !isspecial )
864  	      {
865  	         SCIP_CALL( overEstimatePower(scip, nlhdlrexprdata, oversign, overmultiplier, sol, rowprep, success) );
866  	      }
867  	   }
868  	
869  	   if( *success )
870  	   {
871  	      reformRowprep(scip, nlhdlrexprdata, rowprep, nlhdlrdata->mincutscale, success);
872  	      if( *success )
873  	      {
874  	         SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
875  	      }
876  	   }
877  	
878  	   if( !*success )
879  	      SCIPfreeRowprep(scip, &rowprep);
880  	
881  	   return SCIP_OKAY;
882  	}
883  	
884  	/** callback to detect structure in expression tree */
885  	static
886  	SCIP_DECL_NLHDLRDETECT(nlhdlrDetectSignomial)
887  	{ /*lint --e{715}*/
888  	
889  	   assert(expr != NULL);
890  	   assert(enforcing != NULL);
891  	   assert(participating != NULL);
892  	
893  	   /* for now, we do not get active if separation is already (or does not need to be) provided */
894  	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH )
895  	   {
896  	      return SCIP_OKAY;
897  	   }
898  	
899  	   /* check for product expressions with more than one child */
900  	   if( SCIPisExprProduct(scip, expr) && SCIPexprGetNChildren(expr) >= 2 )
901  	   {
902  	      int c;
903  	      int nf = SCIPexprGetNChildren(expr);
904  	      int nvars = nf + 1;
905  	      SCIP_Bool ismultilinear = TRUE;
906  	
907  	      /* create expression data for the nonlinear handler */
908  	      SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) );
909  	      (*nlhdlrexprdata)->nfactors = nf;
910  	      (*nlhdlrexprdata)->nvars = nvars;
911  	
912  	      /* allocate memory for expression data */
913  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->factors, nf) );
914  	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->exponents, nf) );
915  	
916  	      /* get monomial information */
917  	      SCIP_CALL( SCIPgetExprMonomialData(scip, expr, &((*nlhdlrexprdata)->coef), (*nlhdlrexprdata)->exponents, (*nlhdlrexprdata)->factors) );
918  	
919  	      /* skip multilinear terms, since we wouldn't do better than expr_product */
920  	      for( c = 0; c < nf; c++ )
921  	      {
922  	         if( !SCIPisEQ(scip, (*nlhdlrexprdata)->exponents[c], 1.0) )
923  	         {
924  	            ismultilinear = FALSE;
925  	            break;
926  	         }
927  	      }
928  	
929  	      if( ismultilinear )
930  	      {
931  	         /* if multilinear, we free memory of the expression data and do nothing */
932  	         freeExprDataMem(scip, nlhdlrexprdata, TRUE);
933  	         return SCIP_OKAY;
934  	      }
935  	      else
936  	      {
937  	         SCIP_Real normalize;
938  	         SCIP_Real sumlexponents = 0;
939  	         SCIP_Real sumrexponents = 1;
940  	         int nposvars = 0;
941  	
942  	         /* allocate more memory for expression data */
943  	         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->signs, nvars) );
944  	         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->refexponents, nvars) );
945  	         SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &(*nlhdlrexprdata)->vars, nvars) );
946  	         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->intervals, nvars) );
947  	         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->xstar, nvars) );
948  	         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->facetcoefs, nvars) );
949  	         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlhdlrexprdata)->box, 2 * nvars) );
950  	
951  	         (*nlhdlrexprdata)->isstorecapture = FALSE;
952  	
953  	         /* detect more information for the reformulation: we first compute the sum
954  	          * of positive and negative exponents and determine the sign indicators
955  	          */
956  	         for( c = 0; c < nf; c++ )
957  	         {
958  	            /* capture sub expressions */
959  	            SCIPcaptureExpr((*nlhdlrexprdata)->factors[c]);
960  	            if( (*nlhdlrexprdata)->exponents[c] > 0.0 )
961  	            {
962  	               /* add a positive exponent */
963  	               sumlexponents += (*nlhdlrexprdata)->exponents[c];
964  	               (*nlhdlrexprdata)->signs[c] = TRUE;
965  	               nposvars++;
966  	            }
967  	            else
968  	            {
969  	               /* subtract a negative exponent */
970  	               sumrexponents -= (*nlhdlrexprdata)->exponents[c];
971  	               (*nlhdlrexprdata)->signs[c] = FALSE;
972  	            }
973  	            /* set null to working variables, meaning that they are not stored yet */
974  	         }
975  	         (*nlhdlrexprdata)->signs[nf] = FALSE;
976  	         (*nlhdlrexprdata)->nposvars = nposvars;
977  	         (*nlhdlrexprdata)->nnegvars = nf - nposvars + 1;
978  	
979  	         /* compute the normalization constant */
980  	         normalize = MAX(sumlexponents, sumrexponents);
981  	         /* normalize positive and negative exponents */
982  	         for( c = 0; c < nf; c++ )
983  	         {
984  	            if( (*nlhdlrexprdata)->signs[c] )
985  	               (*nlhdlrexprdata)->refexponents[c] = (*nlhdlrexprdata)->exponents[c] / normalize;
986  	            else
987  	               (*nlhdlrexprdata)->refexponents[c] = -(*nlhdlrexprdata)->exponents[c] / normalize;
988  	         }
989  	         (*nlhdlrexprdata)->refexponents[nf] = 1.0 / normalize;
990  	
991  	         /* tell children that we will use their auxvar and use its activity for estimation */
992  	         for( c = 0; c < nf; c++ )
993  	         {
994  	            SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, (*nlhdlrexprdata)->factors[c], TRUE, FALSE, TRUE, TRUE) );
995  	         }
996  	      }
997  	   }
998  	
999  	   if( *nlhdlrexprdata != NULL )
1000 	   {
1001 	      *participating = SCIP_NLHDLR_METHOD_SEPABOTH;
1002 	   }
1003 	
1004 	#ifdef SCIP_SIGCUT_DEBUG
1005 	   if( *participating )
1006 	   {
1007 	      SCIPdebugMsg(scip, "scip depth: %d, step: %d, expr pointer: %p, expr data pointer: %p, detected expr: total vars (exps) %d ",
1008 	         SCIPgetSubscipDepth(scip), SCIPgetStage(scip), (void *)expr, (void *)nlhdlrexprdata, (*nlhdlrexprdata)->nfactors);
1009 	      SCIPprintExpr(scip, expr, NULL);
1010 	      SCIPinfoMessage(scip, NULL, " participating: %d\n", *participating);
1011 	   }
1012 	#endif
1013 	
1014 	   return SCIP_OKAY;
1015 	}
1016 	
1017 	/** auxiliary evaluation callback of nonlinear handler */
1018 	static SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxSignomial)
1019 	{ /*lint --e{715}*/
1020 	   int c;
1021 	   SCIP_Real val;
1022 	   SCIP_VAR* var;
1023 	
1024 	   *auxvalue = nlhdlrexprdata->coef;
1025 	   for( c = 0; c < nlhdlrexprdata->nfactors; ++c )
1026 	   {
1027 	      assert(nlhdlrexprdata->factors[c] != NULL);
1028 	
1029 	      var = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->factors[c]);
1030 	      assert(var != NULL);
1031 	
1032 	      val = SCIPgetSolVal(scip, sol, var);
1033 	      if( SCIPisPositive(scip, val) )
1034 	      {
1035 	         *auxvalue *= pow(val, nlhdlrexprdata->exponents[c]);
1036 	      }
1037 	      else
1038 	      {
1039 	         *auxvalue = SCIP_INVALID;
1040 	         break;
1041 	      }
1042 	   }
1043 	
1044 	   return SCIP_OKAY;
1045 	}
1046 	
1047 	/** nonlinear handler copy callback */
1048 	static
1049 	SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrSignomial)
1050 	{ /*lint --e{715}*/
1051 	   assert(targetscip != NULL);
1052 	   assert(sourcenlhdlr != NULL);
1053 	   assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
1054 	
1055 	   SCIP_CALL( SCIPincludeNlhdlrSignomial(targetscip) );
1056 	
1057 	   return SCIP_OKAY;
1058 	}
1059 	
1060 	/** callback to free data of handler */
1061 	static
1062 	SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrDataSignomial)
1063 	{ /*lint --e{715}*/
1064 	   assert(nlhdlrdata != NULL);
1065 	
1066 	   SCIPfreeBlockMemory(scip, nlhdlrdata);
1067 	
1068 	   return SCIP_OKAY;
1069 	}
1070 	
1071 	/** callback to free expression specific data */
1072 	static
1073 	SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataSignomial)
1074 	{ /*lint --e{715}*/
1075 	   int c;
1076 	
1077 	   /* release expressions */
1078 	   for( c = 0; c < (*nlhdlrexprdata)->nfactors; c++ )
1079 	   {
1080 	      SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->factors[c]) );
1081 	   }
1082 	
1083 	   /* release variables */
1084 	   if( (*nlhdlrexprdata)->isstorecapture )
1085 	   {
1086 	      for( c = 0; c < (*nlhdlrexprdata)->nvars; c++ )
1087 	      {
1088 	         if( (*nlhdlrexprdata)->vars[c] != NULL )
1089 	         {
1090 	            SCIP_CALL( SCIPreleaseVar(scip, &(*nlhdlrexprdata)->vars[c]) );
1091 	         }
1092 	      }
1093 	   }
1094 	
1095 	   /* free memory */
1096 	   freeExprDataMem(scip, nlhdlrexprdata, FALSE);
1097 	
1098 	   return SCIP_OKAY;
1099 	}
1100 	
1101 	
1102 	/*
1103 	 * nonlinear handler specific interface methods
1104 	 */
1105 	
1106 	/** includes signomial nonlinear handler in nonlinear constraint handler */
1107 	SCIP_RETCODE
1108 	SCIPincludeNlhdlrSignomial(
1109 	   SCIP*                 scip                /**< SCIP data structure */
1110 	   )
1111 	{
1112 	   SCIP_NLHDLRDATA* nlhdlrdata;
1113 	   SCIP_NLHDLR* nlhdlr;
1114 	
1115 	   assert(scip != NULL);
1116 	
1117 	   /* create nonlinear handler specific data */
1118 	   SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata));
1119 	   BMSclearMemory(nlhdlrdata);
1120 	
1121 	   SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY,
1122 	      NLHDLR_ENFOPRIORITY, nlhdlrDetectSignomial, nlhdlrEvalauxSignomial, nlhdlrdata) );
1123 	   assert(nlhdlr != NULL);
1124 	
1125 	   SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrSignomial);
1126 	   SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrDataSignomial);
1127 	   SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataSignomial);
1128 	   SCIPnlhdlrSetSepa(nlhdlr, NULL, NULL, nlhdlrEstimateSignomial, NULL);
1129 	
1130 	   /* parameters */
1131 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxnundervars",
1132 	         "maximum number of variables when underestimating a concave power function",
1133 	         &nlhdlrdata->maxnundervars, TRUE, NLHDLR_MAXNUNDERVARS, 2, 14, NULL, NULL) );
1134 	   SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutscale",
1135 	         "minimum scale factor when scaling a cut",
1136 	         &nlhdlrdata->mincutscale, TRUE, NLHDLR_MINCUTSCALE, 1e-6, 1e6, NULL, NULL) );
1137 	
1138 	   return SCIP_OKAY;
1139 	}
1140