1 /* $Id$
2 *
3 * Name: writeLP.cpp
4 * Author: Pietro Belotti
5 * Purpose: save a problem in LP format (MIQCQP only)
6 *
7 * This file is licensed under the Eclipse Public License (EPL)
8 */
9
10 #include <fstream>
11 #include <iomanip> // to use the setprecision manipulator
12
13 #include "CouenneProblem.hpp"
14 #include "CouenneProblemElem.hpp"
15 #include "CouenneExprVar.hpp"
16
17 #include "CouenneExprGroup.hpp"
18 #include "CouenneExprQuad.hpp"
19
20 using namespace Couenne;
21
22 // print linear or quadratic expression (only used here for LP file
23 // saving). Returns (unprintable) constant terms
24
printLPquad(std::ofstream & f,const expression * ex,double mult)25 double printLPquad (std::ofstream &f, const expression *ex, double mult) {
26
27 expression *mutex = const_cast <expression *> (ex);
28
29 exprGroup *e = dynamic_cast <exprGroup *> (mutex);
30
31 double sign = 1., retval = 0;
32
33 if (!e) {
34 if (mutex -> code () == COU_EXPROPP) {
35 ex = ex -> Argument () -> Original ();
36 mutex = const_cast <expression *> (ex);
37 sign = -1.;
38 e = dynamic_cast <exprGroup *> (mutex);
39 }
40 }
41
42 int wrapCount = 0;
43
44 if (e) {
45
46 // print linear part
47
48 if (e -> getc0 () != 0)
49 retval -= sign * e -> getc0 ();
50
51 for (std::vector <std::pair <exprVar *, CouNumber> >::iterator i = e -> lcoeff () . begin ();
52 i != e -> lcoeff (). end (); ++i) {
53
54 CouNumber c = sign * i -> second;
55
56 if ((c > 0) &&
57 ((i != e -> lcoeff (). begin ()) || (e -> getc0 () != 0)))
58 f << '+';
59
60 f << c << ' ';
61
62 i -> first -> print (f);
63
64 if (!(++wrapCount % 10))
65 f << std::endl;
66 else f << ' ';
67 }
68 } else {
69
70 // assume the expression is a constant, take its value and print
71 // it
72
73 retval -= ex -> Value ();
74 }
75
76 // print quadratic part, if any
77 exprQuad *q = dynamic_cast <exprQuad *> (mutex);
78
79 // todo: an exprGroup has nonlinear components in argList_
80
81 if (q) {
82
83 std::vector <std::pair <exprVar *, exprQuad::sparseQcol> > &Q = q -> getQ ();
84
85 f << " + [ ";
86
87 for (std::vector <std::pair <exprVar *, exprQuad::sparseQcol> >::iterator i = Q. begin ();
88 i != Q . end (); ++i) {
89
90 for (std::vector <std::pair <exprVar *, CouNumber> >::iterator j = i -> second. begin ();
91 j != i -> second . end (); ++j) {
92
93 CouNumber c = mult * sign * j -> second;
94
95 if (c > 0)
96 f << " +";
97
98 f << c << " ";
99
100 if (i -> first -> Index () ==
101 j -> first -> Index ()) {
102 i -> first -> print (f);
103 f << "^2";
104 } else {
105 i -> first -> print (f); f << " * ";
106 j -> first -> print (f);
107 }
108
109 if (!(++wrapCount % 10))
110 f << std::endl;
111 else f << ' ';
112 }
113 }
114
115 f << " ] ";
116 if (mult != 1)
117 f << "/ " << mult;
118 }
119
120 if ((mutex -> nArgs () == 1) &&
121 (mutex -> ArgList () [0] -> Type () == CONST)) {
122
123 retval -= sign * mutex -> ArgList () [0] -> Value ();
124
125 } else {
126
127 f << " + [ ";
128
129 expression *arg = NULL, *arg2 = NULL;
130 double signI;
131
132 for (int i=0; i<mutex -> nArgs ();) {
133
134 if (arg2) {
135
136 arg = arg2;
137 arg2 = NULL;
138 signI = -signI;
139
140 } else {
141
142 signI = sign;
143
144 arg = mutex -> ArgList () [i];
145
146 if (arg -> code () == COU_EXPROPP) {
147
148 arg = arg -> Argument ();
149 signI = - sign;
150 }
151
152 if (arg -> code () == COU_EXPRSUB) {
153
154 arg2 = arg -> ArgList () [1];
155 arg = arg -> ArgList () [0];
156 }
157 }
158
159 if (arg -> code () == COU_EXPROPP) {
160
161 arg = arg -> Argument ();
162 signI = - signI;
163 }
164
165 if (arg -> Type () == CONST) retval -= signI * arg -> Value (); // constant
166
167 else if ((arg -> code () == COU_EXPRPOW) &&
168 (arg -> ArgList () [0] -> Type () == VAR) &&
169 (arg -> ArgList () [1] -> Type () == CONST)) { // xi^2
170
171 f << ((signI > 0) ? " + " : " - ");
172 if (mult != 1) f << mult << " ";
173 arg -> ArgList () [0] -> print (f);
174 f << "^2";
175 }
176
177 else if ((arg -> code () == COU_EXPRMUL) &&
178 (arg -> ArgList () [0] -> Type () == CONST) &&
179 (arg -> ArgList () [1] -> code () == COU_EXPRPOW) &&
180 (arg -> ArgList () [1] -> ArgList () [0] -> Type () == VAR) &&
181 (arg -> ArgList () [1] -> ArgList () [1] -> Type () == CONST)) { // alpha xi^2
182
183 double c = mult * signI * arg -> ArgList () [0] -> Value ();
184 f << ((c > 0) ? '+' : ' ') << c << " ";
185 arg -> ArgList () [1] -> ArgList () [0] -> print (f);
186 f << "^2";
187 }
188
189 else if ((arg -> code () == COU_EXPRMUL) &&
190 (arg -> ArgList () [0] -> Type () == VAR) &&
191 (arg -> ArgList () [1] -> Type () == VAR) &&
192 (arg -> ArgList () [0] -> Index () !=
193 arg -> ArgList () [1] -> Index ())) { // x1 * x2
194
195 f << ((signI > 0) ? " + " : " - ");
196 if (mult != 1) f << mult << " ";
197 arg -> ArgList () [0] -> print (f); f << " * ";
198 arg -> ArgList () [1] -> print (f);
199 }
200
201 else if ((arg -> code () == COU_EXPRMUL) &&
202 (arg -> ArgList () [0] -> Type () == VAR) &&
203 (arg -> ArgList () [1] -> Type () == VAR) &&
204 (arg -> ArgList () [0] -> Index () ==
205 arg -> ArgList () [1] -> Index ())) { // x1 * x1
206
207 f << ((signI > 0) ? " + " : " - ");
208 if (mult != 1) f << mult << " ";
209 arg -> ArgList () [1] -> print (f);
210 f << "^2 ";
211 }
212
213 else if ((arg -> code () == COU_EXPRMUL) &&
214 (arg -> ArgList () [0] -> Type () == CONST) &&
215 (arg -> ArgList () [1] -> code () == COU_EXPRMUL) &&
216 (arg -> ArgList () [1] -> ArgList () [0] -> Type () == VAR) &&
217 (arg -> ArgList () [1] -> ArgList () [1] -> Type () == VAR) &&
218 (arg -> ArgList () [1] -> ArgList () [0] -> Index () !=
219 arg -> ArgList () [1] -> ArgList () [1] -> Index ())) { // alpha x1 * x2
220
221 double c = mult * signI * arg -> ArgList () [0] -> Value ();
222 f << ((c > 0) ? '+' : ' ') << c << " ";
223 arg -> ArgList () [1] -> ArgList () [0] -> print (f); f << " * ";
224 arg -> ArgList () [1] -> ArgList () [1] -> print (f);
225 }
226
227 else if ((arg -> code () == COU_EXPRMUL) &&
228 (arg -> ArgList () [0] -> Type () == CONST) &&
229 (arg -> ArgList () [1] -> code () == COU_EXPRMUL) &&
230 (arg -> ArgList () [1] -> ArgList () [0] -> Type () == VAR) &&
231 (arg -> ArgList () [1] -> ArgList () [1] -> Type () == VAR) &&
232 (arg -> ArgList () [1] -> ArgList () [0] -> Index () !=
233 arg -> ArgList () [1] -> ArgList () [1] -> Index ())) { // alpha x1 * x1
234
235 double c = mult * signI * arg -> ArgList () [0] -> Value ();
236 f << ((c > 0) ? '+' : ' ') << c << " ";
237 arg -> ArgList () [1] -> ArgList () [0] -> print (f);
238 f << "^2 ";
239
240 } else {
241
242 printf ("Can't recognize expression (code: %d), exiting: ", arg -> code ());
243 arg -> print ();
244 printf ("\nExpression: ");
245 ex -> print ();
246 printf ("\n");
247 exit (-1);
248 }
249
250 if (!(++wrapCount % 10))
251 f << std::endl;
252 else f << ' ';
253
254 if (!arg2) ++i;
255 else if (arg -> code () == COU_EXPROPP)
256 signI = - signI;
257 }
258
259 f << " ] ";
260 if (mult != 1)
261 f << "/ " << mult;
262 }
263
264 return retval;
265 }
266
267
268 // store problem in a .mod file (AMPL)
269
writeLP(const std::string & fname)270 void CouenneProblem::writeLP (const std::string &fname) { /// name of the mod file
271
272 // only write LP problems from original.
273
274 for (int i=0; i < nVars (); i++)
275 if (variables_ [i] -> Type () == AUX) {
276 printf ("Auxiliary variables not supported in writeLP yet, bailing out\n");
277 return;
278 }
279
280 if (objectives_ [0] -> Body () -> Linearity () > QUADRATIC) {
281 printf ("Objective is nonlinear and not quadratic, bailing out\n");
282 return;
283 }
284
285 for (int i=0; i < nCons (); i++) {
286 if (constraints_ [i] -> Body () -> Linearity () > QUADRATIC) {
287 printf ("Constraint %d is nonlinear and not quadratic, bailing out\n", i);
288 return;
289 }
290 }
291
292 std::ofstream f (fname.c_str ());
293
294 f << std::setprecision (15);
295
296 // original variables, integer and non //////////////////////////////////////////////
297
298 f << "\\ Problem name (saved using Couenne): " << problemName_ << std::endl << std::endl;
299
300 // objective function /////////////////////////////////////////////////////////////
301
302 //if (objectives_ [0] -> Sense () == MINIMIZE)
303 f << "minimize obj: ";
304
305 double objConst = printLPquad (f, objectives_ [0] -> Body (), 2);
306 if (objConst != 0.) f << ((objConst > 0) ? " + " : " ") << objConst;
307 f << std::endl << std::endl << "Subject To" << std::endl << std::endl;
308
309 // // defined (aux) variables, with formula ///////////////////////////////////////////
310
311 // if (aux) {
312
313 // f << std::endl << "# aux. variables defined" << std::endl << std::endl;
314
315 // for (int i=0; i < nVars (); i++)
316
317 // if ((variables_ [i] -> Type () == AUX) &&
318 // (variables_ [i] -> Multiplicity () > 0)) {
319
320 // f << "aux_" << i << ": "; variables_ [i] -> print (f, false);
321 // f << " = ";
322
323 // variables_ [i] -> Image () -> print (f, false);
324 // f << std::endl;
325 // }
326 // }
327
328 // write constraints //////////////////////////////////////////////////////////////
329
330 // f << std::endl << "# constraints" << std::endl << std::endl;
331
332 // if (!aux) // print h_i(x,y) <= ub, >= lb
333 // for (std::vector <exprVar *>::iterator i = variables_.begin ();
334 // i != variables_.end ();
335 // ++i)
336
337 // if (((*i) -> Type () == AUX) &&
338 // ((*i) -> Multiplicity () > 0)) {
339
340 // CouNumber bound;
341
342 // if ((bound = (*i) -> lb ()) > - COUENNE_INFINITY) {
343 // f << "conAuxLb" << (*i) -> Index () << ": ";
344 // (*i) -> print (f, true);
345 // f << ">= " << bound << std::endl;
346 // }
347
348 // if ((bound = (*i) -> ub ()) < COUENNE_INFINITY) {
349 // f << "conAuxUb" << (*i) -> Index () << ": ";
350 // (*i) -> print (f, true);
351 // f << "<= " << bound << std::endl;
352 // }
353 // }
354
355 for (int i=0; i < nCons (); i++) {
356
357 // get numerical value of lower, upper bound
358 CouNumber lb = (constraints_ [i] -> Lb () -> Value ()),
359 ub = (constraints_ [i] -> Ub () -> Value ());
360
361 f << "con_" << i << ": ";
362 double extraTerm = printLPquad (f, constraints_ [i] -> Body (), 1);
363
364 lb += extraTerm;
365 ub += extraTerm;
366
367 if (lb > - COUENNE_INFINITY) {
368 f << ' ';
369 if (fabs (ub-lb) > COUENNE_EPS)
370 f << '>';
371 f << "= " << lb << std::endl;
372 }
373 else f << " <= " << ub << std::endl;
374
375 // if range constraint, print it once again
376
377 if (( lb > - COUENNE_INFINITY + 1)
378 && (ub < COUENNE_INFINITY - 1)
379 && (fabs (ub-lb) > COUENNE_EPS)) {
380
381 f << "con_" << i << "_rng: ";
382 printLPquad (f, constraints_ [i] -> Body (), 1);
383 f << " <= " << ub << std::endl;
384 }
385 }
386
387 f << "Bounds" << std::endl << std::endl;
388
389 for (int i=0; i < nVars (); i++) {
390
391 if ((Lb (i) == 0) && (Ub (i) >= COUENNE_INFINITY/2))
392 continue;
393
394 if (Lb (i) != 0) f << Lb (i) << " <= ";
395 variables_ [i] -> print (f);
396 if (Ub (i) < + COUENNE_INFINITY/2) f << " <= " << Ub (i);
397 f << std::endl;
398 }
399
400 f << "Generals" << std::endl << std::endl;
401
402 int wrapcount = 0;
403
404 for (int i=0; i < nVars (); i++)
405
406 if (variables_ [i] -> isInteger ()) {
407 variables_ [i] -> print (f);
408 if (!(++wrapcount % 10))
409 f << std::endl;
410 else
411 f << " ";
412 }
413
414 f << std::endl << std::endl << "End" << std::endl;
415
416 f.close ();
417 }
418