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

Last change on this file since 112 was 112, checked in by pbelotti, 13 years ago

isolate restoreUnused to when there is at least one unused -- takes care of reformulated stockcycle problem (thanks to C. D'Ambrosio for the bug report). Check linear terms before adding a full exprGroup -- need a constructor pilot. Quicker solution feasibility check (with initial integrality check).

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