BonCouenneInterface.cpp

Go to the documentation of this file.

/* $Id: BonCouenneInterface.cpp 997 2013-08-11 23:21:38Z pbelotti $ */
// (C) Copyright International Business Machines Corporation (IBM) 2006, 2007
// All Rights Reserved.
// This code is published under the Eclipse Public License (EPL).
//
// Authors :
// Pietro Belotti, Carnegie Mellon University
// Pierre Bonami, International Business Machines Corporation
//
// Date : 12/19/2006


#include "BonCouenneInterface.hpp"
#include "CoinHelperFunctions.hpp"

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

using namespace Couenne;

CouenneInterface::CouenneInterface():
  AmplInterface(),
  have_nlp_solution_ (false)
{}

CouenneInterface::CouenneInterface(const CouenneInterface &other):
  AmplInterface(other),
  have_nlp_solution_ (false)
{}

CouenneInterface * CouenneInterface::clone(bool CopyData){
  return new CouenneInterface(*this);
}

CouenneInterface::~CouenneInterface(){
}

#ifdef COUENNEINTERFACE_FROM_ASL
void 
CouenneInterface::readAmplNlFile(char **& argv, Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
                                 Ipopt::SmartPtr<Ipopt::OptionsList> options,
                                 Ipopt::SmartPtr<Ipopt::Journalist> journalist){
  //  if (!IsValid (app_))
  //createApplication (roptions, options, journalist, "couenne.");
  AmplInterface::readAmplNlFile(argv, roptions, options, journalist);
}
#endif

void
CouenneInterface::extractLinearRelaxation 
(OsiSolverInterface &si, CouenneCutGenerator & couenneCg, bool getObj, bool solveNlp) {

  {
    int nlpLogLevel;
    options () -> GetIntegerValue ("nlp_log_level", nlpLogLevel, "couenne.");
    messageHandler () -> setLogLevel (nlpLogLevel);
  }

  CouenneProblem *p = couenneCg.Problem ();
  bool is_feasible = true;

  if (solveNlp) {

    int nvars = p -> nVars();

    if (p -> doFBBT ()) {

      // include the rhs of auxiliary-based constraints into the FBBT
      // (should be useful with Vielma's problems, for example)

      for (int i=0; i < p -> nCons (); i++) {

        // for each constraint
        CouenneConstraint *con = p -> Con (i);

        // (which has an aux as its body)
        int index = con -> Body () -> Index ();

        if ((index >= 0) && (con -> Body () -> Type () == AUX)) {

          // if there exists violation, add constraint
          CouNumber 
            l = con -> Lb () -> Value (),       
            u = con -> Ub () -> Value ();

          // tighten bounds in Couenne's problem representation
          p -> Lb (index) = CoinMax (l, p -> Lb (index));
          p -> Ub (index) = CoinMin (u, p -> Ub (index));
        }
      }

      t_chg_bounds *chg_bds = new t_chg_bounds [nvars];

      for (int i=0; i<nvars; i++) {
        chg_bds [i].setLower(t_chg_bounds::CHANGED);
        chg_bds [i].setUpper(t_chg_bounds::CHANGED);
      }

      if (!(p -> boundTightening (chg_bds, CglTreeInfo (), NULL))) {
        is_feasible = false;
        //*messageHandler() << "Couenne: Warning, tightened NLP is infeasible" << CoinMessageEol;
      }

      delete [] chg_bds;

      const double 
        *nlb = getColLower (),
        *nub = getColUpper ();

      for (int i=0; i < p -> nOrigVars () - p -> nDefVars (); i++) 
        if (p -> Var (i) -> Multiplicity () > 0) {
          /*printf ("---- %4d [%g,%g] [%g,%g]\n", i,
                  nlb [i], nub [i],
                  p -> Lb (i), p -> Ub (i));*/

          double 
            lower = nlb [i],
            upper = nub [i];

          if (lower       > upper)       CoinSwap (lower,       upper);
          if (p -> Lb (i) > p -> Ub (i)) CoinSwap (p -> Lb (i), p -> Ub (i));

          if (lower < p -> Lb (i) - COUENNE_EPS) setColLower (i, CoinMin (nub[i], p -> Lb (i)));
          if (upper > p -> Ub (i) + COUENNE_EPS) setColUpper (i, CoinMax (nlb[i], p -> Ub (i)));

        } else { 
          // if not enabled, fix them in the NLP solver
          setColLower (i, -COIN_DBL_MAX);
          setColUpper (i,  COIN_DBL_MAX);
        }
    } // ends FBBT part

    if (is_feasible) {
      try {
        options () -> SetNumericValue ("max_cpu_time", CoinMax (0.1, (p -> getMaxCpuTime () - CoinCpuTime ()) / 2));

        initialSolve ();

        if (isDualObjectiveLimitReached() &&
            (getNumIntegers () == 0))
          *messageHandler () << "Couenne: Warning, NLP is unbounded" << CoinMessageEol;
      }
      catch (Bonmin::TNLPSolver::UnsolvedError *E) {
        // wrong, if NLP has problems this is not necessarily true...
        //is_feasible = false;
      }
    }

    if (!is_feasible) {
      OsiAuxInfo * auxInfo = si.getAuxiliaryInfo ();
      Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (auxInfo);

      if (babInfo) 
        babInfo -> setInfeasibleNode ();
    }
    
    if (is_feasible && isProvenOptimal ()) {

      CouNumber obj             = getObjValue    ();
      const CouNumber *solution = getColSolution ();

      if (getNumIntegers () > 0) {

        int
          norig = p -> nOrigVars () - p -> nDefVars (),
          nvars = p -> nVars ();

        bool fractional = false;

        // if problem is integer, check if any integral variable is
        // fractional. If so, round them and re-optimize

        for (int i=0; i<norig; i++)
          if ((p -> Var (i) -> Multiplicity () > 0) &&
              p  -> Var (i) -> isDefinedInteger () &&
              (!::isInteger (solution [i]))) {
            fractional = true;
            break;
          }

        if (fractional) { // try again if solution found by Ipopt is fractional

          double 
            *lbSave = new double [norig],
            *ubSave = new double [norig],

            *lbCur  = new double [nvars],
            *ubCur  = new double [nvars],

            *Y      = new double [nvars];

          CoinCopyN (getColLower (), norig, lbSave);
          CoinCopyN (getColUpper (), norig, ubSave);

          CoinFillN (Y,     nvars, 0.);
          CoinFillN (lbCur, nvars, -COUENNE_INFINITY);
          CoinFillN (ubCur, nvars,  COUENNE_INFINITY);

          CoinCopyN (getColLower (), norig, lbCur);
          CoinCopyN (getColUpper (), norig, ubCur);

          if (p -> getIntegerCandidate (solution, Y, lbCur, ubCur) >= 0) {

            for (int i = getNumCols (); i--;) {

              if (lbCur [i] > ubCur [i])
                CoinSwap (lbCur [i], ubCur [i]);

              if      (Y [i] < lbCur [i]) Y [i] = lbCur [i];
              else if (Y [i] > ubCur [i]) Y [i] = ubCur [i];
            }

            for (int i=0; i<norig; i++)
              if ((p -> Var (i) -> Multiplicity () > 0) &&
                  p  -> Var (i) -> isDefinedInteger ()) {

                setColLower (i, lbCur [i]);
                setColUpper (i, ubCur [i]);
              }

            setColSolution (Y); // use initial solution given 

            try {
              options () -> SetNumericValue ("max_cpu_time", CoinMax (0.1, p -> getMaxCpuTime () - CoinCpuTime ()));

              resolve (); // solve with integer variables fixed

              if (isDualObjectiveLimitReached() &&
                  (getNumIntegers () == 0))
                *messageHandler () << "Couenne: Warning, NLP is is unbounded" << CoinMessageEol;
            }
            catch (Bonmin::TNLPSolver::UnsolvedError *E) {
            }

            //resolve (); 

            obj      = getObjValue ();
            solution = getColSolution ();

            // restore previous bounds on integer variables
            for (int i=0; i<norig; i++)
              if ((p -> Var (i) -> Multiplicity () > 0) &&
                  p  -> Var (i) -> isDefinedInteger ()) {

                if (lbSave [i] > ubSave [i])
                  CoinSwap (lbSave [i], ubSave [i]);

                setColLower (i, lbSave [i]);
                setColUpper (i, ubSave [i]);
              }
          }

          delete [] Y;
          delete [] lbSave;
          delete [] ubSave;
          delete [] lbCur;
          delete [] ubCur;
        }
      }

      // re-check optimality in case resolve () was called
      if (isProvenOptimal () && 
          //      (obj < p -> getCutOff ()) && // check #1 (before re-computing) -- BUG. What if real object is actually better?

#ifdef FM_CHECKNLP2
          (p->checkNLP2(solution, 0, false, true, false, p->getFeasTol())) &&
          (p->getRecordBestSol()->getModSolVal() < p->getCutOff())
#else
          p -> checkNLP (solution, obj, true) && // true for recomputing obj
          (obj < p -> getCutOff ())
#endif
          ) {           // check #2 (real object might be different)

        // tell caller there is an initial solution to be fed to the initHeuristic
        have_nlp_solution_ = true;

        // set cutoff to take advantage of bound tightening

#ifdef FM_CHECKNLP2
        obj = p->getRecordBestSol()->getModSolVal();
#endif

        p -> setCutOff (obj, solution);

        OsiAuxInfo * auxInfo = si.getAuxiliaryInfo ();
        Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (auxInfo);

        if (babInfo) {

#ifdef FM_CHECKNLP2
          babInfo -> setNlpSolution (p->getRecordBestSol()->modSol, 
                                     getNumCols(), obj);
#else
          babInfo -> setNlpSolution (solution, getNumCols (), obj);
#endif
          babInfo -> setHasNlpSolution (true);
        }

#ifdef FM_TRACE_OPTSOL
#ifdef FM_CHECKNLP2
        p->getRecordBestSol()->update();
#else
        p->getRecordBestSol()->update(solution, getNumCols(), 
                                      obj, p->getFeasTol());
#endif
#endif

      }
    } else {


    }
  }

  if (!is_feasible) // nothing else to do, problem infeasible
    return;

  int 
    numcols     = p -> nOrigVars (), // # original               variables
    numcolsconv = p -> nVars     (); // # original + # auxiliary variables

  const double
    *lb = p -> Lb (), //getColLower (),
    *ub = p -> Ub (); //getColUpper ();

   // add original and auxiliary variables to the new problem
   for (int i=0; i<numcols; i++) 
     if (p -> Var (i) -> Multiplicity () > 0) si.addCol (0, NULL,NULL, lb [i],       ub [i],      0);
     else                                     si.addCol (0, NULL,NULL, -COIN_DBL_MAX,COIN_DBL_MAX,0);
   for (int i=numcols; i<numcolsconv; i++)    si.addCol (0, NULL,NULL, -COIN_DBL_MAX,COIN_DBL_MAX,0);

   // get initial relaxation
   OsiCuts cs;
   couenneCg.generateCuts (si, cs);

   // store all (original + auxiliary) bounds in the relaxation
   CouNumber * colLower = new CouNumber [numcolsconv];
   CouNumber * colUpper = new CouNumber [numcolsconv];

   CouNumber *ll = p -> Lb ();
   CouNumber *uu = p -> Ub ();

   // overwrite original bounds, could have improved within generateCuts
   for (int i = numcolsconv; i--;) 
     if (p -> Var (i) -> Multiplicity () > 0) {
       colLower [i] = (ll [i] > - COUENNE_INFINITY) ? ll [i] : -COIN_DBL_MAX;
       colUpper [i] = (uu [i] <   COUENNE_INFINITY) ? uu [i] :  COIN_DBL_MAX;
     } else {
       colLower [i] = -COIN_DBL_MAX;
       colUpper [i] =  COIN_DBL_MAX;
     }

   int numrowsconv = cs.sizeRowCuts ();

   // create matrix and other stuff
   CoinBigIndex * start = new CoinBigIndex [numrowsconv + 1];

   int    * length   = new int    [numrowsconv];
   double * rowLower = new double [numrowsconv];
   double * rowUpper = new double [numrowsconv];

   start[0] = 0;
   int nnz = 0;
   /* fill the four arrays. */
   for(int i = 0 ; i < numrowsconv ; i++)
   {
     OsiRowCut * cut = cs.rowCutPtr (i);

     const CoinPackedVector &v = cut->row();
     start[i+1] = start[i] + v.getNumElements();
     nnz += v.getNumElements();
     length[i] = v.getNumElements();

     rowLower[i] = cut->lb();
     rowUpper[i] = cut->ub();
   }

   assert (nnz == start [numrowsconv]);
   /* Now fill the elements arrays. */
   int * ind = new int[start[numrowsconv]];
   double * elem = new double[start[numrowsconv]];
   for(int i = 0 ; i < numrowsconv ; i++) {

     OsiRowCut * cut = cs.rowCutPtr (i);

     const CoinPackedVector &v = cut->row();

     if(v.getNumElements() != length[i])
       std::cout<<"Empty row"<<std::endl;
     //     cut->print();
     CoinCopyN (v.getIndices(),  length[i], ind  + start[i]);
     CoinCopyN (v.getElements(), length[i], elem + start[i]);
   }

   // Ok everything done now create interface
   CoinPackedMatrix A;
   A.assignMatrix(false, numcolsconv, numrowsconv,
                  start[numrowsconv], elem, ind,
                  start, length);
   if(A.getNumCols() != numcolsconv || A.getNumRows() != numrowsconv){
     std::cout<<"Error in row number"<<std::endl;
   }
   assert(A.getNumElements() == nnz);
   // Objective function
   double * obj = new double[numcolsconv];
   CoinFillN(obj,numcolsconv,0.);

   // some instances have no (or null) objective function, check it here
   if (p -> nObjs () > 0)
     p -> fillObjCoeff (obj);

   // Finally, load interface si with the initial LP relaxation
   si.loadProblem (A, colLower, colUpper, obj, rowLower, rowUpper);

   delete [] rowLower; 
   delete [] rowUpper;
   delete [] colLower;
   delete [] colUpper;
   delete [] obj;

   for (int i=0; i<numcolsconv; i++)
     if ((p -> Var (i) -> Multiplicity () > 0) &&
         (p -> Var (i) -> isDefinedInteger ()))
       si.setInteger (i);
 
   //si.writeMpsNative("toto",NULL,NULL,1);
   //si.writeLp ("toto");
   app_ -> enableWarmStart();

   // restored check. With "TOO FEW DEGREES OF FREEDOM" exception, x_sol() is null
   if (problem () -> x_sol ()) {
     setColSolution (problem () -> x_sol     ());
     setRowPrice    (problem () -> duals_sol ());
   }
}


void CouenneInterface::setAppDefaultOptions(Ipopt::SmartPtr<Ipopt::OptionsList> Options){
  Options->SetStringValue("bonmin.algorithm", "B-Couenne", true, true);
  Options->SetIntegerValue("bonmin.filmint_ecp_cuts", 1, true, true);
}