source: stable/0.1/Couenne/src/readnl/readnl.cpp @ 81

Last change on this file since 81 was 81, checked in by pbelotti, 11 years ago

merged changes from Couenne-trunk:80 (fix casts for AMPL function pointers)

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