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   expr_sum.c
26   	 * @ingroup DEFPLUGINS_EXPR
27   	 * @brief  sum expression handler
28   	 * @author Stefan Vigerske
29   	 * @author Benjamin Mueller
30   	 * @author Felipe Serrano
31   	 */
32   	
33   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34   	
35   	#include <string.h>
36   	#include <stddef.h>
37   	
38   	#include "scip/expr_sum.h"
39   	#include "scip/expr_value.h"
40   	#include "scip/expr_product.h"
41   	#include "scip/expr_exp.h"
42   	#include "scip/expr_pow.h"
43   	#include "symmetry/struct_symmetry.h"
44   	
45   	#define EXPRHDLR_NAME         "sum"
46   	#define EXPRHDLR_DESC         "summation with coefficients and a constant"
47   	#define EXPRHDLR_PRECEDENCE   40000
48   	#define EXPRHDLR_HASHKEY      SCIPcalcFibHash(47161.0)
49   	
50   	/** macro to activate/deactivate debugging information of simplify method */
51   	/*lint -emacro(774,debugSimplify) */
52   	#ifdef SIMPLIFY_DEBUG
53   	#define debugSimplify                   printf
54   	#else
55   	#define debugSimplify                   while( FALSE ) printf
56   	#endif
57   	
58   	/*
59   	 * Data structures
60   	 */
61   	
62   	/** expression data */
63   	struct SCIP_ExprData
64   	{
65   	   SCIP_Real             constant;           /**< constant coefficient */
66   	   SCIP_Real*            coefficients;       /**< coefficients of children */
67   	   int                   coefssize;          /**< size of the coefficients array */
68   	};
69   	
70   	/*
71   	 * Local methods
72   	 */
73   	
74   	/** creates expression data */
75   	static
76   	SCIP_RETCODE createData(
77   	   SCIP*                 scip,               /**< SCIP data structure */
78   	   SCIP_EXPRDATA**       exprdata,           /**< pointer where to store expression data */
79   	   int                   ncoefficients,      /**< number of coefficients (i.e., number of children) */
80   	   SCIP_Real*            coefficients,       /**< array with coefficients for all children (or NULL if all 1.0) */
81   	   SCIP_Real             constant            /**< constant term of sum */
82   	   )
83   	{
84   	   assert(exprdata != NULL);
85   	   assert(ncoefficients >= 0);
86   	
87   	   SCIP_CALL( SCIPallocBlockMemory(scip, exprdata) );
88   	
89   	   if( coefficients != NULL )
90   	   {
91   	      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*exprdata)->coefficients, coefficients, ncoefficients) );
92   	   }
93   	   else
94   	   {
95   	      int i;
96   	
97   	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*exprdata)->coefficients, ncoefficients) );
98   	      for( i = 0; i < ncoefficients; ++i )
99   	         (*exprdata)->coefficients[i] = 1.0;
100  	   }
101  	
102  	   (*exprdata)->coefssize = ncoefficients;
103  	   (*exprdata)->constant  = constant;
104  	
105  	   return SCIP_OKAY;
106  	}
107  	
108  	/** simplifies the `idx`-th child of the sum expression `duplicate` in order for it to be able to be a child of a simplified sum
109  	 *
110  	 * for example, this means that the `idx`-th child cannot be itself a sum
111  	 * if it is, we have to flatten it, i.e., take all its children and make them children of `duplicate`
112  	 */
113  	static
114  	SCIP_RETCODE simplifyTerm(
115  	   SCIP*                 scip,               /**< SCIP data structure */
116  	   SCIP_EXPR*            duplicate,          /**< expression to be simplified */
117  	   int                   idx,                /**< idx of children to be simplified */
118  	   SCIP_Bool*            changed,            /**< pointer to store if some term actually got simplified */
119  	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
120  	   void*                 ownercreatedata     /**< data to pass to ownercreate */
121  	   )
122  	{
123  	   SCIP_EXPR** children;
124  	   SCIP_EXPR* expr;
125  	   SCIP_Real* coefs;
126  	   SCIP_Real constant ;
127  	   SCIP_Real coef;
128  	
129  	   assert(duplicate != NULL);
130  	   assert(idx >= 0);
131  	   assert(idx < SCIPexprGetNChildren(duplicate));
132  	   assert(changed != NULL);
133  	
134  	   children  = SCIPexprGetChildren(duplicate);
135  	   coefs     = SCIPgetCoefsExprSum(duplicate);
136  	   constant  = SCIPgetConstantExprSum(duplicate);
137  	
138  	   coef = coefs[idx];
139  	   expr = children[idx];
140  	   assert(expr != NULL);
141  	
142  	   /* enforces SS3 */
143  	   if( SCIPisExprValue(scip, expr) )
144  	   {
145  	      *changed = TRUE;
146  	      constant += coef * SCIPgetValueExprValue(expr);
147  	      SCIPsetConstantExprSum(duplicate, constant);
148  	
149  	      /* TODO: remove child? */
150  	      coefs[idx] = 0.0;
151  	
152  	      return SCIP_OKAY;
153  	   }
154  	
155  	   /* enforces SS2 */
156  	   if( SCIPisExprSum(scip, expr) )
157  	   {
158  	      *changed = TRUE;
159  	
160  	      /* pass constant to parent */
161  	      constant += coef * SCIPgetConstantExprSum(expr);
162  	      SCIPsetConstantExprSum(duplicate, constant);
163  	
164  	      /* append all children of expr on parent except the first one */
165  	      if( SCIPexprGetNChildren(expr) > 1 )
166  	      {
167  	         int i;
168  	
169  	         for( i = 1; i < SCIPexprGetNChildren(expr); ++i )
170  	         {
171  	            assert(!SCIPisExprSum(scip, SCIPexprGetChildren(expr)[i]));
172  	            SCIP_CALL( SCIPappendExprSumExpr(scip, duplicate, SCIPexprGetChildren(expr)[i],
173  	                  coef * SCIPgetCoefsExprSum(expr)[i]) );
174  	         }
175  	      }
176  	
177  	      /* replace expr with first child; need to get data again since it might be re-allocated */
178  	      assert(!SCIPisExprSum(scip, SCIPexprGetChildren(expr)[0]));
179  	
180  	      coefs = SCIPgetCoefsExprSum(duplicate);
181  	
182  	      coefs[idx] = coef * SCIPgetCoefsExprSum(expr)[0];
183  	      SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, SCIPexprGetChildren(expr)[0]) );
184  	
185  	      return SCIP_OKAY;
186  	   }
187  	
188  	   /* enforce SS9 */
189  	   if( REALABS(coef) != 1.0 && SCIPisExprProduct(scip, expr) )
190  	   {
191  	      SCIP_EXPR* expchild = NULL;
192  	      int i;
193  	
194  	      for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
195  	      {
196  	         SCIP_EXPR* child = SCIPexprGetChildren(expr)[i];
197  	         assert(child != NULL);
198  	
199  	         if( SCIPisExprExp(scip, child) )
200  	         {
201  	            expchild = child;
202  	            break;
203  	         }
204  	      }
205  	
206  	      /* coef != +- 1, term is product and one factor is an exponential -> enforce SS9 */
207  	      if( expchild != NULL )
208  	      {
209  	         SCIP_EXPR* sum;
210  	         SCIP_EXPR* prod;
211  	         SCIP_EXPR* simplifiedprod;
212  	         SCIP_EXPR* simplifiedsum;
213  	         SCIP_EXPR* exponential;
214  	         SCIP_EXPR* simplifiedexp;
215  	         SCIP_Real expconstant;
216  	
217  	         /* inform that expression will change */
218  	         *changed = TRUE;
219  	
220  	         /* compute expchild's coefficient as +- 1.0 * exp(log(abs(coef))) */
221  	         if( coef > 0.0 )
222  	         {
223  	            expconstant = log(coef);
224  	            coefs[idx] = 1.0;
225  	         }
226  	         else
227  	         {
228  	            expconstant = log(-coef);
229  	            coefs[idx] = -1.0;
230  	         }
231  	
232  	         /* add constant to exponential's child */
233  	         SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, SCIPexprGetChildren(expchild), NULL, expconstant, ownercreate,
234  	               ownercreatedata) );
235  	
236  	         /* simplify sum */
237  	         SCIP_CALL( SCIPcallExprSimplify(scip, sum, &simplifiedsum, ownercreate, ownercreatedata) );
238  	         SCIP_CALL( SCIPreleaseExpr(scip, &sum) );
239  	
240  	         /* create exponential with new child */
241  	         SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedsum, ownercreate, ownercreatedata) );
242  	         SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedsum) );
243  	
244  	         /* simplify exponential */
245  	         SCIP_CALL( SCIPcallExprSimplify(scip, exponential, &simplifiedexp, ownercreate, ownercreatedata) );
246  	         SCIP_CALL( SCIPreleaseExpr(scip, &exponential) );
247  	
248  	         /* create product with new child */
249  	         SCIP_CALL( SCIPcreateExprProduct(scip, &prod, 0, NULL, 1.0, ownercreate, ownercreatedata) );
250  	
251  	         for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
252  	         {
253  	            if( SCIPexprGetChildren(expr)[i] == expchild )
254  	            {
255  	               SCIP_CALL( SCIPappendExprChild(scip, prod, simplifiedexp) );
256  	            }
257  	            else
258  	            {
259  	               SCIP_CALL( SCIPappendExprChild(scip, prod, SCIPexprGetChildren(expr)[i]) );
260  	            }
261  	         }
262  	         SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedexp) );
263  	
264  	         /* simplify product */
265  	         SCIP_CALL( SCIPcallExprSimplify(scip, prod, &simplifiedprod, ownercreate, ownercreatedata) );
266  	         SCIP_CALL( SCIPreleaseExpr(scip, &prod) );
267  	
268  	         /* replace current child with simplified product */
269  	         SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, simplifiedprod) );
270  	         SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedprod) );
271  	
272  	         /* since the simplified product can be a sum ( exp(-1)*exp(log(x+y)+1) -> x+y ),
273  	          * we call the function we are in again
274  	          * this is no endless recursion, since the coef is now +- 1
275  	          */
276  	         SCIP_CALL( simplifyTerm(scip, duplicate, idx, changed, ownercreate, ownercreatedata) );
277  	
278  	         return SCIP_OKAY;
279  	      }
280  	   }
281  	
282  	   /* enforce SS10 */
283  	   if( REALABS(coef) != 1.0 && SCIPisExprExp(scip, expr) )
284  	   {
285  	      /* coef != +- 1, term is exponential -> enforce SS10 by moving |coef| into argument of exponential */
286  	
287  	      SCIP_EXPR* sum;
288  	      SCIP_EXPR* simplifiedsum;
289  	      SCIP_EXPR* exponential;
290  	      SCIP_EXPR* simplifiedexp;
291  	      SCIP_Real expconstant;
292  	
293  	      /* inform that expression will change */
294  	      *changed = TRUE;
295  	
296  	      /* compute expchild's coefficient as +- 1.0 * exp(log(abs(coef))) */
297  	      if( coef > 0.0 )
298  	      {
299  	         expconstant = log(coef);
300  	         coefs[idx] = 1.0;
301  	      }
302  	      else
303  	      {
304  	         expconstant = log(-coef);
305  	         coefs[idx] = -1.0;
306  	      }
307  	
308  	      /* add constant to exponential's child */
309  	      SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, SCIPexprGetChildren(expr), NULL, expconstant, ownercreate,
310  	            ownercreatedata) );  /* expconstant+expchild */
311  	
312  	      /* simplify sum */
313  	      SCIP_CALL( SCIPcallExprSimplify(scip, sum, &simplifiedsum, ownercreate, ownercreatedata) );
314  	      SCIP_CALL( SCIPreleaseExpr(scip, &sum) );
315  	
316  	      /* create exponential with new child */
317  	      SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedsum, ownercreate, ownercreatedata) );
318  	      SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedsum) );
319  	
320  	      /* simplify exponential */
321  	      SCIP_CALL( SCIPcallExprSimplify(scip, exponential, &simplifiedexp, ownercreate, ownercreatedata) );
322  	      SCIP_CALL( SCIPreleaseExpr(scip, &exponential) );
323  	
324  	      /* replace current child with simplified exponential */
325  	      SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, simplifiedexp) );
326  	      SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedexp) );
327  	
328  	      return SCIP_OKAY;
329  	   }
330  	
331  	   /* other types of (simplified) expressions can be a child of a simplified sum */
332  	   assert(!SCIPisExprSum(scip, expr));
333  	   assert(!SCIPisExprValue(scip, expr));
334  	
335  	   return SCIP_OKAY;
336  	}
337  	
338  	/** helper struct for expressions sort */
339  	typedef struct
340  	{
341  	   SCIP*                 scip;               /**< SCIP data structure */
342  	   SCIP_EXPR**           exprs;              /**< expressions */
343  	} SORTEXPRDATA;
344  	
345  	static
346  	SCIP_DECL_SORTINDCOMP(sortExprComp)
347  	{
348  	   SORTEXPRDATA* data = (SORTEXPRDATA*)dataptr;
349  	
350  	   return SCIPcompareExpr(data->scip, data->exprs[ind1], data->exprs[ind2]);
351  	}
352  	
353  	/*
354  	 * Callback methods of expression handler
355  	 */
356  	
357  	/** simplifies a sum expression
358  	 *
359  	 * goes through each child and simplifies it; then sorts the simplified children; then sum the children that are equal;
360  	 * finally creates a sum expression with all the children that do not have a 0 coefficient and post-process so that SS6
361  	 * and SS7 are satisfied
362  	 */
363  	static
364  	SCIP_DECL_EXPRSIMPLIFY(simplifySum)
365  	{  /*lint --e{715}*/
366  	   SCIP_EXPR** children;
367  	   SCIP_EXPR* duplicate = NULL;
368  	   SCIP_EXPR** newchildren = NULL;
369  	   SCIP_Real* newcoefs = NULL;
370  	   int nnewchildren;
371  	   SCIP_Real newconstant;
372  	   SCIP_Real* coefs;
373  	   int i;
374  	   int nchildren;
375  	   SCIP_Bool changed;
376  	   SORTEXPRDATA sortdata;
377  	   int* order = NULL;
378  	
379  	   assert(expr != NULL);
380  	   assert(simplifiedexpr != NULL);
381  	   assert(SCIPexprGetHdlr(expr) == SCIPgetExprhdlrSum(scip));
382  	
383  	   changed = FALSE;
384  	
385  	   /* TODO: maybe have a flag to know if it is simplified ? */
386  	   /* TODO: can we do this with a shallow duplicate + copy of children pointer? currently simplifyTerm may modify children,
387  	    * so one would need to be careful
388  	    */
389  	   SCIP_CALL( SCIPduplicateExpr(scip, expr, &duplicate, NULL, NULL, ownercreate, ownercreatedata) );
390  	   assert(duplicate != NULL);
391  	
392  	   nchildren = SCIPexprGetNChildren(duplicate);
393  	   for( i = 0; i < nchildren; i++ )
394  	   {
395  	      /* enforces SS8 TODO: remove child? */
396  	      /* we have to ask for the coefs everytime, since it might get realloced in simpifyTerm */
397  	      if( SCIPgetCoefsExprSum(duplicate)[i] == 0.0 )
398  	      {
399  	         changed = TRUE;
400  	         continue;
401  	      }
402  	
403  	      /* enforces SS2, SS3, SS9, and SS10 */
404  	      SCIP_CALL( simplifyTerm(scip, duplicate, i, &changed, ownercreate, ownercreatedata) );
405  	   }
406  	
407  	   /* simplifyTerm can add new children to duplicate and realloc them; so get them again */
408  	   nchildren = SCIPexprGetNChildren(duplicate);
409  	   children  = SCIPexprGetChildren(duplicate);
410  	   coefs     = SCIPgetCoefsExprSum(duplicate);
411  	
412  	   /* treat zero term case */
413  	   if( nchildren == 0 )
414  	   {
415  	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, SCIPgetConstantExprSum(duplicate), ownercreate, ownercreatedata) );
416  	      goto CLEANUP;
417  	   }
418  	
419  	   /* treat one term case */
420  	   if( nchildren == 1 )
421  	   {
422  	      if( coefs[0] == 0.0 )
423  	      {
424  	         SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, SCIPgetConstantExprSum(duplicate), ownercreate, ownercreatedata) );
425  	         goto CLEANUP;
426  	      }
427  	
428  	      if( coefs[0] == 1.0 && SCIPgetConstantExprSum(duplicate) == 0.0 )
429  	         *simplifiedexpr = children[0]; /* SS7 */
430  	      else
431  	         *simplifiedexpr = changed ? duplicate : expr;
432  	
433  	      SCIPcaptureExpr(*simplifiedexpr);
434  	
435  	      goto CLEANUP;
436  	   }
437  	
438  	   /* enforces SS5: sort children */
439  	   SCIP_CALL( SCIPallocBufferArray(scip, &order, nchildren) );
440  	   for( i = 0; i < nchildren; i++ )
441  	      order[i] = i;
442  	   sortdata.scip = scip;
443  	   sortdata.exprs = children;
444  	   SCIPsortInd(order, sortExprComp, (void*)&sortdata, nchildren);
445  	
446  	   /* create sorted variant of children and coefs */
447  	   SCIP_CALL( SCIPallocBufferArray(scip, &newchildren, nchildren) );
448  	   SCIP_CALL( SCIPallocBufferArray(scip, &newcoefs, nchildren) );
449  	   for( i = 0; i < nchildren; ++i )
450  	   {
451  	      newchildren[i] = children[order[i]];
452  	      newcoefs[i] = coefs[order[i]];
453  	      if( order[i] != i )
454  	         changed = TRUE;
455  	   }
456  	
457  	   /* post-process */
458  	
459  	   /* enforces SS4 */
460  	   nnewchildren = 0;
461  	   for( i = 0; i < nchildren; i++ )
462  	   {
463  	      /* eliminate zero-coefficients */
464  	      if( newcoefs[i] == 0.0 )
465  	      {
466  	         changed = TRUE;
467  	         continue;
468  	      }
469  	
470  	      /* sum equal expressions */
471  	      if( i < nchildren-1 && SCIPcompareExpr(scip, newchildren[i], newchildren[i+1]) == 0 )
472  	      {
473  	         changed = TRUE;
474  	         /* if we substract two almost equal not-so-small numbers, then set new coefficient to 0.0
475  	          * instead of some tiny value that is likely the result of some random round-off error
476  	          * E.g., on instance ex1221, we have x1^2 + b3 = 1.25.
477  	          *   Probing finds an aggregation x1 = 1.11803 - 0.618034 b3.
478  	          *   Simplify would then produce 1.25 + 1e-16 x1 = 1.25.
479  	          */
480  	         if( SCIPisEQ(scip, newcoefs[i], -newcoefs[i+1]) && REALABS(newcoefs[i]) >= 1.0 )
481  	            newcoefs[i+1] = 0.0;
482  	         else
483  	            newcoefs[i+1] += newcoefs[i];
484  	         continue;
485  	      }
486  	
487  	      /* move i-th child to new position */
488  	      newchildren[nnewchildren] = newchildren[i];
489  	      newcoefs[nnewchildren] = newcoefs[i];
490  	      nnewchildren++;
491  	   }
492  	
493  	   /* build sum expression from finalchildren and post-simplify */
494  	   newconstant = SCIPgetConstantExprSum(duplicate);
495  	
496  	   debugSimplify("what to do? finalchildren has length %d\n", nnewchildren); /*lint !e506 !e681*/
497  	
498  	   /* enforces SS6: if they are no children, return value */
499  	   if( nnewchildren == 0 )
500  	   {
501  	      debugSimplify("[sum] got empty list, return value %g\n", newconstant); /*lint !e506 !e681*/
502  	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, newconstant, ownercreate, ownercreatedata) );
503  	
504  	      goto CLEANUP;
505  	   }
506  	
507  	   /* enforces SS7: if list consists of one expr with coef 1.0 and constant is 0, return that expr */
508  	   if( nnewchildren == 1 && newcoefs[0] == 1.0 && newconstant == 0.0 )
509  	   {
510  	      *simplifiedexpr = newchildren[0];
511  	      SCIPcaptureExpr(*simplifiedexpr);
512  	
513  	      goto CLEANUP;
514  	   }
515  	
516  	   /* build sum expression from children */
517  	   if( changed )
518  	   {
519  	      SCIP_CALL( SCIPcreateExprSum(scip, simplifiedexpr, nnewchildren, newchildren, newcoefs, newconstant,
520  	            ownercreate, ownercreatedata) );
521  	
522  	      goto CLEANUP;
523  	   }
524  	
525  	   *simplifiedexpr = expr;
526  	
527  	   /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
528  	   SCIPcaptureExpr(*simplifiedexpr);
529  	
530  	   /* free memory */
531  	 CLEANUP:
532  	   SCIPfreeBufferArrayNull(scip, &newcoefs);
533  	   SCIPfreeBufferArrayNull(scip, &newchildren);
534  	   SCIPfreeBufferArrayNull(scip, &order);
535  	   SCIP_CALL( SCIPreleaseExpr(scip, &duplicate) );
536  	
537  	   return SCIP_OKAY;
538  	}
539  	
540  	/** expression callback to get information for symmetry detection */
541  	static
542  	SCIP_DECL_EXPRGETSYMDATA(getSymDataSum)
543  	{  /*lint --e{715}*/
544  	   SCIP_EXPRDATA* exprdata;
545  	   int i;
546  	
547  	   assert(scip != NULL);
548  	   assert(expr != NULL);
549  	
550  	   exprdata = SCIPexprGetData(expr);
551  	   assert(exprdata != NULL);
552  	
553  	   SCIP_CALL( SCIPallocBlockMemory(scip, symdata) );
554  	
555  	   (*symdata)->nconstants = 1;
556  	   (*symdata)->ncoefficients = exprdata->coefssize;
557  	
558  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->constants, 1) );
559  	   (*symdata)->constants[0] = exprdata->constant;
560  	
561  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->coefficients, exprdata->coefssize) );
562  	   for( i = 0; i < exprdata->coefssize; ++i )
563  	      (*symdata)->coefficients[i] = exprdata->coefficients[i];
564  	
565  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->children, exprdata->coefssize) );
566  	   for( i = 0; i < exprdata->coefssize; ++i )
567  	      (*symdata)->children[i] = SCIPexprGetChildren(expr)[i];
568  	
569  	   return SCIP_OKAY;
570  	}
571  	
572  	/** compares two sum expressions
573  	 *
574  	 *  The order of two sum expressions is a lexicographical order on the terms.
575  	 *
576  	 *  Starting from the *last*, we find the first child where they differ, say, the i-th.
577  	 *  Then u < v <=> u_i < v_i.
578  	 *  If there are no such children and they have different number of children, then u < v <=> nchildren(u) < nchildren(v).
579  	 *  If there are no such children and they have the same number of children, then u < v <=> const(u) < const(v).
580  	 *  Otherwise, they are the same.
581  	 *
582  	 *  Note: we are assuming expression are simplified, so within u, we have u_1 < u_2, etc
583  	 *
584  	 *  Example: y + z < x + y + z, 2*x + 3*y < 3*x + 3*y
585  	 */
586  	static
587  	SCIP_DECL_EXPRCOMPARE(compareSum)
588  	{  /*lint --e{715}*/
589  	   SCIP_Real const1;
590  	   SCIP_Real* coefs1;
591  	   SCIP_EXPR** children1;
592  	   int nchildren1;
593  	   SCIP_Real const2;
594  	   SCIP_Real* coefs2;
595  	   SCIP_EXPR** children2;
596  	   int nchildren2;
597  	   int compareresult;
598  	   int i;
599  	   int j;
600  	
601  	   nchildren1 = SCIPexprGetNChildren(expr1);
602  	   nchildren2 = SCIPexprGetNChildren(expr2);
603  	   children1 = SCIPexprGetChildren(expr1);
604  	   children2 = SCIPexprGetChildren(expr2);
605  	   coefs1 = SCIPgetCoefsExprSum(expr1);
606  	   coefs2 = SCIPgetCoefsExprSum(expr2);
607  	   const1 = SCIPgetConstantExprSum(expr1);
608  	   const2 = SCIPgetConstantExprSum(expr2);
609  	
610  	   for( i = nchildren1 - 1, j = nchildren2 - 1; i >= 0 && j >= 0; --i, --j )
611  	   {
612  	      compareresult = SCIPcompareExpr(scip, children1[i], children2[j]);
613  	      if( compareresult != 0 )
614  	         return compareresult;
615  	      else
616  	      {
617  	         /* expressions are equal, compare coefficient */
618  	         if( (coefs1 ? coefs1[i] : 1.0) < (coefs2 ? coefs2[j] : 1.0) )
619  	            return -1;
620  	         if( (coefs1 ? coefs1[i] : 1.0) > (coefs2 ? coefs2[j] : 1.0) )
621  	            return 1;
622  	
623  	         /* coefficients are equal, continue */
624  	      }
625  	   }
626  	
627  	   /* all children of one expression are children of the other expression, use number of children as a tie-breaker */
628  	   if( i < j )
629  	   {
630  	      assert(i == -1);
631  	      /* expr1 has less elements, hence expr1 < expr2 */
632  	      return -1;
633  	   }
634  	   if( i > j )
635  	   {
636  	      assert(j == -1);
637  	      /* expr1 has more elements, hence expr1 > expr2 */
638  	      return 1;
639  	   }
640  	
641  	   /* everything is equal, use constant/coefficient as tie-breaker */
642  	   assert(i == -1 && j == -1);
643  	   if( const1 < const2 )
644  	      return -1;
645  	   if( const1 > const2 )
646  	      return 1;
647  	
648  	   /* they are equal */
649  	   return 0;
650  	}
651  	
652  	/** expression handler copy callback */
653  	static
654  	SCIP_DECL_EXPRCOPYHDLR(copyhdlrSum)
655  	{  /*lint --e{715}*/
656  	   SCIP_CALL( SCIPincludeExprhdlrSum(scip) );
657  	
658  	   return SCIP_OKAY;
659  	}
660  	
661  	/** expression data copy callback */
662  	static
663  	SCIP_DECL_EXPRCOPYDATA(copydataSum)
664  	{  /*lint --e{715}*/
665  	   SCIP_EXPRDATA* sourceexprdata;
666  	
667  	   assert(targetexprdata != NULL);
668  	   assert(sourceexpr != NULL);
669  	
670  	   sourceexprdata = SCIPexprGetData(sourceexpr);
671  	   assert(sourceexprdata != NULL);
672  	
673  	   SCIP_CALL( createData(targetscip, targetexprdata, SCIPexprGetNChildren(sourceexpr),
674  	            sourceexprdata->coefficients, sourceexprdata->constant) );
675  	
676  	   return SCIP_OKAY;
677  	}
678  	
679  	/** expression data free callback */
680  	static
681  	SCIP_DECL_EXPRFREEDATA(freedataSum)
682  	{  /*lint --e{715}*/
683  	   SCIP_EXPRDATA* exprdata;
684  	
685  	   assert(expr != NULL);
686  	
687  	   exprdata = SCIPexprGetData(expr);
688  	   assert(exprdata != NULL);
689  	
690  	   SCIPfreeBlockMemoryArray(scip, &(exprdata->coefficients), exprdata->coefssize);
691  	   SCIPfreeBlockMemory(scip, &exprdata);
692  	
693  	   SCIPexprSetData(expr, NULL);
694  	
695  	   return SCIP_OKAY;
696  	}
697  	
698  	/** expression print callback */
699  	static
700  	SCIP_DECL_EXPRPRINT(printSum)
701  	{  /*lint --e{715}*/
702  	   SCIP_EXPRDATA* exprdata;
703  	
704  	   assert(expr != NULL);
705  	
706  	   exprdata = SCIPexprGetData(expr);
707  	   assert(exprdata != NULL);
708  	
709  	   /**! [SnippetExprPrintSum] */
710  	   switch( stage )
711  	   {
712  	      case SCIP_EXPRITER_ENTEREXPR :
713  	      {
714  	         /* print opening parenthesis, if necessary */
715  	         if( EXPRHDLR_PRECEDENCE <= parentprecedence )
716  	         {
717  	            SCIPinfoMessage(scip, file, "(");
718  	         }
719  	
720  	         /* print constant, if nonzero */
721  	         if( exprdata->constant != 0.0 )
722  	         {
723  	            SCIPinfoMessage(scip, file, "%g", exprdata->constant);
724  	         }
725  	         break;
726  	      }
727  	
728  	      case SCIP_EXPRITER_VISITINGCHILD :
729  	      {
730  	         SCIP_Real coef;
731  	
732  	         coef = exprdata->coefficients[currentchild];
733  	
734  	         /* print coefficient, if necessary */
735  	         if( coef == 1.0 )
736  	         {
737  	            /* if coefficient is 1.0, then print only "+" if not the first term */
738  	            if( exprdata->constant != 0.0 || currentchild > 0 )
739  	            {
740  	               SCIPinfoMessage(scip, file, "+");
741  	            }
742  	         }
743  	         else if( coef == -1.0 )
744  	         {
745  	            /* if coefficient is -1.0, then print only "-" */
746  	            SCIPinfoMessage(scip, file, "-");
747  	         }
748  	         else
749  	         {
750  	            /* force "+" sign on positive coefficient if not the first term */
751  	            SCIPinfoMessage(scip, file, (exprdata->constant != 0.0 || currentchild > 0) ? "%+g*" : "%g*", coef);
752  	         }
753  	
754  	         break;
755  	      }
756  	
757  	      case SCIP_EXPRITER_LEAVEEXPR :
758  	      {
759  	         /* print closing parenthesis, if necessary */
760  	         if( EXPRHDLR_PRECEDENCE <= parentprecedence )
761  	         {
762  	            SCIPinfoMessage(scip, file, ")");
763  	         }
764  	         break;
765  	      }
766  	
767  	      case SCIP_EXPRITER_VISITEDCHILD :
768  	      default: ;
769  	   }
770  	   /**! [SnippetExprPrintSum] */
771  	
772  	   return SCIP_OKAY;
773  	}
774  	
775  	/** expression point evaluation callback */
776  	static
777  	SCIP_DECL_EXPREVAL(evalSum)
778  	{  /*lint --e{715}*/
779  	   SCIP_EXPRDATA* exprdata;
780  	   int c;
781  	
782  	   assert(expr != NULL);
783  	
784  	   exprdata = SCIPexprGetData(expr);
785  	   assert(exprdata != NULL);
786  	
787  	   /**! [SnippetExprEvalSum] */
788  	   *val = exprdata->constant;
789  	   for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
790  	   {
791  	      assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[c]) != SCIP_INVALID); /*lint !e777*/
792  	
793  	      *val += exprdata->coefficients[c] * SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[c]);
794  	   }
795  	   /**! [SnippetExprEvalSum] */
796  	
797  	   return SCIP_OKAY;
798  	}
799  	
800  	/** expression forward derivative evaluation callback */
801  	static
802  	SCIP_DECL_EXPRFWDIFF(fwdiffSum)
803  	{  /*lint --e{715}*/
804  	   SCIP_EXPRDATA* exprdata;
805  	   int c;
806  	
807  	   assert(expr != NULL);
808  	   assert(dot != NULL);
809  	
810  	   exprdata = SCIPexprGetData(expr);
811  	   assert(exprdata != NULL);
812  	
813  	   *dot = 0.0;
814  	   for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
815  	   {
816  	      assert(SCIPexprGetDot(SCIPexprGetChildren(expr)[c]) != SCIP_INVALID); /*lint !e777*/
817  	
818  	      *dot += exprdata->coefficients[c] * SCIPexprGetDot(SCIPexprGetChildren(expr)[c]);
819  	   }
820  	
821  	   return SCIP_OKAY;
822  	}
823  	
824  	/** expression derivative evaluation callback */
825  	static
826  	SCIP_DECL_EXPRBWDIFF(bwdiffSum)
827  	{  /*lint --e{715}*/
828  	   assert(expr != NULL);
829  	   assert(SCIPexprGetData(expr) != NULL);
830  	   assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr));
831  	   assert(SCIPexprGetChildren(expr)[childidx] != NULL);
832  	   assert(!SCIPisExprValue(scip, SCIPexprGetChildren(expr)[childidx]));
833  	
834  	   *val = SCIPgetCoefsExprSum(expr)[childidx];
835  	
836  	   return SCIP_OKAY;
837  	}
838  	
839  	/** expression backward forward derivative evaluation callback */
840  	static
841  	SCIP_DECL_EXPRBWFWDIFF(bwfwdiffSum)
842  	{  /*lint --e{715}*/
843  	   assert(bardot != NULL);
844  	
845  	   *bardot = 0.0;
846  	
847  	   return SCIP_OKAY;
848  	}
849  	
850  	/** expression interval evaluation callback */
851  	static
852  	SCIP_DECL_EXPRINTEVAL(intevalSum)
853  	{  /*lint --e{715}*/
854  	   SCIP_EXPRDATA* exprdata;
855  	   SCIP_INTERVAL suminterval;
856  	   int c;
857  	
858  	   assert(expr != NULL);
859  	
860  	   exprdata = SCIPexprGetData(expr);
861  	   assert(exprdata != NULL);
862  	
863  	   SCIPintervalSet(interval, exprdata->constant);
864  	
865  	   SCIPdebugMsg(scip, "inteval %p with %d children: %.20g", (void*)expr, SCIPexprGetNChildren(expr), exprdata->constant);
866  	
867  	   for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
868  	   {
869  	      SCIP_INTERVAL childinterval;
870  	
871  	      childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[c]);
872  	      if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
873  	      {
874  	         SCIPintervalSetEmpty(interval);
875  	         break;
876  	      }
877  	
878  	      /* compute coefficients[c] * childinterval and add the result to the so far computed interval */
879  	      if( exprdata->coefficients[c] == 1.0 )
880  	      {
881  	         SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, *interval, childinterval);
882  	      }
883  	      else
884  	      {
885  	         SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &suminterval, childinterval, exprdata->coefficients[c]);
886  	         SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, *interval, suminterval);
887  	      }
888  	
889  	      SCIPdebugMsgPrint(scip, " %+.20g*[%.20g,%.20g]", exprdata->coefficients[c], childinterval.inf, childinterval.sup);
890  	   }
891  	   SCIPdebugMsgPrint(scip, " = [%.20g,%.20g]\n", interval->inf, interval->sup);
892  	
893  	   return SCIP_OKAY;
894  	}
895  	
896  	/** initial estimators callback */
897  	static
898  	SCIP_DECL_EXPRINITESTIMATES(initEstimatesSum)
899  	{  /*lint --e{715}*/
900  	   SCIP_EXPRDATA* exprdata;
901  	
902  	#ifdef SCIP_DEBUG
903  	   SCIPinfoMessage(scip, NULL, "initEstimatesSum %d children: ", SCIPexprGetNChildren(expr));
904  	   SCIPprintExpr(scip, expr, NULL);
905  	   SCIPinfoMessage(scip, NULL, "\n");
906  	#endif
907  	   assert(scip != NULL);
908  	   assert(expr != NULL);
909  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
910  	
911  	   assert(coefs[0] != NULL);
912  	   assert(constant != NULL);
913  	   assert(nreturned != NULL);
914  	
915  	   exprdata = SCIPexprGetData(expr);
916  	   assert(exprdata != NULL);
917  	
918  	   BMScopyMemoryArray(coefs[0], exprdata->coefficients, SCIPexprGetNChildren(expr));
919  	   *constant = exprdata->constant;
920  	   *nreturned = 1;
921  	
922  	   return SCIP_OKAY;
923  	}
924  	
925  	/** expression estimate callback */
926  	static
927  	SCIP_DECL_EXPRESTIMATE(estimateSum)
928  	{  /*lint --e{715}*/
929  	   SCIP_EXPRDATA* exprdata;
930  	
931  	   assert(scip != NULL);
932  	   assert(expr != NULL);
933  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
934  	   assert(islocal != NULL);
935  	   assert(success != NULL);
936  	   assert(branchcand != NULL);
937  	
938  	   exprdata = SCIPexprGetData(expr);
939  	   assert(exprdata != NULL);
940  	
941  	   /* NOTE: nlhdlr_default assumes in nlhdlrInitSepaDefault that this estimator can be used for both under- and overestimation */
942  	
943  	   BMScopyMemoryArray(coefs, exprdata->coefficients, SCIPexprGetNChildren(expr));
944  	   *constant = exprdata->constant;
945  	   *islocal = FALSE;
946  	   *success = TRUE;
947  	
948  	   /* for none of our children, branching would improve the underestimator, so set branchcand[i]=FALSE everywhere
949  	    * if we branch for numerical reasons, then cons-expr-core should figure out what the candidates are
950  	    */
951  	   BMSclearMemoryArray(branchcand, SCIPexprGetNChildren(expr));
952  	
953  	   return SCIP_OKAY;
954  	}
955  	
956  	/** expression reverse propagation callback */
957  	static
958  	SCIP_DECL_EXPRREVERSEPROP(reversepropSum)
959  	{  /*lint --e{715}*/
960  	   SCIP_EXPRDATA* exprdata;
961  	   SCIP_INTERVAL* newbounds;
962  	   int nchildren;
963  	   int nreductions;
964  	
965  	   assert(scip != NULL);
966  	   assert(expr != NULL);
967  	   assert(infeasible != NULL);
968  	
969  	   nchildren = SCIPexprGetNChildren(expr);
970  	   assert(nchildren > 0);
971  	
972  	   exprdata = SCIPexprGetData(expr);
973  	   assert(exprdata != NULL);
974  	
975  	   SCIP_CALL( SCIPallocBufferArray(scip, &newbounds, nchildren) );
976  	
977  	   nreductions = SCIPintervalPropagateWeightedSum(SCIP_INTERVAL_INFINITY, nchildren, childrenbounds,
978  	         exprdata->coefficients, exprdata->constant, bounds, newbounds, infeasible);
979  	
980  	   if( !*infeasible && nreductions > 0 )
981  	      BMScopyMemoryArray(childrenbounds, newbounds, nchildren);
982  	
983  	   SCIPfreeBufferArray(scip, &newbounds);
984  	
985  	   return SCIP_OKAY;
986  	}
987  	
988  	/** sum hash callback */
989  	static
990  	SCIP_DECL_EXPRHASH(hashSum)
991  	{  /*lint --e{715}*/
992  	   SCIP_EXPRDATA* exprdata;
993  	   int c;
994  	
995  	   assert(scip != NULL);
996  	   assert(expr != NULL);
997  	   assert(hashkey != NULL);
998  	   assert(childrenhashes != NULL);
999  	
1000 	   exprdata = SCIPexprGetData(expr);
1001 	   assert(exprdata != NULL);
1002 	
1003 	   /**! [SnippetExprHashSum] */
1004 	   *hashkey = EXPRHDLR_HASHKEY;
1005 	   *hashkey ^= SCIPcalcFibHash(exprdata->constant);
1006 	
1007 	   for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
1008 	      *hashkey ^= SCIPcalcFibHash(exprdata->coefficients[c]) ^ childrenhashes[c];
1009 	   /**! [SnippetExprHashSum] */
1010 	
1011 	   return SCIP_OKAY;
1012 	}
1013 	
1014 	/** expression curvature detection callback */
1015 	static
1016 	SCIP_DECL_EXPRCURVATURE(curvatureSum)
1017 	{  /*lint --e{715}*/
1018 	   SCIP_EXPRDATA* exprdata;
1019 	   int i;
1020 	
1021 	   assert(scip != NULL);
1022 	   assert(expr != NULL);
1023 	   assert(childcurv != NULL);
1024 	   assert(success != NULL);
1025 	
1026 	   exprdata = SCIPexprGetData(expr);
1027 	   assert(exprdata != NULL);
1028 	
1029 	   for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
1030 	      childcurv[i] = SCIPexprcurvMultiply(exprdata->coefficients[i], exprcurvature);
1031 	
1032 	   *success = TRUE;
1033 	
1034 	   return SCIP_OKAY;
1035 	}
1036 	
1037 	/** expression monotonicity detection callback */
1038 	static
1039 	SCIP_DECL_EXPRMONOTONICITY(monotonicitySum)
1040 	{  /*lint --e{715}*/
1041 	   SCIP_EXPRDATA* exprdata;
1042 	
1043 	   assert(scip != NULL);
1044 	   assert(expr != NULL);
1045 	   assert(result != NULL);
1046 	   assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr));
1047 	
1048 	   exprdata = SCIPexprGetData(expr);
1049 	   assert(exprdata != NULL);
1050 	
1051 	   *result = exprdata->coefficients[childidx] >= 0.0 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC;
1052 	
1053 	   return SCIP_OKAY;
1054 	}
1055 	
1056 	/** expression integrality detection callback */
1057 	static
1058 	SCIP_DECL_EXPRINTEGRALITY(integralitySum)
1059 	{  /*lint --e{715}*/
1060 	   SCIP_EXPRDATA* exprdata;
1061 	   int i;
1062 	
1063 	   assert(scip != NULL);
1064 	   assert(expr != NULL);
1065 	   assert(isintegral != NULL);
1066 	
1067 	   exprdata = SCIPexprGetData(expr);
1068 	   assert(exprdata != NULL);
1069 	
1070 	   /**! [SnippetExprIntegralitySum] */
1071 	   *isintegral = EPSISINT(exprdata->constant, 0.0); /*lint !e835*/
1072 	
1073 	   for( i = 0; i < SCIPexprGetNChildren(expr) && *isintegral; ++i )
1074 	   {
1075 	      SCIP_EXPR* child = SCIPexprGetChildren(expr)[i];
1076 	      assert(child != NULL);
1077 	
1078 	      *isintegral = EPSISINT(exprdata->coefficients[i], 0.0) && SCIPexprIsIntegral(child); /*lint !e835*/
1079 	   }
1080 	   /**! [SnippetExprIntegralitySum] */
1081 	
1082 	   return SCIP_OKAY;
1083 	}
1084 	
1085 	/** creates the handler for sum expressions and includes it into SCIP */
1086 	SCIP_RETCODE SCIPincludeExprhdlrSum(
1087 	   SCIP*                 scip                /**< SCIP data structure */
1088 	   )
1089 	{
1090 	   SCIP_EXPRHDLR* exprhdlr;
1091 	
1092 	   SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, EXPRHDLR_PRECEDENCE, evalSum, NULL) );
1093 	   assert(exprhdlr != NULL);
1094 	
1095 	   SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrSum, NULL);
1096 	   SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataSum, freedataSum);
1097 	   SCIPexprhdlrSetSimplify(exprhdlr, simplifySum);
1098 	   SCIPexprhdlrSetCompare(exprhdlr, compareSum);
1099 	   SCIPexprhdlrSetPrint(exprhdlr, printSum);
1100 	   SCIPexprhdlrSetIntEval(exprhdlr, intevalSum);
1101 	   SCIPexprhdlrSetEstimate(exprhdlr, initEstimatesSum, estimateSum);
1102 	   SCIPexprhdlrSetReverseProp(exprhdlr, reversepropSum);
1103 	   SCIPexprhdlrSetHash(exprhdlr, hashSum);
1104 	   SCIPexprhdlrSetDiff(exprhdlr, bwdiffSum, fwdiffSum, bwfwdiffSum);
1105 	   SCIPexprhdlrSetCurvature(exprhdlr, curvatureSum);
1106 	   SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicitySum);
1107 	   SCIPexprhdlrSetIntegrality(exprhdlr, integralitySum);
1108 	   SCIPexprhdlrSetGetSymdata(exprhdlr, getSymDataSum);
1109 	
1110 	   return SCIP_OKAY;
1111 	}
1112 	
1113 	/** creates a sum expression */
1114 	SCIP_RETCODE SCIPcreateExprSum(
1115 	   SCIP*                 scip,               /**< SCIP data structure */
1116 	   SCIP_EXPR**           expr,               /**< pointer where to store expression */
1117 	   int                   nchildren,          /**< number of children */
1118 	   SCIP_EXPR**           children,           /**< children */
1119 	   SCIP_Real*            coefficients,       /**< array with coefficients for all children (or NULL if all 1.0) */
1120 	   SCIP_Real             constant,           /**< constant term of sum */
1121 	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
1122 	   void*                 ownercreatedata     /**< data to pass to ownercreate */
1123 	   )
1124 	{
1125 	   SCIP_EXPRDATA* exprdata;
1126 	
1127 	   SCIP_CALL( createData(scip, &exprdata, nchildren, coefficients, constant) );
1128 	
1129 	   SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPgetExprhdlrSum(scip), exprdata, nchildren, children, ownercreate, ownercreatedata) );
1130 	
1131 	   return SCIP_OKAY;
1132 	}
1133 	
1134 	/** sets the constant of a summation expression */
1135 	void SCIPsetConstantExprSum(
1136 	   SCIP_EXPR*            expr,               /**< sum expression */
1137 	   SCIP_Real             constant            /**< constant */
1138 	   )
1139 	{
1140 	   SCIP_EXPRDATA* exprdata;
1141 	
1142 	   assert(expr != NULL);
1143 	
1144 	   exprdata = SCIPexprGetData(expr);
1145 	   assert(exprdata != NULL);
1146 	
1147 	   exprdata->constant = constant;
1148 	}
1149 	
1150 	/** appends an expression to a sum expression */
1151 	SCIP_RETCODE SCIPappendExprSumExpr(
1152 	   SCIP*                 scip,               /**< SCIP data structure */
1153 	   SCIP_EXPR*            expr,               /**< sum expression */
1154 	   SCIP_EXPR*            child,              /**< expression to be appended */
1155 	   SCIP_Real             childcoef           /**< child's coefficient */
1156 	   )
1157 	{
1158 	   SCIP_EXPRDATA* exprdata;
1159 	   int nchildren;
1160 	
1161 	   assert(expr != NULL);
1162 	   assert(SCIPisExprSum(scip, expr));
1163 	
1164 	   exprdata = SCIPexprGetData(expr);
1165 	   assert(exprdata != NULL);
1166 	
1167 	   nchildren = SCIPexprGetNChildren(expr);
1168 	
1169 	   SCIP_CALL( SCIPensureBlockMemoryArray(scip, &exprdata->coefficients, &exprdata->coefssize, nchildren + 1) );
1170 	
1171 	   assert(exprdata->coefssize > nchildren);
1172 	   exprdata->coefficients[nchildren] = childcoef;
1173 	
1174 	   SCIP_CALL( SCIPappendExprChild(scip, expr, child) );
1175 	
1176 	   return SCIP_OKAY;
1177 	}
1178 	
1179 	/** multiplies given sum expression by a constant */
1180 	void SCIPmultiplyByConstantExprSum(
1181 	   SCIP_EXPR*            expr,               /**< sum expression */
1182 	   SCIP_Real             constant            /**< constant that multiplies sum expression */
1183 	   )
1184 	{
1185 	   int i;
1186 	   SCIP_EXPRDATA* exprdata;
1187 	
1188 	   assert(expr != NULL);
1189 	
1190 	   exprdata = SCIPexprGetData(expr);
1191 	   assert(exprdata != NULL);
1192 	
1193 	   for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
1194 	      exprdata->coefficients[i] *= constant;
1195 	   exprdata->constant *= constant;
1196 	}
1197 	
1198 	/** constructs the expanded product of two sum expressions */
1199 	SCIP_RETCODE SCIPmultiplyBySumExprSum(
1200 	   SCIP*                 scip,               /**< SCIP data structure */
1201 	   SCIP_EXPR**           product,            /**< buffer where to store multiplied sums (expanded as sum) */
1202 	   SCIP_EXPR*            factor1,            /**< first sum */
1203 	   SCIP_EXPR*            factor2,            /**< second sum */
1204 	   SCIP_Bool             simplify,           /**< whether to simplify created terms and sum */
1205 	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
1206 	   void*                 ownercreatedata     /**< data to pass to ownercreate */
1207 	   )
1208 	{
1209 	   SCIP_Real constant1;
1210 	   SCIP_Real constant2;
1211 	   int nchildren1;
1212 	   int nchildren2;
1213 	   int i1;
1214 	   int i2;
1215 	
1216 	   assert(scip != NULL);
1217 	   assert(product != NULL);
1218 	   assert(factor1 != NULL);
1219 	   assert(SCIPisExprSum(scip, factor1));
1220 	   assert(factor2 != NULL);
1221 	   assert(SCIPisExprSum(scip, factor2));
1222 	
1223 	   constant1 = SCIPgetConstantExprSum(factor1);
1224 	   constant2 = SCIPgetConstantExprSum(factor2);
1225 	   nchildren1 = SCIPexprGetNChildren(factor1);
1226 	   nchildren2 = SCIPexprGetNChildren(factor2);
1227 	
1228 	   /* TODO might be nice to integrate more with simplify and construct a simplified sum right away */
1229 	
1230 	   SCIP_CALL( SCIPcreateExprSum(scip, product, 0, NULL, NULL, constant1 * constant2, ownercreate, ownercreatedata) );
1231 	
1232 	   /* first add constant1 * factor2
1233 	    * constant * constant2 already added above
1234 	    */
1235 	   if( constant1 != 0.0 )
1236 	   {
1237 	      for( i2 = 0; i2 < nchildren2; ++i2 )
1238 	      {
1239 	         SCIP_CALL( SCIPappendExprSumExpr(scip, *product, SCIPexprGetChildren(factor2)[i2], constant1 * SCIPgetCoefsExprSum(factor2)[i2]) );
1240 	      }
1241 	   }
1242 	
1243 	   for( i1 = 0; i1 < nchildren1; ++i1 )
1244 	   {
1245 	      SCIP_EXPR* child1;
1246 	      SCIP_EXPR* child2;
1247 	      SCIP_Real coef1;
1248 	      SCIP_Real coef2;
1249 	
1250 	      coef1 = SCIPgetCoefsExprSum(factor1)[i1];
1251 	      child1 = SCIPexprGetChildren(factor1)[i1];
1252 	
1253 	      if( constant2 != 0.0 )
1254 	      {
1255 	         /* add coef1 * child1 * constant2 */
1256 	         SCIP_CALL( SCIPappendExprSumExpr(scip, *product, child1, coef1 * constant2) );
1257 	      }
1258 	
1259 	      for( i2 = 0; i2 < nchildren2; ++i2 )
1260 	      {
1261 	         /* add coef1 * child1 * coef2 * child2 */
1262 	         SCIP_EXPR* termprod;
1263 	         SCIP_EXPR* termprodsimplified;
1264 	         SCIP_EXPR* termfactors[2];
1265 	
1266 	         coef2 = SCIPgetCoefsExprSum(factor2)[i2];
1267 	         child2 = SCIPexprGetChildren(factor2)[i2];
1268 	
1269 	         /* create child1 * child2 expr */
1270 	         termfactors[0] = child1;
1271 	         termfactors[1] = child2;
1272 	         SCIP_CALL( SCIPcreateExprProduct(scip, &termprod, 2, termfactors, 1.0, NULL, NULL) );
1273 	
1274 	         if( simplify )
1275 	         {
1276 	            SCIP_Bool changed;
1277 	            SCIP_Bool infeasible;
1278 	
1279 	            /* simplify child1*child2 */
1280 	            SCIP_CALL( SCIPsimplifyExpr(scip, termprod, &termprodsimplified, &changed, &infeasible, ownercreate, ownercreatedata) );
1281 	            assert(!infeasible);  /* simplify of products should never be infeasible */
1282 	            SCIP_CALL( SCIPreleaseExpr(scip, &termprod) );
1283 	         }
1284 	         else
1285 	            termprodsimplified = termprod;
1286 	
1287 	         /* append to product */
1288 	         SCIP_CALL( SCIPappendExprSumExpr(scip, *product, termprodsimplified, coef1 * coef2) );
1289 	         SCIP_CALL( SCIPreleaseExpr(scip, &termprodsimplified) );
1290 	      }
1291 	   }
1292 	
1293 	   if( simplify )
1294 	   {
1295 	      SCIP_EXPR* prodsimplified;
1296 	      SCIP_Bool changed;
1297 	      SCIP_Bool infeasible;
1298 	
1299 	      SCIP_CALL( SCIPsimplifyExpr(scip, *product, &prodsimplified, &changed, &infeasible, ownercreate, ownercreatedata) );
1300 	      assert(!infeasible);
1301 	      SCIP_CALL( SCIPreleaseExpr(scip, product) );
1302 	      *product = prodsimplified;
1303 	   }
1304 	
1305 	   return SCIP_OKAY;
1306 	}
1307 	
1308 	/** constructs the expanded power of a sum expression
1309 	 *
1310 	 * @attention The number of terms in the expansion grows exponential with the exponent. Be aware of what you wish for.
1311 	 */
1312 	SCIP_EXPORT
1313 	SCIP_RETCODE SCIPpowerExprSum(
1314 	   SCIP*                 scip,               /**< SCIP data structure */
1315 	   SCIP_EXPR**           result,             /**< buffer where to store expanded power of sum */
1316 	   SCIP_EXPR*            base,               /**< sum */
1317 	   int                   exponent,           /**< exponent > 1 */
1318 	   SCIP_Bool             simplify,           /**< whether to simplify created terms and sum */
1319 	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
1320 	   void*                 ownercreatedata     /**< data to pass to ownercreate */
1321 	   )
1322 	{
1323 	   /* for sum_i alpha_i expr_i  and exponent p, this constructs
1324 	   *   sum_{beta} (p over beta) prod_i (alpha_i expr_i)^beta_i   over all multiindex beta such that sum_i beta_i = p
1325 	   * See also https://en.wikipedia.org/wiki/Multinomial_theorem
1326 	   * Calculation of multinomials adapted from https://github.com/m-j-w/MultinomialSeries.jl
1327 	   *
1328 	   * TODO might be nice to integrate more with simplify and construct a simplified sum right away
1329 	   */
1330 	
1331 	   SCIP_EXPR** children;
1332 	   SCIP_EXPR*** childrenpower;
1333 	   SCIP_Real* coefs;
1334 	   SCIP_Real constant;
1335 	   SCIP_Bool haveconst;
1336 	   int nchildren;
1337 	   int nterms;
1338 	   int* factorials;
1339 	   int* beta;
1340 	   int betapos;
1341 	   int restsum;
1342 	   int multinomialcoef;
1343 	   int i;
1344 	
1345 	   SCIP_EXPR* newterm;
1346 	   SCIP_Real newtermcoef;
1347 	
1348 	   assert(scip != NULL);
1349 	   assert(result != NULL);
1350 	   assert(base != NULL);
1351 	   assert(exponent > 1);
1352 	   assert(SCIPisExprSum(scip, base));
1353 	   assert(exponent > 1.0);
1354 	
1355 	   children = SCIPexprGetChildren(base);
1356 	   nchildren = SCIPexprGetNChildren(base);
1357 	   coefs = SCIPgetCoefsExprSum(base);
1358 	   constant = SCIPgetConstantExprSum(base);
1359 	   haveconst = constant != 0.0;
1360 	   nterms = nchildren + (haveconst ? 1 : 0);
1361 	
1362 	#ifdef SCIP_DEBUG
1363 	   SCIPinfoMessage(scip, NULL, "expanding (");
1364 	   SCIPprintExpr(scip, base, NULL);
1365 	   SCIPinfoMessage(scip, NULL, ")^%d\n", exponent);
1366 	#endif
1367 	
1368 	   SCIP_CALL( SCIPcreateExprSum(scip, result, 0, NULL, NULL, 0.0, ownercreate, ownercreatedata) );
1369 	
1370 	   SCIP_CALL( SCIPallocClearBufferArray(scip, &beta, nterms) );
1371 	   SCIP_CALL( SCIPallocBufferArray(scip, &factorials, exponent+1) );
1372 	
1373 	   /* precompute factorials k!, k = 0...exponent */
1374 	   factorials[0] = 1;
1375 	   for( i = 1; i <= exponent; ++i )
1376 	     factorials[i] = factorials[i-1] * i;
1377 	
1378 	   /* precompute children^k, k=1..exponent */
1379 	   SCIP_CALL( SCIPallocBufferArray(scip, &childrenpower, nchildren) );
1380 	   for( i = 0; i < nchildren; ++i )
1381 	   {
1382 	      int k;
1383 	      SCIP_CALL( SCIPallocBufferArray(scip, &childrenpower[i], exponent+1) );
1384 	      childrenpower[i][1] = children[i];
1385 	      for( k = 2; k <= exponent; ++k )
1386 	      {
1387 	         SCIP_CALL( SCIPcreateExprPow(scip, &childrenpower[i][k], children[i], (SCIP_Real)k, NULL, NULL) );
1388 	         if( simplify )
1389 	         {
1390 	            SCIP_Bool changed;
1391 	            SCIP_Bool infeasible;
1392 	            SCIP_EXPR* simplified;
1393 	
1394 	            SCIP_CALL( SCIPsimplifyExpr(scip, childrenpower[i][k], &simplified, &changed, &infeasible, ownercreate, ownercreatedata) );
1395 	            assert(!infeasible);
1396 	            SCIP_CALL( SCIPreleaseExpr(scip, &childrenpower[i][k]) );
1397 	            childrenpower[i][k] = simplified;
1398 	         }
1399 	      }
1400 	   }
1401 	
1402 	   /* first multinomial is (exponent,0,0,...) */
1403 	   beta[0] = exponent;
1404 	   betapos = 0;
1405 	   do
1406 	   {
1407 	      /* compute multinomial coef exponent! /  (beta[0]! * ... * beta[nterms-1]!) */
1408 	      multinomialcoef = factorials[exponent];
1409 	      for( i = 0; i < nterms; ++i )
1410 	      {
1411 	         assert(beta[i] >= 0);
1412 	         assert(beta[i] <= exponent);
1413 	         multinomialcoef /= factorials[beta[i]];
1414 	      }
1415 	
1416 	      SCIPdebugMsg(scip, "multinomial (");
1417 	      for( i = 0; i < nterms; ++i )
1418 	         SCIPdebugPrintf("%d ", beta[i]);
1419 	      SCIPdebugPrintf(")  with coef %d\n", multinomialcoef);
1420 	
1421 	      /* construct new term for expanded sum */
1422 	      SCIP_CALL( SCIPcreateExprProduct(scip, &newterm, 0, NULL, 1.0, ownercreate, ownercreatedata) );
1423 	      newtermcoef = multinomialcoef;
1424 	      for( i = 0; i < nterms; ++i )
1425 	      {
1426 	         if( beta[i] == 0 )
1427 	            continue;
1428 	
1429 	         if( i == nterms-1 && haveconst )
1430 	         {
1431 	            /* if constant term, then update newtermcoef */
1432 	            newtermcoef *= pow(constant, (double)beta[i]);
1433 	            continue;
1434 	         }
1435 	
1436 	         /* alpha_i^beta_i */
1437 	         newtermcoef *= pow(coefs[i], (double)beta[i]);
1438 	
1439 	         /* expr_i^beta_i*/
1440 	         SCIP_CALL( SCIPappendExprChild(scip, newterm, childrenpower[i][beta[i]]) );
1441 	      }
1442 	
1443 	      /* append newterm to sum */
1444 	      switch( SCIPexprGetNChildren(newterm) )
1445 	      {
1446 	         case 0:
1447 	         {
1448 	            /* no factor in product, so it is a constant -> update constant in sum */
1449 	            SCIPsetConstantExprSum(*result, SCIPgetConstantExprSum(*result) + newtermcoef);
1450 	            break;
1451 	         }
1452 	
1453 	         case 1:
1454 	         {
1455 	            /* only one factor in product -> add this factor itself to sum */
1456 	            SCIP_CALL( SCIPappendExprSumExpr(scip, *result, SCIPexprGetChildren(newterm)[0], newtermcoef) );
1457 	            break;
1458 	         }
1459 	
1460 	         default:
1461 	         {
1462 	            if( simplify )
1463 	            {
1464 	               /* simplify product */
1465 	               SCIP_Bool changed;
1466 	               SCIP_Bool infeasible;
1467 	               SCIP_EXPR* simplified;
1468 	
1469 	               SCIP_CALL( SCIPsimplifyExpr(scip, newterm, &simplified, &changed, &infeasible, ownercreate, ownercreatedata) );
1470 	               assert(!infeasible);
1471 	               SCIP_CALL( SCIPreleaseExpr(scip, &newterm) );
1472 	               newterm = simplified;
1473 	            }
1474 	
1475 	            /* append new term to sum */
1476 	            SCIP_CALL( SCIPappendExprSumExpr(scip, *result, newterm, newtermcoef) );
1477 	            break;
1478 	         }
1479 	      }
1480 	      SCIP_CALL( SCIPreleaseExpr(scip, &newterm) );
1481 	
1482 	      /* determine next beta */
1483 	      while( beta[betapos] == 0 )
1484 	         --betapos;
1485 	
1486 	      if( betapos == nterms-1 )
1487 	      {
1488 	         do
1489 	         {
1490 	            if( betapos == 0 )
1491 	               goto TERMINATE;
1492 	            --betapos;
1493 	         }
1494 	         while( beta[betapos] == 0 );
1495 	
1496 	         restsum = 0;
1497 	         for( i = betapos+1; i < nterms; ++i )
1498 	         {
1499 	            restsum += beta[i];
1500 	            beta[i] = 0;
1501 	         }
1502 	         beta[betapos+1] = restsum;
1503 	      }
1504 	
1505 	      if( beta[betapos] > 0 )
1506 	      {
1507 	         --beta[betapos];
1508 	         ++beta[++betapos];
1509 	      }
1510 	   }
1511 	   while( TRUE );  /*lint !e506*/
1512 	
1513 	 TERMINATE:
1514 	
1515 	   if( simplify )
1516 	   {
1517 	      SCIP_Bool changed;
1518 	      SCIP_Bool infeasible;
1519 	      SCIP_EXPR* simplified;
1520 	
1521 	      SCIP_CALL( SCIPsimplifyExpr(scip, *result, &simplified, &changed, &infeasible, ownercreate, ownercreatedata) );
1522 	      assert(!infeasible);
1523 	      SCIP_CALL( SCIPreleaseExpr(scip, result) );
1524 	      *result = simplified;
1525 	   }
1526 	
1527 	   for( i = nchildren-1; i >= 0; --i )
1528 	   {
1529 	      int k;
1530 	      for( k = exponent; k >= 2; --k )
1531 	      {
1532 	         SCIP_CALL( SCIPreleaseExpr(scip, &childrenpower[i][k]) );
1533 	      }
1534 	      SCIPfreeBufferArray(scip, &childrenpower[i]);
1535 	   }
1536 	   SCIPfreeBufferArray(scip, &childrenpower);
1537 	   SCIPfreeBufferArray(scip, &factorials);
1538 	   SCIPfreeBufferArray(scip, &beta);
1539 	
1540 	#ifdef SCIP_DEBUG
1541 	   SCIPinfoMessage(scip, NULL, "-> ");
1542 	   SCIPprintExpr(scip, *result, NULL);
1543 	   SCIPinfoMessage(scip, NULL, "\n");
1544 	#endif
1545 	
1546 	   return SCIP_OKAY;
1547 	}
1548 	
1549 	/* from pub_expr.h */
1550 	
1551 	/** gets the coefficients of a summation expression */
1552 	SCIP_Real* SCIPgetCoefsExprSum(
1553 	   SCIP_EXPR*            expr                /**< sum expression */
1554 	   )
1555 	{
1556 	   SCIP_EXPRDATA* exprdata;
1557 	
1558 	   assert(expr != NULL);
1559 	
1560 	   exprdata = SCIPexprGetData(expr);
1561 	   assert(exprdata != NULL);
1562 	
1563 	   return exprdata->coefficients;
1564 	}
1565 	
1566 	/** gets the constant of a summation expression */
1567 	SCIP_Real SCIPgetConstantExprSum(
1568 	   SCIP_EXPR*            expr                /**< sum expression */
1569 	   )
1570 	{
1571 	   SCIP_EXPRDATA* exprdata;
1572 	
1573 	   assert(expr != NULL);
1574 	
1575 	   exprdata = SCIPexprGetData(expr);
1576 	   assert(exprdata != NULL);
1577 	
1578 	   return exprdata->constant;
1579 	}
1580