Changeset 883


Ignore:
Timestamp:
Aug 3, 2012 9:39:53 AM (7 years ago)
Author:
pbelotti
Message:

adding features to feasibility pump. Checking variables with multiplicity 0

Location:
trunk/Couenne/src
Files:
1 added
13 edited

Legend:

Unmodified
Added
Removed
  • trunk/Couenne/src/expression/partial/CouenneExprHess.cpp

    r716 r883  
    189189  /// for each variable, fill a row of the hessian
    190190  for (int i=0; i < nVars; i++) {
     191
     192    if (p -> Var (i) -> Multiplicity () <= 0)
     193      continue;
    191194
    192195    // create dense row. These will become numL later
  • trunk/Couenne/src/expression/partial/CouenneExprJac.cpp

    r716 r883  
    141141    for (std::set <int>::iterator k = deplist.begin (); k != deplist.end (); ++k) {
    142142
     143      // highly unlikely that x_k shows up in c's body, but you never know...
     144      if (p -> Var (*k) -> Multiplicity () <= 0)
     145        continue;
     146
    143147      expression
    144148        *J = c -> Body () -> differentiate (*k), // derivative of the
     
    201205
    202206    for (std::set <int>::iterator k = deplist.begin (); k != deplist.end (); ++k) {
     207
     208      if (p -> Var (*k) -> Multiplicity () <= 0)
     209        continue;
    203210
    204211      expression
  • trunk/Couenne/src/heuristics/CouenneFPFindSolution.cpp

    r733 r883  
    1414#include "CouenneFPpool.hpp"
    1515#include "CouenneProblem.hpp"
     16#include "CouenneExprVar.hpp"
     17
    1618#include "cons_rowcuts.h"
    1719
    1820#ifdef COIN_HAS_SCIP
    19 /* general SCIP includes */
    2021#include "scip/scip.h"
    21 #include "scip/cons_linear.h"
    22 #include "scip/scipdefplugins.h"
    2322#endif
    2423
     
    2827double CouenneFeasPump::findSolution (double* &sol, int niter, int* nsuciter) {
    2928
    30   /// as found on the notes, these methods can be used, from the most
     29  /// As found on the notes, these methods can be used, from the most
    3130  /// expensive and accurate (exact) method to a cheap, inexact one:
    3231  ///
     
    6261#ifdef COIN_HAS_SCIP
    6362
    64   static int currentmilpmethod = 0;
    65 
    66   int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
    67 
    6863  if (useSCIP_ && problem_ -> nIntVars () > 0) { // if LP, use Clp below
    6964
    70      SCIP* scip;
     65    SCIP_RETCODE retcode = ScipSolve (sol, niter, nsuciter);
    7166
    72      SCIP_VAR** vars;
    73      const SCIP_Real* lbs;
    74      const SCIP_Real* ubs;
    75      const SCIP_Real* objs;
    76      const char* vartypes;
    77      const CoinPackedMatrix * matrix;
    78      const CoinBigIndex* rowstarts;
    79      const int* rowlengths;
    80      const SCIP_Real* coeffs;
    81      const SCIP_Real* lhss;
    82      const SCIP_Real* rhss;
    83      const int* indices;
     67    if (retcode != SCIP_OKAY) {
    8468
    85      SCIP_Real timelimit;
     69      printf ("SCIP did not return successfully\n");
     70      return 1.e40;
     71    }
     72  } else
    8673
    87      double infinity;
    88      int nvars;
    89      int nconss;
    90      int nscipsols;
    91 
    92      bool solveagain;
    93 
    94      // COUENNE_INFINITY , getInfinity()
    95 
    96      // get problem data
    97      nvars    = milp_ -> getNumCols ();
    98      nconss   = milp_ -> getNumRows ();
    99      infinity = milp_ -> getInfinity ();
    100 
    101      // get variable data
    102      lbs =      milp_ -> getColLower ();
    103      ubs =      milp_ -> getColUpper ();
    104      objs =     milp_ -> getObjCoefficients ();
    105      vartypes = milp_ -> getColType ();
    106 
    107      // get row data
    108      lhss = milp_ -> getRowLower ();
    109      rhss = milp_ -> getRowUpper ();
    110 
    111      // get matrix data
    112      matrix     = milp_ -> getMatrixByRow();
    113      rowstarts  = matrix -> getVectorStarts();
    114      rowlengths = matrix -> getVectorLengths();
    115      coeffs     = matrix -> getElements();
    116      indices    = matrix -> getIndices();
    117      
    118      if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC)) {
    119         SCIPdebugMessage("create SCIP problem instance with %d variables and %d constraints.\n", nvars, nconss);
    120      }
    121 
    122      // initialize SCIP
    123      SCIP_CALL_ABORT( SCIPcreate(&scip) );
    124      assert(scip != NULL);
    125    
    126      // include default SCIP plugins
    127      SCIP_CALL_ABORT( SCIPincludeDefaultPlugins(scip) );
    128 
    129      // include row cut constraint hanlder
    130      if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
    131      {
    132         SCIP_CALL_ABORT( SCIPincludeConshdlrRowcuts(scip, couenneCG_, milp_) );
    133      }
    134 
    135      // create problem instance in SCIP
    136      SCIP_CALL_ABORT( SCIPcreateProb(scip, "auxiliary FeasPump MILP", NULL, NULL, NULL, NULL, NULL, NULL, NULL) );
    137 
    138      // allocate local memory for SCIP variable array
    139      SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &vars, nvars) );
    140    
    141      // one variable for objective !!!!!!!!!
    142 
    143      // create variables
    144      for (int i=0; i<nvars; i++) {
    145         char varname[SCIP_MAXSTRLEN]; 
    146 
    147         // check that all data is in valid ranges
    148         assert( 0 <= vartypes[i] && vartypes[i] <= 2);
    149         checkInfinity(scip, lbs[i], infinity);
    150         checkInfinity(scip, ubs[i], infinity);
    151        
    152         // all variables are named x_i
    153         (void) SCIPsnprintf(varname, SCIP_MAXSTRLEN, "x_%d", i);
    154         SCIP_CALL_ABORT( SCIPcreateVar(scip, &vars[i], varname, lbs[i], ubs[i], objs[i],
    155               vartypes[i] == 0 ? SCIP_VARTYPE_CONTINUOUS : (vartypes[i] == 1 ? SCIP_VARTYPE_BINARY : SCIP_VARTYPE_INTEGER),
    156               TRUE, FALSE, NULL, NULL, NULL, NULL, NULL) );
    157 
    158         // add the variable to SCIP
    159         SCIP_CALL_ABORT( SCIPaddVar(scip, vars[i]) );
    160      }
    161 
    162      // create constraints
    163      for (int i=0; i<nconss; i++) {
    164 
    165         SCIP_CONS* cons;
    166        
    167         char consname[SCIP_MAXSTRLEN]; 
    168         (void) SCIPsnprintf(consname, SCIP_MAXSTRLEN, "row_%d", i);
    169 
    170         // check that all data is in valid ranges
    171         checkInfinity(scip, lhss[i], infinity);
    172         checkInfinity(scip, rhss[i], infinity);
    173 
    174         // create an empty linear constraint
    175         SCIP_CALL_ABORT( SCIPcreateConsLinear(scip, &cons, consname, 0, NULL, NULL, lhss[i], rhss[i],
    176               TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
    177              
    178         // add variables to constraint
    179         for(int j=rowstarts[i]; j<rowstarts[i]+rowlengths[i]; j++)       
    180         {
    181            checkInfinity(scip, coeffs[j], infinity);
    182            SCIP_CALL_ABORT( SCIPaddCoefLinear(scip, cons, vars[indices[j]], coeffs[j]) );
    183         }
    184 
    185         // add constraint to SCIP
    186         SCIP_CALL_ABORT( SCIPaddCons(scip, cons) );
    187         SCIP_CALL_ABORT( SCIPreleaseCons(scip, &cons) );       
    188      }
    189  
    190      // determine the method to solve the MILP
    191      if (milpMethod_ == 0 && niter == 0)
    192      {
    193         // initialize currentmilpmethod: at the root node we solve the MILP with SCIP default
    194         // deeper in the tree, we will solve the MILP by a FP heuristic
    195         if( depth == 0 )
    196            currentmilpmethod = 2;
    197         else
    198            currentmilpmethod = 4;
    199      }
    200      else if (milpMethod_ != 0)
    201         currentmilpmethod = milpMethod_; // use a fixed method to solve the MILP
    202      
    203      
    204      // MILP solving loop. If the MILP terminates without a solution, it might get resolved with a more expensive atrategy
    205      do {
    206         solveagain = false;
    207        
    208         // reset parameters if MILP is solved agian
    209         SCIP_CALL( SCIPresetParams(scip) );
    210        
    211         // deactivate SCIP output
    212         if (!(problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC))) {
    213            SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 0) );
    214         }
    215        
    216         // do not abort subproblem on CTRL-C
    217         SCIP_CALL( SCIPsetBoolParam(scip, "misc/catchctrlc", FALSE) );
    218        
    219         // set time limit
    220         timelimit = problem_ -> getMaxCpuTime () - CoinCpuTime ();
    221         SCIP_CALL( SCIPsetRealParam(scip, "limits/time", timelimit) );       
    222 
    223 
    224         if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC)) {
    225            SCIPinfoMessage(scip, NULL, "using MILP method: %d\n",currentmilpmethod);
    226         }
    227        
    228         // tune SCIP differently, depending on the chosen method to solve the MILP
    229         /// -1. Solve the MILP relaxation to proven optimality
    230         ///  0. Let Couenne choose
    231         ///  1. Partially solve the MILP with emphasis on good solutions
    232         ///  2. Solve the MILP relaxation partially, up to a certain node limit
    233         ///  3. Apply RENS to 1
    234         ///  4. Use Objective FP 2.0 for MILPs
    235         switch(currentmilpmethod)
    236         {
    237         case -1: // solve the MILP completely. SCIP's default setting should be best for this
    238            if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
    239            {
    240               SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
    241            }
    242            break;
    243 
    244         case 1: // Be aggressive in finding feasible solutions, but lazy about the dual bound.
    245            // Set limits on overall nodes and stall nodes (nodes without incumbent improvement).
    246            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/stallnodes", 1000) );
    247            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 10000) );
    248 
    249            // disable cutting plane separation
    250            SCIP_CALL_ABORT( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
    251        
    252            // disable expensive presolving
    253            SCIP_CALL_ABORT( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
    254 
    255            // use aggressive primal heuristics
    256            SCIP_CALL_ABORT( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_AGGRESSIVE, TRUE) );
    257        
    258            // use best estimate node selection
    259            if( SCIPfindNodesel(scip, "estimate") != NULL )
    260            {
    261               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "nodeselection/estimate/stdpriority", INT_MAX/4) );
    262            }
    263        
    264            // use inference branching
    265            if( SCIPfindBranchrule(scip, "inference") != NULL )
    266            {
    267               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
    268            }
    269        
    270            // disable conflict analysis
    271            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/useprop", FALSE) );
    272            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/useinflp", FALSE) );
    273            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/useboundlp", FALSE) );
    274            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/usesb", FALSE) );
    275            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/usepseudo", FALSE) );
    276 
    277            break;
    278 
    279         case 2: // default SCIP with node limits
    280            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/stallnodes", 500) );
    281            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 5000) );
    282            break;
    283         case 3: // solve the MILP with RENS. Disable most other features, enable RENS
    284            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 1) );
    285 
    286            // disable cutting plane separation
    287            SCIP_CALL_ABORT( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
    288        
    289            // disable expensive presolving
    290            SCIP_CALL_ABORT( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
    291 
    292            // besides RENS, only use cheap heuristics
    293            SCIP_CALL_ABORT( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
    294 
    295            // use inference branching
    296            if( SCIPfindBranchrule(scip, "inference") != NULL )
    297            {
    298               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
    299            }
    300 
    301            // ensure that RENS is called
    302            if( SCIPfindHeur(scip, "rens") != NULL )
    303            {
    304               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "heuristics/rens/freq", 0) );
    305               SCIP_CALL_ABORT( SCIPsetRealParam(scip, "heuristics/rens/minfixingrate", 0.0) );
    306            }
    307            if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
    308            {
    309               SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
    310            }
    311            break;
    312 
    313         case 4: // solve the MILP with Feasibility Pump. Disable most other features, enable stage 3 for feaspump
    314            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 1) );
    315 
    316            // disable cutting plane separation
    317            SCIP_CALL_ABORT( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
    318        
    319            // disable expensive presolving
    320            SCIP_CALL_ABORT( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
    321 
    322            // besides feaspump, only use cheap heuristics
    323            SCIP_CALL_ABORT( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
    324 
    325            // use inference branching
    326            if( SCIPfindBranchrule(scip, "inference") != NULL )
    327            {
    328               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
    329            }
    330 
    331            // ensure that feasibility pump is called
    332            if( SCIPfindHeur(scip, "feaspump") != NULL )
    333            {
    334               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "heuristics/feaspump/freq", 0) );
    335               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "heuristics/feaspump/maxsols", -1) );
    336 
    337               if( SCIPversion() <= 2.01 )
    338               {
    339                  SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "heuristics/feaspump2/stage3", TRUE) );
    340               }
    341               else
    342               {
    343                  SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "heuristics/feaspump/stage3", TRUE) );
    344               }
    345            }
    346            if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
    347            {
    348               SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
    349            }
    350            break;
    351 
    352         default:
    353            printf("invalid MILP method: %d\n", currentmilpmethod);
    354            assert(false);
    355            break;
    356         }
    357 
    358 #if 0
    359      // writes MILP problem and SCIP settings into a file
    360      SCIP_CALL_ABORT( SCIPwriteOrigProblem(scip, "debug.lp", NULL, FALSE) );
    361      SCIP_CALL_ABORT( SCIPwriteParams(scip, "debug.set", FALSE,TRUE) );
    362 #endif
    363          
    364         // solve the MILP
    365         SCIP_CALL_ABORT( SCIPsolve(scip) );
    366 
    367         nscipsols =  SCIPgetNSols(scip);
    368      
    369         // copy the solution
    370         if( nscipsols)
    371         {
    372            SCIP_SOL** scipsols;
    373            SCIP_SOL* bestsol;
    374            SCIP_Real cutoffbound;
    375 
    376            int nstoredsols;
    377 
    378            /* get incumbent solution */
    379            bestsol = SCIPgetBestSol(scip);
    380            assert(bestsol != NULL);
    381 
    382            /* get SCIP solution pool */
    383            scipsols = SCIPgetSols(scip);
    384            assert(scipsols != NULL);
    385 
    386            if (!sol)
    387              sol = new CouNumber [nvars];
    388 
    389            // get solution values and objective of incumbent
    390            SCIP_CALL_ABORT( SCIPgetSolVals(scip, bestsol, nvars, vars, sol) );
    391            obj = SCIPgetSolOrigObj(scip, bestsol);
    392 
    393            nstoredsols = 0;
    394 
    395            // if we use an objective feaspump, the objective function might be negative
    396            cutoffbound = obj > 0 ? 2*obj : obj / 2.0;
    397 
    398            // insert other SCIP solutions into solution pool
    399            // do not store too many or too poor solutions
    400            for(int i=1; i<nscipsols && nstoredsols < 10 &&
    401                   SCIPgetSolOrigObj(scip,scipsols[i]) <= cutoffbound; i++){
    402               double* tmpsol;
    403 
    404               tmpsol = new CouNumber [nvars];
    405            
    406               // get solution values
    407               SCIP_CALL_ABORT( SCIPgetSolVals(scip, scipsols[i], nvars, vars, tmpsol) );
    408               CouenneFPsolution couennesol = CouenneFPsolution (problem_, tmpsol);
    409 
    410               // add solutions to the pool if they are not in the tabu list
    411               if (tabuPool_. find (couennesol) == tabuPool_ . end ()
    412                   && pool_ -> Set (). find (couennesol) == pool_ -> Set(). end ()
    413                  ){
    414                  pool_ -> Set (). insert (couennesol);
    415 
    416                  ++nstoredsols;
    417               }
    418            }
    419 
    420            ++(*nsuciter);
    421 
    422            // if we succeeded five times in a row, try a cheaper MILP_ solving method next time
    423            // TODO: if we want to use time limits, hitting the time limit would be another good reason to switch
    424            if( *nsuciter >= 3 && currentmilpmethod < 4 )
    425            {
    426               ++currentmilpmethod;
    427               *nsuciter = 0;
    428            }         
    429         }
    430         //try to use a more aggressive, more expensive way to solve the MILP
    431         else if( milpMethod_ == 0 && currentmilpmethod > 1 )
    432         {
    433            --currentmilpmethod;
    434            solveagain = true;
    435            *nsuciter = 0;
    436 
    437            // throw away the current solution process
    438            SCIP_CALL( SCIPfreeTransform(scip) );
    439         }
    440         else
    441            obj = COIN_DBL_MAX; 
    442 
    443      } while (solveagain);
    444    
    445      // release variables before freeing them
    446      for (int i=0; i<nvars; i++) {
    447         SCIP_CALL_ABORT( SCIPreleaseVar(scip, &vars[i]) );
    448      }
    449 
    450      // free memory
    451      SCIPfreeMemoryArray(scip, &vars);
    452      SCIP_CALL_ABORT( SCIPfree(&scip) );
    453    
    454      BMScheckEmptyMemory();     
    455   }
    456   else
    45774#endif     
    45875  {
     
    48097
    48198/// initialize MILP solvers if needed
    482 void CouenneFeasPump::init_MILP () {
    483 
    484 }
     99void CouenneFeasPump::init_MILP () {}
  • trunk/Couenne/src/heuristics/CouenneFPSolveMILP.cpp

    r720 r883  
    159159
    160160    for (int i = 0; i < problem_ -> nVars (); ++i) {
     161
     162      if (problem_ -> Var (i) -> Multiplicity () <= 0)
     163        continue;
    161164
    162165      if (problem_ -> Var (i) -> isInteger () &&
     
    229232
    230233        for (int i = problem_ -> nVars (); i--;)
    231           if (milp_ -> isInteger (i))
     234
     235          if ((problem_ -> Var (i) -> Multiplicity () > 0) &&
     236              (milp_ -> isInteger (i)))
     237
    232238            newLB [i] = newUB [i] = iSol [i];
    233239
  • trunk/Couenne/src/heuristics/CouenneFPSolveNLP.cpp

    r720 r883  
    7373  nlp_     -> setObjective (newObj);
    7474
    75   if (problem_ -> Jnlst () -> ProduceOutput (J_STRONGWARNING, J_NLPHEURISTIC)) {
     75  if (problem_ -> Jnlst () -> ProduceOutput (J_ALL, J_NLPHEURISTIC)) {
    7676    printf ("now solving NLP:\n");
    7777    problem_ -> print ();
  • trunk/Couenne/src/heuristics/CouenneFPcreateMILP.cpp

    r733 r883  
    1818#include "CouenneProblem.hpp"
    1919#include "CouenneProblemElem.hpp"
     20#include "CouenneExprVar.hpp"
    2021
    2122#define COUENNE_EIG_RATIO .1 // how much smaller than the largest eigenvalue should the minimum be set at?
     
    5152        (!isMILP && !intVar))
    5253      // (empty) coeff col vector, lb = 0, ub = inf, obj coeff
    53       lp -> addCol (vec, 0., COIN_DBL_MAX, 1.);
     54      lp -> addCol (vec, 0., (fp -> Problem () -> Var (j) -> Multiplicity () <= 0) ? 0. : COIN_DBL_MAX, 1.);
    5455  }
    5556
     
    130131
    131132    for (int i=0; i<n; i++)
    132       P[i].insert (i, 1.);
     133      if (fp -> Problem () -> Var (i) -> Multiplicity () > 0)
     134        P[i].insert (i, 1.);
    133135  }
    134136
    135137  // Add 2q inequalities
    136138
    137   for (int i = 0, j = n, k = j; k--; ++i) {
     139  for (int i = 0, j = n, k = n; k--; ++i) {
     140
     141    if (fp -> Problem () -> Var (i) -> Multiplicity () <= 0)
     142      continue;
    138143
    139144    // two rows have to be added if:
     
    233238  // Add distance part
    234239  for (int i=0; i<n; ++i)
    235     A [i * (n+1)] += fp -> multDistMILP ();
     240    if (fp -> Problem () -> Var (i) -> Multiplicity () > 0)
     241      A [i * (n+1)] += fp -> multDistMILP ();
    236242
    237243  // Add gradient-parallel term to the hessian, (x_z - x_z_0)^2. This
     
    331337    }
    332338
    333   if (fp -> Problem () -> Jnlst () -> ProduceOutput (Ipopt::J_STRONGWARNING, J_NLPHEURISTIC)) {
    334 
    335     printf ("P:\n");
    336 
    337 
    338     printf ("P^{1/2}:\n");
    339 
    340   }
     339  // if (fp -> Problem () -> Jnlst () -> ProduceOutput (Ipopt::J_STRONGWARNING, J_NLPHEURISTIC)) {
     340  //   printf ("P:\n");
     341  //   printf ("P^{1/2}:\n");
     342  // }
    341343
    342344  free (eigenval);
  • trunk/Couenne/src/heuristics/CouenneFPpool.cpp

    r698 r883  
    4949       ++i) {
    5050
     51    if ((*i) -> Multiplicity () <= 0)
     52      continue;
     53
    5154    CouNumber
    5255      vval = (**i) ();
    53 
    54     if ((*i) -> Multiplicity () <= 0)
    55       continue;
    5656
    5757    if ((*i) -> isInteger ()) {
     
    149149  switch (comparedTerm) {
    150150
    151   case SUM_NINF: return (nNLinf_   + nIinf_   < other.nNLinf_   + other.nIinf_);
    152   case SUM_INF:  return (maxNLinf_ + maxIinf_ < other.maxNLinf_ + other.maxIinf_);
    153   case OBJVAL:   return (objVal_              < other.objVal_);
     151  case SUM_NINF:     return (nNLinf_   + nIinf_   < other.nNLinf_   + other.nIinf_);
     152  case SUM_INF:      return (maxNLinf_ + maxIinf_ < other.maxNLinf_ + other.maxIinf_);
     153  case OBJVAL:       return (objVal_              < other.objVal_);
     154  case INTEGER_VARS: {
     155
     156    // lexicographical comparison
     157
     158    for (int i=0; i<n_; ++i)
     159      if (x_ [i] < other.x_ [i])
     160        return true;
     161    return false;
     162  }
    154163  }
    155164
     
    217226          delta = *x++ - *s++;
    218227
     228          // TODO: check multiplicity () > 0)
     229
    219230          dist += delta * delta;
    220231
     
    229240          continue;
    230241
    231          //update best solution
    232          if( dist < bestdist )
     242        //update best solution
     243        if( dist < bestdist )
    233244         {
    234245            bestdist = dist;
  • trunk/Couenne/src/heuristics/CouenneFPpool.hpp

    r691 r883  
    2626  /// what term to compare: the sum of infeasibilities, the sum of
    2727  /// numbers of infeasible terms, or the objective function
    28   static enum what_to_compare {SUM_NINF = 0, SUM_INF, OBJVAL} comparedTerm_;
     28  static enum what_to_compare {SUM_NINF = 0, SUM_INF, OBJVAL, INTEGER_VARS} comparedTerm_;
    2929
    3030  /// Class containing a solution with infeasibility evaluation
     
    3434
    3535    CouNumber *x_;        ///< solution
    36     int        n_;        ///< number of variables (for independence from CouenneProblem
     36    int        n_;        ///< number of variables (for independence from CouenneProblem)
    3737    int        nNLinf_;   ///< number of NL      infeasibilities
    3838    int        nIinf_;    ///< number of integer infeasibilities
  • trunk/Couenne/src/heuristics/CouenneFeasPump.cpp

    r733 r883  
    3333using namespace Couenne;
    3434
    35 void printDist (CouenneProblem *p, double *iSol, double *nSol);
    36 void printCmpSol (int n, double *iSol, double *nSol, int direction);
     35void printDist   (CouenneProblem *p, double *iSol, double *nSol);
     36void printCmpSol (CouenneProblem *p, double *iSol, double *nSol, int direction);
    3737
    3838// Solve
     
    7979
    8080  if (problem_ -> Jnlst () -> ProduceOutput 
    81       (J_WARNING, J_NLPHEURISTIC)) {
     81      (J_ALL, J_NLPHEURISTIC)) {
    8282
    8383    printf ("initial NLP:\n");
     
    175175
    176176    if (!iSol || z >= COIN_DBL_MAX/2) {
    177 
    178177      problem_ -> Jnlst () -> Printf
    179178        (Ipopt::J_ERROR, J_NLPHEURISTIC,
    180          "FP: breaking out of loop upon %p==NULL or %.3f large\n", iSol, z);
    181 
     179         "FP: breaking out of loop upon %p==NULL or z large (z=%.3f)\n", iSol, z);
    182180      break;
    183181    }
     
    256254        int startIndex = (int) floor (CoinDrand48 () * problem_ -> nOrigVars ());
    257255
    258         for (int ii=problem_ -> nOrigVars (); ii--; lb++, ub++) {
    259 
    260           // make perturbation start from random points
     256        for (int ii = problem_ -> nOrigVars (); ii--; lb++, ub++) {
     257
     258          if (problem_ -> Var (ii) -> Multiplicity () <= 0)
     259            continue;
     260
     261          // make perturbation start from pseudo-random points
    261262
    262263          int i = (startIndex + ii) % problem_ -> nOrigVars ();
     
    283284        }
    284285      }
    285 
    286286    } else tabuPool_. insert (CouenneFPsolution (problem_, iSol));
    287287
     
    291291                                       true, // checkALL
    292292                                       problem_->getFeasTol());
    293     if(isChecked) {
     293    if (isChecked) {
    294294      z = problem_->getRecordBestSol()->getModSolVal();
    295295    }
     
    299299#endif  /* not FM_CHECKNLP2 */
    300300
     301    // Possible improvement: if IP-infeasible or if z does not improve
     302    // the best solution found so far, then try the following:
     303    //
     304    // 1) Fix integer variables to IP solution's corresponding components
     305    // 2) Solve restriction with original obj
     306    // 3) If still MINLP infeasible or z not better
     307    // 4)   Repeat
     308    // 5)     Get solution x* from pool
     309    // 6)     Solve with original object
     310    // 7)   Until MINLP feasible and z better or no more solutions in the pool
     311    // 8) If found solution, set it, otherwise keep previously found (IP-feasible) one
     312
     313    if (!isChecked || z > problem_ -> getCutOff ()) {
     314
     315     
     316
     317    }
     318
     319    // Whatever the previous block yielded, we are now going to check
     320    // if we have a new solution, and if so we save it.
     321
    301322    if (isChecked) {
    302323
    303       problem_ -> Jnlst () -> Printf
    304         (Ipopt::J_ERROR, J_NLPHEURISTIC,
    305          "FP: IP solution is MINLP feasible\n");
     324      problem_ -> Jnlst () -> Printf (Ipopt::J_ERROR, J_NLPHEURISTIC, "FP: IP solution is MINLP feasible\n");
    306325
    307326      // solution is MINLP feasible! Save it.
     
    310329      objVal = z;
    311330
    312 #ifdef FM_CHECKNLP2
    313 #ifdef FM_TRACE_OPTSOL
    314       problem_->getRecordBestSol()->update();
    315       best = problem_->getRecordBestSol()->getSol();
    316       objVal = problem_->getRecordBestSol()->getVal();
    317 #else /* not FM_TRACE_OPTSOL */
    318       best = problem_->getRecordBestSol()->getModSol(problem_->nVars());
    319       objVal = z;
    320 #endif /* not FM_TRACE_OPTSOL */
    321 #else /* not FM_CHECKNLP2 */
    322 #ifdef FM_TRACE_OPTSOL
    323       problem_->getRecordBestSol()->update(iSol, problem_->nVars(),
    324                                            z, problem_->getFeasTol());
    325       best = problem_->getRecordBestSol()->getSol();
    326       objVal = problem_->getRecordBestSol()->getVal();
    327 #else /* not FM_TRACE_OPTSOL */
     331      // Found a MINLP-feasible solution, but to keep diversity do not
     332      // use the best available. Just use this.
     333      //
     334      // #ifdef FM_CHECKNLP2
     335      // #  ifdef FM_TRACE_OPTSOL
     336      //       problem_->getRecordBestSol()->update();
     337      //       best = problem_->getRecordBestSol()->getSol();
     338      //       objVal = problem_->getRecordBestSol()->getVal();
     339      // #  else /* not FM_TRACE_OPTSOL */
     340      //       best = problem_->getRecordBestSol()->getModSol(problem_->nVars());
     341      //       objVal = z;
     342      // #  endif /* not FM_TRACE_OPTSOL */
     343      // #else /* not FM_CHECKNLP2 */
     344      // #  ifdef FM_TRACE_OPTSOL
     345      //       problem_->getRecordBestSol()->update(iSol, problem_->nVars(),
     346      //                                           z, problem_->getFeasTol());
     347      //       best = problem_->getRecordBestSol()->getSol();
     348      //       objVal = problem_->getRecordBestSol()->getVal();
     349      // #  else /* not FM_TRACE_OPTSOL */
     350
    328351      best   = iSol;
    329352      objVal = z;
    330 #endif /* not FM_TRACE_OPTSOL */
    331 #endif /* not FM_CHECKNLP2 */
     353
     354      // #  ifdef FM_TRACE_OPTSOL
     355      //       problem_->getRecordBestSol()->update();
     356      //       best = problem_->getRecordBestSol()->getSol();
     357      //       objVal = problem_->getRecordBestSol()->getVal();
     358
     359
     360      // #  endif /* not FM_TRACE_OPTSOL */
     361      // #endif /* not FM_CHECKNLP2 */
    332362
    333363      if (z < problem_ -> getCutOff ()) {
     
    338368
    339369        if (objInd >= 0) {
    340        
    341370          chg_bds = new t_chg_bounds [problem_ -> nVars ()];
    342371          chg_bds [objInd].setUpper (t_chg_bounds::CHANGED);
     
    346375        bool is_still_feas = problem_ -> btCore (chg_bds);
    347376
     377        if (chg_bds)
     378          delete [] chg_bds;
     379
    348380        // don't tighten MILP if BT says it's infeasible
    349 
    350         if (chg_bds)
    351           delete [] chg_bds;
    352381
    353382        if (!is_still_feas)
     
    362391
    363392        for (int i=problem_ -> nVars (), j=0; i--; ++j, ++plb, ++pub) {
    364            
    365           if (*plb > *mlb++) milp_ -> setColLower (j, *plb);
    366           if (*pub < *mub++) milp_ -> setColUpper (j, *pub);
     393
     394          bool neglect = problem_ -> Var (j) -> Multiplicity () <= 0;
     395
     396          if (*plb > *mlb++) milp_ -> setColLower (j, neglect ? 0. : *plb);
     397          if (*pub < *mub++) milp_ -> setColUpper (j, neglect ? 0. : *pub);
    367398        }
    368399      }
     
    373404               milpCuttingPlane_ == FP_CUT_POST) {
    374405
    375       // solution is not MINLP feasible, it might get cut by
     406      // Solution is IP- but not MINLP feasible: it might get cut by
    376407      // linearization cuts. If so, add a round of cuts and repeat.
    377408
     
    416447      for (int i = 0; i < problem_ -> nVars (); ++i) {
    417448
    418         if (problem_ -> Var (i) -> isInteger () &&
     449        exprVar *e = problem_ -> Var (i);
     450
     451        if (e -> Multiplicity () <= 0)
     452          continue;
     453
     454        if (e -> isInteger () &&
    419455            (fabs (iSol [i] - ceil (iSol [i] - .5)) > 1e-4))
    420456          ++nNonint;
     
    431467
    432468      printDist   (problem_,             iSol, nSol);
    433       printCmpSol (problem_ -> nVars (), iSol, nSol, 0);
     469      printCmpSol (problem_, iSol, nSol, 0);
    434470    }
    435471
     
    571607
    572608#else /* not FM_CHECKNLP2 */
    573     isChecked = problem_ -> checkNLP (nSol, z, true);
     609      isChecked = problem_ -> checkNLP (nSol, z, true);
    574610#endif  /* not FM_CHECKNLP2 */
    575611    }
     
    605641
    606642  if (retval > 0) {
    607 
    608643    if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_NLPHEURISTIC)) {
    609 
    610644      printf ("FP: returning MINLP feasible solution:\n");
    611645      printDist (problem_, best, nSol);
    612646    }
    613 
    614647    CoinCopyN (best, problem_ -> nVars (), newSolution);
    615648  }
     
    661694       ++i) {
    662695
     696    if ((*i) -> Multiplicity () <= 0)
     697      continue;
     698
    663699    CouNumber
    664700      vval = (**i) ();
    665 
    666     if ((*i) -> Multiplicity () <= 0)
    667       continue;
    668701
    669702    if ((*i) -> isInteger ()) {
     
    716749    dist = 0.;
    717750
    718     for (int i = p -> nVars (); i--;)
     751    for (int i = p -> nVars (); i--;) {
     752
     753      if (p -> Var (i) -> Multiplicity () <= 0)
     754        continue;
     755
    719756      dist +=
    720757        (iSol [i] - nSol [i]) *
    721758        (iSol [i] - nSol [i]);
     759    }
    722760
    723761    dist = sqrt (dist);
     
    736774#define WRAP 3
    737775
    738 void printCmpSol (int n, double *iSol, double *nSol, int direction) {
    739  
     776void printCmpSol (CouenneProblem *p, double *iSol, double *nSol, int direction) {
     777
     778  int n = p -> nVars ();
     779
    740780  printf ("i:%p n:%p\nFP # ",
    741781          (void *) iSol, (void *) nSol);
     
    751791  for (int i=0; i<n; i++) {
    752792
     793    if (p -> Var (i) -> Multiplicity () <= 0)
     794      continue;
     795
    753796    if (i && !(i % WRAP))
    754797      printf ("\nFP # ");
    755798
    756799    double
    757       iS = iSol ? iSol [i] : 12345.,
     800      iS = iSol ? iSol [i] : 12345., // flag values, to be seen if something is wrong
    758801      nS = nSol ? nSol [i] : 54321.;
    759802
  • trunk/Couenne/src/heuristics/CouenneFeasPump.hpp

    r713 r883  
    2424#include "IpOptionsList.hpp"
    2525
     26#ifdef COIN_HAS_SCIP
     27#include "scip/scip.h"
     28#endif
     29
    2630struct Scip;
    2731class OsiSolverInterface;
     
    9599    /// (according to the l-1 norm of the hessian) to the current
    96100    /// NLP-feasible (but fractional) solution nsol
    97     CouNumber solveMILP (CouNumber *nSol, CouNumber *&iSol, int niter, int* nsuciter);
     101    virtual CouNumber solveMILP (CouNumber *nSol, CouNumber *&iSol, int niter, int* nsuciter);
    98102
    99103    /// obtain solution to NLP
    100     CouNumber solveNLP  (CouNumber *nSol, CouNumber *&iSol);
     104    virtual CouNumber solveNLP  (CouNumber *nSol, CouNumber *&iSol);
    101105
    102106    /// set new expression as the NLP objective function using
     
    145149    {return nlp_;}
    146150
     151#ifdef COIN_HAS_SCIP
     152    SCIP_RETCODE ScipSolve (double* &sol, int niter, int* nsuciter);
     153#endif
     154
    147155  private:
    148156
  • trunk/Couenne/src/heuristics/CouenneFeasPumpConstructors.cpp

    r720 r883  
    4747  app_ -> Options () -> SetStringValue ("fixed_variable_treatment", "make_parameter");
    4848
     49  // Suppress iteration output from nonlinear solver
     50  app_ -> Options () -> SetStringValue ("sb", "yes", false, true);
     51
    4952  if (status != Solve_Succeeded)
    5053    printf ("FP: Error in initialization\n");
     
    141144
    142145  } else
    143     pool_ = new CouenneFPpool (SUM_NINF);
     146    //pool_ = new CouenneFPpool (SUM_NINF);
     147    pool_ = new CouenneFPpool (INTEGER_VARS);
    144148
    145149  setHeuristicName ("Couenne Feasibility Pump");
     
    281285    for (int i=0; i<problem_ -> nVars (); i++) {
    282286
     287      if (problem_ -> Var (i) -> Multiplicity () <= 0)
     288        continue;
     289
    283290      if (compDistInt_ == FP_DIST_INT &&
    284291          !(problem_ -> Var (i) -> isInteger ()))
     
    329336    for (int i=0; i<num; ++i, ++row, ++col, ++val) {
    330337
     338      if ((problem_ -> Var (*row) -> Multiplicity () <= 0) ||
     339          (problem_ -> Var (*col) -> Multiplicity () <= 0))
     340        continue;
     341
    331342      // check if necessary given options
    332343
     
    380391      // create the argument list (x_i - x_i^0)^2 for all i's
    381392      for (int i=0; i<problem_ -> nVars (); i++) {
     393
     394        if (problem_ -> Var (i) -> Multiplicity () <= 0)
     395          continue;
    382396         
    383397        if ((compDistInt_ == FP_DIST_INT &&
     
    425439  for (int i = problem_ -> nVars (); i--;)
    426440
    427     if (problem_ -> Var (i) -> isInteger ()) {
     441    if ((problem_ -> Var (i) -> isInteger ()) &&
     442        (problem_ -> Var (i) -> Multiplicity () > 0)) {
    428443
    429444      double
  • trunk/Couenne/src/heuristics/Makefile.am

    r765 r883  
    2020        CouenneFPFindSolution.cpp \
    2121        CouenneFPpool.cpp \
    22         CouenneIterativeRounding.cpp
     22        CouenneFPscipSolve.cpp \
     23        CouenneIterativeRounding.cpp
    2324
    2425if COIN_HAS_SCIP
  • trunk/Couenne/src/heuristics/Makefile.in

    r765 r883  
    6363        CouenneFPSolveMILP.cpp CouenneFPSolveNLP.cpp \
    6464        CouenneFPFindSolution.cpp CouenneFPpool.cpp \
    65         CouenneIterativeRounding.cpp cons_rowcuts.cpp
     65        CouenneFPscipSolve.cpp CouenneIterativeRounding.cpp \
     66        cons_rowcuts.cpp
    6667@COIN_HAS_SCIP_TRUE@am__objects_1 = cons_rowcuts.lo
    6768am_libCouenneHeuristics_la_OBJECTS = BonInitHeuristic.lo \
     
    7071        CouenneFPSolveMILP.lo CouenneFPSolveNLP.lo \
    7172        CouenneFPFindSolution.lo CouenneFPpool.lo \
    72         CouenneIterativeRounding.lo $(am__objects_1)
     73        CouenneFPscipSolve.lo CouenneIterativeRounding.lo \
     74        $(am__objects_1)
    7375libCouenneHeuristics_la_OBJECTS =  \
    7476        $(am_libCouenneHeuristics_la_OBJECTS)
     
    324326        CouenneFPSolveMILP.cpp CouenneFPSolveNLP.cpp \
    325327        CouenneFPFindSolution.cpp CouenneFPpool.cpp \
    326         CouenneIterativeRounding.cpp $(am__append_1)
     328        CouenneFPscipSolve.cpp CouenneIterativeRounding.cpp \
     329        $(am__append_1)
    327330
    328331# This is for libtool
     
    400403@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CouenneFPcreateMILP.Plo@am__quote@
    401404@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CouenneFPpool.Plo@am__quote@
     405@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CouenneFPscipSolve.Plo@am__quote@
    402406@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CouenneFeasPump.Plo@am__quote@
    403407@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CouenneFeasPumpConstructors.Plo@am__quote@
Note: See TracChangeset for help on using the changeset viewer.