1 // $Id$
2 //
3 // (C) Copyright XXX 2009
4 // All Rights Reserved.
5 // This code is published under the Eclipse Public License (EPL).
6 //
7 // Authors :
8 // Pietro Belotti, Lehigh University
9 // Stefan Vigerske, Humboldt University
10 //
11 // Date : 07/18/2009
12
13 #include "CouenneConfig.h"
14 #include "CouenneAmplInterface.hpp"
15 #include "CoinPragma.hpp"
16
17 #include <cstdlib>
18
19 #if defined HAVE_CSTDINT
20 #include <cstdint>
21 #elif defined HAVE_STDINT_H
22 #include <stdint.h>
23 #endif
24
25 #include <string>
26
27 #include "BonAmplTMINLP.hpp"
28 #include "BonCbc.hpp"
29
30 #include "CouenneProblem.hpp"
31 #include "CouenneTypes.hpp"
32
33 #include "CouenneExprClone.hpp"
34 #include "CouenneExprGroup.hpp"
35 #include "CouenneExprAbs.hpp"
36 #include "CouenneExprSum.hpp"
37 #include "CouenneExprSub.hpp"
38 #include "CouenneExprMul.hpp"
39 #include "CouenneExprDiv.hpp"
40 #include "CouenneExprInv.hpp"
41 #include "CouenneExprSin.hpp"
42 #include "CouenneExprPow.hpp"
43 #include "CouenneExprLog.hpp"
44 #include "CouenneExprOpp.hpp"
45 #include "CouenneExprCos.hpp"
46 #include "CouenneExprExp.hpp"
47
48 #include "asl.h"
49 #include "nlp.h"
50 #include "getstub.h"
51 #include "opcode.hd"
52
53 // get ASL op. code relative to function pointer passed as parameter
54 size_t getOperator (efunc *);
55
56 #define OBJ_DE ((const ASL_fg *) asl) -> I.obj_de_
57 #define VAR_E ((const ASL_fg *) asl) -> I.var_e_
58 #define CON_DE ((const ASL_fg *) asl) -> I.con_de_
59 #define OBJ_sense ((const ASL_fg *) asl) -> i.objtype_
60
61 #include "r_opn.hd" /* for N_OPS */
62
63 static fint timing = 0;
64
65 static
66 keyword keywds[] = { // must be alphabetical
67 KW(const_cast<char*>("timing"), L_val, &timing, const_cast<char*>("display timings for the run")),
68 };
69
70 static
71 Option_Info Oinfo = { const_cast<char*>("testampl"), const_cast<char*>("ANALYSIS TEST"),
72 const_cast<char*>("concert_options"), keywds, nkeywds, 0, const_cast<char*>("ANALYSIS TEST") };
73
74
75 // (C++) code starts here ///////////////////////////////////////////////////////////////////////////
76
77 using namespace Couenne;
78
registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)79 void CouenneAmplInterface::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions) {
80 roptions->AddStringOption1("nlfile", "name of an ampl .nl file to get the problem from", "", "*", "name of .nl file");
81 }
82
~CouenneAmplInterface()83 CouenneAmplInterface::~CouenneAmplInterface() {
84 delete problem;
85
86 if (asl) {
87 delete[] X0;
88 delete[] havex0;
89 delete[] pi0;
90 delete[] havepi0;
91 ASL_free(&asl);
92 }
93 }
94
95 // create an AMPL problem by using ASL interface to the .nl file
getCouenneProblem()96 CouenneProblem* CouenneAmplInterface::getCouenneProblem() {
97 if (problem)
98 return problem;
99
100 if (!readASLfg())
101 return NULL;
102
103 problem = new CouenneProblem;
104
105 if (!readnl()) {
106 delete problem;
107 problem = NULL;
108 return NULL;
109 }
110
111 return problem;
112 }
113
getTMINLP()114 Ipopt::SmartPtr<Bonmin::TMINLP> CouenneAmplInterface::getTMINLP() {
115 if (IsValid(tminlp))
116 return tminlp;
117
118 if (IsNull(roptions)) {
119 jnlst->Printf(Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "Error: Need registered options to create AmplTMINLP object!\n");
120 return NULL;
121 }
122
123 std::string nlfile;
124 options->GetStringValue("nlfile", nlfile, "");
125 char** argv = new char*[3];
126 argv[0] = const_cast<char*>("dummy");
127 argv[1] = strdup(nlfile.c_str());
128 argv[2] = NULL;
129 tminlp = new Bonmin::AmplTMINLP(GetRawPtr(jnlst), roptions, options, argv);
130
131 free(argv[1]);
132 delete[] argv;
133
134 return tminlp;
135 }
136
writeSolution(Bonmin::Bab & bab)137 bool CouenneAmplInterface::writeSolution(Bonmin::Bab& bab) {
138 const char* message;
139
140 //TODO setup a nicer message
141 if (bab.bestSolution()) {
142 message = "Couenne found a solution.\n";
143 } else {
144 message = "Couenne could not found a solution.\n";
145 }
146
147 write_sol(const_cast<char*>(message), const_cast<double*>(bab.bestSolution()), NULL, NULL);
148
149 return true;
150 }
151
readASLfg()152 bool CouenneAmplInterface::readASLfg() {
153 assert(asl == NULL);
154
155 std::string nlfile;
156 options->GetStringValue("nlfile", nlfile, "");
157
158 if (nlfile == "")
159 return false;
160
161 char** argv = new char*[3];
162 argv[0] = const_cast<char*>("dummy");
163 argv[1] = strdup(nlfile.c_str());
164 argv[2] = NULL;
165
166 // Create the ASL structure
167 asl = (ASL*) ASL_alloc (ASL_read_fg);
168
169 char* stub = getstub (&argv, &Oinfo);
170
171 // Although very intuitive, we shall explain why the second argument
172 // is passed with a minus sign: it is to tell the ASL to retrieve
173 // the nonlinear information too.
174 FILE* nl = jac0dim (stub, - (fint) strlen (stub));
175
176 // Set options in the asl structure
177 want_xpi0 = 1 | 2; // allocate initial values for primal and dual if available
178 obj_no = 0; // always want to work with the first (and only?) objective
179
180 // allocate space for initial values
181 X0 = new real [n_var];
182 havex0 = new char [n_var];
183 pi0 = new real [n_con];
184 havepi0 = new char [n_con];
185
186 // read the rest of the nl file
187 fg_read (nl, ASL_return_read_err | ASL_findgroups);
188
189 //FIXME freeing argv and argv[1] gives segfault !!!
190 // free(argv[1]);
191 // delete[] argv;
192
193 return true;
194 }
195
196 // check if an expression is a null pointer or equals zero
197 //static inline bool is_expr_zero (expr* e) {
198 // return ((e==NULL) || ((((Intcast (e->op)) == OPNUM) &&
199 // (fabs (((expr_n *)e) -> v) < COUENNE_EPS)
200 // // && (fabs (e -> dL) < COUENNE_EPS)
201 // // *** CHECK THIS! dL is the derivative
202 // )));
203 //}
204
205
206 // Reads a MINLP from an AMPL .nl file through the ASL methods
readnl()207 bool CouenneAmplInterface::readnl() {
208
209 std::string nlfile;
210 options->GetStringValue("nlfile", nlfile, "");
211 problem -> setProblemName (nlfile);
212
213 // number of defined variables (aka common expressions)
214 problem -> setNDefVars(como + comc + comb + como1 + comc1);
215
216 // see "hooking your solver to AMPL", by David M. Gay, tables 3, 4, and 5
217
218 // nonlinear in both objectives and constraints
219 if (nlvb >= 0) {
220 for (int i = 0; i < nlvb - nlvbi; i++) problem -> addVariable (false, problem -> domain ());
221 for (int i = 0; i < nlvbi; i++) problem -> addVariable (true, problem -> domain ());
222 }
223
224 // nonlinear in either objectives or constraints
225 if (nlvo > nlvc) {
226 for (int i = 0; i < nlvc - (nlvb + nlvci); i++) problem -> addVariable (false, problem -> domain ());
227 for (int i = 0; i < nlvci; i++) problem -> addVariable (true, problem -> domain ());
228 for (int i = 0; i < nlvo - (nlvc + nlvoi); i++) problem -> addVariable (false, problem -> domain ());
229 for (int i = 0; i < nlvoi; i++) problem -> addVariable (true, problem -> domain ());
230 } else {
231 for (int i = 0; i < nlvo - (nlvb + nlvoi); i++) problem -> addVariable (false, problem -> domain ());
232 for (int i = 0; i < nlvoi; i++) problem -> addVariable (true, problem -> domain ());
233 for (int i = 0; i < nlvc - (nlvo + nlvci); i++) problem -> addVariable (false, problem -> domain ());
234 for (int i = 0; i < nlvci; i++) problem -> addVariable (true, problem -> domain ());
235 }
236
237 for (int i = 0; i < nwv; i++) problem -> addVariable(false, problem -> domain ());//arc
238 for (int i = n_var - (CoinMax (nlvc,nlvo) +niv+nbv+nwv); i--;) problem -> addVariable(false, problem -> domain ());//other
239 for (int i = 0; i < nbv; i++) problem -> addVariable(true, problem -> domain ());//binary
240 for (int i = 0; i < niv; i++) problem -> addVariable(true, problem -> domain ());//int.
241
242 // add space for common expressions
243 for (int i = problem->nDefVars(); i--;) problem -> addVariable(false, problem -> domain ());
244
245 // common expressions (or defined variables) ///////////////////////////////////////
246
247 #ifdef DEBUG
248 printf ("tot var = %d\n", variables_ . size ());
249 printf ("c_vars_ = %d\n", ((const ASL_fg *) asl) -> i.c_vars_ );
250 printf ("comb_ = %d\n", ((const ASL_fg *) asl) -> i.comb_ );
251 printf ("combc_ = %d\n", ((const ASL_fg *) asl) -> i.combc_ );
252 printf ("comc1_ = %d\n", ((const ASL_fg *) asl) -> i.comc1_ );
253 printf ("comc_ = %d\n", ((const ASL_fg *) asl) -> i.comc_ );
254 printf ("como1_ = %d\n", ((const ASL_fg *) asl) -> i.como1_ );
255 printf ("como_ = %d\n", ((const ASL_fg *) asl) -> i.como_ );
256 #endif
257
258 // Each has a linear and a nonlinear part (thanks to Dominique
259 // Orban: http://www.gerad.ca/~orban/drampl/def-vars.html)
260
261 try {
262 for (int i = 0; i < como + comc + comb; i++) {
263
264 struct cexp *common = ((const ASL_fg *) asl) -> I.cexps_ + i;
265 expression *nle = nl2e (common -> e);
266
267 #ifdef DEBUG
268 printf ("cexp %d [%d]: ", i, problem -> nVars()); nle -> print (); printf (" ||| ");
269 #endif
270
271 int nlin = common -> nlin; // Number of linear terms
272
273 if (nlin > 0) {
274
275 int *indexL = new int [nlin+1];
276 CouNumber *coeff = new CouNumber [nlin];
277
278 linpart *L = common -> L;
279
280 for (int j = 0; j < nlin; j++) {
281 //vp = (expr_v *)((char *)L->v.rp - ((char *)&ev.v - (char *)&ev));
282 //Printf( " %-g x[%-d]", L->fac, (int)(vp - VAR_E) );
283 coeff [j] = L [j]. fac;
284 indexL [j] = ((uintptr_t) (L [j].v.rp) - (uintptr_t) VAR_E) / sizeof (expr_v);
285 #ifdef DEBUG
286 Printf( " %+g x_%-3d", L [j]. fac,
287 (expr_v *) (L [j].v.rp) - VAR_E //((const ASL_fg *) asl) -> I.cexps_
288 //L [j]. v.i
289 );
290 #endif
291 }
292
293 indexL [nlin] = -1;
294
295 expression **al = new expression * [1];
296 *al = nle;
297
298 std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
299 problem -> indcoe2vector (indexL, coeff, lcoeff);
300
301 expression *eg = exprGroup::genExprGroup (0, lcoeff, al, 1);
302 problem -> commonExprs (). push_back (eg);
303 }
304 else problem -> commonExprs () . push_back (nle);
305 #ifdef DEBUG
306 printf ("\n");
307 #endif
308 }
309
310 for (int i = 0; i < como1 + comc1; i++) {
311
312 struct cexp1 *common = ((const ASL_fg *) asl) -> I.cexps1_ + i;
313 expression *nle = nl2e (common -> e);
314
315 #ifdef DEBUG
316 printf ("cexp1 %d [%d]: ", i, variables_ . size ()); nle -> print (); printf (" ||| ");
317 #endif
318
319 int nlin = common -> nlin; // Number of linear terms
320
321 if (nlin > 0) {
322
323 int *indexL = new int [nlin+1];
324 CouNumber *coeff = new CouNumber [nlin];
325
326 linpart *L = common -> L;
327
328 for (int j = 0; j < nlin; j++) {
329 //vp = (expr_v *)((char *)L->v.rp - ((char *)&ev.v - (char *)&ev));
330 coeff [j] = L [j]. fac;
331 indexL [j] = ((uintptr_t) (L [j].v.rp) - (uintptr_t) VAR_E) / sizeof (expr_v);
332 #ifdef DEBUG
333 Printf( " %+g x_%-3d", L [j]. fac,
334 (expr_v *) (L [j].v.rp) - VAR_E //((const ASL_fg *) asl) -> I.cexps_
335 //L [j]. v.i
336 );
337 #endif
338 }
339
340 indexL [nlin] = -1;
341
342 expression **al = new expression * [1];
343 *al = nle;
344
345 std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
346 problem -> indcoe2vector (indexL, coeff, lcoeff);
347
348 expression *eg = exprGroup::genExprGroup (0, lcoeff, al, 1);
349 problem -> commonExprs () . push_back (eg);
350 }
351 else problem -> commonExprs () . push_back (nle);
352 #ifdef DEBUG
353 printf ("\n");
354 #endif
355 // addAuxiliary (nl2e (((const ASL_fg *) asl) -> I.cexps1_ [i] . e));
356 }
357
358 // objective functions /////////////////////////////////////////////////////////////
359
360 for (int i = 0; i < n_obj; i++) {
361
362 ////////////////////////////////////////////////
363 int nterms = 0;
364
365 // count nonzero terms in linear part
366
367 for (ograd *objgrad = Ograd [i];
368 objgrad;
369 objgrad = objgrad -> next)
370 if (fabs (objgrad -> coef) > COUENNE_EPS)
371 nterms++;
372
373 expression
374 *body,
375 *nl = nl2e (OBJ_DE [i] . e);
376
377 if (nterms) { // have linear terms
378
379 int *indexL = new int [nterms+1];
380 CouNumber *coeff = new CouNumber [nterms];
381
382 for (ograd *objgrad = Ograd [i]; objgrad; objgrad = objgrad -> next)
383 if (fabs (objgrad -> coef) > COUENNE_EPS) {
384
385 *indexL++ = objgrad -> varno;
386 *coeff++ = objgrad -> coef;
387 }
388
389 *indexL = -1;
390
391 indexL -= nterms;
392 coeff -= nterms;
393
394 std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
395 problem -> indcoe2vector (indexL, coeff, lcoeff);
396
397 if (nl -> code () == COU_EXPRSUM) {
398 body = exprGroup::genExprGroup (0., lcoeff, nl -> ArgList (), nl -> nArgs ());
399 // delete node without deleting children (they are now in body)
400 nl -> ArgList (NULL);
401 delete nl;
402 }
403 else {
404
405 expression **nll = new expression * [1];
406
407 *nll = nl;
408
409 // apparently, objconst (i) is included in the obj expression
410 body = exprGroup::genExprGroup (0., lcoeff, nll, 1);
411 //body = new exprGroup (objconst (i), indexL, coeff, nll, 1);
412 }
413
414 delete [] indexL;
415 delete [] coeff;
416
417 } else
418 // apparently, objconst (i) is included in the obj expression
419 body = nl;
420 //if (fabs (objconst (i) > COUENNE_EPS))
421 //body = new exprSum (nl, new exprConst (objconst (i)));
422 //else
423
424 ///////////////////////////////////////////////////
425
426 expression *subst = Simplified (body);// -> simplify ();
427
428 // if (subst) {
429 // delete body; // VALGRIND
430 // body = subst;
431 // }
432
433 // ThirdParty/ASL/solvers/asl.h, line 336: 0 is minimization, 1 is maximization
434 problem -> addObjective (body, (OBJ_sense [i] == 0) ? "min" : "max");
435 }
436
437 // constraints ///////////////////////////////////////////////////////////////////
438
439 int *nterms = new int [n_con];
440
441 // allocate space for argument list of all constraints' summations
442 // of linear and nonlinear terms
443
444 // init array with # terms of each constraint
445 for (int i = n_con; i--;)
446 *nterms++ = 0;
447 nterms -= n_con;
448
449 cgrad *congrad;
450
451 // count all linear terms
452 if (A_colstarts && A_vals) // Constraints' linear info is stored in A_vals
453 for (register int j = A_colstarts [n_var]; j--;) {
454
455 real coeff = A_vals [j];
456
457 if (fabs (coeff) > COUENNE_EPS)
458 nterms [A_rownos [j]] ++;
459 }
460 else { // Constraints' linear info is stored in Cgrad
461 for (register int i = 0; i < n_con; i++)
462 for (congrad = Cgrad [i];
463 congrad;
464 congrad = congrad -> next)
465 if (fabs (congrad -> coef) > COUENNE_EPS)
466 nterms [i] ++;
467 }
468
469
470 // vectors of the linear part
471 CouNumber **coeff = new CouNumber * [n_con];
472 int **indexL = new int * [n_con];
473
474 for (register int i = n_con; i--;)
475 *indexL++ = NULL;
476
477 indexL -= n_con;
478
479
480 // set linear terms
481
482 if (A_colstarts && A_vals) // Constraints' linear info is stored in A_vals
483 for (int j = 0; j < n_var; j++)
484 for (register int i = A_colstarts [j], k = A_colstarts [j+1] - i; k--; i++) {
485
486 int rowno = A_rownos [i],
487 nt = nterms [rowno] --;
488
489 CouNumber **cline = coeff + rowno;
490 int **iline = indexL + rowno;
491
492 if (*iline==NULL) {
493 *cline = new CouNumber [nt];
494 *iline = new int [nt+1];
495 (*iline) [nt] = -1;
496 }
497
498 (*cline) [--nt] = A_vals [i];
499 (*iline) [nt] = j;
500
501 }
502 else { // Constraints' linear info is stored in Cgrad
503 for (int i=0; i < n_con; i++) {
504
505 int nt = nterms [i];
506
507 CouNumber **cline = coeff + i;
508 int **iline = indexL + i;
509
510 *cline = new CouNumber [nt];
511 *iline = new int [nt+1];
512 (*iline) [nt] = -1;
513
514 for (congrad = Cgrad [i]; congrad; congrad = congrad -> next)
515 if (fabs (congrad -> coef) > COUENNE_EPS) {
516 (*cline) [--nt] = congrad -> coef;
517 (*iline) [nt] = congrad -> varno;
518 }
519 }
520 }
521
522 // set constraints' bound and sign and store nonlinear part ///////////////////////////////
523
524 for (int i = 0; i < n_con; i++) {
525
526 enum con_sign sign;
527 double lb, ub;
528
529 if (Urhsx) {
530 lb = LUrhs [i];
531 ub = Urhsx [i];
532 } else {
533 int j = 2*i;
534 lb = LUrhs [j];
535 ub = LUrhs [j+1];
536 }
537
538 // set constraint sign
539 if (lb > negInfinity)
540 if (ub < Infinity) sign = COUENNE_RNG;
541 else sign = COUENNE_GE;
542 else sign = COUENNE_LE;
543
544 // this is an equality constraint
545 if (fabs (lb - ub) < COUENNE_EPS)
546 sign = COUENNE_EQ;
547
548 expression *body;
549
550 expression **nll = new expression * [1];
551 *nll = nl2e (CON_DE [i] . e);
552
553 if (indexL [i] && (*(indexL [i]) >= 0)) {
554
555 int code = (*nll) -> code ();
556
557 std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
558 problem -> indcoe2vector (indexL [i], coeff [i], lcoeff);
559
560 /*std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
561 for (int i=0, *ind = indexL; *ind >= 0; *ind++, i++)
562 lcoeff.push_back (std::pair <exprVar *, CouNumber> (Var (*ind), coeff [i]));*/
563
564 if ((code == COU_EXPRSUM) ||
565 (code == COU_EXPRGROUP)) {
566
567 body = exprGroup::genExprGroup (0., lcoeff, (*nll) -> ArgList (), (*nll) -> nArgs ());
568 // delete node without deleting children (they are now in body)
569 (*nll) -> ArgList (NULL);
570 delete *nll;
571 delete [] nll;
572 }
573 else body = exprGroup::genExprGroup (0., lcoeff, nll, 1);
574 }
575 else {
576 body = *nll;
577 delete [] nll;
578 }
579
580 expression *subst = Simplified (body);
581
582 // -> simplify ();
583 // if (subst) {
584 // delete body; // VALGRIND
585 // body = subst;
586 // }
587
588 // add them (and set lower-upper bound)
589 switch (sign) {
590
591 case COUENNE_EQ: problem -> addEQConstraint (body, new exprConst (ub)); break;
592 case COUENNE_LE: problem -> addLEConstraint (body, new exprConst (ub)); break;
593 case COUENNE_GE: problem -> addGEConstraint (body, new exprConst (lb)); break;
594 case COUENNE_RNG: problem -> addRNGConstraint (body, new exprConst (lb),
595 new exprConst (ub)); break;
596 default: jnlst->Printf(Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "Error: could not recognize constraint\n"); return false;
597 }
598
599 delete [] indexL [i];
600 delete [] coeff [i];
601 }
602
603 delete [] indexL;
604 delete [] coeff;
605 delete [] nterms;
606
607 } catch (...) {
608 return false;
609 }
610
611 // create room for problem's variables and bounds
612 CouNumber
613 *x = (CouNumber *) malloc ((n_var + problem -> nDefVars() ) * sizeof (CouNumber)),
614 *lb = (CouNumber *) malloc ((n_var + problem -> nDefVars() ) * sizeof (CouNumber)),
615 *ub = (CouNumber *) malloc ((n_var + problem -> nDefVars() ) * sizeof (CouNumber));
616
617 for (int i = n_var + problem -> nDefVars(); i--;) {
618 x [i] = 0.;
619 lb [i] = -COUENNE_INFINITY;
620 ub [i] = COUENNE_INFINITY;
621 }
622
623 problem -> domain () -> push (n_var + problem -> nDefVars(), x, lb, ub);
624 free (x); free (lb); free (ub);
625
626 // suggested:
627 // problem -> domain () -> push (n_var + problem -> nDefVars(), x, lb, ub, false);
628 // //free (x); free (lb); free (ub);
629 // saves three allocations (default last parameter is true, which copies x,l,b)
630
631 // lower and upper bounds ///////////////////////////////////////////////////////////////
632
633 if (LUv) {
634
635 real *Uvx_copy = Uvx;
636
637 if (!Uvx_copy)
638 for (register int i=0; i<n_var; i++) {
639
640 register int j = 2*i;
641
642 problem -> Lb (i) = LUv[j] <= -COUENNE_INFINITY ? -COUENNE_INFINITY : LUv[j] ;
643 problem -> Ub (i) = LUv[j+1] >= COUENNE_INFINITY ? COUENNE_INFINITY : LUv[j+1];
644 }
645 else
646 for (register int i=n_var; i--;) {
647 problem -> Lb (i) = LUv [i] <= -COUENNE_INFINITY ? -COUENNE_INFINITY : LUv[i];
648 problem -> Ub (i) = Uvx_copy [i] >= COUENNE_INFINITY ? COUENNE_INFINITY : Uvx_copy[i];
649 }
650
651 } else
652 for (register int i=n_var; i--;) {
653 problem -> Lb (i) = - COUENNE_INFINITY;
654 problem -> Ub (i) = COUENNE_INFINITY;
655 }
656
657 // initial x ////////////////////////////////////////////////////////////////////
658
659 for (register int i=n_var; i--;)
660
661 if (X0 && havex0 [i]) problem -> X (i) = X0 [i];
662
663 else {
664
665 CouNumber x, l = problem -> Lb (i), u = problem -> Ub (i);
666
667 if (l < - COUENNE_INFINITY)
668 if (u > COUENNE_INFINITY) x = 0.;
669 else x = u;
670 else if (u > COUENNE_INFINITY) x = l;
671 else x = 0.5 * (l+u);
672
673 problem -> X (i) = x;
674 }
675
676 for (register int i=n_var; i < problem -> nDefVars() ; i++) { //FIXME: shouldn't this loop go until n_var + problem -> nDefVars() ?
677
678 problem -> X (i) = 0.;
679 problem -> Lb (i) = -COUENNE_INFINITY;
680 problem -> Ub (i) = COUENNE_INFINITY;
681 }
682
683 return true;
684 }
685
686
687 // warning for non-implemented functions -- return 0 constant expression
688 //expression *notimpl (const std::string &fname) {
689 //static void notimpl (const std::string &fname) {
690 // std::cerr << "*** Error: " << fname << " not implemented" << std::endl;
691 // exit (-1);
692 //}
693
694 // converts an AMPL expression (sub)tree into an expression* (sub)tree
nl2e(expr * e)695 expression *CouenneAmplInterface::nl2e(expr *e) {
696
697 switch (getOperator (e -> op)) {
698
699 case OPPLUS: return new exprSum (nl2e (e -> L.e), nl2e (e -> R.e));
700 case OPMINUS: return new exprSub (nl2e (e -> L.e), nl2e (e -> R.e));
701 case OPMULT: return new exprMul (nl2e (e -> L.e), nl2e (e -> R.e));
702 case OPDIV: return new exprDiv (nl2e (e -> L.e), nl2e (e -> R.e));
703 //case OPREM: notimpl ("remainder");
704 case OPPOW: return new exprPow (nl2e (e -> L.e), nl2e (e -> R.e));
705 //case OPLESS: notimpl ("less");
706 //case MINLIST: notimpl ("min");
707 //case MAXLIST: notimpl ("max");
708 //case FLOOR: notimpl ("floor");
709 //case CEIL: notimpl ("ceil");
710 case ABS: return new exprAbs (nl2e (e -> L.e));
711 case OPUMINUS:return new exprOpp (nl2e (e -> L.e));
712 // return new exprOpp (nl2e (e -> L.e -> L.e));
713 //case OPIFnl: { notimpl ("ifnl");
714
715 // see ASL/solvers/rops.c, IfNL
716 //}
717
718 case OP_tanh: return new exprDiv
719 (new exprSub (new exprExp (nl2e (e -> L.e)),
720 new exprExp (new exprOpp (nl2e (e->L.e)))),
721 new exprSum (new exprExp (nl2e (e -> L.e)),
722 new exprExp (new exprOpp (nl2e (e->L.e)))));
723
724 case OP_tan: {
725 expression *arg;
726 arg = nl2e (e -> L.e);
727 return new exprDiv (new exprSin (arg), new exprCos (new exprClone (arg)));
728 }
729 case OP_sqrt: return new exprPow (nl2e (e -> L.e), new exprConst (0.5));
730 case OP_sinh: return new exprMul (new exprConst (0.5),
731 new exprSub (new exprExp (nl2e (e -> L.e)),
732 new exprExp (new exprOpp (nl2e (e->L.e)))));
733 case OP_sin: return new exprSin (nl2e (e -> L.e));
734 case OP_log10: return new exprMul (new exprConst (1.0 / log (10.0)),
735 new exprLog (nl2e (e -> L.e)));
736 case OP_log: return new exprLog (nl2e (e -> L.e));
737 case OP_exp: return new exprExp (nl2e (e -> L.e));
738 case OP_cosh: return new exprMul (new exprConst (0.5),
739 new exprSum (new exprExp (nl2e (e -> L.e)),
740 new exprExp (new exprOpp (nl2e (e->L.e)))));
741
742 case OP_cos: return new exprCos (nl2e (e -> L.e));
743 //case OP_atanh: notimpl ("atanh");
744 //case OP_atan2: notimpl ("atan2");
745 //case OP_atan: notimpl ("atan");
746 //case OP_asinh: notimpl ("asinh");
747 //case OP_asin: notimpl ("asin");
748 //case OP_acosh: notimpl ("acosh");
749 //case OP_acos: notimpl ("acos");
750
751 case OPSUMLIST: {
752 int i=0;
753 expression **al = new expression * [(e->R.ep - e->L.ep)];
754 for (expr **ep = e->L.ep; ep < e->R.ep; ep++)
755 al [i++] = nl2e (*ep);
756 return new exprSum (al, i);
757 }
758 //case OPintDIV: notimpl ("intdiv");
759 //case OPprecision: notimpl ("precision");
760 //case OPround: notimpl ("round");
761 //case OPtrunc: notimpl ("trunc");
762
763 case OP1POW: return new exprPow (nl2e (e -> L.e),
764 new exprConst (((expr_n *)e->R.e)->v));
765 case OP2POW: return new exprPow (nl2e (e -> L.e),
766 new exprConst (2.));
767 case OPCPOW: return new exprPow (new exprConst (((expr_n *)e->L.e)->v),
768 nl2e (e -> R.e));
769 //case OPFUNCALL: notimpl ("function call");
770 case OPNUM: return new exprConst (((expr_n *)e)->v);
771 //case OPPLTERM: notimpl ("plterm");
772 //case OPIFSYM: notimpl ("ifsym");
773 //case OPHOL: notimpl ("hol");
774 case OPVARVAL: {
775
776 int j = ((expr_v *) e) -> a;
777
778 if (j >= problem -> nOrigVars()) // common expression
779 // use base pointer otherwise the .a field returns an awkward, out-of-bound index
780 j = ((expr_v *) e) - ((const ASL_fg *) asl) -> I.var_e_;
781
782 if (j >= problem -> nOrigVars() + problem -> nDefVars()) {
783 jnlst -> Printf (Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "Error: unknown variable x_%d\n", j);
784 throw -1;
785 }
786
787 return new exprClone (problem -> Variables() [j]);
788 }
789
790 default:
791 jnlst -> Printf (Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "ERROR: unknown operator (address %p), aborting.\n", Intcast (e -> op));
792 throw -2;
793 }
794
795 return new exprConst (0.);
796 }
797