/home/coin/SVN-release/OS-2.3.3/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 

Generated on Sun Mar 20 03:06:15 2011 by  doxygen 1.4.7