source: trunk/Couenne/src/standardize/standardize.cpp @ 306

Last change on this file since 306 was 306, checked in by pbelotti, 10 years ago

eliminate redundant (empty) constraints. Thanks to Leo Liberti and Giuliano Casale for pointing that out.

  • Property svn:keywords set to Id
File size: 12.7 KB
Line 
1/*
2 * $Id: standardize.cpp 306 2010-03-05 04:22:05Z pbelotti $
3 *
4 * Name:    standardize.cpp
5 * Author:  Pietro Belotti
6 * Purpose: standardize all expressions in a problem
7 *
8 * (C) Carnegie-Mellon University, 2006-09.
9 * This file is licensed under the Common Public License (CPL)
10 */
11
12#include <vector>
13
14#include "CoinHelperFunctions.hpp"
15
16#include "CouenneTypes.hpp"
17#include "expression.hpp"
18#include "exprIVar.hpp"
19#include "exprClone.hpp"
20#include "CouenneProblem.hpp"
21#include "CouenneProblemElem.hpp"
22#include "depGraph.hpp"
23
24//#define DEBUG
25
26/// standardize (nonlinear) common expressions, objectives, and constraints
27
28bool CouenneProblem::standardize () {
29
30#ifdef DEBUG
31  printf ("START: current point: %d vars -------------------\n", 
32          domain_.current () -> Dimension ());
33  for (int i=0; i<domain_.current () -> Dimension (); i++)
34    printf ("%3d %20g [%20g %20g]\n", i, domain_.x (i), domain_.lb (i), domain_.ub (i));
35#endif
36
37  bool retval = true;
38
39  // create dependence graph to assign an order to the evaluation (and
40  // bound propagation, and -- in reverse direction -- implied bounds)
41  graph_ = new DepGraph;
42
43  for (std::vector <exprVar *>::iterator i = variables_ . begin ();
44       i != variables_ . end (); ++i)
45    graph_ -> insert (*i);
46
47  // allocate space in auxiliaries_ from commonexprs_
48
49  int initVar = variables_ . size () - commonexprs_ . size ();
50
51  // Defined variables ///////////////////////////////////////////
52
53  // standardize initial aux variables (aka defined variables, aka
54  // common expression)
55
56#ifdef DEBUG
57  if (commonexprs_.size ()) printf ("%d common exprs, initVar = %d = %d - %d\n", 
58                                    commonexprs_.size (), 
59                                    initVar, 
60                                    variables_ . size (), 
61                                    commonexprs_ . size ());
62#endif
63
64  for (std::vector <expression *>::iterator i = commonexprs_ . begin ();
65       i != commonexprs_ . end (); ++i) {
66
67#ifdef DEBUG
68    printf ("-- stdz common expr [%d] :=", initVar); fflush (stdout);
69    (*i) -> print (); printf ("\n"); fflush (stdout);
70#endif
71
72    exprAux *naux = (*i) -> standardize (this, false);
73
74    expression *img = naux -> Image ();
75
76    // trick to obtain same index as common expression: create exprAux
77    // with index initVar and replace occurrences with address of
78    // newly created exprAux through auxiliarize()
79
80    exprAux *newvar = new exprAux (img, initVar, 1 + img -> rank (), exprAux::Unset, &domain_);
81    //img -> isInteger () ? exprAux::Integer : exprAux::Continuous);
82
83    auxiliarize (newvar); // takes care of putting newvar at right position in variables_
84
85    graph_ -> insert (newvar);
86    graph_ -> erase (naux);
87
88    //variables_ . erase (variables_ . end () - 1);
89
90#ifdef DEBUG
91    if (naux) {
92      printf ("done: "); fflush (stdout);
93      naux -> print ();
94      printf (" := "); fflush (stdout);
95      naux -> Image () -> print (); printf ("\n..."); fflush (stdout);
96    } else if (*i) {
97      (*i) -> print ();
98      //printf (" := "); fflush (stdout);
99      //aux -> Image () -> print ();
100      printf ("\n");
101    } else printf ("[n]aux NULL!\n");
102#endif
103
104    initVar++;
105  }
106
107  // OBJECTIVES //////////////////////////////////////////////////////////////////////////////
108
109  for (std::vector <CouenneObjective *>::iterator i = objectives_.begin ();
110       i != objectives_.end (); ++i) {
111
112#ifdef DEBUG
113    printf ("Objective [code: %d]", (*i) -> Body () -> code ()); (*i) -> print ();
114#endif
115
116    exprAux *aux = (*i) -> standardize (this);
117
118#ifdef DEBUG
119    printf ("      --> "); (*i) -> print (); 
120    if (aux) {printf (" -- aux is "); aux -> print (); printf ("\n");}
121#endif
122
123    if (aux) {
124      //delete ((*i) -> Body ());
125      (*i) -> Body (new exprClone (aux));
126    }
127
128#ifdef DEBUG
129    printf ("      --> "); (*i) -> print (); printf ("...................\n");
130#endif
131  }
132
133  // commuted_ is an array with a flag for each original variable,
134  // which is true at position i if initially original variable x_i
135  // went auxiliary
136
137  commuted_ = new bool [nVars ()];
138  for (int i = nVars (); i--;)
139    *commuted_++ = false;
140  commuted_ -= nVars ();
141
142  //  std::vector <CouenneConstraint *> con2;
143  std::vector <std::vector <CouenneConstraint *>::iterator> iters2erase;
144
145  // CONSTRAINTS /////////////////////////////////////////////////////////////////////////////
146
147  for (std::vector <CouenneConstraint *>::iterator i = constraints_.begin (); 
148       i != constraints_.end (); ++i) {
149
150#ifdef DEBUG
151    printf ("############# Constraint "); 
152    (*i) -> print ();
153#endif
154
155    CouNumber
156      conLb = (*((*i) -> Lb ())) (),
157      conUb = (*((*i) -> Ub ())) ();
158
159    // sanity check: if (due to bad model or redundancies) the
160    // constraint's body is constant, check if it's within bounds.
161
162    expression *eBody = (*i) -> Body ();
163
164    if (eBody -> Linearity () <= CONSTANT) {
165
166      CouNumber bodyVal = (*eBody)();
167      if ((bodyVal < conLb - COUENNE_EPS) ||
168          (bodyVal > conUb + COUENNE_EPS)) { // all variables eliminated, but out of bounds
169       
170        jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, 
171                          "Constraint %d: all variables eliminated, but value %g out of bounds [%g,%g]\n", 
172                          bodyVal, conLb, conUb);
173        retval = false;
174        break;
175
176      } else {
177        iters2erase.push_back (i);
178        continue; // all variables eliminated and constraint is redundant
179      }
180    }
181
182    exprAux *aux = (*i) -> standardize (this);
183
184#ifdef DEBUG
185    printf (" ==> [%d] ", aux ? (aux -> Index ()) : -1); 
186    (*i) -> print ();
187#endif
188
189    if (aux) { // save if standardized
190
191      // this is a top level auxiliary, i.e., an auxiliary whose node
192      // in the DAG stands at the maximum level -- no other variable
193      // depends on it as it is the lhs of a constraint.
194      aux -> top_level () = true;
195
196      //printf ("delete %x: ", ((*i) -> Body ())); ((*i) -> Body ()) -> print ();
197      //printf ("\n");
198      //delete ((*i) -> Body ());
199      (*i) -> Body (new exprClone (aux));
200      //      con2.push_back (*i);
201    }
202    else {
203      CouNumber lb, ub;
204      (*i) -> Body () -> getBounds (lb, ub);
205      if ((((*((*i) -> Lb ())) ()) > ub) ||
206          (((*((*i) -> Ub ())) ()) < lb)) {
207        jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, "found infeasible constraint [%g,%g]\n", lb, ub);
208#ifdef DEBUG
209        (*i) -> print ();
210#endif
211        retval = false;
212      }
213      iters2erase.push_back (i);
214    }
215
216    //(*i) -> Body () -> realign (this);
217
218#ifdef DEBUG
219    printf (" --> ");
220    (*i) -> print ();
221    printf ("..............................................................\n");
222    //print ();
223#endif
224
225    /*printf ("=== "); fflush (stdout);
226    aux -> print (); printf (" := "); fflush (stdout);
227    aux -> Image () -> print (); printf ("\n");*/
228  }
229
230  for (unsigned int i = iters2erase.size (); i--;)
231    constraints_. erase (iters2erase [i]);
232
233#ifdef DEBUG
234  // Use with caution. Bounds on auxs are not defined yet, so valgrind complains
235  printf ("done with standardization: (careful, bounds below can be nonsense)\n");
236  print (); 
237#endif
238
239  // Create evaluation order ////////////////////////////////////////////////////
240
241  delete auxSet_;
242
243  // reallocate space for enlarged set of variables
244  domain_.current () -> resize (nVars ());
245
246  //graph_ -> print ();
247  graph_ -> createOrder ();
248  //graph_ -> print ();
249
250  assert (graph_ -> checkCycles () == false);
251
252  // Fill numbering structure /////////////////////////////////////////////////
253
254  int n = nVars ();
255  numbering_ = new int [n];
256  std::set <DepNode *, compNode> vertices = graph_ -> Vertices ();
257
258  for (std::set <DepNode *, compNode>::iterator i = vertices.begin ();
259       i != vertices.end (); ++i)
260
261    numbering_ [(*i) -> Order ()] = (*i) -> Index (); 
262
263  //////////////////////////////////////////////////////////////////////////////
264
265  for (int i = 0; i < nVars (); i++) {
266
267    int ord = numbering_ [i];
268
269    if (variables_ [ord] -> Type () == AUX) {
270
271      // initial auxiliary bounds are infinite (they are later changed
272      // through branching)
273
274      if (variables_ [ord] -> Index () >= nOrigVars_) { // and one that was not an original, originally...
275
276        domain_.lb (ord) = -COIN_DBL_MAX;
277        domain_.ub (ord) =  COIN_DBL_MAX;
278      }
279
280      //printf ("from "); variables_ [ord] -> Lb    () -> print ();
281
282      // tighten them with propagated bounds
283      variables_ [ord] -> crossBounds ();
284
285      //printf ("to "); variables_ [ord] -> Lb    () -> print (); printf (", now eval\n");
286
287#ifdef DEBUG
288      printf (":::: %3d %10g [%10g, %10g]", 
289              ord, domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
290#endif
291
292      // and evaluate them
293      domain_.(ord) = (*(variables_ [ord] -> Image ())) ();
294      domain_.lb (ord) = (*(variables_ [ord] -> Lb    ())) ();
295      domain_.ub (ord) = (*(variables_ [ord] -> Ub    ())) ();
296
297#ifdef DEBUG
298      printf (" --> %10g [%10g, %10g] [", 
299              domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
300      variables_ [ord] -> Lb () -> print (); printf (",");
301      variables_ [ord] -> Ub () -> print (); printf ("]\n");
302#endif
303
304      bool integer = variables_ [ord] -> isInteger ();
305
306      if (integer) {
307        domain_.lb (ord) = ceil  (domain_.lb (ord) - COUENNE_EPS);
308        domain_.ub (ord) = floor (domain_.ub (ord) + COUENNE_EPS);
309      }
310    }
311  }
312
313  // TODO: resolve duplicate index in exprQuad before restoring this
314
315  std::string delete_redund;
316  if (bonBase_)
317    bonBase_ -> options () -> GetStringValue ("delete_redundant", delete_redund, "couenne."); 
318  else
319    delete_redund = "yes";
320
321  if (delete_redund == "yes")
322
323    // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
324    for (std::vector <exprVar *>::iterator i = variables_.begin (); 
325       i != variables_.end (); ++i)
326
327    if ((*i) -> Type () == AUX) {
328
329      int type = (*i) -> Image () -> Type ();
330
331      if ((type == VAR) || (type == AUX)) {
332
333        // found w_k = x_h.
334        //
335        // Check if either is integer, the survivor will be integer too
336        // Replace all occurrences of w_k with x_h
337
338        /*printf ("redundancy: ");
339        (*i)             -> print (); printf (" := ");
340        (*i) -> Image () -> print (); printf ("\n");*/
341
342        // use the info on the variable to be discarded: tighter
343        // bounds and integrality that the replacement might not have.
344
345        int 
346          indStays  = (*i) -> Image () -> Index (), // index h
347          indLeaves = (*i)             -> Index (); // index k
348
349        if (indStays == indLeaves)  // very strange case, w_h = x_h
350          continue;
351
352        // do not swap! x_h could be in turn an auxiliary...
353        //
354        //if (indStays > indLeaves)
355        //{int swap = indStays; indStays = indLeaves; indLeaves = swap;} // swap
356
357        exprVar
358          *varStays  = variables_ [indStays],
359          *varLeaves = variables_ [indLeaves];
360
361        // intersect features of the two variables (integrality, bounds)
362
363        varStays -> lb () = varLeaves -> lb () = CoinMax (varStays -> lb (), varLeaves -> lb ());
364        varStays -> ub () = varLeaves -> ub () = CoinMin (varStays -> ub (), varLeaves -> ub ());
365
366        if (varStays  -> isInteger () ||
367            varLeaves -> isInteger ()) {
368
369          varStays -> lb () = ceil  (varStays -> lb ());
370          varStays -> ub () = floor (varStays -> ub ());
371
372          if (varStays -> Type () == AUX)
373            varStays -> setInteger (true);
374          else {
375            //expression *old = varStays; // !!! leak
376            variables_ [indStays] = varStays = new exprIVar (indStays, &domain_);
377            auxiliarize (varStays); // replace it everywhere in the problem
378            //delete old;
379          }
380        }
381
382        auxiliarize (varLeaves, varStays); // now replace occurrences of w_k with x_h
383
384        //if (varLeaves -> Index () >= nOrigVars_) // why check? It's not there anymore.
385        varLeaves -> zeroMult (); // disable this variable
386      }
387    }
388
389  /// re-check integrality. This is necessary as the initial
390  /// integrality check is done on some continuous variables, which
391  /// may turn out to be identical to other, integer, variables. See
392  /// minlplib/ex1223.nl, where x_29 = x_4^2 and x_4=x_9, with x_4
393  /// declared continuous and x_9 integer
394  for (int ii=0; ii < nVars (); ii++) {
395
396    int i = numbering_ [ii];
397
398    if ((Var (i) -> Multiplicity () > 0) &&
399        (Var (i) -> Type () == AUX) &&
400        (Var (i) -> Image () -> isInteger ()))
401      Var (i) -> setInteger (true);
402
403    //if (Var (i) -> Multiplicity () == 0)
404    //Lb (i) = Ub (i) = 0.;
405  }
406
407  // check how many multiplications there are
408
409//   int nmul = 0;
410//   // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
411//   for (std::vector <exprVar *>::iterator i = variables_.begin ();
412//        i != variables_.end (); ++i) {
413
414//     if ((*i) -> Type () != AUX ||
415//      (*i) -> Multiplicity () <= 0)
416//       continue;
417
418//     expression *img = (*i) -> Image ();
419
420//     if (img -> code () != COU_EXPRMUL)
421//       continue;
422
423//     expression **args = img -> ArgList ();
424
425//     if ((args [0] -> Type () == AUX ||
426//       args [0] -> Type () == VAR) &&
427//      (args [1] -> Type () == AUX ||
428//       args [1] -> Type () == VAR))
429//       nmul++;
430//   }
431
432//   printf ("MULS: %d/%d\n", nmul, variables_.size());
433
434//   exit (-1);
435
436  // TODO: re-compute ranks
437
438  delete [] commuted_;  commuted_ = NULL;
439  delete    graph_;     graph_    = NULL;
440
441  return retval;
442}
Note: See TracBrowser for help on using the repository browser.