/home/coin/SVN-release/OS-2.4.2/Couenne/src/problem/getIntegerCandidate.cpp

Go to the documentation of this file.
00001 /* $Id: getIntegerCandidate.cpp 752 2011-08-08 03:45:07Z pbelotti $
00002  *
00003  * Name:    getIntegerCandidate.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: generate integer NLP point Y starting from fractional 
00006  *          solution using bound tightening
00007  *
00008  * (C) Carnegie-Mellon University, 2007-08.
00009  * This file is licensed under the Eclipse Public License (EPL)
00010  */
00011 
00012 #include "CoinTime.hpp"
00013 #include "CoinHelperFunctions.hpp"
00014 
00015 #include "CouenneProblem.hpp"
00016 #include "CouenneProblemElem.hpp"
00017 #include "CouenneExprVar.hpp"
00018 
00019 #include "CouenneRecordBestSol.hpp"
00020 
00021 using namespace Couenne;
00022 
00023 // lose patience after this many iterations of non-improving valid
00024 // tightening (step 1)
00025 #define VALID_ONLY_THRESHOLD 5 
00026 
00033 
00034 int CouenneProblem::getIntegerCandidate (const double *xFrac, double *xInt, 
00035                                          double *lb, double *ub) const {
00036   fillIntegerRank ();
00037 
00038   if (numberInRank_.size () == 0) // there is no original integer to fix
00039     return 0;
00040 
00041   CouNumber *store_optimum = optimum_; // temporary place for debug
00042                                        // optimal solution --- we
00043                                        // don't want to debug while
00044                                        // faking bounds, or we'd cut
00045                                        // the optimum all the time
00046   optimum_ = NULL;
00047 
00048   int
00049     ncols  = nVars (), 
00050     retval = 0;
00051 
00052   double
00053     *olb   = new double [ncols],       *oub   = new double [ncols],      // outer bounds
00054     *dualL = new double [nOrigVars_],  *dualR = new double [nOrigVars_]; // lb[objind] per fix index/direction
00055 
00056   // copy in current bounds
00057   CoinCopyN (Lb (), ncols, olb);
00058   CoinCopyN (Ub (), ncols, oub);
00059 
00060   // for now save fractional point into integer point
00061   CoinCopyN (xFrac, nOrigVars_ - ndefined_, xInt);
00062 
00063   domain_.push (nVars (), xInt, lb, ub);
00064 
00065   enum fixType *fixed = new enum fixType [nOrigVars_ - ndefined_]; // integer variables that were fixed
00066 
00067   for (int i=0; i<nOrigVars_ - ndefined_; i++) 
00068     fixed [i] = (Var (i) -> isDefinedInteger () && // work on integer variables only
00069                  Var (i) -> Multiplicity () > 0) ? // don't care if unused variable
00070       UNFIXED : CONTINUOUS;
00071 
00072   // A more sophisticated rounding procedure
00073   //
00074   // Call aggressive tightening for all integer variables, setting
00075   // each first at floor(x_i) and then at ceil(x_i).
00076   //
00077   //
00078   // for each value r of the rank: 1 .. maxrank {
00079   //
00080   //   restrict search of unfixed original variables with rank r to
00081   //   [floor,ceil] box around LP value
00082   //
00083   //   do {
00084   //
00085   //     threshold = nVarWithThisRank / k
00086   //     niter = 0;
00087   //     for each such (unfixed) variable {
00088   //
00089   //       try left
00090   //       try right
00091   //
00092   //       if      both infeasible, abort
00093   //       else if one infeasible, fix
00094   //       else if both feasible
00095   //         if niter < threshold skip
00096   //         else                 fix based on lb(left<>right)
00097   //     }
00098   //
00099   //     if no variables fixed, fix first unfixed
00100   //     ++niter
00101   //
00102   //   } while there are unfixed variables at this rank
00103   // }   
00104   //
00105   // check value of objective -> set cutoff and run bound tightening
00106   //
00107 
00108   const int infeasible = 1;
00109 
00110   try {
00111 
00112     // for all values of the integer rank
00113 
00114     jnlst_ -> Printf (Ipopt::J_ERROR, J_NLPHEURISTIC,
00115                       "Heuristic: looking for an initial point\n");
00116 
00117     if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
00118       printf ("=       ===========================================\n");
00119       printf ("= BEGIN ===========================================\n");
00120       printf ("=       ===========================================\n");
00121       for (int i=0; i<nOrigVars_; i++)
00122         if (variables_ [i] -> Multiplicity () > 0)
00123           printf ("#### %4d: %d %c %2d frac %20g  [%20g,%20g]\n", 
00124                   i, fixed [i], 
00125                   variables_ [i] -> isInteger () ? 'I' : ' ',
00126                   integerRank_ ? integerRank_ [i] : -1,
00127                   xFrac [i], Lb (i), Ub (i));
00128       printf ("---\n");
00129       for (int i=nOrigVars_; i<nVars (); i++)
00130         if (variables_ [i] -> Multiplicity () > 0)
00131           printf ("#### %4d:   %c    frac %20g   [%20g,%20g]\n", 
00132                   i, variables_ [i] -> isInteger () ? 'I' : ' ',
00133                   //(integerRank_ && integerRank_ [i]) ? 'F' : ' ',
00134                   X (i), Lb (i), Ub (i));
00135       printf ("===================================================\n");
00136       printf ("===================================================\n");
00137       printf ("===================================================\n");
00138     }
00139 
00140     const int 
00141       N_VARS_HUGE   = 10000,
00142       N_VARS_LARGE  = 1000,
00143       N_VARS_MEDIUM = 100,
00144       N_VARS_SMALL  = 10,
00145       N_VARS_TINY   = 5;
00146 
00147     int 
00148       ntrials   = 0, 
00149       nvars     = nVars (),
00150       maxtrials =
00151       (nvars >= N_VARS_HUGE)   ? 0 :
00152       (nvars >= N_VARS_LARGE)  ? 1 :
00153       (nvars >= N_VARS_MEDIUM) ? 2 :
00154       (nvars >= N_VARS_SMALL)  ? 4 :
00155       (nvars >= N_VARS_TINY)   ? 8 : 16;
00156 
00157     int rank = 1;
00158 
00159     for (std::vector <int>::iterator rNum = numberInRank_.begin(); 
00160          ++rNum != numberInRank_.end(); rank++) 
00161 
00162       if (*rNum > 0) {
00163 
00164         jnlst_ -> Printf (Ipopt::J_STRONGWARNING, J_NLPHEURISTIC,
00165                           "Analyzing %d variables with rank %d\n", *rNum, rank);
00166 
00167         if (CoinCpuTime () > maxCpuTime_)
00168           break;
00169 
00170         // *rNum is the number of variable with integer rank equal to rank
00171 
00172         // start restricting around current integer box
00173         for (int i=0; i<nOrigVars_ - ndefined_; i++) 
00174           if ((Var (i) -> Multiplicity () > 0) && // alive variable
00175               (Var (i) -> isDefinedInteger ()) && // integer, may fix if independent of other integers
00176               (integerRank_ [i] == rank)) {
00177 
00178             Lb (i) = CoinMax (Lb (i), floor (xFrac [i] - COUENNE_EPS)); 
00179             Ub (i) = CoinMin (Ub (i), ceil  (xFrac [i] + COUENNE_EPS));
00180           }
00181 
00182         // translate current NLP point+bounds into higher-dimensional space
00183         initAuxs ();
00184 
00185         if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
00186           printf ("= RANK LEVEL = %d [%d] ==================================\n", rank, *rNum);
00187           for (int i=0; i<nOrigVars_; i++)
00188             if (Var (i) -> Multiplicity () > 0) // alive variable
00189               printf ("#### %4d: %d %c %2d frac %20g -> int %20g [%20g,%20g]\n", 
00190                       i, fixed [i], 
00191                       variables_ [i] -> isInteger () ? 'I' : ' ',
00192                       integerRank_ ? integerRank_ [i] : -1,
00193                       xFrac [i], xInt [i], Lb (i), Ub (i));
00194           printf ("--------------------\n");
00195           for (int i=nOrigVars_; i<nVars (); i++)
00196             if (Var (i) -> Multiplicity () > 0) // alive variable
00197               printf ("#### %4d:   %c    frac %20g   [%20g,%20g]\n", 
00198                       i, variables_ [i] -> isInteger () ? 'I' : ' ',
00199                       //(integerRank_ && integerRank_ [i]) ? 'F' : ' ',
00200                       X (i), Lb (i), Ub (i));
00201           printf ("=================================================\n");
00202         }
00203 
00204         //CoinCopyN (xFrac, nOrigVars_, xInt);// TODO: re-copy first nOrigVars_ variables into xInt?
00205 
00206         int remaining = *rNum;
00207 
00208         do {
00209 
00210           bool one_fixed = false;
00211 
00212           for (int i=0; i<nOrigVars_ - ndefined_; i++) 
00213 
00214             if ((Var (i) -> Multiplicity () > 0) && // alive 
00215                 (integerRank_ [i] == rank)       && // at this rank
00216                 (fixed [i] == UNFIXED)           && // still to be fixed
00217                 Var (i) -> isDefinedInteger ()) {   // and integer
00218 
00219               if (CoinCpuTime () > maxCpuTime_)
00220                 break;
00221 
00222               // check if initAuxs() closed any bound; if so, fix variable
00223               //if (Lb (i) == Ub (i)) { // doesn't work
00224               if (ceil (Lb (i) - COUENNE_EPS) + COUENNE_EPS >= floor (Ub (i) + COUENNE_EPS)) {
00225 
00226                 X (i) = Lb (i) = Ub (i) = xInt [i] = ceil (Lb (i) - COUENNE_EPS);
00227                 fixed [i] = FIXED;
00228                 one_fixed = true;
00229                 --remaining;
00230                 continue;
00231               }
00232 
00233               // if variable already integer valued, fix it (TODO: do
00234               // this only if impatient, e.g. if #variables huge)
00235               if (ceil (xInt [i] - COUENNE_EPS) - COUENNE_EPS <= xInt [i]) {
00236 
00237                 X (i) = Lb (i) = Ub (i) = xInt [i] = ceil (xInt [i] - COUENNE_EPS);
00238                 fixed [i] = FIXED;
00239                 one_fixed = true;
00240                 --remaining;
00241                 continue;
00242               }
00243 
00244               // otherwise, test rounding up and down
00245               int result = testIntFix (i, xFrac [i], fixed, xInt,
00246                                        dualL, dualR, olb, oub, ntrials < maxtrials);
00247 
00248               jnlst_ -> Printf (Ipopt::J_STRONGWARNING, J_NLPHEURISTIC, 
00249                                 "testing %d [%g -> %g], res = %d\n", i, xFrac [i], xInt [i], result);
00250 
00251               if (result > 0) {
00252                 one_fixed = true;
00253                 --remaining;
00254               } else if (result < 0) 
00255                 throw infeasible;
00256             }
00257 
00258           // if none fixed, fix first unfixed variable with this rank
00259 
00260           if (!one_fixed) {
00261 
00262             int index = 0;
00263 
00264             // find first unfixed integer at this rank
00265             while ((index < nOrigVars_ - ndefined_) && 
00266                    (!(Var (index) -> isInteger ()) ||
00267                     (integerRank_ [index] != rank) ||
00268                     (fixed [index] != UNFIXED)))
00269               ++index;
00270 
00271             assert (index < nOrigVars_ - ndefined_);
00272 
00273             jnlst_ -> Printf (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC, 
00274                               "none fixed, fix %d from %g [%g,%g] [L=%g, R=%g]", 
00275                               index, xFrac [index], Lb (index), Ub (index), 
00276                               dualL [index], dualR [index]);
00277 
00278             Lb (index) = Ub (index) = X (index) = xInt [index] = 
00279               ((dualL [index] < dualR [index] - COUENNE_EPS) ? floor (xFrac [index]) :
00280                (dualL [index] > dualR [index] + COUENNE_EPS) ? ceil  (xFrac [index]) :
00281                ((CoinDrand48 () > xFrac [index] - floor (xFrac [index])) ? 
00282                 floor (xFrac [index]) : ceil (xFrac [index])));
00283 
00284             jnlst_ -> Printf (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC, " to %g\n", xInt [index]);
00285 
00286             fixed [index] = FIXED;
00287 
00288             --remaining;
00289           }
00290 
00291           ntrials++;
00292 
00293           if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
00294             printf ("--- remaining = %d --------------------------- \n", remaining);
00295             for (int i=0; i<nOrigVars_; i++)
00296               if (variables_ [i] -> Multiplicity () > 0)
00297                 printf ("#### %4d: %d %c %2d frac %20g -> int %20g  [%20g,%20g]\n", 
00298                         i, fixed [i], 
00299                         variables_ [i] -> isInteger () ? 'I' : ' ',
00300                         integerRank_ ? integerRank_ [i] : -1,
00301                         xFrac [i], xInt [i], Lb (i), Ub (i));
00302             printf ("---------------------------\n");
00303 
00304           }
00305 
00306           if (CoinCpuTime () > maxCpuTime_)
00307             break;
00308 
00309         } while (remaining > 0);
00310       } // for
00311 
00312     // save tightened bounds in NLP space. Sanity check
00313     for (int i = nOrigVars_ - ndefined_; i--;)
00314       if (Var (i) -> Multiplicity () > 0) {
00315 
00316         if (fixed [i] == FIXED)       // integer point, fixed
00317           lb [i] = ub [i] = X (i) = xInt [i];
00318 
00319         else if (Lb (i) > Ub (i)) {    // non-sense bounds, fix them
00320           xInt [i] = X (i) = lb [i] = ub [i] = 
00321             (fixed [i] == CONTINUOUS) ?
00322                           (0.5 * (Lb (i) + Ub (i))) :
00323             COUENNE_round (0.5 * (Lb (i) + Ub (i)));
00324 
00325           if (fixed [i] != CONTINUOUS)
00326             fixed [i] = FIXED;
00327 
00328         } else {                        // normal case
00329           lb [i] = Lb (i);
00330           ub [i] = Ub (i);
00331           if      (xInt [i] < lb [i]) X (i) = xInt [i] = lb [i];
00332           else if (xInt [i] > ub [i]) X (i) = xInt [i] = ub [i];
00333         }
00334       }
00335 
00336     restoreUnusedOriginals (xInt);
00337 
00338     // if initial point is feasible, compute corresponding objective
00339     // and update if upper bound improves
00340     initAuxs ();
00341     int objind = Obj (0) -> Body () -> Index ();
00342 
00343     if (X (objind) < getCutOff ()) {
00344 
00345       const CouNumber *x = X ();
00346       CouNumber xp = x [objind];
00347 
00348       if 
00349 #ifdef FM_CHECKNLP2
00350         (checkNLP2(x, 0, false, true, true, feas_tolerance_)) 
00351         // false for not caring about objective value, 
00352         // true for stopping at first violation and true for checkAll
00353 #else
00354         (checkNLP (x, xp, true)) // true for recomputing xp
00355 #endif   
00356           { 
00357             
00358 #ifdef FM_TRACE_OPTSOL
00359 #ifdef FM_CHECKNLP2
00360             recBSol->update();
00361             xp = recBSol->getVal();
00362 #else
00363             recBSol->update(x, nVars(), xp, feas_tolerance_);
00364 #endif
00365 #endif
00366             setCutOff (xp, x);
00367             jnlst_ -> Printf (Ipopt::J_DETAILED, J_NLPHEURISTIC, 
00368                               "new cutoff from getIntCand: %g\n", xp);
00369           }
00370     }
00371   } // try
00372 
00373   catch (int i) {
00374 
00375     if (i == infeasible)
00376       retval = -1;
00377   }
00378 
00380 
00381   if (jnlst_->ProduceOutput(Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
00382     if (retval >= 0) {
00383       printf ("- Done: retval %d ----------------------------------------------------------------\n", 
00384               retval);
00385       for (int i=0; i<nOrigVars_; i++)
00386         if (variables_ [i] -> Multiplicity () > 0)
00387           printf ("#### %4d: %d %c frac %20g -> int %20g [%20g,%20g]\n", 
00388                   i, fixed [i], variables_ [i] -> isInteger () ? 'I' : ' ',
00389                   xFrac [i], xInt [i], lb [i], ub [i]);
00390     } else printf ("no good point was found\n");
00391   }
00392 
00393   jnlst_ -> Printf (Ipopt::J_ERROR, J_NLPHEURISTIC, "Heuristic done\n");
00394 
00395 
00396   delete [] fixed;
00397 
00398   delete [] olb;   delete [] oub;
00399   delete [] dualL; delete [] dualR;
00400 
00401   domain_.pop ();
00402 
00403   jnlst_ -> Printf (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC, "Done with GetIntegerCandidate\n");
00404 
00405   optimum_ = store_optimum; // restore 
00406 
00407   return retval;
00408 }

Generated on Wed Nov 30 03:04:08 2011 by  doxygen 1.4.7