source: trunk/Couenne/src/standardize/splitAux.cpp @ 91

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

fixed serious bug. Failed to find optimum due to duplicate variable...

File size: 13.0 KB
Line 
1/*
2 * Name:    splitAux.cpp
3 * Author:  Pietro Belotti
4 * Purpose: extract auxiliary variable from implicit equality constraint
5 *
6 * (C) Carnegie-Mellon University, 2007.
7 * This file is licensed under the Common Public License (CPL)
8 */
9
10#include "CouenneProblem.hpp"
11
12#include "CoinHelperFunctions.hpp"
13
14#include "exprClone.hpp"
15#include "exprSum.hpp"
16#include "exprMul.hpp"
17#include "exprGroup.hpp"
18#include "exprQuad.hpp"
19#include "lqelems.hpp"
20
21//#define DEBUG
22
23/// given an element of a sum, check if it is a variable (possibly
24/// with a coefficient) and return its index (and the coefficient) if
25/// it has not been spotted as an auxiliary
26void elementBreak (expression *, int &, CouNumber &);
27
28
29/// split a constraint w - f(x) = c into w's index (returned)
30/// and rest = f(x) + c
31
32int CouenneProblem::splitAux (CouNumber rhs, expression *body, expression *&rest, bool *wentAux) {
33
34  int auxInd = -1,          // index of the auxiliary to be extracted
35    code = body -> code (); // type of expression
36
37  expression **alist = body -> ArgList ();
38
39#ifdef DEBUG
40  printf ("|||||||||| Splitting "); body -> print (); printf ("\n");
41#endif
42
43  switch (code) { // constraint h(x) = 0 may be in different forms:
44                  // subtraction, sum, linear group
45
46    /////////////////////////////////////////////////////////////////////////////
47
48  case COU_EXPRSUB: {
49
50    // simplest case, we have w-f(x)=rhs or f(x)-w=rhs or f(x)-g(x)=rhs
51
52    int pos = 0;
53    CouNumber coeff = 1;
54
55    auxInd = (*alist) -> Index ();
56
57    if (auxInd < 0)
58      elementBreak (*alist, auxInd, coeff); // check first element
59
60    if ((auxInd < 0) || 
61        wentAux [auxInd] || 
62        (alist [1] -> dependsOn (auxInd, TAG_AND_RECURSIVE) >= 1)) {
63
64      auxInd = (alist [1]) -> Index ();
65
66      if (auxInd < 0)
67        elementBreak (alist [1], auxInd, coeff); // check second element
68      else coeff = 1;
69
70      if ((auxInd < 0) ||       // no index found
71          wentAux [auxInd] ||   // or, found but invalid
72          (alist [0] -> dependsOn (auxInd, TAG_AND_RECURSIVE) >= 1)) 
73                                // or, variable depends upon itself
74        return -1;
75
76      pos = 1;
77    }
78
79    //////////////////////
80
81    // what remains is the "independent" expression
82    expression *clone = alist [1 - pos] -> clone (&domain_); 
83
84    expression *auxdef = 
85      (fabs (coeff - 1) < COUENNE_EPS) ?          // if coefficient is 1
86      (clone) :                                   //     do nothing
87      new exprMul (new exprConst (coeff), clone); //     else multiply it by coefficient
88
89    rest = (fabs (rhs) < COUENNE_EPS) ?                              // no extra constant?
90      (auxdef) :                                                     // just put other argument
91      (new exprSum (auxdef, new exprConst ((pos==1) ? -rhs : rhs))); // otherwise sum it with \pm rhs
92
93  } break;
94
95  ////////////////////////////////////////////////////////////////////////////
96
97  case COU_EXPRQUAD:
98  case COU_EXPRGROUP:
99  case COU_EXPRSUM: {
100
101    // an expr{Sum,Group,Quad}. Decompose the expression and
102    // re-assemble (to look for right variable)
103
104    // data structure to be used below if there is a linear term.
105    // which specifies position within arrays (negative for linear
106    // part of exprGroup, positive for all elements of exprSum)
107
108    int maxindex = -1, nlin = 0, which = 1;
109    CouNumber c0 = 0., auxcoe = 1;
110    bool which_was_set = false;
111
112    // check indices of linear part /////////////////////////////////
113
114    if (code != COU_EXPRSUM) {
115
116      exprGroup *egBody = dynamic_cast <exprGroup *> (body);
117      exprGroup::lincoeff &lcoe = egBody -> lcoeff ();
118
119      // import exprGroup linear structure
120
121      c0 = egBody -> getc0 ();
122
123      for (int i=0, n = lcoe.size (); n--; i++, nlin++) {
124
125        int j = lcoe [i]. first -> Index ();
126
127        //lincoe [i] = lcoe [i]. second;
128
129        if ((j > maxindex) && 
130            !(wentAux [j]) && 
131            (fabs (lcoe [i]. second) > COUENNE_EPS)) {
132
133          // fake cut in linind and check dependence. Only mark if
134          // dependsOn() gives 0
135
136          exprVar *saveVar = lcoe [i].first;
137          //lcoe [i].first = new exprVar (nVars ()); // better use index -1
138          lcoe [i].first = new exprVar (-1);
139
140          if (body -> dependsOn (j, TAG_AND_RECURSIVE) == 0) {
141
142#ifdef DEBUG
143            printf ("body does not depend on x%d: ", j);
144            body -> print ();
145            printf ("\n");
146#endif
147            // mark which with negative number
148            which    = - nlin - 1;
149            auxcoe   = lcoe [i]. second;
150            maxindex = j;
151          }
152
153          delete lcoe [i].first;
154          lcoe [i].first = saveVar;
155        }
156      }
157    }
158
159    if (which != 1) 
160      which_was_set = true;
161    else which = -1;
162
163    // check indices of elements of (nonlinear) sum /////////////////////////////////
164
165    for (int i = body -> nArgs (); i--;) {
166
167      CouNumber coeff = 1;
168      int index = alist [i] -> Index ();
169
170      if (index < 0)
171        elementBreak (alist [i], index, coeff);
172
173      if ((index > maxindex) &&
174          !(wentAux [index]) &&
175          (fabs (coeff) > COUENNE_EPS)) {
176
177        // fake a cut in the arglist and check
178
179        expression *cut = alist [i];
180        alist [i] = new exprConst (0.);
181
182        // not enough... check now linear (and quadratic!) terms
183
184        if (body -> dependsOn (index, TAG_AND_RECURSIVE) == 0) {
185
186          maxindex = index;
187          which    = i;
188          auxcoe   = coeff;
189        }
190
191        delete alist [i];
192        alist [i] = cut;
193      }
194    }
195
196    if (!which_was_set && (which == -1)) // which has not been set
197      return -1;
198
199    ///////////////////////////////////////////////////////////////////
200
201    if (maxindex < 0) break; // no substitute found ==> no hidden auxiliary
202
203    // create a new exprGroup, exprQuad, or exprSum with all elements but the
204    // extracted auxiliary
205
206    // start with exprSum
207    int nargs = body -> nArgs ();
208    expression **newarglist;
209
210#ifdef DEBUG
211    printf (" [[ind %d, coe %.1g, wh %d, nargs %d, nlin %d]] ", 
212            maxindex, auxcoe, which, nargs, nlin);
213#endif
214
215    if (nargs > 0) { // there is an element in the nonlinear sum to be drawn
216
217      int j, mid = (which < 0) ? nargs : which;
218
219      newarglist = new expression * [nargs + 1];
220
221      for (j=0; j<mid;   j++) newarglist [j]   = alist [j] -> clone (&domain_);
222      for (j++; j<nargs; j++) newarglist [j-1] = alist [j] -> clone (&domain_);
223
224      // nl arglist is done, later decide whether to incorporate it as
225      // it is or with a coefficient
226
227    } else { // no nonlinear arguments, or the only one is the new aux
228
229      nargs++; // !!!!!!!!!!!!!!!!!!!!!!!!!
230      newarglist  = new expression *;
231      *newarglist = new exprConst (0.);
232    }
233
234    // form rhs linear part ////////////////////////////////////////////////////
235
236    int       *linind2 = NULL;
237    CouNumber *lincoe2 = NULL;
238
239    // in case this was (and will be) an exprQuad
240    int *qindI = NULL, 
241        *qindJ = NULL;
242    CouNumber *qcoe = NULL;
243
244    if (nlin > 0) { // there is an element in the linear sum to be drawn
245
246      exprGroup *egBody = dynamic_cast <exprGroup *> (body);
247      exprGroup::lincoeff &lcoe = egBody -> lcoeff ();
248
249      int mid = (which >= 0) ? nlin : - which - 1;
250
251      linind2 = new int       [nlin + 1];
252      lincoe2 = new CouNumber [nlin + 1];
253
254      register int j;
255
256#ifdef DEBUG
257      //for (j=0; j<mid;  j++) printf ("{%g x%d} ", lincoe [j], linind [j]);
258      //for (j++; j<nlin; j++) printf ("{%g x%d} ", lincoe [j], linind [j]);
259#endif
260
261      CouNumber divider = -1. / auxcoe;
262
263      for (j=0; j<mid;  j++){linind2[j]  =lcoe[j].first->Index();lincoe2[j]  =divider*lcoe[j].second;}
264      for (j++; j<nlin; j++){linind2[j-1]=lcoe[j].first->Index();lincoe2[j-1]=divider*lcoe[j].second;}
265
266      linind2 [j-1] = -1; // terminate list of indices
267
268#ifdef DEBUG
269      for (j=0; j<mid;  j++) printf ("<%g x%d> ", lincoe2 [j],   linind2 [j]);
270      for (j++; j<nlin; j++) printf ("<%g x%d> ", lincoe2 [j-1], linind2 [j-1]);
271#endif
272
273      if (code == COU_EXPRQUAD) { // copy quadratic elements
274
275        exprQuad *eq = dynamic_cast <exprQuad *> (body);
276
277        int nqt = eq -> getnQTerms (), j=0;
278
279        qindI = new int [nqt];
280        qindJ = new int [nqt];
281        qcoe  = new CouNumber [nqt];
282
283        exprQuad::sparseQ &M = eq -> getQ ();
284
285        for (exprQuad::sparseQ::iterator row = M.begin (); 
286             row != M.end (); ++row) {
287
288          int xind = row -> first -> Index ();
289
290          for (exprQuad::sparseQcol::iterator col = row -> second.begin (); 
291               col != row -> second.end (); ++col, ++j) {
292
293            qindI [j] = xind;
294            qindJ [j] = col -> first -> Index ();
295            qcoe  [j] = divider * col -> second;
296          }
297        }
298      }
299
300      // nl arglist is done, later decide whether to incorporate it as
301      // it is or with a coefficient
302    }
303
304    // the extracted index is one term of...
305    if (which >= 0) --nargs; // ...the nonlinear sum
306    else            --nlin;  // ...the linear part
307
308#ifdef DEBUG
309    printf ("\n::: auxcoe %g, rhs %g, lin %d, nl %d\n", auxcoe, rhs, nlin, nargs);
310#endif
311
312    // all is ready to take the independent stuff to the other side of
313    // the inequality.
314
315    if ((code == COU_EXPRQUAD) || 
316        (code == COU_EXPRGROUP) && (nlin > 0)) { 
317
318      // an exprGroup with at least one linear term left
319      //
320      // build new vectors for index and coeff. Two cases:
321      //
322      // 1)  f(x) + c0 -  w = rhs   =====>   w =       f(x) + c0 - rhs
323      // 2)  f(x) + c0 + aw = rhs   =====>   w = -1/a (f(x) + c0 - rhs), a != -1
324
325      if (fabs (auxcoe + 1) < COUENNE_EPS) {
326
327        //std::vector <std::pair <exprVar *, CouNumber> >
328        exprGroup::lincoeff lcoeff;
329        indcoe2vector (linind2, lincoe2, lcoeff);
330
331        if (code == COU_EXPRGROUP)
332          rest = new exprGroup (c0-rhs, lcoeff, newarglist, nargs);
333
334        else {
335          std::vector <quadElem> qcoeff;
336          indcoe2vector (qindI, qindJ, qcoe, qcoeff);
337          rest = new exprQuad  (c0-rhs, lcoeff, qcoeff, newarglist, nargs);
338        }
339      }
340      else {
341
342        expression **mullist = new expression * [1];
343
344        // only nl term is constant
345        if ((nargs <= 1) && ((*newarglist) -> Linearity () <= CONSTANT)) {
346          *mullist = new exprConst ((*newarglist) -> Value ());
347          //delete *newarglist;
348          delete [] newarglist;
349        }
350        else if ((nargs <= 1) && 
351                 ((*newarglist) -> code () == COU_EXPROPP) &&
352                 (fabs (auxcoe - 1) < COUENNE_EPS))
353          //*mullist = new exprSum (newarglist, nargs);
354          *mullist = (*newarglist) -> Argument ();//, nargs);
355        else { // multiple nonlinear terms, multiply them by -1/auxcoe
356
357          if (nargs <= 1)
358            *mullist = new exprMul (new exprConst (-1. / auxcoe), 
359                                    *newarglist); // !!! one more leak?
360          else 
361            *mullist = new exprMul (new exprConst (-1. / auxcoe), 
362                                    new exprSum (newarglist, nargs));
363        }
364
365        std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
366        indcoe2vector (linind2, lincoe2, lcoeff);
367
368        // final outcome: -1/a (f(x) + c0 - rhs)
369        if (code == COU_EXPRGROUP)
370          rest = new exprGroup ((rhs - c0) / auxcoe, lcoeff,        mullist,1);
371
372        else {
373          std::vector <quadElem> qcoeff;
374          indcoe2vector (qindI, qindJ, qcoe, qcoeff);
375          rest = new exprQuad  ((rhs - c0) / auxcoe, lcoeff, qcoeff, mullist,1);
376        }
377      }
378    }
379    else { // simple exprSum
380
381      if (fabs (rhs) > COUENNE_EPS) { // have to add constant to exprSum
382
383        if ((nargs == 1) && ((*newarglist) -> Type () == CONST)) {
384
385          CouNumber val = (*newarglist) -> Value () - rhs;
386          delete *newarglist;
387          *newarglist = new exprConst (val);
388        } 
389        else newarglist [nargs++] = new exprConst (-rhs);
390      }
391
392      // now exprSum is complete with -rhs. Send it to right hand side
393
394      expression *auxDef;
395
396      if (nargs==1) {
397        auxDef = *newarglist;
398        delete [] newarglist;
399      } else auxDef = new exprSum (newarglist, nargs);
400
401      if (fabs (auxcoe + 1) < COUENNE_EPS)
402        rest = auxDef;
403      else if ((fabs (auxcoe - 1) < COUENNE_EPS) && 
404               (auxDef -> code () == COU_EXPROPP))
405        rest = auxDef -> Argument ();
406      else // TODO: check if auxdef is an exprOpp or an exprMul
407           // (k*f(x)) and -1/auxcoe simplifies
408        rest = new exprMul (new exprConst (-1./auxcoe), 
409                            ((auxDef -> Type () == AUX) ||
410                             (auxDef -> Type () == VAR)) ? new exprClone (auxDef) : auxDef);
411    }
412
413    if (linind2) {
414      delete [] linind2;
415      delete [] lincoe2;
416    }
417
418#ifdef DEBUG
419    printf ("gotten "); rest -> print (); printf ("\n");
420#endif
421
422    auxInd = maxindex;
423
424  } break;
425
426  default: break;
427  } // end switch () ////////////////////////////////////////////////////////
428
429  if ((auxInd < 0) || (wentAux [auxInd]))
430    return -1;
431
432  // we have a variable, meaning the constraint body is of the form
433  // w +/- f(x) or f(x) +/- w
434
435  // standardize remaining of the expression
436
437#ifdef DEBUG
438  printf ("standardize rest (2nd level) "); fflush (stdout);
439  rest -> print (); printf (" [body = ");
440  body -> print (); printf ("]");
441#endif
442
443  // second argument is false to tell standardize not to create a new
444  // auxiliary variable (we are creating it ourselves, it's aux)
445  exprAux *aux = rest -> standardize (this, false);
446
447#ifdef DEBUG
448  printf (" {body = "); body -> print (); printf ("} ");
449  if (aux) {
450    printf (" --- aux = "); 
451    aux -> print ();
452    printf (" := ");
453    aux -> Image () -> print ();
454  }
455#endif
456
457  if (aux) {
458    //delete rest;
459    //rest = aux -> Image () -> clone (&domain_);
460    rest = aux -> Image ();// -> clone (&domain_);
461    aux -> Image (NULL);
462    delete aux;
463  }
464
465#ifdef DEBUG
466  printf (" ==> "); rest -> print ();
467  printf (" and "); fflush (stdout);
468  rest -> print (); printf ("\n");
469#endif
470
471  return auxInd;
472}
Note: See TracBrowser for help on using the repository browser.