/home/coin/SVN-release/OS-2.1.1/Couenne/src/main/BonCouenneInterface.cpp

Go to the documentation of this file.
00001 /* $Id: BonCouenneInterface.cpp 242 2009-07-20 12:41:03Z stefan $ */
00002 // (C) Copyright International Business Machines Corporation (IBM) 2006, 2007
00003 // All Rights Reserved.
00004 // This code is published under the Common Public License.
00005 //
00006 // Authors :
00007 // Pietro Belotti, Carnegie Mellon University
00008 // Pierre Bonami, International Business Machines Corporation
00009 //
00010 // Date : 12/19/2006
00011 
00012 
00013 #include "BonCouenneInterface.hpp"
00014 #include "CoinHelperFunctions.hpp"
00015 #include "CouenneProblem.hpp"
00016 
00017 namespace Bonmin {
00018 
00020 CouenneInterface::CouenneInterface():
00021   AmplInterface(),
00022   have_nlp_solution_ (false)
00023 {}
00024 
00026 CouenneInterface::CouenneInterface(const CouenneInterface &other):
00027   AmplInterface(other),
00028   have_nlp_solution_ (false)
00029 {}
00030 
00032 CouenneInterface * CouenneInterface::clone(bool CopyData){
00033   return new CouenneInterface(*this);
00034 }
00035 
00037 CouenneInterface::~CouenneInterface(){
00038 }
00039 
00040 #ifdef COIN_HAS_ASL    
00041 void 
00042 CouenneInterface::readAmplNlFile(char **& argv, Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
00043                                  Ipopt::SmartPtr<Ipopt::OptionsList> options,
00044                                  Ipopt::SmartPtr<Ipopt::Journalist> journalist){
00045   //  if (!IsValid (app_))
00046   //createApplication (roptions, options, journalist, "couenne.");
00047   AmplInterface::readAmplNlFile(argv, roptions, options, journalist);
00048 }
00049 #endif
00050 
00063 void
00064 CouenneInterface::extractLinearRelaxation 
00065 (OsiSolverInterface &si, CouenneCutGenerator & couenneCg, bool getObj, bool solveNlp) {
00066 
00067   CouenneProblem *p = couenneCg.Problem ();
00068   bool is_feasible = true;
00069 
00070   if (solveNlp) {
00071 
00072     int nvars = p -> nVars();
00073 
00074     if (p -> doFBBT ()) {
00075 
00076       // include the rhs of auxiliary-based constraints into the FBBT
00077       // (should be useful with Vielma's problems, for example)
00078 
00079       for (int i=0; i < p -> nCons (); i++) {
00080 
00081         // for each constraint
00082         CouenneConstraint *con = p -> Con (i);
00083 
00084         // (which has an aux as its body)
00085         int index = con -> Body () -> Index ();
00086 
00087         if ((index >= 0) && (con -> Body () -> Type () == AUX)) {
00088 
00089           // if there exists violation, add constraint
00090           CouNumber 
00091             l = con -> Lb () -> Value (),       
00092             u = con -> Ub () -> Value ();
00093 
00094           // tighten bounds in Couenne's problem representation
00095           p -> Lb (index) = CoinMax (l, p -> Lb (index));
00096           p -> Ub (index) = CoinMin (u, p -> Ub (index));
00097         }
00098       }
00099 
00100       t_chg_bounds *chg_bds = new t_chg_bounds [nvars];
00101 
00102       for (int i=0; i<nvars; i++) {
00103         chg_bds [i].setLower(t_chg_bounds::CHANGED);
00104         chg_bds [i].setUpper(t_chg_bounds::CHANGED);
00105       }
00106 
00107       if (!(p -> boundTightening (chg_bds, NULL))) {
00108         is_feasible = false;
00109         *messageHandler() << "warning, tightened NLP is infeasible" << CoinMessageEol;
00110       }
00111 
00112       delete [] chg_bds;
00113 
00114       const double 
00115         *nlb = getColLower (),
00116         *nub = getColUpper ();
00117 
00118       for (int i=0; i < p -> nOrigVars () - p -> nDefVars (); i++) 
00119         if (p -> Var (i) -> Multiplicity () > 0) {
00120           /*printf ("---- %4d [%g,%g] [%g,%g]\n", i,
00121                   nlb [i], nub [i],
00122                   p -> Lb (i), p -> Ub (i));*/
00123           if (nlb [i] < p -> Lb (i) - COUENNE_EPS) setColLower (i, p -> Lb (i));
00124           if (nub [i] > p -> Ub (i) + COUENNE_EPS) setColUpper (i, p -> Ub (i));
00125         } else { 
00126           // if not enabled, fix them in the NLP solver
00127           setColLower (i, -COIN_DBL_MAX);
00128           setColUpper (i,  COIN_DBL_MAX);
00129         }
00130     } // ends FBBT part
00131 
00132     if (is_feasible) {
00133       try {
00134         initialSolve ();
00135       }
00136       catch (TNLPSolver::UnsolvedError *E) {
00137         // wrong, if NLP has problems this is not necessarily true...
00138         //is_feasible = false;
00139       }
00140     }
00141 
00142     if (!is_feasible) {
00143       OsiAuxInfo * auxInfo = si.getAuxiliaryInfo ();
00144       BabInfo * babInfo = dynamic_cast <BabInfo *> (auxInfo);
00145 
00146       if (babInfo) 
00147         babInfo -> setInfeasibleNode ();
00148     }
00149     
00150     if (is_feasible && isProvenOptimal ()) {
00151 
00152       CouNumber obj             = getObjValue    ();
00153       const CouNumber *solution = getColSolution ();
00154 
00155       if (getNumIntegers () > 0) {
00156 
00157         int
00158           norig = p -> nOrigVars (),
00159           nvars = p -> nVars ();
00160 
00161         bool fractional = false;
00162 
00163         // if problem is integer, check if any integral variable is
00164         // fractional. If so, round them and re-optimize
00165 
00166         for (int i=0; i<norig; i++)
00167           if (p  -> Var (i) -> isInteger () &&
00168               (p -> Var (i) -> Multiplicity () > 0) &&
00169               (!::isInteger (solution [i]))) {
00170             fractional = true;
00171             break;
00172           }
00173 
00174         if (fractional) { // try again if solution found by Ipopt is fractional
00175 
00176           double 
00177             *lbSave = new double [norig],
00178             *ubSave = new double [norig],
00179 
00180             *lbCur  = new double [nvars],
00181             *ubCur  = new double [nvars],
00182 
00183             *Y      = new double [nvars];
00184 
00185           CoinCopyN (getColLower (), norig, lbSave);
00186           CoinCopyN (getColUpper (), norig, ubSave);
00187 
00188           CoinFillN (Y,     nvars, 0.);
00189           CoinFillN (lbCur, nvars, -COUENNE_INFINITY);
00190           CoinFillN (ubCur, nvars,  COUENNE_INFINITY);
00191 
00192           CoinCopyN (getColLower (), norig, lbCur);
00193           CoinCopyN (getColUpper (), norig, ubCur);
00194 
00195           if (p -> getIntegerCandidate (solution, Y, lbCur, ubCur) >= 0) {
00196 
00197             for (int i = getNumCols (); i--;) {
00198 
00199               if (lbCur [i] > ubCur [i]) {
00200                 double swap = lbCur [i];
00201                 lbCur [i] = ubCur [i];
00202                 ubCur [i] = swap;
00203               }
00204 
00205               if      (Y [i] < lbCur [i]) Y [i] = lbCur [i];
00206               else if (Y [i] > ubCur [i]) Y [i] = ubCur [i];
00207             }
00208 
00209             for (int i=0; i<norig; i++)
00210               if ((p -> Var (i) -> Multiplicity () > 0) &&
00211                   p  -> Var (i) -> isInteger ()) {
00212                 setColLower (i, lbCur [i]);
00213                 setColUpper (i, ubCur [i]);
00214               }
00215 
00216             setColSolution (Y); // use initial solution given 
00217 
00218             try {
00219               resolve (); // solve with integer variables fixed
00220             }
00221             catch(TNLPSolver::UnsolvedError *E) {
00222             }
00223 
00224             //resolve (); 
00225 
00226             obj      = getObjValue ();
00227             solution = getColSolution ();
00228 
00229             // restore previous bounds on integer variables
00230             for (int i=0; i<norig; i++)
00231               if ((p -> Var (i) -> Multiplicity () > 0) &&
00232                   p  -> Var (i) -> isInteger ()) {
00233                 setColLower (i, lbSave [i]);
00234                 setColUpper (i, ubSave [i]);
00235               }
00236           }
00237 
00238           delete [] Y;
00239           delete [] lbSave;
00240           delete [] ubSave;
00241           delete [] lbCur;
00242           delete [] ubCur;
00243         } 
00244       }
00245 
00246       // re-check optimality in case resolve () was called
00247       if (isProvenOptimal () && 
00248           (obj < p -> getCutOff ()) && // first check (before re-computing)
00249           p -> checkNLP (solution, obj, true) && // true for recomputing obj
00250           (obj < p -> getCutOff ())) { // second check (real object might be different)
00251 
00252         // tell caller there is an initial solution to be fed to the initHeuristic
00253         have_nlp_solution_ = true;
00254 
00255         // set cutoff to take advantage of bound tightening
00256         //printf ("new cutoff %g from BonCouenneInterface\n", obj);
00257         p -> setCutOff (obj);
00258 
00259         OsiAuxInfo * auxInfo = si.getAuxiliaryInfo ();
00260         BabInfo * babInfo = dynamic_cast <BabInfo *> (auxInfo);
00261 
00262         if (babInfo) {
00263           babInfo -> setNlpSolution (solution, getNumCols (), obj);
00264           babInfo -> setHasNlpSolution (true);
00265         }
00266       }
00267     }
00268   }
00269 
00270   if (!is_feasible) // nothing else to do, problem infeasible
00271     return;
00272 
00273   int 
00274     numcols     = getNumCols (), // # original               variables
00275     numcolsconv = p -> nVars (); // # original + # auxiliary variables
00276 
00277   const double
00278     *lb = getColLower (),
00279     *ub = getColUpper ();
00280 
00281    // add original and auxiliary variables to the new problem
00282    for (int i=0; i<numcols; i++) 
00283      if (p -> Var (i) -> Multiplicity () > 0) si.addCol (0, NULL,NULL, lb [i],       ub [i],      0);
00284      else                                     si.addCol (0, NULL,NULL, -COIN_DBL_MAX,COIN_DBL_MAX,0);
00285    for (int i=numcols; i<numcolsconv; i++)    si.addCol (0, NULL,NULL, -COIN_DBL_MAX,COIN_DBL_MAX,0);
00286 
00287    // get initial relaxation
00288    OsiCuts cs;
00289    couenneCg.generateCuts (si, cs);
00290 
00291    // store all (original + auxiliary) bounds in the relaxation
00292    CouNumber * colLower = new CouNumber [numcolsconv];
00293    CouNumber * colUpper = new CouNumber [numcolsconv];
00294 
00295    CouNumber *ll = p -> Lb ();
00296    CouNumber *uu = p -> Ub ();
00297 
00298    // overwrite original bounds, could have improved within generateCuts
00299    for (int i = numcolsconv; i--;) 
00300      if (p -> Var (i) -> Multiplicity () > 0) {
00301        colLower [i] = (ll [i] > - COUENNE_INFINITY) ? ll [i] : -COIN_DBL_MAX;
00302        colUpper [i] = (uu [i] <   COUENNE_INFINITY) ? uu [i] :  COIN_DBL_MAX;
00303      } else {
00304        colLower [i] = -COIN_DBL_MAX;
00305        colUpper [i] =  COIN_DBL_MAX;
00306      }
00307 
00308    int numrowsconv = cs.sizeRowCuts ();
00309 
00310    // create matrix and other stuff
00311    CoinBigIndex * start = new CoinBigIndex [numrowsconv + 1];
00312 
00313    int    * length   = new int    [numrowsconv];
00314    double * rowLower = new double [numrowsconv];
00315    double * rowUpper = new double [numrowsconv];
00316 
00317    start[0] = 0;
00318    int nnz = 0;
00319    /* fill the four arrays. */
00320    for(int i = 0 ; i < numrowsconv ; i++)
00321    {
00322      OsiRowCut * cut = cs.rowCutPtr (i);
00323 
00324      const CoinPackedVector &v = cut->row();
00325      start[i+1] = start[i] + v.getNumElements();
00326      nnz += v.getNumElements();
00327      length[i] = v.getNumElements();
00328 
00329      rowLower[i] = cut->lb();
00330      rowUpper[i] = cut->ub();
00331    }
00332 
00333    assert (nnz == start [numrowsconv]);
00334    /* Now fill the elements arrays. */
00335    int * ind = new int[start[numrowsconv]];
00336    double * elem = new double[start[numrowsconv]];
00337    for(int i = 0 ; i < numrowsconv ; i++) {
00338 
00339      OsiRowCut * cut = cs.rowCutPtr (i);
00340 
00341      const CoinPackedVector &v = cut->row();
00342 
00343      if(v.getNumElements() != length[i])
00344        std::cout<<"Empty row"<<std::endl;
00345      //     cut->print();
00346      CoinCopyN (v.getIndices(),  length[i], ind  + start[i]);
00347      CoinCopyN (v.getElements(), length[i], elem + start[i]);
00348    }
00349 
00350    // Ok everything done now create interface
00351    CoinPackedMatrix A;
00352    A.assignMatrix(false, numcolsconv, numrowsconv,
00353                   start[numrowsconv], elem, ind,
00354                   start, length);
00355    if(A.getNumCols() != numcolsconv || A.getNumRows() != numrowsconv){
00356      std::cout<<"Error in row number"<<std::endl;
00357    }
00358    assert(A.getNumElements() == nnz);
00359    // Objective function
00360    double * obj = new double[numcolsconv];
00361    CoinFillN(obj,numcolsconv,0.);
00362 
00363    // some instances have no (or null) objective function, check it here
00364    if (p -> nObjs () > 0)
00365      p -> fillObjCoeff (obj);
00366 
00367    // Finally, load interface si with the initial LP relaxation
00368    si.loadProblem (A, colLower, colUpper, obj, rowLower, rowUpper);
00369 
00370    delete [] rowLower; 
00371    delete [] rowUpper;
00372    delete [] colLower;
00373    delete [] colUpper;
00374    delete [] obj;
00375 
00376    for (int i=0; i<numcolsconv; i++)
00377      if ((p -> Var (i) -> Multiplicity () > 0) &&
00378          (p -> Var (i) -> isDefinedInteger ()))
00379        si.setInteger (i);
00380  
00381    //si.writeMpsNative("toto",NULL,NULL,1);
00382    //si.writeLp ("toto");
00383    app_ -> enableWarmStart();
00384 
00385    // restored check. With "TOO FEW DEGREES OF FREEDOM" exception, x_sol() is null
00386    if (problem () -> x_sol ()) {
00387      setColSolution (problem () -> x_sol     ());
00388      setRowPrice    (problem () -> duals_sol ());
00389    }
00390 }
00391 
00392 
00394 void CouenneInterface::setAppDefaultOptions(Ipopt::SmartPtr<Ipopt::OptionsList> Options){
00395   Options->SetStringValue("bonmin.algorithm", "B-Couenne", true, true);
00396   Options->SetIntegerValue("bonmin.filmint_ecp_cuts", 1, true, true);
00397 }
00398 } 

Generated on Mon May 3 03:05:20 2010 by  doxygen 1.4.7