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
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;
91 std::cout << "FINISH SET FACTORY OPTION " << std::endl;
92 //m_osrouteSolver = new OSBearcatSolverXij( osoption);
93
94 //share the common parameters
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
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
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
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;
308
309 //pass options
312
313
314 //get the solver interface
316
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
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
432
433 //go into branch and bound
434 m_message = "";
435 std::cout << "START BRANCH AND BOUND = " << std::endl;
437
438 //demand values
439 //m_osrouteSolver->m_demand;
440
441 std::cout << "FINISH BRANCH AND BOUND = " << std::endl;
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
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 ****";
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
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
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) {
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
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
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] >=
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;
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;
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;
1695
1697
1698 //get the solver interface
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 = " <<
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.
virtual void buildSolverInstance()
The implementation of the corresponding virtual function.
OsiSolverInterface * osiSolver
osiSolver is the osi solver object – in this case clp, glpk, cbc, cplex, symphony or dylp
virtual void solve()
The implementation of the corresponding virtual function.
std::string sSolverName
sSolverName is the name of the Coin solver used, e.g.
OSInstance * osinstance
osinstance holds the problem instance in-memory as an OSInstance object
OSOption * osoption
osoption holds the solver options in-memory as an OSOption object
OSResult * osresult
osresult holds the solution or result of the model in-memory as an OSResult object
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
int masterColumnResetValue
when the number of columns in the master hits masterColumnResetValue do a column purge
double artVarCoeff
artVarCoeff is the coefficient of the artificial variable in the objective function
int columnLimit
columnLimit is the limit on the number of columns that can be generated in a single call to solveRest...
double optTolPerCent
we fathom a node if UB*(1 - optTolPerCent) <= LB
double zeroTol
we terminate column generation when the reduced costs are not smaller than zeroTol
int nodeLimit
nodeLimit is the limit on the number of nodes that are allowed in the branch and bound tree
static std::map< std::string, OSDecompSolverFactory * > factories
OSDecompParam m_osDecompParam
share the parameters with the decomposition solver
virtual void getCutsMultiCommod(const double *thetaVar, const int numThetaVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)=0
This is the routine that generates the multi-item cuts.
int m_maxBmatrixCon
m_maxBmatrixCon is the maximum number of B matrix constraints it is the number of tour breaking const...
int m_numNodes
m_numNodes is the number of nodes (both pickup and hub) in the model
int m_maxMasterRows
m_maxMasterColumns is the maximumn number of rows we allow in the master, in this application it is e...
int m_numBmatrixCon
m_numBmatrixCon is the number of constraints in B - 1, we have the -1 because: m_pntBmatrix[ k] point...
int m_numHubs
m_numHubs is the number of hubs/routes
std::string * m_variableNames
virtual void initializeDataStructures()=0
allocate memory and initialize arrays
std::set< std::pair< int, double > > intVarSet
intVarSet holds and std::pair where the first element is the index of an integer variable and the sec...
OSInstance * m_osinstanceMaster
virtual 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)=0
RETURN VALUES:
OSOption * m_osoption
int m_maxMasterColumns
m_maxMasterColumns is the maximumn number of columns we allow in the master
virtual void getCutsTheta(const double *thetaVar, const int numThetaVar, int &numNewRows, int *&numNonz, int **&colIdx, double **&values, double *&rowLB, double *&rowUB)=0
RETURN VALUES:
virtual void resetMaster(std::map< int, int > &inVars, OsiSolverInterface *si)=0
INPUT:
virtual void getBranchingCut(const double *thetaVar, const int numThetaVar, const std::map< int, int > &varConMap, int &varIdx, int &numNonz, int *&indexes, double *&values)=0
Dense Version.
virtual void pauHana(std::vector< int > &m_zOptIndexes, std::vector< double > &m_zRootLPx_vals, int numNodes, int numColsGen, std::string message)=0
virtual OSInstance * getInitialRestrictedMaster()=0
int getConstraintNumber()
Get number of constraints.
int getVariableNumber()
Get number of variables.
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
std::vector< SolverOption * > getSolverOptions(std::string solver_name)
Get the options associated with a given solver.
std::string getSolutionStatusType(int solIdx)
Get the [i]th optimization solution status type, where i equals the given solution index.
#define OSDBL_MAX