00001
00014 #include "OSBearcatSolverXij.h"
00015
00016 #include "OSErrorClass.h"
00017 #include "OSDataStructures.h"
00018
00019 #include "OSInstance.h"
00020 #include "OSCoinSolver.h"
00021 #include "OSConfig.h"
00022 #include "OSResult.h"
00023
00024 #include "OSiLReader.h"
00025 #include "OSiLWriter.h"
00026 #include "OSoLReader.h"
00027 #include "OSoLWriter.h"
00028 #include "OSrLReader.h"
00029 #include "OSrLWriter.h"
00030 #include "OSInstance.h"
00031 #include "OSFileUtil.h"
00032 #include "OSStringUtil.h"
00033
00034 #include "CoinTime.hpp"
00035
00036 #include "ClpFactorization.hpp"
00037 #include "ClpNetworkMatrix.hpp"
00038 #include "OsiClpSolverInterface.hpp"
00039
00040 #include <iostream>
00041
00042 #include <algorithm>
00043
00044
00045 #ifdef HAVE_CMATH
00046 # include <cmath>
00047 #else
00048 # ifdef HAVE_MATH_H
00049 # include <math.h>
00050 # else
00051 # error "don't have header file for math"
00052 # endif
00053 #endif
00054
00055 #ifdef HAVE_CTIME
00056 # include <ctime>
00057 #else
00058 # ifdef HAVE_TIME_H
00059 # include <time.h>
00060 # else
00061 # error "don't have header file for time"
00062 # endif
00063 #endif
00064
00065 using std::ostringstream;
00066
00067
00068 OSBearcatSolverXij::OSBearcatSolverXij() {
00069 std::cout << "INSIDE OSBearcatSolverXij CONSTRUCTOR with OSOption argument" << std::endl;
00070 }
00071
00072 OSBearcatSolverXij::OSBearcatSolverXij(OSOption *osoption) {
00073 std::cout << "INSIDE OSBearcatSolverXij CONSTRUCTOR with OSOption argument" << std::endl;
00074
00075 m_hubPoint = NULL;
00076
00077 m_bestIPValue = OSDBL_MAX;
00078 m_bestLPValue = -OSDBL_MAX;
00079
00080 m_upperBoundL = NULL;
00081 m_upperBoundLMax = -OSINT_MAX;
00082 m_minDemand = OSINT_MAX;
00083
00084 m_u = NULL;
00085 m_v = NULL;
00086 m_g = NULL;
00087 m_px = NULL;
00088 m_tx =NULL;
00089 m_varIdx = NULL;
00090
00091 m_optL = NULL;
00092 m_optD = NULL;
00093 m_vv = NULL;
00094 m_vvpnt = NULL;
00095
00096 m_demand = NULL;
00097 m_cost = NULL;
00098
00099 m_rc = NULL;
00100
00101 m_routeCapacity = NULL;
00102 m_routeMinPickup = NULL;
00103
00104 m_osoption = osoption;
00105
00106 m_multiCommodCutLimit = 250;
00107 m_numMultCuts = 0;
00108
00109
00110 }
00111
00112 void OSBearcatSolverXij::initializeDataStructures(){
00113
00114 int k;
00115 int i;
00116 int l;
00117 int tmpVal;
00118
00119 try{
00120
00121
00122
00123 getOptions( m_osoption);
00124
00125 if(m_cost == NULL) throw ErrorClass("Option file did not contain cost data");
00126
00127 getVariableIndexMap();
00128
00129
00130 m_maxMasterRows = m_maxBmatrixCon + m_numNodes;
00131
00132 m_upperBoundL = new int[ m_numHubs];
00133 m_lowerBoundL = new int[ m_numHubs];
00134
00135 m_hubPoint = new int[ m_numHubs];
00136
00137 permuteHubs();
00138
00139
00140 for(k = 0; k < m_numHubs; k++){
00141
00142
00143 tmpVal = 0;
00144 for(i = 0; i < m_numHubs; i++)
00145 if( i != k) tmpVal += m_routeCapacity[ i];
00146
00147 m_lowerBoundL[ k] = std::max( m_routeMinPickup[ k], (m_totalDemand - tmpVal) ) ;
00148
00149 m_upperBoundL[ k] = m_routeCapacity[ k];
00150
00151
00152 if(m_upperBoundL[ k] > m_totalDemand) m_upperBoundL[ k] = m_totalDemand;
00153 if( m_upperBoundL[ k] > m_upperBoundLMax) m_upperBoundLMax = m_upperBoundL[ k];
00154
00155 }
00156
00157
00158
00159
00160 m_varIdx = new int[ m_upperBoundLMax + 1];
00161
00162
00163 m_u = new double*[ m_numNodes];
00164 m_v = new double*[ m_numNodes];
00165 m_g = new double*[ m_numNodes];
00166
00167 m_px = new int*[ m_numNodes];
00168 m_tx = new int*[ m_numNodes];
00169
00170
00171
00181 for (i = 0; i < m_numNodes; i++) {
00182
00183
00184 m_u[ i] = new double[ m_upperBoundLMax + 1];
00185 m_v[ i] = new double[ m_upperBoundLMax + 1];
00186
00187
00188
00189 for(l = 0; l <= m_upperBoundLMax; l++){
00190
00191 m_u[ i][l] = OSDBL_MAX;
00192 m_v[ i][l] = OSDBL_MAX;
00193 }
00194
00195 m_g[ i] = new double[ m_numNodes];
00196 m_px[ i] = new int[ m_upperBoundLMax + 1];
00197 m_tx[ i] = new int[m_upperBoundLMax + 1];
00198
00199
00200 }
00201
00202
00203
00204 m_optL = new int[ m_numHubs];
00205 m_optD = new int[ m_numHubs];
00206
00207 m_vv = new double*[ m_numHubs];
00208 m_vvpnt = new int*[ m_numHubs];
00209 m_rc = new double*[ m_numHubs];
00210
00211 for (k = 0; k < m_numHubs; k++) {
00212
00213
00214 m_vv[ k] = new double[ m_totalDemand + 1];
00215 m_vvpnt[ k] = new int[ m_totalDemand + 1];
00216
00217
00218
00219
00220 m_rc[ k] = new double[ m_upperBoundL[ k]*(m_numNodes*m_numNodes - m_numNodes)];
00221
00222
00223 }
00224
00225 m_optValHub = new double[ m_numHubs];
00226
00227 m_variableNames = new std::string[ m_numNodes*(m_numNodes - 1)];
00228
00229 createVariableNames();
00230
00231
00232
00233 m_pntAmatrix = new int[ m_numNodes - m_numHubs + 1];
00234
00235
00236
00237
00238
00239 m_Amatrix = new int[ (m_numNodes - m_numHubs)*(m_numNodes - 1) ];
00240 createAmatrix();
00241
00242
00243
00244 int numVar = m_numNodes*m_numNodes - m_numNodes ;
00245 m_tmpScatterArray = new int[ numVar ];
00246 for(i = 0; i < numVar; i++){
00247
00248 m_tmpScatterArray[ i] = 0;
00249
00250 }
00251
00252
00253 m_newColumnNonz = new int[ m_numHubs] ;
00254 m_costVec = new double[ m_numHubs];
00255 m_newColumnRowIdx = new int*[ m_numHubs];
00256 m_newColumnRowValue = new double*[ m_numHubs];
00257
00258 for (k = 0; k < m_numHubs; k++) {
00259
00260
00261 m_newColumnRowValue[ k] = new double[ m_maxBmatrixCon + m_numNodes];
00262 m_newColumnRowIdx[ k] = new int[ m_maxBmatrixCon + m_numNodes];
00263
00264 }
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276 m_newRowNonz = new int[ 100] ;
00277 m_newRowColumnIdx = new int*[ 100] ;
00278 m_newRowColumnValue = new double*[ 100] ;
00279 m_newRowUB = new double[ 100];
00280 m_newRowLB = new double[ 100];
00281
00282
00283
00284
00285 for (k = 0; k < 100; k++) {
00286
00287
00288 m_newRowColumnValue[ k] = new double[ m_maxMasterColumns];
00289 m_newRowColumnIdx[ k] = new int[ m_maxMasterColumns];
00290
00291 }
00292
00293
00294 m_convexityRowIndex = new int[ m_maxMasterColumns];
00295
00296
00297 branchCutIndexes = new int[ m_maxMasterColumns];
00298 branchCutValues = new double[ m_maxMasterColumns];
00299
00300 m_thetaPnt = new int[ m_maxMasterColumns + 1];
00301 for(i = 0; i <= m_maxMasterColumns; i++){
00302 m_thetaPnt[ i] = 0;
00303 }
00304 m_thetaCost = new double[ m_maxMasterColumns];
00305 m_thetaIndex = new int[ m_maxThetaNonz];
00306 m_numThetaVar = 0;
00307 m_numThetaNonz = 0;
00308 m_thetaPnt[ m_numThetaVar] = 0;
00309
00310
00311
00312
00313
00314 m_pntBmatrix = new int[ m_maxBmatrixCon];
00315
00316 m_BmatrixIdx = new int[ m_maxBmatrixNonz];
00317 m_BmatrixRowIndex = new int[ m_maxBmatrixCon];
00318 for(i = 0; i < m_maxBmatrixCon; i++) m_BmatrixRowIndex[ i] = -1;
00319
00320
00321 m_BmatrixVal = new double[ m_maxBmatrixNonz];
00322
00323
00324 m_numBmatrixCon = 0;
00325 m_numBmatrixNonz = 0;
00326 m_pntBmatrix[ m_numBmatrixCon] = 0;
00327
00328 ;
00329
00330
00331 m_separationIndexMap = new int[m_numNodes*(m_numNodes - 1)];
00332
00333 for(i = 0; i < m_numNodes*(m_numNodes - 1); i++){
00334
00335 m_separationIndexMap[ i] = OSINT_MAX;
00336
00337 }
00338
00339
00340
00341
00342
00343 getSeparationInstance();
00344 for(k = 0; k < m_numHubs; k++){
00345
00346 m_multCommodCutSolvers.push_back( getMultiCommodInstance( k) );
00347
00348 }
00349
00350 } catch (const ErrorClass& eclass) {
00351
00352 throw ErrorClass(eclass.errormsg);
00353
00354 }
00355
00356
00357 }
00358
00359
00360 OSBearcatSolverXij::~OSBearcatSolverXij(){
00361
00362 std::cout << "INSIDE ~OSBearcatSolverXij DESTRUCTOR" << std::endl;
00363
00364
00365
00366
00367 int i;
00368
00369 delete[] m_routeCapacity;
00370 m_routeCapacity = NULL;
00371
00372
00373 delete[] m_routeMinPickup;
00374 m_routeMinPickup = NULL;
00375
00376 for(i = 0; i < m_numNodes; i++){
00377
00378
00379
00380 delete[] m_v[i];
00381 delete[] m_g[i];
00382 delete[] m_px[i];
00383 delete[] m_tx[i];
00384 delete[] m_u[i];
00385
00386
00387 }
00388
00389 delete[] m_u;
00390 m_u= NULL;
00391
00392 delete[] m_v;
00393 m_v= NULL;
00394
00395 delete[] m_g;
00396 m_g= NULL;
00397
00398 delete[] m_px;
00399 m_px= NULL;
00400
00401 delete[] m_tx;
00402 m_tx= NULL;
00403
00404
00405
00406 if(m_demand != NULL){
00407
00408 delete[] m_demand;
00409 }
00410
00411
00412 if(m_varIdx != NULL){
00413 delete[] m_varIdx;
00414 m_varIdx = NULL;
00415
00416 }
00417
00418 if(m_cost != NULL ){
00419 delete[] m_cost;
00420 m_cost = NULL;
00421 }
00422
00423 for(i = 0; i < m_numHubs; i++){
00424
00425 delete[] m_vv[i];
00426 delete[] m_vvpnt[i];
00427 delete[] m_rc[ i];
00428
00429
00430 }
00431 delete[] m_optL;
00432 m_optL = NULL;
00433 delete[] m_optD;
00434 m_optD = NULL;
00435 delete[] m_vv;
00436 m_vv = NULL;
00437 delete[] m_vvpnt;
00438 m_vvpnt = NULL;
00439
00440 delete[] m_rc;
00441 m_rc = NULL;
00442
00443 delete[] m_upperBoundL;
00444 m_upperBoundL = NULL;
00445
00446 delete[] m_lowerBoundL;
00447 m_lowerBoundL = NULL;
00448
00449 delete[] m_hubPoint;
00450 m_hubPoint = NULL;
00451
00452 delete[] m_optValHub;
00453 m_optValHub = NULL;
00454
00455 delete[] m_variableNames;
00456 m_variableNames = NULL;
00457
00458 delete[] m_pntAmatrix;
00459 m_pntAmatrix = NULL;
00460
00461 delete[] m_Amatrix;
00462 m_Amatrix = NULL;
00463
00464 delete[] m_tmpScatterArray;
00465 m_tmpScatterArray = NULL;
00466
00467 delete[] m_newColumnNonz ;
00468 m_newColumnNonz = NULL;
00469 delete[] m_costVec ;
00470 m_costVec = NULL;
00471
00472 for(i = 0; i < m_numHubs; i++){
00473
00474 delete[] m_newColumnRowIdx[i];
00475 delete[] m_newColumnRowValue[i];
00476 }
00477
00478 delete[] m_newColumnRowIdx;
00479 m_newColumnRowIdx = NULL;
00480
00481 delete[] m_newColumnRowValue;
00482 m_newColumnRowValue = NULL;
00483
00484
00485 delete[] m_convexityRowIndex;
00486 m_convexityRowIndex = NULL;
00487
00488
00489
00490 delete[] m_newRowNonz;
00491 m_newRowNonz = NULL;
00492
00493 delete[] m_newRowUB ;
00494 m_newRowUB = NULL;
00495
00496 delete[] m_newRowLB ;
00497 m_newRowLB = NULL;
00498
00499
00500 for (i = 0; i < m_numHubs; i++) {
00501
00502 delete[] m_newRowColumnValue[ i];
00503 delete[] m_newRowColumnIdx[ i];
00504
00505 }
00506
00507 delete[] m_newRowColumnIdx;
00508 m_newRowColumnIdx = NULL;
00509
00510 delete[] m_newRowColumnValue;
00511 m_newRowColumnValue = NULL;
00512
00513
00514 delete[] branchCutIndexes ;
00515 branchCutIndexes = NULL;
00516
00517 delete[] branchCutValues ;
00518 branchCutValues = NULL;
00519
00520
00521 delete[] m_thetaPnt;
00522 m_thetaPnt = NULL;
00523
00524 delete[] m_thetaIndex;
00525 m_thetaIndex = NULL;
00526
00527
00528 delete[] m_thetaCost;
00529 m_thetaCost = NULL;
00530
00531
00532 delete[] m_pntBmatrix ;
00533 m_pntBmatrix = NULL;
00534
00535 delete[] m_BmatrixIdx ;
00536 m_BmatrixIdx = NULL;
00537
00538 delete[] m_BmatrixVal ;
00539 m_BmatrixVal = NULL;
00540
00541 delete[] m_BmatrixRowIndex ;
00542 m_BmatrixRowIndex = NULL;
00543
00544
00545
00546
00547 delete[] m_separationIndexMap;
00548 m_separationIndexMap = NULL;
00549
00550 delete m_separationClpModel;
00551 m_separationClpModel = NULL;
00552
00553 delete m_osinstanceSeparation;
00554 m_osinstanceSeparation = NULL;
00555
00556 std::vector<CoinSolver*>::iterator vit;
00557 for(vit = m_multCommodCutSolvers.begin(); vit < m_multCommodCutSolvers.end(); vit++){
00558
00559 delete *vit;
00560 }
00561
00562 }
00563
00564
00565
00566
00567
00568
00569
00570 void OSBearcatSolverXij::getOptL( double** c) {
00571
00572
00573
00574 int k;
00575 int d;
00576 int d1;
00577 int kountVar;
00578 double testVal;
00579 int l;
00580
00581 double trueMin;
00582
00583 bool isFeasible;
00584 isFeasible = false;
00585
00586 kountVar = 0;
00587
00588
00589
00590 m_vv[ m_hubPoint[0] ][ 0] = 0;
00591 for(d = 1; d <= m_totalDemand; d++){
00592
00593 m_vv[ m_hubPoint[0] ][ d] = OSDBL_MAX;
00594
00595 }
00596
00597 for(k = 1; k < m_numHubs - 1; k++){
00598 for(d = 0; d <= m_totalDemand; d++){
00599
00600 m_vv[ m_hubPoint[ k] ][ d] = OSDBL_MAX;
00601
00602 }
00603 }
00604
00605
00606
00607
00608 int dlower;
00609 dlower = 0;
00610
00611 for(k = 1; k < m_numHubs; k++){
00612
00613 dlower += m_lowerBoundL[ m_hubPoint[k - 1] ];
00614
00615
00616 for(d = dlower; d <= m_totalDemand; d++){
00617
00618 m_vv[ m_hubPoint[ k] ][ d] = OSDBL_MAX;
00619
00620
00621
00622 for(d1 = 0; d1 <= d; d1++){
00623
00624 l = d - d1;
00625
00626
00627
00628 if( (m_vv[ m_hubPoint[ k - 1] ][ d1] < OSDBL_MAX) && (l <= m_upperBoundL[ m_hubPoint[ k - 1] ]) && (l >= m_lowerBoundL[ m_hubPoint[k - 1] ]) ){
00629
00630
00631
00632
00633
00634
00635 testVal = qrouteCost( m_hubPoint[ k - 1], l, c[ m_hubPoint[ k - 1] ], &kountVar);
00636
00637
00638
00639
00640 if( m_vv[ m_hubPoint[ k-1] ][ d1] + testVal < m_vv[ m_hubPoint[ k] ][ d] ){
00641
00642 m_vv[ m_hubPoint[ k] ][ d] = m_vv[ m_hubPoint[ k-1] ][ d1] + testVal;
00643
00644 m_vvpnt[ m_hubPoint[ k] ][ d] = d1;
00645
00646 }
00647
00648
00649 }
00650
00651 }
00652
00653 }
00654
00655
00656
00657 }
00658
00659 trueMin = OSDBL_MAX;
00660
00661
00662
00663
00664
00665
00666
00667
00668 for(d = dlower; d < m_totalDemand; d++){
00669
00670
00671 l = m_totalDemand - d;
00672
00673 if(m_vv[ m_hubPoint[ m_numHubs - 1] ][ d] < OSDBL_MAX && l <= m_upperBoundL[ m_hubPoint[ m_numHubs - 1] ] && l >= m_lowerBoundL[ m_hubPoint[ m_numHubs - 1] ]){
00674
00675
00676
00677
00678
00679 isFeasible = true;
00680
00681
00682 testVal = qrouteCost(m_hubPoint[m_numHubs -1] , l, c[ m_hubPoint[ m_numHubs -1] ], &kountVar);
00683
00684
00685
00686
00687 if(m_vv[ m_hubPoint[ m_numHubs - 1] ][ d] + testVal < trueMin){
00688
00689 trueMin = m_vv[ m_hubPoint[ m_numHubs -1] ][ d] + testVal;
00690 m_optD[ m_hubPoint[ m_numHubs -1] ] = d;
00691 m_optL[ m_hubPoint[ m_numHubs -1] ] = l;
00692
00693 }
00694
00695
00696 }
00697 }
00698
00699
00700
00701 if( isFeasible == false){
00702
00703 std::cout << "NOT ENOUGH CAPACITY " << std::endl;
00704 for(k = 0; k < m_numHubs; k++) std::cout << " k perm = " << m_hubPoint[ k ]<< std::endl;
00705 throw ErrorClass( "NOT ENOUGH CAPACITY ");
00706 }
00707
00708 k = m_numHubs - 1;
00709
00710 while( k - 1 >= 0) {
00711
00712 m_optD[ m_hubPoint[ k - 1] ] = m_vvpnt[ m_hubPoint[ k] ][ m_optD[ m_hubPoint[ k] ] ];
00713
00714 m_optL[ m_hubPoint[ k - 1] ] = m_optD[ m_hubPoint[ k ] ] - m_optD[ m_hubPoint[ k - 1] ] ;
00715
00716 k--;
00717
00718 }
00719
00720 }
00721
00722
00723
00724
00725
00726
00727 double OSBearcatSolverXij::qrouteCost(const int& k, const int& l, const double* c, int* kountVar){
00728
00729
00730
00731
00732 double rcost;
00733 rcost = OSDBL_MAX;
00734
00735
00736
00737 if(l <= 0){
00738
00739 std::cout << "LVALUE NEGATIVE OR ZERO " << l << std::endl;
00740 exit( 1);
00741 }
00742
00743
00744
00745 if(l > m_upperBoundL[ k]){
00746
00747 std::cout << "LVALUE BIGGER THAN UPPER BOUND " << l << std::endl;
00748 exit( 1);
00749 }
00750
00751
00752
00753
00754
00755 int startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
00756 c+= startPnt ;
00757
00758
00759
00760 *kountVar = 0;
00761 int bestLastNode;
00762
00763 bestLastNode = OSINT_MAX;
00764 int i;
00765 int j;
00766 int l1;
00767 int l2;
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782 for(i = m_numHubs; i < m_numNodes; i++){
00783
00784
00785 for(l1 = m_minDemand; l1 <= l; l1++){
00786
00787
00788 m_u[i][l1] = OSDBL_MAX;
00789
00790 m_v[i][l1] = OSDBL_MAX;
00791 m_px[i][l1] = -1;
00792
00793
00794 if(l1 == *(m_demand + i) ){
00795
00796 m_px[i][l1] = k;
00797
00798
00799 m_u[i][l1] = *(c + k*(m_numNodes - 1) + i - 1);
00800
00801
00802
00803 }
00804
00805 }
00806 }
00807
00808
00809
00810
00811
00812
00813
00814 if(l == m_minDemand){
00815
00816 for(i = m_numHubs; i < m_numNodes; i++){
00817
00818
00819 if( m_u[i][l] + *(c + i*(m_numNodes-1) + k ) < rcost){
00820
00821 rcost = m_u[i][l] + *(c + i*(m_numNodes-1) + k );
00822
00823
00824
00825
00826
00827 bestLastNode = i;
00828 }
00829
00830 }
00831
00832
00833
00834 *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
00835 *(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1) + bestLastNode - 1;
00836
00837
00838 return rcost;
00839 }
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850 for(l2 = m_minDemand + 1; l2 <= l; l2++){
00851
00852 for(i = m_numHubs; i < m_numNodes; i++) {
00853
00854
00855 if( *(m_demand + i) < l2 ){
00856
00857 for(j = m_numHubs; j < i; j++){
00858
00859
00860
00861
00862 l1 = l2 - *(m_demand + i);
00863
00864 if( m_px[j][ l1 ] != i ){
00865
00866
00867 m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i - 1) ;
00868
00869
00870
00871
00872 }else{
00873
00874 m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i - 1) ;
00875
00876
00877
00878 }
00879
00880
00881
00882 if(m_g[j][i] < m_u[i][l2] ){
00883
00884 m_u[i][l2] = m_g[j][i];
00885 m_px[i][l2] = j;
00886
00887 }
00888
00889
00890
00891 }
00892
00893
00894 for(j = i + 1; j < m_numNodes; j++){
00895
00896
00897
00898 l1 = l2 - *(m_demand + i);
00899
00900 if( m_px[j][ l1 ] != i ){
00901
00902
00903 m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i ) ;
00904
00905
00906 }else{
00907
00908 m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i ) ;
00909
00910 }
00911
00912
00913
00914 if(m_g[j][i] < m_u[i][l2] ){
00915
00916 m_u[i][l2] = m_g[j][i];
00917 m_px[i][l2] = j;
00918
00919 }
00920
00921
00922 }
00923
00924
00925
00926
00927
00928 for(j =m_numHubs; j < m_numNodes; j++){
00929
00930 if(j != i){
00931
00932 if( (m_g[j][i] < m_v[i][l2] ) && (m_px[i][l2] != j) ){
00933
00934
00935
00936 m_v[i][l2] = m_g[j][i];
00937 m_tx[i][l2] = j;
00938
00939
00940 }
00941
00942 }
00943
00944
00945 }
00946
00947
00948 if(l2 == l ){
00949
00950 if( m_u[i][l2] + *(c + i*(m_numNodes-1) + k ) < rcost){
00951
00952 rcost = m_u[i][l2] + *(c + i*(m_numNodes-1) + k );
00953
00954 bestLastNode = i;
00955 }
00956
00957 }
00958
00959
00960 }
00961
00962
00963 }
00964
00965
00966 }
00967
00968
00969
00970
00971
00972
00973
00974
00975 int currentNode;
00976 int successorNode;
00977 int lvalue;
00978
00979
00980
00981
00982
00983
00984
00985 *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1) + k ;
00986
00987
00988
00989
00990 if( bestLastNode == OSINT_MAX) return OSDBL_MAX;
00991
00992
00993 successorNode = k;
00994 currentNode = bestLastNode;
00995
00996 lvalue = l ;
00997
00998
00999
01000
01001 while(currentNode != k){
01002
01003 if(*kountVar == m_upperBoundLMax + 1) return OSDBL_MAX;
01004 if( m_px[ currentNode][ lvalue ] != successorNode){
01005
01006
01007
01008
01009 successorNode = currentNode;
01010 currentNode = m_px[ currentNode][ lvalue ];
01011
01012
01013 if(currentNode - successorNode > 0){
01014
01015
01016
01017
01018 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
01019
01020
01021 }else{
01022
01023
01024
01025
01026 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
01027
01028 }
01029
01030
01031 }else{
01032
01033
01034
01035 successorNode = currentNode;
01036 currentNode = m_tx[ currentNode][ lvalue ];
01037
01038 if(currentNode - successorNode > 0){
01039
01040
01041
01042 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode;
01043
01044 }else{
01045
01046
01047 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1) + successorNode - 1 ;
01048
01049 }
01050
01051 }
01052
01053
01054
01055
01056 lvalue = lvalue - *(m_demand + successorNode);
01057
01058
01059
01060 }
01061
01062
01063
01064
01065
01066
01067 return rcost;
01068
01069 }
01070
01071
01072
01073
01074 void OSBearcatSolverXij::getColumns(const double* yA, const int numARows,
01075 const double* yB, const int numBRows,
01076 int &numNewColumns, int* &numNonzVec, double* &costVec,
01077 int** &rowIdxVec, double** &valuesVec, double &lowerBound)
01078 {
01079
01080
01081
01082 int i;
01083 int j;
01084 int numCoulpingConstraints;
01085 numCoulpingConstraints = m_numNodes - m_numHubs;
01086
01087 int numVar;
01088 numVar = m_numNodes*m_numNodes - m_numHubs;
01089 int numNonz;
01090
01091 try{
01092
01093
01094
01095 if(numARows != m_numNodes) throw ErrorClass("inconsistent row count in getColumns");
01096
01097
01098
01099
01100 calcReducedCost( yA, yB );
01101
01102 int kountVar;
01103 double testVal;
01104 testVal = 0;
01105 int k;
01106 int startPntInc;
01107 int rowCount;
01108 double rowValue;
01109
01111
01112 double cpuTime;
01113 double start = CoinCpuTime();
01114 getOptL( m_rc);
01115 cpuTime = CoinCpuTime() - start;
01116 std::cout << "DYNAMIC PROGRAMMING CPU TIME " << cpuTime << std::endl;
01117 m_lowerBnd = 0.0;
01118 for(k = 0; k < m_numHubs; k++){
01119
01120
01121
01122 startPntInc = (m_optL[ k] - 1)*(m_numNodes*m_numNodes - m_numNodes);
01123
01124 std::cout << " whichBlock = " << k << " L = " << m_optL[ k] << std::endl;
01125
01126 testVal += m_optL[ k];
01127
01128 kountVar = 0;
01129
01131 m_optValHub[ k] = qrouteCost(k, m_optL[ k], m_rc[ k], &kountVar);
01132
01133 m_optValHub[ k] -= yA[ k + numCoulpingConstraints];
01134
01135 std::cout << "Best Reduced Cost Hub " << k << " = " << m_optValHub[ k] << std::endl;
01136 m_lowerBnd += m_optValHub[ k];
01137
01138
01139
01140
01141
01142 m_costVec[ k] = 0.0;
01143
01144 for(j = 0; j < kountVar; j++){
01145
01146
01147
01148
01149 m_tmpScatterArray[ m_varIdx[ j] - startPntInc ] += 1;
01150
01151
01152
01153
01154 m_costVec[ k] += m_cost[ m_varIdx[ j] - startPntInc ];
01155
01156 }
01157
01158
01159
01160 numNonz = 0;
01161
01162 for(i = 0; i < numCoulpingConstraints; i++){
01163
01164 rowCount = 0;
01165
01166 for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
01167
01168
01169
01170
01171 rowCount += m_tmpScatterArray[ m_Amatrix[ j] ];
01172
01173
01174 }
01175
01176 if(rowCount > 0){
01177
01178
01179 m_newColumnRowIdx[ k][ numNonz] = i;
01180 m_newColumnRowValue[ k][ numNonz] = rowCount;
01181 numNonz++;
01182 }
01183
01184
01185 }
01186
01187
01188 m_newColumnRowIdx[ k][ numNonz] = m_numNodes - m_numHubs + k;
01189 m_newColumnRowValue[ k][ numNonz++] = 1.0;
01190
01191
01192
01193
01194
01195 for(i = 0; i < m_numBmatrixCon; i++){
01196
01197 if(m_BmatrixRowIndex[ i] == -1 || m_BmatrixRowIndex[ i] == k ){
01198
01199
01200 rowValue = 0;
01201
01202 for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
01203
01204
01205
01206
01207
01208
01209
01210
01211 rowValue += m_tmpScatterArray[ m_BmatrixIdx[ j] ]*m_BmatrixVal[ j];
01212
01213
01214 }
01215
01216 if( rowValue > m_osDecompParam.zeroTol || rowValue < -m_osDecompParam.zeroTol){
01217
01218
01219 m_newColumnRowIdx[ k][ numNonz] = i + m_numNodes;
01220 m_newColumnRowValue[ k][ numNonz++] = rowValue;
01221
01222 }
01223
01224 }
01225
01226
01227 }
01228
01229 m_newColumnNonz[ k] = numNonz;
01230
01231
01232
01233
01234
01235
01236 for(j = 0; j < kountVar; j++){
01237
01238
01239
01240 m_thetaIndex[ m_numThetaNonz++ ] = m_varIdx[ j] - startPntInc ;
01241 m_tmpScatterArray[ m_varIdx[ j] - startPntInc ] = 0;
01242
01243
01244
01245 }
01246
01247
01248
01249
01250 intVarSet.insert ( std::pair<int,double>( m_numThetaVar, 1.0) );
01251 m_convexityRowIndex[ m_numThetaVar] = k;
01252 m_costVec[ k] = m_optL[ k]*m_costVec[ k];
01253 m_thetaCost[ m_numThetaVar++ ] = m_costVec[ k];
01254 m_thetaPnt[ m_numThetaVar ] = m_numThetaNonz;
01255
01256
01257
01258
01259
01260
01261
01294
01295
01296
01297
01298 }
01299
01300
01301
01302 numNonzVec = m_newColumnNonz;
01303 costVec = m_costVec;
01304 rowIdxVec = m_newColumnRowIdx;
01305 valuesVec = m_newColumnRowValue;
01306 std::cout << "Lower Bound = " << m_lowerBnd << std::endl;
01307
01308
01309 if(testVal != m_totalDemand) {
01310
01311 std::cout << "TOTAL DEMAND = " << m_totalDemand << std::endl;
01312 std::cout << "Test Value = " << testVal << std::endl;
01313 throw ErrorClass( "inconsistent demand calculation" );
01314 }
01315
01316
01317
01318
01319
01320
01321 } catch (const ErrorClass& eclass) {
01322
01323 throw ErrorClass(eclass.errormsg);
01324
01325 }
01326
01327
01328
01329 numNewColumns = m_numHubs;
01330 lowerBound = m_lowerBnd;
01331
01332 std::cout << "LEAVING GET COLUMNS" << std::endl;
01333 return;
01334
01335 }
01336
01337
01338
01632 OSInstance* OSBearcatSolverXij::getInitialRestrictedMaster( ){
01633
01634 getInitialSolution();
01635
01636 std::cout << "Executing OSBearcatSolverXij::getInitialRestrictedMaster( )" << std::endl;
01637
01638
01639 int numVarArt;
01640
01641 numVarArt = m_numNodes;
01642
01643 int numXVar;
01644 numXVar = m_numNodes*m_numNodes - m_numNodes;
01645
01646 double* xVar = NULL;
01647 xVar = new double[ numXVar];
01648
01649 int i;
01650 int k;
01651 std::string testFileName;
01652 std::string osil;
01653
01654 std::map<int, std::map<int, std::vector<int> > >::iterator mit;
01655 std::map<int, std::vector<int> >::iterator mit2;
01656 std::vector<int>::iterator vit;
01657
01658 m_osinstanceMaster = NULL;
01659
01660
01661
01662 int kountNonz;
01663 int kount;
01664 m_numberOfSolutions = 1;
01665 int numThetaVar = m_numberOfSolutions*m_numHubs;
01666
01667 double *values = new double[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar + numVarArt];
01668 int *indexes = new int[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar + numVarArt] ;
01669 int *starts = new int[ numThetaVar + 1 + numVarArt];
01670 kount = 0;
01671 starts[ 0] = 0;
01672 int startsIdx;
01673 startsIdx = 0;
01674 startsIdx++;
01675
01676 double tspObjValue;
01677 double demandSum;
01678
01679 SparseVector *objcoeff = NULL;
01680
01681 for(i = 0; i < numVarArt; i++){
01682 m_convexityRowIndex[ m_numThetaVar] = -1;
01683 m_thetaPnt[ m_numThetaVar++] = 0;
01684
01685 }
01686 try {
01687
01688
01689 m_osinstanceMaster = new OSInstance();
01690 m_osinstanceMaster->setInstanceDescription("The Initial Restricted Master");
01691
01692
01693
01694 m_osinstanceMaster->setVariableNumber( m_numberOfSolutions*m_numHubs + numVarArt);
01695
01696
01697 m_osinstanceMaster->setObjectiveNumber( 1);
01698
01699 objcoeff = new SparseVector( m_numberOfSolutions*m_numHubs + numVarArt);
01700
01701
01702 m_osinstanceMaster->setConstraintNumber( m_numNodes);
01703
01704
01705
01706 int varNumber;
01707 varNumber = 0;
01708 std::string masterVarName;
01709 kountNonz = 0;
01710
01711 for(i = 0; i < m_numNodes; i++){
01712
01713 objcoeff->indexes[ varNumber ] = varNumber ;
01714
01715
01716
01717
01718
01719 objcoeff->values[ varNumber ] = m_osDecompParam.artVarCoeff;
01720
01721 m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt("AP", i ) ,
01722 0, 1.0, 'C');
01723
01724 values[ kountNonz] = 1;
01725 indexes[ kountNonz++] = i ;
01726 starts[ startsIdx++] = kountNonz;
01727
01728 }
01729
01730 m_bestIPValue = 0;
01731 mit = m_initSolMap.find( 0);
01732 if(mit == m_initSolMap.end() ) throw ErrorClass("Problem finding first solution in m_initSolMap");
01733
01734 CoinSolver *solver = NULL;
01735 solver = getTSP(m_numNodes, m_cost);
01736
01737
01738 for(k = 0; k < m_numHubs; k++){
01739
01740 mit2 = mit->second.find( k);
01741 if(mit2 == mit->second.end() ) throw ErrorClass("Problem finding hub k in the solution map");
01742
01743 tspObjValue = getRouteDistance(m_numNodes, mit2->first, solver,
01744 mit2->second, xVar);
01745 demandSum = 0;
01746
01747 for ( vit = mit2->second.begin() ; vit != mit2->second.end(); vit++ ) {
01748
01749 demandSum += m_demand[ *vit];
01750 values[ kountNonz] = 1;
01751 indexes[ kountNonz++] = *vit - m_numHubs ;
01752
01753 }
01754
01755
01756
01757 for(i = 0; i < numXVar; i++){
01758
01759
01760 if(xVar[ i] > .1) {
01761 m_thetaIndex[ m_numThetaNonz++ ] = i;
01762
01763 }
01764 }
01765
01766
01767
01768 values[ kountNonz] = 1;
01769 indexes[ kountNonz++] = m_numNodes - m_numHubs + k ;
01770
01771
01772
01773
01774 m_convexityRowIndex[ m_numThetaVar] = k;
01775 m_thetaCost[ m_numThetaVar++ ] = demandSum*tspObjValue;
01776 m_thetaPnt[ m_numThetaVar ] = m_numThetaNonz;
01777
01778 masterVarName = makeStringFromInt("theta(", k);
01779 masterVarName += makeStringFromInt(",", 0);
01780 masterVarName += ")";
01781 intVarSet.insert ( std::pair<int,double>(varNumber, 1.0) );
01782 m_osinstanceMaster->addVariable(varNumber++, masterVarName, 0, 1, 'C');
01783
01784 std::cout << "Optimal Objective Value = " << demandSum*tspObjValue << std::endl;
01785
01786
01787 objcoeff->indexes[ k + numVarArt ] = k + numVarArt ;
01788 objcoeff->values[ k + numVarArt ] = demandSum*tspObjValue;
01789
01790 m_bestIPValue += demandSum*tspObjValue;
01791
01792 std::cout << "m_bestIPValue = " << m_bestIPValue << std::endl;
01793 starts[ startsIdx++] = kountNonz;
01794
01795 }
01796 std::cout << " m_numThetaVar " << m_numThetaVar << std::endl;
01797
01798
01799
01800 for( i = 0; i < m_numNodes - m_numHubs ; i++){
01801
01802 m_osinstanceMaster->addConstraint(i, makeStringFromInt("visitNode_", i + m_numHubs) , 1.0, 1.0, 0);
01803 }
01804
01805 kount = 0;
01806
01807
01808 for( i = m_numNodes - m_numHubs; i < m_numNodes ; i++){
01809
01810 m_osinstanceMaster->addConstraint(i, makeStringFromInt("convexityRowRoute_", kount++ ) , 1.0, 1.0, 0);
01811 }
01812
01813 m_osinstanceMaster->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
01814
01815
01816
01817 m_osinstanceMaster->setLinearConstraintCoefficients(kountNonz , true,
01818 values, 0, kountNonz - 1, indexes, 0, kountNonz - 1, starts, 0, startsIdx);
01819
01820
01821 std::cout << m_osinstanceMaster->printModel( ) << std::endl;
01822 std::cout << "NONZ = " << kountNonz << std::endl;
01823
01824
01825 delete objcoeff;
01826 objcoeff = NULL;
01827 delete[] xVar;
01828 xVar = NULL;
01829 delete solver->osinstance;
01830 delete solver;
01831
01832
01833
01834
01835 } catch (const ErrorClass& eclass) {
01836 std::cout << std::endl << std::endl << std::endl;
01837
01838 if(objcoeff != NULL) delete objcoeff;
01839 if(xVar != NULL) delete xVar;
01840
01841 throw ErrorClass(eclass.errormsg);
01842 }
01843
01844 return m_osinstanceMaster;
01845 }
01846
01847
01848
01849 void OSBearcatSolverXij::getOptions(OSOption *osoption) {
01850
01851
01852 std::cout << "Executing getOptions(OSOption *osoption)" << std::endl;
01853
01854 try{
01855
01856 int i;
01857
01858
01859 std::vector<SolverOption*> solverOptions;
01860 std::vector<SolverOption*>::iterator vit;
01861 std::vector<int>::iterator vit2;
01862 std::vector<int >demand;
01863 std::vector<std::string >nodeName;
01864 std::vector<int >routeCapacity;
01865 std::vector<int >routeMinPickup;
01866 std::vector<double >arcCost;
01867 std::vector<double >::iterator vit3;
01868 std::vector<std::string>::iterator vit4;
01869
01870 m_numberOfSolutions = 0;
01871 solverOptions = osoption->getSolverOptions("routeSolver");
01872 if (solverOptions.size() == 0) throw ErrorClass( "options for routeSolver not available");
01873
01874
01875 int tmpVal;
01876 double tmpCost;
01877
01878
01879 std::string routeString;
01880 std::string solutionString;
01881 std::string::size_type pos1;
01882 std::string::size_type pos2;
01883 std::string::size_type pos3;
01884
01885
01886 std::map<int, std::map<int, std::vector<int> > >::iterator mit;
01887 std::map<int, std::vector<int> >::iterator mit2;
01888 int solutionNumber;
01889 int routeNumber;
01890
01891
01892 for (vit = solverOptions.begin(); vit != solverOptions.end(); vit++) {
01893
01894
01895
01896
01897
01898
01899
01900 if( (*vit)->name.find("numHubs") != std::string::npos){
01901
01902
01903 std::istringstream hubBuffer( (*vit)->value);
01904 hubBuffer >> m_numHubs;
01905 std::cout << "numHubs = " << m_numHubs << std::endl;
01906
01907 }else{
01908
01909 if((*vit)->name.find("numNodes") != std::string::npos){
01910
01911
01912 std::istringstream numNodesBuffer( (*vit)->value);
01913 numNodesBuffer >> m_numNodes;
01914 std::cout << "numNodes = " << m_numNodes << std::endl;
01915
01916 }else{
01917 if((*vit)->name.find("totalDemand") != std::string::npos){
01918
01919
01920 std::istringstream totalDemandBuffer( (*vit)->value);
01921 totalDemandBuffer >> m_totalDemand;
01922 std::cout << "m_totalDemand = " << m_totalDemand << std::endl;
01923
01924 }else{
01925 if((*vit)->name.find("routeMinPickup") != std::string::npos){
01926
01927
01928 std::istringstream routeMinPickupBuffer( (*vit)->value);
01929 routeMinPickupBuffer >> tmpVal;
01930 routeMinPickup.push_back( tmpVal);
01931
01932
01933 }else{
01934 if( (*vit)->name.find("demand") != std::string::npos ){
01935
01936
01937 std::istringstream demandBuffer( (*vit)->value);
01938 demandBuffer >> tmpVal;
01939 if(tmpVal <= 0 && demand.size() > (unsigned int) m_numHubs) throw ErrorClass("must have strictly positive demand");
01940 if(tmpVal < m_minDemand && demand.size() > (unsigned int) m_numHubs ) m_minDemand = tmpVal;
01941 demand.push_back( tmpVal);
01942 nodeName.push_back( (*vit)->description);
01943
01944
01945 }else{
01946 if((*vit)->name.find("routeCapacity") != std::string::npos ){
01947 std::istringstream routeCapacityBuffer( (*vit)->value);
01948 routeCapacityBuffer >> tmpVal;
01949 routeCapacity.push_back( tmpVal);
01950 std::cout << "m_routeCapacity = " << tmpVal << std::endl;
01951
01952 }else{
01953
01954 if((*vit)->name.find("osilFile") != std::string::npos ){
01955 m_initOSiLFile = (*vit)->value;
01956 std::cout << "m_initOSiLFile = " << m_initOSiLFile << std::endl;
01957
01958 }else{
01959
01960 if( (*vit)->name.find("restrictedMasterSolution") != std::string::npos ){
01961
01962
01963
01964
01965
01966
01967 pos1 = (*vit)->category.find( ":");
01968 if(pos1 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
01969
01970
01971 solutionString = (*vit)->category.substr( 0, pos1);
01972 routeString = (*vit)->category.substr( pos1 + 1);
01973
01974 pos2 = solutionString.find( "n");
01975 if(pos2 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
01976
01977 std::istringstream solutionBuffer( solutionString.substr( pos2 + 1) );
01978 solutionBuffer >> solutionNumber;
01979
01980
01981
01982 pos3 = routeString.find( "e");
01983 if(pos3 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
01984 std::istringstream routeBuffer( routeString.substr( pos3 + 1) );
01985 routeBuffer >> routeNumber;
01986
01987 std::istringstream nodeBuffer( (*vit)->value);
01988 nodeBuffer >> tmpVal;
01989
01990 mit = m_initSolMap.find( solutionNumber );
01991
01992 if( mit != m_initSolMap.end() ){
01993
01994
01995
01996 mit2 = mit->second.find( routeNumber);
01997
01998 if(mit2 != mit->second.end() ){
01999
02000
02001
02002 mit2->second.push_back( tmpVal);
02003
02004
02005 }else{
02006
02007
02008 std::vector<int> tmpVec;
02009 tmpVec.push_back( tmpVal) ;
02010 mit->second.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
02011
02012
02013 }
02014
02015 }else{
02016
02017 std::vector<int> tmpVec;
02018 tmpVec.push_back( tmpVal) ;
02019
02020 std::map<int, std::vector<int> > tmpMap;
02021 tmpMap.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
02022 m_initSolMap.insert( std::pair<int, std::map<int, std::vector<int> > >(solutionNumber, tmpMap) ) ;
02023
02024 }
02025 }
02026 else{
02027 if( (*vit)->name.find("maxMasterColumns") != std::string::npos){
02028
02029
02030 std::istringstream maxMasterColumns( (*vit)->value);
02031 maxMasterColumns >> m_maxMasterColumns;
02032 std::cout << "m_maxMasterColumn = " << m_maxMasterColumns << std::endl;
02033
02034 }else{
02035 if( (*vit)->name.find("maxThetaNonz") != std::string::npos){
02036
02037 std::istringstream maxThetaNonz( (*vit)->value);
02038 maxThetaNonz >> m_maxThetaNonz;
02039 std::cout << "m_maxThetaNonz = " << m_maxThetaNonz << std::endl;
02040
02041 }else{
02042 if( (*vit)->name.find("use1OPTstart") != std::string::npos){
02043 m_use1OPTstart = false;
02044 if ( (*vit)->value.find("true") != std::string::npos ) m_use1OPTstart = true;
02045 std::cout << "m_use1OPTstart = " << m_use1OPTstart << std::endl;
02046
02047 }else{
02048 if( (*vit)->name.find("maxBmatrixCon") != std::string::npos ){
02049
02050 std::istringstream maxBmatrixCon( (*vit)->value);
02051 maxBmatrixCon >> m_maxBmatrixCon;
02052 std::cout << "m_maxBmatrixCon = " << m_maxBmatrixCon << std::endl;
02053
02054 }else{
02055 if( (*vit)->name.find("maxBmatrixNonz") != std::string::npos ){
02056
02057 std::istringstream maxBmatrixNonz( (*vit)->value);
02058 maxBmatrixNonz >> m_maxBmatrixNonz;
02059 std::cout << "m_maxBmatrixNonz = " << m_maxBmatrixNonz << std::endl;
02060
02061
02062 }else{
02063
02064 if( (*vit)->name.find("cost") != std::string::npos ){
02065
02066
02067 std::istringstream costBuffer( (*vit)->value);
02068 costBuffer >> tmpCost;
02069 arcCost.push_back( tmpCost);
02070
02071 }
02072 }
02073 }
02074 }
02075 }
02076 }
02077 }
02078 }
02079 }
02080 }
02081 }
02082 }
02083 }
02084 }
02085
02086 }
02087
02088
02089
02090 i = 0;
02091 m_routeCapacity = new int[ m_numHubs];
02092 if( (unsigned int)m_numHubs != routeCapacity.size( ) ) throw ErrorClass("inconsistent number of HUBS");
02093 for (vit2 = routeCapacity.begin(); vit2 != routeCapacity.end(); vit2++) {
02094
02095 *(m_routeCapacity + i++) = *vit2;
02096
02097 std::cout << "ROUTE CAP = " << *vit2 << std::endl;
02098
02099 }
02100 routeCapacity.clear();
02101
02102
02103
02104 i = 0;
02105 m_routeMinPickup = new int[ m_numHubs];
02106 if( (unsigned int)m_numHubs != routeMinPickup.size( ) ) throw ErrorClass("inconsistent number of HUBS");
02107
02108 for(int k = 0; k < m_numHubs; k++){
02109 m_routeMinPickup[ k] = 1;
02110 }
02111 for (vit2 = routeMinPickup.begin(); vit2 != routeMinPickup.end(); vit2++) {
02112
02113 *(m_routeMinPickup + i++) = *vit2;
02114
02115 }
02116 routeMinPickup.clear();
02117
02118
02119
02120 i = 0;
02121 m_demand = new int[ m_numNodes];
02122 if( (unsigned int)m_numNodes != demand.size( ) )
02123 throw ErrorClass("inconsistent number of demand nodes");
02124 for (vit2 = demand.begin(); vit2 != demand.end(); vit2++) {
02125
02126 *(m_demand + i++) = *vit2;
02127
02128 }
02129 demand.clear();
02130
02131
02132 i = 0;
02133 m_nodeName = new std::string[ m_numNodes];
02134
02135 for (vit4 = nodeName.begin(); vit4 != nodeName.end(); vit4++) {
02136
02137 *(m_nodeName + i++) = *vit4;
02138
02139 }
02140 nodeName.clear();
02141
02142
02143
02144
02145 m_cost = NULL;
02146 m_costSetInOption = false;
02147 if(arcCost.size() != (unsigned int)(m_numNodes*m_numNodes - m_numNodes) )
02148 throw ErrorClass("input cost vector not consistent with number of nodes");
02149 if(arcCost.size() >= 1){
02150 m_costSetInOption = true;
02151 i = 0;
02152 m_cost = new double[ m_numNodes*m_numNodes - m_numNodes ];
02153 for (vit3 = arcCost.begin(); vit3 != arcCost.end(); vit3++) {
02154
02155 *(m_cost + i++) = *vit3;
02156
02157 }
02158 arcCost.clear();
02159 }
02160
02161
02162 m_numberOfSolutions = m_initSolMap.size();
02163
02164
02165 } catch (const ErrorClass& eclass) {
02166
02167 throw ErrorClass(eclass.errormsg);
02168
02169 }
02170
02171 }
02172
02173
02174
02175 void OSBearcatSolverXij::getCutsTheta(const double* theta, const int numTheta,
02176 int &numNewRows, int* &numNonz, int** &colIdx,
02177 double** &values, double* &rowLB, double* &rowUB) {
02178
02179
02180
02181
02182 int i;
02183 int j;
02184 int k;
02185 int index;
02186 int rowKount;
02187 int tmpKount;
02188 int indexAdjust = m_numNodes - m_numHubs;
02189 double* tmpRhs;
02190 int numSepRows = m_osinstanceSeparation->getConstraintNumber() ;
02191
02192 tmpRhs = new double[ numSepRows ];
02193
02194 for(i = 0; i < numSepRows; i++){
02195
02196 tmpRhs[ i] = 0;
02197 }
02198
02199 try{
02200 m_osinstanceSeparation->bConstraintsModified = true;
02201
02202 if(numTheta != m_numThetaVar ) throw
02203 ErrorClass("number of master varibles in OSBearcatSolverXij::getCuts inconsistent");
02204
02205
02206
02207
02208
02209
02210
02211
02212 for(i = 0; i < numTheta; i++){
02213
02214
02215 if(theta[ i] > m_osDecompParam.zeroTol){
02216
02217
02218 for(j = m_thetaPnt[ i ]; j < m_thetaPnt[ i + 1 ]; j++ ){
02219
02220
02221
02222
02223
02224 rowKount = m_separationIndexMap[ m_thetaIndex[ j] ];
02225
02226
02227
02228 if(rowKount < OSINT_MAX ){
02229
02230 tmpRhs[ rowKount] -= theta[ i];
02231
02232 }
02233
02234 }
02235 }
02236 }
02237
02238
02239
02240
02241 for(i = indexAdjust; i < numSepRows - 1; i++){
02242
02243 if(-tmpRhs[ i] > 1 + m_osDecompParam.zeroTol ){
02244
02245
02246
02247
02248 int tmpKount = indexAdjust;
02249 for(int i1 = m_numHubs; i1 < m_numNodes; i1++){
02250
02251
02252
02253 for(int j1 = i1+1; j1 < m_numNodes; j1++){
02254
02255 if(tmpKount == i){
02256
02257
02258
02259
02260
02261 m_BmatrixVal[ m_numBmatrixNonz ] = 1 ;
02262
02263 m_BmatrixIdx[ m_numBmatrixNonz++ ] = i1*(m_numNodes - 1) + j1 - 1 ;
02264
02265 m_BmatrixVal[ m_numBmatrixNonz ] = 1 ;
02266
02267 m_BmatrixIdx[ m_numBmatrixNonz++ ] = j1*(m_numNodes - 1) + i1 ;
02268 m_numBmatrixCon++;
02269 m_pntBmatrix[ m_numBmatrixCon ] = m_numBmatrixNonz;
02270
02271 numNewRows = 1;
02272
02273 m_newRowNonz[ 0] = 0;
02274 m_newRowUB[ 0] = 1;
02275 m_newRowLB[ 0] = 0;
02276
02277
02278
02279
02280 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
02281 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
02282
02283
02284 std::cout << " m_BmatrixIdx[ j] " << m_BmatrixIdx[ j] << std::endl ;
02285
02286 m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 1;
02287
02288 }
02289
02290
02291
02292
02293 for(k = 0; k < m_numThetaVar ; k++){
02294
02295
02296 tmpKount = 0;
02297
02298
02299 for(j = m_thetaPnt[k]; j < m_thetaPnt[k + 1] ; j++){
02300
02301 if(m_tmpScatterArray[ m_thetaIndex[ j] ] > 0 ){
02302
02303 std::cout << " Variable " << m_variableNames[ m_thetaIndex[ j] ] << std::endl;
02304
02305 tmpKount++;
02306
02307 }
02308
02309 }
02310
02311 if(tmpKount > 0){
02312
02313
02314 m_newRowColumnIdx[0][ m_newRowNonz[ 0] ] = k ;
02315
02316 m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
02317
02318
02319 }
02320
02321 }
02322
02323
02324
02325
02326
02327
02328
02329 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
02330 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
02331
02332 m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 0;
02333
02334 }
02335
02336 numNonz = m_newRowNonz;
02337 colIdx = m_newRowColumnIdx;
02338 values = m_newRowColumnValue;
02339 rowUB = m_newRowUB;
02340 rowLB = m_newRowLB;
02341
02342 delete[] tmpRhs;
02343 tmpRhs = NULL;
02344
02345
02346
02347 m_numThetaVar++;
02348 m_convexityRowIndex[ m_numThetaVar] = -1;
02349 m_thetaPnt[ m_numThetaVar] = m_numThetaNonz;
02350
02351
02352
02353 return;
02354
02355 }
02356
02357 tmpKount++;
02358
02359 }
02360
02361 }
02362
02363
02364 }
02365
02366 m_separationClpModel->setRowUpper(i, tmpRhs[ i] );
02367 m_separationClpModel->setRowLower(i, tmpRhs[ i] );
02368
02369 }
02370
02371
02372
02373
02374
02375 std::vector<int> dualIdx;
02376 std::vector<int>::iterator vit1;
02377 std::vector<int>::iterator vit2;
02378
02379
02380
02381
02382 for(k = 0; k < indexAdjust; k++){
02383
02384
02385
02386 m_separationClpModel->setRowUpper(k, 0.0);
02387
02388
02389 m_separationClpModel->setLogLevel( 0);
02390
02391 m_separationClpModel->primal();
02392
02393 if(m_separationClpModel->getObjValue() > m_osDecompParam.zeroTol){
02394 std::cout << "DOING SEPARATION FOR NODE " << k + m_numHubs << std::endl;
02395 std::cout << "SEPERATION OBJ VALUE = " << m_separationClpModel->getObjValue() << std::endl;
02396 numNewRows = 1;
02397
02398 for(i = 0; i < m_numNodes - m_numHubs ; i++){
02399
02400 if( m_separationClpModel->getRowPrice()[ i] - m_osDecompParam.zeroTol <= -1) dualIdx.push_back( i) ;
02401 }
02402
02403 for (vit1 = dualIdx.begin(); vit1 != dualIdx.end(); vit1++) {
02404
02405 i = *vit1 + m_numHubs;
02406
02407 for (vit2 = dualIdx.begin(); vit2 != dualIdx.end(); vit2++) {
02408
02409 j = *vit2 + m_numHubs;
02410
02411 if( i > j ){
02412
02413 index = i*(m_numNodes -1) + j;
02414 std::cout << "CUT VARIABLE = " << m_variableNames[ index ] <<std::endl;
02415 m_BmatrixVal[ m_numBmatrixNonz ] = 1 ;
02416 m_BmatrixIdx[ m_numBmatrixNonz++ ] = index ;
02417
02418 }else{
02419
02420 if( i < j ){
02421
02422 index = i*(m_numNodes -1) + j - 1;
02423 std::cout << "CUT VARIABLE = " << m_variableNames[ index ] <<std::endl;
02424 m_BmatrixVal[ m_numBmatrixNonz ] = 1 ;
02425 m_BmatrixIdx[ m_numBmatrixNonz++ ] = index ;
02426
02427 }
02428 }
02429
02430 }
02431 }
02432
02433
02434 m_numBmatrixCon++;
02435 m_pntBmatrix[ m_numBmatrixCon ] = m_numBmatrixNonz;
02436
02437
02438
02439
02440
02441
02442 for(i = indexAdjust; i < numSepRows - 1; i++){
02443
02444 m_separationClpModel->setRowUpper(i, 0.0 );
02445 m_separationClpModel->setRowLower(i, 0.0 );
02446
02447
02448 }
02449 m_separationClpModel->setRowUpper(k, 1.0);
02450 delete[] tmpRhs;
02451 tmpRhs = NULL;
02452
02453
02454 m_newRowNonz[ 0] = 0;
02455 m_newRowUB[ 0] = dualIdx.size() - 1;
02456 m_newRowLB[ 0] = 0;
02457
02458 dualIdx.clear();
02459
02460
02461
02462
02463 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
02464 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
02465
02466 m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 1;
02467
02468 }
02469
02470
02471
02472
02473 for(i = 0; i < m_numThetaVar ; i++){
02474
02475
02476 tmpKount = 0;
02477 for(j = m_thetaPnt[i]; j < m_thetaPnt[i + 1] ; j++){
02478
02479 if(m_tmpScatterArray[ m_thetaIndex[ j] ] > 0 ){
02480
02481 tmpKount++;
02482
02483 }
02484
02485 }
02486
02487 if(tmpKount > 0){
02488
02489
02490 m_newRowColumnIdx[0][ m_newRowNonz[ 0] ] = i ;
02491
02492 m_newRowColumnValue[0][ m_newRowNonz[ 0]++ ] = tmpKount;
02493
02494
02495 }
02496
02497 }
02498
02499
02500
02501
02502 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
02503 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
02504
02505 m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 0;
02506
02507 }
02508
02509
02510
02511 numNonz = m_newRowNonz;
02512 colIdx = m_newRowColumnIdx;
02513 values = m_newRowColumnValue;
02514 rowUB = m_newRowUB;
02515 rowLB = m_newRowLB;
02516 m_numThetaVar++;
02517 m_convexityRowIndex[ m_numThetaVar] = -1;
02518 m_thetaPnt[ m_numThetaVar] = m_numThetaNonz;
02519
02520
02521
02522 return;
02523
02524
02525
02526 }
02527 m_separationClpModel->setRowUpper(k, 1.0);
02528 dualIdx.clear();
02529
02530 }
02531
02532
02533
02534
02535 for(i = indexAdjust; i < numSepRows - 1; i++){
02536
02537 m_separationClpModel->setRowUpper(i, 0.0 );
02538 m_separationClpModel->setRowLower(i, 0.0 );
02539
02540
02541 }
02542
02543 delete[] tmpRhs;
02544 tmpRhs = NULL;
02545
02546 } catch (const ErrorClass& eclass) {
02547
02548 throw ErrorClass(eclass.errormsg);
02549
02550 }
02551
02552
02553
02554 }
02555
02556
02557 void OSBearcatSolverXij::getCutsMultiCommod(const double* theta, const int numTheta,
02558 int &numNewRows, int* &numNonz, int** &colIdx,
02559 double** &values, double* &rowLB, double* &rowUB) {
02560
02561
02562
02563
02564 numNewRows = 0;
02565
02566
02567
02568 int i;
02569 int j;
02570 int k;
02571 int ivalue;
02572 int jvalue;
02573 int thetaIdx;
02574 int inodenum;
02575 int jnodenum;
02576 int j1;
02577 int j2;
02578 int kount;
02579 double wVal;
02580 double uVal;
02581 bool foundCut;
02582 double objVal;
02583
02584 int ckijIndex;
02585 int numNonHubs;
02586 numNonHubs = m_numNodes - m_numHubs;
02587
02588 int numVar;
02589 double rowValue;
02590
02591 double* scatterValues;
02592 int numXijVar = m_numNodes*m_numNodes - m_numNodes;
02593 scatterValues = new double[ numXijVar ];
02594 for(i = 0; i < numXijVar; i++ )scatterValues[ i] = 0;
02595
02596 double tmpCoeff;
02597
02598 double *wcoeff = NULL;
02599 wcoeff = new double[ numNonHubs];
02600 CoinSolver *solver;
02601
02602 std::cout << std::endl << std::endl;
02603 std::cout << "INSIDE getCutsMultiCommod " << std::endl;
02604 std::cout << std::endl << std::endl;
02605
02606 std::map<int, std::vector<std::pair<int,double> > > xVarMap;
02607 std::map<int, std::vector<std::pair<int,double> > >::iterator mit;
02608 std::vector<std::pair<int,double> >::iterator vit;
02609
02610 std::map<std::pair<int, int>,int>xVarIndexMap;
02611 std::pair<int,int> indexPair;
02612 ostringstream cutString;
02613
02614
02615
02616
02617
02618
02619
02647
02648
02649 try{
02650
02651 for(i = 0; i < numTheta; i++){
02652 xVarMap[ m_convexityRowIndex[ i] ] ;
02653
02654 if(theta[ i] > m_osDecompParam.zeroTol){
02655
02656
02657 for(j = m_thetaPnt[ i ]; j < m_thetaPnt[ i + 1 ]; j++ ){
02658
02659 mit = xVarMap.find( m_convexityRowIndex[ i]) ;
02660
02661 if(mit != xVarMap.end() ){
02662
02663 mit->second.push_back( std::pair<int, double>( m_thetaIndex[ j], theta[ i]) );
02664
02665 }
02666 }
02667 }
02668 }
02669
02670
02671
02672
02673 for(k = 0; k < m_numHubs; k++){
02674
02675
02676 mit = xVarMap.find( k) ;
02677
02678 solver = m_multCommodCutSolvers[ k];
02679
02680 numVar = solver->osiSolver->getNumCols();
02681 for(i = 0; i < numVar; i++ ) solver->osiSolver->setObjCoeff( i, 0);
02682
02683 for(i = 0; i < numNonHubs; i++) wcoeff[ i ] = 0;
02684
02685 if(mit != xVarMap.end() ){
02686
02687 std::cout << "CONVEXITY ROW " << " = " << mit->first << std::endl;
02688
02689 for(vit = mit->second.begin(); vit < mit->second.end(); vit++){
02690
02691
02692
02693 ivalue = vit->first /(m_numNodes - 1);
02694
02695 jvalue = vit->first - ivalue*(m_numNodes - 1);
02696
02697 if( jvalue >= ivalue ){
02698
02699
02700 inodenum = ivalue;
02701 jnodenum = jvalue + 1;
02702
02703
02704 }else{
02705
02706
02707
02708 inodenum = ivalue;
02709 jnodenum = jvalue;
02710 }
02711
02712 if(jnodenum != k){
02713
02714
02715 wcoeff[ jnodenum - m_numHubs ] += vit->second;
02716
02717 if( inodenum == k ){
02718 ckijIndex = jnodenum - m_numHubs;
02719 } else{
02720
02721 inodenum -= m_numHubs;
02722 jnodenum -= m_numHubs;
02723
02724
02725 if( jnodenum > inodenum) ckijIndex = inodenum*(numNonHubs - 1) + jnodenum - 1 ;
02726 else ckijIndex = inodenum*(numNonHubs - 1) + jnodenum ;
02727
02728
02729 ckijIndex += numNonHubs ;
02730
02731 }
02732
02733
02734 ckijIndex += numNonHubs;
02735
02736 tmpCoeff = solver->osiSolver->getObjCoefficients()[ ckijIndex] ;
02737
02738
02739 tmpCoeff = tmpCoeff - vit->second;
02740 if( ckijIndex > numVar - 1) throw ErrorClass( "problem with ckijIndex calculation");
02741
02742 solver->osiSolver->setObjCoeff( ckijIndex, tmpCoeff );
02743 }
02744
02745 }
02746
02747 }
02748 foundCut = false;
02749
02750 for(i = 0; i < numNonHubs; i++){
02751
02752
02753 solver->osiSolver->setObjCoeff( i, wcoeff[ i ] );
02754
02755
02756 solver->solve();
02757
02758
02759
02760 if(solver->osiSolver->getObjValue() > .1){
02761 objVal = 0;
02762 std::cout << "Separation Obj Value = " <<
02763 solver->osiSolver->getObjValue() << " Node " << m_nodeName[i + m_numHubs] << std::endl;
02764
02765 foundCut = true;
02766 m_newRowNonz[ numNewRows ] = 0;
02767
02768
02769 cutString.str("");
02770 cutString << "";
02771
02772
02773 kount = numNonHubs;
02774
02775 wVal = solver->osiSolver->getColSolution()[ i];
02776 objVal += wVal*solver->osiSolver->getObjCoefficients()[ i];
02777
02778 if(wVal < m_osDecompParam.zeroTol )throw ErrorClass("problem with wVal in generating a multicommodity cut");
02779
02780
02781 for(j2 = m_numHubs; j2 < m_numNodes; j2++){
02782
02783 indexPair.first = k;
02784 indexPair.second = j2;
02785
02786 uVal = solver->osiSolver->getColSolution()[ kount];
02787 objVal += uVal*solver->osiSolver->getObjCoefficients()[ kount];
02788
02789
02790
02791
02792 if (j2 == (i + m_numHubs) ){
02793
02794 if( (wVal - uVal) > m_osDecompParam.zeroTol || (uVal - wVal) > m_osDecompParam.zeroTol ){
02795
02796 cutString << " =";
02797 cutString << (wVal - uVal);
02798 cutString << "*" ;
02799
02800
02801
02802
02803
02804 m_BmatrixVal[ m_numBmatrixNonz ] = (wVal - uVal)/wVal ;
02805
02806
02807 if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
02808 std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
02809 throw ErrorClass( "index mapping problem in generating multicommodity cut");
02810 }else{
02811 m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
02812 cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
02813 }
02814
02815 }
02816
02817 }else{
02818
02819 if( (-uVal > m_osDecompParam.zeroTol) || (uVal > m_osDecompParam.zeroTol) ){
02820
02821 cutString << " ";
02822 cutString << - uVal;
02823 cutString << "*" ;
02824
02825
02826
02827
02828
02829
02830 m_BmatrixVal[ m_numBmatrixNonz ] = -uVal/wVal ;
02831
02832
02833
02834 if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
02835 std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
02836 throw ErrorClass( "index mapping problem in generating multicommodity cut");
02837 }else{
02838 m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
02839 cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
02840 }
02841 }
02842
02843 }
02844
02845 kount++;
02846 }
02847
02848
02849 for(j1 = m_numHubs; j1 < m_numNodes; j1++){
02850
02851
02852
02853
02854 for(j2 = m_numHubs; j2 < j1; j2++){
02855
02856 indexPair.first = j1;
02857 indexPair.second = j2;
02858
02859 uVal = solver->osiSolver->getColSolution()[ kount];
02860 objVal += uVal*solver->osiSolver->getObjCoefficients()[ kount];
02861
02862
02863
02864 if (j2 == (i + m_numHubs) ){
02865
02866
02867 if( (wVal - uVal) > m_osDecompParam.zeroTol || (uVal - wVal) > m_osDecompParam.zeroTol ){
02868
02869 cutString << " +";
02870 cutString << (wVal - uVal) ;
02871 cutString << "*" ;
02872
02873
02874
02875
02876 m_BmatrixVal[ m_numBmatrixNonz ] = (wVal - uVal)/wVal ;
02877
02878
02879
02880 if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
02881 std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
02882 throw ErrorClass( "index mapping problem in generating multicommodity cut");
02883 }else{
02884 m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
02885 cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
02886 }
02887
02888 }
02889
02890 }else{
02891
02892 if( (-uVal > m_osDecompParam.zeroTol) || (uVal > m_osDecompParam.zeroTol) ){
02893
02894 cutString << " ";
02895 cutString << - uVal;
02896 cutString << "*" ;
02897
02898
02899
02900
02901
02902 m_BmatrixVal[ m_numBmatrixNonz ] = -uVal/wVal ;
02903
02904
02905
02906 if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
02907 std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
02908 throw ErrorClass( "index mapping problem in generating multicommodity cut");
02909 }else{
02910 m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
02911 cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
02912 }
02913 }
02914 }
02915 kount++;
02916 }
02917
02918
02919 for(j2 = j1 + 1; j2 < m_numNodes; j2++){
02920
02921 indexPair.first = j1;
02922 indexPair.second = j2;
02923
02924 uVal = solver->osiSolver->getColSolution()[ kount];
02925 objVal += uVal*solver->osiSolver->getObjCoefficients()[ kount];
02926
02927
02928
02929 if (j2 == (i + m_numHubs) ){
02930
02931 if( (wVal - uVal) > m_osDecompParam.zeroTol || (uVal - wVal) > m_osDecompParam.zeroTol ){
02932
02933 cutString << " +";
02934 cutString << (wVal - uVal);
02935 cutString << "*" ;
02936
02937
02938
02939
02940 m_BmatrixVal[ m_numBmatrixNonz ] = (wVal - uVal)/wVal ;
02941
02942
02943
02944 if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
02945 std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
02946 throw ErrorClass( "index mapping problem in generating multicommodity cut");
02947 }else{
02948 m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
02949 cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
02950 }
02951 }
02952
02953 }else{
02954
02955 if( (-uVal > m_osDecompParam.zeroTol) || (uVal > m_osDecompParam.zeroTol) ){
02956
02957 cutString << " ";
02958 cutString << - uVal;
02959 cutString << "*" ;
02960
02961
02962
02963
02964 m_BmatrixVal[ m_numBmatrixNonz ] = -uVal/wVal ;
02965
02966
02967 if(m_xVarIndexMap.find( indexPair) == m_xVarIndexMap.end() ){
02968 std::cout << " Getting ready to call ErrorClass: kount = " << kount << std::endl;
02969 throw ErrorClass( "index mapping problem in generating multicommodity cut");
02970 }else{
02971 m_BmatrixIdx[ m_numBmatrixNonz++ ] = m_xVarIndexMap[ indexPair];
02972 cutString << m_variableNames[ m_xVarIndexMap[ indexPair] ];
02973 }
02974 }
02975 }
02976 kount++;
02977 }
02978 }
02979 cutString << std::endl;
02980
02981
02982
02983
02984
02985
02986 m_numBmatrixCon++;
02987 m_pntBmatrix[ m_numBmatrixCon ] = m_numBmatrixNonz;
02988
02989
02990
02991
02992
02993
02994 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
02995 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
02996
02997 m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 1;
02998 scatterValues[ m_BmatrixIdx[ j] ] = m_BmatrixVal[ j];
02999
03000
03001 }
03002
03003
03004
03005
03006
03007 if(numTheta != m_numThetaVar)throw
03008 ErrorClass( "Inconsistent Number of theta variables in multicommondity cut separation problem" );
03009
03010
03011 for(thetaIdx = 0; thetaIdx < m_numThetaVar ; thetaIdx++){
03012
03013
03014 if(m_convexityRowIndex[ thetaIdx] == k){
03015
03016
03017 rowValue = 0;
03018 for(j = m_thetaPnt[ thetaIdx]; j < m_thetaPnt[ thetaIdx + 1] ; j++){
03019
03020
03021
03022 rowValue += m_tmpScatterArray[ m_thetaIndex[ j] ]*scatterValues[ m_thetaIndex[ j] ];
03023
03024 }
03025
03026 if(rowValue > m_osDecompParam.zeroTol || -rowValue > m_osDecompParam.zeroTol){
03027
03028 m_newRowColumnIdx[ numNewRows][ m_newRowNonz[ numNewRows ] ] = thetaIdx ;
03029 m_newRowColumnValue[ numNewRows][ m_newRowNonz[ numNewRows]++ ] = rowValue;
03030 }
03031
03032
03033 }
03034
03035
03036 }
03037
03038 m_newRowLB[ numNewRows] = -OSDBL_MAX;
03039 m_newRowUB[ numNewRows] = 0;
03040 numNewRows++;
03041
03042
03043
03044 for(j = m_pntBmatrix[ m_numBmatrixCon - 1] ;
03045 j < m_pntBmatrix[ m_numBmatrixCon ] ; j++){
03046
03047 m_tmpScatterArray[ m_BmatrixIdx[ j] ] = 0;
03048 scatterValues[ m_BmatrixIdx[ j] ] = 0;
03049
03050 }
03051
03052
03053
03054
03055
03056
03057
03058
03059 m_BmatrixRowIndex[ m_numBmatrixCon - 1] = k;
03061
03062 }
03063
03064 solver->osiSolver->setObjCoeff( i, 0 );
03065
03066
03067
03068 }
03069 std::cout << std::endl << std::endl;
03070
03071
03072
03073 for(i = 0; i < numVar; i++ ) solver->osiSolver->setObjCoeff( i, 0 );
03074
03075
03076 }
03077
03078 delete[] wcoeff;
03079 wcoeff = NULL;
03080
03081 delete[] scatterValues;
03082 scatterValues = NULL;
03083
03084 numNonz = m_newRowNonz;
03085 colIdx = m_newRowColumnIdx;
03086 values = m_newRowColumnValue;
03087 rowUB = m_newRowUB;
03088 rowLB = m_newRowLB;
03089
03090 for(i = 0; i < numNewRows; i++){
03091
03092
03093
03094 m_numThetaVar++;
03095 m_convexityRowIndex[ m_numThetaVar] = -1;
03096 m_thetaPnt[ m_numThetaVar] = m_numThetaNonz;
03097 }
03098
03099 return;
03100
03101
03102 } catch (const ErrorClass& eclass) {
03103 if(wcoeff != NULL){
03104 delete[] wcoeff;
03105 wcoeff = NULL;
03106 }
03107
03108
03109
03110 if(scatterValues != NULL) {
03111 delete[] scatterValues;
03112 scatterValues = NULL;
03113 }
03114
03115
03116 throw ErrorClass(eclass.errormsg);
03117
03118 }
03119
03120
03121 }
03122
03123 void OSBearcatSolverXij::getCutsX(const double* x, const int numX,
03124 int &numNewRows, int* &numNonz, int** &colIdx,
03125 double** &values, double* &rowLB, double* &rowUB) {
03126
03127
03128
03129 int i;
03130 int j;
03131 int k;
03132 int index;
03133 int rowKount;
03134
03135
03136 int indexAdjust = m_numNodes - m_numHubs;
03137 double* tmpRhs;
03138 int numSepRows = m_osinstanceSeparation->getConstraintNumber() ;
03139
03140 tmpRhs = new double[ numSepRows ];
03141
03142 for(i = 0; i < numSepRows; i++){
03143
03144 tmpRhs[ i] = 0;
03145 }
03146
03147 try{
03148 m_osinstanceSeparation->bConstraintsModified = true;
03149
03150 for(i = 0; i < numX; i++){
03151
03152
03153 if(x[ i] > m_osDecompParam.zeroTol){
03154
03155
03156 rowKount = m_separationIndexMap[ i ];
03157
03158 if(rowKount < OSINT_MAX ){
03159
03160 tmpRhs[ rowKount] -= x[ i];
03161
03162 }
03163
03164 }
03165 }
03166
03167 for(i = indexAdjust; i < numSepRows - 1; i++){
03168
03169 if(-tmpRhs[ i] > 1 + m_osDecompParam.zeroTol ){
03170
03171
03172
03173
03174
03175 int tmpKount = indexAdjust;
03176 for(int i1 = m_numHubs; i1 < m_numNodes; i1++){
03177
03178 for(int j1 = i1+1; j1 < m_numNodes; j1++){
03179
03180 if(tmpKount == i){
03181
03182 numNewRows = 1;
03183
03184 m_newRowNonz[ 0] = 2;
03185 m_newRowUB[ 0] = 1;
03186 m_newRowLB[ 0] = 0;
03187
03188 m_newRowColumnIdx[ 0][ 0 ] = i1*(m_numNodes - 1) + j1 - 1;
03189 m_newRowColumnIdx[ 0][ 1 ] = j1*(m_numNodes - 1) + i1;
03190 m_newRowColumnValue[ 0][ 0] = 1;
03191 m_newRowColumnValue[ 0][ 1] = 1;
03192
03193 numNonz = m_newRowNonz;
03194 colIdx = m_newRowColumnIdx;
03195 values = m_newRowColumnValue;
03196 rowUB = m_newRowUB;
03197 rowLB = m_newRowLB;
03198
03199 delete[] tmpRhs;
03200 tmpRhs = NULL;
03201 return;
03202
03203
03204
03205 }
03206
03207 tmpKount++;
03208
03209 }
03210
03211 }
03212
03213
03214 }
03215
03216 m_separationClpModel->setRowUpper(i, tmpRhs[ i] );
03217 m_separationClpModel->setRowLower(i, tmpRhs[ i] );
03218
03219 }
03220
03221
03222
03223
03224
03225 std::vector<int> dualIdx;
03226 std::vector<int>::iterator vit1;
03227 std::vector<int>::iterator vit2;
03228
03229
03230
03231
03232 for(k = 0; k < indexAdjust; k++){
03233 std::cout << std::endl << std::endl;
03234
03235
03236 m_separationClpModel->setRowUpper(k, 0.0);
03237
03238
03239 m_separationClpModel->primal();
03240
03241 if(m_separationClpModel->getObjValue() > m_osDecompParam.zeroTol){
03242 std::cout << "DOING SEPARATION FOR NODE " << k + m_numHubs << std::endl;
03243 std::cout << "SEPERATION OBJ = " << m_separationClpModel->getObjValue() << std::endl;
03244 numNewRows = 1;
03245 m_newRowNonz[ 0] = 0;
03246 m_newRowLB[ 0] = 0;
03247
03248 for(i = 0; i < m_numNodes - m_numHubs ; i++){
03249
03250 if( m_separationClpModel->getRowPrice()[ i] - m_osDecompParam.zeroTol <= -1) dualIdx.push_back( i) ;
03251 }
03252
03253 for (vit1 = dualIdx.begin(); vit1 != dualIdx.end(); vit1++) {
03254
03255 i = *vit1 + m_numHubs;
03256
03257 for (vit2 = dualIdx.begin(); vit2 != dualIdx.end(); vit2++) {
03258
03259 j = *vit2 + m_numHubs;
03260
03261 if( i > j ){
03262
03263 index = i*(m_numNodes -1) + j;
03264 std::cout << "CUT VARIABLE = " << m_variableNames[ index] <<std::endl;
03265 m_newRowColumnValue[ 0][ m_newRowNonz[ 0] ] = 1.0;
03266 m_newRowColumnIdx[ 0][ m_newRowNonz[ 0]++ ] = index;
03267
03268 }else{
03269
03270 if( i < j ){
03271
03272 index = i*(m_numNodes -1) + j - 1;
03273 std::cout << "CUT VARIABLE = " << m_variableNames[ index] <<std::endl;
03274 m_newRowColumnValue[ 0][ m_newRowNonz[ 0] ] = 1.0;
03275 m_newRowColumnIdx[ 0][ m_newRowNonz[ 0]++ ] = index;
03276
03277 }
03278 }
03279
03280 }
03281 }
03282
03283
03284 m_newRowUB[ 0] = dualIdx.size() - 1;
03285
03286 dualIdx.clear();
03287
03288
03289 for(i = indexAdjust; i < numSepRows - 1; i++){
03290
03291 m_separationClpModel->setRowUpper(i, 0.0 );
03292 m_separationClpModel->setRowLower(i, 0.0 );
03293
03294
03295 }
03296 m_separationClpModel->setRowUpper(k, 1.0);
03297 delete[] tmpRhs;
03298 tmpRhs = NULL;
03299
03300
03301 numNonz = m_newRowNonz;
03302 colIdx = m_newRowColumnIdx;
03303 values = m_newRowColumnValue;
03304 rowUB = m_newRowUB;
03305 rowLB = m_newRowLB;
03306
03307 return;
03308
03309
03310
03311 }
03312 m_separationClpModel->setRowUpper(k, 1.0);
03313 dualIdx.clear();
03314
03315 }
03316
03317
03318
03319
03320 for(i = indexAdjust; i < numSepRows - 1; i++){
03321
03322 m_separationClpModel->setRowUpper(i, 0.0 );
03323 m_separationClpModel->setRowLower(i, 0.0 );
03324
03325
03326 }
03327
03328 delete[] tmpRhs;
03329 tmpRhs = NULL;
03330
03331 } catch (const ErrorClass& eclass) {
03332
03333 throw ErrorClass(eclass.errormsg);
03334
03335 }
03336
03337
03338 }
03339
03340
03341 void OSBearcatSolverXij::calcReducedCost( const double* yA, const double* yB){
03342
03343 int k;
03344 int i;
03345 int j;
03346 int l;
03347 int kount;
03348
03349 int tmpVal;
03350 tmpVal = m_numNodes - 1;
03351
03352 for(k = 0; k < m_numHubs; k++){
03353 kount = 0;
03354
03355 for(l = 1; l <= m_upperBoundL[ k]; l++){
03356
03357
03358 for(i = 0; i< m_numNodes; i++){
03359
03360
03361
03362 for(j = 0; j < i; j++){
03363
03364 if(j < m_numHubs){
03365
03366 m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j ] ;
03367
03368 }else{
03369
03370 m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j ] - yA[ j - m_numHubs] ;
03371 }
03372
03373
03374 }
03375
03376
03377
03378 for(j = i + 1; j < m_numNodes; j++){
03379
03380
03381 if(j < m_numHubs){
03382
03383 m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j - 1 ];
03384
03385 } else {
03386
03387
03388 m_rc[k][ kount++] = l*m_cost[ i*tmpVal + j - 1 ] - yA[ j - m_numHubs ];
03389
03390 }
03391
03392 }
03393
03394
03395 }
03396
03397
03398 }
03399
03400
03401 }
03402
03403
03404
03405
03406 int startPnt ;
03407
03408 for(i = 0; i < m_numBmatrixCon; i++){
03409
03410
03411
03412 for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
03413
03414
03415
03416 for(k = 0; k < m_numHubs; k++){
03417
03418
03419 for(l = 1; l <= m_upperBoundL[ k]; l++){
03420
03421
03422 startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
03423
03424 if(m_BmatrixRowIndex[ i] == -1 || m_BmatrixRowIndex[ i] == k ) m_rc[ k][ startPnt + m_BmatrixIdx[ j] ] -= yB[ i]*m_BmatrixVal[ j];
03425
03426 }
03427
03428 }
03429
03430
03431 }
03432
03433 }
03434
03435 }
03436
03437
03438 void OSBearcatSolverXij::createVariableNames( ){
03439
03440 int i;
03441 int j;
03442 int kount;
03443
03444 kount = 0;
03445
03446 for(i = 0; i< m_numNodes; i++){
03447
03448
03449 for(j = 0; j < i; j++){
03450
03451 if(m_nodeName[ i] == "") m_variableNames[ kount] = makeStringFromInt("x(" , i);
03452 else m_variableNames[ kount] = "x(" + m_nodeName[ i];
03453 if(m_nodeName[ i] == "") m_variableNames[ kount] += makeStringFromInt( "," , j);
03454 else m_variableNames[ kount] += "," + m_nodeName[ j];
03455 m_variableNames[ kount] += ")";
03456
03457
03458 kount++;
03459
03460 }
03461
03462 for(j = i + 1; j < m_numNodes; j++){
03463
03464 if(m_nodeName[ i] == "") m_variableNames[ kount] = makeStringFromInt("x(" , i);
03465 else m_variableNames[ kount] = "x(" + m_nodeName[ i];
03466 if(m_nodeName[ i] == "") m_variableNames[ kount] += makeStringFromInt( "," , j);
03467 else m_variableNames[ kount] += "," + m_nodeName[ j];
03468 m_variableNames[ kount] += ")";
03469
03470
03471 kount++;
03472
03473 }
03474
03475
03476 }
03477 }
03478
03479 void OSBearcatSolverXij::createAmatrix(){
03480
03481
03482
03483
03484
03485
03486
03487 int i;
03488 int j;
03489 int numNonz;
03490
03491
03492 m_pntAmatrix[ 0] = 0;
03493 numNonz = 0;
03494
03495 for(j = m_numHubs; j < m_numNodes; j++){
03496
03497
03498 for(i = 0; i < j; i++){
03499
03500 m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j - 1 ;
03501
03502 }
03503
03504 for(i = j + 1; i < m_numNodes; i++){
03505
03506 m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j ;
03507
03508 }
03509
03510 m_pntAmatrix[ j - m_numHubs + 1] = numNonz;
03511
03512 }
03513
03514
03515
03516
03517
03518
03519
03520
03521
03522
03523
03524
03525
03526 }
03527
03528 void OSBearcatSolverXij::pauHana( std::vector<int> &m_zOptIndexes,
03529 std::vector<double> &m_zRootLPx_vals,
03530 int numNodes, int numColsGen, std::string message){
03531
03532 std::cout << std::endl;
03533 std::cout << " PAU HANA TIME! " << std::endl;
03534
03535 double cost;
03536 cost = 0;
03537 std::vector<int>::iterator vit;
03538 try{
03539 int i;
03540 int j;
03541
03542
03543
03544
03545
03546
03547
03548
03549
03550
03551
03552
03553
03554
03555
03556 double routeDemand;
03557
03558
03559 int ivalue;
03560 int jvalue;
03561
03562 for(vit = m_zOptIndexes.begin() ; vit != m_zOptIndexes.end(); vit++){
03563
03564 i = *vit;
03565 std::cout << "x variables for column " << i << " CONVEXITY ROW = "<< m_convexityRowIndex[ i] << std::endl;
03566
03567
03568
03569 routeDemand = 0;
03570
03571 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
03572
03573
03574
03575 std::cout << m_variableNames[ m_thetaIndex[ j] ] << " = " << 1 << " DISTANCE = " << m_cost[ m_thetaIndex[ j] ] << std::endl;
03576
03577 ivalue = m_thetaIndex[ j] /(m_numNodes - 1);
03578
03579 jvalue = m_thetaIndex[ j] - ivalue*(m_numNodes - 1);
03580
03581 if( jvalue >= ivalue ){
03582
03583
03584
03585 routeDemand += m_demand[ jvalue + 1];
03586
03587
03588 }else{
03589
03590
03591
03592 routeDemand += m_demand[ jvalue ];
03593 }
03594 }
03595
03596 std::cout << "route demand = " << routeDemand << std::endl << std::endl;
03597 std::cout << "distance for this route " << m_thetaCost[ i ] / routeDemand << std::endl;
03598
03599 }
03600
03601
03602
03603
03604 std::cout << std::endl << std::endl;
03605 std::cout << message << std::endl;
03606 std::cout << "LINEAR PROGRAMMING RELAXATION VALUE = " << m_rootLPValue << std::endl;
03607 std::cout << "NONLINEAR RELAXATION VALUE = " << calcNonlinearRelax( m_zRootLPx_vals) << std::endl;
03608 std::cout << "LOWER BOUND VALUE = " << m_bestLPValue << std::endl;
03609 std::cout << "FINAL BEST IP SOLUTION VALUE = " << m_bestIPValue << std::endl;
03610 std::cout << "NUMBER OF COLUMNS IN FINAL MASTER = " << m_numThetaVar << std::endl;
03611
03612
03613 std::cout << "NUMBER OF GENERATED COLUMNS = " << numColsGen << std::endl;
03614 std::cout << "NUMBER OF GENERATED CUTS = " << m_numBmatrixCon << std::endl;
03615 std::cout << "NUMBER OF NODES = " << numNodes << std::endl;
03616 std::cout << " PAU!!!" << std::endl;
03617
03618 std::cout << std::endl << std::endl;
03619
03620 std::cout << std::endl << std::endl;
03621 }catch (const ErrorClass& eclass) {
03622
03623 throw ErrorClass(eclass.errormsg);
03624
03625 }
03626
03627 }
03628
03629 double OSBearcatSolverXij::calcNonlinearRelax( std::vector<double> &m_zRootLPx_vals){
03630
03631
03632 std::vector<double>::iterator dvit;
03633 std::vector<double>::iterator dvit2;
03634 int index = 0;
03635 int i;
03636 int j;
03637 int ivalue;
03638 int jvalue;
03639 double *hubDemand = NULL;
03640 double *hubCost = NULL;
03641 hubDemand = new double[m_numHubs ];
03642 hubCost = new double[ m_numHubs];
03643
03644 double objVal = 0.0;
03645 double objValAux = 0.0;
03646
03647 std::map<int, std::vector<double> > extPointDemands;
03648 std::map<int, std::vector<double> > extPointCosts;
03649 std::map<int, std::vector<double> > extPointValues;
03650
03651 std::map<int, std::vector<double> >::iterator mit;
03652
03653 double tmpDemand;
03654 double tmpCost;
03655
03656 for(i = 0; i < m_numHubs; i++){
03657 hubDemand[ i] = 0;
03658 hubCost[ i] = 0;
03659
03660 }
03661
03662 try{
03663 for (dvit = m_zRootLPx_vals.begin(); dvit < m_zRootLPx_vals.end(); dvit++ ){
03664
03665 if(*dvit > m_osDecompParam.zeroTol){
03666
03667
03668
03669
03670 tmpDemand = 0;
03671 tmpCost = 0;
03672
03673 for(j = m_thetaPnt[ index]; j < m_thetaPnt[ index + 1] ; j++){
03674
03675
03676
03677
03678 hubCost[ m_hubPoint[ m_convexityRowIndex[ index] ] ] += m_cost[ m_thetaIndex[ j] ]*( *dvit);
03679 tmpCost += m_cost[ m_thetaIndex[ j] ];
03680
03681 ivalue = m_thetaIndex[ j] /(m_numNodes - 1);
03682 jvalue = m_thetaIndex[ j] - ivalue*(m_numNodes - 1);
03683
03684 if( jvalue >= ivalue ){
03685
03686
03687
03688 hubDemand[ m_hubPoint[ m_convexityRowIndex[ index] ] ] += m_demand[ jvalue + 1]*( *dvit);
03689 tmpDemand += m_demand[ jvalue + 1];
03690
03691
03692
03693 }else{
03694
03695
03696
03697 hubDemand[ m_hubPoint[ m_convexityRowIndex[ index] ] ] += m_demand[ jvalue ]*( *dvit);
03698 tmpDemand += m_demand[ jvalue];
03699
03700 }
03701 }
03702
03703
03704 extPointDemands[ m_convexityRowIndex[ index] ].push_back( tmpDemand);
03705 extPointCosts[ m_convexityRowIndex[ index] ].push_back( tmpCost);
03706 extPointValues[ m_convexityRowIndex[ index] ].push_back(*dvit);
03707
03708 }
03709
03710
03711
03712 index++;
03713
03714 }
03715
03716 int mapSize;
03717 objValAux = 0;
03718 for (i = 0; i < m_numHubs; i++){
03719
03720
03721 objVal += hubDemand[ m_hubPoint[ i] ]*hubCost[ m_hubPoint[ i]];
03722 tmpDemand = 0;
03723 tmpCost = 0;
03724
03725 mapSize = extPointDemands[ m_hubPoint[ i] ].size();
03726
03727
03728
03729 for (j = 0; j < mapSize; j++){
03730
03731 tmpDemand += extPointDemands[ m_hubPoint[ i] ][j]*extPointValues[ m_hubPoint[ i] ][j];
03732 tmpCost += extPointCosts[ m_hubPoint[ i] ][j]*extPointValues[ m_hubPoint[ i] ][j];
03733
03734
03735
03736 objValAux +=
03737 extPointCosts[ m_hubPoint[ i] ][j]*extPointDemands[ m_hubPoint[ i] ][j]*extPointValues[ m_hubPoint[ i] ][j];
03738
03739 }
03740
03741
03742
03743
03744
03745 }
03746
03747 std::cout << "AUXILIARY VALUE COST = " << objValAux << std::endl;
03748
03749
03750 delete[] hubDemand;
03751 hubDemand = NULL;
03752 delete[] hubCost;
03753 hubCost = NULL;
03754
03755
03756 return objVal;
03757
03758 }catch (const ErrorClass& eclass) {
03759
03760
03761 delete[] hubDemand;
03762 hubDemand = NULL;
03763 delete[] hubCost;
03764 hubCost = NULL;
03765
03766
03767 throw ErrorClass(eclass.errormsg);
03768
03769 }
03770
03771
03772 }
03773
03774
03775 OSInstance* OSBearcatSolverXij::getSeparationInstance(){
03776
03777
03778
03779
03780 m_osinstanceSeparation = NULL;
03781
03782
03783
03784
03785 int kountNonz;
03786 int kount;
03787 int startsIdx;
03788
03789 int numYvar = (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1);
03790 int numVvar = m_numNodes - m_numHubs;
03791
03792 int numCon = (m_numNodes - m_numHubs) + (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1)/2 + 1;
03793 double *values = new double[ 2*numYvar + 2*numVvar];
03794 int *indexes = new int[ 2*numYvar + 2*numVvar];
03795 int *starts = new int[ numYvar + numVvar + 1];
03796 starts[ 0] = 0;
03797 startsIdx = 0;
03798 startsIdx++;
03799 kountNonz = 0;
03800 int i;
03801 int j;
03802
03803
03804 std::string separationVarName;
03805 std::string separationConName;
03806
03807 try {
03808
03809 m_osinstanceSeparation = new OSInstance();
03810
03811
03812
03813 m_osinstanceSeparation->setInstanceDescription("The Tour Breaking Separation Problem");
03814
03815
03816
03817 m_osinstanceSeparation->setConstraintNumber( numCon);
03818
03819
03820
03821 for( i = 0; i < m_numNodes - m_numHubs ; i++){
03822
03823 m_osinstanceSeparation->addConstraint(i, makeStringFromInt("nodeRow_", i+ m_numHubs ) , 0.0, 1.0, 0);
03824
03825 }
03826
03827
03828
03829 int rowKounter;
03830 rowKounter = m_numNodes - m_numHubs;
03831
03832 for(i = m_numHubs; i < m_numNodes; i++){
03833
03834
03835
03836 for(j = i+1; j < m_numNodes; j++){
03837 separationConName = makeStringFromInt("Row_(", i);
03838 separationConName += makeStringFromInt(",", j);
03839 separationConName += ")";
03840
03841 m_osinstanceSeparation->addConstraint(rowKounter++, separationConName , 0, 0, 0);
03842 }
03843
03844 }
03845
03846
03847 m_osinstanceSeparation->addConstraint(rowKounter++, "kludgeRow" , 0, m_numNodes, 0);
03848
03849
03850 m_osinstanceSeparation->setVariableNumber( numYvar + numVvar);
03851
03852
03853
03854 std::cout << "NUMBER OF VARIABLES SET = " << numYvar + numVvar << std::endl;
03855
03856 for(i = 0; i < numVvar; i++){
03857
03858 separationVarName = makeStringFromInt("v", i + m_numHubs);
03859
03860 m_osinstanceSeparation->addVariable(i, separationVarName, 0, 1, 'C');
03861
03862 values[ kountNonz ] = -1.0;
03863 indexes[ kountNonz ] = i;
03864 kountNonz++;
03865
03866 values[ kountNonz ] = 1.0;
03867 indexes[ kountNonz ] = rowKounter - 1;
03868 kountNonz++;
03869
03870
03871
03872 starts[ startsIdx++ ] = kountNonz;
03873
03874
03875 }
03876
03877 kount = numVvar;
03878
03879 int i1;
03880 int j1;
03881 int kountCon;
03882 kountCon = m_numNodes - m_numHubs;
03883
03884 for(i1 = 0; i1 < m_numNodes - m_numHubs; i1++){
03885
03886
03887 i = i1 + m_numHubs;
03888
03889 for(j1 = i1 + 1; j1 < m_numNodes - m_numHubs; j1++){
03890
03891
03892 j = j1 + m_numHubs;
03893
03894 separationVarName = makeStringFromInt("y(", i);
03895 separationVarName += makeStringFromInt(",", j);
03896 separationVarName += ")";
03897 m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
03898
03899
03900
03901
03902 m_separationIndexMap[ i*(m_numNodes - 1) + (j - 1) ] = kountCon;
03903
03904 values[ kountNonz ] = 1.0;
03905 indexes[ kountNonz ] = i1;
03906 kountNonz++;
03907
03908 values[ kountNonz ] = -1.0;
03909 indexes[ kountNonz ] = kountCon ;
03910 kountNonz++;
03911
03912 starts[ startsIdx++ ] = kountNonz;
03913
03914
03915
03916
03917 separationVarName = makeStringFromInt("y(", j );
03918 separationVarName += makeStringFromInt(",", i);
03919 separationVarName += ")";
03920 m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
03921
03922 values[ kountNonz ] = 1.0;
03923 indexes[ kountNonz ] = j1;
03924 kountNonz++;
03925
03926
03927 m_separationIndexMap[ j*(m_numNodes - 1) + i ] = kountCon;
03928
03929 values[ kountNonz ] = -1.0;
03930 indexes[ kountNonz ] = kountCon ;
03931 kountNonz++;
03932
03933 starts[ startsIdx++ ] = kountNonz;
03934
03935
03936 kountCon++;
03937
03938
03939 }
03940
03941 }
03942
03943 std::cout << "NUMBER OF VARIABLES ADDED = " << kount << std::endl;
03944
03945
03946 m_osinstanceSeparation->setObjectiveNumber( 1);
03947 SparseVector *objcoeff = new SparseVector( numVvar);
03948
03949
03950 for(i = 0; i < numVvar; i++){
03951
03952 objcoeff->indexes[ i] = i;
03953 objcoeff->values[ i] = 1.0;
03954
03955 }
03956
03957
03958
03959
03960
03961 m_osinstanceSeparation->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
03962
03963
03964 m_osinstanceSeparation->setLinearConstraintCoefficients(kountNonz , true,
03965 values, 0, kountNonz - 1, indexes, 0, kountNonz - 1, starts, 0, startsIdx);
03966
03967
03968
03969
03970
03971 delete objcoeff;
03972
03973
03974
03975
03976
03977
03978 CoinPackedMatrix* matrix;
03979 bool columnMajor = true;
03980 double maxGap = 0;
03981 matrix = new CoinPackedMatrix(
03982 columnMajor,
03983 columnMajor? m_osinstanceSeparation->getConstraintNumber() : m_osinstanceSeparation->getVariableNumber(),
03984 columnMajor? m_osinstanceSeparation->getVariableNumber() : m_osinstanceSeparation->getConstraintNumber(),
03985 m_osinstanceSeparation->getLinearConstraintCoefficientNumber(),
03986 columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->values : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->values,
03987 columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->indexes : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->indexes,
03988 columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->starts : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->starts,
03989 0, 0, maxGap );
03990
03991 ClpNetworkMatrix network( *matrix);
03992
03993 m_separationClpModel = new ClpSimplex();
03994
03995
03996 m_separationClpModel->setOptimizationDirection( 1);
03997 m_separationClpModel->loadProblem( network, m_osinstanceSeparation->getVariableLowerBounds(),
03998 m_osinstanceSeparation->getVariableUpperBounds(),
03999 m_osinstanceSeparation->getDenseObjectiveCoefficients()[0],
04000 m_osinstanceSeparation->getConstraintLowerBounds(), m_osinstanceSeparation->getConstraintUpperBounds()
04001 );
04002
04003 m_separationClpModel->factorization()->maximumPivots(200 + m_separationClpModel->numberRows() / 100);
04004
04005
04006 delete matrix;
04007
04008 }catch (const ErrorClass& eclass) {
04009
04010 throw ErrorClass(eclass.errormsg);
04011
04012 }
04013
04014 return NULL;
04015 }
04016
04017
04018
04019 int OSBearcatSolverXij::getBranchingVar(const double* theta, const int numThetaVar ) {
04020
04021 int varIdx;
04022 varIdx = -1;
04023 int i;
04024 int j;
04025 int numVar = m_numNodes*m_numNodes - m_numHubs ;
04026
04027 double from1Distance;
04028 double from0Distance;
04029 double fraction;
04030 double minFraction;
04031
04032 double *xvalues;
04033
04034
04035 xvalues = new double[ numVar];
04036 for(i = 0; i < numVar; i++){
04037 xvalues[ i] = 0;
04038 }
04039
04040 try{
04041 if(numThetaVar != m_numThetaVar) throw ErrorClass("inconsistent number of variables in getBranchingVar");
04042
04043 for(i = 0; i < m_numThetaVar; i++){
04044
04045 if( ( theta[ i ] > m_osDecompParam.zeroTol ) && ( theta[ i ] < 1 - m_osDecompParam.zeroTol ) ){
04046
04047 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
04048
04049 xvalues[ m_thetaIndex[ j] ] += theta[ i ] ;
04050
04051 }
04052
04053 }
04054
04055
04056 }
04057
04058
04059 minFraction = 1.0;
04060
04061
04062 for(i = 0; i < m_numHubs; i++){
04063
04064 for( j = 0; j < i; j++){
04065
04066
04067 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
04068 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
04069 fraction = std::max(from1Distance, from0Distance);
04070
04071 if(fraction < minFraction){
04072
04073 minFraction = fraction;
04074 varIdx = i*(m_numNodes - 1) + j;
04075 }
04076
04077 }
04078
04079 for(j = i + 1; j < m_numNodes; j++){
04080
04081
04082
04083 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
04084 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
04085 fraction = std::max(from1Distance, from0Distance);
04086
04087 if(fraction < minFraction) {
04088
04089 minFraction = fraction;
04090 varIdx = i*(m_numNodes - 1) + j - 1;
04091 }
04092
04093
04094 }
04095
04096 }
04097
04098
04099
04100 if(minFraction > 1 - m_osDecompParam.zeroTol){
04101
04102 for(i = m_numHubs; i < m_numNodes; i++){
04103
04104 for( j = 0; j < i; j++){
04105
04106
04107 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
04108 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
04109 fraction = std::max(from1Distance, from0Distance);
04110
04111 if(fraction < minFraction) {
04112
04113 minFraction = fraction;
04114 varIdx = i*(m_numNodes - 1) + j ;
04115 }
04116
04117 }
04118
04119 for(j = i + 1; j < m_numNodes; j++){
04120
04121
04122
04123 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
04124 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
04125 fraction = std::max(from1Distance, from0Distance);
04126
04127 if(fraction < minFraction) {
04128
04129 minFraction = fraction;
04130 varIdx = i*(m_numNodes - 1) + j - 1;
04131 }
04132
04133 }
04134
04135 }
04136
04137 }
04138
04139
04140 delete[] xvalues;
04141
04142 xvalues = NULL;
04143
04144 return varIdx;
04145
04146 }catch (const ErrorClass& eclass) {
04147
04148 delete[] xvalues;
04149 xvalues = NULL;
04150
04151 throw ErrorClass(eclass.errormsg);
04152
04153 }
04154
04155
04156 }
04157
04158
04159
04160 int OSBearcatSolverXij::getBranchingVar(const int* thetaIdx, const double* theta,
04161 const int numThetaVar) {
04162
04163 int varIdx;
04164 varIdx = -1;
04165 int i;
04166 int j;
04167 int numVar = m_numNodes*m_numNodes - m_numHubs ;
04168 double from1Distance;
04169 double from0Distance;
04170 double fraction;
04171 double minFraction;
04172
04173 double *xvalues;
04174
04175
04176 xvalues = new double[ numVar];
04177 for(i = 0; i < numVar; i++){
04178 xvalues[ i] = 0;
04179 }
04180
04181 try{
04182
04183
04184 for(i = 0; i < numThetaVar; i++){
04185
04186 if( ( theta[ i ] > m_osDecompParam.zeroTol ) && ( theta[ i ] < 1 - m_osDecompParam.zeroTol ) ){
04187
04188 for(j = m_thetaPnt[ thetaIdx[ i] ]; j < m_thetaPnt[ thetaIdx[ i] + 1] ; j++){
04189
04190 xvalues[ m_thetaIndex[ j] ] += theta[ i ] ;
04191
04192 }
04193
04194 }
04195
04196
04197 }
04198
04199
04200
04201 minFraction = 1.0;
04202
04203
04204 for(i = 0; i < m_numHubs; i++){
04205
04206 for( j = 0; j < i; j++){
04207
04208
04209 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
04210 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
04211 fraction = std::max(from1Distance, from0Distance);
04212
04213 if(fraction < minFraction){
04214
04215 minFraction = fraction;
04216 varIdx = i*(m_numNodes - 1) + j;
04217 }
04218
04219 }
04220
04221 for(j = i + 1; j < m_numNodes; j++){
04222
04223
04224
04225 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
04226 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
04227 fraction = std::max(from1Distance, from0Distance);
04228
04229 if(fraction < minFraction) {
04230
04231 minFraction = fraction;
04232 varIdx = i*(m_numNodes - 1) + j - 1;
04233 }
04234
04235
04236 }
04237
04238 }
04239
04240
04241
04242 std::cout << "MIN FRACTION = " << minFraction << std::endl;
04243
04244 if(minFraction > 1 - m_osDecompParam.zeroTol){
04245
04246 for(i = m_numHubs; i < m_numNodes; i++){
04247
04248
04249
04250 for( j = 0; j < i; j++){
04251
04252
04253 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j ];
04254 from0Distance = xvalues[ i*(m_numNodes - 1) + j ];
04255 fraction = std::max(from1Distance, from0Distance);
04256
04257 if(fraction < minFraction) {
04258
04259 minFraction = fraction;
04260 varIdx = i*(m_numNodes - 1) + j ;
04261 }
04262
04263 }
04264
04265 for(j = i + 1; j < m_numNodes; j++){
04266
04267
04268
04269 from1Distance = 1 - xvalues[ i*(m_numNodes - 1) + j - 1 ];
04270 from0Distance = xvalues[ i*(m_numNodes - 1) + j - 1 ];
04271 fraction = std::max(from1Distance, from0Distance);
04272
04273 if(fraction < minFraction) {
04274
04275 minFraction = fraction;
04276 varIdx = i*(m_numNodes - 1) + j - 1;
04277 }
04278
04279 }
04280
04281 }
04282
04283 }
04284
04285
04286
04287 delete[] xvalues;
04288 xvalues = NULL;
04289
04290 return varIdx;
04291
04292 }catch (const ErrorClass& eclass) {
04293
04294 delete[] xvalues;
04295 xvalues = NULL;
04296
04297 throw ErrorClass(eclass.errormsg);
04298
04299 }
04300
04301
04302 }
04303
04304
04305 void OSBearcatSolverXij::getBranchingCut(const double* thetaVar, const int numThetaVar,
04306 const std::map<int, int> &varConMap, int &varIdx, int &numNonz,
04307 int* &indexes, double* &values) {
04308
04309
04310 int i;
04311 int j;
04312 int kount;
04313 numNonz = 0;
04314
04315
04316
04317 try{
04318
04319 if(numThetaVar != m_numThetaVar) throw ErrorClass("inconsistent number of variables in getBranchingCut");
04320
04321
04322 varIdx = getBranchingVar(thetaVar, numThetaVar );
04323 std::cout << "VARIABLE INDEX = " << varIdx << std::endl;
04324 std::cout << "Branching on Variable: " << m_variableNames[ varIdx] << std::endl;
04325
04326
04327
04328
04329 if( varConMap.find( varIdx) == varConMap.end() ){
04330
04331 for(i = 0; i < m_numThetaVar; i++){
04332
04333 kount = 0;
04334
04335 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
04336
04337 if ( m_thetaIndex[ j] == varIdx) kount++ ;
04338
04339 }
04340
04341
04342
04343 if(kount > 0){
04344
04345 branchCutIndexes[ numNonz] = i;
04346 branchCutValues[ numNonz++] = kount ;
04347
04348 }
04349
04350 }
04351
04352
04353
04354 m_BmatrixVal[ m_numBmatrixNonz] = 1;
04355 m_BmatrixIdx[ m_numBmatrixNonz++] = varIdx;
04356 m_numBmatrixCon++;
04357 m_pntBmatrix[ m_numBmatrixCon] = m_numBmatrixNonz;
04358
04359
04360
04361
04362 m_numThetaVar++;
04363 m_convexityRowIndex[ m_numThetaVar] = -1;
04364 m_thetaPnt[ m_numThetaVar] = m_numThetaNonz;
04365
04366
04367
04368
04369 }
04370
04371
04372
04373 indexes = branchCutIndexes;
04374 values = branchCutValues;
04375
04376 return;
04377
04378 }catch (const ErrorClass& eclass) {
04379
04380 throw ErrorClass(eclass.errormsg);
04381
04382 }
04383
04384 }
04385
04386
04387 void OSBearcatSolverXij::getBranchingCut(const int* thetaIdx, const double* thetaVar,
04388 const int numThetaVar, const std::map<int, int> &varConMap,
04389 int &varIdx, int &numNonz, int* &indexes, double* &values) {
04390
04391
04392 int i;
04393 int j;
04394 int kount;
04395 numNonz = 0;
04396
04397
04398
04399 try{
04400
04401
04402
04403 varIdx = getBranchingVar(thetaIdx, thetaVar, numThetaVar );
04404
04405 std::cout << "Branching on Variable: " << m_variableNames[ varIdx] << std::endl;
04406
04407
04408
04409
04410 if( varConMap.find( varIdx) == varConMap.end() ){
04411
04412
04413
04414
04415
04416
04417 for(i = 0; i < m_numThetaVar; i++){
04418
04419 kount = 0;
04420
04421 for(j = m_thetaPnt[ i]; j < m_thetaPnt[ i + 1] ; j++){
04422
04423 if ( m_thetaIndex[ j] == varIdx) kount++ ;
04424
04425 }
04426
04427
04428
04429 if(kount > 0){
04430
04431 branchCutIndexes[ numNonz] = i;
04432 branchCutValues[ numNonz++] = kount ;
04433
04434 }
04435
04436 }
04437
04438
04439
04440 m_BmatrixVal[ m_numBmatrixNonz] = 1.0;
04441 m_BmatrixIdx[ m_numBmatrixNonz++] = varIdx;
04442 m_numBmatrixCon++;
04443 m_pntBmatrix[ m_numBmatrixCon] = m_numBmatrixNonz;
04444
04445
04446
04447
04448 m_numThetaVar++;
04449 m_convexityRowIndex[ m_numThetaVar] = -1;
04450 m_thetaPnt[ m_numThetaVar] = m_numThetaNonz;
04451
04452
04453
04454
04455 }
04456
04457
04458
04459 indexes = branchCutIndexes;
04460 values = branchCutValues;
04461
04462 return;
04463
04464 }catch (const ErrorClass& eclass) {
04465
04466 throw ErrorClass(eclass.errormsg);
04467
04468 }
04469
04470 }
04471
04472
04473 void OSBearcatSolverXij::getInitialSolution(){
04474
04475 int k;
04476 double *routeDemand;
04477 routeDemand = NULL;
04478 routeDemand = new double[ m_numHubs];
04479 std::map<int, std::vector<int> > routeMap;
04480 std::vector<int>::iterator vit;
04481 bool isFeasible;
04482
04483
04484 isFeasible = true;
04485 try{
04486
04487 routeMap = m_initSolMap[ 0];
04488
04489 for(k = 0; k < m_numHubs; k++){
04490
04491 routeDemand[ k] = 0;
04492 for(vit = routeMap[k].begin(); vit!=routeMap[k].end(); ++vit){
04493
04494
04495 routeDemand[ k] += m_demand[ *vit];
04496 }
04497
04498 if(routeDemand[k] > m_routeCapacity[ k] ){
04499
04500 isFeasible = false;
04501 break;
04502
04503 }
04504 }
04505
04506 delete[] routeDemand;
04507 routeDemand = NULL;
04508
04509 if(isFeasible == false) getFeasibleSolution();
04510
04511 OneOPT();
04512
04513
04514 }catch (const ErrorClass& eclass) {
04515
04516 if(routeDemand == NULL){
04517
04518 delete[] routeDemand;
04519 routeDemand = NULL;
04520 }
04521
04522 throw ErrorClass(eclass.errormsg);
04523
04524 }
04525
04526
04527 }
04528
04529
04530 void OSBearcatSolverXij::resetMaster( std::map<int, int> &inVars, OsiSolverInterface *si){
04531
04532 int i;
04533 int j;
04534
04535 int kount;
04536 int numNonz;
04537
04538 std::map<int, int>::iterator mit;
04539
04540 int numVars = inVars.size();
04541 int numVarArt;
04542
04543
04544
04545 numVarArt = m_numNodes + m_numBmatrixCon;
04546
04547
04548 std::vector<double> valuesVec;
04549 double *values = NULL;
04550
04551 std::vector<int> indexesVec;
04552 int *indexes = NULL ;
04553
04554 int *starts = new int[ numVars + 1 + numVarArt];
04555
04556 int startsIdx;
04557
04558
04559
04560
04561 int* thetaPntTmp;
04562 int* thetaIndexTmp;
04563 int* tmpConvexity = new int[ m_numThetaVar];
04564
04565
04566 numNonz = 0;
04567
04568 for(mit = inVars.begin(); mit != inVars.end(); mit++){
04569
04570 numNonz += m_thetaPnt[mit->first + 1 ] - m_thetaPnt[ mit->first ];
04571 }
04572
04573
04574 thetaPntTmp = new int[ numVars + 1];
04575 thetaIndexTmp = new int[ numNonz];
04576
04577
04578
04579 for(mit = inVars.begin(); mit != inVars.end(); mit++){
04580
04581
04582
04583 if( m_convexityRowIndex[ mit->first] == -1) throw ErrorClass( "we have an artificial variable in reset master");
04584
04585
04586 }
04587
04588
04589 kount = 0;
04590 numNonz = 0;
04591 thetaPntTmp[ kount] = 0;
04592
04593 for(mit = inVars.begin(); mit != inVars.end(); mit++){
04594
04595
04596
04597 kount++;
04598
04599 for(i = m_thetaPnt[ mit->first ]; i < m_thetaPnt[mit->first + 1 ]; i++){
04600
04601 thetaIndexTmp[ numNonz++] = m_thetaIndex[ i];
04602
04603
04604
04605 }
04606
04607 thetaPntTmp[ kount] = numNonz;
04608
04609
04610
04611
04612
04613 inVars[ mit->first] = numVarArt + kount - 1 ;
04614
04615 }
04616
04617
04618
04619
04620
04621
04622
04623 for(i = 0; i < m_numThetaVar; i++) tmpConvexity[ i] = m_convexityRowIndex[ i];
04624
04625
04626
04627 m_numThetaVar = 0;
04628 m_numThetaNonz = 0;
04629 for(i = 0; i < numVarArt; i++){
04630
04631 m_convexityRowIndex[ m_numThetaVar] = -1;
04632 m_thetaPnt[ m_numThetaVar++] = 0;
04633
04634
04635 }
04636
04637
04638 intVarSet.clear();
04639 for(mit = inVars.begin(); mit != inVars.end(); mit++){
04640
04641
04642 intVarSet.insert ( std::pair<int,double>(mit->second, 1.0) );
04643
04644
04645
04646
04647 m_convexityRowIndex[ m_numThetaVar] = tmpConvexity[ mit->first];
04648
04649 m_thetaPnt[ m_numThetaVar++ ] = m_numThetaNonz;
04650
04651 for(j = thetaPntTmp[ mit->second - numVarArt]; j < thetaPntTmp[ mit->second - numVarArt + 1 ]; j++){
04652
04653
04654 m_thetaIndex[ m_numThetaNonz ] = thetaIndexTmp[ m_numThetaNonz] ;
04655 m_numThetaNonz++;
04656
04657 }
04658
04659 }
04660
04661 m_thetaPnt[ m_numThetaVar ] = m_numThetaNonz;
04662
04663
04664
04665
04666
04667
04668
04669
04670
04671 si->writeLp( "gailTest" );
04672
04673
04674
04675
04676
04677 numNonz = 0;
04678 startsIdx = 0;
04679 starts[ startsIdx++] = numNonz;
04680
04681 for(i = 0; i < numVarArt; i++){
04682 numNonz++;
04683 starts[ startsIdx++] = numNonz;
04684 valuesVec.push_back( 1.0);
04685 indexesVec.push_back( i);
04686
04687 }
04688
04689
04690 int rowCount;
04691
04692 int numAmatrixRows;
04693 numAmatrixRows = m_numNodes - m_numHubs;
04694
04695 for(mit = inVars.begin(); mit != inVars.end(); mit++){
04696
04697
04698 valuesVec.push_back( 1.0);
04699 indexesVec.push_back( numAmatrixRows + m_convexityRowIndex[ mit->second] );
04700
04701 numNonz++;
04702
04703 for(j = m_thetaPnt[ mit->second ]; j < m_thetaPnt[ mit->second + 1 ]; j++){
04704
04705 m_tmpScatterArray[ m_thetaIndex[ j] ]++;
04706
04707
04708
04709 }
04710
04711
04712
04713
04714 for(i = 0; i < numAmatrixRows; i++){
04715
04716 rowCount = 0;
04717
04718 for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
04719
04720
04721
04722
04723 rowCount += m_tmpScatterArray[ m_Amatrix[ j] ];
04724
04725
04726 }
04727
04728 if(rowCount > 0){
04729
04730 numNonz++;
04731
04732
04733 valuesVec.push_back( rowCount);
04734 indexesVec.push_back( i);
04735
04736
04737 }
04738
04739
04740 }
04741
04742
04743
04744 for(i = 0; i < m_numBmatrixCon; i++){
04745
04746 rowCount = 0;
04747
04748 for(j = m_pntBmatrix[ i]; j < m_pntBmatrix[ i + 1]; j++){
04749
04750
04751
04752
04753 rowCount += m_tmpScatterArray[ m_BmatrixIdx[ j] ];
04754
04755
04756 }
04757
04758 if(rowCount > 0){
04759 numNonz++;
04760
04761
04762
04763 valuesVec.push_back( rowCount);
04764 indexesVec.push_back( i + m_numNodes);
04765 }
04766
04767
04768 }
04769
04770
04771
04772
04773 for(j = m_thetaPnt[ mit->second ]; j < m_thetaPnt[ mit->second + 1 ]; j++){
04774
04775 m_tmpScatterArray[ m_thetaIndex[ j] ] = 0;
04776
04777 }
04778
04779 starts[ startsIdx++] = numNonz;
04780
04781 }
04782
04783
04784
04785 values = new double[ numNonz];
04786 indexes = new int[ numNonz];
04787
04788 if( (unsigned int)numNonz != valuesVec.size() ) throw ErrorClass("dimension problem in reset");
04789 if( (unsigned int)numNonz != indexesVec.size() ) throw ErrorClass("dimension problem in reset");
04790
04791 for(i = 0; i < numNonz; i++){
04792
04793 values[ i] = valuesVec[i];
04794 indexes[ i] = indexesVec[i];
04795
04796 }
04797
04798
04799
04800
04801
04802
04803
04804 delete m_osinstanceMaster;
04805 m_osinstanceMaster = NULL;
04806
04807
04808 m_osinstanceMaster = new OSInstance();
04809 m_osinstanceMaster->setInstanceDescription("The Restricted Master");
04810
04811
04812
04813 m_osinstanceMaster->setVariableNumber( numVars + numVarArt );
04814
04815
04816
04817
04818 SparseVector *objcoeff = new SparseVector( numVars + numVarArt);
04819
04820
04821
04822
04823 m_osinstanceMaster->setConstraintNumber( m_numNodes + m_numBmatrixCon);
04824
04825
04826
04827
04828 int varNumber;
04829 varNumber = 0;
04830
04831
04832
04833 for(i = 0; i < numVarArt; i++){
04834
04835 objcoeff->indexes[ varNumber ] = varNumber ;
04836
04837 objcoeff->values[ varNumber ] = m_osDecompParam.artVarCoeff;
04838
04839 m_thetaCost[ varNumber] = m_osDecompParam.artVarCoeff;
04840
04841 m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt("x", i ) ,
04842 0, 1.0, 'C');
04843
04844
04845 }
04846
04847
04848 kount = 0;
04849 for(mit = inVars.begin(); mit != inVars.end(); mit++){
04850
04851 objcoeff->indexes[ varNumber ] = varNumber ;
04852
04853 objcoeff->values[ varNumber ] = si->getObjCoefficients()[ mit->first] ;
04854
04855 m_thetaCost[ varNumber] = si->getObjCoefficients()[ mit->first];
04856
04857 m_osinstanceMaster->addVariable(varNumber++, makeStringFromInt("x", kount + numVarArt ) ,
04858 0, 1.0, 'C');
04859
04860 kount++;
04861
04862
04863
04864 }
04865
04866
04867
04868 for(i = 0; i < m_numNodes; i++){
04869
04870 m_osinstanceMaster->addConstraint(i, makeStringFromInt("con", i),
04871 1.0, 1.0, 0);
04872
04873 }
04874
04875
04876 for(i = m_numNodes; i < m_numBmatrixCon + m_numNodes; i++){
04877
04878 m_osinstanceMaster->addConstraint(i, makeStringFromInt("con", i),
04879 si->getRowLower()[ i], si->getRowUpper()[ i], 0);
04880
04881
04882 }
04883
04884
04885
04886 m_osinstanceMaster->setObjectiveNumber( 1);
04887 m_osinstanceMaster->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
04888
04889
04890 m_osinstanceMaster->setLinearConstraintCoefficients(numNonz , true,
04891 values, 0, numNonz - 1, indexes, 0, numNonz - 1, starts, 0, startsIdx);
04892
04893
04894
04895
04896
04897 delete[] tmpConvexity;
04898 tmpConvexity = NULL;
04899 delete[] thetaPntTmp;
04900 thetaPntTmp = NULL;
04901 delete[] thetaIndexTmp;
04902 thetaIndexTmp = NULL;
04903 delete objcoeff;
04904 objcoeff = NULL;
04905 }
04906
04907
04908
04909 double OSBearcatSolverXij::getRouteDistance(int numNodes, int hubIndex,
04910 CoinSolver* solver, std::vector<int> zk, double* xVar){
04911
04912 int i;
04913 int j;
04914 int numVar;
04915 numVar = numNodes*numNodes - numNodes;
04916
04917 int kount;
04918 double objValue;
04919 std::vector<int>::iterator vit;
04920 std::map<std::pair<int, int>, int > cutMap;
04921 std::map<std::pair<int, int>, int >::iterator mit;
04922 std::vector<IndexValuePair*> primalValPair;
04923
04924
04925
04926
04927
04928 for(i = 0; i < numVar; i++) xVar[ i] = 0;
04929
04930 int numNewRows;
04931 int* numRowNonz = NULL;
04932 int** colIdx = NULL;
04933 double** rowValues = NULL ;
04934 double* rowLB;
04935 double* rowUB;
04936
04937
04938
04939 OsiSolverInterface *si = solver->osiSolver;
04940 try{
04941
04942
04943 si->setRowUpper(hubIndex, 1.0);
04944 si->setRowUpper(hubIndex + numNodes, 1.0);
04945 si->setRowLower(hubIndex, 1.0);
04946 si->setRowLower(hubIndex + numNodes, 1.0);
04947
04948
04949
04950 kount = zk.size();
04951
04952 int idx1;
04953 int idx2;
04954
04955 for(i = 0; i < kount ; i++){
04956
04957
04958 si->setRowUpper(zk[ i], 1.0);
04959 si->setRowUpper(zk[ i] + numNodes, 1.0);
04960 si->setRowLower(zk[ i], 1.0);
04961 si->setRowLower(zk[ i] + numNodes, 1.0);
04962
04963 idx1 = hubIndex;
04964 idx2 = zk[i ];
04965
04966
04967 if( idx1 > idx2 ){
04968 si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
04969 }else{
04970
04971 if( idx1 < idx2 ){
04972
04973 si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
04974
04975 }
04976 }
04977
04978 idx1 = zk[ i];
04979 idx2 = hubIndex;
04980
04981
04982 if( idx1 > idx2 ){
04983
04984 si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
04985
04986 }else{
04987
04988 if( idx1 < idx2 ){
04989
04990 si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
04991
04992 }
04993 }
04994
04995 for(j = i + 1; j < kount; j++){
04996
04997 idx1 = zk[i];
04998 idx2 = zk[j];
04999
05000
05001 if( idx1 > idx2 ){
05002
05003 si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
05004
05005 }else{
05006
05007 if( idx1 < idx2 ){
05008
05009 si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
05010 }
05011 }
05012
05013
05014 idx1 = zk[j];
05015 idx2 = zk[i];
05016
05017
05018 if( idx1 > idx2 ){
05019 si->setColUpper(idx1*(numNodes -1) + idx2, 1.0);
05020 }else{
05021
05022 if( idx1 < idx2 ){
05023
05024 si->setColUpper(idx1*(numNodes -1) + idx2 - 1, 1.0);
05025 }
05026 }
05027 }
05028 }
05029
05030 solver->solve();
05031
05032
05033
05034 if(solver->osresult->getSolutionStatusType( 0 ) != "optimal" ) throw ErrorClass("Trouble with integer program used to construct initial master");
05035 objValue = solver->osresult->getObjValue(0, 0) ;
05036
05037 primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
05038
05039
05040 bool isCutAdded;
05041 isCutAdded = true;
05042 while(isCutAdded == true){
05043
05044 isCutAdded = false;
05045
05046 for(i = 0; i < numVar; i++) xVar[ primalValPair[ i]->idx ] = primalValPair[ i]->value;
05047
05048 numNewRows = 0;
05049 getCutsX(xVar, numVar, numNewRows, numRowNonz,
05050 colIdx,rowValues, rowLB, rowUB);
05051
05052 if(numNewRows >= 1){
05053 isCutAdded = true;
05054 std::cout << "WE HAVE A CUT " << std::endl;
05055 std::cout << "EXPRESS CUT IN X(I, J) SPACE" << std::endl;
05056
05057
05058 for(i = 0; i < numNewRows; i++){
05059
05060 std::cout << "CUT UPPER BOUND = " << rowUB[ i] << std::endl;
05061 si->addRow(numRowNonz[ i], colIdx[ i], rowValues[ i], rowLB[ i], rowUB[ i] ) ;
05062 }
05063
05064 std::cout << "A CUT WAS ADDED, CALL SOLVE AGAIN" << std::endl;
05065 solver->solve();
05066 if(solver->osresult->getSolutionStatusType( 0 ) != "optimal" ) throw ErrorClass("Trouble with integer program used to construct initial master");
05067 primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
05068 std::cout << "New Solution Status = " << solver->osresult->getSolutionStatusType( 0 ) << std::endl;
05069 std::cout << "Optimal Objective Value = " << solver->osresult->getObjValue(0, 0) << std::endl;
05070
05071
05072
05073
05074 }
05075
05076 }
05077
05078 objValue = solver->osresult->getObjValue(0, 0) ;
05079
05080 primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
05081
05082
05083
05084 std::cout << "Objective Function Value = " << objValue << std::endl;
05085
05086
05087 for(i = 0; i < numVar; i++) si->setColUpper(i, 0);
05088 for(i = 0; i < 2*numNodes; i++) si->setRowUpper(i, 0);
05089 for(i = 0; i < 2*numNodes; i++) si->setRowLower(i, 0);
05090
05091 return objValue;
05092
05093
05094 } catch (const ErrorClass& eclass) {
05095
05096 std::cout << std::endl << std::endl << std::endl;
05097
05098 if( xVar != NULL)
05099 delete xVar;
05100
05101 throw ErrorClass(eclass.errormsg);
05102 }
05103
05104 }
05105
05106
05107
05108 CoinSolver* OSBearcatSolverXij::getTSP(int numNodes, double* cost){
05109
05110
05111 int i;
05112 int j;
05113 int numVar;
05114 numVar = numNodes*numNodes - numNodes;
05115 int numNonz;
05116 int kount;
05117
05118 std::vector<int>::iterator vit;
05119 double* rhsVec;
05120 CoinSolver *solver = NULL;
05121 std::map<std::pair<int, int>, int > cutMap;
05122 std::map<std::pair<int, int>, int >::iterator mit;
05123
05124
05125 int numCuts;
05126
05127 rhsVec = new double[ 2*numNodes];
05128
05129 for(i = 0; i < 2*numNodes; i++){
05130
05131
05132
05133 rhsVec[ i] = 0;
05134 }
05135
05136
05137
05138
05139 numCuts = 0;
05140 std::pair <int,int> pairIJ;
05141 std::pair <int,int> pairJI;
05142
05143 for(i = 0; i < numNodes - 1; i++){
05144
05145
05146 for(j = i + 1; j < numNodes; j++){
05147
05148
05149 pairIJ = std::make_pair(i, j);
05150 pairJI = std::make_pair(j, i);
05151
05152 cutMap[pairIJ ] = 2*numNodes + numCuts;
05153 cutMap[pairJI ] = 2*numNodes + numCuts;
05154 numCuts++;
05155
05156 }
05157 }
05158
05159
05160
05161 OSInstance *osinstance = new OSInstance();
05162 try{
05163
05164 osinstance->setInstanceSource("generated from OSBearcatSoleverXij");
05165 osinstance->setInstanceDescription("Find the minimum tour for a route");
05166 osinstance->setVariableNumber( numVar);
05167
05168 for(i = 0; i < numVar; i++){
05169
05170 osinstance->addVariable(i, m_variableNames[ i], 0, 0, 'B');
05171
05172 }
05173
05174 osinstance->setObjectiveNumber( 1);
05175
05176
05177 SparseVector *objcoeff = new SparseVector( numVar);
05178
05179 for(i = 0; i < numVar; i++){
05180
05181 objcoeff->indexes[ i] = i;
05182 objcoeff->values[ i] = cost[ i];
05183
05184 }
05185
05186 osinstance->addObjective(-1, "minimizeDistance", "min", 0.0, 1.0, objcoeff);
05187 objcoeff->bDeleteArrays = true;
05188 delete objcoeff;
05189
05190 osinstance->setConstraintNumber( 2*numNodes + numCuts);
05191
05192
05193 for(i = 0; i < numNodes; i++){
05194
05195 osinstance->addConstraint(i, makeStringFromInt("enter_node_", i) , rhsVec[i], rhsVec[i], 0);
05196
05197 }
05198
05199 for(i = numNodes; i < 2*numNodes; i++){
05200
05201 osinstance->addConstraint(i, makeStringFromInt("leave_node_", i - numNodes) , rhsVec[i], rhsVec[i], 0);
05202
05203 }
05204
05205
05206
05207 for(i = 0; i < numCuts; i++){
05208
05209 osinstance->addConstraint(2*numNodes + i, makeStringFromInt("XijCut_", i) , 0, 1, 0);
05210
05211 }
05212
05213
05214
05215 numNonz = numVar*numVar - numVar + 2*numCuts;
05216
05217
05218 double *values = new double[ numNonz];
05219 int *rowIndexes = new int[ numNonz];
05220 int *starts = new int[ numVar + 1];
05221
05222
05223 kount = 0;
05224 numNonz = 0;
05225 starts[ kount++] = 0;
05226 for(i = 0; i < numNodes; i++){
05227
05228 for(j = 0; j < i; j++){
05229
05230
05231 rowIndexes[ numNonz] = j;
05232 values[ numNonz++] = 1.0;
05233
05234 rowIndexes[ numNonz] = numNodes + i;
05235 values[ numNonz++] = 1.0;
05236
05237
05238 pairIJ = std::make_pair(i, j);
05239 mit = cutMap.find( pairIJ);
05240 if(mit != cutMap.end() ){
05241
05242 rowIndexes[ numNonz] = mit->second;
05243 values[ numNonz++] = 1.0;
05244 }
05245
05246
05247 starts[ kount++] = numNonz;
05248
05249
05250 }
05251
05252 for(j = i+1; j < numNodes; j++){
05253
05254
05255 rowIndexes[ numNonz] = j;
05256 values[ numNonz++] = 1.0;
05257
05258 rowIndexes[ numNonz] = numNodes + i;
05259 values[ numNonz++] = 1.0;
05260
05261
05262 pairIJ = std::make_pair(i, j);
05263 mit = cutMap.find( pairIJ);
05264 if(mit != cutMap.end() ){
05265
05266 rowIndexes[ numNonz] = mit->second;
05267 values[ numNonz++] = 1.0;
05268 }
05269
05270
05271 starts[ kount++] = numNonz;
05272
05273 }
05274
05275
05276 }
05277
05278 osinstance->setLinearConstraintCoefficients(numNonz, true, values, 0, numNonz - 1,
05279 rowIndexes, 0, numNonz - 1, starts, 0, numVar);
05280
05281
05282
05283
05284 solver = new CoinSolver();
05285 solver->sSolverName ="cbc";
05286 solver->osinstance = osinstance;
05287 solver->osoption = m_osoption;
05288 solver->buildSolverInstance();
05289 solver->osoption = m_osoption;
05290
05291 delete[] rhsVec;
05292 rhsVec = NULL;
05293
05294 return solver;
05295
05296
05297
05298 } catch (const ErrorClass& eclass) {
05299
05300 std::cout << std::endl << std::endl << std::endl;
05301
05302 if(rhsVec != NULL){
05303 delete[] rhsVec;
05304 rhsVec = NULL;
05305 }
05306
05307
05308 throw ErrorClass(eclass.errormsg);
05309 }
05310
05311
05312 }
05313
05314
05315 bool OSBearcatSolverXij::OneOPT(){
05316
05317 int k;
05318 int kprime;
05319 double *routeCost = NULL;
05320 int *routeDemand = NULL;
05321 double *xVar = NULL;
05322 int numXVar = m_numNodes*m_numNodes - m_numNodes;
05323 double routeCostInf;
05324
05325 routeCostInf = OSINT_MAX;
05326 double feasMult = 10000;
05327
05328 routeCost = new double[m_numHubs];
05329 routeDemand = new int[m_numHubs];
05330 xVar = new double[ numXVar];
05331
05332
05333
05334 std::map<int, std::vector<int> > routeMap;
05335 std::vector<int> tmpVec;
05336 std::vector<int>::iterator vit;
05337 std::vector<int>::iterator vit2;
05338
05339 routeMap = m_initSolMap[ 0];
05340 CoinSolver *solver = NULL;
05341
05342 double totalCost;
05343 bool foundMoBetta;
05344 bool foundLocalImprovement;
05345
05346 try{
05347
05348 solver = getTSP(m_numNodes, m_cost);
05349
05350 totalCost = 0;
05351 for(k = 0; k < m_numHubs; k++){
05352
05353 routeDemand[ k] = 0;
05354 for(vit = routeMap[k].begin(); vit!=routeMap[k].end(); ++vit){
05355
05356 std::cout << "node i = " << *vit << " demand " << m_demand[ *vit] << std::endl;
05357 routeDemand[ k] += m_demand[ *vit];
05358 }
05359
05360 if(routeDemand[k] <= m_routeCapacity[ k] ){
05361
05362 routeCost[ k] = getRouteDistance(m_numNodes, k, solver,
05363 routeMap[k], xVar)*routeDemand[ k];
05364
05365
05366 std::cout << "rout = " << k << " cost " << routeCost[ k] << std::endl;
05367 totalCost += routeCost[ k];
05368 }else{
05369 std::cout << "rout demand " << routeDemand[k] << " route capacity = " << m_routeCapacity[ k] << std::endl;
05370 routeCost[ k] = routeCostInf;
05371 totalCost += routeCost[ k];
05372
05373 }
05374
05375
05376 }
05377
05378
05379 foundMoBetta = true;
05380 double improvement;
05381
05382 double tmpCostk;
05383 double tmpCostkPrime;
05384 double optCostk;
05385 double optCostkPrime;
05386 double tmpVal;
05387 int tmpIdx;
05388 std::vector<int>::iterator vitIdx;
05389
05390 while(foundMoBetta == true){
05391 foundMoBetta = false;
05392
05393 for(k = 0; k < m_numHubs; k++){
05394
05395 foundLocalImprovement = false;
05396
05397 for(kprime = 0; kprime < m_numHubs; kprime++){
05398
05399 if(kprime != k && routeCost[ kprime] < routeCostInf){
05400
05401
05402 improvement = 0;
05403 optCostk = routeCostInf;
05404 optCostkPrime = routeCostInf;
05405
05406 for(vit = routeMap[k].begin(); vit!=routeMap[k].end(); ++vit){
05407
05408 foundLocalImprovement = false;
05409
05410
05411 if( m_demand[ *vit] + routeDemand[ kprime] <= m_routeCapacity[ kprime] ){
05412
05413 tmpCostk = routeCostInf;
05414 tmpCostkPrime = routeCostInf;
05415
05416
05417
05418
05419 routeMap[ kprime].push_back( *vit);
05420
05421 tmpCostkPrime = getRouteDistance(m_numNodes, kprime, solver,
05422 routeMap[kprime], xVar)*(routeDemand[ kprime] + m_demand[ *vit] );
05423
05424
05425 routeMap[ kprime].pop_back( );
05426
05427
05428
05429 if(routeCost[k] == routeCostInf ){
05430
05431
05432
05433
05434
05435 if( routeDemand[ k] - m_demand[ *vit] <= m_routeCapacity[ k]) {
05436
05437 for(vit2 = routeMap[k].begin(); vit2 != routeMap[k].end(); vit2++){
05438
05439 if(vit != vit2) tmpVec.push_back( *vit2);
05440
05441 }
05442
05443 tmpCostk = getRouteDistance(m_numNodes, k, solver,
05444 tmpVec, xVar)*(routeDemand[ k] - m_demand[ *vit] );
05445
05446 tmpVec.clear();
05447
05448 }
05449
05450 tmpVal = std::max(*vit, routeDemand[ k] - m_demand[ *vit])*feasMult
05451 + ( routeCost[kprime] - tmpCostkPrime);
05452
05453 if( tmpVal > improvement) {
05454 foundLocalImprovement = true;
05455 improvement = tmpVal;
05456 tmpIdx = *vit;
05457 vitIdx = vit;
05458 optCostk = tmpCostk;
05459 optCostkPrime = tmpCostkPrime;
05460
05461 }
05462
05463
05464 }else{
05465
05466
05467
05468
05469 for(vit2 = routeMap[k].begin(); vit2 != routeMap[k].end(); vit2++){
05470
05471 if(vit != vit2) tmpVec.push_back( *vit2);
05472
05473 }
05474
05475 tmpCostk = getRouteDistance(m_numNodes, k, solver,
05476 tmpVec, xVar)*(routeDemand[ k] - m_demand[ *vit] );
05477
05478 tmpVec.clear();
05479
05480 if( ( (routeCost[k] + routeCost[kprime]) -
05481 ( tmpCostkPrime + tmpCostk ) ) > improvement ) {
05482
05483 improvement = (routeCost[k] + routeCost[kprime]) -
05484 ( tmpCostkPrime + tmpCostk );
05485
05486 foundLocalImprovement = true;
05487 tmpIdx = *vit;
05488 vitIdx = vit;
05489 optCostk = tmpCostk;
05490 optCostkPrime = tmpCostkPrime;
05491
05492
05493
05494 }
05495 }
05496 }
05497
05498 if(foundLocalImprovement == true)break;
05499 }
05500
05501
05502
05503
05504 if( foundLocalImprovement == true) {
05505
05506
05507 std::cout << "index = " << tmpIdx << std::endl;
05508
05509 routeMap[ kprime].push_back( tmpIdx);
05510
05511 routeDemand[ kprime] += m_demand[ tmpIdx];
05512
05513 routeCost[ kprime] = optCostkPrime;
05514
05515
05516
05517
05518
05519
05520 std::cout << "tmpIdx = " << tmpIdx << std::endl;
05521 std::cout << "vitIdx = " << *vitIdx << std::endl;
05522 routeMap[ k].erase( vitIdx);
05523
05524 routeDemand[ k] -= m_demand[ tmpIdx];
05525
05526 routeCost[ k] = optCostk;
05527
05528
05529
05530
05531
05532
05533 foundMoBetta = true;
05534
05535 }
05536 }
05537 }
05538 }
05539
05540
05541 }
05542
05543
05544
05545 totalCost = 0;
05546 for(k = 0; k < m_numHubs; k++){
05547
05548 std::cout << std::endl << std::endl;
05549
05550 std::cout << "ROUTE " << k << " Total Demand = " << routeDemand[ k] << std::endl;
05551 std::cout << "ROUTE " << k << " Total Cost = " << routeCost[ k] << std::endl;
05552 std::cout << "ROUTE " << k << " Nodes " << std::endl;
05553
05554
05555
05556
05557 totalCost += routeCost[ k];
05558 }
05559
05560 std::cout << "Total Cost = " << totalCost << std::endl;
05561
05562 m_initSolMap[ 0] = routeMap;
05563
05564
05565
05566 delete[] routeCost;
05567 routeCost = NULL;
05568 delete[] routeDemand;
05569 routeDemand = NULL;
05570 delete[] xVar;
05571 xVar = NULL;
05572 delete solver->osinstance;
05573 delete solver;
05574
05575 if( totalCost >= routeCostInf )return false;
05576 else return true;
05577
05578
05579 } catch (const ErrorClass& eclass) {
05580
05581 std::cout << std::endl << std::endl << std::endl;
05582
05583 if(routeCost != NULL){
05584 delete[] routeCost;
05585 routeCost = NULL;
05586 }
05587
05588 if(routeDemand != NULL){
05589 delete[] routeDemand;
05590 routeDemand = NULL;
05591 }
05592
05593 if(xVar != NULL){
05594 delete[] xVar;
05595 xVar = NULL;
05596 }
05597
05598
05599
05600
05601 throw ErrorClass(eclass.errormsg);
05602 }
05603
05604 }
05605
05606
05607
05608
05609
05610 CoinSolver* OSBearcatSolverXij::getMultiCommodInstance(int hubIndex){
05611
05612
05613 int i;
05614 int j;
05615 int numVar;
05616 int numNonHubs;
05617 numNonHubs = m_numNodes - m_numHubs;
05618
05619 numVar = numNonHubs;
05620
05621 numVar += numNonHubs;
05622
05623 numVar += numNonHubs*numNonHubs - numNonHubs;
05624 int numNonz;
05625 int kount;
05626 int numCon;
05627 CoinSolver *solver = NULL;
05628 SparseVector *objcoeff = NULL;
05629
05630 numCon = numNonHubs + (numNonHubs*numNonHubs - numNonHubs) + 1;
05631
05632
05633
05634 OSInstance *osinstance = new OSInstance();
05635 try{
05636
05637 osinstance->setInstanceSource("generated from OSBearcatSoleverXij");
05638 osinstance->setInstanceDescription("Instance for finding a multicommodity cut");
05639 osinstance->setVariableNumber( numVar);
05640
05641 numVar = 0;
05642
05643 for(i = m_numHubs; i < m_numNodes; i++){
05644
05645 if(m_nodeName[ i] != "")
05646 osinstance->addVariable(numVar++, "w[" + m_nodeName[ i] +"]", 0, OSDBL_MAX, 'C');
05647 else
05648 osinstance->addVariable(numVar++, makeStringFromInt("w[", i) +"]", 0, OSDBL_MAX, 'C');
05649
05650 m_tmpVarMap.insert( std::pair<int,std::string>(numVar, "w[" + m_nodeName[ i] +"]" ) );
05651
05652 }
05653
05654
05655
05656 for(j = m_numHubs; j < m_numNodes; j++){
05657
05658 if(m_nodeName[ hubIndex ] != "" && m_nodeName[ j] != "")
05659 osinstance->addVariable(numVar++, "u[" + m_nodeName[ hubIndex] + "," + m_nodeName[ j] +"]", 0, OSDBL_MAX, 'C');
05660 else
05661 osinstance->addVariable(numVar++, makeStringFromInt("u[", hubIndex) + makeStringFromInt(",", j) +"]", 0, OSDBL_MAX, 'C');
05662
05663 m_tmpVarMap.insert( std::pair<int,std::string>(numVar, m_nodeName[ hubIndex] + "," + m_nodeName[ j] +"]") );
05664
05665
05666 }
05667
05668 for(i = m_numHubs; i < m_numNodes; i++){
05669
05670
05671 for(j = m_numHubs; j < i; j++){
05672
05673
05674
05675 if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
05676 osinstance->addVariable(numVar++, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", 0, OSDBL_MAX, 'C');
05677 else
05678 osinstance->addVariable(numVar++, makeStringFromInt("u[", i) + makeStringFromInt(",", j) +"]", 0, OSDBL_MAX, 'C');
05679
05680
05681 m_tmpVarMap.insert( std::pair<int,std::string>(numVar, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]") );
05682
05683 }
05684
05685 for(j = i + 1; j < m_numNodes; j++){
05686
05687 if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
05688 osinstance->addVariable(numVar++, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", 0, OSDBL_MAX, 'C');
05689 else
05690 osinstance->addVariable(numVar++, makeStringFromInt("u[", i) + makeStringFromInt(",", j) +"]", 0, OSDBL_MAX, 'C');
05691
05692 m_tmpVarMap.insert( std::pair<int,std::string>(numVar, "u[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]") );
05693
05694
05695 }
05696
05697 }
05698
05699
05700
05701 osinstance->setObjectiveNumber( 1);
05702
05703
05704
05705 objcoeff = new SparseVector( numVar);
05706
05707 for(i = 0; i < numVar; i++){
05708
05709 objcoeff->indexes[ i] = i;
05710 objcoeff->values[ i] = 0;
05711
05712 }
05713
05714 osinstance->addObjective(-1, "cutViolation", "max", 0.0, 1.0, objcoeff);
05715 objcoeff->bDeleteArrays = true;
05716 delete objcoeff;
05717
05718 osinstance->setConstraintNumber( numCon );
05719
05720 numCon = 0;
05721 for(i = m_numHubs; i < m_numNodes; i++){
05722
05723 if(m_nodeName[ hubIndex] != "" && m_nodeName[ i] != "")
05724 osinstance->addConstraint(numCon++, "dualCon[" + m_nodeName[ hubIndex] + "," + m_nodeName[ i] + "]", -OSDBL_MAX, 0, 0);
05725 else
05726 osinstance->addConstraint(numCon++, makeStringFromInt("dualCon[", hubIndex) + makeStringFromInt(",", i) +"]", -OSDBL_MAX, 0, 0);
05727 }
05728
05729
05730 for(i = m_numHubs; i < m_numNodes; i++){
05731
05732
05733 for(j = m_numHubs; j < i; j++){
05734
05735 if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
05736 osinstance->addConstraint(numCon++, "dualCon[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", -OSDBL_MAX, 0, 0);
05737 else
05738 osinstance->addConstraint(numCon++, makeStringFromInt("dualCon[", i) + makeStringFromInt(",", j) +"]", -OSDBL_MAX, 0, 0);
05739
05740
05741 }
05742
05743 for(j = i + 1; j < m_numNodes; j++){
05744
05745 if(m_nodeName[ i] != "" && m_nodeName[ j] != "")
05746 osinstance->addConstraint(numCon++, "dualCon[" + m_nodeName[ i] + "," + m_nodeName[ j] +"]", -OSDBL_MAX, 0, 0);
05747 else
05748 osinstance->addConstraint(numCon++, makeStringFromInt("dualCon[", i) + makeStringFromInt(",", j) +"]", -OSDBL_MAX, 0, 0);
05749
05750
05751 }
05752
05753 }
05754
05755 double upperBound;
05756 upperBound = numVar - numNonHubs ;
05757 upperBound = 100;
05758 osinstance->addConstraint(numCon++, "boundConstraint", -OSDBL_MAX, upperBound, 0);
05759
05760
05761 numCon = numNonHubs + (numNonHubs*numNonHubs - numNonHubs) + 1;
05762 numNonz = 2*numNonHubs;
05763 numNonz += 3*(numNonHubs*numNonHubs - numNonHubs);
05764 numNonz += (numNonHubs*numNonHubs - numNonHubs) + numNonHubs;
05765
05766
05767 double *values = new double[ numNonz];
05768 int *columnIndexes = new int[ numNonz];
05769
05770 int *starts = new int[ numCon + 1];
05771
05772
05773 kount = 0;
05774 numNonz = 0;
05775 starts[ kount++] = 0;
05776
05777
05779
05780
05781 int uijKount;
05782 uijKount = numNonHubs;
05783
05784 for(j = m_numHubs; j < m_numNodes; j++){
05785
05786
05787
05788 columnIndexes[ numNonz] = j - m_numHubs ;
05789 values[ numNonz++] = 1.0;
05790
05791
05792 columnIndexes[ numNonz] = uijKount ;
05793 values[ numNonz++] = -1.0;
05794
05795 starts[ kount++] = numNonz;
05796 uijKount++;
05797 }
05798
05799
05800 for(i = m_numHubs; i < m_numNodes; i++){
05801
05802
05803 for(j = m_numHubs; j < i; j++){
05804
05805
05806
05807 columnIndexes[ numNonz] = i - m_numHubs ;
05808 values[ numNonz++] = -1.0;
05809
05810
05811 columnIndexes[ numNonz] = j - m_numHubs ;
05812 values[ numNonz++] = 1.0;
05813
05814
05815 columnIndexes[ numNonz] = uijKount ;
05816 values[ numNonz++] = -1.0;
05817
05818
05819 starts[ kount++] = numNonz;
05820 uijKount++;
05821
05822
05823 }
05824
05825 for(j = i + 1; j < m_numNodes; j++){
05826
05827
05828
05829 columnIndexes[ numNonz] = i - m_numHubs ;
05830 values[ numNonz++] = -1.0;
05831
05832
05833 columnIndexes[ numNonz] = j - m_numHubs ;
05834 values[ numNonz++] = 1.0;
05835
05836
05837 columnIndexes[ numNonz] = uijKount ;
05838 values[ numNonz++] = -1.0;
05839
05840
05841 starts[ kount++] = numNonz;
05842 uijKount++;
05843
05844
05845 }
05846
05847 }
05848
05849
05850 for(i = numNonHubs; i < numVar; i++ ){
05851
05852
05853 columnIndexes[ numNonz] = i ;
05854 values[ numNonz++] = 1.0;
05855
05856 }
05857
05858 starts[ kount++] = numNonz;
05859 osinstance->setLinearConstraintCoefficients(numNonz, false, values, 0, numNonz - 1,
05860 columnIndexes, 0, numNonz - 1, starts, 0, numVar);
05861
05862
05863
05864
05865 solver = new CoinSolver();
05866 solver->sSolverName ="clp";
05867 solver->osinstance = osinstance;
05868 solver->buildSolverInstance();
05869 solver->osoption = m_osoption;
05870
05871
05872 return solver;
05873
05874
05875
05876 } catch (const ErrorClass& eclass) {
05877
05878 if( objcoeff != NULL ){
05879 delete objcoeff;
05880 objcoeff = NULL;
05881 }
05882
05883 throw ErrorClass(eclass.errormsg);
05884 }
05885
05886
05887 }
05888
05889
05890
05891 void OSBearcatSolverXij::getFeasibleSolution(){
05892
05893 int i;
05894 int k;
05895 int numVar;
05896 int numNonHubs;
05897 numNonHubs = m_numNodes - m_numHubs;
05898
05899 numVar = m_numHubs*numNonHubs;
05900
05901 int numNonz;
05902 int kount;
05903
05904 CoinSolver *solver = NULL;
05905 SparseVector *objcoeff = NULL;
05906
05907 std::pair<int,int> indexPair1;
05908 std::pair<int,int> indexPair2;
05909
05910 std::map<int, std::vector<int> > routeMap;
05911
05912 OSInstance *osinstance = new OSInstance();
05913 try{
05914
05915 osinstance->setInstanceSource("generated from OSBearcatSoleverXij");
05916 osinstance->setInstanceDescription("Generalized Assignment Instance for finding a feasible solution");
05917 osinstance->setVariableNumber( numVar);
05918
05919 numVar = 0;
05920
05921 for(k = 0; k < m_numHubs; k++){
05922
05923 for(i = m_numHubs; i < m_numNodes; i++){
05924
05925 if( m_nodeName[ k] != "" && m_nodeName[ i] != "")
05926 osinstance->addVariable(numVar++, "x(" + m_nodeName[ k] + "," + m_nodeName[ i] +")", 0, 1, 'B');
05927 else
05928 osinstance->addVariable(numVar++, makeStringFromInt("x(" ,k) + makeStringFromInt(",", i) +")", 0, 1, 'B');
05929
05930 }
05931
05932 }
05933
05934
05935 osinstance->setObjectiveNumber( 1);
05936
05937
05938
05939 objcoeff = new SparseVector( numVar);
05940
05941 kount = 0;
05942 for(k = 0; k < m_numHubs; k++){
05943
05944 indexPair1.first = k;
05945 indexPair2.second = k;
05946
05947 for(i = m_numHubs; i < m_numNodes; i++){
05948
05949 indexPair1.second = i;
05950 indexPair2.first = i;
05951 objcoeff->indexes[ kount] = kount;
05952
05953 if( (m_xVarIndexMap.find( indexPair1) == m_xVarIndexMap.end() ) ||
05954 (m_xVarIndexMap.find( indexPair2) == m_xVarIndexMap.end() ) ){
05955 throw ErrorClass("Index problem in generalized assignment problem to find feasible solution");
05956 }
05957
05958 objcoeff->values[ kount++] = (m_cost[ m_xVarIndexMap[ indexPair1] ] +
05959 m_cost[ m_xVarIndexMap[ indexPair2] ])*m_demand[i] ;
05960
05961 }
05962
05963 }
05964
05965 osinstance->addObjective(-1, "feasibililtyObj", "min", 0.0, 1.0, objcoeff);
05966 objcoeff->bDeleteArrays = true;
05967 delete objcoeff;
05968 objcoeff = NULL;
05969
05970 osinstance->setConstraintNumber( m_numNodes );
05971
05972
05973 for(k = 0; k < m_numHubs; k++){
05974
05975 if(m_nodeName[ k] != "" )
05976 osinstance->addConstraint(k, "capacityCon[" + m_nodeName[ k] + "]", 0, m_routeCapacity[ k], 0);
05977 else
05978 osinstance->addConstraint(k, makeStringFromInt("dualCon[", k) +"]", 0, m_routeCapacity[ k], 0);
05979 }
05980
05981
05982 for(i = m_numHubs; i < m_numNodes; i++){
05983
05984 if(m_nodeName[ i] != "" )
05985 osinstance->addConstraint(i, "assingCon[" + m_nodeName[ i] +"]", 1, 1, 0);
05986 else
05987 osinstance->addConstraint(i, makeStringFromInt("assignCon[", i) +"]", 1, 1, 0);
05988
05989 }
05990
05991
05992
05993
05994 numNonz = 2*numVar;
05995
05996
05997 double *values = new double[ numNonz];
05998 int *rowIndexes = new int[ numNonz];
05999 int *starts = new int[ numVar + 1];
06000
06001
06002 kount = 0;
06003 numNonz = 0;
06004 starts[ kount++] = 0;
06005
06006
06008
06009
06010 for(k = 0; k < m_numHubs; k++){
06011
06012
06013 for(i = m_numHubs; i < m_numNodes; i++){
06014
06015
06016 rowIndexes[ numNonz] = k ;
06017 values[ numNonz++] = m_demand[ i];
06018
06019 rowIndexes[ numNonz] = i ;
06020 values[ numNonz++] = 1.0;
06021
06022 starts[ kount++] = numNonz;
06023
06024 }
06025 }
06026
06027
06028 osinstance->setLinearConstraintCoefficients(numNonz, true, values, 0, numNonz - 1,
06029 rowIndexes, 0, numNonz - 1, starts, 0, numVar);
06030
06031
06032
06033
06034
06035 solver = new CoinSolver();
06036 solver->sSolverName ="cbc";
06037 solver->osinstance = osinstance;
06038 solver->buildSolverInstance();
06039 solver->osoption = m_osoption;
06040 solver->solve();
06041
06042
06043 OSResult *osresult;
06044 osresult = solver->osresult;
06045 std::string solStatus;
06046 std::vector<IndexValuePair*> primalValPair;
06047 int vecSize;
06048 double optSolValue;
06049
06050 solStatus = osresult->getSolutionStatusType( 0 );
06051
06052 if( solStatus.find("ptimal") != std::string::npos ){
06053
06054 optSolValue = osresult->getOptimalObjValue( -1, 0);
06055 std::cout << "OPTIMAL SOLUTION VALUE " << optSolValue << std::endl;
06056 }else{
06057 throw ErrorClass("There is no feasible solution to this problem!");
06058 }
06059
06060 primalValPair = osresult->getOptimalPrimalVariableValues( 0);
06061 vecSize = primalValPair.size();
06062
06063
06064 kount = 0;
06065 for(k = 0; k < m_numHubs; k++){
06066
06067 indexPair1.first = k;
06068 std::cout << std::endl << std::endl;
06069 for(i = m_numHubs; i < m_numNodes; i++){
06070
06071 indexPair1.second = i;
06072 if(primalValPair[ kount ]->value > m_osDecompParam.zeroTol) {
06073
06074 std::cout << "Variable = " << m_variableNames[ m_xVarIndexMap[ indexPair1] ]
06075 << " value = " << primalValPair[ kount ]->value << std::endl;
06076
06077 routeMap[k].push_back( i);
06078 }
06079
06080 kount++;
06081 }
06082
06083
06084 }
06085
06086 m_initSolMap[ 0] = routeMap;
06087 delete solver;
06088 solver = NULL;
06089
06090
06091 } catch (const ErrorClass& eclass) {
06092
06093 if( objcoeff != NULL ){
06094 delete objcoeff;
06095 objcoeff = NULL;
06096 }
06097
06098 throw ErrorClass(eclass.errormsg);
06099 }
06100
06101
06102 }
06103
06104 void OSBearcatSolverXij::getVariableIndexMap(){
06105
06106 int i;
06107 int j;
06108 int kount;
06109 std::pair<int,int> indexPair;
06110
06111 kount = 0;
06112 for(i = 0; i < m_numNodes; i++){
06113
06114 for(j = 0; j < i; j++){
06115
06116 indexPair.first = i;
06117 indexPair.second = j;
06118 m_xVarIndexMap[ indexPair] = kount;
06119 kount++;
06120 }
06121
06122 for(j = i + 1; j < m_numNodes; j++){
06123
06124 indexPair.first = i;
06125 indexPair.second = j;
06126 m_xVarIndexMap[ indexPair] = kount;
06127 kount++;
06128 }
06129 }
06130
06131
06132 }
06133
06134
06135 void OSBearcatSolverXij::permuteHubs(){
06136
06137 int k1;
06138 int k2;
06139
06140 double tmpVal;
06141 double *tmpCap;
06142 int tmpIdx;
06143 tmpCap = new double[ m_numHubs];
06144
06145 for(k1 = 0; k1 < m_numHubs; k1++) tmpCap[ k1] = m_routeCapacity[ k1];
06146 for(k1 = 0; k1 < m_numHubs; k1++) m_hubPoint[ k1] = k1;
06147
06148 for(k1 = 0; k1 < m_numHubs - 1; k1++){
06149
06150
06151 for(k2 = k1 + 1; k2 < m_numHubs; k2++){
06152
06153 if( tmpCap[ k2 ] < tmpCap[ k1 ] ){
06154
06155
06156 tmpVal = tmpCap[ k1 ];
06157 tmpCap[ k1 ] = tmpCap[ k2 ];
06158 tmpCap[ k2 ] = tmpVal;
06159
06160 tmpIdx = m_hubPoint[ k1];
06161 m_hubPoint[ k1] = m_hubPoint[ k2];
06162 m_hubPoint[ k2] = tmpIdx;
06163
06164 }
06165
06166 }
06167 }
06168
06169 for(k1 = 0; k1 < m_numHubs; k1++) std::cout << "m_hubPoint = " << m_hubPoint[ k1] << std::endl;
06170 for(k1 = 0; k1 < m_numHubs; k1++) std::cout << "tmp Cap = " << tmpCap[ k1] << std::endl;
06171 for(k1 = 0; k1 < m_numHubs; k1++) std::cout << "hub capacity = " << m_routeCapacity[ m_hubPoint[ k1] ]<< std::endl;
06172
06173 delete[] tmpCap;
06174 tmpCap = NULL;
06175
06176 }
06177
06178
06179
06180
06181
06182
06183