source: trunk/Couenne/src/readnl/readnl.cpp @ 79

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

name fixing. Trying to fix a crippling cast conversion for function pointer comparison (dangerous...)

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