1 /* $Id$
2 *
3 * Name: conv-exprMul-reformulate.cpp
4 * Author: Pietro Belotti
5 * Purpose: utility methods to reformulate multiplications into bilinear and trilinear terms
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include <queue>
12
13 #include "CouenneTypes.hpp"
14 #include "CouenneExprMul.hpp"
15 #include "CouenneExprTrilinear.hpp"
16 #include "CouenneExprBMul.hpp"
17 #include "CouenneExprConst.hpp"
18 #include "CouenneExprPow.hpp"
19 #include "CouenneExprAux.hpp"
20 #include "CouenneExprClone.hpp"
21 #include "CouenneProblem.hpp"
22
23 using namespace Couenne;
24
25 /// check if two arguments point to the same variable
26
areSameVariables(expression * v1,expression * v2)27 inline bool areSameVariables (expression *v1, expression *v2) {
28
29 register int t1 = v1 -> Type (), t2;
30 return (((t1 == VAR) || (t1 == AUX)) &&
31 (((t2 = v2 -> Type ()) == VAR) || (t2 == AUX)) &&
32 (v1 -> Index () == v2 -> Index ()));
33 }
34
35
36 /// Create standard formulation of this expression
37
standardize(CouenneProblem * p,bool addAux)38 exprAux *exprMul::standardize (CouenneProblem *p, bool addAux) {
39
40 //printf ("standardizing exprMul (addaux=%d) ", addAux); fflush (stdout); print (); printf ("\n");
41
42 CouNumber coeff = 1.;
43 std::map <int, CouNumber> indCoe;
44
45 p -> flattenMul (this, coeff, indCoe);
46
47 int nArgs = 0;
48
49 expression **arglist = new expression * [indCoe.size ()];
50
51 for (std::map <int, CouNumber>::iterator i = indCoe.begin (); i != indCoe.end (); ++i)
52 if (i -> second == 1.) arglist [nArgs++] = new exprClone (p -> Var (i -> first));
53 else arglist [nArgs++] = new exprPow (new exprClone (p -> Var (i -> first)), new exprConst (i -> second));
54
55 //while (nargs_--)
56 //delete arglist_ [nargs_];
57
58 //delete [] arglist_;
59
60 arglist_ = arglist;
61 nargs_ = (int) indCoe.size();
62
63 //printf ("new mul [%d]: %g * ", nargs_, coeff); fflush (stdout); print (); printf (" -- ");
64
65 exprOp::standardize (p);
66
67 //printf ("standardized: "); fflush (stdout); print (); printf ("\n");
68
69 if (nargs_ == 1) {
70
71 if (coeff != 1.) {
72
73 // leak: arglist_ should be deleted
74
75 expression *simMul = new exprMul (new exprConst (coeff), new exprClone (arglist_ [0]));
76 return (addAux ? (p -> addAuxiliary (simMul)) : new exprAux (simMul, p -> domain ()));
77
78 } else {
79
80 // quick fix for nvs08: expression x_0^.5 * x_0^3 is quite
81 // unusual and is translated into an exprMul with a single
82 // argument.
83
84 exprAux *var = dynamic_cast <exprAux *> (arglist_ [0] -> Copy ());
85 return var;
86 }
87 }
88
89 //if (nargs_ == 1) // TODO: what happens really?
90 //return NULL;
91 /* {
92 exprAux *aux = arglist_ [0];
93 arglist_ [0] = NULL;
94 return aux;
95 } */
96
97 // enable this when binary products are available
98
99 #if 0
100 // check if it is a product of binary variables
101 bool isBinProd = true;
102 for (int i=nargs_; i--;) {
103
104 expression *arg = arglist_ [i];
105 if (arg -> isInteger ()) {
106
107 CouNumber lb, ub;
108 arg -> getBounds (lb, ub);
109 if ((fabs (ceil (lb - COUENNE_EPS)) > 0.) ||
110 (fabs (floor (ub + COUENNE_EPS) - 1.) > 0.)) { // if not all conditions hold,
111 isBinProd = false;
112 break;
113 }
114 } else {
115 isBinProd = false;
116 break;
117 }
118 }
119
120 if (isBinProd) {
121 //printf ("found a BinProd!\n");
122 }
123 #endif
124
125 enum Couenne::TrilinDecompType type = p -> getTrilinDecompType ();
126
127 if (nargs_ <= 2)
128 type = Couenne::rAI;
129
130 expression *retExpr;
131
132 switch (type) {
133
134 case Couenne::treeDecomp: {
135
136 //printf ("trying treeDecomp on "); print (); printf ("\n");
137
138 // A more hierarchical decomposition. Example:
139 //
140 // (x1 x2 x3 x4 x5 x6 x7 x8 x9 x10)
141 //
142 // becomes
143 //
144 // ((((x1 x2) (x3 x4)) ((x5 x6) (x7 x8))) (x9 x10))
145 //
146 // so that x1 .. x10 are the leaves of a binary tree whose
147 // non-leaf nodes are auxiliaries
148 //
149 // use a queue to parse the leaf level first (by pushing all
150 // original members) and push the newly added auxs into the queue
151 // while extracting _pairs_ of elements
152
153 std::queue <expression *> queue;
154
155 for (int i=0; i<nargs_; i++)
156 queue.push (arglist_ [i]);
157
158 expression *aux;
159
160 while (queue.size() > 1) {
161
162 expression *arg1 = queue.front (); queue.pop ();
163 expression *arg2 = queue.front (); queue.pop ();
164
165 //printf ("Coupling "); arg1 -> print ();
166 //printf (" with "); arg2 -> print ();
167
168 if (areSameVariables (arg1, arg2)) aux = new exprPow (new exprClone (arg1), new exprConst (2.));
169 else aux = new exprMul (new exprClone (arg1), new exprClone (arg2));
170
171 //printf (" --> "); aux -> print (); printf ("\n");
172
173 if (!(queue.empty ()))
174 aux = p -> addAuxiliary (aux);
175
176 queue.push (aux);
177 }
178
179 aux = queue.front (); queue.pop ();
180
181 retExpr = aux; //(addAux ? (p -> addAuxiliary (aux)) : new exprAux (this, p -> domain ()));
182 } break;
183
184 // ----------------------------------------------------------------------------------------------
185
186 case Couenne::tri_bi:
187 case Couenne::bi_tri: { // the two cases are very similar
188
189 //printf ("trying %s on ", type==tri_bi ? "tri+bi" : "bi+tri"); print (); printf (": "); fflush (stdout);
190
191 // rAI-tre (ok, funny name, but you get the meaning): same as rAI,
192 // but with trilinear terms. A product
193 //
194 // x1 x2 x3... xk
195 //
196 // is decomposed as
197 //
198 // (...((x1 x2 x3) x4 x5) x6 x7) ... ) x[k-1] xk)
199
200 expression *aux = arglist_ [0];
201
202 for (int i = 1; i < nargs_;) {
203
204 if (i < nargs_ - 1 && ((type==tri_bi) || (i!=1))) { // this is the only point of departure between tri_bi and bi_tri
205
206 // there are at least two remaining arguments: can create trilinear
207 if (areSameVariables (aux, arglist_ [i])) {
208
209 if (areSameVariables (aux, arglist_ [i+1])) // this is a term (x_i x_i x_i)
210
211 if (i == nargs_ - 2) aux = new exprPow (new exprClone (aux), new exprConst (3.));
212 else aux = p -> addAuxiliary (new exprPow (new exprClone (aux), new exprConst (3.)));
213
214 else {
215
216 aux = p -> addAuxiliary (new exprPow (new exprClone (aux), new exprConst (2.)));
217 //aux -> print (); printf (" (tri0) := "); fflush (stdout); aux -> Image () -> print (); printf ("; "); fflush (stdout);
218 if (i == nargs_ - 2) aux = new exprMul (new exprClone (aux), new exprClone (arglist_ [i+1]));
219 else aux = p -> addAuxiliary (new exprMul (new exprClone (aux), new exprClone (arglist_ [i+1])));
220 }
221
222 } else
223 if (areSameVariables (aux, arglist_ [i+1])) {
224
225 printf ("Couenne, exprTrilinear: bad ordering of factors in product, aborting\n");
226 exit (-1);
227
228 } else {
229
230 aux = new exprTrilinear (new exprClone (aux),
231 new exprClone (arglist_ [i]),
232 new exprClone (arglist_ [i+1]));
233
234 if (i != nargs_ - 2)
235 aux = p -> addAuxiliary (aux);
236 }
237
238 //aux -> print (); if (i != nargs_ - 2) {printf (" (tri) := "); aux -> Image () -> print (); printf ("; "); fflush (stdout);}
239
240 i += 2; // covered two variables
241
242 } else {
243
244 if (areSameVariables (aux, arglist_ [i])) aux = new exprPow (new exprClone (aux), new exprConst (2.));
245 else aux = new exprMul (new exprClone (aux), new exprClone (arglist_ [i]));
246
247 if (i==1) // must be bi+tri
248 aux = p -> addAuxiliary (aux); // necessary to introduce the auxiliary variable
249
250 //aux -> print (); if (i==1) {printf (" (bi) := "); fflush (stdout); aux -> Image () -> print (); printf ("; "); fflush (stdout);}
251
252 i++; // covered the last variable
253 }
254 }
255
256 //printf ("\n");
257
258 retExpr = aux; //(addAux ? (p -> addAuxiliary (aux)) : new exprAux (this, p -> domain ()));
259 } break;
260
261 // ----------------------------------------------------------------------------------------------
262
263 case Couenne::rAI:
264 default:
265
266 //printf ("trying good ol' rAI on "); print (); printf ("\n");
267
268 // rAI -- recursive Arithmetic Intervals (see Ryoo and Sahinidis,
269 // JOGO 19 (2001):403-424):
270 //
271 // All multilinear expressions are decomposed as
272 //
273 // (x_1 (x_2 (x_3... (x_{k-1} x_k))...)
274
275 expression *aux = arglist_ [0];
276
277 for (int i = 1; i < nargs_ - 1; i++)
278 aux = (areSameVariables (aux, arglist_ [i])) ?
279 (p -> addAuxiliary (new exprPow (new exprClone (aux), new exprConst (2.)))) :
280 (p -> addAuxiliary (new exprMul (new exprClone (aux), new exprClone (arglist_ [i]))));
281
282 if (areSameVariables (aux, arglist_ [nargs_ - 1]))
283 aux = new exprPow (new exprClone (aux), new exprConst (2.));
284 else aux = new exprMul (new exprClone (aux), new exprClone (arglist_ [nargs_ - 1]));
285
286 retExpr = aux;
287 }
288
289 if (coeff != 1.) {
290
291 retExpr = p -> addAuxiliary (retExpr);
292 retExpr = new exprMul (new exprConst (coeff), new exprClone (retExpr));
293 }
294
295 return (addAux ? (p -> addAuxiliary (retExpr)) : new exprAux (retExpr, p -> domain ()));
296 }
297