/home/coin/SVN-release/OS-2.4.0/OS/applications/columnGen/code/OSColGenApp.cpp

Go to the documentation of this file.
00001 /* $Id: OSColGenApp.cpp 3038 2009-11-07 11:43:44Z kmartin $ */
00012 #include "OSColGenApp.h"
00013 #include "OSErrorClass.h" 
00014 #include "OSDataStructures.h"
00015 #include "OSBearcatSolverXij.h"
00016 #include "OSConfig.h"
00017 #include "OSResult.h" 
00018 #include "OSOption.h"
00019 #include "OSiLReader.h"        
00020 #include "OSiLWriter.h" 
00021 #include "OSoLReader.h"        
00022 #include "OSoLWriter.h" 
00023 #include "OSrLReader.h"          
00024 #include "OSrLWriter.h"      
00025 #include "OSInstance.h"  
00026 #include "OSFileUtil.h"  
00027 
00028 #include  "OSDecompSolverFactory.h"
00029 
00030 
00031 
00032 #ifdef COIN_HAS_COUENNE    
00033 #include "OSCouenneSolver.h"
00034 #endif
00035 
00036 #ifdef COIN_HAS_IPOPT    
00037 #include "OSIpoptSolver.h"
00038 #endif
00039 
00040 
00041 
00042 
00043 #include<vector>
00044 
00045 #include<map>
00046 #include <sstream>
00047 using std::ostringstream;
00048 
00049 
00050 
00051 
00052 OSColGenApp::OSColGenApp():
00053         m_osinstanceMaster(NULL) {
00054         std::cout << "INSIDE OSColGenApp CONSTRUCTOR" << std::endl;
00055         
00056 }//end OSColGenApp Constructor
00057 
00058 
00059 OSColGenApp::OSColGenApp(   OSOption *osoption) {
00060           std::cout << "INSIDE OSColGenApp CONSTRUCTOR" << std::endl;
00061           //std::cout << "the contructor things whichBlock = " << m_whichBlock<< std::endl;
00062           
00063           //get parameters-options
00064           //set default values:
00065           
00066           m_calledBranchAndBound = false;
00067           
00068         
00069           
00070           m_osDecompParam.nodeLimit = 1000;
00071           m_osDecompParam.columnLimit = 20000;
00072           m_osDecompParam.masterColumnResetValue = 5000;
00073           m_osDecompParam.zeroTol = .0001;
00074           m_osDecompParam.artVarCoeff = 1000000;
00075           m_osDecompParam.optTolPerCent = 0;
00076           
00077           m_osoption = osoption;
00078           //get the options for the OSDecompSolver
00079           getOptions( m_osoption);
00080           
00081           
00082           m_osinstanceMaster = NULL;
00083           m_osrouteSolver = NULL;
00084           
00085           m_osrouteSolver = NULL;
00086           m_factoryInit  = new OSDecompFactoryInitializer();
00087           std::cout << "CREATE THE FACTORY " << std::endl;
00088           m_osrouteSolver = OSDecompSolverFactory::factories[ "OSBearcatSolverXij"]->create();
00089           std::cout << "FINISH FACTORY CREATION " << std::endl;
00090           std::cout << "SET FACTORY OPTION " << std::endl;
00091           m_osrouteSolver->m_osoption =  m_osoption;
00092           std::cout << "FINISH SET FACTORY OPTION " << std::endl;
00093           //m_osrouteSolver = new OSBearcatSolverXij( osoption);
00094 
00095           //share the common parameters
00096           m_osrouteSolver->m_osDecompParam = m_osDecompParam;
00097           m_osrouteSolver->initializeDataStructures();
00098          
00099           
00100           //initialize the bounds
00101           m_zUB = OSDBL_MAX;
00102           m_zLB = -OSDBL_MAX;
00103           
00104           //set the column number
00105           m_numColumnsOld = 0;
00106             
00107 
00108 
00109 } //end OSColGenApp Constructor
00110 
00111 
00112 OSColGenApp::~OSColGenApp(){
00113         
00114         std::cout << "INSIDE ~OSColGenApp DESTRUCTOR" << std::endl;
00115 
00116         //kipp -- why doesn't m_osrouteSolver delete the 
00117         //m_osinstanceMaster object
00118         
00119         if( m_osinstanceMaster != NULL) delete  m_osinstanceMaster;
00120 
00121         if( m_osrouteSolver != NULL) delete  m_osrouteSolver;
00122 
00123         //finally delete the factories
00124         
00125         delete m_factoryInit;
00126         
00127 }//end ~OSColGenApp() destructor
00128 
00129 
00130 void OSColGenApp::getInitialRestrictedMaster( ){
00131 
00132         m_osinstanceMaster = m_osrouteSolver->getInitialRestrictedMaster( );
00133         //std::cout << m_osinstanceMaster->printModel(  ) << std::endl;
00134         
00135 }//end generateInitialRestrictedMaster
00136 
00137 
00138 void OSColGenApp::getCuts(const  double* thetaVar, const int numThetaVar,
00139                 int &numNewRows, int*  &numNonz, int** &colIdx,
00140                 double** &values, double* &rowLB, double* &rowUB) {
00141         
00142         m_osrouteSolver->getCutsTheta( thetaVar, numThetaVar,
00143                         numNewRows, numNonz, colIdx, values, rowLB, rowUB);
00144         
00145         m_calledBranchAndBound = false;
00146         //for now let's always get these cuts, hence the default false
00147         if(numNewRows == 0 && m_calledBranchAndBound == false
00148                         && m_osrouteSolver->m_numMultCuts <= m_osrouteSolver->m_multiCommodCutLimit) {
00149                 m_osrouteSolver->getCutsMultiCommod( thetaVar, numThetaVar,
00150                                 numNewRows, numNonz, colIdx, values, rowLB, rowUB);     
00151                 
00152                 
00153                 m_osrouteSolver->m_numMultCuts += numNewRows;
00154                 
00155                 //double lhs;
00156                 //for(int i = 0; i < numNewRows; i++){
00157                 //      lhs = 0;
00158                         //for(int j = 0; j < numNonz[ i]; j++){ 
00159                                         
00160                         //      lhs += m_si->getColSolution()[ colIdx[i][j] ]*values[i][j];
00161                         //      std::cout << " cut coefficient = " << values[i][j] << " theta value =  " << m_si->getColSolution()[ colIdx[i][j] ] << std::endl;
00162                                 
00163                         //}
00164                         
00165                         //std::cout << "LHS = " << lhs << std::endl;
00166                         
00167                 //}// loop over number of new  rows
00168                 
00169                 //exit( 1);
00170         }//end on if
00171         
00172         
00173 }//end getCuts
00174 
00175 
00176 void OSColGenApp::getColumns(const  double* yA, const int numARows,
00177                         const  double* yB, const int numBRows,
00178                         int &numNewColumns, int* &numNonz, double* &cost, 
00179                         int** &rowIdx, double** &values, double &lowerBound) {
00180 
00181         m_osrouteSolver->getColumns(yA, numARows,
00182                         yB, numBRows, numNewColumns, numNonz, 
00183                         cost, rowIdx, values,  lowerBound);
00184         
00185 
00186         
00187         
00188 }//end generateColumns
00189 
00190 void OSColGenApp::getOptions(OSOption *osoption) {
00191         
00192         
00193         //get any options relevant to OSColGenApp
00194         
00195         try{
00196 
00197                 std::vector<SolverOption*> solverOptions;
00198                 std::vector<SolverOption*>::iterator vit;
00199                 
00200                 solverOptions = osoption->getSolverOptions("OSDecompSolver");
00201                 if (solverOptions.size() == 0) throw ErrorClass( "options for OSDecompSolver not available");
00202                 
00203                 
00204                 for (vit = solverOptions.begin(); vit != solverOptions.end(); vit++) {
00205                         
00206                         if((*vit)->name.find("columnLimit") !=  std::string::npos){
00207                                 
00208                                 
00209                                 std::istringstream columnLimitBuffer( (*vit)->value);
00210                                 columnLimitBuffer >> m_osDecompParam.columnLimit;
00211                                 std::cout << "columnLimit = " <<  m_osDecompParam.columnLimit <<  std::endl;
00212                                 
00213                         }else{
00214                                 
00215                                 if( (*vit)->name.find("artVarCoeff") !=  std::string::npos ){
00216                                         
00217                                         std::istringstream artVarCoeffBuffer( (*vit)->value);
00218                                         artVarCoeffBuffer >> m_osDecompParam.artVarCoeff;
00219                                         std::cout << "artVarCoeff = " <<  m_osDecompParam.artVarCoeff <<  std::endl;
00220                                         
00221                                 }else{
00222                                         
00223                                         if( (*vit)->name.find("zeroTol") !=  std::string::npos){
00224                                                 
00225                                                 std::istringstream zeroTolBuffer( (*vit)->value);
00226                                                 zeroTolBuffer >> m_osDecompParam.zeroTol;
00227                                                 std::cout << "zeroTol = " <<  m_osDecompParam.zeroTol <<  std::endl;
00228                                                 
00229                                         }else{
00230                                                 
00231                                                 
00232                                                 if( (*vit)->name.find("nodeLimit") !=  std::string::npos){
00233                                                         
00234                                                         std::istringstream nodeLimitBuffer( (*vit)->value);
00235                                                         nodeLimitBuffer >> m_osDecompParam.nodeLimit;
00236                                                         std::cout << "nodeLimit = " <<  m_osDecompParam.nodeLimit <<  std::endl;
00237                                                         
00238                                                 }else{
00239                                                         
00240                                                         if( (*vit)->name.find("masterColumnResetValue") !=  std::string::npos){
00241                                                                 
00242                                                                 std::istringstream masterColumnResetValueBuffer( (*vit)->value);
00243                                                                 masterColumnResetValueBuffer >> m_osDecompParam.masterColumnResetValue;
00244                                                                 std::cout << "masterColumnResetValue = " <<  m_osDecompParam.masterColumnResetValue <<  std::endl;
00245                                                         }else{
00246                                                                 
00247                                                                 if( (*vit)->name.find("optTolPerCent") !=  std::string::npos){
00248                                                                         
00249                                                                         std::istringstream optTolPerCentBuffer( (*vit)->value);
00250                                                                         optTolPerCentBuffer >> m_osDecompParam.optTolPerCent;
00251                                                                         std::cout << "masterColumnResetValue = " <<  m_osDecompParam.optTolPerCent<<  std::endl;
00252                                                                 }
00253                                                         }
00254                                                 }
00255                                         }
00256                                 }
00257                         }
00258                 }
00259         
00260         } catch (const ErrorClass& eclass) {
00261 
00262                 throw ErrorClass(eclass.errormsg);
00263 
00264         }               
00265         
00266 }//end getOptions
00267 
00268 
00269 void OSColGenApp::solve(){
00270         
00271         
00272 
00273         int *cbasis = NULL;
00274         int *rbasis = NULL;
00275         int *new_cbasis = NULL;
00276         
00277         
00278         
00279         std::set<std::pair<int, double> >::iterator sit;
00280         std::vector<int>::iterator vit;
00281         std::map<int, int>::iterator mit;
00282         
00283         int numCols;
00284         int numRows;
00285         int i;
00286 
00287         //initialize upper bound
00288         m_zUB = m_osrouteSolver->m_bestIPValue;
00289         
00290         //initialize number of columns and nodes generated
00291         
00292         m_numColumnsGenerated = 0;
00293         m_numNodesGenerated = 0;
00294         std::cout << " m_zUB  " << m_zUB  << std::endl;
00295         
00296         try{
00297                   
00298                 // the solver
00299                   
00300                 m_solver = new CoinSolver();
00301                 
00302                 // the solver interface
00303                 
00304                 //kipp -- later have clp be an option
00305                 //I guess for now it must be an Osi solver
00306                 m_solver->sSolverName ="clp";
00307                 //std::cout << m_osinstanceMaster->printModel(  ) << std::endl;
00308                 m_solver->osinstance = m_osinstanceMaster;
00309                 
00310                 //pass options
00311                 m_solver->osoption = m_osoption;        
00312                 m_solver->buildSolverInstance();
00313 
00314                 
00315                 //get the solver interface
00316                 m_si = m_solver->osiSolver;
00317                 
00318                 m_yA = new double[m_osinstanceMaster->getConstraintNumber() ];
00319                 //kipp -- hard coding, come back and fix with option
00320                 //kipp -- do all of the newing in a separate routine
00321                 m_yB = new double[  m_osrouteSolver->m_maxBmatrixCon];
00322                 
00323                 
00324                 m_maxCols = m_osrouteSolver->m_maxMasterColumns;
00325                 m_maxRows = m_osrouteSolver->m_maxMasterRows;
00326                 
00327                 m_theta = new double[ m_maxCols];
00328 
00329                 
00330                 //get initial LP relaxation of master
00331                 solveRestrictedMasterRelaxation();
00332                 //exit( 1);
00333                 //get the solution vector
00334                 numCols = m_si->getNumCols();
00335                 numRows = m_si->getNumRows();
00336                 
00337                 //kipp -- just testing
00338                 cbasis = new int[ numCols];
00339                 rbasis = new int[ numRows ];
00340                 m_si->getBasisStatus( cbasis, rbasis);
00341                 
00342                 for(i = 0; i < numCols; i++){
00343                         //get the basic primal variables
00344                         if(cbasis[ i] == 1) m_zOptRootLP.push_back( i);
00345                         //get the LP relaxation
00346                         *(m_theta + i) = m_si->getColSolution()[i];     
00347                         
00348                         m_zRootLPx_vals.push_back( *(m_theta + i) );
00349                         
00351                         int j;
00352                         if( *(m_theta + i) > m_osDecompParam.zeroTol){
00353                                 std::cout <<  "x variables for column "  << i  << std::endl;
00354                                 for(j = m_osrouteSolver->m_thetaPnt[ i];  j < m_osrouteSolver->m_thetaPnt[ i + 1] ;  j++){      
00355                                         std::cout <<  m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[  j] ]  << " = "  <<  *(m_theta + i)  << std::endl;       
00356                                 }       
00357                         }                       
00359                         
00360                         
00361                 }
00362                 m_zLB =  m_si->getObjValue();
00363                 m_zRootLP = m_si->getObjValue();
00364                 //print LP value at node
00365                 std::cout <<  "optimal LP value at root node = "  <<  m_zLB << std::endl;
00366                 //get the optimal LP root solution
00367                 
00368                 
00369                 //exit( 1);
00370 
00371                 for ( sit = m_osrouteSolver->intVarSet.begin() ; 
00372                                 sit != m_osrouteSolver->intVarSet.end(); sit++ ){
00373                         
00374                         m_si->setInteger( sit->first);
00375 
00376                 }
00377                 
00378                 CbcModel model(  *m_si);
00379                 OsiSolverInterface *ipSolver = model.solver();
00380                 std::cout <<  "start solving master as integer program  "   << std::endl;
00381                 ipSolver->branchAndBound();
00382                 std::cout <<  "done solving master as integer program  "   << std::endl;
00383                 //CbcMain0(  model);    
00384                 //CbcMain1( 0, 0, model);
00385                 //kipp -- put in check to make sure we get an integer solution
00386                 if( ipSolver->getObjValue() < m_zUB) m_zUB = ipSolver->getObjValue() ;
00387                 
00388                 //get the solution
00389                 numCols = m_si->getNumCols();
00390                 
00391                 for(i = 0; i < numCols; i++){   
00392                         
00393                         //get the indexes of integer variables
00394                         if( model.getColSolution()[i] > m_osDecompParam.zeroTol){
00395                                 
00396                                 m_zOptIndexes.push_back( i) ;
00397 
00398                         }
00399                         
00400                 }
00401 
00402                 for ( sit = m_osrouteSolver->intVarSet.begin() ; 
00403                                 sit != m_osrouteSolver->intVarSet.end(); sit++ ){
00404                         
00405                         m_si->setContinuous( sit->first);
00406                         m_si->setColUpper( sit->first, sit->second);
00407                         
00408                 }
00409         
00410                 std::cout << "OPTIMAL LP VALUE = " << m_zLB << std::endl;
00411                 std::cout << "CURRENT BEST IP VALUE = " << m_zUB << std::endl;
00412                 
00414                 //start reset
00415                 /*
00416                 int tmpCols =  m_numColumnsGenerated;
00417                 resetMaster();
00418                 numCols = m_si->getNumCols();
00419                 new_cbasis = new int[ numCols  ];
00420                 for(i = 0; i < numCols; i++) new_cbasis[ i] = 3;
00421                 for (vit = m_zOptRootLP.begin(); vit != m_zOptRootLP.end(); vit++ ) new_cbasis[  *vit ] = 1;
00422                 m_si->setBasisStatus( new_cbasis, rbasis );
00423                 solveRestrictedMasterRelaxation();
00424                 //m_si->initialSolve();
00425                 //exit( 1);
00426                 for(i = 0; i < numCols; i++) *(m_theta + i) = m_si->getColSolution()[i];        
00427                 std::cout << "NUMBER OF NEW GENERATED COLUMNS =  "  << m_numColumnsGenerated - tmpCols << std::endl;
00428                 */
00429                 //end reset
00431                 
00432                 m_osrouteSolver->m_rootLPValue = m_zRootLP;
00433                 
00434                 //go into branch and bound
00435                 m_message = "";
00436                 std::cout << "START BRANCH AND BOUND =  "   << std::endl;
00437                 if(m_zLB + m_osDecompParam.zeroTol <  m_zUB) branchAndBound();
00438                 
00439                 //demand values
00440                 //m_osrouteSolver->m_demand;
00441                 
00442                 std::cout << "FINISH BRANCH AND BOUND =  "   << std::endl;
00443                 printTreeInfo();
00444                 m_osrouteSolver->m_bestLPValue = m_zLB;
00445                 m_osrouteSolver->m_bestIPValue = m_zUB; 
00446                 if(m_message == "") m_message = "********  WE ARE OPTIMAL  *******";
00447                 m_osrouteSolver->pauHana( m_zOptIndexes, m_zRootLPx_vals, 
00448                                 m_numNodesGenerated, m_numColumnsGenerated, m_message);
00449                 
00450                 
00451                 delete m_solver;
00452                 
00453                 delete[] m_yA;
00454                 m_yA = NULL;
00455                 
00456                 delete[] m_yB;
00457                 m_yB = NULL;
00458                 
00459                 delete[] m_theta;
00460                 m_theta = NULL;
00461 
00462                 
00463                 delete[] cbasis;
00464                 cbasis = NULL;
00465                 if(new_cbasis != NULL) delete[] new_cbasis;
00466                 new_cbasis = NULL;
00467                 delete[] rbasis;
00468                 rbasis = NULL;
00469         
00470 
00471         } catch (const ErrorClass& eclass) {
00472                 
00473                 delete m_solver;
00474                 
00475                 delete[] m_yA;
00476                 m_yA = NULL;
00477                 
00478                 delete[] m_yB;
00479                 m_yB = NULL;
00480                 
00481                 delete[] m_theta;
00482                 m_theta = NULL; 
00483                 
00484 
00485         
00486                 delete[] cbasis;
00487                 cbasis = NULL;
00488                 if(new_cbasis != NULL) delete[] new_cbasis;
00489                 new_cbasis = NULL;
00490                 delete[] rbasis;
00491                 rbasis = NULL;
00492 
00493                 throw ErrorClass(eclass.errormsg);
00494 
00495 }
00496         
00497 }//end solve
00498 
00499 
00500 void OSColGenApp::solveRestrictedMasterRelaxation( ){
00501         
00502         int i;
00503         int k;
00504         //we include convexity constraints in this number
00505         int numARows;
00506         // dimension y to number of nodes
00507         int numBRows;
00508         int numCols;
00509         
00510         //getColumns function call return parameters
00511         int numNewColumns;
00512         int* numNonz = NULL;
00513         double* cost = NULL; 
00514         int** rowIdx = NULL; 
00515         double** values = NULL ; 
00516         double lowerBound;
00517         //end of getColumns function call return parameters
00518         
00519         double collb; // kipp -- put in getColumns
00520         double colub; // kipp -- put in getColumns
00521         //all of our m_theta columns have a lower bound of 0 and upper bound of 1
00522         collb = 0.0;
00523         colub = 1.0;
00524         //kipp -- I would like to use OSDBL_MAX but Clp likes this better
00525         //double bigM = 1.0e24;
00526         double bigM  = m_osDecompParam.artVarCoeff;
00527         //getRows function call return parameters
00528         int numNewRows;
00529         int* numRowNonz = NULL;
00530         int** colIdx = NULL; 
00531         double** rowValues = NULL ; 
00532         double* rowLB;
00533         double* rowUB;
00534         //end of getRows function call return parameters        
00535         //art variables
00536 
00537         int rowArtIdx;
00538         double rowArtVal;
00539         
00540         bool isCutAdded;
00541         
00542         
00543 
00544         try{
00545                 numARows = m_osrouteSolver->m_numNodes;
00546                 
00547                 
00548                 isCutAdded = true;
00549                 
00550                 while(isCutAdded == true ){
00551                         
00552                         isCutAdded = false;
00553                         //start out loop on if cuts found
00554                         std::cout << "CALL Solve  " << " Number of columns =  " <<  m_si->getNumCols() <<  std::endl;
00555                         //we are going through OS here, m_solver is a CoinSolver object
00556                         //now solve
00557                         m_solver->solve();
00558                         //m_si->initialSolve();
00559                         std::cout << "Solution Status =  " << m_solver->osresult->getSolutionStatusType( 0 ) << std::endl;
00560                         //std::cout <<  m_solver->osrl << std::endl;
00561         
00562                         if(m_si->getNumRows() != numARows + m_osrouteSolver->m_numBmatrixCon ) {
00563                                 std::cout << "m_si->getNumRows() = " << m_si->getNumRows() << std::endl;
00564                                 std::cout << "numARows = " << numARows << std::endl;
00565                                 std::cout << "m_numBmatrixCon = " << m_osrouteSolver->m_numBmatrixCon << std::endl;
00566                                 throw ErrorClass("detect a row number inconsistency in solveRestrictedMasterRelaxation");
00567                         }
00568                         
00569                         
00570                         
00571                         if(m_si->getRowPrice() == NULL  ) 
00572                                 throw ErrorClass("problem getting dual values in solveRestrictedMasterRelaxation");
00573                         
00574                         
00575                         numBRows = m_osrouteSolver->m_numBmatrixCon;
00576                         
00577                         for(i = 0; i <  numARows; i++){
00578                                 
00579                                 *(m_yA + i) = m_si->getRowPrice()[ i];
00580                                 
00581                         }
00582                         
00583                         for(i = numARows; i <  numARows + numBRows; i++){
00584                                 
00585                                 *(m_yB + i - numARows) = m_si->getRowPrice()[ i];
00586                                 
00587                         }
00588                         
00589                         lowerBound = -1;
00590                         int loopKount = 0;
00593                         while(lowerBound < -m_osDecompParam.zeroTol && loopKount < m_osDecompParam.columnLimit){
00594                                 loopKount++;
00595                                 
00596                                 //kipp here is where the while loop goes
00597                                 //start while loop
00598                                 getColumns(m_yA, numARows, m_yB, numBRows, numNewColumns, 
00599                                                 numNonz, cost,  rowIdx, values,  lowerBound);
00600                                 
00601                                 std::cout << "Lower Bound = " <<  lowerBound   << std::endl;
00602 
00603                                 
00604                         
00605                                 for(k = 0; k < numNewColumns; k++){
00606                                         
00607                                         m_si->addCol( numNonz[ k], rowIdx[k], values[k],
00608                                                         collb, colub,  cost[ k]) ;
00609                                         
00610                                         
00611                                         m_numColumnsGenerated++;
00612                                         
00613                                 }
00614                         
00615                                 std::cout << " OBJ VALUE = " << m_si->getObjValue() << std::endl;
00616                                 
00617                                 std::cout << "m_zUB  = " <<  m_zUB << std::endl;
00618                                 
00619                                 if(lowerBound + m_si->getObjValue() > (1 - m_osDecompParam.optTolPerCent)*m_zUB) break;
00620                         
00621                                 std::cout << std::endl  << std::endl << std::endl;
00622                                 std::cout << "CALL Solve  " << " Number of columns =  " <<  m_si->getNumCols() <<  std::endl;
00623                                 m_solver->solve();
00624                                 //m_si->initialSolve();
00625                                 std::cout << "Solution Status =  " << m_solver->osresult->getSolutionStatusType( 0 ) << std::endl;
00626                                 std::cout << "Number of solver interface columns =  " <<  m_si->getNumCols()  << std::endl;
00627                                 //m_numNodes is number of artificial variables
00628                                 
00629                                 numCols = m_si->getNumCols();
00630                                 
00631                                 if( numCols != m_osrouteSolver->m_numThetaVar ) throw ErrorClass("number variables in solver not consistent with master");
00632                                 if( numCols + m_osrouteSolver->m_numHubs >= m_maxCols) {
00633 
00634                                         m_message = " *****  COLUMN LIMIT EXCEEDED -- INSIDE solveRestrictedMasterRelaxation ****";
00635                                         printTreeInfo();
00636                                         m_osrouteSolver->m_bestLPValue = m_zLB;
00637                                         m_osrouteSolver->m_bestIPValue = m_zUB; 
00638                                         m_osrouteSolver->pauHana( m_zOptIndexes, m_zRootLPx_vals,
00639                                                         m_numNodesGenerated, m_numColumnsGenerated, m_message);
00640                                         throw ErrorClass("we ran out of columns");
00641                                 }
00642                                 
00643                                 for(i = 0; i <  numARows; i++){
00644                                         
00645                                         *(m_yA + i) = m_si->getRowPrice()[ i];
00646                                         
00647                                 }
00648                                                         
00649                                 for(i = numARows; i <  numARows + numBRows; i++){
00650                                         
00651                                         *(m_yB + i - numARows) = m_si->getRowPrice()[ i];
00652                                         
00653                                 }
00654                                 
00655                         }//end while on column generation
00657                         
00658                         if( loopKount >=  m_osDecompParam.columnLimit)throw ErrorClass("we exceeded loop kount in column generation");
00659                         
00660                         //get a primal solution
00661                         numCols = m_si->getNumCols();
00662                         for(i=0; i < numCols; i++){
00663                                 *(m_theta + i) = m_si->getColSolution()[i];
00664                         }
00665         
00666                         
00667                         numNewRows = 0;
00668                         
00669                         //do not get cuts if LP relaxation worse than upper bound
00670                         if(m_si->getObjValue() < (1 - m_osDecompParam.optTolPerCent)*m_zUB) 
00671                                 getCuts(m_theta, numCols, numNewRows, numRowNonz, 
00672                                         colIdx,rowValues, rowLB, rowUB);
00673                 
00674                         
00675                         if( numNewRows >= 1 ){
00676                                 
00677                                 isCutAdded = true;
00678                                 
00679                                 for(i = 0; i < numNewRows; i++){
00680                                         
00681                                         m_si->addRow(numRowNonz[ i], colIdx[ i], rowValues[ i], rowLB[ i], rowUB[ i] ) ;
00682                                         
00683                                         //add two  variables for this row so we can never be infeasible
00684                                         //m_si->addCol( numNonz, rowIdx[k], values[k],
00685                                         //              collb, colub,  cost[ k]) ;
00686                                         
00687                                         //add the artificial variable for the UB                                        
00688                                         rowArtVal = -1.0;
00689                                         rowArtIdx = m_si->getNumRows() - 1;
00690                                         //m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
00691                                         //m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, 1, bigM);
00692                                         //add the artificial variable for the LB                                        
00693                                         rowArtVal = 1.0;
00694                                         //m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
00695                                         m_si->addCol(1, &rowArtIdx, &rowArtVal, 0, 1, bigM);
00696                                         m_numColumnsGenerated++;
00697                                                         
00698                                 }
00699                                 
00700                                 std::cout << std::endl;
00701                                 std::cout << "CUTS WERE ADDED CALL SOLVE" << std::endl;
00702                                 m_solver->solve();
00703                                 
00704                         } // end if on whether or not we added cuts
00705                         
00706 
00707 
00708                 
00709                 }//end while on isCutAdded
00710 
00711                 
00712                 
00713         } catch (const ErrorClass& eclass) {
00714 
00715                 throw ErrorClass(eclass.errormsg);
00716 
00717         }               
00718         
00719         
00720 }//end solveRestrictedMasterRelaxation
00721 
00722 
00723 bool OSColGenApp::isInteger( const double *thetaVar, const int numThetaVar, 
00724                 const double tol){
00725         
00726         
00727         bool isInt;
00728         isInt = true;
00729         int i;
00730         
00731         try{    
00732                 
00733                 for(i = 0; i < numThetaVar; i++){
00734                         
00735                         if( (thetaVar[ i] > tol) && (thetaVar[ i] < 1 -  tol) ){
00736                                 
00737                                 isInt = false;
00738                                 break;
00739                         }
00740                         
00741                 }
00742                 
00743                 return isInt;
00744                 
00745         } catch (const ErrorClass& eclass) {
00746         
00747                 throw ErrorClass(eclass.errormsg);
00748         
00749         }               
00750 
00751         
00752         
00753 }//end isInteger
00754 
00755 
00756 void OSColGenApp::printDebugInfo( ){
00757         
00758         int numCols;
00759         int numRows;
00760         std::set<std::pair<int, double> >::iterator sit;
00761         int i;
00762         
00763         numCols = m_si->getNumCols();   
00764         
00765         
00766         for(i = 0;  i < numCols; i++){
00767                 
00768                 std::cout << "PROCESSING THETA COLUMN " << i <<  "  value =  " <<  m_si->getColSolution()[i] << std::endl;
00769                 
00770                 for(int j = m_osrouteSolver->m_thetaPnt[ i]; j <  m_osrouteSolver->m_thetaPnt[ i + 1]; j++ ){
00771                         
00772                         //std::cout << m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[ j] ] << std::endl;
00773                         
00774                 }
00775         }
00776 
00777         numRows = m_si->getNumRows();
00778         
00779         for(i = m_osrouteSolver->m_numNodes;  i < numRows; i++){
00780                 
00781                 std::cout << "PROCESSING ROW " << i  << std::endl;
00782                 
00783                 for(int j = m_osrouteSolver->m_pntBmatrix[ i  -   m_osrouteSolver->m_numNodes]; j <  m_osrouteSolver->m_pntBmatrix[ i + 1 -  m_osrouteSolver->m_numNodes]; j++ ){
00784                         
00785                         //std::cout << m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_Bmatrix[ j] ] << std::endl;
00786                         
00787                 }
00788         }
00789         //check integer variables and upper bounds -- loop over integer variable set
00790 
00791         //
00792         for ( sit = m_osrouteSolver->intVarSet.begin() ; 
00793                         sit != m_osrouteSolver->intVarSet.end(); sit++ ){
00794                 
00795                 //std::cout << "Integer variable  " << sit->first << " Upper Bound = "  << sit->second  << std::endl;   
00796                 
00797         }       
00798 }//end printDebugInfo
00799 
00800 
00801 bool OSColGenApp::branchAndBound( ){
00802         
00803         m_calledBranchAndBound = true;
00804         
00809         std::map<int, int> varConMap;
00810         
00811         std::vector<OSNode*> nodeVec;
00812         std::vector<OSNode*>::iterator vit;
00813         
00814         
00815         std::map<int, OSNode*>::iterator mit;
00816         int bestNodeID;
00817         double bestNodeBound;
00818 
00819         
00820         OSNode *osnode = NULL;
00821         OSNode *osnodeLeftChild = NULL;
00822         OSNode *osnodeRightChild = NULL;
00823         
00824         bool bandbWorked;
00825         bandbWorked = true;
00826         int numCols;
00827         int rowIdx;
00828         rowIdx = 0;
00829         
00830         bool leftNodeCreated = false;
00831         bool rightNodeCreated = false;
00832         
00833 
00834 
00835         try{
00836                 
00837                 //get the solution
00838                 numCols = m_si->getNumCols();   
00839                 //kipp -- imporant this is now found earliear.
00840                 //for(i = 0; i < numCols; i++){ 
00841                 //      //get the LP relaxation
00842                 //      std::cout  << "theta = " <<   *(m_theta + i) << std::endl; 
00843                 //}
00844                 
00845                 //NOTE -- we must know theta here
00846                 
00847                 //create a branching cut 
00848                 createBranchingCut(m_theta, numCols, varConMap, rowIdx);
00849         
00850 
00851 
00853                         
00854                 osnodeLeftChild = createChild(osnode, varConMap, rowIdx, 1, 1);
00855                 if(osnodeLeftChild != NULL){
00856                         //finally set the nodeID
00857                         //and record parent ID
00858                         //m_numNodesGenerated++;
00859                         osnodeLeftChild->nodeID = m_numNodesGenerated;
00860                         osnodeLeftChild->parentID = 0;
00861                         //nodeVec.push_back( osnodeLeftChild);
00862                         m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeLeftChild->nodeID, osnodeLeftChild) );
00863 
00864                 }
00865                 
00867                 
00869                 
00870                 osnodeRightChild = createChild(osnode, varConMap, rowIdx, 0, 0);
00871                 if(osnodeRightChild != NULL){
00872                         //finally set the nodeID
00873                         //and record parent ID
00874                         //m_numNodesGenerated++;
00875                         osnodeRightChild->nodeID = m_numNodesGenerated;
00876                         osnodeRightChild->parentID = 0;
00877                         //nodeVec.push_back( osnodeRightChild);
00878                         m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeRightChild->nodeID, osnodeRightChild) );
00879                 }
00880                         
00882                 
00883                 // now loop
00884                 //kipp -- make this an option
00885 
00886                 std::cout << "ENTERING THE WHILE IN BRANCH AND BOUND" << std::endl;
00887                 std::cout << "m_numNodesGenerated = " <<  m_numNodesGenerated  << std::endl;
00888                 //while( (nodeVec.size() > 0) && (m_numNodesGenerated <= nodeLimit) ){
00889                 while( m_nodeMap.size() > 0 ){
00890                         
00891                         if(m_numNodesGenerated > m_osDecompParam.nodeLimit ){
00892                                 m_message = "******* NODE LIMIT EXCEEDED *******";
00893                                 return false;
00894                         }
00895                         
00896                         
00897                         if( m_si->getNumCols() > m_maxCols ){
00898                                 m_message = "******* COLUMN LIMIT EXCEEDED *******";
00899                                 return false;
00900                         }
00901                         
00902                         //kipp -- experimental
00903                         //m_osDecompParam.masterColumnResetValue = 3000;
00904                         //if( m_si->getNumCols() > 200000) {
00905                         if( (m_numColumnsGenerated - m_numColumnsOld)  > 
00906                                 m_osDecompParam.masterColumnResetValue)  {
00907                                 m_numColumnsOld = m_numColumnsGenerated;
00908                                 std::cout << "DOING A MASTER RESET IN BRANCH AND BOUND" << std::endl;
00909                                 std::cout << "NUMBER OF COLUMNS BEFORE RESET = " << m_si->getNumCols() << std::endl;
00910                                 resetMaster();
00911                                 std::cout << "NUMBER OF COLUMNS AFTER RESET = " << m_si->getNumCols() << std::endl;
00912                                 //int tmpCols =  m_numColumnsGenerated;
00913                                 //solveRestrictedMasterRelaxation();
00914                                 //std::cout << "NUMBER OF NEW GENERATED COLUMNS =  "  << m_numColumnsGenerated - tmpCols << std::endl;
00915                                 //exit( 1);
00916                         }
00917                         
00918                         leftNodeCreated = false;
00919                         rightNodeCreated = false;
00920                         //grab a node -- for now the last node, we do FIFO
00921                         //osnode =  nodeVec.back();
00922                         
00923                         //let's loop and find node with the largest nodeID -- this will
00924                         //corespond to fifo
00925                         
00926                         bestNodeID = -1;
00927                         bestNodeBound = OSDBL_MAX;
00928                         //mit->first is the the OSNode nodeID
00929                         //mit->second is an OSNode
00930                         for (mit = m_nodeMap.begin(); mit != m_nodeMap.end(); mit++ ){
00931                                 
00932                                 //FIFO criterions
00933                                 //if( mit->second->nodeID > bestNodeID) bestNodeID =  mit->second->nodeID;
00934                                 
00935                                 //Best Bound criterion
00936                                 if( mit->second->lpValue < bestNodeBound) {
00937                                         
00938                                         bestNodeBound = mit->second->lpValue;
00939                                         bestNodeID =  mit->first;
00940                                         //note same as:bestNodeID =  mit->second->nodeID;
00941                                         
00942                                         
00943                                 }
00944                                 
00945                         }
00946 
00947                         //get the node
00948                         mit = m_nodeMap.find( bestNodeID );
00949                         if(mit == m_nodeMap.end() ) throw ErrorClass("a node selection problem in branch and bound");
00950                         osnode = mit->second;
00951                         
00952                         if( osnode->lpValue < (1 - m_osDecompParam.optTolPerCent)*m_zUB - m_osDecompParam.zeroTol){
00953                         
00954                         
00955                                 //create a branching cut 
00956                                 std::cout << "CREATE A BRANCHING CUT " << std::endl;
00957                                 createBranchingCut(osnode->thetaIdx, osnode->theta, osnode->thetaNumNonz, 
00958                                                 varConMap, rowIdx);
00959                                 
00961                                 /*
00962                                 std::map<int, int>::iterator tmpit;
00963                                 for (tmpit = varConMap.begin() ; tmpit != varConMap.end(); tmpit++ ){
00964                                         
00965                                         std::cout << std::endl;
00966                                         std::cout << "LOOPING OVER VARIABLE " << m_osrouteSolver->m_variableNames[ tmpit->first ] << std::endl;
00967                                         std::cout << "ROW UB  = " << osnode->rowUB[  tmpit->second] << std::endl;
00968                                         std::cout << "ROW LB  = " << osnode->rowLB[  tmpit->second] << std::endl;
00969                                 kippster
00970                                 }
00971                                 */
00973                                 
00974                                 std::cout << "BEST NODE ID " << bestNodeID << std::endl;
00975                                 std::cout << "NODE LP VALUE =  " << osnode->lpValue << std::endl;
00976                                 //check for node consistency
00977                                 checkNodeConsistency( rowIdx, osnode);
00978                                 // create children
00979                                 //create the left node
00980                                 
00981                                 osnodeLeftChild = createChild(osnode, varConMap, rowIdx, 1, 1);
00982                                 if(osnodeLeftChild != NULL){
00983                                         //finally set the nodeID
00984                                         //and record parent ID
00985                                         //m_numNodesGenerated++;
00986                                         osnodeLeftChild->nodeID = m_numNodesGenerated;
00987                                         osnodeLeftChild->parentID = osnode->nodeID;
00988                                         leftNodeCreated = true;
00989                                 }
00990                                 
00991                                 //create the right node
00992                                 osnodeRightChild = createChild(osnode, varConMap, rowIdx, 0, 0);
00993                                 if(osnodeRightChild != NULL){
00994                                         //finally set the nodeID
00995                                         //and record parent ID
00996                                         //m_numNodesGenerated++;
00997                                         osnodeRightChild->nodeID = m_numNodesGenerated;
00998                                         osnodeRightChild->parentID = osnode->nodeID;
00999                                         rightNodeCreated = true;
01000                                 }
01001                                 
01002                                 
01003                                 //nodeVec.erase( nodeVec.end() - 1) ;
01004                                 m_nodeMap.erase( mit);
01005                                 delete osnode;
01006                                 
01007                                 //if( leftNodeCreated == true) nodeVec.push_back( osnodeLeftChild) ;
01008                                 //if( rightNodeCreated == true) nodeVec.push_back( osnodeRightChild) ;
01009                                 
01010                                 if( leftNodeCreated == true) 
01011                                         m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeLeftChild->nodeID, osnodeLeftChild) ) ;
01012                                 
01013                                 if( rightNodeCreated == true) 
01014                                         m_nodeMap.insert ( std::pair<int, OSNode*>(osnodeRightChild->nodeID, osnodeRightChild) ) ;
01015 
01016                                  
01017                         }else{
01018                                 
01019                                 //fathom node by virtue of the upper bound
01020                                 std::cout << "FATHAM BY UPPER BOUND " << std::endl;
01021                                 //nodeVec.erase( nodeVec.end() - 1) ;
01022                                 m_nodeMap.erase( mit);
01023                                 delete osnode;
01024 
01025                         }//end if on lp bound check
01026                         
01027                         //kipp -- critical reset upper and lower bounds 
01028                         //kipp don't forget to erase the parent node
01029                         
01030                 }//end the while
01031                 
01032 
01033                 
01034                 if(m_numNodesGenerated > 0){
01035                         
01036                         m_zLB = (1 - m_osDecompParam.optTolPerCent)*m_zUB;
01037                 }else{
01038                         
01039                         m_zLB = m_zUB;
01040                 }
01041                 
01042                 
01043                 //exit( 1);
01044         
01045                 return bandbWorked;
01046 
01047         } catch (const ErrorClass& eclass) {
01048                 
01049                 throw ErrorClass(eclass.errormsg);
01050 
01051         }               
01052 
01053 }// end branchAndBound
01054 
01055 OSNode* OSColGenApp::createChild(const OSNode *osnodeParent, std::map<int, int> &varConMap,
01056                 const int rowIdx, const double rowLB, const double rowUB){
01057         
01058         m_numNodesGenerated++;
01059         
01060         OSNode *osnodeChild;
01061         osnodeChild = NULL;
01062         int numRows;
01063         int numCols;
01064         
01065         int tmpColNum ;
01066         int tmpRowNum  ;
01067         
01068         std::map<int, int>::iterator mit;
01069 
01070 
01071         
01072         int i;
01073         int k;
01074         int childRowIdxNumNonz;
01075         childRowIdxNumNonz = 0;
01076 
01077         //we want to store the solution vector (theta space)
01078         //in sparse format
01079         int thetaNumNonz;
01080 
01081         
01082         try{
01083 
01084                 if(osnodeParent != NULL) childRowIdxNumNonz = osnodeParent->rowIdxNumNonz + 1;
01085                 else childRowIdxNumNonz = 1;
01086                         
01087                 //set upper and lower bounds correctly
01088                 //set the parent values 
01089                 if(osnodeParent != NULL){
01090                         for(i = 0; i < osnodeParent->rowIdxNumNonz; i++){
01091                                 
01092                                 
01093                                 m_si->setRowLower( osnodeParent->rowIdx[ i], osnodeParent->rowLB[ i]);
01094                                 m_si->setRowUpper( osnodeParent->rowIdx[ i], osnodeParent->rowUB[ i]);
01095                                 
01096                                 
01097                         }
01098                 }
01099                 //set the new value
01100                 m_si->setRowLower( rowIdx, rowLB);
01101                 m_si->setRowUpper( rowIdx, rowUB);
01102                 //now solve
01103                 
01104                 //print out the restricted master
01105 
01106                 //if(rowUB == 0) m_si->writeLp( "gailTest2" );
01107         
01108                 //exit( 1);
01109                 std::cout << "CALL SOLVE FROM CREATE CHILD "  << std::endl;
01110                 //kipp -- important, you really need to verify that an optimal solution was obtained!!!
01111                 if(osnodeParent != NULL){
01112                         
01113                         tmpColNum = m_si->getNumCols() ;
01114                         tmpRowNum = m_si->getNumRows() ;
01115                         int *tmpColParent = new int[ tmpColNum];
01116                         int *tmpRowParent = new int[ tmpRowNum ];
01117                         
01118                         for(k = 0; k < tmpColNum; k++){
01119                                 
01120                                 if( m_si->getObjCoefficients()[k]  >=  
01121                                         m_osDecompParam.artVarCoeff -  m_osDecompParam.zeroTol)
01122                                         tmpColParent[ k ] = 3;
01123                                 //else if( osnodeParent->reducedCostIdx.find(k) == osnodeParent->reducedCostIdx.end() )
01124                                 //      tmpColParent[ k ] = 3;
01125                                 else tmpColParent[ k] = 0;
01126                                 
01127                         }
01128                         
01129                         for(k = 0; k < tmpRowNum; k++){
01130                                 
01131                                 tmpRowParent[ k] = 0;
01132                         }
01133                         
01134                         
01135                 
01136                         //for(k = 0; k < osnodeParent->colBasisStatus.size(); k++){
01137                                 
01138                                 //make the basis status of artificial variables 3
01139                                 //that is, nonbasic at lower bound
01140                                 //if( m_si->getObjCoefficients()[ osnodeParent->colBasisStatus[k].first]  
01141                                 //    >=  m_osDecompParam.artVarCoeff -  m_osDecompParam.zeroTol)
01142                                 //      tmpColParent[osnodeParent->colBasisStatus[k].first  ] = 3;
01143                                 //else
01144                                         //tmpColParent[osnodeParent->colBasisStatus[k].first  ] = 
01145                                         //              osnodeParent->colBasisStatus[k].second;
01146                         //}
01147                         
01148                         m_si->setBasisStatus(tmpColParent, tmpRowParent);
01149                         solveRestrictedMasterRelaxation( );
01150                         
01151                         //kippster extra error checking
01152                         //kippster check on upper and lower bound
01153                         //we are in rowIdx and the theta here should correspond to the same xijk
01154                         //getRowActivity()
01155                         
01156                         /*
01157                         m_si->initialSolve();
01158                         if(m_si->getRowActivity()[ rowIdx] > rowLB ||
01159                                         m_si->getRowActivity()[ rowIdx] < rowUB ) {
01160                                         
01161                                         std::cout << "Row lower bound = " << rowLB << std::endl;
01162                                         std::cout << "Row upper bound = " << rowUB << std::endl;
01163                                         std::cout << "Row activity  = " << m_si->getRowActivity()[ rowIdx] << std::endl;
01164                                         throw ErrorClass( "Violating a branching cut UB and LB");
01165                                 
01166                         }
01167                         */              
01168                         
01169                         //first get the column index nonzero elements in row rowIdx
01170                         
01171                         
01172                         
01176                         
01177                         delete[] tmpColParent;
01178                         tmpColParent = NULL;
01179                         delete[] tmpRowParent;
01180                         tmpRowParent = NULL;
01181                         
01182                         //solveRestrictedMasterRelaxation( osnodeParent->colBasisStatus, 
01183                         //              osnodeParent->rowBasisStatus);
01184                 } else { 
01185                         
01186                         solveRestrictedMasterRelaxation( );
01187                 }
01188                 
01189                 
01190                 std::cout << std::endl << std::endl;
01191                 std::cout << "FINISH SOLVING THE MASTER "  << std::endl;
01192                 
01193 
01194                 //
01195                 //now reset the upper and lower bound
01196                 //set the parent values 
01197                 if( osnodeParent != NULL){
01198                         for(i = 0; i < osnodeParent->rowIdxNumNonz; i++){
01199                                 
01200                                 
01201                                 m_si->setRowLower( osnodeParent->rowIdx[ i], 0);
01202                                 m_si->setRowUpper( osnodeParent->rowIdx[ i], 1);
01203                                 
01204                                 
01205                         }
01206                 }
01207                 //reset the new value
01208                 m_si->setRowLower( rowIdx, 0);
01209                 m_si->setRowUpper( rowIdx, 1);
01210                 
01211                 // let's try and fathom the node
01212                 // if we are not as good a upper bound
01213                 // we fathom, if we are integer we fathom
01214                 std::cout << std::endl << std::endl;
01215                 std::cout << "MESSAGE: START CREATION OF A CHILD NODE" << std::endl;
01216                 std::cout << "LB " << rowLB  <<  " UB = " << rowUB << std::endl;
01217                 std::cout << "MESSAGE: LP RELAXATION VALUE OF POTENTIAL CHILD NODE  " << m_si->getObjValue() << std::endl;
01218                 std::cout << "MESSAGE: OPTIMALITY STATUS OF NODE IS " << m_si->isProvenOptimal() << std::endl;
01219                 
01220                 if( m_si->getObjValue() < (1 - m_osDecompParam.optTolPerCent)*m_zUB - m_osDecompParam.zeroTol && m_si->isProvenOptimal() == 1) {
01221                         // okay cannot fathom based on bound try integrality
01222                         std::cout << "MESSAGE: WE CANNOT FATHOM THE CHILD BASED ON UPPER BOUND " << std::endl;
01223                         numCols = m_si->getNumCols();
01224                         numRows = m_si->getNumRows();
01225                         thetaNumNonz = 0;
01226                         
01227                         for(i = 0; i < numCols; i++){   
01228                                 //get the LP relaxation
01229                                 *(m_theta + i) = m_si->getColSolution()[i];     
01230                                 if( *(m_theta + i) > m_osDecompParam.zeroTol) thetaNumNonz++;
01231                                 
01232                         }
01233                         if( isInteger( m_theta, numCols, m_osDecompParam.zeroTol) == true){
01234                                 //fathom by integrality
01235                                 std::cout << "MESSAGE:  WE HAVE AN INTEGRALITY FATHOM " << m_zUB << std::endl;
01236                                 if( m_zUB > m_si->getObjValue() ){
01237                                         
01238                                         m_zUB = m_si->getObjValue() ;
01239                                         //clear out out solution vector
01240                                         if( m_zOptIndexes.size() > 0) m_zOptIndexes.clear();
01241                                         
01242                                         for(i = 0; i < numCols; i++){
01243                                                 
01244                                                 if( *(m_theta + i) > m_osDecompParam.zeroTol) m_zOptIndexes.push_back( i) ;
01245                                                 
01246                                         }
01247                                 }
01248                                 
01249                         }else{
01250                                 //create the child node
01251                                 std::cout << "MESSAGE:  WE ARE CREATING A CHILD NODE WITH NUMBER COLUMNS = "<<  numCols << std::endl;
01252                                 osnodeChild = new OSNode(childRowIdxNumNonz,  thetaNumNonz );
01253                                 
01254                                 
01255                                 //set the basis
01256                                 /*
01257                                 tmpColNum = m_si->getNumCols() ;
01258                                 tmpRowNum = m_si->getNumRows() ;
01259                                 int *tmpColChild = new int[ tmpColNum];
01260                                 int *tmpRowChild = new int[ tmpRowNum ];
01261                                 
01262                                 m_si->getBasisStatus(tmpColChild, tmpRowChild);
01263                                 
01264                                 for(k = 0; k < tmpColNum; k++){
01265                                         
01266                                         osnodeChild->colBasisStatus.push_back(std::make_pair(k, tmpColChild[ k]) );
01267                                         
01268                                 }
01269                                 
01270                                 for(k = 0; k < tmpRowNum; k++){
01271                                         
01272                                         osnodeChild->rowBasisStatus.push_back(std::make_pair(k, tmpRowChild[ k]) );
01273                                         
01274                                 }
01275                                 
01276                                 delete[] tmpColChild;
01277                                 tmpColChild = NULL;
01278 
01279                                 delete[] tmpRowChild;
01280                                 tmpRowChild = NULL;
01281                                 */
01282                                 
01283                                 //now set bound arrays 
01284                                 if(osnodeParent == NULL){
01285                                         osnodeChild->rowIdx[ 0] = rowIdx;
01286                                         if(rowLB <=  m_osDecompParam.zeroTol) osnodeChild->rowLB[ 0] = 0;
01287                                                 else osnodeChild->rowLB[ 0] = 1;
01288                                         
01289                                         if(rowUB <=  m_osDecompParam.zeroTol) osnodeChild->rowUB[ 0] = 0;
01290                                                 else osnodeChild->rowUB[ 0] = 1;
01291                                         
01292                                         
01293                                 }else{
01294                                         //set old values
01295                                         for(i = 0; i < osnodeParent->rowIdxNumNonz; i++){
01296                                                 
01297                                                 osnodeChild->rowIdx[ i] = osnodeParent->rowIdx[ i];
01298                                                 osnodeChild->rowLB[ i] = osnodeParent->rowLB[ i];
01299                                                 osnodeChild->rowUB[ i] = osnodeParent->rowUB[ i];                               
01300                                                 
01301                                         }
01302                                         //set new value
01303                                         
01304                                         osnodeChild->rowIdx[ childRowIdxNumNonz - 1] = rowIdx;
01305                                         
01306                                         
01307                                         
01308                                         if(rowLB <=  m_osDecompParam.zeroTol) osnodeChild->rowLB[ childRowIdxNumNonz - 1 ] = 0;
01309                                                 else osnodeChild->rowLB[ childRowIdxNumNonz - 1 ] = 1;
01310                                                                                 
01311                                         if(rowUB <=  m_osDecompParam.zeroTol) osnodeChild->rowUB[ childRowIdxNumNonz - 1 ] = 0;
01312                                                 else osnodeChild->rowUB[ childRowIdxNumNonz - 1 ] = 1;
01313                                                                                 
01314                                         
01315                                         
01316                                 }
01317                                 //set child lp value
01318                                 osnodeChild->lpValue = m_si->getObjValue();
01319                                 //set theta
01320                                 thetaNumNonz = 0;
01321                                 for(i = 0; i < numCols; i++){
01322                                         
01323                                         if( *(m_theta + i) > m_osDecompParam.zeroTol){
01324                                                 
01325                                                 osnodeChild->thetaIdx[ thetaNumNonz] = i;
01326                                                 osnodeChild->theta[ thetaNumNonz] = *(m_theta + i);
01327                                                 
01328                                                 thetaNumNonz++;
01329                                                 //std::cout <<  "x variables for column "  << i  << std::endl;
01330                                                 //for(int j = m_osrouteSolver->m_thetaPnt[ i ];  j < m_osrouteSolver->m_thetaPnt[ i + 1] ;  j++){
01331                                                 //      std::cout <<  m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[  j] ]  << " = "  <<  *(m_theta + i) << std::endl;        
01332                                                 //}
01333                                         }
01334                                         
01335                                         //add the reduced costs
01336                                         
01337                                         if(m_si->getReducedCost()[i] < (m_zUB - osnodeChild->lpValue) ) osnodeChild->reducedCostIdx.insert( i);
01338                                         
01339                                 }
01340                         }//end else on isInteger
01341                 }// end if on upper bound test
01342 
01343                 std::cout << std::endl << std::endl;
01344                 return osnodeChild;
01345         
01346         } catch (const ErrorClass& eclass) {
01347         
01348 
01349                 throw ErrorClass(eclass.errormsg);
01350 
01351         }               
01352         
01353 }//end createChild
01354 
01355 
01356 void OSColGenApp::createBranchingCut(const int* thetaIdx, const double* theta, 
01357                 const int numThetaVar, std::map<int, int> &varConMap, int &rowIdx){
01358         
01359         int varIdx;
01360         int numNonz;
01361         int* indexes;
01362         double* values; 
01363         
01364 
01365         //for(int i = 0; i < numThetaVar; i++){
01366         //      std::cout <<  "x variables for column "  << thetaIdx[i]  << std::endl;
01367         //      for(int j = m_osrouteSolver->m_thetaPnt[ thetaIdx[ i] ];  j < m_osrouteSolver->m_thetaPnt[ thetaIdx[ i] + 1] ;  j++){
01368         //              std::cout <<  m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[  j] ]  << " = "  << theta[ i]  << std::endl;
01369         //      }       
01370         //}
01371         
01372         //kipp -- I would like to use OSDBL_MAX but Clp likes this better
01373         //double bigM  = 1.0e24;
01374         double bigM  = m_osDecompParam.artVarCoeff;
01375         double rowArtVal ;
01376         
01377         std::map<int, int>::iterator mit;
01378         
01379         //get the branching cut information
01380         
01381         //
01382         //for(int i = 0; i < numThetaVar; i++ ) std::cout << "theta idx = " << thetaIdx[ i] << "  theta = " << theta[ i] << std::endl;
01383 
01384         
01385         m_osrouteSolver->getBranchingCut(thetaIdx, theta, numThetaVar, 
01386                         varConMap, varIdx, numNonz, indexes,  values);
01387         
01388         
01389                 
01390         //std::cout << "varIDX = " << varIdx << std::endl;
01391         //std::cout << "numNonz1 = " << numNonz << std::endl;   
01392         
01393         
01394 
01395 
01396         //if numNonz is greater than zero:
01397         // 1) add add new variable to map -- at this point varConMap is empty
01398         // 2) add constraint then add to the formulation
01399         // 3) add  variables
01400         
01401 
01402         
01403         if( numNonz >0){
01404 
01405                 
01406                 //add the row
01407                 //make upper and lower bound 0 and 1 first 
01408                 m_si->addRow(numNonz, indexes, values, 0, 1) ;
01409                 
01410                 //add the artificial variables
01411                 //add the artificial variable for the UB                                        
01412                 rowArtVal = -1.0;
01413                 rowIdx = m_si->getNumRows() - 1;
01414                 
01415                 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
01416                 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
01417                 //add the artificial variable for the LB                                        
01418                 rowArtVal = 1.0;
01419                 m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
01420                 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
01421                 m_numColumnsGenerated++;
01422                 
01423                 //insert into map -- this is the first variable
01424                 varConMap.insert( std::pair<int, int>(varIdx , rowIdx) );
01425                 
01426                 m_rowIdxVarMap.insert( std::pair<int, int>(rowIdx , varIdx) );
01427                 
01428                 
01429                 
01430         } else{
01431                 //the variable varIdx is in the map
01432                 //get the constraint associated with this variable
01433                 //throw and exception if varIdx not a key
01434         
01435                 mit = varConMap.find( varIdx);
01436                 if( mit == varConMap.end() ) throw ErrorClass("in branchAndBound getBranchingCut() returned inconsistent value for varIdx");
01437                                 else rowIdx = mit->second;
01438 
01439 
01440         }//end if on numNonz
01441         
01442         
01443         
01444 }//end createBranchingCut Sparse
01445 
01446 
01447 
01448 void OSColGenApp::createBranchingCut(const double* theta, 
01449                 const int numThetaVar, std::map<int, int> &varConMap, int &rowIdx){
01450         
01451         int varIdx;
01452         int numNonz;
01453         int* indexes;
01454         double* values; 
01455         
01456         //kipp -- I would like to use OSDBL_MAX but Clp likes this better
01457         //double bigM  = 1.0e24;
01458         double bigM  = m_osDecompParam.artVarCoeff;
01459         double rowArtVal ;
01460         
01461         std::map<int, int>::iterator mit;
01462         
01463         //get the branching cut information
01464         m_osrouteSolver->getBranchingCut( theta, numThetaVar, 
01465                         varConMap, varIdx, numNonz, indexes,  values);
01466 
01467                 
01468         std::cout << "varIDX2 = " << varIdx << std::endl;
01469         std::cout << "numNonz2 = " << numNonz << std::endl;     
01470         
01471         
01472         //for(int i = 0; i < numNonz; i++){
01473         //      std::cout <<  "x variables for column "  << indexes[i]  << std::endl;
01474         //      for(int j = m_osrouteSolver->m_thetaPnt[ indexes[ i] ];  j < m_osrouteSolver->m_thetaPnt[ indexes[ i] + 1] ;  j++){
01476         //      }       
01477         //}
01478 
01479         //if numNonz is greater than zero:
01480         // 1) add add new variable to map -- at this point varConMap is empty
01481         // 2) add constraint then add to the formulation
01482         // 3) add artificial variables
01483         
01484         if( numNonz >0){
01485                 
01486                 //add the row
01487                 //make upper and lower bound 0 and 1 first 
01488                 m_si->addRow(numNonz, indexes, values, 0, 1) ;
01489                 
01490                 //add the artificial variables
01491                 //add the artificial variable for the UB                                        
01492                 rowArtVal = -1.0;
01493                 rowIdx = m_si->getNumRows() - 1;
01494                 
01495                 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
01496                 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
01497                 //add the artificial variable for the LB                                        
01498                 rowArtVal = 1.0;
01499                 
01500                 m_si->addCol(1, &rowIdx, &rowArtVal, 0, 1.0, bigM);
01501                 //m_si->addCol(1, &rowIdx, &rowArtVal, 0, OSDBL_MAX, bigM);
01502                 m_numColumnsGenerated++;
01503                 
01504                 //insert into map -- this is the first variable
01505                 varConMap.insert ( std::pair<int,int>(varIdx , rowIdx) );
01506                 m_rowIdxVarMap.insert( std::pair<int, int>(rowIdx , varIdx) );
01507                 
01508                 
01509                 
01510         } else{
01511                 //the variable varIdx is in the map
01512                 //get the constraint associated with this variable
01513                 //throw and exception if varIdx not a key
01514         
01515                 mit = varConMap.find( varIdx);
01516                 if( mit == varConMap.end() ) throw ErrorClass("in branchAndBound getBranchingCut() returned inconsistent value for varIdx");
01517                                 else rowIdx = mit->second;
01518 
01519 
01520         }//end if on numNonz
01521         
01522         
01523         
01524 }//end createBranchingCut Dense
01525 
01526 
01527 void  OSColGenApp::resetMaster(){
01528         
01529         //kipp -- temp stuff here delete later
01531         
01532         //std::map<int, int> inVars;
01533         std::map<int, int>::iterator mit;
01534         std::set<int>::iterator sit;
01535         std::vector<int>::iterator vit;
01536         std::vector<std::pair<int, int> >::iterator vit2;
01537         std::map<int, OSNode*>::iterator mit2;
01538         int i;
01539         int kount = 0;
01540 
01541         
01542         inVars.clear();
01543         
01544         try{
01545 
01546                 //first add the columns corresponding to the root node solution
01547                 
01548                 //for(i = 0; i < m_si->getNumCols(); i++){
01549                         
01550                         //if( m_si->getColSolution()[i] > m_osDecompParam.zeroTol) inVars.insert( std::pair<int, int>(i, kount++) );
01551                         //if( m_si->getObjCoefficients()[i] < 10000) inVars.insert( std::pair<int, int>(i, kount++) );
01552                 //}     
01553 
01554                 for(vit = m_zOptRootLP.begin() ; vit != m_zOptRootLP.end(); vit++){
01555                         
01556                         inVars.insert( std::pair<int, int>(*vit, kount++) );
01557                         
01558                 }
01559                 
01560                 
01561 
01562                 //next add the integer varaibles in the best known integer solution
01563                 for(vit = m_zOptIndexes.begin() ; vit != m_zOptIndexes.end(); vit++){
01564                         
01565                         if( inVars.find( *vit ) == inVars.end() ) inVars.insert( std::pair<int, int>(*vit, kount++) );
01566                         
01567                         //for(int k = m_osrouteSolver->m_thetaPnt[*vit]; k <  m_osrouteSolver->m_thetaPnt[*vit + 1]; k++){
01568                                 
01569                                 //std::cout << m_osrouteSolver-> m_variableNames[ m_osrouteSolver->m_thetaIndex[ k] ] <<  std::endl;
01570                                 
01571                         //}
01572                 }
01573                 
01574 
01575                 
01576                 //now loop over the nodes in the branch and bound tree
01577                 //kipp -- this is hardcoded play with it later
01578                 double tmpEps = OSDBL_MAX;
01579                 for(mit2 = m_nodeMap.begin(); mit2 !=m_nodeMap.end(); mit2++){
01580                         
01581                         std::cout << "NUMBER OF REDUCED COSTS = " << mit2->second->reducedCostIdx.size() << std::endl;
01582                         for(sit = mit2->second->reducedCostIdx.begin(); 
01583                                         sit != mit2->second->reducedCostIdx.end(); sit++){
01584                                 
01585                                 if( ( inVars.find( *sit ) == inVars.end() )   
01586                                                 && (m_si->getObjCoefficients()[*sit] < m_osDecompParam.artVarCoeff) 
01587                                                 && (m_si->getReducedCost()[*sit] < tmpEps )
01588                                 ) 
01589                                         inVars.insert( std::pair<int, int>(*sit, kount++) );
01590                         
01591                                 
01592                         }
01593                         
01594                         //insert the thetat variables
01595                         for(i = 0; i < mit2->second->thetaNumNonz; i++){
01596                                 
01597                                 if( inVars.find( mit2->second->thetaIdx[ i] ) == inVars.end() )
01598                                         
01599                                         inVars.insert( std::pair<int, int>(mit2->second->thetaIdx[ i], kount++) );
01600                                 
01601                         }
01602 
01603                 }
01604                 
01605         
01606                 
01607                 //for(mit = inVars.begin();  mit != inVars.end(); mit++)  std::cout << "mit->first " <<   mit->first << "  mit->second   " << mit->second   << std::endl;       
01608                 std::cout << "NUMBER OF COLUMNS =  " <<  inVars.size()  <<  std::endl;
01609                 std::cout << "CALLING osroute solver reset " << std::endl;
01610                 m_osrouteSolver->resetMaster( inVars, m_si );
01611                 //std::cout << "THE MAPPING AFTER A RESET: "   <<  std::endl;
01612                 //for(mit = inVars.begin();  mit != inVars.end(); mit++)  std::cout << "mit->first " <<   mit->first << "  mit->second   " << mit->second   << std::endl;
01613                         
01614                 //int numVars =   m_osrouteSolver->m_osinstanceMaster->getVariableNumber();
01615                 int numVars =   m_si->getNumCols();
01616                 double *tmpVals = NULL;
01617                 tmpVals = new double[ numVars];
01618                 
01619                 for(i = 0;  i < numVars; i++)  tmpVals[ i ] = 0;
01620                         
01621                 for(mit = inVars.begin(); mit != inVars.end(); mit++){
01622                         //tmpVals now point to old index values
01623                         tmpVals[ mit->second] = m_theta[ mit->first] ;
01624                 }
01625                 
01626                 
01627                 for(i = 0; i  < numVars; i++) m_theta[ i] = tmpVals[ i] ;
01628                 
01629                 //reset the nodes in the branch and bound tree
01630                 
01631                 for(mit2 = m_nodeMap.begin(); mit2 !=m_nodeMap.end(); mit2++){
01632 
01633                         
01634                         for(i = 0; i < mit2->second->thetaNumNonz; i++){
01635                                 
01636                                 //inVars[mit2->second->thetaIdx[ i] ]
01637                                 if( inVars.find( mit2->second->thetaIdx[ i] ) == inVars.end() ) throw ErrorClass("index problem in resetMaster");
01638 
01639                                 //kipp check to make sure we do not index an aritifical variable
01640                                 mit2->second->thetaIdx[ i] =  inVars[ mit2->second->thetaIdx[ i] ] ;
01641                                 //if( inVars.find( mit2->second->thetaIdx[ i] ) != inVars.end() ) 
01642                                 //      inVars.insert( std::pair<int, int>(mit2->second->thetaIdx[ i], kount++) );
01643                                 
01644                         }
01645                         
01646                         
01647                         for(vit2 = mit2->second->colBasisStatus.begin(); 
01648                                         vit2 != mit2->second->colBasisStatus.end(); vit2++){
01649                                 
01650                                 (*vit2).first =  inVars[ (*vit2).first ] ;
01651                                 
01652                                         
01653                         }
01654                         
01655                         //reset reduced cost indexes
01656                         std::set<int> tmpSet;
01657                         for(sit = mit2->second->reducedCostIdx.begin(); 
01658                                         sit != mit2->second->reducedCostIdx.end(); sit++){
01659                                 
01660                                 tmpSet.insert( inVars[ *sit ] );
01661                         }
01662                         
01663                         mit2->second->reducedCostIdx.clear();
01664                         
01665                         for(sit = tmpSet.begin(); sit != tmpSet.end(); sit++){
01666                                 
01667                                 //make sure that variable *sit is in the new reset master
01668                                 
01669                                 if( inVars.find( *sit) !=  inVars.end()  ) 
01670                                         mit2->second->reducedCostIdx.insert( *sit );
01671                         }
01672                         tmpSet.clear();
01673                         //end reset reduced cost indexes
01674 
01675                 }//end loop over nodes in tree -- mit2
01676                 
01677                 //reset the indexes of variables in the current integer incumbent
01678                 for(vit = m_zOptIndexes.begin() ; vit != m_zOptIndexes.end(); vit++) *vit = inVars[ *vit ];
01679                 
01680                 
01681                 //reset the indexes of variables in the root LP
01682                 for(vit = m_zOptRootLP.begin() ; vit != m_zOptRootLP.end(); vit++) *vit = inVars[ *vit ];
01683 
01685                 
01686                 delete m_solver;
01687                 m_osinstanceMaster = m_osrouteSolver->m_osinstanceMaster;
01688                 m_solver = new CoinSolver();
01689 
01690                 //kipp -- later have clp be an option
01691                 //I guess for now it must be an Osi solver
01692                 m_solver->sSolverName ="cbc";
01693                 //std::cout << m_osinstanceMaster->printModel(  ) << std::endl;
01694                 m_solver->osinstance = m_osrouteSolver->m_osinstanceMaster;
01695                 m_solver->osoption = m_osoption;
01696                 
01697                 m_solver->buildSolverInstance();
01698         
01699                 //get the solver interface
01700                 m_si = m_solver->osiSolver;
01701                 
01702                 if(m_si->getNumCols() !=  m_osrouteSolver->m_osinstanceMaster->getVariableNumber() ) 
01703                         throw ErrorClass("there is an inconsistency in the the model rebuid in resetMaster");
01704                 
01705                 std::cout << "OSINTANCE NUMBER OF COLUMNS = "  << m_osrouteSolver->m_osinstanceMaster->getVariableNumber() << std::endl;
01706                 std::cout << "OSINTANCE NUMBER OF ROWS = "  << m_osrouteSolver->m_osinstanceMaster->getConstraintNumber() << std::endl;
01707                 std::cout << "SOLVER INTERFACE NUMBER OF COLUMNS = "  << m_si->getNumCols() << std::endl;
01708                 std::cout << "SOLVER INTERFACE NUMBER OF ROWS = "  <<m_si->getNumRows() << std::endl;
01709                 
01710                 
01711                 //kipp this is a check, DO NOT do in production run
01712                 
01713                 
01714                 
01715                 double lpVal;
01716                 
01717                 for(mit2 = m_nodeMap.begin(); mit2 !=m_nodeMap.end(); mit2++){
01718 
01719                         lpVal = 0;
01720                         
01721                         for(i = 0; i < mit2->second->thetaNumNonz; i++){
01722                                 
01723                                 lpVal +=  m_si->getObjCoefficients()[  mit2->second->thetaIdx[ i] ]*mit2->second->theta[ i];
01724                                 
01725                                 
01726                         }
01727                         
01728                         if( ( lpVal - mit2->second->lpValue > m_osDecompParam.zeroTol )  ||   
01729                                         (mit2->second->lpValue - lpVal > m_osDecompParam.zeroTol )      ) throw ErrorClass( "uh oh, problem with node lp value" );
01730                         
01731                         //std::cout << "lpVal = " << lpVal <<  "  lpValue = " << mit2->second->lpValue << std::endl ;
01732                 }               
01733                 
01734                 
01735                 //end check
01736                 
01737                 //m_si->writeLp( "gailTest2" );
01738                 
01739                 delete[] tmpVals;
01740                 tmpVals = NULL;
01741         
01742         } catch (const ErrorClass& eclass) {
01743 
01744                 throw ErrorClass(eclass.errormsg);
01745 
01746         }       
01747 
01748         
01749 }//end resetMaster
01750 
01751 void OSColGenApp::printTreeInfo(){
01752         
01753         std::map<int, OSNode*>::iterator mit;
01754         int i;
01755         
01756         std::cout <<  std::endl << std::endl;
01757         
01758         std::cout << "NUMBER OF REMAINING DANGLING NODES  = " << m_nodeMap.size() << std::endl;
01759 
01760         if( m_nodeMap.size() > 0)  m_zLB = OSDBL_MAX; //find best LP value over dangling nodes
01761         
01762         for ( mit = m_nodeMap.begin() ; 
01763                         mit != m_nodeMap.end(); mit++ ){
01764                 
01765                 std::cout << "NODE ID VALUE = " << mit->second->nodeID << " " ;
01766                 std::cout << "  NODE LP VALUE = " << mit->second->lpValue << std::endl;
01767                 
01768                 for(i = 0; i < mit->second->rowIdxNumNonz; i++){
01769                         
01770                         std::cout << "CONSTRAINT =  " << mit->second->rowIdx[ i]  ;
01771                         std::cout << "  CONSTRAINT LB = " <<  mit->second->rowLB[ i]  ;
01772                         std::cout << "  CONSTRAINT UB = " <<  mit->second->rowUB[ i]  << std::endl;
01773                 }
01774                         
01775                 if(  mit->second->lpValue < m_zLB) m_zLB =  mit->second->lpValue;
01776 
01777                 
01778         }
01779         m_nodeMap.clear();
01780         
01781         
01782 }//end printTreeInfo
01783 
01784 
01785 void OSColGenApp::checkNodeConsistency( const int rowIdx, const OSNode *osnode){
01786         try{
01787                 if( osnode == NULL) return;
01788                 //we are going to throw an exception if we try to add a constraint to a node that is already there
01789                 std::set<int> indexSet;
01790                 int i;
01791                 int j;
01792                 int rowIdxNumNonz = 0;
01793                 int thetaNumNonz = 0;
01794                 rowIdxNumNonz = osnode->rowIdxNumNonz;
01795                 thetaNumNonz = osnode->thetaNumNonz;
01796                 std::map<int, double> varSumMap;
01797                 
01798                 std::cout << "MESSAGE: CHECKING FOR NODE CONSISTENCY CONSTRAINT" << std::endl;
01799                 
01800                 for(i = 0; i < thetaNumNonz; i++){
01801                         
01802                         
01803                         //loop over theta variables
01804                         std::cout << "theta idx " << osnode->thetaIdx[ i] << " theta value " << osnode->theta[ i] << std::endl;
01805                         
01806                         for(j = m_osrouteSolver->m_thetaPnt[ osnode->thetaIdx[ i]  ]; j <  m_osrouteSolver->m_thetaPnt[ osnode->thetaIdx[ i] + 1 ]; j++ ){
01807                                 
01808                                 if( varSumMap.find( m_osrouteSolver->m_thetaIndex[ j] )  == varSumMap.end() ){
01809                                         
01810                                         varSumMap[ m_osrouteSolver->m_thetaIndex[ j] ] = osnode->theta[ i];
01811                                         
01812                                 }else{
01813                                         
01814                                         varSumMap[ m_osrouteSolver->m_thetaIndex[ j] ] += osnode->theta[ i];
01815                                 }
01816                                 
01817                                 std::cout << "xijk idx " << m_osrouteSolver->m_thetaIndex[ j] << " variable name =  " <<  
01818                                                 m_osrouteSolver->m_variableNames[ m_osrouteSolver->m_thetaIndex[ j] ] << std::endl;
01819                                 
01820                         }
01821                 
01822                 }
01823                 
01824                 
01825 
01826                 for(i = 0; i < rowIdxNumNonz; i++){
01827                         
01828                         std::cout << " row number " << osnode->rowIdx[ i] << " LB = " << osnode->rowLB[ i] << " UB = " 
01829                                         << osnode->rowUB[ i] ;
01830                         
01831                         std::cout << " variable index = " << m_rowIdxVarMap[ osnode->rowIdx[ i] ] ;
01832                         
01833                         std::cout << " variable name = " << m_osrouteSolver->m_variableNames[ m_rowIdxVarMap[ osnode->rowIdx[ i] ] ] ;
01834                         
01835                         std::cout << " variable sum = " <<  varSumMap[ m_rowIdxVarMap[ osnode->rowIdx[ i] ]]    << std::endl ;
01836                                 
01837                         if(indexSet.find( osnode->rowIdx[ i] ) == indexSet.end() ){
01838                                 
01839                                 indexSet.insert(  osnode->rowIdx[ i] );
01840                                 
01841                         }else{
01842                                 
01843                                 
01844                                 throw ErrorClass( "We are trying to add an existing constraint to a node" );
01845                         }
01846                         
01847                         
01848                 
01849                 }
01850 
01851         } catch (const ErrorClass& eclass) {
01852 
01853         throw ErrorClass(eclass.errormsg);
01854 
01855         }       
01856 }//end checkNodeConsistency

Generated on Thu Sep 22 03:06:00 2011 by  doxygen 1.4.7