source: stable/0.3/Couenne/src/readnl/readnl.cpp @ 578

Last change on this file since 578 was 578, checked in by pbelotti, 9 years ago

added fictitious objective when none is declared in the .nl. Fixing branch on large value with very negative lower and very large upper bounds

  • Property svn:keywords set to Id
File size: 14.4 KB
Line 
1/* $Id: readnl.cpp 578 2011-05-21 21:03:04Z pbelotti $
2 *
3 * Name:    readnl.cpp
4 * Author:  Pietro Belotti
5 * Purpose: define a reader for .nl files. Adapted from ampl2ev3 by L. Liberti and S. Galli
6 *
7 * (C) Carnegie-Mellon University, 2006.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include "CouenneProblem.hpp"
12
13#if   defined HAVE_CSTDINT
14#include "cstdint"
15#elif defined HAVE_STDINT_H
16#include "stdint.h"
17#endif
18
19#include "CouenneTypes.hpp"
20
21#include "exprSum.hpp"
22#include "exprMul.hpp"
23#include "exprClone.hpp"
24#include "exprGroup.hpp"
25
26#include "asl.h"
27#include "nlp.h"
28#include "getstub.h"
29#include "opcode.hd"
30
31#define OBJ_DE    ((const ASL_fg *) asl) -> I.obj_de_
32#define VAR_E     ((const ASL_fg *) asl) -> I.var_e_
33#define CON_DE    ((const ASL_fg *) asl) -> I.con_de_
34#define OBJ_sense ((const ASL_fg *) asl) -> i.objtype_
35
36//#define DEBUG
37
38// check if an expression is a null pointer or equals zero
39inline bool is_expr_zero (expr* e)
40{return ((e==NULL) || ((((Intcast (e->op)) == OPNUM) && 
41                          (fabs (((expr_n *)e) -> v)  < COUENNE_EPS) 
42                          //  && (fabs (e -> dL) < COUENNE_EPS)
43                          // *** CHECK THIS! dL is the derivative
44                          )));} 
45
46// Reads a MINLP from an AMPL .nl file through the ASL methods
47int CouenneProblem::readnl (const ASL *asl) {
48
49  problemName_ = filename;
50
51  // number of defined variables (aka common expressions)
52  ndefined_ = como + comc + comb + como1 + comc1; 
53
54  // see "hooking your solver to AMPL", by David M. Gay, tables 3, 4, and 5
55
56  // nonlinear in both objectives and constraints
57  if (nlvb >= 0) {
58    for (int i = 0; i < nlvb - nlvbi; i++) addVariable (false, &domain_);
59    for (int i = 0; i < nlvbi;        i++) addVariable (true,  &domain_);
60  }
61
62  // nonlinear in either objectives or constraints
63  if (nlvo > nlvc) {
64    for (int i = 0; i < nlvc - (nlvb + nlvci); i++) addVariable (false, &domain_);
65    for (int i = 0; i < nlvci;                 i++) addVariable (true,  &domain_);
66    for (int i = 0; i < nlvo - (nlvc + nlvoi); i++) addVariable (false, &domain_);
67    for (int i = 0; i < nlvoi;                 i++) addVariable (true,  &domain_);
68  } else {
69    for (int i = 0; i < nlvo - (nlvb + nlvoi); i++) addVariable (false, &domain_);
70    for (int i = 0; i < nlvoi;                 i++) addVariable (true,  &domain_);
71    for (int i = 0; i < nlvc - (nlvo + nlvci); i++) addVariable (false, &domain_);
72    for (int i = 0; i < nlvci;                 i++) addVariable (true,  &domain_);
73  }
74
75  for (int i = 0; i < nwv; i++)                                  addVariable(false, &domain_);//arc
76  for (int i = n_var - (CoinMax (nlvc,nlvo) +niv+nbv+nwv); i--;) addVariable(false, &domain_);//other
77  for (int i = 0; i < nbv; i++)                                  addVariable(true,  &domain_);//binary
78  for (int i = 0; i < niv; i++)                                  addVariable(true,  &domain_);//int.
79
80  // add space for common expressions
81  for (int i = ndefined_; i--;)                                  addVariable(false, &domain_);
82
83  // common expressions (or defined variables) ///////////////////////////////////////
84
85#ifdef DEBUG
86  printf ("tot var = %d\n", variables_ . size ());
87  printf ("c_vars_ = %d\n", ((const ASL_fg *) asl) -> i.c_vars_ );
88  printf ("comb_ = %d\n",   ((const ASL_fg *) asl) -> i.comb_  );
89  printf ("combc_ = %d\n",  ((const ASL_fg *) asl) -> i.combc_ );
90  printf ("comc1_ = %d\n",  ((const ASL_fg *) asl) -> i.comc1_ );
91  printf ("comc_ = %d\n",   ((const ASL_fg *) asl) -> i.comc_  );
92  printf ("como1_ = %d\n",  ((const ASL_fg *) asl) -> i.como1_ );
93  printf ("como_ = %d\n",   ((const ASL_fg *) asl) -> i.como_  );
94#endif
95
96  // Each has a linear and a nonlinear part (thanks to Dominique
97  // Orban: http://www.gerad.ca/~orban/drampl/def-vars.html)
98
99  for (int i = 0; i < como + comc + comb; i++) {
100
101    struct cexp *common = ((const ASL_fg *) asl) -> I.cexps_ + i;
102    expression *nle = nl2e (common -> e, asl);
103
104#ifdef DEBUG
105    printf ("cexp  %d [%d]: ", i, variables_ . size ()); nle -> print ();  printf (" ||| ");
106#endif
107
108    int nlin = common -> nlin;  // Number of linear terms
109
110    if (nlin > 0) {
111
112      int       *indexL = new int       [nlin+1];
113      CouNumber *coeff  = new CouNumber [nlin];
114
115      linpart *L = common -> L;
116
117      for (int j = 0; j < nlin; j++) {
118        //vp = (expr_v *)((char *)L->v.rp - ((char *)&ev.v - (char *)&ev));
119        //Printf( " %-g x[%-d]", L->fac, (int)(vp - VAR_E) );   
120        coeff [j] = L [j]. fac;
121        indexL [j] = ((uintptr_t) (L [j].v.rp) - (uintptr_t) VAR_E) / sizeof (expr_v);
122#ifdef DEBUG
123        Printf( " %+g x_%-3d", L [j]. fac, 
124                (expr_v *) (L [j].v.rp) - VAR_E //((const ASL_fg *) asl) -> I.cexps_
125                //L [j]. v.i
126                );
127#endif
128      }
129
130      indexL [nlin] = -1;
131
132      expression **al = new expression * [1];
133      *al = nle;
134
135      std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
136      indcoe2vector (indexL, coeff, lcoeff);
137
138      expression *eg = exprGroup::genExprGroup (0, lcoeff, al, 1);
139      commonexprs_ . push_back (eg);
140    } 
141    else commonexprs_ . push_back (nle);
142#ifdef DEBUG
143    printf ("\n");
144#endif
145  }
146
147  for (int i = 0; i < como1 + comc1; i++) {
148
149    struct cexp1 *common = ((const ASL_fg *) asl) -> I.cexps1_ + i;
150    expression *nle = nl2e (common -> e, asl);
151
152#ifdef DEBUG
153    printf ("cexp1 %d [%d]: ", i, variables_ . size ()); nle -> print ();  printf (" ||| ");
154#endif
155
156    int nlin = common -> nlin;  // Number of linear terms
157
158    if (nlin > 0) {
159
160      int       *indexL = new int       [nlin+1];
161      CouNumber *coeff  = new CouNumber [nlin];
162
163      linpart *L = common -> L;
164
165      for (int j = 0; j < nlin; j++) {
166        //vp = (expr_v *)((char *)L->v.rp - ((char *)&ev.v - (char *)&ev));
167        coeff  [j] = L [j]. fac;
168        indexL [j] = ((uintptr_t) (L [j].v.rp) - (uintptr_t) VAR_E) / sizeof (expr_v);
169#ifdef DEBUG
170        Printf( " %+g x_%-3d", L [j]. fac, 
171                (expr_v *) (L [j].v.rp) - VAR_E //((const ASL_fg *) asl) -> I.cexps_
172                //L [j]. v.i
173                );
174#endif
175      }
176
177      indexL [nlin] = -1;
178
179      expression **al = new expression * [1];
180      *al = nle;
181
182      std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
183      indcoe2vector (indexL, coeff, lcoeff);
184
185      expression *eg = exprGroup::genExprGroup (0, lcoeff, al, 1);
186      commonexprs_ . push_back (eg);
187    } 
188    else commonexprs_ . push_back (nle);
189#ifdef DEBUG
190    printf ("\n");
191#endif
192    //    addAuxiliary (nl2e (((const ASL_fg *) asl) -> I.cexps1_ [i] . e, asl));
193  }
194
195  // objective functions /////////////////////////////////////////////////////////////
196
197  if (n_obj == 0) {
198
199    // strange, no objective function. Add one equal to zero
200
201    jnlst_ -> Printf (Ipopt::J_ERROR, J_PROBLEM, "Couenne: warning, no objective function found\nAdded fictitious function f(x)=0\n");     
202    addObjective (new exprConst (0.), "min");
203  }
204
205  for (int i = 0; i < n_obj; i++) {
206
207    ////////////////////////////////////////////////
208    int nterms = 0;
209
210    // count nonzero terms in linear part
211 
212    for (ograd *objgrad = Ograd [i];
213         objgrad;
214         objgrad = objgrad -> next)
215      if (fabs (objgrad -> coef) > COUENNE_EPS)
216        nterms++;
217
218    expression
219      *body,
220      *nl = nl2e (OBJ_DE [i] . e, asl);
221
222    if (nterms) { // have linear terms
223
224      int       *indexL = new int       [nterms+1];
225      CouNumber *coeff  = new CouNumber [nterms];
226
227      for (ograd *objgrad = Ograd [i]; objgrad; objgrad = objgrad -> next)
228        if (fabs (objgrad -> coef) > COUENNE_EPS) {
229
230          *indexL++ = objgrad -> varno;
231          *coeff++  = objgrad -> coef;
232        }
233
234      *indexL = -1;
235
236      indexL -= nterms;
237      coeff  -= nterms;
238
239      std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
240      indcoe2vector (indexL, coeff, lcoeff);
241
242      if (nl -> code () == COU_EXPRSUM) {
243        body = exprGroup::genExprGroup (0., lcoeff, nl -> ArgList (), nl -> nArgs ());
244        // delete node without deleting children (they are now in body)
245        nl -> ArgList (NULL);
246        delete nl;
247      }
248      else {
249
250        expression **nll = new expression * [1];
251
252        *nll = nl;
253
254        // apparently, objconst (i) is included in the obj expression
255        body = exprGroup::genExprGroup (0., lcoeff, nll, 1);
256        //body = new exprGroup (objconst (i), indexL, coeff, nll, 1);
257      }
258
259      delete [] indexL;
260      delete [] coeff;
261
262    } else
263      // apparently, objconst (i) is included in the obj expression
264      body = nl;
265      //if (fabs (objconst (i) > COUENNE_EPS))
266      //body = new exprSum (nl, new exprConst (objconst (i)));
267      //else
268
269    ///////////////////////////////////////////////////
270
271    expression *subst = body -> simplify ();
272
273    if (subst) {
274      delete body; // VALGRIND
275      body = subst;
276    }
277
278    // ThirdParty/ASL/solvers/asl.h, line 336: 0 is minimization, 1 is maximization
279    addObjective (body, (OBJ_sense [i] == 0) ? "min" : "max");
280  }
281
282  // constraints ///////////////////////////////////////////////////////////////////
283
284  int *nterms = new int [n_con];
285
286  // allocate space for argument list of all constraints' summations
287  // of linear and nonlinear terms
288
289  // init array with # terms of each constraint
290  for (int i = n_con; i--;) 
291    *nterms++ = 0;
292  nterms -= n_con;
293
294  cgrad *congrad;
295
296  // count all linear terms
297  if (A_colstarts && A_vals)         // Constraints' linear info is stored in A_vals
298    for (register int j = A_colstarts [n_var]; j--;) {
299
300      real coeff = A_vals [j];
301
302      if (fabs (coeff) > COUENNE_EPS)
303        nterms [A_rownos [j]] ++;
304    }
305  else {                             // Constraints' linear info is stored in Cgrad
306    for (register int i = 0; i < n_con; i++)
307      for (congrad = Cgrad [i]; 
308           congrad; 
309           congrad = congrad -> next) 
310        if (fabs (congrad -> coef) > COUENNE_EPS) 
311          nterms [i] ++;
312  }
313
314
315  // vectors of the linear part
316  CouNumber **coeff  = new CouNumber * [n_con];
317  int       **indexL = new int       * [n_con];
318
319  for (register int i = n_con; i--;) 
320    *indexL++ = NULL;
321
322  indexL -= n_con;
323
324
325  // set linear terms
326
327  if (A_colstarts && A_vals)         // Constraints' linear info is stored in A_vals
328    for (int j = 0; j < n_var; j++)
329      for (register int i = A_colstarts [j], k = A_colstarts [j+1] - i; k--; i++) {
330
331        int rowno = A_rownos [i],
332            nt    = nterms [rowno] --;
333
334        CouNumber **cline = coeff  + rowno;
335        int       **iline = indexL + rowno;
336
337        if (*iline==NULL) {
338          *cline = new CouNumber [nt];
339          *iline = new int       [nt+1];
340          (*iline) [nt] = -1;
341        }
342
343        (*cline) [--nt] = A_vals [i];
344        (*iline)   [nt] = j;
345
346      }
347  else {                             // Constraints' linear info is stored in Cgrad
348    for (int i=0; i < n_con; i++) {
349
350      int nt = nterms [i];
351
352      CouNumber **cline = coeff + i;
353      int       **iline = indexL + i;
354
355      *cline = new CouNumber [nt];
356      *iline = new int       [nt+1];
357      (*iline) [nt] = -1;
358
359      for (congrad = Cgrad [i]; congrad; congrad = congrad -> next) 
360        if (fabs (congrad -> coef) > COUENNE_EPS) {
361          (*cline) [--nt] = congrad -> coef;
362          (*iline)   [nt] = congrad -> varno;
363        }
364    }
365  }
366
367  // set constraints' bound and sign and store nonlinear part ///////////////////////////////
368
369  for (int i = 0; i < n_con; i++) {
370
371    enum con_sign sign;
372    double lb, ub;
373
374    if (Urhsx) {
375      lb = LUrhs [i];
376      ub = Urhsx [i];
377    } else {
378      int j = 2*i;
379      lb = LUrhs [j];
380      ub = LUrhs [j+1];
381    }
382
383    // set constraint sign
384    if (lb > negInfinity)
385      if (ub < Infinity) sign = COUENNE_RNG;
386      else               sign = COUENNE_GE;
387    else                 sign = COUENNE_LE;
388
389    // this is an equality constraint 
390    if (fabs (lb - ub) < COUENNE_EPS)
391      sign = COUENNE_EQ;
392
393    expression *body;
394
395    expression **nll = new expression * [1];
396    *nll = nl2e (CON_DE [i] . e, asl);
397
398    if (indexL [i] && (*(indexL [i]) >= 0)) {
399
400      int code = (*nll) -> code ();
401
402      std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
403      indcoe2vector (indexL [i], coeff [i], lcoeff);
404
405      /*std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
406      for (int i=0, *ind = indexL; *ind >= 0; *ind++, i++)
407      lcoeff.push_back (std::pair <exprVar *, CouNumber> (Var (*ind), coeff [i]));*/
408
409      if ((code == COU_EXPRSUM) || 
410          (code == COU_EXPRGROUP)) {
411
412        body    = exprGroup::genExprGroup (0., lcoeff, (*nll) -> ArgList (), (*nll) -> nArgs ());
413        // delete node without deleting children (they are now in body)
414        (*nll) -> ArgList (NULL);
415        delete *nll;
416        delete [] nll;
417      }
418      else body = exprGroup::genExprGroup (0., lcoeff, nll, 1);
419    }
420    else {
421      body = *nll;
422      delete [] nll;
423    }
424
425    expression *subst = body -> simplify ();
426    if (subst) {
427      delete body; // VALGRIND
428      body = subst;
429    }
430
431    // add them (and set lower-upper bound)
432    switch (sign) {
433
434    case COUENNE_EQ:  addEQConstraint  (body, new exprConst (ub)); break;
435    case COUENNE_LE:  addLEConstraint  (body, new exprConst (ub)); break;
436    case COUENNE_GE:  addGEConstraint  (body, new exprConst (lb)); break;
437    case COUENNE_RNG: addRNGConstraint (body, new exprConst (lb), 
438                                              new exprConst (ub)); break;
439    default: printf ("could not recognize constraint\n"); return -1;
440    }
441
442    delete [] indexL [i];
443    delete [] coeff  [i];
444  }
445
446  delete [] indexL;
447  delete [] coeff;
448
449  // create room for problem's variables and bounds
450  CouNumber
451    *= (CouNumber *) malloc ((n_var + ndefined_) * sizeof (CouNumber)),
452    *lb = (CouNumber *) malloc ((n_var + ndefined_) * sizeof (CouNumber)),
453    *ub = (CouNumber *) malloc ((n_var + ndefined_) * sizeof (CouNumber));
454
455  for (int i = n_var + ndefined_; i--;) {
456    x  [i] =  0.;
457    lb [i] = -COUENNE_INFINITY;
458    ub [i] =  COUENNE_INFINITY;
459  }
460
461  domain_.push (n_var + ndefined_, x, lb, ub);
462
463  free (x); free (lb); free (ub);
464
465  // lower and upper bounds ///////////////////////////////////////////////////////////////
466
467  if (LUv) {
468
469    real *Uvx_copy = Uvx;
470
471    if (!Uvx_copy)
472      for (register int i=0; i<n_var; i++) {
473
474        register int j = 2*i;
475
476        Lb (i) = LUv [j];
477        Ub (i) = LUv [j+1];
478      }
479    else
480      for (register int i=n_var; i--;) {
481        Lb (i) = LUv [i];
482        Ub (i) = Uvx_copy [i];
483      }
484
485  } else
486    for (register int i=n_var; i--;) {
487      Lb (i) = - COUENNE_INFINITY;
488      Ub (i) =   COUENNE_INFINITY;
489    }
490
491  // initial x ////////////////////////////////////////////////////////////////////
492
493  for (register int i=n_var; i--;) 
494
495    if (X0 && havex0 [i]) X (i) = X0 [i]; 
496
497    else {
498
499      CouNumber x, l = Lb (i), u = Ub (i);
500
501      if      (l < - COUENNE_INFINITY)
502        if    (u >   COUENNE_INFINITY)  x = 0.;
503        else                            x = u;
504      else if (u >   COUENNE_INFINITY)  x = l;
505      else                              x = 0.5 * (l+u);
506
507      X (i) = x;
508    }
509
510  for (register int i=n_var; i<ndefined_; i++) {
511
512    X  (i) =  0.;
513    Lb (i) = -COUENNE_INFINITY;
514    Ub (i) =  COUENNE_INFINITY;
515  }
516
517  delete [] nterms;
518
519  return 0;
520}
Note: See TracBrowser for help on using the repository browser.