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

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