- Main Page
- Namespaces
- Classes
- Files
- Related Pages

- File List
- File Members

# /home/coin/SVN-release/OS-2.3.5/OS/applications/columnGen/code/OSRouteSolver.cpp

Go to the documentation of this file.

    00001 /* $Id: OSRouteSolver.cpp 3038 2009-11-07 11:43:44Z kmartin $ */
    00019 #include "OSRouteSolver.h"
    00020
    00021 #include "OSErrorClass.h"
    00022 #include "OSDataStructures.h"
    00023
    00024 #include "OSInstance.h"
    00025 #include "OSCoinSolver.h"
    00026 #include "OSConfig.h"
    00027 #include "OSResult.h"
    00028
    00029 #include "OSiLReader.h"
    00030 #include "OSiLWriter.h"
    00031 #include "OSoLReader.h"
    00032 #include "OSoLWriter.h"
    00033 #include "OSrLReader.h"
    00034 #include "OSrLWriter.h"
    00035 #include "OSInstance.h"
    00036 #include "OSFileUtil.h"
    00037
    00038
    00039 #include "ClpSimplex.hpp"
    00040 #include "ClpFactorization.hpp"
    00041 #include "ClpNetworkMatrix.hpp"
    00042 #include "OsiClpSolverInterface.hpp"
    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
    00056
    00057 std::string makeStringFromInt(std::string theString, int theInt);
    00058
    00059
    00060 OSRouteSolver::OSRouteSolver() {
    00061         std::cout << "INSIDE OSRouteSolver CONSTRUCTOR with OSOption argument" << std::endl;
    00062 }//end default OSRouteSolver constructor
    00063
    00064 OSRouteSolver::OSRouteSolver(OSOption *osoption) {
    00065         std::cout << "INSIDE OSRouteSolver CONSTRUCTOR with OSOption argument" << std::endl;
    00066
    00067         m_u = NULL;
    00068         m_v = NULL;
    00069         m_g = NULL;
    00070         m_px = NULL;
    00071         m_tx =NULL;
    00072         m_varIdx = NULL;
    00073
    00074         m_optL = NULL;
    00075         m_optD = NULL;
    00076         m_vv = NULL;
    00077         m_vvpnt = NULL;
    00078
    00079         m_demand = NULL;
    00080         m_cost = NULL;
    00081
    00082         m_rc = NULL;
    00083
    00084         try{
    00085
    00086                 getOptions( osoption);
    00087
    00088                 if(m_routeCapacity > m_totalDemand){
    00089                         m_upperBoundL = m_totalDemand;
    00090                 }else{
    00091                         m_upperBoundL = m_routeCapacity;
    00092                 }
    00093
    00094
    00095
    00096                 //m_varIdx = new int[ m_numNodes];
    00097                 m_varIdx = new int[ m_upperBoundL + 1];
    00098
    00099
    00100                 m_u = new double*[ m_numNodes];
    00101                 m_v = new double*[ m_numNodes];
    00102                 m_g = new double*[ m_numNodes];
    00103
    00104                 m_px = new int*[ m_numNodes];
    00105                 m_tx = new int*[ m_numNodes];
    00106
    00107                 int k;
    00108                 int i;
    00109
    00110                 int l;
    00111
    00112
    00113                 for (i = 0; i < m_numNodes; i++) {
    00114
    00115
    00116                         m_u[ i] = new double[ m_upperBoundL + 1];
    00117                         m_v[ i] = new double[ m_upperBoundL + 1];
    00118
    00119                         for(l = 0; l <= m_upperBoundL; l++){
    00120                                 m_u[ i][l] = OSDBL_MAX;
    00121                                 m_v[ i][l] = OSDBL_MAX;
    00122                         }
    00123
    00124                         m_g[ i] = new double[ m_numNodes];
    00125                         m_px[ i] = new int[ m_upperBoundL + 1];
    00126                         m_tx[ i] = new int[m_upperBoundL + 1];
    00127
    00128
    00129                 }
    00130
    00131
    00132                 //outer dynamic programming arrays
    00133                 m_optL = new int[ m_numHubs];
    00134                 m_optD = new int[ m_numHubs];
    00135
    00136                 m_vv = new double*[ m_numHubs];
    00137                 m_vvpnt = new int*[ m_numHubs];
    00138                 m_cost = new double*[ m_numHubs];
    00139
    00140                 for (k = 0; k < m_numHubs; k++) {
    00141
    00142
    00143                         m_vv[ k] = new double[ m_totalDemand + 1];
    00144                         m_vvpnt[ k] = new int[ m_totalDemand + 1];
    00145                         m_cost[ k] = new double[ m_numNodes*m_numNodes - m_numNodes];
    00146
    00147
    00148                 }
    00149
    00150                 //allocate memory for the reduced cost vector.
    00151                 //assume order is k, l, i, j
    00152                 m_rc = new double[ m_numHubs*m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes)];
    00153
    00154
    00155                 //allocate memory for convexity row
    00156                 m_psi = new double[ m_numHubs];
    00157
    00158                 //allocate memory for node assignment
    00159                 m_phi = new double[ m_numNodes];
    00160
    00161                 m_optValHub = new double[ m_numHubs];
    00162
    00163                 m_variableNames = new string[ m_numNodes*(m_numNodes - 1)];
    00164
    00165                 createVariableNames();
    00166
    00167                 //these are constraints that say we must be incident to each
    00168                 //non-hub node -- there are  m_numNodes - m_numHubs of these
    00169                 m_pntAmatrix = new int[ m_numNodes - m_numHubs + 1];
    00170                 //the variables -- the variable space we are in is x_{ij} NOT
    00171                 // x_{ijk} -- a bit tricky
    00172                 //m_Amatrix[ j] is a variable index -- this logic works
    00173                 //since the Amatrix coefficient is 1 -- we don't need a value
    00174                 //it indexes variable that points into the node
    00175                 m_Amatrix = new  int[ (m_numNodes - m_numHubs)*(m_numNodes - 1) ];
    00176                 createAmatrix();
    00177
    00178                 //this has size of the number of x variables
    00179                 int numVar = m_numNodes*m_numNodes - m_numHubs ;
    00180                 m_tmpScatterArray = new int[ numVar ];
    00181                 for(i = 0; i < numVar; i++){
    00182
    00183                         m_tmpScatterArray[ i] = 0;
    00184
    00185                 }
    00186
    00187                 m_nonzVec = new int[ m_numHubs] ;
    00188                 m_costVec = new double[ m_numHubs];
    00189
    00190
    00191                 m_newColumnRowIdx = new int*[ m_numHubs];
    00192                 m_newColumnRowValue = new double*[ m_numHubs];
    00193
    00194 //kipp change -- put the 1000 in as an option
    00195 //hardcoding
    00197                 for (k = 0; k < m_numHubs; k++) {
    00198
    00199                         m_newColumnRowValue[ k] = new double[ 1000];
    00200                         m_newColumnRowIdx[ k] = new int[ 1000];
    00201
    00202                 }
    00203                 //for now, the number of columns will be 10000
    00204                 //for now number of nonzeros will be 500000
    00205                 m_thetaPnt = new int[ 10000];
    00206                 m_thetaCost = new double[ 10000];
    00207                 m_thetaIndex = new int[ 500000];
    00208                 m_numThetaVar = 0;
    00209                 m_numThetaNonz = 0;
    00210                 m_thetaPnt[ m_numThetaVar++ ] = 0;
    00211
    00212                 //kipp -- move this later
    00213                 getSeparationInstance();
    00214         } catch (const ErrorClass& eclass) {
    00215
    00216                 throw ErrorClass(eclass.errormsg);
    00217
    00218         }
    00219
    00220
    00221 }//end OSRouteSolver Constructor
    00222
    00223
    00224 OSRouteSolver::~OSRouteSolver(){
    00225
    00226         std::cout << "INSIDE ~OSRouteSolver DESTRUCTOR" << std::endl;
    00227
    00228
    00229
    00230         //delete data structures for arrays used in calculating minimum reduced cost
    00231         int i;
    00232
    00233         for(i = 0; i < m_numNodes; i++){
    00234
    00235
    00236
    00237             delete[] m_v[i];
    00238             delete[] m_g[i];
    00239             delete[] m_px[i];
    00240             delete[] m_tx[i];
    00241                 delete[] m_u[i];
    00242
    00243
    00244         }
    00245
    00246         delete[] m_u;
    00247         m_u= NULL;
    00248
    00249         delete[] m_v;
    00250         m_v= NULL;
    00251
    00252         delete[] m_g;
    00253         m_g= NULL;
    00254
    00255         delete[] m_px;
    00256         m_px= NULL;
    00257
    00258         delete[] m_tx;
    00259         m_tx= NULL;
    00260
    00261
    00262
    00263         if(m_demand != NULL){
    00264                 //std::cout << "I AM DELETING m_demand" << std::endl;
    00265                 delete[] m_demand;
    00266         }
    00267
    00268
    00269         if(m_varIdx != NULL) delete[] m_varIdx;
    00270
    00271         for(i = 0; i < m_numHubs; i++){
    00272
    00273                 delete[] m_vv[i];
    00274                 delete[] m_vvpnt[i];
    00275                 delete[] m_cost[ i];
    00276
    00277
    00278         }
    00279         delete[] m_optL;
    00280         m_optL = NULL;
    00281         delete[] m_optD;
    00282         m_optD = NULL;
    00283         delete[] m_vv;
    00284         m_vv = NULL;
    00285         delete[] m_vvpnt;
    00286         m_vvpnt = NULL;
    00287
    00288         delete[] m_cost;
    00289         m_cost = NULL;
    00290
    00291         delete[] m_rc;
    00292         m_rc = NULL;
    00293
    00294         delete[] m_phi;
    00295         m_phi = NULL;
    00296
    00297         delete[] m_psi;
    00298         m_psi = NULL;
    00299
    00300         delete[] m_optValHub;
    00301         m_optValHub = NULL;
    00302
    00303         delete[] m_variableNames;
    00304         m_variableNames = NULL;
    00305
    00306         delete[] m_pntAmatrix;
    00307         m_pntAmatrix = NULL;
    00308
    00309         delete[] m_Amatrix;
    00310         m_Amatrix = NULL;
    00311
    00312         delete[] m_tmpScatterArray;
    00313         m_tmpScatterArray = NULL;
    00314
    00315         delete[] m_nonzVec  ;
    00316         m_nonzVec = NULL;
    00317         delete[] m_costVec ;
    00318         m_costVec = NULL;
    00319
    00320         for(i = 0; i < m_numHubs; i++){
    00321
    00322                 delete[] m_newColumnRowIdx[i];
    00323                 delete[] m_newColumnRowValue[i];
    00324         }
    00325
    00326         delete[] m_newColumnRowIdx;
    00327         m_newColumnRowIdx = NULL;
    00328
    00329         delete[] m_newColumnRowValue;
    00330         m_newColumnRowValue = NULL;
    00331
    00332         delete[] m_thetaPnt;
    00333         m_thetaPnt = NULL;
    00334
    00335         delete[] m_thetaIndex;
    00336         m_thetaIndex = 0;
    00337
    00338
    00339         delete[] m_thetaCost;
    00340         m_thetaCost = 0;
    00341
    00342
    00343 }//end ~OSRouteSolver
    00344
    00345
    00346
    00347
    00348
    00349
    00350
    00351 void OSRouteSolver::getOptL(const  double* c) {
    00352
    00353         //initialize the first HUB
    00354
    00355         int k;
    00356         int d;
    00357         int d1;
    00358         int kountVar;
    00359         double testVal;
    00360         int l;
    00361         //int startPntInc;
    00362         double trueMin;
    00363
    00364         kountVar = 0;
    00365         //startPntInc = m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes);
    00366
    00367         m_vv[ 0][ 0] = 0;
    00368         for(d = 1; d <=  m_totalDemand; d++){
    00369
    00370                 m_vv[ 0][ d] = OSDBL_MAX;
    00371
    00372         }
    00373
    00374         //now loop over the other HUBS
    00375
    00376
    00377
    00378         for(k = 1; k < m_numHubs; k++){
    00379
    00380                 for(d = 1; d <= m_totalDemand; d++){
    00381
    00382                         m_vv[ k][ d] = OSDBL_MAX;
    00383
    00384                         //d1 is the state variable at stage k -1
    00385                         for(d1 = 0; d1 <= m_totalDemand; d1++){
    00386
    00387                                 l = d - d1;
    00388
    00389                                 if( (m_vv[ k - 1][ d1] < OSDBL_MAX) &&  (l <= m_upperBoundL) && (l >= 1) ){
    00390
    00391
    00392                                         // l was the decision at state d1 in stage k-1
    00393                                         // l + d1 brings us to state d at stage k
    00394                                         // d is the total carried on routes 0 -- k-1
    00395
    00396                                         testVal = qrouteCost(k - 1,  l,  c,  &kountVar);
    00397
    00398                                         //std::cout << "L = " << l << std::endl;
    00399                                         //std::cout << "testVal " << testVal << std::endl;
    00400
    00401                                         if( m_vv[ k-1][ d1]  +  testVal < m_vv[  k][ d] ){
    00402
    00403                                                 m_vv[ k][ d] =  m_vv[ k-1][ d1]  +  testVal;
    00404                                                 //now point to the best way to get to d
    00405                                                 m_vvpnt[ k][ d]  = d1;
    00406
    00407                                         }
    00408
    00409
    00410                                 }
    00411
    00412                         }
    00413
    00414                 }
    00415
    00416                 //c+=  startPntInc ;
    00417
    00418         }//  end for on k
    00419
    00420         trueMin = OSDBL_MAX;
    00421         //we now enter the last stage through the other hubs
    00422         // have satisfied total demand d
    00423
    00424         int  dlower = 0;
    00425         if (m_numHubs > 1) dlower = 1;
    00426
    00427         for(d = dlower; d < m_totalDemand; d++){
    00428
    00429                 //std::cout << "m_vv[ m_numHubs - 1 ][ d]  " << m_vv[ m_numHubs - 1 ][ d]  << std::endl;
    00430                 l = m_totalDemand - d;
    00431
    00432                 if(m_vv[ m_numHubs - 1 ][ d]  < OSDBL_MAX  && l <= m_upperBoundL && l >= 1){
    00433
    00434
    00435
    00436                         //std::cout << "LL = " <<  l  << std::endl;
    00437
    00438
    00439                         testVal = qrouteCost(m_numHubs -1 ,  l,  c,  &kountVar);
    00440
    00441                         //std::cout << "l = " << l << std::endl;
    00442                         //std::cout << "testVal = " << testVal << std::endl;
    00443
    00444                         if(m_vv[ m_numHubs - 1][ d] + testVal < trueMin){
    00445
    00446                                 trueMin = m_vv[ m_numHubs -1][ d] + testVal;
    00447                                 m_optD[  m_numHubs -1 ] = d;
    00448                                 m_optL[  m_numHubs -1 ] = l;
    00449
    00450                         }
    00451
    00452
    00453                 }
    00454         }
    00455
    00456         std::cout << "TRUE MIN = " <<  trueMin << std::endl;
    00457
    00458         k = m_numHubs -1;
    00459
    00460         while( k - 1 >= 0) {
    00461
    00462                 m_optD[  k - 1 ] = m_vvpnt[ k][ m_optD[  k  ] ];
    00463
    00464                 m_optL[ k - 1 ] =  m_optD[  k  ] - m_optD[  k - 1 ] ;
    00465
    00466                 //std::cout << "k = " <<  k << std::endl;
    00467                 //std::cout << "m_optD[  k  ]  = " <<  m_optD[  k  ] << std::endl;
    00468                 //std::cout << "m_optD[  k -1 ] " << m_optD[  k - 1 ]  << std::endl;
    00469
    00470                 k--;
    00471
    00472
    00473         }
    00474
    00475 }//end getOptL
    00476
    00477
    00478
    00479
    00480
    00481
    00482 double OSRouteSolver::qrouteCost(const int& k, const int& l, const double* c, int* kountVar){
    00483
    00484         //critical -- nodes 0, ..., m_numNodes - 1 are the hub nodes
    00485         // we are doing the calculation for hub k, k <= m_numNodes - 1
    00486         //l  should be the total demand on the network -- we are NOT using 0 based counting
    00487         double rcost;
    00488         rcost = OSDBL_MAX;
    00489
    00490
    00491
    00492         if(l < 0){
    00493
    00494                 std::cout  << "LVALUE  NEGATIVE " << l  << std::endl;
    00495                 exit( 1);
    00496         }
    00497
    00498
    00499
    00500         if(l > m_upperBoundL){
    00501
    00502                 std::cout  << "LVALUE  BIGGER THAN UPPER BOUND " << l  << std::endl;
    00503                 exit( 1);
    00504         }
    00505
    00506         //start of the cost vector for hub k plus move to start of l demands
    00507         // now move the pointer up
    00508         int startPnt = k*m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes) + (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
    00509
    00510         //int startPnt = (l - 1)*(m_numNodes*m_numNodes - m_numNodes);
    00511         c+=  startPnt ;
    00512
    00513
    00514
    00515         *kountVar = 0;
    00516         int bestLastNode;
    00517         int i;
    00518         int j;
    00519         int l1;
    00520         int l2;
    00521         //for(i = 0; i < 20; i++){
    00522         //      std::cout << "i =  " << i  <<  " c[i] = " << *(c + i) << std::endl ;
    00523         //}
    00524
    00525
    00526
    00527         // l  is the total demand on this network
    00528         //address of node (j, i) is j*(m_numNodes-1) + i when i < j
    00529         //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
    00530
    00531         //
    00532         // initialize
    00533
    00534         for(i = m_numHubs; i < m_numNodes; i++){
    00535
    00536
    00537                 for(l1 = m_minDemand; l1 <= l; l1++){  //l-1  is total demand on network
    00538
    00539                         m_u[i][l1] = OSDBL_MAX;
    00540                         m_v[i][l1] = OSDBL_MAX;
    00541                         m_px[i][l1] = -1; //a node we don't have
    00542                         if(l1 == *(m_demand + i) ){
    00543
    00544                                 m_px[i][l1] = k;
    00545                                 // want the cost for arc (k, i)
    00546                                 // note: k < i so use that formula
    00547                                 m_u[i][l1] = *(c + k*(m_numNodes - 1) + i - 1);  // put in correct cost
    00548
    00549
    00550                         }
    00551                 }
    00552         }
    00553         //end initialize
    00554
    00555
    00556         //
    00557
    00558         //if l = minDemand we visit only one node, let's get the reduced cost in this case
    00559         if(l == m_minDemand){
    00560
    00561                 for(i = m_numHubs; i < m_numNodes; i++){
    00562
    00563
    00564                         if(  m_u[i][l] + *(c + i*(m_numNodes-1) + k )  < rcost){
    00565
    00566                                 rcost = m_u[i][l] + *(c + i*(m_numNodes-1) + k );
    00567
    00568                                 //std::cout << " m_u[i][l2] = "  << m_u[i][l2] << std::endl;
    00569
    00570                                 //std::cout << " *(c + i*(m_numNodes-1) + k ) = "  << *(c + i*(m_numNodes-1) + k ) << std::endl;
    00571                                 //push node back
    00572                                 bestLastNode = i;
    00573                         }
    00574
    00575                 }
    00576
    00577                 //go from node bestLastNode to node k
    00578                 //node bestLastNode is a higher number than k
    00579                 *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1)  +  k ;
    00580                 *(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1)  + bestLastNode - 1;
    00581
    00582
    00583                 return rcost;
    00584         }//end if on l == minDemand
    00585
    00586
    00587
    00588
    00589 // now calculate values for demand 2 or greater
    00590         //address of node (j, i) is j*(m_numNodes-1) + i when i < j
    00591         //address of node (j, i) is j*(m_numNodes-1) + i - 1 when i > j
    00592         // we start l2 at 2 since demand must be at least 1
    00593         // change to min demand + 1
    00594         int lowerVal = m_minDemand + 1;
    00595         for(l2 = lowerVal; l2 <= l; l2++){// loop over possible demand values assuming we have already gone to at least one node
    00596
    00597                 for(i = m_numHubs; i < m_numNodes; i++) { //we are finding least cost to node i
    00598
    00599
    00600                         if( *(m_demand + i) < l2 ){ // kipp < OR <= ????
    00601
    00602                                 for(j = m_numHubs; j < i; j++){ //we are coming from node j
    00603
    00604
    00605
    00606                                         //calculate g  -- l2 - d[ i]  is the demand trough and including node j
    00607                                         l1 = l2 - *(m_demand + i);
    00608
    00609                                         if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
    00610
    00611
    00612                                                 m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i - 1) ;
    00613
    00614
    00615
    00616
    00617                                         }else{
    00618
    00619                                                 m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i - 1) ;
    00620
    00621
    00622
    00623                                         }
    00624
    00625                                         // update u and the pointer for p
    00626
    00627                                         if(m_g[j][i] < m_u[i][l2] ){
    00628
    00629                                                 m_u[i][l2] = m_g[j][i];
    00630                                                 m_px[i][l2] =  j;
    00631
    00632                                         }
    00633
    00634
    00635
    00636                                 }//end of first for loop on j, j < i
    00637
    00638
    00639                                 for(j = i + 1; j < m_numNodes; j++){ //we are coming from node j
    00640
    00641
    00642                                         //calculate g  -- l2 - d[ i]  is the demand trough and including node j
    00643                                         l1 = l2 - *(m_demand + i);
    00644
    00645                                         if( m_px[j][ l1 ] != i ){//check to make sure we did not come into j from i
    00646
    00647
    00648                                                 m_g[j][i] = m_u[ j][ l1 ] + *(c + j*(m_numNodes-1) + i ) ;
    00649
    00650
    00651                                         }else{
    00652
    00653                                                 m_g[j][i] = m_v[ j][ l1] + *(c + j*(m_numNodes-1) + i ) ;
    00654
    00655                                         }
    00656
    00657                                         // update u and the pointer for p
    00658
    00659                                         if(m_g[j][i] < m_u[i][l2] ){
    00660
    00661                                                 m_u[i][l2] = m_g[j][i];
    00662                                                 m_px[i][l2] =  j;
    00663
    00664                                         }
    00665
    00666
    00667                                 }//end of second for loop on j, j > i
    00668
    00669
    00670                                 //now calculate the second best solution and point
    00671                                 //right now m_px[i][l2] points to the best j node
    00672
    00673                                 for(j =m_numHubs; j < m_numNodes; j++){ //we are coming from node j
    00674
    00675                                         if(j != i){
    00676
    00677                                                 if( (m_g[j][i] < m_v[i][l2] ) && (m_px[i][l2] != j)  ){ // kipp the && gives the second best
    00678
    00679                                                         //if(g(j, i) < v(i, l2),
    00680
    00681                                                         m_v[i][l2] = m_g[j][i];
    00682                                                         m_tx[i][l2] = j;
    00683
    00684
    00685                                                 }
    00686
    00687                                         }
    00688
    00689
    00690                                 }//end second best calculation, anothe for loop on j
    00691
    00692                                 //now if l2 = l  we are done
    00693                                 if(l2 == l ){
    00694
    00695                                         if(  m_u[i][l2] + *(c + i*(m_numNodes-1) + k )  < rcost){
    00696
    00697                                                 rcost = m_u[i][l2] + *(c + i*(m_numNodes-1) + k );
    00698
    00699                                                 bestLastNode = i;
    00700                                         }
    00701
    00702                                 }
    00703
    00704
    00705                         }//end if on demand less than l2
    00706
    00707
    00708                 }//i loop
    00709
    00710
    00711         }//l2 loop
    00712
    00713
    00714         //std::cout << "best Last Node = "  << bestLastNode << std::endl;
    00715
    00716         // now get the path that gives the reduced cost
    00717
    00718
    00719         int currentNode;
    00720         int successorNode;
    00721         int lvalue;
    00722
    00723         //initialize
    00724         // we work backwords from bestLastNode
    00725         //in our recursion we recurse on the currentNode and assume
    00726         //it is in the optimal path
    00727
    00728         //node bestLastNode is a higher number than k
    00729         *(m_varIdx + (*kountVar)++) = startPnt + bestLastNode*(m_numNodes - 1)  +  k ;
    00730
    00731
    00732
    00733
    00734         //by successor, I mean node after current node in the path
    00735         successorNode = k;
    00736         currentNode = bestLastNode;
    00737         //the lvalue is the demand through the currentNode
    00738         lvalue = l ;
    00739
    00740
    00741         while(currentNode != k){
    00742
    00743                 if( m_px[ currentNode][ lvalue ] != successorNode){
    00744
    00745
    00746
    00747                         //update nodes
    00748                         successorNode = currentNode;
    00749                         currentNode = m_px[ currentNode][ lvalue ];
    00750
    00751
    00752                         if(currentNode - successorNode > 0){
    00753                                  //below diagonal
    00754
    00755                                 *(m_varIdx  + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1)  +  successorNode;
    00756
    00757
    00758                         }else{
    00759                                  //above diagonal
    00760
    00761                                 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1)  +  successorNode  - 1 ;
    00762
    00763                         }
    00764
    00765
    00766                 }else{ //take second best
    00767
    00768
    00769                         //update nodes
    00770                         successorNode = currentNode;
    00771                         currentNode = m_tx[ currentNode][ lvalue ];
    00772
    00773                         if(currentNode - successorNode > 0){
    00774                                  //below diagonal
    00775                                 *(m_varIdx  + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1)  +  successorNode;
    00776
    00777                         }else{
    00778                                  //above diagonal
    00779                                 *(m_varIdx + (*kountVar)++) = startPnt + currentNode*(m_numNodes - 1)  +  successorNode  - 1 ;
    00780
    00781                         }
    00782
    00783                 }
    00784
    00785                 //update lvalue
    00786                 lvalue = lvalue - *(m_demand + successorNode);
    00787
    00788
    00789
    00790         }//end while
    00791
    00792
    00793         //go from node k to bestLastNode -- already done in loop above
    00794         //*(m_varIdx + (*kountVar)++) = startPnt + k*(m_numNodes - 1)  +  currentNode - 1;
    00795
    00796
    00797         return rcost;
    00798
    00799 }//end qroute
    00800
    00801
    00802
    00803
    00804 void OSRouteSolver::getColumns(const  double* y, const int numRows,
    00805                 int &numColumns, int* numNonzVec, double* costVec, double* rcostVec,
    00806                 int** rowIdxVec, double** valuesVec, double &lowerBound)
    00807 {
    00808
    00809 //first strip of the phi dual values and then the convexity row costs
    00810
    00811         int i;
    00812         int j;
    00813         int numCoulpingConstraints;
    00814         numCoulpingConstraints = m_numNodes - m_numHubs;
    00815
    00816         int numVar;
    00817         numVar = m_numNodes*m_numNodes - m_numHubs;
    00818         int numNonz;
    00819
    00820         try{
    00821
    00822
    00823
    00824                 if(numRows != m_numNodes) throw ErrorClass("inconsistent row count in getColumns");
    00825
    00826                 for(i = 0; i < m_numHubs; i++){
    00827
    00828                         m_phi[i] = 0.0;
    00829                 }
    00830
    00831                 for(i = 0; i < numCoulpingConstraints; i++){
    00832
    00833                         m_phi[ i + m_numHubs] = y[ i];
    00834
    00835                 }
    00836
    00837                 for(i = 0; i < m_numHubs; i++){
    00838
    00839                         m_psi[ i ] = y[ i + numCoulpingConstraints];
    00840
    00841                 }
    00842
    00843                 //get the reduced costs
    00844                 calcReducedCost( m_cost,  m_phi,  m_rc);
    00845
    00846
    00847
    00848                 int kountVar;
    00849                 double testVal;
    00850                 testVal = 0;
    00851                 int k;
    00852                 int startPntInc;
    00853                 int rowCount;
    00854
    00855
    00856
    00857                 getOptL( m_rc);
    00858                 m_lowerBnd = 0.0;
    00859                 for(k = 0; k < m_numHubs; k++){
    00860
    00861                         startPntInc  =  k*m_upperBoundL*(m_numNodes*m_numNodes - m_numNodes) + (m_optL[ k] - 1)*(m_numNodes*m_numNodes - m_numNodes);
    00862
    00863                         std::cout << " whichBlock =  " << k << "  L = " << m_optL[ k] << std::endl;
    00864
    00865                         testVal += m_optL[ k];
    00866
    00867                         kountVar = 0;
    00868
    00869                         m_optValHub[ k] = qrouteCost(k,  m_optL[ k], m_rc,  &kountVar);
    00870
    00871                         m_optValHub[ k] -= m_psi[ k ];
    00872
    00873
    00874
    00875                         std::cout << "Best Reduced Cost Hub " << k << " =  "  << m_optValHub[ k] << std::endl;
    00876                         m_lowerBnd += m_optValHub[ k];
    00877
    00878                         //loop over the rows, scatter each row and figure
    00879                         //out the column coefficients in this row
    00880                         //first scatter the sparse array m_varIdx[ j]
    00881
    00882                         m_costVec[ k] = 0.0;
    00883
    00884                         for(j = 0; j < kountVar; j++){
    00885
    00886
    00887                                 //we are counting the NUMBER of times the variable used
    00888                                 //the same variable can appear more than once in m_varIdx
    00889                                 m_tmpScatterArray[ m_varIdx[ j] - startPntInc  ] += 1;
    00890
    00891                                 // is variable m_varIdx[ j] - startPntInc in this row
    00892
    00893                                 m_costVec[ k] += m_cost[k][  m_varIdx[ j] - startPntInc  ];
    00894
    00895                         }
    00896
    00897
    00898
    00899                         numNonz = 0;
    00900                         //multiply the sparse array by each constraint
    00901                         for(i = 0; i < numCoulpingConstraints; i++){
    00902
    00903                                 rowCount = 0;
    00904
    00905                                 for(j = m_pntAmatrix[ i]; j < m_pntAmatrix[ i + 1]; j++){
    00906
    00907                                         //m_Amatrix[ j] is a variable index -- this logic works
    00908                                         //since the Amatrix coefficient is 1 -- we don't need a value
    00909                                         //it indexes variable that points into the node
    00910                                         rowCount += m_tmpScatterArray[  m_Amatrix[ j] ];
    00911
    00912
    00913                                 }
    00914
    00915                                 if(rowCount > 0){
    00916
    00917                                         //std::cout << "GAIL NODE " << i + m_numHubs <<  " COUNT = " << rowCount << std::endl;
    00918                                         m_newColumnRowIdx[ k][ numNonz] = i;
    00919                                         m_newColumnRowValue[ k][ numNonz] = rowCount;
    00920                                         numNonz++;
    00921                                 }
    00922
    00923
    00924                         }//end loop on coupling constraints
    00925
    00926                         //add a 1 in the convexity row
    00927                         m_newColumnRowIdx[ k][ numNonz] = m_numNodes - m_numHubs + k;
    00928                         m_newColumnRowValue[ k][ numNonz] = 1.0;
    00929                         numNonz++;
    00930
    00931                         m_nonzVec[ k] = numNonz;
    00932
    00933                         //zero out the scatter vector and store the generated column
    00934                         for(j = 0; j < kountVar; j++){
    00935
    00936
    00937                                 m_thetaIndex[ m_numThetaNonz++ ] =  m_varIdx[ j] - startPntInc ;
    00938
    00939                                 m_tmpScatterArray[ m_varIdx[ j] - startPntInc  ]  = 0;
    00940
    00941                                 // is variable m_varIdx[ j] - startPntInc in this row
    00942
    00943                         }
    00944
    00945                         m_costVec[ k] =  m_optL[ k]*m_costVec[ k];
    00946                         m_thetaCost[ m_numThetaVar ] = m_costVec[ k];
    00947                         m_thetaPnt[ m_numThetaVar++ ]  = m_numThetaNonz;
    00948
    00949                         //std::cout << "GAIL COST VECTOR " <<  m_costVec[ k] << std::endl;
    00950
    00951                         //stuff for debugging
    00952                         //*****************************//
    00985                         //**************************//
    00986                         //end debugging stuff
    00987
    00988
    00989                 }//end of loop on hubs
    00990
    00991                 std::cout << "Lower Bound = " << m_lowerBnd << std::endl;
    00992
    00993
    00994                 if(testVal != m_totalDemand) {
    00995
    00996                         std::cout  << "TOTAL DEMAND = " << m_totalDemand << std::endl;
    00997                         std::cout  << "Test Value = " << testVal << std::endl;
    00998                         throw  ErrorClass( "inconsistent demand calculation" );
    00999                 }
    01000
    01001
    01002
    01003
    01004
    01005
    01006         } catch (const ErrorClass& eclass) {
    01007
    01008                 throw ErrorClass(eclass.errormsg);
    01009
    01010         }
    01011
    01012
    01013         //set method parameters
    01014         numColumns = m_numHubs;
    01015         lowerBound =  m_lowerBnd;
    01016
    01017         std::cout << "LEAVING GET COLUMNS" << std::endl;
    01018         return;
    01019
    01020 }//end getColumns
    01021
    01022
    01023
    01024
    01025 OSInstance* OSRouteSolver::getInitialRestrictedMaster( ){
    01026
    01027
    01028         std::cout << "Executing OSRouteSolver::getInitialRestrictedMaster( )" << std::endl;
    01029
    01030         // define the classes
    01031         FileUtil *fileUtil = NULL;
    01032         OSiLReader *osilreader = NULL;
    01033         DefaultSolver *solver  = NULL;
    01034         OSInstance *osinstance = NULL;
    01035
    01036
    01037         // end classes
    01038
    01039         std::string testFileName;
    01040         std::string osil;
    01041
    01042         //std::vector< int> indexes;
    01043         fileUtil = new FileUtil();
    01044
    01045
    01046         std::map<int, std::map<int, std::vector<int> > >::iterator  mit;
    01047         std::map<int, std::vector<int> >::iterator  mit2;
    01048         std::vector<int>::iterator  vit;
    01049
    01050         m_osinstanceMaster = NULL;
    01051
    01052         //add linear constraint coefficients
    01053         //number of values will nodes.size() the coefficients in the node constraints
    01054         //plus coefficients in convexity constraints which is the number of varaibles
    01055         int kountNonz;
    01056         int kount;
    01057         int numThetaVar = m_numberOfSolutions*m_numHubs;
    01058         double *values = new double[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar];
    01059         int *indexes = new int[ m_numberOfSolutions*(m_numNodes-m_numHubs) + numThetaVar];
    01060         int *starts = new int[ numThetaVar + 1];
    01061         kount = 0;
    01062
    01063         starts[ 0] = 0;
    01064
    01065         int startsIdx;
    01066         startsIdx = 0;
    01067
    01068         std::vector<IndexValuePair*> primalValPair;
    01069
    01070         try {
    01071
    01072                 if(m_initOSiLFile.size() == 0) throw ErrorClass("OSiL file to generate restricted master missing");
    01073                 osil = fileUtil->getFileAsString( m_initOSiLFile.c_str());
    01074
    01075                 osilreader = new OSiLReader();
    01076                 osinstance = osilreader->readOSiL(osil);
    01077
    01078
    01079
    01080                 int i;
    01081                 int j;
    01082                 int k;
    01083                 //fill in the cost vector first
    01084                 //the x vector starts at 2*m_numHubs
    01085
    01086                 int idx1;
    01087                 int idx2;
    01088
    01089
    01090                 idx2 = 0;  //zRouteDemand have 0 coefficients in obj
    01091                 //fill in m_cost from the cost coefficient in osil
    01092                 for(k = 0; k < m_numHubs; k++){
    01093
    01094                         idx1 = 0;
    01095
    01096                         for(i = 0; i < m_numNodes; i++){
    01097
    01098                                 for(j = 0; j < i; j++){
    01099
    01100                                         m_cost[k][idx1++ ] = osinstance->instanceData->objectives->obj[0]->coef[ idx2++ ]->value;
    01101                                 }
    01102
    01103                                 for(j = i + 1; j < m_numNodes; j++){
    01104
    01105                                         m_cost[k][idx1++ ] = osinstance->instanceData->objectives->obj[0]->coef[ idx2++ ]->value;
    01106
    01107                                 }
    01108                         }
    01109                 }
    01110
    01111
    01112
    01113                 //get variable names for checking purposes
    01114                 std::string* varNames;
    01115                 varNames =  osinstance->getVariableNames();
    01116
    01117
    01118                 //start building the restricted master here
    01119
    01120
    01121                 m_osinstanceMaster = new OSInstance();
    01122                 m_osinstanceMaster->setInstanceDescription("The Initial Restricted Master");
    01123
    01124                 // first the variables
    01125                 m_osinstanceMaster->setVariableNumber( m_numberOfSolutions*m_numHubs);
    01126
    01127                 // now add the objective function
    01128                 m_osinstanceMaster->setObjectiveNumber( 1);
    01129                 SparseVector *objcoeff = new SparseVector( m_numberOfSolutions*m_numHubs);
    01130
    01131
    01132                 // now the constraints
    01133                 m_osinstanceMaster->setConstraintNumber( m_numNodes);
    01134
    01135
    01136                 //addVariable(int index, string name, double lowerBound, double upperBound, char type, double init, string initString);
    01137                 // we could use setVariables() and add all the variable with one method call -- below is easier
    01138
    01139                 //
    01140
    01141                 int varNumber;
    01142                 varNumber = 0;
    01143                 std::string masterVarName;
    01144                 kountNonz = 0;
    01145                 // now get the primal solution
    01146                 //solve the model for solution in the osoption object
    01147                 for ( mit = m_initSolMap.begin() ; mit != m_initSolMap.end(); mit++ ){
    01148
    01149                         //kipp change upper and lower bounds here on z variables
    01150                         //loop over nodes and routes and set bound
    01151                         //set kount to the start of the z variables
    01152                         //go past the x variables
    01153                         kount  =  2*m_numHubs + m_numHubs*(m_numNodes*m_numNodes - m_numNodes);
    01154                         osinstance->bVariablesModified = true;
    01155                         for ( mit2 = mit->second.begin() ; mit2 != mit->second.end(); mit2++ ){ //we are looping over routes in solution mit
    01156
    01157
    01158
    01159                                 startsIdx++;
    01160                                 starts[ startsIdx ] = kountNonz + mit2->second.size() + 1; //the +1 comes from the convexity row
    01161
    01162
    01163                                 //make sure all lower bounds on z variables on this route back to 0.0
    01164                                 for(i = 0; i < m_numNodes; i++){
    01165                                         osinstance->instanceData->variables->var[ kount + mit2->first*m_numNodes + i]->lb = 0.0;
    01166                                 }
    01167
    01168
    01169                                 for ( vit = mit2->second.begin() ; vit != mit2->second.end(); vit++ ){
    01170
    01171
    01172                                         osinstance->instanceData->variables->var[ kount + mit2->first*m_numNodes + *vit]->lb = 1.0;
    01173                                         //std::cout << "FIXING LOWER BOUND ON VARIABLE " << osinstance->getVariableNames()[ kount + mit2->first*m_numNodes + *vit ] << std::endl;
    01174
    01175                                         values[ kountNonz] = 1.0;
    01176                                         indexes[ kountNonz ] = *vit - m_numHubs ;  //0 based counting
    01177                                         kountNonz++;
    01178
    01179                                 }
    01180
    01181
    01182                                  //now for the convexity row coefficient
    01183                                 values[ kountNonz] = 1;
    01184                                 indexes[ kountNonz ] = m_numNodes - m_numHubs +  mit2->first  ;
    01185                                 kountNonz++;
    01186
    01187
    01188                         }
    01189
    01190                         solver = new CoinSolver();
    01191                         solver->sSolverName ="cbc";
    01192                         solver->osinstance = osinstance;
    01193
    01194                         solver->buildSolverInstance();
    01195                         solver->solve();
    01196
    01197                         //get the solver solution status
    01198
    01199                         std::cout << "Solution Status =  " << solver->osresult->getSolutionStatusType( 0 ) << std::endl;
    01200
    01201                         //get the optimal objective function value
    01202
    01203
    01204                         primalValPair = solver->osresult->getOptimalPrimalVariableValues( 0);
    01205
    01206                         //loop over routes again to create master objective coefficients
    01207
    01208                         for(k = 0; k < m_numHubs; k++){
    01209
    01210
    01211                                 //lets get the x variables
    01212                                 //the variables for this route should be from 2*numHubs + k*(numNodes - 1*)*(numNodes - 1)
    01213                                 idx1 = 2*m_numHubs + k*m_numNodes*(m_numNodes - 1);
    01214                                 idx2 = idx1 + m_numNodes*(m_numNodes - 1);
    01215                                 //end of x variables
    01216
    01217                                 //std::cout << "HUB " <<  k  << " VARIABLES" << std::endl;
    01218
    01219
    01220
    01221                                 for(i = idx1; i < idx2; i++){
    01222                                         if(  primalValPair[ i]->value > .1 ){
    01223                                                 //std::cout <<  osinstance->getVariableNames()[  primalValPair[ i]->idx  ] << std::endl;
    01224                                                 //std::cout <<  m_variableNames[  primalValPair[ i]->idx  -  k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs  ] << std::endl;
    01225                                                 m_thetaIndex[ m_numThetaNonz++ ] = primalValPair[ i]->idx  -  k*(m_numNodes - 1)*m_numNodes - 2*m_numHubs;
    01226                                         }
    01227
    01228                                 }
    01229
    01230                                 m_thetaCost[ m_numThetaVar ] = primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value;
    01231                                 m_thetaPnt[ m_numThetaVar++ ] = m_numThetaNonz;
    01232
    01233                                 masterVarName = makeStringFromInt("theta(", mit->first);
    01234                                 masterVarName += makeStringFromInt(",", k);
    01235                                 masterVarName += ")";
    01236                                 std::cout << masterVarName << std::endl;
    01237                                 m_osinstanceMaster->addVariable(varNumber++, masterVarName, 0, 1, 'C');
    01238
    01239                                 std::cout << "Optimal Objective Value = " << primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value << std::endl;
    01240
    01241                                 objcoeff->indexes[ k + (mit->first)*m_numHubs] = k +  (mit->first)*m_numHubs;
    01242                                 objcoeff->values[ k + (mit->first)*m_numHubs] = primalValPair[ k]->value*primalValPair[ k + m_numHubs]->value;
    01243
    01244
    01245                                 std::cout <<  osinstance->getVariableNames()[ k ] << std::endl;
    01246                                 std::cout <<  osinstance->getVariableNames()[ k + m_numHubs ] << std::endl;
    01247
    01248
    01249                         }//end for on k -- hubs
    01250
    01251
    01252                         primalValPair.clear();
    01253                         delete solver;
    01254                         solver = NULL;
    01255                 }//end for on number of solutions
    01256
    01257                 //add the constraints
    01258                 //add the row saying we must visit each node
    01259                 for( i =  0; i < m_numNodes - m_numHubs ; i++){
    01260
    01261                         m_osinstanceMaster->addConstraint(i,  makeStringFromInt("visitNode_", i + m_numHubs) , 1.0, 1.0, 0);
    01262                 }
    01263
    01264                 kount = 0;
    01265
    01266                 //add the convexity row
    01267                 for( i =  m_numNodes - m_numHubs; i < m_numNodes ; i++){
    01268
    01269                         m_osinstanceMaster->addConstraint(i,  makeStringFromInt("convexityRowRoute_", kount++ ) , 1.0, 1.0, 0);
    01270                 }
    01271
    01272                 m_osinstanceMaster->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
    01273
    01274                 std::cout << "kountNonz = " << kountNonz << std::endl;
    01275
    01276                 //add the linear constraints coefficients
    01277                 m_osinstanceMaster->setLinearConstraintCoefficients(kountNonz , true,
    01278                                 values, 0, kountNonz - 1,  indexes, 0, kountNonz - 1, starts, 0, startsIdx);
    01279
    01280
    01281
    01282                 delete objcoeff;
    01283
    01284                 //delete[] values;
    01285                 //delete[] starts;
    01286                 //delete[] indexes;
    01287                 delete osilreader;
    01288                 osilreader = NULL;
    01289
    01290
    01291
    01292         } catch (const ErrorClass& eclass) {
    01293                 std::cout << std::endl << std::endl << std::endl;
    01294                 if (osilreader != NULL)
    01295                         delete osilreader;
    01296                 if (solver != NULL)
    01297                         delete solver;
    01298
    01299
    01300                 //  Problem with the parser
    01301                 throw ErrorClass(eclass.errormsg);
    01302         }
    01303
    01304         delete fileUtil;
    01305         fileUtil = NULL;
    01306
    01307         return m_osinstanceMaster;
    01308 }//end generateInitialRestrictedMaster
    01309
    01310
    01311
    01312 void OSRouteSolver::getOptions(OSOption *osoption) {
    01313
    01314
    01315         std::cout << "Executing getOptions(OSOption *osoption)" << std::endl;
    01316         //get any options relevant to OSColGenApp
    01317         try{
    01318                 std::vector<SolverOption*> solverOptions;
    01319                 std::vector<SolverOption*>::iterator vit;
    01320                 std::vector<int >demand;
    01321
    01322                 m_numberOfSolutions = 0;
    01323                 solverOptions = osoption->getSolverOptions("decomp");
    01324                 //iterate over the vector
    01325
    01326                 int tmpVal;
    01327
    01328                 std::string routeString; //variable for parsing a category option
    01329                 std::string solutionString; //variable for parsing a category option
    01330                 string::size_type pos1; //variable for parsing a category option
    01331                 string::size_type pos2; //variable for parsing a category option
    01332                 string::size_type pos3; //variable for parsing a category option
    01333
    01334
    01335                 std::map<int, std::map<int, std::vector<int> > >::iterator  mit;
    01336                 std::map<int, std::vector<int> >::iterator  mit2;
    01337                 int solutionNumber;
    01338                 int routeNumber;
    01339
    01340
    01341                 for (vit = solverOptions.begin(); vit != solverOptions.end(); vit++) {
    01342
    01343
    01344                         //std::cout << (*vit)->name << std::endl;
    01345
    01346                         //(*vit)->name.compare("initialCol") == 0
    01347                         //if(rowNames[ i3].find("routeCapacity(1)") == string::npos )
    01348
    01349                         if( (*vit)->name.find("numHubs") !=  std::string::npos){
    01350
    01351
    01352                                 std::istringstream hubBuffer( (*vit)->value);
    01353                                 hubBuffer >> m_numHubs;
    01354                                 std::cout << "numHubs = " << m_numHubs <<  std::endl;
    01355
    01356                         }else{
    01357
    01358                                 if((*vit)->name.find("numNodes") !=  std::string::npos){
    01359
    01360
    01361                                         std::istringstream numNodesBuffer( (*vit)->value);
    01362                                         numNodesBuffer >> m_numNodes;
    01363                                         std::cout << "numNodes = " <<  m_numNodes <<  std::endl;
    01364
    01365                                 }else{
    01366                                         if((*vit)->name.find("totalDemand") !=  std::string::npos){
    01367
    01368
    01369                                                 std::istringstream totalDemandBuffer( (*vit)->value);
    01370                                                 totalDemandBuffer >> m_totalDemand;
    01371                                                 std::cout << "m_totalDemand = " << m_totalDemand <<  std::endl;
    01372
    01373                                         }else{
    01374                                                 if((*vit)->name.find("minDemand") !=  std::string::npos){
    01375
    01376
    01377                                                         std::istringstream minDemandBuffer( (*vit)->value);
    01378                                                         minDemandBuffer >> m_minDemand;
    01379                                                         std::cout << "m_minDemand = " << m_minDemand <<  std::endl;
    01380
    01381                                                 }else{
    01382                                                         if( (*vit)->name.find("demand") !=  std::string::npos ){
    01383
    01384
    01385                                                                 std::istringstream demandBuffer( (*vit)->value);
    01386                                                                 demandBuffer >> tmpVal;
    01387                                                                 demand.push_back( tmpVal);
    01388                                                                 //std::cout << "demand = " << tmpVal <<  std::endl;
    01389
    01390                                                         }else{
    01391                                                                 if((*vit)->name.find("routeCapacity") !=  std::string::npos ){
    01392                                                                         std::istringstream routeCapacityBuffer( (*vit)->value);
    01393                                                                         routeCapacityBuffer >> m_routeCapacity;
    01394                                                                         std::cout << "m_routeCapacity = " << m_routeCapacity <<  std::endl;
    01395
    01396                                                                 }else{
    01397
    01398                                                                         if((*vit)->name.find("osilFile") !=  std::string::npos ){
    01399                                                                                 m_initOSiLFile = (*vit)->value;
    01400                                                                                 std::cout << "m_initOSiLFile = " << m_initOSiLFile <<  std::endl;
    01401
    01402                                                                         }else{
    01403
    01404                                                                                 if( (*vit)->name.find("restrictedMasterSolution") !=  std::string::npos  ){
    01405                                                                                         //std::istringstream buffer( (*vit)->value);
    01406                                                                                         //buffer >> m_numberOfSolutions;
    01407
    01408                                                                                         //get the block number and solution number
    01409                                                                                         //first get routeString and soluionString
    01410                                                                                         //parse the category string base on :
    01411                                                                                         pos1 = (*vit)->category.find( ":");
    01412                                                                                         if(pos1 == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
    01413
    01414                                                                                         //solutionString = (*vit)->category.substr( pos1 + 1, pos2 - pos1 - 1);
    01415                                                                                         solutionString = (*vit)->category.substr( 0,  pos1);
    01416                                                                                         routeString = (*vit)->category.substr( pos1 + 1);
    01417
    01418                                                                                         pos2 = solutionString.find( "n");
    01419                                                                                         if(pos2  == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
    01420
    01421                                                                                         std::istringstream solutionBuffer( solutionString.substr( pos2 + 1)  );
    01422                                                                                         solutionBuffer >> solutionNumber;
    01423                                                                                         //std::cout << "solution number = " << solutionNumber  << std::endl;
    01424
    01425
    01426                                                                                         pos3 = routeString.find( "e");
    01427                                                                                         if(pos3  == std::string::npos ) throw ErrorClass("OSoL category attribute not properly defined");
    01428                                                                                         std::istringstream routeBuffer( routeString.substr( pos3 + 1)  );
    01429                                                                                         routeBuffer >> routeNumber;
    01430                                                                                         //std::cout << "route number = " <<  routeNumber << std::endl;
    01431                                                                                         std::istringstream nodeBuffer( (*vit)->value);
    01432                                                                                         nodeBuffer >> tmpVal;
    01433
    01434                                                                                         mit = m_initSolMap.find(  solutionNumber  );
    01435
    01436                                                                                         if( mit  != m_initSolMap.end() ){
    01437                                                                                                 // we have solution from before
    01438                                                                                                 // do we have a new route?
    01439
    01440                                                                                                 mit2 = mit->second.find( routeNumber);
    01441
    01442                                                                                                 if(mit2 != mit->second.end() ){
    01443                                                                                                 // we have a route from before and solution from before
    01444
    01445
    01446                                                                                                         mit2->second.push_back( tmpVal);
    01447
    01448
    01449                                                                                                 }else{
    01450
    01451                                                                                                         //we have a new route, but old solution
    01452                                                                                                         std::vector<int> tmpVec;
    01453                                                                                                         tmpVec.push_back( tmpVal) ;
    01454                                                                                                         mit->second.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
    01455
    01456
    01457                                                                                                 }
    01458
    01459                                                                                         }else{
    01460                                                                                                 // we have a new solution
    01461                                                                                                 std::vector<int> tmpVec;
    01462                                                                                                 tmpVec.push_back( tmpVal) ;
    01463
    01464                                                                                                 std::map<int, std::vector<int> > tmpMap;
    01465                                                                                                 tmpMap.insert( std::pair<int,std::vector<int> >(routeNumber, tmpVec) );
    01466                                                                                                 m_initSolMap.insert( std::pair<int, std::map<int, std::vector<int> > >(solutionNumber, tmpMap) )  ;
    01467
    01468                                                                                         }
    01469                                                                                 }
    01470                                                                         }
    01471                                                                 }
    01472                                                         }
    01473                                                 }
    01474                                         }
    01475                                 }
    01476                         }//end if on solver options
    01477
    01478                 }//end for loop on options
    01479
    01480                 //now fill in demand
    01481                 m_demand = new int[ m_numNodes];
    01482                 std::vector<int>::iterator vit2;
    01483                 if(m_numNodes != demand.size( ) ) throw ErrorClass("inconsistent number of demand nodes");
    01484                 int i;
    01485                 i = 0;
    01486                 for (vit2 = demand.begin(); vit2 != demand.end(); vit2++) {
    01487
    01488                         *(m_demand + i++) = *vit2;
    01489
    01490                 }
    01491                 demand.clear();
    01492
    01493                 //kipp -- fill in numberOfRestricedMasterSolutions from map size
    01494                  m_numberOfSolutions = m_initSolMap.size();
    01495
    01496
    01497         } catch (const ErrorClass& eclass) {
    01498
    01499                 throw ErrorClass(eclass.errormsg);
    01500
    01501         }
    01502
    01503 }//end getOptions
    01504
    01505
    01506 void OSRouteSolver::getCuts(const  double* x){
    01507
    01508
    01509 }//end getCuts
    01510
    01511 void OSRouteSolver::calcReducedCost( double** c, double* phi, double* d){
    01512
    01513         int k;
    01514         int i;
    01515         int j;
    01516         int l;
    01517         int kount;
    01518         kount = 0;
    01519         int tmpVal;
    01520         tmpVal = m_numNodes - 1;
    01521
    01522         for(k = 0; k < m_numHubs; k++){
    01523
    01524                 for(l = 0; l < m_upperBoundL; l++){
    01525
    01526
    01527                         for(i = 0; i< m_numNodes; i++){
    01528
    01529                                 //if we have (i, j) where j is hub then do not subtract off phi[ j]
    01530                                 for(j = 0; j < i; j++){
    01531
    01532                                         m_rc[ kount++] = (l + 1)*c[k][ i*tmpVal + j ] - phi[ j];
    01533
    01534                                 }
    01535
    01536                                 for(j = i + 1; j < m_numNodes; j++){
    01537
    01538                                         m_rc[ kount++] = (l + 1)*c[k][ i*tmpVal + j - 1 ] - phi[ j];
    01539
    01540                                 }
    01541
    01542
    01543                         }
    01544
    01545
    01546                 }
    01547
    01548
    01549         }
    01550
    01551 }//end calcReducedCost
    01552
    01553
    01554 void OSRouteSolver::createVariableNames( ){
    01555
    01556         int i;
    01557         int j;
    01558         int kount;
    01559
    01560         kount = 0;
    01561
    01562         for(i = 0; i< m_numNodes; i++){
    01563
    01564                 //if we have (i, j) where j is hub then do not subtract off phi[ j]
    01565                 for(j = 0; j < i; j++){
    01566
    01567                         m_variableNames[ kount] = makeStringFromInt("x(" , i);
    01568                         m_variableNames[ kount] += makeStringFromInt( "," , j);
    01569                         m_variableNames[ kount] +=  ")";
    01570                         //std::cout << "GAIL VARIABLE NAME " << m_variableNames[ kount] << std::endl;
    01571
    01572                         kount++;
    01573
    01574                 }
    01575
    01576                 for(j = i + 1; j < m_numNodes; j++){
    01577
    01578                         m_variableNames[ kount] = makeStringFromInt("x(" , i);
    01579                         m_variableNames[ kount] += makeStringFromInt( "," , j);
    01580                         m_variableNames[ kount] +=  ")";
    01581
    01582                         //std::cout << "GAIL VARIABLE NAME " << m_variableNames[ kount] << std::endl;
    01583                         kount++;
    01584
    01585                 }
    01586
    01587
    01588         }
    01589 }//end createVariableNames
    01590
    01591 void OSRouteSolver::createAmatrix(){
    01592
    01593         //arrays for the coupling constraint matrix
    01594         //int* m_pntAmatrix;
    01595         //int* m_Amatrix;
    01596
    01597
    01598         int i;
    01599         int j;
    01600
    01601         int numNonz;
    01602
    01603         //loop over nodes
    01604         m_pntAmatrix[ 0] = 0;
    01605         numNonz = 0;
    01606
    01607         for(j = m_numHubs; j < m_numNodes; j++){
    01608
    01609
    01610                 for(i = 0; i < j; i++){
    01611
    01612                         m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j - 1 ;
    01613
    01614                 }
    01615
    01616                 for(i = j + 1; i < m_numNodes; i++){
    01617
    01618                         m_Amatrix[ numNonz++] = i*(m_numNodes - 1) + j ;
    01619
    01620                 }
    01621
    01622                 m_pntAmatrix[ j - m_numHubs + 1]  = numNonz;
    01623
    01624         }
    01625
    01626         /*
    01627         for(i = 0; i < m_numNodes - m_numHubs; i++){
    01628
    01629                 for(j = m_pntAmatrix[ i]; j <  m_pntAmatrix[ i + 1]; j++){
    01630
    01631                         std::cout << m_variableNames[  m_Amatrix[ j ] ] << std::endl;
    01632
    01633                 }
    01634
    01635         }
    01636         */
    01637
    01638 }//end createAmatrix
    01639
    01640 void OSRouteSolver::pauHana(const double* theta){
    01641
    01642         std::cout <<  std::endl;
    01643         std::cout << "     PAU HANA TIME! " << std::endl;
    01644
    01645         int i;
    01646         int j;
    01647         int k;
    01648         double cost = 0;
    01649         for(i = 0; i < m_numThetaVar - 1  ; i++){
    01650
    01651                 cost += theta[ i]*m_thetaCost[ i + 1];
    01652                 //std::cout << "COLUMN VALUE = " << theta[ i] << std::endl;
    01653         }
    01654
    01655
    01656         std::cout << "FINAL SOLUTION VALUE = " << cost << std::endl;
    01657         std::cout << "NUMBER OF GENERATED COLUMNS = " << m_numThetaVar - 1 << std::endl;
    01658         std::cout << std::endl <<  std::endl;
    01659
    01660         float numSets;
    01661         int kount;
    01662
    01663         numSets = floor( double((m_numThetaVar - 1 ) / m_numHubs));
    01664
    01665
    01666         //kipp throw exception if number of columns not an even multiple of m_numHubs
    01667         kount = 0;
    01668
    01669
    01670         for(i = 0; i < numSets  ; i++){
    01671
    01672                 for(k = 0; k < m_numHubs; k++){
    01673
    01674                         if( theta[ kount ] > .001){
    01675
    01676                                 std::cout << "HUB = "  <<  k << "  THETA = " << kount << " = "  << theta[ kount] << std::endl;
    01677
    01678                                 for(j = m_thetaPnt[ kount ];  j <  m_thetaPnt[ kount + 1];  j++){
    01679
    01680                                         std::cout << "VARIABLE "  <<  m_variableNames[ m_thetaIndex[ j] ]   << std::endl;
    01681
    01682                                 }
    01683
    01684                         }//loop on if positive
    01685
    01686                         kount++;
    01687
    01688                 }//loop on hubs
    01689
    01690         }//loop on sets
    01691
    01692 }//end pauHana -- no pun intended
    01693
    01694
    01695 OSInstance* OSRouteSolver::getSeparationInstance(){
    01696
    01697
    01698
    01699
    01700         m_osinstanceSeparation = NULL;
    01701
    01702         //add linear constraint coefficients
    01703         //number of values will nodes.size() the coefficients in the node constraints
    01704         //plus coefficients in convexity constraints which is the number of varaibles
    01705         int kountNonz;
    01706         int kount;
    01707         int startsIdx;
    01708         //we build these on nodes that do not include the hubs
    01709         int numYvar = (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1);
    01710         int numVvar = m_numNodes - m_numHubs;
    01711         //the plus 1 is for the kludge row
    01712         int numCon = (m_numNodes - m_numHubs) + (m_numNodes - m_numHubs)*(m_numNodes - m_numHubs - 1)/2 + 1;
    01713         double *values = new double[ 2*numYvar + 2*numVvar];
    01714         int *indexes = new int[ 2*numYvar + 2*numVvar];
    01715         int *starts = new int[ numYvar + numVvar + 1];
    01716         starts[ 0] = 0;
    01717         startsIdx = 0;
    01718         startsIdx++;
    01719         kountNonz = 0;
    01720         int i;
    01721         int j;
    01722
    01723
    01724         std::string separationVarName;
    01725         std::string separationConName;
    01726
    01727         try {
    01728
    01729                 m_osinstanceSeparation = new OSInstance();
    01730
    01731                 //start building the separation instance
    01732
    01733                 m_osinstanceSeparation->setInstanceDescription("The Tour Breaking Separation Problem");
    01734
    01735
    01736                 // now the constraints
    01737                 m_osinstanceSeparation->setConstraintNumber( numCon);
    01738
    01739
    01740                 //add the node rows
    01741                 for( i =  0; i < m_numNodes - m_numHubs ; i++){
    01742
    01743                         m_osinstanceSeparation->addConstraint(i,  makeStringFromInt("nodeRow_", i+  m_numHubs ) , 0.0, 1.0, 0);
    01744
    01745                 }
    01746
    01747                 //add the variable rows rows
    01748
    01749                 int rowKounter;
    01750                 rowKounter = m_numNodes - m_numHubs;
    01751
    01752                 for(i = m_numHubs; i < m_numNodes; i++){
    01753
    01754
    01755
    01756                         for(j = i+1; j < m_numNodes; j++){
    01757                                 separationConName = makeStringFromInt("Row_(", i);
    01758                                 separationConName += makeStringFromInt(",", j);
    01759                                 separationConName += ")";
    01760
    01761                                 m_osinstanceSeparation->addConstraint(rowKounter++,  separationConName , 0, 0, 0);
    01762                         }
    01763
    01764                 }
    01765
    01766                 // the klude row so we have +/-1 in every column
    01767                 m_osinstanceSeparation->addConstraint(rowKounter++,  "kludgeRow" , 0, m_numNodes, 0);
    01768
    01769                 //  the variables
    01770                 m_osinstanceSeparation->setVariableNumber(  numYvar + numVvar);
    01771
    01772
    01773
    01774                 std::cout << "NUMBER OF VARIABLES SET = " << numYvar + numVvar << std::endl;
    01775                 //add the v variables
    01776                 for(i = 0; i < numVvar; i++){
    01777
    01778                         separationVarName = makeStringFromInt("v", i + m_numHubs);
    01779
    01780                         m_osinstanceSeparation->addVariable(i, separationVarName, 0, 1, 'C');
    01781
    01782                         values[ kountNonz ] = -1.0;
    01783                         indexes[ kountNonz ] = i;
    01784                         kountNonz++;
    01785
    01786                         values[ kountNonz ] = 1.0;
    01787                         indexes[ kountNonz ] = rowKounter - 1;
    01788                         kountNonz++;
    01789
    01790
    01791
    01792                         starts[ startsIdx++ ] = kountNonz;
    01793
    01794
    01795                 }
    01796                 //add the y variables
    01797                 kount = numVvar;
    01798
    01799                 int i1;
    01800                 int j1;
    01801                 int kountCon;
    01802
    01803                 //adjust for fact we don't use hub nodes
    01804                 i1 = i - m_numHubs;
    01805                 j1 = j - m_numHubs;
    01806
    01807                 kountCon = m_numNodes - m_numHubs;
    01808
    01809                 for(i1 = 0; i1 < m_numNodes - m_numHubs; i1++){
    01810
    01811
    01812
    01813                         for(j1 = i1 + 1; j1 < m_numNodes - m_numHubs; j1++){
    01814
    01815
    01816
    01817                                 separationVarName = makeStringFromInt("y(", i1 + m_numHubs);
    01818                                 separationVarName += makeStringFromInt(",", j1 + m_numHubs);
    01819                                 separationVarName += ")";
    01820                                 m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
    01821
    01822                                 values[ kountNonz ] = 1.0;
    01823                                 indexes[ kountNonz ] = i1;
    01824                                 kountNonz++;
    01825
    01826                                 values[ kountNonz ] = -1.0;
    01827                                 indexes[ kountNonz ] = kountCon ;
    01828                                 kountNonz++;
    01829
    01830                                 starts[ startsIdx++ ] = kountNonz;
    01831
    01832
    01833
    01834
    01835                                 separationVarName = makeStringFromInt("y(", j1 + m_numHubs);
    01836                                 separationVarName += makeStringFromInt(",", i1 + m_numHubs);
    01837                                 separationVarName += ")";
    01838                                 m_osinstanceSeparation->addVariable(kount++, separationVarName, 0, 1, 'C');
    01839
    01840                                 values[ kountNonz ] = 1.0;
    01841                                 indexes[ kountNonz ] = j1;
    01842                                 kountNonz++;
    01843
    01844                                 values[ kountNonz ] = -1.0;
    01845                                 indexes[ kountNonz ] = kountCon ;
    01846                                 kountNonz++;
    01847
    01848                                 starts[ startsIdx++ ] = kountNonz;
    01849
    01850
    01851                                 kountCon++;
    01852
    01853
    01854                         }
    01855
    01856                 }
    01857
    01858                 std::cout << "NUMBER OF VARIABLES ADDED = " << kount << std::endl;
    01859
    01860                 // now add the objective function
    01861                 m_osinstanceSeparation->setObjectiveNumber( 1);
    01862                 SparseVector *objcoeff = new SparseVector( numVvar);
    01863
    01864
    01865                 for(i = 0; i < numVvar; i++){
    01866
    01867                         objcoeff->indexes[ i] = i;
    01868                         objcoeff->values[ i] = 1.0;
    01869
    01870                 }
    01871
    01872
    01873
    01874
    01875
    01876                 m_osinstanceSeparation->addObjective(-1, "objfunction", "min", 0.0, 1.0, objcoeff);
    01877                 //now for the nonzeros
    01878                 //add the linear constraints coefficients
    01879                 m_osinstanceSeparation->setLinearConstraintCoefficients(kountNonz , true,
    01880                                 values, 0, kountNonz - 1,  indexes, 0, kountNonz - 1, starts, 0, startsIdx);
    01881
    01882
    01883
    01884                 std::cout << m_osinstanceSeparation->printModel(  ) << std::endl;
    01885                 //below is temporty see if we can setup as a Clp network problem
    01886             CoinPackedMatrix* matrix;
    01887             bool columnMajor = true;
    01888             double maxGap = 0;
    01889                 matrix = new CoinPackedMatrix(
    01890                 columnMajor, //Column or Row Major
    01891                 columnMajor? m_osinstanceSeparation->getConstraintNumber() : m_osinstanceSeparation->getVariableNumber(), //Minor Dimension
    01892                 columnMajor? m_osinstanceSeparation->getVariableNumber() : m_osinstanceSeparation->getConstraintNumber(), //Major Dimension
    01893                 m_osinstanceSeparation->getLinearConstraintCoefficientNumber(), //Number of nonzeroes
    01894                 columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->values : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->values, //Pointer to matrix nonzeroes
    01895                 columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->indexes : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->indexes, //Pointer to start of minor dimension indexes -- change to allow for row storage
    01896                 columnMajor? m_osinstanceSeparation->getLinearConstraintCoefficientsInColumnMajor()->starts : m_osinstanceSeparation->getLinearConstraintCoefficientsInRowMajor()->starts, //Pointers to start of columns.
    01897                 0,   0, maxGap );
    01898
    01899                 ClpNetworkMatrix network( *matrix);
    01900
    01901
    01902                 ClpSimplex  model;
    01903
    01904                 //if( m_osinstanceSeparation->getObjectiveMaxOrMins()[0] == "min") osiSolver->setObjSense(1.0);
    01905
    01906
    01907             model.setOptimizationDirection( 1);
    01908                 model.loadProblem( network, m_osinstanceSeparation->getVariableLowerBounds(),
    01909                                 m_osinstanceSeparation->getVariableUpperBounds(),
    01910                                 m_osinstanceSeparation->getDenseObjectiveCoefficients()[0],
    01911                                 m_osinstanceSeparation->getConstraintLowerBounds(), m_osinstanceSeparation->getConstraintUpperBounds()
    01912                 );
    01913
    01914
    01915
    01916                 //model.loadProblem(network, lowerColumn, upperColumn, objective,
    01917             //lower, upper);
    01918
    01919              model.factorization()->maximumPivots(200 + model.numberRows() / 100);
    01920              model.factorization()->maximumPivots(1000);
    01921              //model.factorization()->maximumPivots(1);
    01922              if (model.numberRows() < 50)
    01923                   model.messageHandler()->setLogLevel( 10);
    01924              //model.dual();
    01925                 model.primal();
    01926
    01927
    01928              if (model.numberRows() < 50)
    01929                   model.messageHandler()->setLogLevel(10);
    01930              //model.dual();
    01931                 model.primal();
    01932                 for(i = 0; i < numYvar + numVvar ; i++){
    01933                         std::cout <<   m_osinstanceSeparation->getVariableNames()[ i]   << " = " << model.getColSolution()[ i] << std::endl;
    01934                 }
    01935
    01936                 for(i = 0; i < rowKounter ; i++){
    01937                         std::cout <<   m_osinstanceSeparation->getConstraintNames()[ i]   << " = " << model.getRowPrice()[ i] << std::endl;
    01938                 }
    01939
    01940
    01941
    01942
    01943
    01944                 //exit( 1);
    01945
    01946                 //
    01947                 delete objcoeff;
    01948
    01949         }catch (const ErrorClass& eclass) {
    01950
    01951                 throw ErrorClass(eclass.errormsg);
    01952
    01953         }
    01954
    01955         return NULL;
    01956 }//end getSeparationInstance
    01957
    01958
    01959 std::string makeStringFromInt(std::string theString, int theInt){
    01960         ostringstream outStr;
    01961         outStr << theString;
    01962         outStr << theInt;
    01963         return outStr.str();
    01964 }//end makeStringFromInt
    01965
    01966
    01967