1 /* $Id$
2 *
3 * Name: decomposeTerm.cpp
4 * Author: Pietro Belotti
5 * Purpose: decompose sums and products
6 *
7 * (C) Carnegie-Mellon University, 2007-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include <stdio.h>
12
13 #include "CouenneExprQuad.hpp"
14 #include "CouenneJournalist.hpp"
15
16 #include "CouenneProblem.hpp"
17 #include "CouenneExprAux.hpp"
18 #include "CouenneExprMul.hpp"
19 #include "CouenneExprPow.hpp"
20 #include "CouenneLQelems.hpp"
21
22 using namespace Couenne;
23
24 /// given (expression *) element of sum, returns (coe,ind0,ind1)
25 /// depending on element:
26 ///
27 /// 1) a * x_i ^ 2 ---> (a,i,?) return COU_EXPRPOW
28 /// 2) a * x_i ---> (a,i,?) return COU_EXPRVAR
29 /// 3) a * x_i * x_j ---> (a,i,j) return COU_EXPRMUL
30 /// 4) a ---> (a,?,?) return COU_EXPRCONST
31 ///
32 /// x_i and/or x_j may come from standardizing other (linear or
33 /// quadratic operator) sub-expressions
34
decomposeTerm(expression * term,CouNumber initCoe,CouNumber & c0,LinMap & lmap,QuadMap & qmap)35 void CouenneProblem::decomposeTerm (expression *term,
36 CouNumber initCoe,
37 CouNumber &c0,
38 LinMap &lmap,
39 QuadMap &qmap) {
40
41 switch (term -> code ()) {
42
43 // easy cases ////////////////////////////////////////////////////////////////////////
44
45 case COU_EXPRCONST: /// a constant
46 c0 += initCoe * term -> Value ();
47 break;
48
49 case COU_EXPRVAR: /// a variable
50 lmap.insert (term -> Index (), initCoe);
51 break;
52
53 case COU_EXPROPP: /// the opposite of a term
54 decomposeTerm (term -> Argument (), -initCoe, c0, lmap, qmap);
55 break;
56
57 case COU_EXPRSUB: /// a subtraction
58 decomposeTerm (term -> ArgList () [0], initCoe, c0, lmap, qmap);
59 decomposeTerm (term -> ArgList () [1], -initCoe, c0, lmap, qmap);
60 break;
61
62 case COU_EXPRQUAD: { /// a quadratic form
63
64 exprQuad *t = dynamic_cast <exprQuad *> (term -> isaCopy () ?
65 term -> Copy () :
66 term);
67 exprQuad::sparseQ &M = t -> getQ ();
68
69 for (exprQuad::sparseQ::iterator row = M.begin ();
70 row != M.end (); ++row) {
71
72 int xind = row -> first -> Index ();
73
74 for (exprQuad::sparseQcol::iterator col = row -> second.begin ();
75 col != row -> second.end (); ++col) {
76 qmap.insert (xind, col -> first -> Index (), initCoe * col -> second);
77 }
78 }
79 } // NO break here, exprQuad generalizes exprGroup
80
81 case COU_EXPRGROUP: { /// a linear term
82
83 exprGroup *t = dynamic_cast <exprGroup *> (term -> isaCopy () ?
84 term -> Copy () :
85 term);
86 exprGroup::lincoeff &lcoe = t -> lcoeff ();
87
88 // for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
89 for (int n = lcoe.size (), i=0; n--; i++)
90 lmap.insert (lcoe [i].first -> Index (), initCoe * lcoe [i].second);
91
92 c0 += initCoe * t -> getc0 ();
93 } // NO break here, exprGroup generalizes exprSum
94
95 case COU_EXPRSUM: { /// a sum of (possibly) nonlinear elements
96
97 expression **al = term -> ArgList ();
98 for (int i = term -> nArgs (); i--;)
99 decomposeTerm (*al++, initCoe, c0, lmap, qmap);
100
101 } break;
102
103 // not-so-easy cases /////////////////////////////////////////////////////////////////
104 //
105 // cannot add terms as it may fill up the triplet
106
107 case COU_EXPRMUL: { /// a product of n factors /////////////////////////////////////////
108
109 std::map <int, CouNumber> indices;
110 CouNumber coe = initCoe;
111
112 // return list of variables (some of which auxiliary)
113 flattenMul (term, coe, indices);
114
115 if (jnlst_ -> ProduceOutput (Ipopt::J_ALL, J_REFORMULATE)) {
116 printf ("from flattenmul: [%g] ", coe);
117 for (std::map <int, CouNumber>::iterator itt = indices.begin ();
118 itt != indices.end(); ++itt)
119 printf (" %d,%g",
120 itt -> first,
121 itt -> second);
122 printf ("\n");
123 }
124
125 // based on number of factors, decide what to return
126 switch (indices.size ()) {
127
128 case 0: // no variables in multiplication (hmmm...)
129 c0 += coe;
130 break;
131
132 case 1: { // only one term (may be with exponent != 1)
133
134 std::map <int, CouNumber>::iterator one = indices.begin ();
135 int index = one -> first;
136 CouNumber expon = one -> second;
137
138 if (fabs (expon - 1) < COUENNE_EPS) lmap.insert (index, coe);
139 else if (fabs (expon - 2) < COUENNE_EPS) qmap.insert (index, index, coe);
140 else {
141 exprAux *aux = addAuxiliary
142 (new exprPow (new exprClone (Var (index)),
143 new exprConst (expon)));
144
145 //linsert (lmap, aux -> Index (), initCoe); // which of these three is correct?
146 //linsert (lmap, aux -> Index (), initCoe * coe);
147 lmap.insert (aux -> Index (), coe);
148 }
149 } break;
150
151 case 2: { // two terms
152
153 int ind0, ind1;
154
155 std::map <int, CouNumber>::iterator one = indices.begin (),
156 two = one;
157 ++two; // now "two" points to the other variable
158
159 // first variable
160 if (fabs (one -> second - 1) > COUENNE_EPS) {
161 exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (one -> first)),
162 new exprConst (one -> second)));
163 ind0 = aux -> Index ();
164 } else ind0 = one -> first;
165
166 // second variable
167 if (fabs (two -> second - 1) > COUENNE_EPS) {
168 exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (two -> first)),
169 new exprConst (two -> second)));
170 ind1 = aux -> Index ();
171 } else ind1 = two -> first;
172
173 qmap.insert (ind0, ind1, coe);
174 } break;
175
176 default: {
177
178 // create new auxiliary variable containing product of 3+ factors
179
180 expression **al = new expression * [indices.size ()];
181 std::map <int, CouNumber>::iterator one = indices.begin ();
182
183 for (int i=0; one != indices.end (); ++one, i++)
184 if (fabs (one -> second - 1) > COUENNE_EPS) {
185 exprAux *aux = addAuxiliary (new exprPow (new exprClone (Var (one -> first)),
186 new exprConst (one -> second)));
187 al [i] = new exprClone (aux);
188 } else al [i] = new exprClone (Var (one -> first));
189
190 // TODO: when we have a convexification for \prod_{i \in I}...
191 // exprAux *aux = addAuxiliary (new exprMul (al, indices.size ()));
192
193 exprMul *mul = new exprMul (al, indices.size ());
194 exprAux *aux = mul -> standardize (this);
195 lmap.insert (aux -> Index (), coe);
196
197 } break;
198 }
199
200 } break; // end of case COU_EXPRMUL
201
202 case COU_EXPRPOW: { // expression = f(x)^g(x) ////////////////////////////////////////////////
203
204 expression **al = term -> ArgList ();
205
206 if (al [1] -> Type () != CONST) {
207
208 // non-constant exponent, standardize the whole term and add
209 // linear component (single aux)
210
211 expression *aux = term -> standardize (this);
212 if (!aux) aux = term;
213 lmap.insert (aux -> Index (), initCoe);
214
215 } else { // this is of the form f(x)^k. If k=2, return square. If
216 // k=1, return var. Otherwise, generate new auxiliary.
217
218 expression *aux = (*al) -> standardize (this);
219
220 if (!aux)
221 aux = *al; // it was a simple variable, and was not standardized.
222
223 CouNumber expon = al [1] -> Value ();
224 int ind = aux -> Index ();
225
226 if (fabs (expon - 1.) == 0.) lmap.insert (ind, initCoe);
227 else if (fabs (expon - 2.) == 0.) qmap.insert (ind, ind, initCoe);
228 else {
229 exprAux *aux2 = addAuxiliary
230 (new exprPow (new exprClone (aux), new exprConst (expon))); // TODO: FIX!
231 lmap.insert (aux2 -> Index (), initCoe);
232 }
233 }
234 } break;
235
236 default: { /// otherwise, simply standardize expression
237
238 expression *aux = term -> standardize (this);
239 if (!aux)
240 aux = term;
241 lmap.insert (aux -> Index (), initCoe);
242 } break;
243 }
244 }
245