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

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

fixed reformulation bug affecting st_glmp_fp3

  • Property svn:keywords set to Id
File size: 11.9 KB
Line 
1/*
2 * $Id: standardize.cpp 294 2010-02-13 14:04:23Z 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    exprAux *aux = (*i) -> standardize (this);
156
157#ifdef DEBUG
158    printf (" ==> [%d] ", aux ? (aux -> Index ()) : -1); 
159    (*i) -> print ();
160#endif
161
162    if (aux) { // save if standardized
163
164      // this is a top level auxiliary, i.e., an auxiliary whose node
165      // in the DAG stands at the maximum level -- no other variable
166      // depends on it as it is the lhs of a constraint.
167      aux -> top_level () = true;
168
169      //printf ("delete %x: ", ((*i) -> Body ())); ((*i) -> Body ()) -> print ();
170      //printf ("\n");
171      //delete ((*i) -> Body ());
172      (*i) -> Body (new exprClone (aux));
173      //      con2.push_back (*i);
174    }
175    else {
176      CouNumber lb, ub;
177      (*i) -> Body () -> getBounds (lb, ub);
178      if ((((*((*i) -> Lb ())) ()) > ub) ||
179          (((*((*i) -> Ub ())) ()) < lb)) {
180        jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, "found infeasible constraint [%g,%g]\n", lb, ub);
181#ifdef DEBUG
182        (*i) -> print ();
183#endif
184        retval = false;
185      }
186      iters2erase.push_back (i);
187    }
188
189    //(*i) -> Body () -> realign (this);
190
191#ifdef DEBUG
192    printf (" --> ");
193    (*i) -> print ();
194    printf ("..............................................................\n");
195    //print ();
196#endif
197
198    /*printf ("=== "); fflush (stdout);
199    aux -> print (); printf (" := "); fflush (stdout);
200    aux -> Image () -> print (); printf ("\n");*/
201  }
202
203  for (unsigned int i = iters2erase.size (); i--;)
204    constraints_. erase (iters2erase [i]);
205
206#ifdef DEBUG
207  // Use with caution. Bounds on auxs are not defined yet, so valgrind complains
208  printf ("done with standardization: (careful, bounds below can be nonsense)\n");
209  print (); 
210#endif
211
212  // Create evaluation order ////////////////////////////////////////////////////
213
214  delete auxSet_;
215
216  // reallocate space for enlarged set of variables
217  domain_.current () -> resize (nVars ());
218
219  //graph_ -> print ();
220  graph_ -> createOrder ();
221  //graph_ -> print ();
222
223  assert (graph_ -> checkCycles () == false);
224
225  // Fill numbering structure /////////////////////////////////////////////////
226
227  int n = nVars ();
228  numbering_ = new int [n];
229  std::set <DepNode *, compNode> vertices = graph_ -> Vertices ();
230
231  for (std::set <DepNode *, compNode>::iterator i = vertices.begin ();
232       i != vertices.end (); ++i)
233
234    numbering_ [(*i) -> Order ()] = (*i) -> Index (); 
235
236  //////////////////////////////////////////////////////////////////////////////
237
238  for (int i = 0; i < nVars (); i++) {
239
240    int ord = numbering_ [i];
241
242    if (variables_ [ord] -> Type () == AUX) {
243
244      // initial auxiliary bounds are infinite (they are later changed
245      // through branching)
246
247      if (variables_ [ord] -> Index () >= nOrigVars_) { // and one that was not an original, originally...
248
249        domain_.lb (ord) = -COIN_DBL_MAX;
250        domain_.ub (ord) =  COIN_DBL_MAX;
251      }
252
253      //printf ("from "); variables_ [ord] -> Lb    () -> print ();
254
255      // tighten them with propagated bounds
256      variables_ [ord] -> crossBounds ();
257
258      //printf ("to "); variables_ [ord] -> Lb    () -> print (); printf (", now eval\n");
259
260#ifdef DEBUG
261      printf (":::: %3d %10g [%10g, %10g]", 
262              ord, domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
263#endif
264
265      // and evaluate them
266      domain_.(ord) = (*(variables_ [ord] -> Image ())) ();
267      domain_.lb (ord) = (*(variables_ [ord] -> Lb    ())) ();
268      domain_.ub (ord) = (*(variables_ [ord] -> Ub    ())) ();
269
270#ifdef DEBUG
271      printf (" --> %10g [%10g, %10g] [", 
272              domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
273      variables_ [ord] -> Lb () -> print (); printf (",");
274      variables_ [ord] -> Ub () -> print (); printf ("]\n");
275#endif
276
277      bool integer = variables_ [ord] -> isInteger ();
278
279      if (integer) {
280        domain_.lb (ord) = ceil  (domain_.lb (ord) - COUENNE_EPS);
281        domain_.ub (ord) = floor (domain_.ub (ord) + COUENNE_EPS);
282      }
283    }
284  }
285
286  // TODO: resolve duplicate index in exprQuad before restoring this
287
288  std::string delete_redund;
289  if (bonBase_)
290    bonBase_ -> options () -> GetStringValue ("delete_redundant", delete_redund, "couenne."); 
291  else
292    delete_redund = "yes";
293
294  if (delete_redund == "yes")
295
296    // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
297    for (std::vector <exprVar *>::iterator i = variables_.begin (); 
298       i != variables_.end (); ++i)
299
300    if ((*i) -> Type () == AUX) {
301
302      int type = (*i) -> Image () -> Type ();
303
304      if ((type == VAR) || (type == AUX)) {
305
306        // found w_k = x_h.
307        //
308        // Check if either is integer, the survivor will be integer too
309        // Replace all occurrences of w_k with x_h
310
311        /*printf ("redundancy: ");
312        (*i)             -> print (); printf (" := ");
313        (*i) -> Image () -> print (); printf ("\n");*/
314
315        // use the info on the variable to be discarded: tighter
316        // bounds and integrality that the replacement might not have.
317
318        int 
319          indStays  = (*i) -> Image () -> Index (), // index h
320          indLeaves = (*i)             -> Index (); // index k
321
322        if (indStays == indLeaves)  // very strange case, w_h = x_h
323          continue;
324
325        // do not swap! x_h could be in turn an auxiliary...
326        //
327        //if (indStays > indLeaves)
328        //{int swap = indStays; indStays = indLeaves; indLeaves = swap;} // swap
329
330        exprVar
331          *varStays  = variables_ [indStays],
332          *varLeaves = variables_ [indLeaves];
333
334        // intersect features of the two variables (integrality, bounds)
335
336        varStays -> lb () = varLeaves -> lb () = CoinMax (varStays -> lb (), varLeaves -> lb ());
337        varStays -> ub () = varLeaves -> ub () = CoinMin (varStays -> ub (), varLeaves -> ub ());
338
339        if (varStays  -> isInteger () ||
340            varLeaves -> isInteger ()) {
341
342          varStays -> lb () = ceil  (varStays -> lb ());
343          varStays -> ub () = floor (varStays -> ub ());
344
345          if (varStays -> Type () == AUX)
346            varStays -> setInteger (true);
347          else {
348            //expression *old = varStays; // !!! leak
349            variables_ [indStays] = varStays = new exprIVar (indStays, &domain_);
350            auxiliarize (varStays); // replace it everywhere in the problem
351            //delete old;
352          }
353        }
354
355        auxiliarize (varLeaves, varStays); // now replace occurrences of w_k with x_h
356
357        //if (varLeaves -> Index () >= nOrigVars_) // why check? It's not there anymore.
358        varLeaves -> zeroMult (); // disable this variable
359      }
360    }
361
362  /// re-check integrality. This is necessary as the initial
363  /// integrality check is done on some continuous variables, which
364  /// may turn out to be identical to other, integer, variables. See
365  /// minlplib/ex1223.nl, where x_29 = x_4^2 and x_4=x_9, with x_4
366  /// declared continuous and x_9 integer
367  for (int ii=0; ii < nVars (); ii++) {
368
369    int i = numbering_ [ii];
370
371    if ((Var (i) -> Multiplicity () > 0) &&
372        (Var (i) -> Type () == AUX) &&
373        (Var (i) -> Image () -> isInteger ()))
374      Var (i) -> setInteger (true);
375
376    //if (Var (i) -> Multiplicity () == 0)
377    //Lb (i) = Ub (i) = 0.;
378  }
379
380  // check how many multiplications there are
381
382//   int nmul = 0;
383//   // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
384//   for (std::vector <exprVar *>::iterator i = variables_.begin ();
385//        i != variables_.end (); ++i) {
386
387//     if ((*i) -> Type () != AUX ||
388//      (*i) -> Multiplicity () <= 0)
389//       continue;
390
391//     expression *img = (*i) -> Image ();
392
393//     if (img -> code () != COU_EXPRMUL)
394//       continue;
395
396//     expression **args = img -> ArgList ();
397
398//     if ((args [0] -> Type () == AUX ||
399//       args [0] -> Type () == VAR) &&
400//      (args [1] -> Type () == AUX ||
401//       args [1] -> Type () == VAR))
402//       nmul++;
403//   }
404
405//   printf ("MULS: %d/%d\n", nmul, variables_.size());
406
407//   exit (-1);
408
409  // TODO: re-compute ranks
410
411  delete [] commuted_;  commuted_ = NULL;
412  delete    graph_;     graph_    = NULL;
413
414  return retval;
415}
Note: See TracBrowser for help on using the repository browser.