getIntegerCandidate.cpp

Go to the documentation of this file.

/* $Id: getIntegerCandidate.cpp 1257 2018-08-28 10:23:00Z pbelotti $
 *
 * Name:    getIntegerCandidate.cpp
 * Author:  Pietro Belotti
 * Purpose: generate integer NLP point Y starting from fractional 
 *          solution using bound tightening
 *
 * (C) Carnegie-Mellon University, 2007-08.
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "CoinTime.hpp"
#include "CoinHelperFunctions.hpp"

#include "CouenneProblem.hpp"
#include "CouenneProblemElem.hpp"
#include "CouenneExprVar.hpp"

#include "CouenneRecordBestSol.hpp"

using namespace Couenne;

// lose patience after this many iterations of non-improving valid
// tightening (step 1)
#define VALID_ONLY_THRESHOLD 5 


int CouenneProblem::getIntegerCandidate (const double *xFrac, double *xInt, 
                                         double *lb, double *ub) const {
  fillIntegerRank ();

  if (numberInRank_.size () == 0) // there is no original integer to fix
    return 0;

  CouNumber *store_optimum = optimum_; // temporary place for debug
                                       // optimal solution --- we
                                       // don't want to debug while
                                       // faking bounds, or we'd cut
                                       // the optimum all the time
  optimum_ = NULL;

  int
    ncols  = nVars (), 
    retval = 0;

  double
    *olb   = new double [ncols],       *oub   = new double [ncols],      // outer bounds
    *dualL = new double [nOrigVars_],  *dualR = new double [nOrigVars_]; // lb[objind] per fix index/direction

  // copy in current bounds
  CoinCopyN (Lb (), ncols, olb);
  CoinCopyN (Ub (), ncols, oub);

  // for now save fractional point into integer point
  CoinCopyN (xFrac, nOrigVars_ - ndefined_, xInt);

  domain_.push (nVars (), xInt, lb, ub);

  enum fixType *fixed = new enum fixType [nOrigVars_ - ndefined_]; // integer variables that were fixed

  for (int i=0; i<nOrigVars_ - ndefined_; i++) 
    fixed [i] = (Var (i) -> isDefinedInteger () && // work on integer variables only
                 Var (i) -> Multiplicity () > 0) ? // don't care if unused variable
      UNFIXED : CONTINUOUS;

  // A more sophisticated rounding procedure
  //
  // Call aggressive tightening for all integer variables, setting
  // each first at floor(x_i) and then at ceil(x_i).
  //
  //
  // for each value r of the rank: 1 .. maxrank {
  //
  //   restrict search of unfixed original variables with rank r to
  //   [floor,ceil] box around LP value
  //
  //   do {
  //
  //     threshold = nVarWithThisRank / k
  //     niter = 0;
  //     for each such (unfixed) variable {
  //
  //       try left
  //       try right
  //
  //       if      both infeasible, abort
  //       else if one infeasible, fix
  //       else if both feasible
  //         if niter < threshold skip
  //         else                 fix based on lb(left<>right)
  //     }
  //
  //     if no variables fixed, fix first unfixed
  //     ++niter
  //
  //   } while there are unfixed variables at this rank
  // }   
  //
  // check value of objective -> set cutoff and run bound tightening
  //

  const int infeasible = 1;

  try {

    // for all values of the integer rank

    jnlst_ -> Printf (Ipopt::J_ERROR, J_NLPHEURISTIC,
                      "Heuristic: looking for an initial point\n");

    if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
      printf ("=       ===========================================\n");
      printf ("= BEGIN ===========================================\n");
      printf ("=       ===========================================\n");
      for (int i=0; i<nOrigVars_; i++)
        if (variables_ [i] -> Multiplicity () > 0)
          printf ("#### %4d: %d %c %2d frac %20g  [%20g,%20g]\n", 
                  i, fixed [i], 
                  variables_ [i] -> isInteger () ? 'I' : ' ',
                  integerRank_ ? integerRank_ [i] : -1,
                  xFrac [i], Lb (i), Ub (i));
      printf ("---\n");
      for (int i=nOrigVars_; i<nVars (); i++)
        if (variables_ [i] -> Multiplicity () > 0)
          printf ("#### %4d:   %c    frac %20g   [%20g,%20g]\n", 
                  i, variables_ [i] -> isInteger () ? 'I' : ' ',
                  //(integerRank_ && integerRank_ [i]) ? 'F' : ' ',
                  X (i), Lb (i), Ub (i));
      printf ("===================================================\n");
      printf ("===================================================\n");
      printf ("===================================================\n");
    }

    const int 
      N_VARS_HUGE   = 10000,
      N_VARS_LARGE  = 1000,
      N_VARS_MEDIUM = 100,
      N_VARS_SMALL  = 10,
      N_VARS_TINY   = 5;

    int 
      ntrials   = 0, 
      nvars     = nVars (),
      maxtrials =
      (nvars >= N_VARS_HUGE)   ? 0 :
      (nvars >= N_VARS_LARGE)  ? 1 :
      (nvars >= N_VARS_MEDIUM) ? 2 :
      (nvars >= N_VARS_SMALL)  ? 4 :
      (nvars >= N_VARS_TINY)   ? 8 : 16;

    int rank = 1;

    for (std::vector <int>::iterator rNum = numberInRank_.begin(); 
         ++rNum != numberInRank_.end(); rank++) 

      if (*rNum > 0) {

        jnlst_ -> Printf (Ipopt::J_STRONGWARNING, J_NLPHEURISTIC,
                          "Analyzing %d variables with rank %d\n", *rNum, rank);

        if (CoinCpuTime () > maxCpuTime_)
          break;

        // *rNum is the number of variable with integer rank equal to rank

        // start restricting around current integer box
        for (int i=0; i<nOrigVars_ - ndefined_; i++) 
          if ((Var (i) -> Multiplicity () > 0) && // alive variable
              (Var (i) -> isDefinedInteger ()) && // integer, may fix if independent of other integers
              (integerRank_ [i] == rank)) {

            Lb (i) = CoinMax (Lb (i), floor (xFrac [i] - COUENNE_EPS)); 
            Ub (i) = CoinMin (Ub (i), ceil  (xFrac [i] + COUENNE_EPS));
          }

        // translate current NLP point+bounds into higher-dimensional space
        initAuxs ();

        if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
          printf ("= RANK LEVEL = %d [%d] ==================================\n", rank, *rNum);
          for (int i=0; i<nOrigVars_; i++)
            if (Var (i) -> Multiplicity () > 0) // alive variable
              printf ("#### %4d: %d %c %2d frac %20g -> int %20g [%20g,%20g]\n", 
                      i, fixed [i], 
                      variables_ [i] -> isInteger () ? 'I' : ' ',
                      integerRank_ ? integerRank_ [i] : -1,
                      xFrac [i], xInt [i], Lb (i), Ub (i));
          printf ("--------------------\n");
          for (int i=nOrigVars_; i<nVars (); i++)
            if (Var (i) -> Multiplicity () > 0) // alive variable
              printf ("#### %4d:   %c    frac %20g   [%20g,%20g]\n", 
                      i, variables_ [i] -> isInteger () ? 'I' : ' ',
                      //(integerRank_ && integerRank_ [i]) ? 'F' : ' ',
                      X (i), Lb (i), Ub (i));
          printf ("=================================================\n");
        }

        //CoinCopyN (xFrac, nOrigVars_, xInt);// TODO: re-copy first nOrigVars_ variables into xInt?

        int remaining = *rNum;

        do {

          bool one_fixed = false;

          for (int i=0; i<nOrigVars_ - ndefined_; i++) 

            if ((Var (i) -> Multiplicity () > 0) && // alive 
                (integerRank_ [i] == rank)       && // at this rank
                (fixed [i] == UNFIXED)           && // still to be fixed
                Var (i) -> isDefinedInteger ()) {   // and integer

              if (CoinCpuTime () > maxCpuTime_)
                break;

              // check if initAuxs() closed any bound; if so, fix variable
              //if (Lb (i) == Ub (i)) { // doesn't work
              if (ceil (Lb (i) - COUENNE_EPS) + COUENNE_EPS >= floor (Ub (i) + COUENNE_EPS)) {

                X (i) = Lb (i) = Ub (i) = xInt [i] = ceil (Lb (i) - COUENNE_EPS);
                fixed [i] = FIXED;
                one_fixed = true;
                --remaining;
                continue;
              }

              // if variable already integer valued, fix it (TODO: do
              // this only if impatient, e.g. if #variables huge)
              if (ceil (xInt [i] - COUENNE_EPS) - COUENNE_EPS <= xInt [i]) {

                X (i) = Lb (i) = Ub (i) = xInt [i] = ceil (xInt [i] - COUENNE_EPS);
                fixed [i] = FIXED;
                one_fixed = true;
                --remaining;
                continue;
              }

              // otherwise, test rounding up and down
              int result = testIntFix (i, xFrac [i], fixed, xInt,
                                       dualL, dualR, olb, oub, ntrials < maxtrials);

              jnlst_ -> Printf (Ipopt::J_STRONGWARNING, J_NLPHEURISTIC, 
                                "testing %d [%g -> %g], res = %d\n", i, xFrac [i], xInt [i], result);

              if (result > 0) {
                one_fixed = true;
                --remaining;
              } else if (result < 0) 
                throw infeasible;
            }

          // if none fixed, fix first unfixed variable with this rank

          if (!one_fixed) {

            int index = 0;

            // find first unfixed integer at this rank
            while ((index < nOrigVars_ - ndefined_) && 
                   (!(Var (index) -> isInteger ()) ||
                    (integerRank_ [index] != rank) ||
                    (fixed [index] != UNFIXED)))
              ++index;

            assert (index < nOrigVars_ - ndefined_);

            jnlst_ -> Printf (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC, 
                              "none fixed, fix %d from %g [%g,%g] [L=%g, R=%g]", 
                              index, xFrac [index], Lb (index), Ub (index), 
                              dualL [index], dualR [index]);

            Lb (index) = Ub (index) = X (index) = xInt [index] = 
              ((dualL [index] < dualR [index] - COUENNE_EPS) ? floor (xFrac [index]) :
               (dualL [index] > dualR [index] + COUENNE_EPS) ? ceil  (xFrac [index]) :
               ((CoinDrand48 () > xFrac [index] - floor (xFrac [index])) ? 
                floor (xFrac [index]) : ceil (xFrac [index])));

            jnlst_ -> Printf (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC, " to %g\n", xInt [index]);

            fixed [index] = FIXED;

            --remaining;
          }

          ntrials++;

          if (jnlst_ -> ProduceOutput (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
            printf ("--- remaining = %d --------------------------- \n", remaining);
            for (int i=0; i<nOrigVars_; i++)
              if (variables_ [i] -> Multiplicity () > 0)
                printf ("#### %4d: %d %c %2d frac %20g -> int %20g  [%20g,%20g]\n", 
                        i, fixed [i], 
                        variables_ [i] -> isInteger () ? 'I' : ' ',
                        integerRank_ ? integerRank_ [i] : -1,
                        xFrac [i], xInt [i], Lb (i), Ub (i));
            printf ("---------------------------\n");

          }

          if (CoinCpuTime () > maxCpuTime_)
            break;

        } while (remaining > 0);
      } // for

    // save tightened bounds in NLP space. Sanity check
    for (int i = nOrigVars_ - ndefined_; i--;)
      if (Var (i) -> Multiplicity () > 0) {

        if (fixed [i] == FIXED)       // integer point, fixed
          lb [i] = ub [i] = X (i) = xInt [i];

        else if (Lb (i) > Ub (i)) {    // non-sense bounds, fix them
          xInt [i] = X (i) = lb [i] = ub [i] = 
            (fixed [i] == CONTINUOUS) ?
                          (0.5 * (Lb (i) + Ub (i))) :
            COUENNE_round (0.5 * (Lb (i) + Ub (i)));

          if (fixed [i] != CONTINUOUS)
            fixed [i] = FIXED;

        } else {                        // normal case
          lb [i] = Lb (i);
          ub [i] = Ub (i);
          if      (xInt [i] < lb [i]) X (i) = xInt [i] = lb [i];
          else if (xInt [i] > ub [i]) X (i) = xInt [i] = ub [i];
        }
      }

    restoreUnusedOriginals (xInt);

    // if initial point is feasible, compute corresponding objective
    // and update if upper bound improves
    initAuxs ();
    int objind = Obj (0) -> Body () -> Index ();

    CouNumber xp = objind >= 0 ? X (objind) : Obj (0) -> Body () -> Value ();

    const CouNumber *x = X ();

    if  (checkNLP0 (x, xp, true, false, true, true)) {
// #ifdef FM_CHECKNLP2
//      (checkNLP2(x, 0, false, true, true, feas_tolerance_)) 
//      // false for not caring about objective value, 
//      // true for stopping at first violation and true for checkAll
// #else
//      (checkNLP (x, xp, true)) // true for recomputing xp
// #endif        
      
#ifdef FM_TRACE_OPTSOL
#ifdef FM_CHECKNLP2
      recBSol->update();
      xp = recBSol->getVal();
#else
      recBSol -> update (x, nVars(), xp, feas_tolerance_);
#endif
#endif
      if (xp < getCutOff ()) {

        setCutOff (xp, x);
        jnlst_ -> Printf (Ipopt::J_DETAILED, J_NLPHEURISTIC, 
                          "new cutoff from getIntCand: %g\n", xp);
      }
    }
  } // try

  catch (int i) {

    if (i == infeasible)
      retval = -1;

    //
    // While Couenne does not use this solution (retval==-1), the FP can still use a rounded solution
    //

    for (int j=nVars(); j--;)
      xInt [j] = ceil (xInt [j] - 0.5);
  }


  if (jnlst_->ProduceOutput(Ipopt::J_MOREVECTOR, J_NLPHEURISTIC)) {
    if (retval >= 0) {
      printf ("- Done: retval %d ----------------------------------------------------------------\n", 
              retval);
      for (int i=0; i<nOrigVars_; i++)
        if (variables_ [i] -> Multiplicity () > 0)
          printf ("#### %4d: %d %c frac %20g -> int %20g [%20g,%20g]\n", 
                  i, fixed [i], variables_ [i] -> isInteger () ? 'I' : ' ',
                  xFrac [i], xInt [i], lb [i], ub [i]);
    } else printf ("no good point was found\n");
  }

  jnlst_ -> Printf (Ipopt::J_ERROR, J_NLPHEURISTIC, "Heuristic done\n");


  delete [] fixed;

  delete [] olb;   delete [] oub;
  delete [] dualL; delete [] dualR;

  domain_.pop ();

  jnlst_ -> Printf (Ipopt::J_MOREVECTOR, J_NLPHEURISTIC, "Done with GetIntegerCandidate\n");

  optimum_ = store_optimum; // restore 

  return retval;
}