My Project
OSColGenApp.cpp
Go to the documentation of this file.
1/* $Id: OSColGenApp.cpp 3038 2009-11-07 11:43:44Z kmartin $ */
12#include "OSColGenApp.h"
13#include "OSErrorClass.h"
14#include "OSDataStructures.h"
15#include "OSBearcatSolverXij.h"
16#include "OSConfig.h"
17#include "OSResult.h"
18#include "OSOption.h"
19#include "OSiLReader.h"
20#include "OSiLWriter.h"
21#include "OSoLReader.h"
22#include "OSoLWriter.h"
23#include "OSrLReader.h"
24#include "OSrLWriter.h"
25#include "OSInstance.h"
26#include "OSFileUtil.h"
28
29
30
31#ifdef COIN_HAS_COUENNE
32#include "OSCouenneSolver.h"
33#endif
34
35#ifdef COIN_HAS_IPOPT
36#include "OSIpoptSolver.h"
37#endif
38
39
40
41
42#include <vector>
43#include <map>
44#include <sstream>
45
46using std::ostringstream;
47
48
49
50
52 m_osinstanceMaster(NULL) {
53 std::cout << "INSIDE OSColGenApp CONSTRUCTOR" << std::endl;
54
55}//end OSColGenApp Constructor
56
57
59 std::cout << "INSIDE OSColGenApp CONSTRUCTOR" << std::endl;
60 //std::cout << "the contructor things whichBlock = " << m_whichBlock<< std::endl;
61
62 //get parameters-options
63 //set default values:
64
66
67
68
69 m_osDecompParam.nodeLimit = 1000;
70 m_osDecompParam.columnLimit = 20000;
71 m_osDecompParam.masterColumnResetValue = 5000;
72 m_osDecompParam.zeroTol = .0001;
73 m_osDecompParam.artVarCoeff = 1000000;
74 m_osDecompParam.optTolPerCent = 0;
75
77 //get the options for the OSDecompSolver
79
80
81 m_osinstanceMaster = NULL;
82 m_osrouteSolver = NULL;
83
84 m_osrouteSolver = NULL;
86 std::cout << "CREATE THE FACTORY " << std::endl;
87 m_osrouteSolver = OSDecompSolverFactory::factories[ "OSBearcatSolverXij"]->create();
88 std::cout << "FINISH FACTORY CREATION " << std::endl;
89 std::cout << "SET FACTORY OPTION " << std::endl;
90 m_osrouteSolver->m_osoption = m_osoption;
91 std::cout << "FINISH SET FACTORY OPTION " << std::endl;
92 //m_osrouteSolver = new OSBearcatSolverXij( osoption);
93
94 //share the common parameters
95 m_osrouteSolver->m_osDecompParam = m_osDecompParam;
96 m_osrouteSolver->initializeDataStructures();
97
98
99 //initialize the bounds
101 m_zLB = -OSDBL_MAX;
102
103 //set the column number
104 m_numColumnsOld = 0;
105
106
107
108} //end OSColGenApp Constructor
109
110
112
113 std::cout << "INSIDE ~OSColGenApp DESTRUCTOR" << std::endl;
114
115 //kipp -- why doesn't m_osrouteSolver delete the
116 //m_osinstanceMaster object
117
118 if( m_osinstanceMaster != NULL) delete m_osinstanceMaster;
119
120 if( m_osrouteSolver != NULL) delete m_osrouteSolver;
121
122 //finally delete the factories
123
124 delete m_factoryInit;
125
126}//end ~OSColGenApp() destructor
127
128
130
131 m_osinstanceMaster = m_osrouteSolver->getInitialRestrictedMaster( );
132 //std::cout << m_osinstanceMaster->printModel( ) << std::endl;
133
134}//end generateInitialRestrictedMaster
135
136
137void OSColGenApp::getCuts(const double* thetaVar, const int numThetaVar,
138 int &numNewRows, int* &numNonz, int** &colIdx,
139 double** &values, double* &rowLB, double* &rowUB) {
140
141 m_osrouteSolver->getCutsTheta( thetaVar, numThetaVar,
142 numNewRows, numNonz, colIdx, values, rowLB, rowUB);
143
145 //for now let's always get these cuts, hence the default false
146 if(numNewRows == 0 && m_calledBranchAndBound == false
147 && m_osrouteSolver->m_numMultCuts <= m_osrouteSolver->m_multiCommodCutLimit) {
148 m_osrouteSolver->getCutsMultiCommod( thetaVar, numThetaVar,
149 numNewRows, numNonz, colIdx, values, rowLB, rowUB);
150
151
152 m_osrouteSolver->m_numMultCuts += numNewRows;
153
154 //double lhs;
155 //for(int i = 0; i < numNewRows; i++){
156 // lhs = 0;
157 //for(int j = 0; j < numNonz[ i]; j++){
158
159 // lhs += m_si->getColSolution()[ colIdx[i][j] ]*values[i][j];
160 // std::cout << " cut coefficient = " << values[i][j] << " theta value = " << m_si->getColSolution()[ colIdx[i][j] ] << std::endl;
161
162 //}
163
164 //std::cout << "LHS = " << lhs << std::endl;
165
166 //}// loop over number of new rows
167
168 //exit( 1);
169 }//end on if
170
171
172}//end getCuts
173
174
175void OSColGenApp::getColumns(const double* yA, const int numARows,
176 const double* yB, const int numBRows,
177 int &numNewColumns, int* &numNonz, double* &cost,
178 int** &rowIdx, double** &values, double &lowerBound) {
179
180 m_osrouteSolver->getColumns(yA, numARows,
181 yB, numBRows, numNewColumns, numNonz,
182 cost, rowIdx, values, lowerBound);
183
184
185
186
187}//end generateColumns
188
190
191
192 //get any options relevant to OSColGenApp
193
194 try{
195
196 std::vector<SolverOption*> solverOptions;
197 std::vector<SolverOption*>::iterator vit;
198
199 solverOptions = osoption->getSolverOptions("OSDecompSolver");
200 if (solverOptions.size() == 0) throw ErrorClass( "options for OSDecompSolver not available");
201
202
203 for (vit = solverOptions.begin(); vit != solverOptions.end(); vit++) {
204
205 if((*vit)->name.find("columnLimit") != std::string::npos){
206
207
208 std::istringstream columnLimitBuffer( (*vit)->value);
209 columnLimitBuffer >> m_osDecompParam.columnLimit;
210 std::cout << "columnLimit = " << m_osDecompParam.columnLimit << std::endl;
211
212 }else{
213
214 if( (*vit)->name.find("artVarCoeff") != std::string::npos ){
215
216 std::istringstream artVarCoeffBuffer( (*vit)->value);
217 artVarCoeffBuffer >> m_osDecompParam.artVarCoeff;
218 std::cout << "artVarCoeff = " << m_osDecompParam.artVarCoeff << std::endl;
219
220 }else{
221
222 if( (*vit)->name.find("zeroTol") != std::string::npos){
223
224 std::istringstream zeroTolBuffer( (*vit)->value);
225 zeroTolBuffer >> m_osDecompParam.zeroTol;
226 std::cout << "zeroTol = " << m_osDecompParam.zeroTol << std::endl;
227
228 }else{
229
230
231 if( (*vit)->name.find("nodeLimit") != std::string::npos){
232
233 std::istringstream nodeLimitBuffer( (*vit)->value);
234 nodeLimitBuffer >> m_osDecompParam.nodeLimit;
235 std::cout << "nodeLimit = " << m_osDecompParam.nodeLimit << std::endl;
236
237 }else{
238
239 if( (*vit)->name.find("masterColumnResetValue") != std::string::npos){
240
241 std::istringstream masterColumnResetValueBuffer( (*vit)->value);
242 masterColumnResetValueBuffer >> m_osDecompParam.masterColumnResetValue;
243 std::cout << "masterColumnResetValue = " << m_osDecompParam.masterColumnResetValue << std::endl;
244 }else{
245
246 if( (*vit)->name.find("optTolPerCent") != std::string::npos){
247
248 std::istringstream optTolPerCentBuffer( (*vit)->value);
249 optTolPerCentBuffer >> m_osDecompParam.optTolPerCent;
250 std::cout << "masterColumnResetValue = " << m_osDecompParam.optTolPerCent<< std::endl;
251 }
252 }
253 }
254 }
255 }
256 }
257 }
258
259 } catch (const ErrorClass& eclass) {
260
261 throw ErrorClass(eclass.errormsg);
262
263 }
264
265}//end getOptions
266
267
269
270
271
272 int *cbasis = NULL;
273 int *rbasis = NULL;
274 int *new_cbasis = NULL;
275
276
277
278 std::set<std::pair<int, double> >::iterator sit;
279 std::vector<int>::iterator vit;
280 std::map<int, int>::iterator mit;
281
282 int numCols;
283 int numRows;
284 int i;
285
286 //initialize upper bound
287 m_zUB = m_osrouteSolver->m_bestIPValue;
288
289 //initialize number of columns and nodes generated
290
293 std::cout << " m_zUB " << m_zUB << std::endl;
294
295 try{
296
297 // the solver
298
299 m_solver = new CoinSolver();
300
301 // the solver interface
302
303 //kipp -- later have clp be an option
304 //I guess for now it must be an Osi solver
305 m_solver->sSolverName ="clp";
306 //std::cout << m_osinstanceMaster->printModel( ) << std::endl;
307 m_solver->osinstance = m_osinstanceMaster;
308
309 //pass options
310 m_solver->osoption = m_osoption;
311 m_solver->buildSolverInstance();
312
313
314 //get the solver interface
315 m_si = m_solver->osiSolver;
316
317 m_yA = new double[m_osinstanceMaster->getConstraintNumber() ];
318 //kipp -- hard coding, come back and fix with option
319 //kipp -- do all of the newing in a separate routine
320 m_yB = new double[ m_osrouteSolver->m_maxBmatrixCon];
321
322
323 m_maxCols = m_osrouteSolver->m_maxMasterColumns;
324 m_maxRows = m_osrouteSolver->m_maxMasterRows;
325
326 m_theta = new double[ m_maxCols];
327
328
329 //get initial LP relaxation of master
331 //exit( 1);
332 //get the solution vector
333 numCols = m_si->getNumCols();
334 numRows = m_si->getNumRows();
335
336 //kipp -- just testing
337 cbasis = new int[ numCols];
338 rbasis = new int[ numRows ];
339 m_si->getBasisStatus( cbasis, rbasis);
340
341 for(i = 0; i < numCols; i++){
342 //get the basic primal variables
343 if(cbasis[ i] == 1) m_zOptRootLP.push_back( i);
344 //get the LP relaxation
345 *(m_theta + i) = m_si->getColSolution()[i];
346
347 m_zRootLPx_vals.push_back( *(m_theta + i) );
348
350 int j;
351 if( *(m_theta + i) > m_osDecompParam.zeroTol){
352 std::cout << "x variables for column " << i << std::endl;
353 for(j = m_osrouteSolver->m_thetaPnt[ i]; j < m_osrouteSolver->m_thetaPnt[ i + 1] ; j++){
354 std::cout << m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[ j] ] << " = " << *(m_theta + i) << std::endl;
355 }
356 }
358
359
360 }
361 m_zLB = m_si->getObjValue();
362 m_zRootLP = m_si->getObjValue();
363 //print LP value at node
364 std::cout << "optimal LP value at root node = " << m_zLB << std::endl;
365 //get the optimal LP root solution
366
367
368 //exit( 1);
369
370 for ( sit = m_osrouteSolver->intVarSet.begin() ;
371 sit != m_osrouteSolver->intVarSet.end(); sit++ ){
372
373 m_si->setInteger( sit->first);
374
375 }
376
377 CbcModel model( *m_si);
378 OsiSolverInterface *ipSolver = model.solver();
379 std::cout << "start solving master as integer program " << std::endl;
380 ipSolver->branchAndBound();
381 std::cout << "done solving master as integer program " << std::endl;
382 //CbcMain0( model);
383 //CbcMain1( 0, 0, model);
384 //kipp -- put in check to make sure we get an integer solution
385 if( ipSolver->getObjValue() < m_zUB) m_zUB = ipSolver->getObjValue() ;
386
387 //get the solution
388 numCols = m_si->getNumCols();
389
390 for(i = 0; i < numCols; i++){
391
392 //get the indexes of integer variables
393 if( model.getColSolution()[i] > m_osDecompParam.zeroTol){
394
395 m_zOptIndexes.push_back( i) ;
396
397 }
398
399 }
400
401 for ( sit = m_osrouteSolver->intVarSet.begin() ;
402 sit != m_osrouteSolver->intVarSet.end(); sit++ ){
403
404 m_si->setContinuous( sit->first);
405 m_si->setColUpper( sit->first, sit->second);
406
407 }
408
409 std::cout << "OPTIMAL LP VALUE = " << m_zLB << std::endl;
410 std::cout << "CURRENT BEST IP VALUE = " << m_zUB << std::endl;
411
413 //start reset
414 /*
415 int tmpCols = m_numColumnsGenerated;
416 resetMaster();
417 numCols = m_si->getNumCols();
418 new_cbasis = new int[ numCols ];
419 for(i = 0; i < numCols; i++) new_cbasis[ i] = 3;
420 for (vit = m_zOptRootLP.begin(); vit != m_zOptRootLP.end(); vit++ ) new_cbasis[ *vit ] = 1;
421 m_si->setBasisStatus( new_cbasis, rbasis );
422 solveRestrictedMasterRelaxation();
423 //m_si->initialSolve();
424 //exit( 1);
425 for(i = 0; i < numCols; i++) *(m_theta + i) = m_si->getColSolution()[i];
426 std::cout << "NUMBER OF NEW GENERATED COLUMNS = " << m_numColumnsGenerated - tmpCols << std::endl;
427 */
428 //end reset
430
431 m_osrouteSolver->m_rootLPValue = m_zRootLP;
432
433 //go into branch and bound
434 m_message = "";
435 std::cout << "START BRANCH AND BOUND = " << std::endl;
436 if(m_zLB + m_osDecompParam.zeroTol < m_zUB) branchAndBound();
437
438 //demand values
439 //m_osrouteSolver->m_demand;
440
441 std::cout << "FINISH BRANCH AND BOUND = " << std::endl;
443 m_osrouteSolver->m_bestLPValue = m_zLB;
444 m_osrouteSolver->m_bestIPValue = m_zUB;
445 if(m_message == "") m_message = "******** WE ARE OPTIMAL *******";
448
449
450 delete m_solver;
451
452 delete[] m_yA;
453 m_yA = NULL;
454
455 delete[] m_yB;
456 m_yB = NULL;
457
458 delete[] m_theta;
459 m_theta = NULL;
460
461
462 delete[] cbasis;
463 cbasis = NULL;
464 if(new_cbasis != NULL) delete[] new_cbasis;
465 new_cbasis = NULL;
466 delete[] rbasis;
467 rbasis = NULL;
468
469
470 } catch (const ErrorClass& eclass) {
471
472 delete m_solver;
473
474 delete[] m_yA;
475 m_yA = NULL;
476
477 delete[] m_yB;
478 m_yB = NULL;
479
480 delete[] m_theta;
481 m_theta = NULL;
482
483
484
485 delete[] cbasis;
486 cbasis = NULL;
487 if(new_cbasis != NULL) delete[] new_cbasis;
488 new_cbasis = NULL;
489 delete[] rbasis;
490 rbasis = NULL;
491
492 throw ErrorClass(eclass.errormsg);
493
494}
495
496}//end solve
497
498
500
501 int i;
502 int k;
503 //we include convexity constraints in this number
504 int numARows;
505 // dimension y to number of nodes
506 int numBRows;
507 int numCols;
508
509 //getColumns function call return parameters
510 int numNewColumns;
511 int* numNonz = NULL;
512 double* cost = NULL;
513 int** rowIdx = NULL;
514 double** values = NULL ;
515 double lowerBound;
516 //end of getColumns function call return parameters
517
518 double collb; // kipp -- put in getColumns
519 double colub; // kipp -- put in getColumns
520 //all of our m_theta columns have a lower bound of 0 and upper bound of 1
521 collb = 0.0;
522 colub = 1.0;
523 //kipp -- I would like to use OSDBL_MAX but Clp likes this better
524 //double bigM = 1.0e24;
525 double bigM = m_osDecompParam.artVarCoeff;
526 //getRows function call return parameters
527 int numNewRows;
528 int* numRowNonz = NULL;
529 int** colIdx = NULL;
530 double** rowValues = NULL ;
531 double* rowLB;
532 double* rowUB;
533 //end of getRows function call return parameters
534 //art variables
535
536 int rowArtIdx;
537 double rowArtVal;
538
539 bool isCutAdded;
540
541
542
543 try{
544 numARows = m_osrouteSolver->m_numNodes;
545
546
547 isCutAdded = true;
548
549 while(isCutAdded == true ){
550
551 isCutAdded = false;
552 //start out loop on if cuts found
553 std::cout << "CALL Solve " << " Number of columns = " << m_si->getNumCols() << std::endl;
554 //we are going through OS here, m_solver is a CoinSolver object
555 //now solve
556 m_solver->solve();
557 //m_si->initialSolve();
558 std::cout << "Solution Status = " << m_solver->osresult->getSolutionStatusType( 0 ) << std::endl;
559 //std::cout << m_solver->osrl << std::endl;
560
561 if(m_si->getNumRows() != numARows + m_osrouteSolver->m_numBmatrixCon ) {
562 std::cout << "m_si->getNumRows() = " << m_si->getNumRows() << std::endl;
563 std::cout << "numARows = " << numARows << std::endl;
564 std::cout << "m_numBmatrixCon = " << m_osrouteSolver->m_numBmatrixCon << std::endl;
565 throw ErrorClass("detect a row number inconsistency in solveRestrictedMasterRelaxation");
566 }
567
568
569
570 if(m_si->getRowPrice() == NULL )
571 throw ErrorClass("problem getting dual values in solveRestrictedMasterRelaxation");
572
573
574 numBRows = m_osrouteSolver->m_numBmatrixCon;
575
576 for(i = 0; i < numARows; i++){
577
578 *(m_yA + i) = m_si->getRowPrice()[ i];
579
580 }
581
582 for(i = numARows; i < numARows + numBRows; i++){
583
584 *(m_yB + i - numARows) = m_si->getRowPrice()[ i];
585
586 }
587
588 lowerBound = -1;
589 int loopKount = 0;
592 while(lowerBound < -m_osDecompParam.zeroTol && loopKount < m_osDecompParam.columnLimit){
593 loopKount++;
594
595 //kipp here is where the while loop goes
596 //start while loop
597 getColumns(m_yA, numARows, m_yB, numBRows, numNewColumns,
598 numNonz, cost, rowIdx, values, lowerBound);
599
600 std::cout << "Lower Bound = " << lowerBound << std::endl;
601
602
603
604 for(k = 0; k < numNewColumns; k++){
605
606 m_si->addCol( numNonz[ k], rowIdx[k], values[k],
607 collb, colub, cost[ k]) ;
608
609
611
612 }
613
614 std::cout << " OBJ VALUE = " << m_si->getObjValue() << std::endl;
615
616 std::cout << "m_zUB = " << m_zUB << std::endl;
617
618 if(lowerBound + m_si->getObjValue() > (1 - m_osDecompParam.optTolPerCent)*m_zUB) break;
619
620 std::cout << std::endl << std::endl << std::endl;
621 std::cout << "CALL Solve " << " Number of columns = " << m_si->getNumCols() << std::endl;
622 m_solver->solve();
623 //m_si->initialSolve();
624 std::cout << "Solution Status = " << m_solver->osresult->getSolutionStatusType( 0 ) << std::endl;
625 std::cout << "Number of solver interface columns = " << m_si->getNumCols() << std::endl;
626 //m_numNodes is number of artificial variables
627
628 numCols = m_si->getNumCols();
629
630 if( numCols != m_osrouteSolver->m_numThetaVar ) throw ErrorClass("number variables in solver not consistent with master");
631 if( numCols + m_osrouteSolver->m_numHubs >= m_maxCols) {
632
633 m_message = " ***** COLUMN LIMIT EXCEEDED -- INSIDE solveRestrictedMasterRelaxation ****";
635 m_osrouteSolver->m_bestLPValue = m_zLB;
636 m_osrouteSolver->m_bestIPValue = m_zUB;
639 throw ErrorClass("we ran out of columns");
640 }
641
642 for(i = 0; i < numARows; i++){
643
644 *(m_yA + i) = m_si->getRowPrice()[ i];
645
646 }
647
648 for(i = numARows; i < numARows + numBRows; i++){
649
650 *(m_yB + i - numARows) = m_si->getRowPrice()[ i];
651
652 }
653
654 }//end while on column generation
656
657 if( loopKount >= m_osDecompParam.columnLimit)throw ErrorClass("we exceeded loop kount in column generation");
658
659 //get a primal solution
660 numCols = m_si->getNumCols();
661 for(i=0; i < numCols; i++){
662 *(m_theta + i) = m_si->getColSolution()[i];
663 }
664
665
666 numNewRows = 0;
667
668 //do not get cuts if LP relaxation worse than upper bound
669 if(m_si->getObjValue() < (1 - m_osDecompParam.optTolPerCent)*m_zUB)
670 getCuts(m_theta, numCols, numNewRows, numRowNonz,
671 colIdx,rowValues, rowLB, rowUB);
672
673
674 if( numNewRows >= 1 ){
675
676 isCutAdded = true;
677
678 for(i = 0; i < numNewRows; i++){
679
680 m_si->addRow(numRowNonz[ i], colIdx[ i], rowValues[ i], rowLB[ i], rowUB[ i] ) ;
681
682 //add two variables for this row so we can never be infeasible
683 //m_si->addCol( numNonz, rowIdx[k], values[k],
684 // collb, colub, cost[ k]) ;
685
686 //add the artificial variable for the UB
687 rowArtVal = -1.0;
688 rowArtIdx = m_si->getNumRows() - 1;
689 //m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
690 //m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, 1, bigM);
691 //add the artificial variable for the LB
692 rowArtVal = 1.0;
693 //m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
694 m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, 1, bigM);
696
697 }
698
699 std::cout << std::endl;
700 std::cout << "CUTS WERE ADDED CALL SOLVE" << std::endl;
701 m_solver->solve();
702
703 } // end if on whether or not we added cuts
704
705
706
707
708 }//end while on isCutAdded
709
710
711
712 } catch (const ErrorClass& eclass) {
713
714 throw ErrorClass(eclass.errormsg);
715
716 }
717
718
719}//end solveRestrictedMasterRelaxation
720
721
722bool OSColGenApp::isInteger( const double *thetaVar, const int numThetaVar,
723 const double tol){
724
725
726 bool isInt;
727 isInt = true;
728 int i;
729
730 try{
731
732 for(i = 0; i < numThetaVar; i++){
733
734 if( (thetaVar[ i] > tol) && (thetaVar[ i] < 1 - tol) ){
735
736 isInt = false;
737 break;
738 }
739
740 }
741
742 return isInt;
743
744 } catch (const ErrorClass& eclass) {
745
746 throw ErrorClass(eclass.errormsg);
747
748 }
749
750
751
752}//end isInteger
753
754
756
757 int numCols;
758 int numRows;
759 std::set<std::pair<int, double> >::iterator sit;
760 int i;
761
762 numCols = m_si->getNumCols();
763
764
765 for(i = 0; i < numCols; i++){
766
767 std::cout << "PROCESSING THETA COLUMN " << i << " value = " << m_si->getColSolution()[i] << std::endl;
768
769 for(int j = m_osrouteSolver->m_thetaPnt[ i]; j < m_osrouteSolver->m_thetaPnt[ i + 1]; j++ ){
770
771 //std::cout << m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[ j] ] << std::endl;
772
773 }
774 }
775
776 numRows = m_si->getNumRows();
777
778 for(i = m_osrouteSolver->m_numNodes; i < numRows; i++){
779
780 std::cout << "PROCESSING ROW " << i << std::endl;
781
782 for(int j = m_osrouteSolver->m_pntBmatrix[ i - m_osrouteSolver->m_numNodes]; j < m_osrouteSolver->m_pntBmatrix[ i + 1 - m_osrouteSolver->m_numNodes]; j++ ){
783
784 //std::cout << m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_Bmatrix[ j] ] << std::endl;
785
786 }
787 }
788 //check integer variables and upper bounds -- loop over integer variable set
789
790 //
791 for ( sit = m_osrouteSolver->intVarSet.begin() ;
792 sit != m_osrouteSolver->intVarSet.end(); sit++ ){
793
794 //std::cout << "Integer variable " << sit->first << " Upper Bound = " << sit->second << std::endl;
795
796 }
797}//end printDebugInfo
798
799
801
803
808 std::map<int, int> varConMap;
809
810 std::vector<OSNode*> nodeVec;
811 std::vector<OSNode*>::iterator vit;
812
813
814 std::map<int, OSNode*>::iterator mit;
815 int bestNodeID;
816 double bestNodeBound;
817
818
819 OSNode *osnode = NULL;
820 OSNode *osnodeLeftChild = NULL;
821 OSNode *osnodeRightChild = NULL;
822
823 bool bandbWorked;
824 bandbWorked = true;
825 int numCols;
826 int rowIdx;
827 rowIdx = 0;
828
829 bool leftNodeCreated = false;
830 bool rightNodeCreated = false;
831
832
833
834 try{
835
836 //get the solution
837 numCols = m_si->getNumCols();
838 //kipp -- imporant this is now found earliear.
839 //for(i = 0; i < numCols; i++){
840 // //get the LP relaxation
841 // std::cout << "theta = " << *(m_theta + i) << std::endl;
842 //}
843
844 //NOTE -- we must know theta here
845
846 //create a branching cut
847 createBranchingCut(m_theta, numCols, varConMap, rowIdx);
848
849
850
852
853 osnodeLeftChild = createChild(osnode, varConMap, rowIdx, 1, 1);
854 if(osnodeLeftChild != NULL){
855 //finally set the nodeID
856 //and record parent ID
857 //m_numNodesGenerated++;
858 osnodeLeftChild->nodeID = m_numNodesGenerated;
859 osnodeLeftChild->parentID = 0;
860 //nodeVec.push_back( osnodeLeftChild);
861 m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeLeftChild->nodeID, osnodeLeftChild) );
862
863 }
864
866
868
869 osnodeRightChild = createChild(osnode, varConMap, rowIdx, 0, 0);
870 if(osnodeRightChild != NULL){
871 //finally set the nodeID
872 //and record parent ID
873 //m_numNodesGenerated++;
874 osnodeRightChild->nodeID = m_numNodesGenerated;
875 osnodeRightChild->parentID = 0;
876 //nodeVec.push_back( osnodeRightChild);
877 m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeRightChild->nodeID, osnodeRightChild) );
878 }
879
881
882 // now loop
883 //kipp -- make this an option
884
885 std::cout << "ENTERING THE WHILE IN BRANCH AND BOUND" << std::endl;
886 std::cout << "m_numNodesGenerated = " << m_numNodesGenerated << std::endl;
887 //while( (nodeVec.size() > 0) && (m_numNodesGenerated <= nodeLimit) ){
888 while( m_nodeMap.size() > 0 ){
889
890 if(m_numNodesGenerated > m_osDecompParam.nodeLimit ){
891 m_message = "******* NODE LIMIT EXCEEDED *******";
892 return false;
893 }
894
895
896 if( m_si->getNumCols() > m_maxCols ){
897 m_message = "******* COLUMN LIMIT EXCEEDED *******";
898 return false;
899 }
900
901 //kipp -- experimental
902 //m_osDecompParam.masterColumnResetValue = 3000;
903 //if( m_si->getNumCols() > 200000) {
905 m_osDecompParam.masterColumnResetValue) {
907 std::cout << "DOING A MASTER RESET IN BRANCH AND BOUND" << std::endl;
908 std::cout << "NUMBER OF COLUMNS BEFORE RESET = " << m_si->getNumCols() << std::endl;
909 resetMaster();
910 std::cout << "NUMBER OF COLUMNS AFTER RESET = " << m_si->getNumCols() << std::endl;
911 //int tmpCols = m_numColumnsGenerated;
912 //solveRestrictedMasterRelaxation();
913 //std::cout << "NUMBER OF NEW GENERATED COLUMNS = " << m_numColumnsGenerated - tmpCols << std::endl;
914 //exit( 1);
915 }
916
917 leftNodeCreated = false;
918 rightNodeCreated = false;
919 //grab a node -- for now the last node, we do FIFO
920 //osnode = nodeVec.back();
921
922 //let's loop and find node with the largest nodeID -- this will
923 //corespond to fifo
924
925 bestNodeID = -1;
926 bestNodeBound = OSDBL_MAX;
927 //mit->first is the the OSNode nodeID
928 //mit->second is an OSNode
929 for (mit = m_nodeMap.begin(); mit != m_nodeMap.end(); mit++ ){
930
931 //FIFO criterions
932 //if( mit->second->nodeID > bestNodeID) bestNodeID = mit->second->nodeID;
933
934 //Best Bound criterion
935 if( mit->second->lpValue < bestNodeBound) {
936
937 bestNodeBound = mit->second->lpValue;
938 bestNodeID = mit->first;
939 //note same as:bestNodeID = mit->second->nodeID;
940
941
942 }
943
944 }
945
946 //get the node
947 mit = m_nodeMap.find( bestNodeID );
948 if(mit == m_nodeMap.end() ) throw ErrorClass("a node selection problem in branch and bound");
949 osnode = mit->second;
950
951 if( osnode->lpValue < (1 - m_osDecompParam.optTolPerCent)*m_zUB - m_osDecompParam.zeroTol){
952
953
954 //create a branching cut
955 std::cout << "CREATE A BRANCHING CUT " << std::endl;
956 createBranchingCut(osnode->thetaIdx, osnode->theta, osnode->thetaNumNonz,
957 varConMap, rowIdx);
958
960 /*
961 std::map<int, int>::iterator tmpit;
962 for (tmpit = varConMap.begin() ; tmpit != varConMap.end(); tmpit++ ){
963
964 std::cout << std::endl;
965 std::cout << "LOOPING OVER VARIABLE " << m_osrouteSolver->m_variableNames[ tmpit->first ] << std::endl;
966 std::cout << "ROW UB = " << osnode->rowUB[ tmpit->second] << std::endl;
967 std::cout << "ROW LB = " << osnode->rowLB[ tmpit->second] << std::endl;
968 kippster
969 }
970 */
972
973 std::cout << "BEST NODE ID " << bestNodeID << std::endl;
974 std::cout << "NODE LP VALUE = " << osnode->lpValue << std::endl;
975 //check for node consistency
976 checkNodeConsistency( rowIdx, osnode);
977 // create children
978 //create the left node
979
980 osnodeLeftChild = createChild(osnode, varConMap, rowIdx, 1, 1);
981 if(osnodeLeftChild != NULL){
982 //finally set the nodeID
983 //and record parent ID
984 //m_numNodesGenerated++;
985 osnodeLeftChild->nodeID = m_numNodesGenerated;
986 osnodeLeftChild->parentID = osnode->nodeID;
987 leftNodeCreated = true;
988 }
989
990 //create the right node
991 osnodeRightChild = createChild(osnode, varConMap, rowIdx, 0, 0);
992 if(osnodeRightChild != NULL){
993 //finally set the nodeID
994 //and record parent ID
995 //m_numNodesGenerated++;
996 osnodeRightChild->nodeID = m_numNodesGenerated;
997 osnodeRightChild->parentID = osnode->nodeID;
998 rightNodeCreated = true;
999 }
1000
1001
1002 //nodeVec.erase( nodeVec.end() - 1) ;
1003 m_nodeMap.erase( mit);
1004 delete osnode;
1005
1006 //if( leftNodeCreated == true) nodeVec.push_back( osnodeLeftChild) ;
1007 //if( rightNodeCreated == true) nodeVec.push_back( osnodeRightChild) ;
1008
1009 if( leftNodeCreated == true)
1010 m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeLeftChild->nodeID, osnodeLeftChild) ) ;
1011
1012 if( rightNodeCreated == true)
1013 m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeRightChild->nodeID, osnodeRightChild) ) ;
1014
1015
1016 }else{
1017
1018 //fathom node by virtue of the upper bound
1019 std::cout << "FATHAM BY UPPER BOUND " << std::endl;
1020 //nodeVec.erase( nodeVec.end() - 1) ;
1021 m_nodeMap.erase( mit);
1022 delete osnode;
1023
1024 }//end if on lp bound check
1025
1026 //kipp -- critical reset upper and lower bounds
1027 //kipp don't forget to erase the parent node
1028
1029 }//end the while
1030
1031
1032
1033 if(m_numNodesGenerated > 0){
1034
1035 m_zLB = (1 - m_osDecompParam.optTolPerCent)*m_zUB;
1036 }else{
1037
1038 m_zLB = m_zUB;
1039 }
1040
1041
1042 //exit( 1);
1043
1044 return bandbWorked;
1045
1046 } catch (const ErrorClass& eclass) {
1047
1048 throw ErrorClass(eclass.errormsg);
1049
1050 }
1051
1052}// end branchAndBound
1053
1054OSNode* OSColGenApp::createChild(const OSNode *osnodeParent, std::map<int, int> &varConMap,
1055 const int rowIdx, const double rowLB, const double rowUB){
1056
1058
1059 OSNode *osnodeChild;
1060 osnodeChild = NULL;
1061 int numRows;
1062 int numCols;
1063
1064 int tmpColNum ;
1065 int tmpRowNum ;
1066
1067 std::map<int, int>::iterator mit;
1068
1069
1070
1071 int i;
1072 int k;
1073 int childRowIdxNumNonz;
1074 childRowIdxNumNonz = 0;
1075
1076 //we want to store the solution vector (theta space)
1077 //in sparse format
1078 int thetaNumNonz;
1079
1080
1081 try{
1082
1083 if(osnodeParent != NULL) childRowIdxNumNonz = osnodeParent->rowIdxNumNonz + 1;
1084 else childRowIdxNumNonz = 1;
1085
1086 //set upper and lower bounds correctly
1087 //set the parent values
1088 if(osnodeParent != NULL){
1089 for(i = 0; i < osnodeParent->rowIdxNumNonz; i++){
1090
1091
1092 m_si->setRowLower( osnodeParent->rowIdx[ i], osnodeParent->rowLB[ i]);
1093 m_si->setRowUpper( osnodeParent->rowIdx[ i], osnodeParent->rowUB[ i]);
1094
1095
1096 }
1097 }
1098 //set the new value
1099 m_si->setRowLower( rowIdx, rowLB);
1100 m_si->setRowUpper( rowIdx, rowUB);
1101 //now solve
1102
1103 //print out the restricted master
1104
1105 //if(rowUB == 0) m_si->writeLp( "gailTest2" );
1106
1107 //exit( 1);
1108 std::cout << "CALL SOLVE FROM CREATE CHILD " << std::endl;
1109 //kipp -- important, you really need to verify that an optimal solution was obtained!!!
1110 if(osnodeParent != NULL){
1111
1112 tmpColNum = m_si->getNumCols() ;
1113 tmpRowNum = m_si->getNumRows() ;
1114 int *tmpColParent = new int[ tmpColNum];
1115 int *tmpRowParent = new int[ tmpRowNum ];
1116
1117 for(k = 0; k < tmpColNum; k++){
1118
1119 if( m_si->getObjCoefficients()[k] >=
1120 m_osDecompParam.artVarCoeff - m_osDecompParam.zeroTol)
1121 tmpColParent[ k ] = 3;
1122 //else if( osnodeParent->reducedCostIdx.find(k) == osnodeParent->reducedCostIdx.end() )
1123 // tmpColParent[ k ] = 3;
1124 else tmpColParent[ k] = 0;
1125
1126 }
1127
1128 for(k = 0; k < tmpRowNum; k++){
1129
1130 tmpRowParent[ k] = 0;
1131 }
1132
1133
1134
1135 //for(k = 0; k < osnodeParent->colBasisStatus.size(); k++){
1136
1137 //make the basis status of artificial variables 3
1138 //that is, nonbasic at lower bound
1139 //if( m_si->getObjCoefficients()[ osnodeParent->colBasisStatus[k].first]
1140 // >= m_osDecompParam.artVarCoeff - m_osDecompParam.zeroTol)
1141 // tmpColParent[osnodeParent->colBasisStatus[k].first ] = 3;
1142 //else
1143 //tmpColParent[osnodeParent->colBasisStatus[k].first ] =
1144 // osnodeParent->colBasisStatus[k].second;
1145 //}
1146
1147 m_si->setBasisStatus(tmpColParent, tmpRowParent);
1149
1150 //kippster extra error checking
1151 //kippster check on upper and lower bound
1152 //we are in rowIdx and the theta here should correspond to the same xijk
1153 //getRowActivity()
1154
1155 /*
1156 m_si->initialSolve();
1157 if(m_si->getRowActivity()[ rowIdx] > rowLB ||
1158 m_si->getRowActivity()[ rowIdx] < rowUB ) {
1159
1160 std::cout << "Row lower bound = " << rowLB << std::endl;
1161 std::cout << "Row upper bound = " << rowUB << std::endl;
1162 std::cout << "Row activity = " << m_si->getRowActivity()[ rowIdx] << std::endl;
1163 throw ErrorClass( "Violating a branching cut UB and LB");
1164
1165 }
1166 */
1167
1168 //first get the column index nonzero elements in row rowIdx
1169
1170
1171
1175
1176 delete[] tmpColParent;
1177 tmpColParent = NULL;
1178 delete[] tmpRowParent;
1179 tmpRowParent = NULL;
1180
1181 //solveRestrictedMasterRelaxation( osnodeParent->colBasisStatus,
1182 // osnodeParent->rowBasisStatus);
1183 } else {
1184
1186 }
1187
1188
1189 std::cout << std::endl << std::endl;
1190 std::cout << "FINISH SOLVING THE MASTER " << std::endl;
1191
1192
1193 //
1194 //now reset the upper and lower bound
1195 //set the parent values
1196 if( osnodeParent != NULL){
1197 for(i = 0; i < osnodeParent->rowIdxNumNonz; i++){
1198
1199
1200 m_si->setRowLower( osnodeParent->rowIdx[ i], 0);
1201 m_si->setRowUpper( osnodeParent->rowIdx[ i], 1);
1202
1203
1204 }
1205 }
1206 //reset the new value
1207 m_si->setRowLower( rowIdx, 0);
1208 m_si->setRowUpper( rowIdx, 1);
1209
1210 // let's try and fathom the node
1211 // if we are not as good a upper bound
1212 // we fathom, if we are integer we fathom
1213 std::cout << std::endl << std::endl;
1214 std::cout << "MESSAGE: START CREATION OF A CHILD NODE" << std::endl;
1215 std::cout << "LB " << rowLB << " UB = " << rowUB << std::endl;
1216 std::cout << "MESSAGE: LP RELAXATION VALUE OF POTENTIAL CHILD NODE " << m_si->getObjValue() << std::endl;
1217 std::cout << "MESSAGE: OPTIMALITY STATUS OF NODE IS " << m_si->isProvenOptimal() << std::endl;
1218
1219 if( m_si->getObjValue() < (1 - m_osDecompParam.optTolPerCent)*m_zUB - m_osDecompParam.zeroTol && m_si->isProvenOptimal() == 1) {
1220 // okay cannot fathom based on bound try integrality
1221 std::cout << "MESSAGE: WE CANNOT FATHOM THE CHILD BASED ON UPPER BOUND " << std::endl;
1222 numCols = m_si->getNumCols();
1223 numRows = m_si->getNumRows();
1224 thetaNumNonz = 0;
1225
1226 for(i = 0; i < numCols; i++){
1227 //get the LP relaxation
1228 *(m_theta + i) = m_si->getColSolution()[i];
1229 if( *(m_theta + i) > m_osDecompParam.zeroTol) thetaNumNonz++;
1230
1231 }
1232 if( isInteger( m_theta, numCols, m_osDecompParam.zeroTol) == true){
1233 //fathom by integrality
1234 std::cout << "MESSAGE: WE HAVE AN INTEGRALITY FATHOM " << m_zUB << std::endl;
1235 if( m_zUB > m_si->getObjValue() ){
1236
1237 m_zUB = m_si->getObjValue() ;
1238 //clear out out solution vector
1239 if( m_zOptIndexes.size() > 0) m_zOptIndexes.clear();
1240
1241 for(i = 0; i < numCols; i++){
1242
1243 if( *(m_theta + i) > m_osDecompParam.zeroTol) m_zOptIndexes.push_back( i) ;
1244
1245 }
1246 }
1247
1248 }else{
1249 //create the child node
1250 std::cout << "MESSAGE: WE ARE CREATING A CHILD NODE WITH NUMBER COLUMNS = "<< numCols << std::endl;
1251 osnodeChild = new OSNode(childRowIdxNumNonz, thetaNumNonz );
1252
1253
1254 //set the basis
1255 /*
1256 tmpColNum = m_si->getNumCols() ;
1257 tmpRowNum = m_si->getNumRows() ;
1258 int *tmpColChild = new int[ tmpColNum];
1259 int *tmpRowChild = new int[ tmpRowNum ];
1260
1261 m_si->getBasisStatus(tmpColChild, tmpRowChild);
1262
1263 for(k = 0; k < tmpColNum; k++){
1264
1265 osnodeChild->colBasisStatus.push_back(std::make_pair(k, tmpColChild[ k]) );
1266
1267 }
1268
1269 for(k = 0; k < tmpRowNum; k++){
1270
1271 osnodeChild->rowBasisStatus.push_back(std::make_pair(k, tmpRowChild[ k]) );
1272
1273 }
1274
1275 delete[] tmpColChild;
1276 tmpColChild = NULL;
1277
1278 delete[] tmpRowChild;
1279 tmpRowChild = NULL;
1280 */
1281
1282 //now set bound arrays
1283 if(osnodeParent == NULL){
1284 osnodeChild->rowIdx[ 0] = rowIdx;
1285 if(rowLB <= m_osDecompParam.zeroTol) osnodeChild->rowLB[ 0] = 0;
1286 else osnodeChild->rowLB[ 0] = 1;
1287
1288 if(rowUB <= m_osDecompParam.zeroTol) osnodeChild->rowUB[ 0] = 0;
1289 else osnodeChild->rowUB[ 0] = 1;
1290
1291
1292 }else{
1293 //set old values
1294 for(i = 0; i < osnodeParent->rowIdxNumNonz; i++){
1295
1296 osnodeChild->rowIdx[ i] = osnodeParent->rowIdx[ i];
1297 osnodeChild->rowLB[ i] = osnodeParent->rowLB[ i];
1298 osnodeChild->rowUB[ i] = osnodeParent->rowUB[ i];
1299
1300 }
1301 //set new value
1302
1303 osnodeChild->rowIdx[ childRowIdxNumNonz - 1] = rowIdx;
1304
1305
1306
1307 if(rowLB <= m_osDecompParam.zeroTol) osnodeChild->rowLB[ childRowIdxNumNonz - 1 ] = 0;
1308 else osnodeChild->rowLB[ childRowIdxNumNonz - 1 ] = 1;
1309
1310 if(rowUB <= m_osDecompParam.zeroTol) osnodeChild->rowUB[ childRowIdxNumNonz - 1 ] = 0;
1311 else osnodeChild->rowUB[ childRowIdxNumNonz - 1 ] = 1;
1312
1313
1314
1315 }
1316 //set child lp value
1317 osnodeChild->lpValue = m_si->getObjValue();
1318 //set theta
1319 thetaNumNonz = 0;
1320 for(i = 0; i < numCols; i++){
1321
1322 if( *(m_theta + i) > m_osDecompParam.zeroTol){
1323
1324 osnodeChild->thetaIdx[ thetaNumNonz] = i;
1325 osnodeChild->theta[ thetaNumNonz] = *(m_theta + i);
1326
1327 thetaNumNonz++;
1328 //std::cout << "x variables for column " << i << std::endl;
1329 //for(int j = m_osrouteSolver->m_thetaPnt[ i ]; j < m_osrouteSolver->m_thetaPnt[ i + 1] ; j++){
1330 // std::cout << m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[ j] ] << " = " << *(m_theta + i) << std::endl;
1331 //}
1332 }
1333
1334 //add the reduced costs
1335
1336 if(m_si->getReducedCost()[i] < (m_zUB - osnodeChild->lpValue) ) osnodeChild->reducedCostIdx.insert( i);
1337
1338 }
1339 }//end else on isInteger
1340 }// end if on upper bound test
1341
1342 std::cout << std::endl << std::endl;
1343 return osnodeChild;
1344
1345 } catch (const ErrorClass& eclass) {
1346
1347
1348 throw ErrorClass(eclass.errormsg);
1349
1350 }
1351
1352}//end createChild
1353
1354
1355void OSColGenApp::createBranchingCut(const int* thetaIdx, const double* theta,
1356 const int numThetaVar, std::map<int, int> &varConMap, int &rowIdx){
1357
1358 int varIdx;
1359 int numNonz;
1360 int* indexes;
1361 double* values;
1362
1363
1364 //for(int i = 0; i < numThetaVar; i++){
1365 // std::cout << "x variables for column " << thetaIdx[i] << std::endl;
1366 // for(int j = m_osrouteSolver->m_thetaPnt[ thetaIdx[ i] ]; j < m_osrouteSolver->m_thetaPnt[ thetaIdx[ i] + 1] ; j++){
1367 // std::cout << m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[ j] ] << " = " << theta[ i] << std::endl;
1368 // }
1369 //}
1370
1371 //kipp -- I would like to use OSDBL_MAX but Clp likes this better
1372 //double bigM = 1.0e24;
1373 double bigM = m_osDecompParam.artVarCoeff;
1374 double rowArtVal ;
1375
1376 std::map<int, int>::iterator mit;
1377
1378 //get the branching cut information
1379
1380 //
1381 //for(int i = 0; i < numThetaVar; i++ ) std::cout << "theta idx = " << thetaIdx[ i] << " theta = " << theta[ i] << std::endl;
1382
1383
1384 m_osrouteSolver->getBranchingCut(thetaIdx, theta, numThetaVar,
1385 varConMap, varIdx, numNonz, indexes, values);
1386
1387
1388
1389 //std::cout << "varIDX = " << varIdx << std::endl;
1390 //std::cout << "numNonz1 = " << numNonz << std::endl;
1391
1392
1393
1394
1395 //if numNonz is greater than zero:
1396 // 1) add add new variable to map -- at this point varConMap is empty
1397 // 2) add constraint then add to the formulation
1398 // 3) add variables
1399
1400
1401
1402 if( numNonz >0){
1403
1404
1405 //add the row
1406 //make upper and lower bound 0 and 1 first
1407 m_si->addRow(numNonz, indexes, values, 0, 1) ;
1408
1409 //add the artificial variables
1410 //add the artificial variable for the UB
1411 rowArtVal = -1.0;
1412 rowIdx = m_si->getNumRows() - 1;
1413
1414 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
1415 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
1416 //add the artificial variable for the LB
1417 rowArtVal = 1.0;
1418 m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
1419 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
1421
1422 //insert into map -- this is the first variable
1423 varConMap.insert( std::pair<int, int>(varIdx , rowIdx) );
1424
1425 m_rowIdxVarMap.insert( std::pair<int, int>(rowIdx , varIdx) );
1426
1427
1428
1429 } else{
1430 //the variable varIdx is in the map
1431 //get the constraint associated with this variable
1432 //throw and exception if varIdx not a key
1433
1434 mit = varConMap.find( varIdx);
1435 if( mit == varConMap.end() ) throw ErrorClass("in branchAndBound getBranchingCut() returned inconsistent value for varIdx");
1436 else rowIdx = mit->second;
1437
1438
1439 }//end if on numNonz
1440
1441
1442
1443}//end createBranchingCut Sparse
1444
1445
1446
1447void OSColGenApp::createBranchingCut(const double* theta,
1448 const int numThetaVar, std::map<int, int> &varConMap, int &rowIdx){
1449
1450 int varIdx;
1451 int numNonz;
1452 int* indexes;
1453 double* values;
1454
1455 //kipp -- I would like to use OSDBL_MAX but Clp likes this better
1456 //double bigM = 1.0e24;
1457 double bigM = m_osDecompParam.artVarCoeff;
1458 double rowArtVal ;
1459
1460 std::map<int, int>::iterator mit;
1461
1462 //get the branching cut information
1463 m_osrouteSolver->getBranchingCut( theta, numThetaVar,
1464 varConMap, varIdx, numNonz, indexes, values);
1465
1466
1467 std::cout << "varIDX2 = " << varIdx << std::endl;
1468 std::cout << "numNonz2 = " << numNonz << std::endl;
1469
1470
1471 //for(int i = 0; i < numNonz; i++){
1472 // std::cout << "x variables for column " << indexes[i] << std::endl;
1473 // for(int j = m_osrouteSolver->m_thetaPnt[ indexes[ i] ]; j < m_osrouteSolver->m_thetaPnt[ indexes[ i] + 1] ; j++){
1475 // }
1476 //}
1477
1478 //if numNonz is greater than zero:
1479 // 1) add add new variable to map -- at this point varConMap is empty
1480 // 2) add constraint then add to the formulation
1481 // 3) add artificial variables
1482
1483 if( numNonz >0){
1484
1485 //add the row
1486 //make upper and lower bound 0 and 1 first
1487 m_si->addRow(numNonz, indexes, values, 0, 1) ;
1488
1489 //add the artificial variables
1490 //add the artificial variable for the UB
1491 rowArtVal = -1.0;
1492 rowIdx = m_si->getNumRows() - 1;
1493
1494 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
1495 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
1496 //add the artificial variable for the LB
1497 rowArtVal = 1.0;
1498
1499 m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
1500 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
1502
1503 //insert into map -- this is the first variable
1504 varConMap.insert ( std::pair<int,int>(varIdx , rowIdx) );
1505 m_rowIdxVarMap.insert( std::pair<int, int>(rowIdx , varIdx) );
1506
1507
1508
1509 } else{
1510 //the variable varIdx is in the map
1511 //get the constraint associated with this variable
1512 //throw and exception if varIdx not a key
1513
1514 mit = varConMap.find( varIdx);
1515 if( mit == varConMap.end() ) throw ErrorClass("in branchAndBound getBranchingCut() returned inconsistent value for varIdx");
1516 else rowIdx = mit->second;
1517
1518
1519 }//end if on numNonz
1520
1521
1522
1523}//end createBranchingCut Dense
1524
1525
1527
1528 //kipp -- temp stuff here delete later
1530
1531 //std::map<int, int> inVars;
1532 std::map<int, int>::iterator mit;
1533 std::set<int>::iterator sit;
1534 std::vector<int>::iterator vit;
1535 std::vector<std::pair<int, int> >::iterator vit2;
1536 std::map<int, OSNode*>::iterator mit2;
1537 int i;
1538 int kount = 0;
1539
1540
1541 inVars.clear();
1542
1543 try{
1544
1545 //first add the columns corresponding to the root node solution
1546
1547 //for(i = 0; i < m_si->getNumCols(); i++){
1548
1549 //if( m_si->getColSolution()[i] > m_osDecompParam.zeroTol) inVars.insert( std::pair<int, int>(i, kount++) );
1550 //if( m_si->getObjCoefficients()[i] < 10000) inVars.insert( std::pair<int, int>(i, kount++) );
1551 //}
1552
1553 for(vit = m_zOptRootLP.begin() ; vit != m_zOptRootLP.end(); vit++){
1554
1555 inVars.insert( std::pair<int, int>(*vit, kount++) );
1556
1557 }
1558
1559
1560
1561 //next add the integer varaibles in the best known integer solution
1562 for(vit = m_zOptIndexes.begin() ; vit != m_zOptIndexes.end(); vit++){
1563
1564 if( inVars.find( *vit ) == inVars.end() ) inVars.insert( std::pair<int, int>(*vit, kount++) );
1565
1566 //for(int k = m_osrouteSolver->m_thetaPnt[*vit]; k < m_osrouteSolver->m_thetaPnt[*vit + 1]; k++){
1567
1568 //std::cout << m_osrouteSolver-> m_variableNames[ m_osrouteSolver->m_thetaIndex[ k] ] << std::endl;
1569
1570 //}
1571 }
1572
1573
1574
1575 //now loop over the nodes in the branch and bound tree
1576 //kipp -- this is hardcoded play with it later
1577 double tmpEps = OSDBL_MAX;
1578 for(mit2 = m_nodeMap.begin(); mit2 !=m_nodeMap.end(); mit2++){
1579
1580 std::cout << "NUMBER OF REDUCED COSTS = " << mit2->second->reducedCostIdx.size() << std::endl;
1581 for(sit = mit2->second->reducedCostIdx.begin();
1582 sit != mit2->second->reducedCostIdx.end(); sit++){
1583
1584 if( ( inVars.find( *sit ) == inVars.end() )
1585 && (m_si->getObjCoefficients()[*sit] < m_osDecompParam.artVarCoeff)
1586 && (m_si->getReducedCost()[*sit] < tmpEps )
1587 )
1588 inVars.insert( std::pair<int, int>(*sit, kount++) );
1589
1590
1591 }
1592
1593 //insert the thetat variables
1594 for(i = 0; i < mit2->second->thetaNumNonz; i++){
1595
1596 if( inVars.find( mit2->second->thetaIdx[ i] ) == inVars.end() )
1597
1598 inVars.insert( std::pair<int, int>(mit2->second->thetaIdx[ i], kount++) );
1599
1600 }
1601
1602 }
1603
1604
1605
1606 //for(mit = inVars.begin(); mit != inVars.end(); mit++) std::cout << "mit->first " << mit->first << " mit->second " << mit->second << std::endl;
1607 std::cout << "NUMBER OF COLUMNS = " << inVars.size() << std::endl;
1608 std::cout << "CALLING osroute solver reset " << std::endl;
1609 m_osrouteSolver->resetMaster( inVars, m_si );
1610 //std::cout << "THE MAPPING AFTER A RESET: " << std::endl;
1611 //for(mit = inVars.begin(); mit != inVars.end(); mit++) std::cout << "mit->first " << mit->first << " mit->second " << mit->second << std::endl;
1612
1613 //int numVars = m_osrouteSolver->m_osinstanceMaster->getVariableNumber();
1614 int numVars = m_si->getNumCols();
1615 double *tmpVals = NULL;
1616 tmpVals = new double[ numVars];
1617
1618 for(i = 0; i < numVars; i++) tmpVals[ i ] = 0;
1619
1620 for(mit = inVars.begin(); mit != inVars.end(); mit++){
1621 //tmpVals now point to old index values
1622 tmpVals[ mit->second] = m_theta[ mit->first] ;
1623 }
1624
1625
1626 for(i = 0; i < numVars; i++) m_theta[ i] = tmpVals[ i] ;
1627
1628 //reset the nodes in the branch and bound tree
1629
1630 for(mit2 = m_nodeMap.begin(); mit2 !=m_nodeMap.end(); mit2++){
1631
1632
1633 for(i = 0; i < mit2->second->thetaNumNonz; i++){
1634
1635 //inVars[mit2->second->thetaIdx[ i] ]
1636 if( inVars.find( mit2->second->thetaIdx[ i] ) == inVars.end() ) throw ErrorClass("index problem in resetMaster");
1637
1638 //kipp check to make sure we do not index an aritifical variable
1639 mit2->second->thetaIdx[ i] = inVars[ mit2->second->thetaIdx[ i] ] ;
1640 //if( inVars.find( mit2->second->thetaIdx[ i] ) != inVars.end() )
1641 // inVars.insert( std::pair<int, int>(mit2->second->thetaIdx[ i], kount++) );
1642
1643 }
1644
1645
1646 for(vit2 = mit2->second->colBasisStatus.begin();
1647 vit2 != mit2->second->colBasisStatus.end(); vit2++){
1648
1649 (*vit2).first = inVars[ (*vit2).first ] ;
1650
1651
1652 }
1653
1654 //reset reduced cost indexes
1655 std::set<int> tmpSet;
1656 for(sit = mit2->second->reducedCostIdx.begin();
1657 sit != mit2->second->reducedCostIdx.end(); sit++){
1658
1659 tmpSet.insert( inVars[ *sit ] );
1660 }
1661
1662 mit2->second->reducedCostIdx.clear();
1663
1664 for(sit = tmpSet.begin(); sit != tmpSet.end(); sit++){
1665
1666 //make sure that variable *sit is in the new reset master
1667
1668 if( inVars.find( *sit) != inVars.end() )
1669 mit2->second->reducedCostIdx.insert( *sit );
1670 }
1671 tmpSet.clear();
1672 //end reset reduced cost indexes
1673
1674 }//end loop over nodes in tree -- mit2
1675
1676 //reset the indexes of variables in the current integer incumbent
1677 for(vit = m_zOptIndexes.begin() ; vit != m_zOptIndexes.end(); vit++) *vit = inVars[ *vit ];
1678
1679
1680 //reset the indexes of variables in the root LP
1681 for(vit = m_zOptRootLP.begin() ; vit != m_zOptRootLP.end(); vit++) *vit = inVars[ *vit ];
1682
1684
1685 delete m_solver;
1686 m_osinstanceMaster = m_osrouteSolver->m_osinstanceMaster;
1687 m_solver = new CoinSolver();
1688
1689 //kipp -- later have clp be an option
1690 //I guess for now it must be an Osi solver
1691 m_solver->sSolverName ="cbc";
1692 //std::cout << m_osinstanceMaster->printModel( ) << std::endl;
1693 m_solver->osinstance = m_osrouteSolver->m_osinstanceMaster;
1694 m_solver->osoption = m_osoption;
1695
1696 m_solver->buildSolverInstance();
1697
1698 //get the solver interface
1699 m_si = m_solver->osiSolver;
1700
1701 if(m_si->getNumCols() != m_osrouteSolver->m_osinstanceMaster->getVariableNumber() )
1702 throw ErrorClass("there is an inconsistency in the the model rebuid in resetMaster");
1703
1704 std::cout << "OSINTANCE NUMBER OF COLUMNS = " << m_osrouteSolver->m_osinstanceMaster->getVariableNumber() << std::endl;
1705 std::cout << "OSINTANCE NUMBER OF ROWS = " << m_osrouteSolver->m_osinstanceMaster->getConstraintNumber() << std::endl;
1706 std::cout << "SOLVER INTERFACE NUMBER OF COLUMNS = " << m_si->getNumCols() << std::endl;
1707 std::cout << "SOLVER INTERFACE NUMBER OF ROWS = " <<m_si->getNumRows() << std::endl;
1708
1709
1710 //kipp this is a check, DO NOT do in production run
1711
1712
1713
1714 double lpVal;
1715
1716 for(mit2 = m_nodeMap.begin(); mit2 !=m_nodeMap.end(); mit2++){
1717
1718 lpVal = 0;
1719
1720 for(i = 0; i < mit2->second->thetaNumNonz; i++){
1721
1722 lpVal += m_si->getObjCoefficients()[ mit2->second->thetaIdx[ i] ]*mit2->second->theta[ i];
1723
1724
1725 }
1726
1727 if( ( lpVal - mit2->second->lpValue > m_osDecompParam.zeroTol ) ||
1728 (mit2->second->lpValue - lpVal > m_osDecompParam.zeroTol ) ) throw ErrorClass( "uh oh, problem with node lp value" );
1729
1730 //std::cout << "lpVal = " << lpVal << " lpValue = " << mit2->second->lpValue << std::endl ;
1731 }
1732
1733
1734 //end check
1735
1736 //m_si->writeLp( "gailTest2" );
1737
1738 delete[] tmpVals;
1739 tmpVals = NULL;
1740
1741 } catch (const ErrorClass& eclass) {
1742
1743 throw ErrorClass(eclass.errormsg);
1744
1745 }
1746
1747
1748}//end resetMaster
1749
1751
1752 std::map<int, OSNode*>::iterator mit;
1753 int i;
1754
1755 std::cout << std::endl << std::endl;
1756
1757 std::cout << "NUMBER OF REMAINING DANGLING NODES = " << m_nodeMap.size() << std::endl;
1758
1759 if( m_nodeMap.size() > 0) m_zLB = OSDBL_MAX; //find best LP value over dangling nodes
1760
1761 for ( mit = m_nodeMap.begin() ;
1762 mit != m_nodeMap.end(); mit++ ){
1763
1764 std::cout << "NODE ID VALUE = " << mit->second->nodeID << " " ;
1765 std::cout << " NODE LP VALUE = " << mit->second->lpValue << std::endl;
1766
1767 for(i = 0; i < mit->second->rowIdxNumNonz; i++){
1768
1769 std::cout << "CONSTRAINT = " << mit->second->rowIdx[ i] ;
1770 std::cout << " CONSTRAINT LB = " << mit->second->rowLB[ i] ;
1771 std::cout << " CONSTRAINT UB = " << mit->second->rowUB[ i] << std::endl;
1772 }
1773
1774 if( mit->second->lpValue < m_zLB) m_zLB = mit->second->lpValue;
1775
1776
1777 }
1778 m_nodeMap.clear();
1779
1780
1781}//end printTreeInfo
1782
1783
1784void OSColGenApp::checkNodeConsistency( const int rowIdx, const OSNode *osnode){
1785 try{
1786 if( osnode == NULL) return;
1787 //we are going to throw an exception if we try to add a constraint to a node that is already there
1788 std::set<int> indexSet;
1789 int i;
1790 int j;
1791 int rowIdxNumNonz = 0;
1792 int thetaNumNonz = 0;
1793 rowIdxNumNonz = osnode->rowIdxNumNonz;
1794 thetaNumNonz = osnode->thetaNumNonz;
1795 std::map<int, double> varSumMap;
1796
1797 std::cout << "MESSAGE: CHECKING FOR NODE CONSISTENCY CONSTRAINT" << std::endl;
1798
1799 for(i = 0; i < thetaNumNonz; i++){
1800
1801
1802 //loop over theta variables
1803 std::cout << "theta idx " << osnode->thetaIdx[ i] << " theta value " << osnode->theta[ i] << std::endl;
1804
1805 for(j = m_osrouteSolver->m_thetaPnt[ osnode->thetaIdx[ i] ]; j < m_osrouteSolver->m_thetaPnt[ osnode->thetaIdx[ i] + 1 ]; j++ ){
1806
1807 if( varSumMap.find( m_osrouteSolver->m_thetaIndex[ j] ) == varSumMap.end() ){
1808
1809 varSumMap[ m_osrouteSolver->m_thetaIndex[ j] ] = osnode->theta[ i];
1810
1811 }else{
1812
1813 varSumMap[ m_osrouteSolver->m_thetaIndex[ j] ] += osnode->theta[ i];
1814 }
1815
1816 std::cout << "xijk idx " << m_osrouteSolver->m_thetaIndex[ j] << " variable name = " <<
1817 m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[ j] ] << std::endl;
1818
1819 }
1820
1821 }
1822
1823
1824
1825 for(i = 0; i < rowIdxNumNonz; i++){
1826
1827 std::cout << " row number " << osnode->rowIdx[ i] << " LB = " << osnode->rowLB[ i] << " UB = "
1828 << osnode->rowUB[ i] ;
1829
1830 std::cout << " variable index = " << m_rowIdxVarMap[ osnode->rowIdx[ i] ] ;
1831
1832 std::cout << " variable name = " << m_osrouteSolver->m_variableNames[ m_rowIdxVarMap[ osnode->rowIdx[ i] ] ] ;
1833
1834 std::cout << " variable sum = " << varSumMap[ m_rowIdxVarMap[ osnode->rowIdx[ i] ]] << std::endl ;
1835
1836 if(indexSet.find( osnode->rowIdx[ i] ) == indexSet.end() ){
1837
1838 indexSet.insert( osnode->rowIdx[ i] );
1839
1840 }else{
1841
1842
1843 throw ErrorClass( "We are trying to add an existing constraint to a node" );
1844 }
1845
1846
1847
1848 }
1849
1850 } catch (const ErrorClass& eclass) {
1851
1852 throw ErrorClass(eclass.errormsg);
1853
1854 }
1855}//end checkNodeConsistency
OSOption * osoption
Implements a solve method for the Coin solvers.
used for throwing exceptions.
std::string errormsg
errormsg is the error that is causing the exception to be thrown
std::vector< int > m_zOptRootLP
m_zOptRootLP is the vector theta indexes corresponding to optimal LP solution at the roor tnode
void printDebugInfo()
OSOption * m_osoption
Definition OSColGenApp.h:54
void getColumns(const double *yA, const int numARows, const double *yB, const int numBRows, int &numNewColumns, int *&numNonz, double *&cost, int **&rowIdx, double **&values, double &lowerBound)
RETURN VALUES: int numNewColumns – number of new columns generated int* numNonz – number of nonzeros ...
~OSColGenApp()
Default destructor.
double m_zUB
m_zUB is the upper bound
OsiSolverInterface * m_si
Definition OSColGenApp.h:84
OSDecompFactoryInitializer * m_factoryInit
Definition OSColGenApp.h:43
std::vector< int > m_zOptIndexes
m_zOptIndexes is the vector theta indexes corresponding to the current m_zUB
double * m_yB
m_yB is the dual for the cuts that get added
double m_zRootLP
m_zRootLP is the value of the root LP relaxation
CoinSolver * m_solver
Definition OSColGenApp.h:82
std::map< int, int > inVars
OSDecompParam m_osDecompParam
Application specific parameters.
Definition OSColGenApp.h:88
std::string m_message
m_message is the message to the pauHana routine
double * m_yA
m_yA is the dual values for the initial restricted constraints
int m_numNodesGenerated
kount the nodes generated
Definition OSColGenApp.h:91
void createBranchingCut(const int *thetaIdx, const double *theta, const int numThetaVar, std::map< int, int > &varConMap, int &rowIdx)
INPUT: – sparse version int* thetaIdx – index vector of nonzero theta variables double* theta – the s...
OSDecompSolver * m_osrouteSolver
Definition OSColGenApp.h:46
int m_maxCols
m_maxCols is the maximum number of columns we can have
void solveRestrictedMasterRelaxation()
kipp – document
OSInstance * m_osinstanceMaster
Definition OSColGenApp.h:53
void getInitialRestrictedMaster()
std::vector< double > m_zRootLPx_vals
m_zRootLPx_vals holds root node optimal LP solution nonzero values
Definition OSColGenApp.h:60
std::map< int, int > m_rowIdxVarMap
map the variable generated at a node with a variable
Definition OSColGenApp.h:78
bool branchAndBound()
the method that invokes and controls branch and bound
void getOptions(OSOption *osoption)
int m_numColumnsOld
when m_numColumnsGenerated - m_numColumnsOld hits masterColumnResetValue we do a reset
void getCuts(const double *thetaVar, const int numThetaVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)
RETURN VALUES: int numNewRows – number of new rows generated int* numNonz – number of nonzeros in eac...
double m_zLB
m_zLB is the lower bound
std::map< int, OSNode * > m_nodeMap
nodeMap is the map of nodes in the branch and bound tree
bool m_calledBranchAndBound
this variable is true if we have called the branchAndBound() method
Definition OSColGenApp.h:73
int m_numColumnsGenerated
kount the columns generated
Definition OSColGenApp.h:94
bool isInteger(const double *thetaVar, const int numThetaVar, const double tol)
INPUT: double* thetaVar – the vector of primal master values int numThetaVar – size of master primal ...
OSColGenApp()
Default Constructor.
int m_maxRows
m_maxRows is the maximum number of rows we can have
void printTreeInfo()
OSNode * createChild(const OSNode *osnode, std::map< int, int > &varConMap, const int rowIdx, const double rowLB, const double rowUB)
INPUT: OSNode* osnode – the parent node for which we create a child std::map<int, int> varConMap – th...
void checkNodeConsistency(const int rowIdx, const OSNode *osnode)
double * m_theta
m_theta is the primal values in the master
static std::map< std::string, OSDecompSolverFactory * > factories
double * rowUB
rowUB is a vector of row upper bounds
Definition OSNode.h:50
int * thetaIdx
theta is an array of primal solution variable indexes
Definition OSNode.h:65
int parentID
parentID is the node ID of the parent
Definition OSNode.h:34
int rowIdxNumNonz
rowIdxNumNonz is the number of non-zero elements in rowIndex
Definition OSNode.h:42
int thetaNumNonz
thetaNumNonz is the number of non-zero elements in the theta variable solution at this node
Definition OSNode.h:60
int nodeID
nodeID is the node ID
Definition OSNode.h:39
int * rowIdx
rowIdx is a vector of row indexes for which we are setting the upper and lower bounds
Definition OSNode.h:47
double lpValue
lpValue is the LP relaxation for the node
Definition OSNode.h:56
double * theta
theta is an array of primal positive values this is used for branching and creating new children node...
Definition OSNode.h:71
std::set< int > reducedCostIdx
reducedCostVec will hold variables within a tolerance on their reduced costs.
Definition OSNode.h:93
double * rowLB
rowLB is a vector of row lower bounds
Definition OSNode.h:53
The Option Class.
Definition OSOption.h:3565
#define OSDBL_MAX