Changeset 684 for branches/devel/Cbc


Ignore:
Timestamp:
Jul 11, 2007 5:05:34 PM (12 years ago)
Author:
forrest
Message:

updated example

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Cbc/examples/link2.cpp

    r681 r684  
    1414#include "CbcModel.hpp"
    1515#include "CbcBranchActual.hpp"
     16#include "CbcCutGenerator.hpp"
    1617#include "OsiChooseVariable.hpp"
    1718#include "CoinModel.hpp"
     
    2930constrained problem.
    3031
    31 Having variables integer would make easier
     32Having variables integer would probably make easier
    3233************************************************************************/
    3334
     
    3940  CoinModel build;
    4041  // size
    41   int numberX=5;
    42   int numberConstraints=20;
     42  int numberX=20;
     43  int numberConstraints=1000;
    4344  // variables are x sub j and w sub ij which is x sub j * x sub i
    4445  char name[20];
     
    4748  // pointers to W (-1 if not existing)
    4849  int ** whichW = new int * [numberX];
    49   int numberW=0;
    50   // build columns
    51   for (i=0;i<numberX;i++) {
    52     sprintf(name,"X_%d",i);
    53     build.addColumn(0,NULL,NULL,0.0,1.0,0.0,name);
    54   }
    55   // now W
    56   for (i=0;i<numberX;i++) {
    57     // and do pointers to W
     50  for (i=0;i<numberX;i++) {
    5851    int * sequenceW = new int[numberX];
    5952    int j;
    60     for (j=0;j<i;j++)
     53    for (j=0;j<numberX;j++)
    6154      sequenceW[j]=-1;
    62     for (;j<numberX;j++) {
    63       sequenceW[j]=numberX+numberW;
    64       assert (sequenceW[j]==build.numberColumns());
    65       sprintf(name,"W_%d_%d",i,j);
    66       // minimize sums of squares
    67       double value = (i==j) ? 1.0 : 0.0;
    68       build.addColumn(0,NULL,NULL,0.0,1.0,value,name);
    69       numberW++;
    70     }
    7155    whichW[i]=sequenceW;
     56  }
     57  /* now decide on W.
     58     When both x and y are continuous then the default method for x*y will branch on x
     59     until it is fixed and x*y can be modeled exactly.  In this case it is better not to
     60     use a triangular version for W but to randomize so each variable is the "x" one
     61     about the same number of times.
     62  */
     63  int neededW=(numberX*(numberX+1))/2;
     64  int * wI = new int[neededW];
     65  int * wJ = new int[neededW];
     66  while (neededW) {
     67    bool gotColumn=false;
     68    int iColumn=-1;
     69    int jColumn=-1;
     70    while(!gotColumn) {
     71      iColumn=(int) floor(CoinDrand48()*numberX);
     72      jColumn=(int) floor(CoinDrand48()*numberX);
     73      if (iColumn==numberX||jColumn==numberX)
     74        continue;
     75      int * sequenceW = whichW[iColumn];
     76      if (sequenceW[jColumn]==-1) {
     77        // available
     78        gotColumn=true;
     79        break;
     80      }
     81    }
     82    int * sequenceW = whichW[iColumn];
     83    sequenceW[jColumn]=1;
     84    neededW--;
     85    if (iColumn!=jColumn) {
     86      int * sequenceW = whichW[jColumn];
     87      sequenceW[iColumn]=-2;
     88    }
     89  }
     90  // build columns
     91  for (i=0;i<numberX;i++) {
     92    sprintf(name,"X_%d",i);
     93    // If we know tighter bounds then better
     94    build.addColumn(0,NULL,NULL,0.0,1.0,0.0,name);
     95  }
     96  // now do W
     97  int numberW=0;
     98  for (i=0;i<numberX;i++) {
     99    int * sequenceW = whichW[i];
     100    for (int j=0;j<numberX;j++) {
     101      if (sequenceW[j]>0) {
     102        sequenceW[j]=numberX+numberW;
     103        assert (sequenceW[j]==build.numberColumns());
     104        sprintf(name,"W_%d_%d",i,j);
     105        // minimize sums of squares of x so ..
     106        double value = (i==j) ? 1.0 : 0.0;
     107        build.addColumn(0,NULL,NULL,0.0,1.0,value,name);
     108        wI[numberW]=i;
     109        wJ[numberW]=j;
     110        numberW++;
     111      }
     112    }
    72113  }
    73114  int numberColumns = build.numberColumns();
    74115  int * column = new int [numberColumns];
    75116  double * element = new double [numberColumns];
    76   // first row says sum of x >= 1.0
     117  // first row says sum of x >= 1.0 *** no make equality?
    77118  for (i=0;i<numberX;i++) {
    78119    column[i]=i;
    79120    element[i]=1.0;
    80121  }
     122  // Should be stronger if we make sum x == 1 rather >= 1
     123#define EQUALITY
     124#ifndef EQUALITY
    81125  build.addRow(numberX,column,element,1.0,COIN_DBL_MAX);
     126#else
     127  build.addRow(numberX,column,element,1.0,1.0);
     128#endif
    82129  // Now define W
    83130  for (i=0;i<numberX;i++) {
     
    98145  }
    99146#if 1
    100   // Optional linear constraints
     147  /* Optional linear constraints
     148     These are theoretically redundant but may tighten relaxation
     149  */
    101150  for (i=0;i<numberX;i++) {
    102151    // We want to add in sum over j W_i_j >= X_i
     
    117166    column[n]=i;
    118167    element[n++]=-1.0;
     168#ifndef EQUALITY
    119169    build.addRow(n,column,element,0.0,COIN_DBL_MAX);
     170#else
     171    build.addRow(n,column,element,0.0,0.0);
     172#endif
    120173  }
    121174#endif
    122175  // Now add in random quadratics - no linear term but that would be easy (X variables)
    123   double scaleFactor=0.5;
     176  double rhs=0.5; // So constraints do something
     177  if (numberX==8)
     178    rhs=0.366;
     179  else if (numberX==20)
     180    rhs=0.25; // so will be quickish!
    124181#define USE_CUTS
    125182#ifndef USE_CUTS
     
    136193      }
    137194    }
    138     build.addRow(n,column,element,-COIN_DBL_MAX,scaleFactor);
     195    // make sure last one feasible at 1.0
     196    if (element[n-1]>rhs)
     197      element[n-1]=0.9*rhs;
     198    build.addRow(n,column,element,-COIN_DBL_MAX,rhs);
    139199  }
    140200#else
     
    155215      }
    156216    }
    157     storedCuts.addCut(-COIN_DBL_MAX,scaleFactor,n,column,element);
    158     build2.addRow(n,column,element,-COIN_DBL_MAX,scaleFactor);
     217    // make sure last one feasible at 1.0
     218    if (element[n-1]>rhs)
     219      element[n-1]=0.9*rhs;
     220    storedCuts.addCut(-COIN_DBL_MAX,rhs,n,column,element);
     221    build2.addRow(n,column,element,-COIN_DBL_MAX,rhs);
    159222  }
    160223  // load from coin model
     
    174237  // load from coin model
    175238  OsiSolverLink solver1;
    176   //OsiSolverInterface * solver2 = solver1.clone();
    177   ///CbcModel model;
    178   //model.assignSolver(solver2,true);
    179   //OsiSolverLink * si =
    180   //dynamic_cast<OsiSolverLink *>(model.solver()) ;
    181   //assert (si != NULL);
    182   //solver1.setSpecialOptions2(16);
    183239  solver1.setDefaultMeshSize(0.001);
    184   // need some relative granularity
     240  // Use coarse grid so we can solve quickly
    185241  solver1.setDefaultMeshSize(0.01);
     242  solver1.setDefaultMeshSize(0.001);
     243  //solver1.setDefaultMeshSize(0.0001);
    186244  solver1.setIntegerPriority(1000);
    187245  solver1.setBiLinearPriority(10000);
    188246  solver1.load(build,true,2);
     247  // Add more objects with even coarser grid
     248  solver1.setBiLinearPriorities(10,0.2);
    189249  CbcModel model(solver1);
    190250  // set more stuff
     
    208268  model.setLogLevel(1);
    209269#ifdef USE_CUTS
    210   model.addCutGenerator(&storedCuts,1,"Stored cuts",false,true);
    211   model.setMaximumCutPasses(4);
    212 #endif
     270  model.addCutGenerator(&storedCuts,1,"Stored cuts",true,false);
     271  model.cutGenerator(0)->setMustCallAgain(true);;
     272  //model.addCutGenerator(&storedCuts,1,"Stored cuts",false,true);
     273  //model.setMaximumCutPasses(4);
     274#endif
     275  //model.setNumberThreads(4);
    213276  model.branchAndBound();
    214277
     
    219282           <<std::endl;
    220283
     284  int numberGenerators = model.numberCutGenerators();
     285  for (int iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
     286    CbcCutGenerator * generator = model.cutGenerator(iGenerator);
     287    std::cout<<generator->cutGeneratorName()<<" was tried "
     288             <<generator->numberTimesEntered()<<" times and created "
     289             <<generator->numberCutsInTotal()<<" cuts of which "
     290             <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
     291    if (generator->timing())
     292      std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
     293    else
     294      std::cout<<std::endl;
     295  }
    221296
    222297  if (model.getMinimizationObjValue()<1.0e50) {
    223298   
    224299    const double * solution = model.bestSolution();
    225     for (i=0;i<numberColumns;i++)
    226       printf("%d %g\n",i,solution[i]);
     300    double obj=0.0;
     301    for (i=0;i<numberX;i++) {
     302      printf("%d %s %g\n",i,build.columnName(i),solution[i]);
     303      obj += solution[i]*solution[i];
     304    }
     305    printf("true objective %g\n",obj);
     306    for (;i<numberColumns;i++) {
     307      int iX = wI[i-numberX];
     308      int jX = wJ[i-numberX];
     309      printf("%d %s %g - should be %g\n",i,build.columnName(i),solution[i],solution[iX]*solution[jX]);
     310    }
     311#ifdef USE_CUTS
     312    double worstViolation=0.0;
     313    double sumViolation=0.0;
     314    for (i=0;i<numberConstraints;i++) {
     315      const OsiRowCut * cut = storedCuts.rowCutPointer(i);
     316      const CoinPackedVector row = cut->row();
     317      int n = row.getNumElements();
     318      const int * column = row.getIndices();
     319      const double * element = row.getElements();
     320      double sum=0.0;
     321      double sumApprox=0.0;
     322      for (int j=0;j<n;j++) {
     323        int iColumn = column[j];
     324        sumApprox += solution[iColumn]*element[j];
     325        int iX = wI[iColumn-numberX];
     326        int jX = wJ[iColumn-numberX];
     327        sum += solution[iX]*solution[jX]*element[j];
     328      }
     329      if (sum>cut->ub()+1.0e-6)
     330        printf("*** %d %g (approx %g) > %g\n",i,sum,sumApprox,cut->ub());
     331      else if (sum>cut->ub()-1.0e-3)
     332        printf("    %d %g (approx %g) =~ %g\n",i,sum,sumApprox,cut->ub());
     333      if (sum>cut->ub()) {
     334        double violation = sum-cut->ub();
     335        sumViolation += violation;
     336        worstViolation = CoinMax(worstViolation,violation);
     337      }
     338    }
     339    printf("Total violation %g, worst was %g\n",sumViolation,worstViolation);
     340#endif
    227341  }
    228342  // deletes
     
    231345  }
    232346  delete [] whichW;
     347  delete [] wI;
     348  delete [] wJ;
    233349  return 0;
    234350}   
Note: See TracChangeset for help on using the changeset viewer.