Changeset 684
 Timestamp:
 Jul 11, 2007 5:05:34 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel/Cbc/examples/link2.cpp
r681 r684 14 14 #include "CbcModel.hpp" 15 15 #include "CbcBranchActual.hpp" 16 #include "CbcCutGenerator.hpp" 16 17 #include "OsiChooseVariable.hpp" 17 18 #include "CoinModel.hpp" … … 29 30 constrained problem. 30 31 31 Having variables integer would make easier32 Having variables integer would probably make easier 32 33 ************************************************************************/ 33 34 … … 39 40 CoinModel build; 40 41 // size 41 int numberX= 5;42 int numberConstraints= 20;42 int numberX=20; 43 int numberConstraints=1000; 43 44 // variables are x sub j and w sub ij which is x sub j * x sub i 44 45 char name[20]; … … 47 48 // pointers to W (1 if not existing) 48 49 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++) { 58 51 int * sequenceW = new int[numberX]; 59 52 int j; 60 for (j=0;j< i;j++)53 for (j=0;j<numberX;j++) 61 54 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 squares67 double value = (i==j) ? 1.0 : 0.0;68 build.addColumn(0,NULL,NULL,0.0,1.0,value,name);69 numberW++;70 }71 55 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==numberXjColumn==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 } 72 113 } 73 114 int numberColumns = build.numberColumns(); 74 115 int * column = new int [numberColumns]; 75 116 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? 77 118 for (i=0;i<numberX;i++) { 78 119 column[i]=i; 79 120 element[i]=1.0; 80 121 } 122 // Should be stronger if we make sum x == 1 rather >= 1 123 #define EQUALITY 124 #ifndef EQUALITY 81 125 build.addRow(numberX,column,element,1.0,COIN_DBL_MAX); 126 #else 127 build.addRow(numberX,column,element,1.0,1.0); 128 #endif 82 129 // Now define W 83 130 for (i=0;i<numberX;i++) { … … 98 145 } 99 146 #if 1 100 // Optional linear constraints 147 /* Optional linear constraints 148 These are theoretically redundant but may tighten relaxation 149 */ 101 150 for (i=0;i<numberX;i++) { 102 151 // We want to add in sum over j W_i_j >= X_i … … 117 166 column[n]=i; 118 167 element[n++]=1.0; 168 #ifndef EQUALITY 119 169 build.addRow(n,column,element,0.0,COIN_DBL_MAX); 170 #else 171 build.addRow(n,column,element,0.0,0.0); 172 #endif 120 173 } 121 174 #endif 122 175 // 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! 124 181 #define USE_CUTS 125 182 #ifndef USE_CUTS … … 136 193 } 137 194 } 138 build.addRow(n,column,element,COIN_DBL_MAX,scaleFactor); 195 // make sure last one feasible at 1.0 196 if (element[n1]>rhs) 197 element[n1]=0.9*rhs; 198 build.addRow(n,column,element,COIN_DBL_MAX,rhs); 139 199 } 140 200 #else … … 155 215 } 156 216 } 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[n1]>rhs) 219 element[n1]=0.9*rhs; 220 storedCuts.addCut(COIN_DBL_MAX,rhs,n,column,element); 221 build2.addRow(n,column,element,COIN_DBL_MAX,rhs); 159 222 } 160 223 // load from coin model … … 174 237 // load from coin model 175 238 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);183 239 solver1.setDefaultMeshSize(0.001); 184 // need some relative granularity240 // Use coarse grid so we can solve quickly 185 241 solver1.setDefaultMeshSize(0.01); 242 solver1.setDefaultMeshSize(0.001); 243 //solver1.setDefaultMeshSize(0.0001); 186 244 solver1.setIntegerPriority(1000); 187 245 solver1.setBiLinearPriority(10000); 188 246 solver1.load(build,true,2); 247 // Add more objects with even coarser grid 248 solver1.setBiLinearPriorities(10,0.2); 189 249 CbcModel model(solver1); 190 250 // set more stuff … … 208 268 model.setLogLevel(1); 209 269 #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); 213 276 model.branchAndBound(); 214 277 … … 219 282 <<std::endl; 220 283 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 } 221 296 222 297 if (model.getMinimizationObjValue()<1.0e50) { 223 298 224 299 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[inumberX]; 308 int jX = wJ[inumberX]; 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[iColumnnumberX]; 326 int jX = wJ[iColumnnumberX]; 327 sum += solution[iX]*solution[jX]*element[j]; 328 } 329 if (sum>cut>ub()+1.0e6) 330 printf("*** %d %g (approx %g) > %g\n",i,sum,sumApprox,cut>ub()); 331 else if (sum>cut>ub()1.0e3) 332 printf(" %d %g (approx %g) =~ %g\n",i,sum,sumApprox,cut>ub()); 333 if (sum>cut>ub()) { 334 double violation = sumcut>ub(); 335 sumViolation += violation; 336 worstViolation = CoinMax(worstViolation,violation); 337 } 338 } 339 printf("Total violation %g, worst was %g\n",sumViolation,worstViolation); 340 #endif 227 341 } 228 342 // deletes … … 231 345 } 232 346 delete [] whichW; 347 delete [] wI; 348 delete [] wJ; 233 349 return 0; 234 350 }
Note: See TracChangeset
for help on using the changeset viewer.