source: trunk/Couenne/src/problem/getIntegerCandidate.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.5 KB
Line 
1/*
2 * Name:    getIntegerCandidate.cpp
3 * Author:  Pietro Belotti
4 * Purpose: generate integer NLP point Y starting from fractional
5 *          solution using bound tightening
6 *
7 * (C) Carnegie-Mellon University, 2007-08.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include "CoinTime.hpp"
12#include "CouenneProblem.hpp"
13
14// lose patience after this many iterations of non-improving valid
15// tightening (step 1)
16#define VALID_ONLY_THRESHOLD 5
17
18/// GRASP for finding integer feasible solutions
19///
20/// Generate integer NLP point xInt starting from fractional solution
21/// xFrac, using (feasibility-based, i.e. cheap) bound tightening
22///
23/// return -1 if the problem is infeasible
24
25int CouenneProblem::getIntegerCandidate (const double *xFrac, double *xInt, 
26                                         double *lb, double *ub) const {
27  fillIntegerRank ();
28
29  if (numberInRank_.size () == 0) // there is no original integer to fix
30    return 0;
31
32  CouNumber *store_optimum = optimum_; // temporary place for debug optimal solution
33  optimum_ = NULL;
34
35  int
36    ncols    = nVars (), 
37    retval   = 0;
38
39  double
40    *olb   = new double [ncols],       *oub   = new double [ncols],      // outer bounds
41    *dualL = new double [nOrigVars_],  *dualR = new double [nOrigVars_]; // lb[objind] per fix index/direction
42
43  // copy in current bounds
44  CoinCopyN (Lb (), ncols, olb);
45  CoinCopyN (Ub (), ncols, oub);
46
47  // for now save fractional point into integer point
48  CoinCopyN (xFrac, nOrigVars_, xInt);
49
50  domain_.push (nVars (), xInt, lb, ub);
51
52  CoinCopyN (lb, nVars (), Lb ());
53  CoinCopyN (ub, nVars (), Ub ());
54
55  enum fixType *fixed = new enum fixType [nOrigVars_]; // integer variables that were fixed
56
57  for (int i=0; i<nOrigVars_; i++) 
58    fixed [i] = (Var (i) -> isInteger () &&        // work on integer variables only
59                 Var (i) -> Multiplicity () > 0) ? // don't care if unused variable
60      UNFIXED : CONTINUOUS;
61
62  // A more sophisticated rounding procedure
63  //
64  // Call aggressive tightening for all integer variables, setting
65  // each first at floor(x_i) and then at ceil(x_i).
66  //
67  //
68  // for each value r of the rank: 1 .. maxrank {
69  //
70  //   restrict search of unfixed original variables with rank r to
71  //   [floor,ceil] box around LP value
72  //
73  //   do {
74  //
75  //     threshold = nVarWithThisRank / k
76  //     niter = 0;
77  //     for each such (unfixed) variable {
78  //
79  //       try left
80  //       try right
81  //
82  //       if      both infeasible, abort
83  //       else if one infeasible, fix
84  //       else if both feasible
85  //         if niter < threshold skip
86  //         else                 fix based on lb(left<>right)
87  //     }
88  //
89  //     if no variables fixed, fix first unfixed
90  //     ++niter
91  //
92  //   } while there are unfixed variables at this rank
93  // }   
94  //
95  // check value of objective -> set cutoff and run bound tightening
96  //
97
98  const int infeasible = 1;
99
100  try {
101
102    // for all values of the integer rank
103
104    int rank = 1;
105
106    if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
107      printf ("=       ===========================================\n");
108      printf ("= BEGIN ===========================================\n");
109      printf ("=       ===========================================\n");
110      for (int i=0; i<nOrigVars_; i++)
111        if (variables_ [i] -> Multiplicity () > 0)
112          printf ("#### %4d: %d %c %2d frac %20g  [%20g,%20g]\n", 
113                  i, fixed [i], 
114                  variables_ [i] -> isInteger () ? 'I' : ' ',
115                  integerRank_ ? integerRank_ [i] : -1,
116                  xFrac [i], Lb (i), Ub (i));
117      printf ("---\n");
118      for (int i=nOrigVars_; i<nVars (); i++)
119        if (variables_ [i] -> Multiplicity () > 0)
120          printf ("#### %4d:   %c    frac %20g   [%20g,%20g]\n", 
121                  i, variables_ [i] -> isInteger () ? 'I' : ' ',
122                  //(integerRank_ && integerRank_ [i]) ? 'F' : ' ',
123                  X (i), Lb (i), Ub (i));
124      printf ("===================================================\n");
125      printf ("===================================================\n");
126      printf ("===================================================\n");
127    }
128
129    for (std::vector <int>::iterator rNum = numberInRank_.begin(); 
130         ++rNum != numberInRank_.end(); rank++) 
131
132      if (*rNum > 0) {
133
134        if (CoinCpuTime () > maxCpuTime_)
135          break;
136
137        // *rNum is the number of variable with integer rank equal to rank
138
139        // start restricting around current integer box
140        for (int i=0; i<nOrigVars_; i++) 
141          if ((Var (i) -> Multiplicity () > 0) && // alive variable
142              (Var (i) -> isInteger    ())     && // integer, may fix if independent of other integers
143              (integerRank_ [i] == rank)) {
144
145            Lb (i) = CoinMax (Lb (i), floor (xFrac [i])); 
146            Ub (i) = CoinMin (Ub (i), ceil  (xFrac [i]));
147          }
148
149        // translate current NLP point+bounds into higher-dimensional space
150        initAuxs ();
151
152        if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
153          printf ("= RANK LEVEL = %d [%d] ==================================\n", rank, *rNum);
154          for (int i=0; i<nOrigVars_; i++)
155            if (Var (i) -> Multiplicity () > 0) // alive variable
156              printf ("#### %4d: %d %c %2d frac %20g -> int %20g [%20g,%20g]\n", 
157                      i, fixed [i], 
158                      variables_ [i] -> isInteger () ? 'I' : ' ',
159                      integerRank_ ? integerRank_ [i] : -1,
160                      xFrac [i], xInt [i], Lb (i), Ub (i));
161          printf ("--------------------\n");
162          for (int i=nOrigVars_; i<nVars (); i++)
163            if (Var (i) -> Multiplicity () > 0) // alive variable
164              printf ("#### %4d:   %c    frac %20g   [%20g,%20g]\n", 
165                      i, variables_ [i] -> isInteger () ? 'I' : ' ',
166                      //(integerRank_ && integerRank_ [i]) ? 'F' : ' ',
167                      X (i), Lb (i), Ub (i));
168          printf ("=================================================\n");
169        }
170
171        //CoinCopyN (xFrac, nOrigVars_, xInt);// TODO: re-copy first nOrigVars_ variables into xInt?
172
173        int 
174          remaining = *rNum,
175          ntrials   = 0, 
176          maxtrials = 3;// rNum / divider;
177
178        do {
179
180          bool one_fixed = false;
181
182          for (int i=0; i<nOrigVars_; i++) 
183
184            if ((Var (i) -> Multiplicity () > 0) && // alive
185                (integerRank_ [i] == rank)       && // at this rank
186                (fixed [i] == UNFIXED)           && // still to be fixed
187                Var (i) -> isInteger ()) {          // and integer
188
189              if (CoinCpuTime () > maxCpuTime_)
190                break;
191
192              // check if initAuxs() closed any bound; if so, fix variable
193              //if (Lb (i) == Ub (i)) { // doesn't work
194              if (ceil (Lb (i) - COUENNE_EPS) + COUENNE_EPS >= floor (Ub (i) + COUENNE_EPS)) {
195
196                X (i) = xInt [i] = ceil (Lb (i) - COUENNE_EPS);
197                fixed [i] = FIXED;
198                one_fixed = true;
199                --remaining;
200                continue;
201              }
202
203              // otherwise, test rounding up and down
204              int result = testIntFix (i, xFrac [i], fixed, xInt,
205                                       dualL, dualR, olb, oub, ntrials < maxtrials);
206
207              jnlst_ -> Printf (J_MOREVECTOR, J_NLPHEURISTIC, 
208                                "testing %d [%g -> %g], res = %d\n", i, xFrac [i], xInt [i], result);
209
210              if (result > 0) {
211                one_fixed = true;
212                --remaining;
213              } else if (result < 0) 
214                throw infeasible;
215            }
216
217          // if none fixed, fix first unfixed variable with this rank
218
219          if (!one_fixed) {
220
221            int index = 0;
222
223            // find first unfixed integer at this rank
224            while ((index < nOrigVars_) && 
225                   (!(Var (index) -> isInteger ()) ||
226                    (integerRank_ [index] != rank) ||
227                    (fixed [index] != UNFIXED)))
228              index++;
229
230            assert (index < nOrigVars_);
231
232            jnlst_ -> Printf (J_MOREVECTOR, J_NLPHEURISTIC, 
233                              "none fixed, fix %d from %g [%g,%g] [L=%g, R=%g]", 
234                              index, xFrac [index], Lb (index), Ub (index), 
235                              dualL [index], dualR [index]);
236
237            Lb (index) = Ub (index) = X (index) = xInt [index] = 
238              ((dualL [index] < dualR [index] - COUENNE_EPS) ? floor (xFrac [index]) :
239               (dualL [index] > dualR [index] + COUENNE_EPS) ? ceil  (xFrac [index]) :
240               ((CoinDrand48 () > xFrac [index] - floor (xFrac [index])) ? 
241                floor (xFrac [index]) : ceil (xFrac [index])));
242
243            jnlst_ -> Printf (J_MOREVECTOR, J_NLPHEURISTIC, " to %g\n", xInt [index]);
244
245            fixed [index] = FIXED;
246
247            --remaining;
248          }
249
250          ntrials++;
251
252          if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
253            printf ("--- remaining = %d --------------------------- \n", remaining);
254            for (int i=0; i<nOrigVars_; i++)
255              if (variables_ [i] -> Multiplicity () > 0)
256                printf ("#### %4d: %d %c %2d frac %20g -> int %20g  [%20g,%20g]\n", 
257                        i, fixed [i], 
258                        variables_ [i] -> isInteger () ? 'I' : ' ',
259                        integerRank_ ? integerRank_ [i] : -1,
260                        xFrac [i], xInt [i], Lb (i), Ub (i));
261            printf ("---------------------------\n");
262
263          }
264        } while (remaining > 0);
265      } // for
266
267    // save tightened bounds in NLP space. Sanity check
268    for (int i = nOrigVars_; i--;)
269      if (Var (i) -> Multiplicity () > 0) {
270
271        if (fixed [i] == FIXED)       // integer point, fixed
272          lb [i] = ub [i] = X (i) = xInt [i];
273
274        else if (Lb (i) > Ub (i))     // non-sense bounds, fix them
275          xInt [i] = X (i) = lb [i] = ub [i] = 
276            (fixed [i] == CONTINUOUS) ?
277                          (0.5 * (Lb (i) + Ub (i))) :
278            COUENNE_round (0.5 * (Lb (i) + Ub (i)));
279
280        else {                        // normal case
281          lb [i] = Lb (i);
282          ub [i] = Ub (i);
283          if      (xInt [i] < lb [i]) X (i) = xInt [i] = lb [i];
284          else if (xInt [i] > ub [i]) X (i) = xInt [i] = ub [i];
285        }
286      }
287
288    restoreUnusedOriginals (xInt);
289
290    // if initial point is feasible, compute corresponding objective
291    // and update if upper bound improves
292    initAuxs ();
293    int objind = Obj (0) -> Body () -> Index ();
294
295    if (X (objind) < getCutOff ()) {
296
297      const CouNumber *x = X ();
298      CouNumber xp = x [objind];
299
300      if (checkNLP (x, xp, true)) { // true for recomputing xp
301        setCutOff (xp);
302        jnlst_ -> Printf (J_DETAILED, J_NLPHEURISTIC, 
303                          "new cutoff from getIntCand: %g\n", xp);
304      }
305    }
306  } // try
307
308  catch (int i) {
309
310    if (i == infeasible)
311      retval = -1;
312  }
313
314  ////////////////////////////////////////////////////////////////////////////////
315
316  if (jnlst_->ProduceOutput(Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
317    if (retval >= 0) {
318      printf ("- Done: retval %d ----------------------------------------------------------------\n", 
319              retval);
320      for (int i=0; i<nOrigVars_; i++)
321        if (variables_ [i] -> Multiplicity () > 0)
322          printf ("#### %4d: %d %c frac %20g -> int %20g [%20g,%20g]\n", 
323                  i, fixed [i], variables_ [i] -> isInteger () ? 'I' : ' ',
324                  xFrac [i], xInt [i], lb [i], ub [i]);
325    } else printf ("no good point was found\n");
326  }
327
328  delete [] fixed;
329
330  delete [] olb;   delete [] oub;
331  delete [] dualL; delete [] dualR;
332
333  domain_.pop ();
334
335  jnlst_ -> Printf (J_MOREVECTOR, J_NLPHEURISTIC, "Done with GetIntegerCandidate\n");
336
337  optimum_ = store_optimum; // restore
338
339  return retval;
340}
Note: See TracBrowser for help on using the repository browser.